add sim energy deposit info and light calo info

Author: lynnt20

sbncode/CAFMaker/CAFMakerParams.h

@@ -465,6 +465,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"
@@ -616,6 +622,12 @@ namespace caf
       Comment("Label of CVN scores."),
       "cvn" 
     };
+
+    Atom<string> LightCaloLabel {
+      Name("LightCaloLabel"),
+      Comment("Label of light calorimetry producer"),
+      "lightcalo"
+    };
 
   };
 }

sbncode/CAFMaker/CAFMaker_module.cc

@@ -97,6 +97,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"
@@ -119,6 +120,7 @@
 #include "sbnobj/Common/Reco/CRUMBSResult.h"
 #include "sbnobj/Common/Reco/OpT0FinderResult.h"
 #include "sbnobj/Common/Reco/CorrectedOpFlashTiming.h"
+#include "sbnobj/Common/Reco/LightCalo.h"
 #include "sbnobj/SBND/Timing/TimingInfo.hh"
 #include "sbnobj/SBND/Timing/FrameShiftInfo.hh"
 
@@ -1422,6 +1424,15 @@ void CAFMaker::produce(art::Event& evt) noexcept {
     art::fill_ptr_vector(simchannels, simchannel_handle);
   }
 
+  // get sim energy deposits if they're there
+  ::art::Handle<std::vector<sim::SimEnergyDeposit>> sed_handle;
+  GetByLabelStrict(evt, fParams.SimEnergyDepositLabel().encode(), sed_handle);
+
+  std::vector<art::Ptr<sim::SimEnergyDeposit>> seds; 
+  if (sed_handle.isValid()){
+    art::fill_ptr_vector(seds, sed_handle);
+  }
+
   art::Handle<std::vector<simb::MCFlux>> mcflux_handle;
   GetByLabelStrict(evt, std::string("generator"), mcflux_handle);
 
@@ -1614,6 +1625,34 @@ void CAFMaker::produce(art::Event& evt) noexcept {
     } // end for fm
   } // end for i (mctruths)
 
+
+  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 (size_t n_dep=0; n_dep < seds.size(); n_dep++){
+      auto sed = seds[n_dep];
+      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()) {
@@ -1995,6 +2034,13 @@ 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 = nullptr;
+    if (foLightCalo.isValid()) {
+      slcLightCalo = foLightCalo.at(0).get();
+    }
 
     art::FindOneP<lcvn::Result> foCVNResult =
       FindOnePStrict<lcvn::Result>(sliceList, evt,
@@ -2250,6 +2296,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

@@ -634,6 +634,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

@@ -47,7 +47,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"
 
@@ -108,6 +108,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);