Implementing SimEnergyDeposit filter for ICARUS variation studies in release tag version v09_89_01_01

.github/pull_request_template.md

@@ -0,0 +1,9 @@
+### Description 
+Please provide a detailed description of the changes this pull request introduces. If available, also link to a docdb link where the issue/change have been presented on/discussed.
+
+- [X] Have you added a label? (bug/enhancement/physics etc.)
+- [X] Have you assigned at least 1 reviewer?
+- [X] Is this PR related to an open issue / project?
+- [X] Does this PR affect CAF data format? If so, please assign a CAF maintainer as additional reviewer.
+- [X] Does this PR require merging another PR in a different repository (such as sbnanobj/sbnobj etc.)? If so, please link it in the description.
+- [X] Are you submitting this PR on behalf of someone else who made the code changes? If so, please mention them in the description.

CMakeLists.txt

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

sbncode/CAFMaker/CAFMakerParams.h

@@ -41,6 +41,10 @@ namespace caf
       Comment("Provide a string to override the automatic filename."), ""
     };
 
+    Atom<bool> SaveGENIEEventRecord { Name("SaveGENIEEventRecord"),
+      Comment("Whether to produce GENIE event record to the output file"), false
+    };
+
     Atom<float> PrescaleFactor { Name("PrescaleFactor"),
 	Comment("Factor by which to prescale unblind events"), 10
     };

sbncode/CAFMaker/CAFMaker_module.cc

@@ -87,6 +87,8 @@
 #include "messagefacility/MessageLogger/MessageLogger.h"
 
 // LArSoft includes
+#include "larcore/Geometry/WireReadout.h"
+#include "larcore/Geometry/Geometry.h"
 #include "lardataobj/RecoBase/PFParticle.h"
 #include "lardataobj/RecoBase/Slice.h"
 #include "lardataobj/RecoBase/Track.h"
@@ -99,8 +101,6 @@
 #include "nusimdata/SimulationBase/MCNeutrino.h"
 #include "nusimdata/SimulationBase/GTruth.h"
 
-#include "fhiclcpp/ParameterSetRegistry.h"
-
 #include "sbnobj/Common/EventGen/MeVPrtl/MeVPrtlTruth.h"
 #include "sbnobj/Common/Reco/RangeP.h"
 #include "sbnobj/Common/SBNEventWeight/EventWeightMap.h"
@@ -116,6 +116,11 @@
 #include "sbnobj/Common/Reco/CRUMBSResult.h"
 #include "sbnobj/Common/Reco/OpT0FinderResult.h"
 
+// GENIE
+#include "Framework/EventGen/EventRecord.h"
+#include "Framework/Ntuple/NtpMCEventRecord.h"
+#include "nugen/EventGeneratorBase/GENIE/GENIE2ART.h"
+
 #include "canvas/Persistency/Provenance/ProcessConfiguration.h"
 #include "larcoreobj/SummaryData/POTSummary.h"
 
