Merge branch 'feature/lynnt_lightcalorimetry'

Author: hgreenlee

sbncode/CAFMaker/CAFMakerParams.h

@@ -483,6 +483,12 @@ namespace caf
       "simdrift"
     };
 
+    Atom<art::InputTag> SimEnergyDepositLabel {
+      Name("SimEnergyDepositLabel"),
+      Comment("Label of input sim::SimEnergyDeposit objects."),
+      art::InputTag("ionandscint", "priorSCE","G4")
+    };
+
     Atom<bool> FillTrueParticles {
       Name("FillTrueParticles"),
       Comment("Whether to fill the rec.true_particles branch. The information on true particles"
@@ -634,6 +640,12 @@ namespace caf
       Comment("Label of CVN scores."),
       "cvn" 
     };
+
+    Atom<string> LightCaloLabel {
+      Name("LightCaloLabel"),
+      Comment("Label of light calorimetry producer"),
+      "lightcalo"
+    };
 
     Atom<std::string> fBlipTag { 
       Name("BlipTag"),

sbncode/CAFMaker/CAFMaker_module.cc

@@ -98,6 +98,7 @@
 #include "lardataobj/RecoBase/MCSFitResult.h"
 #include "lardataobj/RecoBase/Cluster.h"
 #include "lardataobj/AnalysisBase/MVAOutput.h"
+#include "lardataobj/Simulation/SimEnergyDeposit.h"
 
 #include "nusimdata/SimulationBase/MCFlux.h"
 #include "nusimdata/SimulationBase/MCTruth.h"
@@ -121,6 +122,7 @@
 #include "sbnobj/Common/Reco/OpT0FinderResult.h"
 #include "sbnobj/SBND/CRT/CRTVeto.hh"
 #include "sbnobj/Common/Reco/CorrectedOpFlashTiming.h"
+#include "sbnobj/Common/Reco/LightCalo.h"
 #include "sbnobj/SBND/Timing/TimingInfo.hh"
 #include "sbnobj/SBND/Timing/FrameShiftInfo.hh"
 
@@ -1650,6 +1652,36 @@ void CAFMaker::produce(art::Event& evt) noexcept {
     } // end for fm
   } // end for i (mctruths)
 
+  // get sim energy deposits if they're there
+  ::art::Handle<std::vector<sim::SimEnergyDeposit>> sed_handle;
+  GetByLabelStrict(evt, fParams.SimEnergyDepositLabel().encode(), sed_handle);
+
+  if (!isRealData && sed_handle.isValid()){
+    art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
+
+    srtruthbranch.dep.reserve(mctruths.size());
+    for (size_t n=0; n<mctruths.size();n++){
+      SRTrueDeposit init;
+      init.electrons = 0;
+      init.photons = 0;
+      init.energy = 0; 
+      srtruthbranch.dep.push_back(init);
+    }
+
+    for (sim::SimEnergyDeposit const& sed: *sed_handle){
+      const auto trackID = sed.TrackID();
+
+      art::Ptr<simb::MCTruth> mctruth = pi_serv->TrackIdToMCTruth_P(trackID);
+      auto it = std::find(mctruths.begin(), mctruths.end(), mctruth);
+      if (it == mctruths.end()) continue; 
+
+      auto idx = std::distance(mctruths.begin(), it);
+      srtruthbranch.dep.at(idx).energy    += sed.Energy()*1e-3; // GeV
+      srtruthbranch.dep.at(idx).photons   += sed.NumPhotons();
+      srtruthbranch.dep.at(idx).electrons += sed.NumElectrons();
+    }
+  }
+
   // get the number of events generated in the gen stage
   unsigned n_gen_evt = 0;
   for (const art::ProcessConfiguration &process: evt.processHistory()) {
@@ -2072,6 +2104,10 @@ void CAFMaker::produce(art::Event& evt) noexcept {
      if (fmCorrectedOpFlash.isValid())
       slcCorrectedOpFlash = fmCorrectedOpFlash.at(0);
 
+    art::FindOneP<sbn::LightCalo> foLightCalo =
+      FindOnePStrict<sbn::LightCalo>(sliceList,evt,
+        fParams.LightCaloLabel() + slice_tag_suff);
+    const sbn::LightCalo *slcLightCalo = foLightCalo.isValid()? foLightCalo.at(0).get() : nullptr;
 
     art::FindOneP<lcvn::Result> foCVNResult =
       FindOnePStrict<lcvn::Result>(sliceList, evt,
@@ -2334,6 +2370,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
     FillSliceCRUMBS(slcCRUMBS, recslc);
     FillSliceOpT0Finder(slcOpT0, recslc);
     FillSliceBarycenter(slcHits, slcSpacePoints, recslc);
+    FillSliceLightCalo(slcLightCalo, recslc);
     FillTPCPMTBarycenterMatch(barycenterMatch, recslc);
     FillCorrectedOpFlashTiming(slcCorrectedOpFlash, recslc);
     FillCVNScores(cvnResult, recslc);

sbncode/CAFMaker/CMakeLists.txt

@@ -25,6 +25,7 @@ art_make_library( LIBRARY_NAME sbncode_CAFMaker
                   caf_RecoUtils
                   lardataobj::AnalysisBase
                   lardataobj::RecoBase
+                  lardataobj::Simulation
                   larrecodnn::CVN_func
                   larcorealg::Geometry
                   larcore::Geometry_Geometry_service

sbncode/CAFMaker/FillReco.cxx

@@ -699,6 +699,17 @@ namespace caf
     }
   }
 
+  void FillSliceLightCalo(const sbn::LightCalo *lightcalo,
+                          caf::SRSlice &slice)
+  {
+    if (lightcalo != nullptr) {
+      slice.lightcalo.charge    = lightcalo->charge;
+      slice.lightcalo.light     = lightcalo->light;
+      slice.lightcalo.energy    = lightcalo->energy;
+      slice.lightcalo.bestplane = lightcalo->bestplane; 
+    }
+  }
+
   void FillSliceBarycenter(const std::vector<art::Ptr<recob::Hit>> &inputHits,
                            const std::vector<art::Ptr<recob::SpacePoint>> &inputPoints,
                            caf::SRSlice &slice)

sbncode/CAFMaker/FillReco.h

@@ -49,7 +49,7 @@
 #include "sbnobj/Common/PMT/Data/PMTBeamSignal.hh"
 #include "sbnobj/SBND/Timing/TimingInfo.hh"
 #include "sbnobj/SBND/Timing/FrameShiftInfo.hh"
-
+#include "sbnobj/Common/Reco/LightCalo.h"
 #include "nusimdata/SimulationBase/MCParticle.h"
 #include "nusimdata/SimulationBase/MCTruth.h"
 
@@ -111,6 +111,9 @@ namespace caf
   void FillSliceOpT0Finder(const std::vector<art::Ptr<sbn::OpT0Finder>> &opt0_v,
                            caf::SRSlice &slice);
 
+  void FillSliceLightCalo(const sbn::LightCalo *lightcalo,
+                          caf::SRSlice& slice);
+
   void FillSliceBarycenter(const std::vector<art::Ptr<recob::Hit>> &inputHits,
                            const std::vector<art::Ptr<recob::SpacePoint>> &inputPoints,
                            caf::SRSlice &slice);