Merge branch 'release/v10_05_00'

Author: kjplows

CMakeLists.txt

@@ -16,7 +16,7 @@
 cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 
 find_package(cetmodules 3.20.00 REQUIRED)
-project(sbncode VERSION 10.04.08 LANGUAGES CXX)
+project(sbncode VERSION 10.05.00 LANGUAGES CXX)
 
 message(STATUS "\n\n  ==========================   ${PROJECT_NAME}   ==========================")
 

sbncode/CAFMaker/CAFMakerParams.h

@@ -44,6 +44,10 @@ namespace caf
     Atom<bool> SaveGENIEEventRecord { Name("SaveGENIEEventRecord"),
       Comment("Whether to produce GENIE event record to the output file"), false
     };
+    
+    Atom<bool> OverrideRealData { Name("OverrideRealData"),
+	      Comment("when true, some algorithms (e.g. PoT count) treat events as MC rather than real data -- e.g. set it if the event is an overlay"), false
+    };
 
     Atom<float> PrescaleFactor { Name("PrescaleFactor"),
 	Comment("Factor by which to prescale unblind events"), 10

sbncode/CAFMaker/CAFMaker_module.cc

@@ -199,7 +199,8 @@ class CAFMaker : public art::EDProducer {
 
   std::string fSourceFile;
   std::uint32_t fSourceFileHash;
-
+ 
+  bool fOverrideRealData;
   bool fFirstInSubRun;
   unsigned int fIndexInFile = SRHeader::NoSourceIndex;
   bool fFirstBlindInSubRun;
@@ -401,6 +402,7 @@ class CAFMaker : public art::EDProducer {
   // Note: we will define isRealData on a per event basis in produce function [using event.isRealData()], at least for now.
 
   fCafFilename = fParams.CAFFilename();
+  fOverrideRealData = fParams.OverrideRealData();
   fFlatCafFilename = fParams.FlatCAFFilename();
 
   // Normally CAFMaker is run wit no output ART stream, so these go
@@ -847,9 +849,18 @@ void CAFMaker::beginSubRun(art::SubRun& sr) {
     std::cout << std::endl;
     abort();
   }
-
+  if ( fOverrideRealData ) {
+	// Expects a generator POT summary then...
+	if(auto pot_handle = sr.getHandle<sumdata::POTSummary>(fParams.GenLabel())){
+      		fSubRunPOT = pot_handle->totgoodpot;
+      		fTotalPOT += fSubRunPOT;
+    	}else{
+      		std::cout << "Did not find MC POT info under " << fParams.GenLabel() << std::endl;
+      		if(fParams.StrictMode()) abort();
+    	}
+  }else{
   if(bnb_spill){
-    FillExposure(*bnb_spill, fBNBInfo, fSubRunPOT);
+    		FillExposure(*bnb_spill, fBNBInfo, fSubRunPOT);
     fTotalPOT += fSubRunPOT;
 
     // Find the spill for each event and fill the event map:
@@ -908,10 +919,9 @@ void CAFMaker::beginSubRun(art::SubRun& sr) {
                 << std::endl;
       if(fParams.StrictMode()) abort();
     }
-
     // Otherwise, if one label is blank, maybe no POT was the expected result
   }
-
+ }
   std::cout << "POT: " << fSubRunPOT << std::endl;
 
   fFirstInSubRun = true;
@@ -1263,8 +1273,8 @@ void CAFMaker::produce(art::Event& evt) noexcept {
 
   bool const firstInFile = (fIndexInFile++ == 0);
 
-  // is this event real data?
-  bool isRealData = evt.isRealData();
+  // is this event real data? -- BH: if fOverrideRealData, treat it as MC. Otherwise, get the info from the art event.
+  bool isRealData = !fOverrideRealData && evt.isRealData();
 
   std::unique_ptr<std::vector<caf::StandardRecord>> srcol(
       new std::vector<caf::StandardRecord>);
@@ -1459,7 +1469,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
 
 	  int iparticle=0;
 	  genie::GHepParticle * p = 0;
-	  while( genie_rec->Particle(iparticle) != 0 ) {
+	  while( (genie_rec->Particle(iparticle) != 0) && (iparticle < 250) ) {
 	    p = genie_rec->Particle(iparticle);
 	    fGenieEvtRec_brStdHepPdg[iparticle] = p->Pdg();
 	    fGenieEvtRec_brStdHepStatus[iparticle] = (int) p->Status(); 
@@ -1654,26 +1664,50 @@ void CAFMaker::produce(art::Event& evt) noexcept {
 
   // Get all of the OpFlashes
   std::vector<caf::SROpFlash> srflashes;
+  if(fDet == kICARUS)
+  {
+    for (const std::string& pandora_tag_suffix : pandora_tag_suffixes) {
+      art::Handle<std::vector<recob::OpFlash>> flashes_handle;
+      GetByLabelStrict(evt, fParams.OpFlashLabel() + pandora_tag_suffix, flashes_handle);
+      // fill into event
+      if (flashes_handle.isValid()) {
+        const std::vector<recob::OpFlash> &opflashes = *flashes_handle;
+        int cryostat = ( pandora_tag_suffix.find("W") != std::string::npos ) ? 1 : 0;
 
-  for (const std::string& pandora_tag_suffix : pandora_tag_suffixes) {
-    art::Handle<std::vector<recob::OpFlash>> flashes_handle;
-    GetByLabelStrict(evt, fParams.OpFlashLabel() + pandora_tag_suffix, flashes_handle);
-    // fill into event
-    if (flashes_handle.isValid()) {
-      const std::vector<recob::OpFlash> &opflashes = *flashes_handle;
-      int cryostat = ( pandora_tag_suffix.find("W") != std::string::npos ) ? 1 : 0;
+        // get associated OpHits for each OpFlash
+        art::FindMany<recob::OpHit> findManyHits(flashes_handle, evt, fParams.OpFlashLabel() + pandora_tag_suffix);
 
-      // get associated OpHits for each OpFlash
-      art::FindMany<recob::OpHit> findManyHits(flashes_handle, evt, fParams.OpFlashLabel() + pandora_tag_suffix);
+        int iflash=0;
+        for (const recob::OpFlash& flash : opflashes) {
 
-      int iflash=0;
-      for (const recob::OpFlash& flash : opflashes) {
+          std::vector<recob::OpHit const*> const& ophits = findManyHits.at(iflash);
 
-        std::vector<recob::OpHit const*> const& ophits = findManyHits.at(iflash);
+          srflashes.emplace_back();
+          FillICARUSOpFlash(flash, ophits, cryostat, srflashes.back());
+          iflash++;
+        }
+      }
+    }
+  }
+  else if(fDet == kSBND)
+  {
+    std::vector<std::string> tpc_suffixes_sbnd = {"tpc0", "tpc1"};
 
-        srflashes.emplace_back();
-        FillOpFlash(flash, ophits, cryostat, srflashes.back());
-        iflash++;
+    for (size_t tpc=0; tpc<tpc_suffixes_sbnd.size(); tpc++) {
+      art::Handle<std::vector<recob::OpFlash>> flashes_handle;
+      GetByLabelStrict(evt, fParams.OpFlashLabel() + tpc_suffixes_sbnd[tpc], flashes_handle);
+      // fill into event
+      if (flashes_handle.isValid()) {
+        const std::vector<recob::OpFlash> &opflashes = *flashes_handle;
+        // get associated OpHits for each OpFlash
+        art::FindMany<recob::OpHit> findManyHits(flashes_handle, evt, fParams.OpFlashLabel() + tpc_suffixes_sbnd[tpc]);
+        int iflash=0;
+        for (const recob::OpFlash& flash : opflashes) {
+          std::vector<recob::OpHit const*> const& ophits = findManyHits.at(iflash);
+          srflashes.emplace_back();
+          FillSBNDOpFlash(flash, ophits, tpc, srflashes.back());
+          iflash++;
+        }
       }
     }
   }

sbncode/CAFMaker/FillReco.cxx

@@ -179,7 +179,7 @@ namespace caf
   }
 
 
-  void FillOpFlash(const recob::OpFlash &flash,
+  void FillICARUSOpFlash(const recob::OpFlash &flash,
                   std::vector<recob::OpHit const*> const& hits,
                   int cryo, 
                   caf::SROpFlash &srflash,
@@ -227,6 +227,41 @@ namespace caf
     }
   }
 
+  void FillSBNDOpFlash(const recob::OpFlash &flash,
+    std::vector<recob::OpHit const*> const& hits,
+    int tpc, 
+    caf::SROpFlash &srflash,
+    bool allowEmpty) {
+
+    srflash.setDefault();
+
+    srflash.time = flash.Time();
+    srflash.timewidth = flash.TimeWidth();
+
+    double firstTime = std::numeric_limits<double>::max();
+    for(const auto& hit: hits){
+    double const hitTime = hit->HasStartTime()? hit->StartTime(): hit->PeakTime();
+    if (firstTime > hitTime)
+    firstTime = hitTime;
+    }
+    srflash.firsttime = firstTime;
+    srflash.tpc = tpc;
+
+    srflash.totalpe = flash.TotalPE();
+    srflash.fasttototal = flash.FastToTotal();
+    srflash.onbeamtime = flash.OnBeamTime();
+
+    srflash.center.SetXYZ( -9999.f, flash.YCenter(), flash.ZCenter() );
+    srflash.width.SetXYZ( -9999.f, flash.YWidth(), flash.ZWidth() );
+
+    // Checks if ( recob::OpFlash.XCenter() != std::numeric_limits<double>::max() )
+    // See LArSoft OpFlash.h at https://nusoft.fnal.gov/larsoft/doxsvn/html/OpFlash_8h_source.html
+    if ( flash.hasXCenter() ) {
+    srflash.center.SetX( flash.XCenter() );
+    srflash.width.SetX( flash.XWidth() );
+    }
+  }
+
   std::vector<float> double_to_float_vector(const std::vector<double>& v)
   {
     std::vector<float> ret;

sbncode/CAFMaker/FillReco.h

@@ -217,11 +217,18 @@ namespace caf
                         caf::SRSBNDCRTTrack &srsbndcrttrack,
                         bool allowEmpty = false);
 
-  void FillOpFlash(const recob::OpFlash &flash,
+  void FillICARUSOpFlash(const recob::OpFlash &flash,
                   std::vector<recob::OpHit const*> const& hits,
                   int cryo,
                   caf::SROpFlash &srflash,
                   bool allowEmpty = false);
+
+  void FillSBNDOpFlash(const recob::OpFlash &flash,
+                  std::vector<recob::OpHit const*> const& hits,
+                  int tpc,
+                  caf::SROpFlash &srflash,
+                  bool allowEmpty = false);
+                  
   void FillCRTPMTMatch(const sbn::crt::CRTPMTMatching &match,
 		       caf::SRCRTPMTMatch &srmatch,
 		       bool allowEmpty = false);

sbncode/CMakeLists.txt

@@ -20,6 +20,7 @@ add_subdirectory(GeometryTools)
 add_subdirectory(CosmicID)
 add_subdirectory(DetSim)
 add_subdirectory(Cluster3D)
+add_subdirectory(HitFinder)
 
 # Supera
 #

sbncode/Calibration/TrackCaloSkimmer_module.cc

@@ -234,7 +234,7 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e)
 
   // Hits (ICARUS style)
   art::FindManyP<anab::T0> fmT0CRTHit(tracks, e, fCRTHitT0producer);
-  art::FindManyP<sbn::crt::CRTHitT0TaggingInfo> fmCRTHitT0TaggingInfo(PFParticleList, e, fCRTHitT0producer);
+  art::FindManyP<sbn::crt::CRTHitT0TaggingInfo> fmCRTHitT0TaggingInfo(tracks, e, fCRTHitT0producer);
 
   // Track - associated data
   art::FindManyP<recob::Track> fmTracks(PFParticleList, e, fTRKproducer);

sbncode/DetSim/AdjustSimForTrigger_module.cc

@@ -22,6 +22,7 @@
 #include "lardataobj/Simulation/AuxDetSimChannel.h"
 #include "lardataobj/Simulation/BeamGateInfo.h"
 #include "lardataobj/Simulation/SimEnergyDeposit.h"
+#include "lardataobj/Simulation/SimEnergyDepositLite.h"
 #include "lardataobj/Simulation/SimPhotons.h"
 
 #include <lardata/DetectorInfoServices/DetectorClocksService.h>
@@ -49,11 +50,13 @@ class AdjustSimForTrigger : public art::EDProducer {
   art::InputTag fInitAuxDetSimChannelLabel;
   art::InputTag fInitBeamGateInfoLabel;
   art::InputTag fInitSimEnergyDepositLabel;
+  art::InputTag fInitSimEnergyDepositLiteLabel;
   art::InputTag fInitSimPhotonsLabel;
   art::InputTag fInitWaveformLabel;
   bool fShiftAuxDetIDEs;
   bool fShiftBeamGateInfo;
   bool fShiftSimEnergyDeposits;
+  bool fShiftSimEnergyDepositLites;
   bool fShiftSimPhotons;
   bool fShiftWaveforms;
   double fAdditionalOffset;
@@ -66,29 +69,33 @@ AdjustSimForTrigger::AdjustSimForTrigger(fhicl::ParameterSet const& p)
   , fInitAuxDetSimChannelLabel(p.get<art::InputTag>("InitAuxDetSimChannelLabel", "undefined"))
   , fInitBeamGateInfoLabel{p.get<art::InputTag>("InitBeamGateInfoLabel", "undefined")}
   , fInitSimEnergyDepositLabel{p.get<art::InputTag>("InitSimEnergyDepositLabel", "undefined")}
+  , fInitSimEnergyDepositLiteLabel{p.get<art::InputTag>("InitSimEnergyDepositLiteLabel", "undefined")}
   , fInitSimPhotonsLabel{p.get<art::InputTag>("InitSimPhotonsLabel", "undefined")}
   , fInitWaveformLabel(p.get<art::InputTag>("InitWaveformLabel", "undefined"))
   , fShiftAuxDetIDEs{p.get<bool>("ShiftAuxDetIDEs", false)}
   , fShiftBeamGateInfo{p.get<bool>("ShiftBeamGateInfo", false)}
   , fShiftSimEnergyDeposits{p.get<bool>("ShiftSimEnergyDeposits", false)}
+  , fShiftSimEnergyDepositLites{p.get<bool>("ShiftSimEnergyDepositLites", false)}
   , fShiftSimPhotons{p.get<bool>("ShiftSimPhotons", false)}
   , fShiftWaveforms{p.get<bool>("ShiftWaveforms", false)}
   , fAdditionalOffset{p.get<double>("AdditionalOffset", 0.)}
 {
   if (!(fShiftSimEnergyDeposits || fShiftSimPhotons || fShiftWaveforms || fShiftAuxDetIDEs ||
-        fShiftBeamGateInfo)) {
+        fShiftBeamGateInfo || fShiftSimEnergyDepositLites)) {
     throw art::Exception(art::errors::EventProcessorFailure)
       << kModuleName << ": NO SHIFTS ENABLED!\n";
   }
   mf::LogInfo(kModuleName) << std::boolalpha << "SHIFTING AUXDETIDES? " << fShiftAuxDetIDEs << '\n'
                            << "SHIFTING BEAMGATEINFO? " << fShiftBeamGateInfo << '\n'
                            << "SHIFTING SIMENERGYDEPOSITS? " << fShiftSimEnergyDeposits << '\n'
+                           << "SHIFTING SIMENERGYDEPOSITLITES? " << fShiftSimEnergyDepositLites << '\n'
                            << "SHIFTING SIMPHOTONS? " << fShiftSimPhotons << '\n'
                            << "SHIFTING OPDETWAVEFORMS? " << fShiftWaveforms;
 
   if (fShiftAuxDetIDEs) { produces<std::vector<sim::AuxDetSimChannel>>(); }
   if (fShiftBeamGateInfo) { produces<std::vector<sim::BeamGateInfo>>(); }
   if (fShiftSimEnergyDeposits) { produces<std::vector<sim::SimEnergyDeposit>>(); }
+  if (fShiftSimEnergyDepositLites) { produces<std::vector<sim::SimEnergyDepositLite>>(); }
   if (fShiftSimPhotons) { produces<std::vector<sim::SimPhotons>>(); }
   if (fShiftWaveforms) { produces<std::vector<raw::OpDetWaveform>>(); }
 }
@@ -205,6 +212,24 @@ void AdjustSimForTrigger::produce(art::Event& e)
     }
     e.put(std::move(pSimEDeps));
   }
+  // and SimEnergyDepositLite's
+  if (fShiftSimEnergyDepositLites) {
+    auto const& simEDepLites =
+      e.getProduct<std::vector<sim::SimEnergyDepositLite>>(fInitSimEnergyDepositLiteLabel);
+
+    auto pSimEDepLites = std::make_unique<std::vector<sim::SimEnergyDepositLite>>();
+
+    for (auto const& inSimEDepLite : simEDepLites) {
+      double energy = inSimEDepLite.Energy(); 
+      geo::Point_t middlePos = inSimEDepLite.Position();
+      double middleTime = inSimEDepLite.Time() + timeShiftForTrigger_ns;
+      int ID = inSimEDepLite.TrackID();
+      
+      pSimEDepLites->emplace_back(energy, middlePos, middleTime, ID);
+    }
+
+    e.put(std::move(pSimEDepLites));
+  }
 
   // Repeat for sim::SimPhotons
   if (fShiftSimPhotons) {

sbncode/DetSim/FilterSimEnergyDeposits_module.cc

@@ -18,6 +18,7 @@
 #include "larcorealg/Geometry/BoxBoundedGeo.h"
 
 #include "lardataobj/Simulation/SimEnergyDeposit.h"
+#include "lardataobj/Simulation/SimEnergyDepositLite.h"
 
 #include <memory>
 
@@ -43,6 +44,7 @@ class FilterSimEnergyDeposits : public art::EDProducer {
   // Declare member data here.
 
   art::InputTag fInitSimEnergyDepositLabel;
+  art::InputTag fInitSimEnergyDepositLiteLabel;
   geo::BoxBoundedGeo fBox;
   static constexpr auto kModuleName = "FilterSimEnergyDeposit";
   double ShiftX(double z) const;
@@ -56,6 +58,7 @@ class FilterSimEnergyDeposits : public art::EDProducer {
 FilterSimEnergyDeposits::FilterSimEnergyDeposits(fhicl::ParameterSet const& p)
   : EDProducer{p}
   , fInitSimEnergyDepositLabel{p.get<art::InputTag>("InitSimEnergyDepositLabel", "undefined")}
+  , fInitSimEnergyDepositLiteLabel{p.get<art::InputTag>("InitSimEnergyDepositLiteLabel", "undefined")}
   , fBox{p.get<double>("P1_X"), p.get<double>("P2_X"),
 	 p.get<double>("P1_Y"), p.get<double>("P2_Y"),
 	 p.get<double>("P1_Z"), p.get<double>("P2_Z")}
@@ -66,6 +69,9 @@ FilterSimEnergyDeposits::FilterSimEnergyDeposits(fhicl::ParameterSet const& p)
   // Call appropriate produces<>() functions here.
   // Call appropriate consumes<>() for any products to be retrieved by this module.
   produces<std::vector<sim::SimEnergyDeposit>>();
+
+  
+  if (fInitSimEnergyDepositLiteLabel != "undefined") { produces<std::vector<sim::SimEnergyDepositLite>>(); }
 }
 
 
@@ -126,6 +132,27 @@ void FilterSimEnergyDeposits::produce(art::Event& e)
 					       origID));
   }
   e.put(std::move(pSimEDeps));
+
+  // also futz with the SimEnergyDepositLite's if configured to
+  if (fInitSimEnergyDepositLiteLabel == "undefined") return;
+
+  auto const& simEDepLites =
+    e.getProduct<std::vector<sim::SimEnergyDepositLite>>(fInitSimEnergyDepositLiteLabel);
+  auto pSimEDepLites = std::make_unique<std::vector<sim::SimEnergyDepositLite>>();
+  pSimEDepLites->reserve(simEDepLites.size());
+
+  for (auto const& inSimEDepLite : simEDepLites) {
+    geo::Point_t middlePos = inSimEDepLite.Position();
+    if(fBox.ContainsPosition(middlePos))
+      continue;
+    double energy = inSimEDepLite.Energy();
+    double middleTime = inSimEDepLite.Time();
+    int ID = inSimEDepLite.TrackID();
+
+    pSimEDepLites->emplace_back(energy, middlePos, middleTime, ID);
+  }
+
+  e.put(std::move(pSimEDepLites));
 }
 
 DEFINE_ART_MODULE(FilterSimEnergyDeposits)