@@ -190,6 +195,7 @@ class CAFMaker : public art::EDProducer {
   std::string fFlatCafPrescaleFilename;
 
   std::string fSourceFile;
+  std::uint32_t fSourceFileHash;
 
   bool fFirstInSubRun;
   unsigned int fIndexInFile = SRHeader::NoSourceIndex;
@@ -232,6 +238,14 @@ class CAFMaker : public art::EDProducer {
   TTree* fFlatTreeb = 0;
   TTree* fFlatTreep = 0;
 
+  // GENIE EventRecord
+  genie::NtpMCEventRecord * fGenieEvtRec = 0;
+  TTree                   * fGenieTree = 0;
+  genie::NtpMCEventRecord * fFlatGenieEvtRec = 0;
+  TTree                   * fFlatGenieTree = 0;
+  bool fSaveGENIEEventRecord;
+  unsigned int fGenieEventCounter;
+
   flat::Flat<caf::StandardRecord>* fFlatRecord = 0;
   flat::Flat<caf::StandardRecord>* fFlatRecordb = 0;
   flat::Flat<caf::StandardRecord>* fFlatRecordp = 0;
@@ -374,6 +388,8 @@ class CAFMaker : public art::EDProducer {
   // setup volume definitions
   InitVolumes();
 
+  fSaveGENIEEventRecord = fParams.SaveGENIEEventRecord();
+
 }
 
 //......................................................................
@@ -536,9 +552,14 @@ void CAFMaker::InitVolumes() {
 //......................................................................
 CAFMaker::~CAFMaker()
 {
+
+  delete fGenieEvtRec;
+  delete fGenieTree;
   delete fRecTree;
   delete fFile;
 
+  delete fFlatGenieEvtRec;
+  delete fFlatGenieTree;
   delete fFlatRecord;
   delete fFlatTree;
   delete fFlatFile;
@@ -640,6 +661,13 @@ void CAFMaker::respondToOpenInputFile(const art::FileBlock& fb) {
   fFirstBlindInFile = true;
   fFirstPrescaleInFile = true;
   fSourceFile = inputBasename;
+  // Getting full hash
+  size_t fSourceFileHashFull = std::hash<std::string>{}(fSourceFile);
+  // truncate the full hash into a 32-bit integer
+  // This is required to be used as one of the ingredient of TTree::BuildIndex();
+  // it shifts the "major" value by 32-bit (https://root.cern/doc/master/classTTreeIndex.html#a08aac749ab22fd5c8ab792a0061a4b0f),
+  // so should be less than or equal to 32-bit
+  fSourceFileHash = static_cast<std::uint32_t>(fSourceFileHashFull);
 
 }
 
@@ -948,6 +976,7 @@ void CAFMaker::AddMetadataToFile(TFile* outfile, const std::map<std::string, std
 //......................................................................
 void CAFMaker::InitializeOutfiles()
 {
+
   if(fParams.CreateCAF()){
 
     mf::LogInfo("CAFMaker") << "Output filename is " << fCafFilename;
@@ -976,7 +1005,15 @@ void CAFMaker::InitializeOutfiles()
       AddEnvToFile(fFilep);
     }
 
-  }
+    if (fSaveGENIEEventRecord) {
+      fGenieTree = new TTree( "GenieEvtRecTree", "GenieEvtRecTree" );
+      fGenieTree->Branch("GenieEvtRec", &fGenieEvtRec);
+      fGenieTree->Branch("GENIEEntry", &fGenieEventCounter, "GENIEEntry/i");
+      fGenieTree->Branch("SourceFileHash", &fSourceFileHash, "SourceFileHash/i");
+    }
+
+  }     
+
 
   if(fParams.CreateFlatCAF()){
     mf::LogInfo("CAFMaker") << "Output flat filename is " << fFlatCafFilename;
@@ -1011,8 +1048,17 @@ void CAFMaker::InitializeOutfiles()
       AddEnvToFile(fFlatFilep);
     }
 
+    if (fSaveGENIEEventRecord){
+      fFlatGenieTree = new TTree( "GenieEvtRecTree", "GenieEvtRecTree" );
+      fFlatGenieTree->Branch("GenieEvtRec", &fFlatGenieEvtRec);
+      fFlatGenieTree->Branch("GENIEEntry", &fGenieEventCounter, "GENIE/i");
+      fFlatGenieTree->Branch("SourceFileHash", &fSourceFileHash, "SourceFileHash/i");
+    }
+
   }
 
+  fGenieEventCounter = 0;
+
   fFileNumber = -1;
   fTotalPOT = 0;
   fSubRunPOT = 0;
@@ -1256,7 +1302,9 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
   auto const dprop =
     art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clock_data);
-  const geo::GeometryCore *geometry = lar::providerFrom<geo::Geometry>();
+  const geo::GeometryCore* geom = lar::providerFrom<geo::Geometry>();
+  const geo::WireReadoutGeom &wireReadout =
+    art::ServiceHandle<geo::WireReadout>()->Get();
 
   auto const *sce = lar::providerFrom<spacecharge::SpaceChargeService>();
 
@@ -1284,7 +1332,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   if ( !isRealData ) {
     art::ServiceHandle<cheat::BackTrackerService> bt_serv;
 
-    id_to_ide_map = PrepSimChannels(simchannels, *geometry);
+    id_to_ide_map = PrepSimChannels(simchannels, wireReadout);
     id_to_truehit_map = PrepTrueHits(hits, clock_data, *bt_serv);
     id_to_hit_energy_map = SetupIDHitEnergyMap(hits, clock_data, *bt_serv);
   }
@@ -1329,10 +1377,33 @@ void CAFMaker::produce(art::Event& evt) noexcept {
       std::cout << "Failed to get GTruth object!" << std::endl;
     }
 
+    //GENIE EventRecord
+
     srtruthbranch.nu.push_back(SRTrueInteraction());
     srtruthbranch.nnu ++;
 
-    if ( !isRealData ) FillTrueNeutrino(mctruth, mcflux, gtruth, true_particles, id_to_truehit_map, srtruthbranch.nu.back(), i, fActiveVolumes);
+    if ( !isRealData ){
+
+      FillTrueNeutrino(mctruth, mcflux, gtruth, true_particles, id_to_truehit_map, srtruthbranch.nu.back(), i, fActiveVolumes);
+
+      srtruthbranch.nu.back().genie_evtrec_idx = fGenieEventCounter;
+
+      // GENIE event record
+      if(fSaveGENIEEventRecord){
+        genie::EventRecord* genie_rec = evgb::RetrieveGHEP(*mctruth, gtruth);
+        if(fGenieTree){
+          fGenieEvtRec->Fill(fGenieEventCounter, genie_rec);
+          fGenieTree->Fill();
+        }
+        if(fFlatGenieTree){
+          fFlatGenieEvtRec->Fill(fGenieEventCounter, genie_rec);
+          fFlatGenieTree->Fill();
+        }
+      }
+
+      fGenieEventCounter++;
+
+    }
 
     // Don't check for syst weight assocations until we have something (MCTruth
     // corresponding to a neutrino) that could plausibly be reweighted. This
@@ -1481,22 +1552,21 @@ void CAFMaker::produce(art::Event& evt) noexcept {
       }
     }
 
-  // Get all of the CRTPMT Matches ..
+  // Get all of the CRTPMT Matches
   std::vector<caf::SRCRTPMTMatch> srcrtpmtmatches;
-  std::cout << "srcrtpmtmatches.size = " << srcrtpmtmatches.size() << "\n";
   art::Handle<std::vector<sbn::crt::CRTPMTMatching>> crtpmtmatch_handle;
   GetByLabelStrict(evt, fParams.CRTPMTLabel(), crtpmtmatch_handle);
   if(crtpmtmatch_handle.isValid()){
-    std::cout << "valid handle! label: " << fParams.CRTPMTLabel() << "\n";
     const std::vector<sbn::crt::CRTPMTMatching> &crtpmtmatches = *crtpmtmatch_handle;
     for (unsigned i = 0; i < crtpmtmatches.size(); i++) {
+      int topen = 0, topex = 0, sideen = 0, sideex = 0; // bottomen = 0, bottomex = 0;
       srcrtpmtmatches.emplace_back();
-      FillCRTPMTMatch(crtpmtmatches[i],srcrtpmtmatches.back());
+      FillCRTPMTMatch(crtpmtmatches[i],
+		      topen, topex, sideen, sideex, 
+		      //bottomen, bottomex,
+		      srcrtpmtmatches.back());
     }
   }
-  else{
-    std::cout << "crtpmtmatch_handle.isNOTValid!\n";
-  }
 
   // Get all of the OpFlashes
   std::vector<caf::SROpFlash> srflashes;
@@ -1979,7 +2049,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
         FillTrackRangeP(*thisTrack[0], rangePs, trk);
 
         if (fmChi2PID.isValid()) {
-           FillTrackChi2PID(fmChi2PID.at(iPart), lar::providerFrom<geo::Geometry>(), trk);
+           FillTrackChi2PID(fmChi2PID.at(iPart), trk);
         }
         if (fmScatterClosestApproach.isValid() && fmScatterClosestApproach.at(iPart).size()==1) {
            FillTrackScatterClosestApproach(fmScatterClosestApproach.at(iPart).front(), trk);
@@ -1994,7 +2064,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
           FillTrackCalo(fmCalo.at(iPart), fmTrackHit.at(iPart),
               (fParams.FillHitsNeutrinoSlices() && NeutrinoSlice) || fParams.FillHitsAllSlices(),
               fParams.TrackHitFillRRStartCut(), fParams.TrackHitFillRREndCut(),
-              lar::providerFrom<geo::Geometry>(), dprop, trk);
+              dprop, trk);
         }
         if (fmCRTHitMatch.isValid() && fDet == kICARUS) {
           art::FindManyP<sbn::crt::CRTHit> CRTT02Hit = FindManyPStrict<sbn::crt::CRTHit>
@@ -2032,7 +2102,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
             FillTrackTruth(fmTrackHit.at(iPart), id_to_hit_energy_map, true_particles, clock_data, trk);
             // Hit truth information corresponding to Calo-Points
             // Assumes truth matching and calo-points are filled
-            if (mc_particles.isValid() && fParams.FillTrackCaloTruth()) FillTrackCaloTruth(id_to_ide_map, *mc_particles, geometry, clock_data, sce, trk);
+            if (mc_particles.isValid() && fParams.FillTrackCaloTruth()) FillTrackCaloTruth(id_to_ide_map, *mc_particles, *geom, wireReadout, clock_data, sce, trk);
           }
         }
       } // thisTrack exists
