Updates to Enable ICARUS Spring Production

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(); 

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)

sbncode/DetSim/fcl/icarus_simedepfilter.fcl

@@ -3,6 +3,7 @@ BEGIN_PROLOG
 simedepfilter_ind1gap: {
   module_type: "FilterSimEnergyDeposits"
   InitSimEnergyDepositLabel: "ionization"
+  InitSimEnergyDepositLiteLabel: "sedlite"
   P1_X: -400 #cm, create a box that spans the full x,y range. Not super-precise.
   P1_Y: -200 #cm, create a box that spans the full x,y range. Not super-precise.
   P1_Z: -2 #cm, one edge of the filtered gap based on the 37 mm total size of the wire gap
@@ -16,4 +17,4 @@ simedepfilter_ind1gap: {
 }
 
 
-END_PROLOG
+END_PROLOG