@@ -2041,7 +2111,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
         assert(thisShower.size() == 1);
 
         SRShower& shw = pfp.shw;
-        FillShowerVars(*thisShower[0], vertex, fmShowerHit.at(iPart), lar::providerFrom<geo::Geometry>(), producer, shw);
+        FillShowerVars(*thisShower[0], vertex, fmShowerHit.at(iPart), wireReadout, producer, shw);
 
         // We may have many residuals per shower depending on how many showers ar in the slice
         if (fmShowerRazzle.isValid() && fmShowerRazzle.at(iPart).size()==1) {
@@ -2189,6 +2259,8 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   rec.hdr.ngenevt = n_gen_evt;
   rec.hdr.mctype  = mctype;
   rec.hdr.sourceName = fSourceFile;
+  rec.hdr.sourceNameHash = fSourceFileHash;
+
   rec.hdr.sourceIndex = fIndexInFile;
   rec.hdr.first_in_file = firstInFile;
   rec.hdr.first_in_subrun = fFirstInSubRun;
@@ -2349,6 +2421,10 @@ void CAFMaker::endJob() {
 
     AddHistogramsToFile(fFile);
     fRecTree->SetDirectory(fFile);
+    if(fGenieTree){
+      fGenieTree->BuildIndex("SourceFileHash", "GENIEEntry");
+      fGenieTree->SetDirectory(fFile);
+    }
     if (fParams.CreateBlindedCAF()) {
       fRecTreeb->SetDirectory(fFileb);
       fRecTreep->SetDirectory(fFilep);
@@ -2370,6 +2446,10 @@ void CAFMaker::endJob() {
 
     AddHistogramsToFile(fFlatFile);
     fFlatTree->SetDirectory(fFlatFile);
+    if(fFlatGenieTree){
+      fFlatGenieTree->BuildIndex("SourceFileHash", "GENIEEntry");
+      fFlatGenieTree->SetDirectory(fFlatFile);
+    }
     if (fParams.CreateBlindedCAF() && fFlatFileb) {
       fFlatTreeb->SetDirectory(fFlatFileb);
       fFlatTreep->SetDirectory(fFlatFilep);
@@ -2384,6 +2464,7 @@ void CAFMaker::endJob() {
       fFlatFilep->cd();
       fFlatFilep->Write();
     }
+
   }
 
   std::map<std::string, std::string> metamap;

sbncode/CAFMaker/CMakeLists.txt

@@ -43,6 +43,8 @@ art_make_library( LIBRARY_NAME sbncode_CAFMaker
                   systematicstools::interface
                   systematicstools::interpreters
                   systematicstools::utility
+                  ${GENIE_LIB_LIST}
+                  nugen::EventGeneratorBase_GENIE
                   )
 
 cet_build_plugin ( CAFMaker art::module