Merge pull request #411 from SBNSoftware/feature/hlay_razzled_pr

Razzled
SBNSoftware · Jan 26, 2024 · 35d7a87 · 35d7a87
2 parents 1be9863 + f5700b7
commit 35d7a87
Show file tree

Hide file tree

Showing 11 changed files with 1,212 additions and 95 deletions.
diff --git a/sbncode/CAFMaker/CAFMakerParams.h b/sbncode/CAFMaker/CAFMakerParams.h
@@ -204,6 +204,12 @@ namespace caf
       "pandoraShowerRazzle"
     };
 
+    Atom<string> PFPRazzledLabel {
+      Name("PFPRazzledLabel"),
+      Comment("Base label of pfp mva particle-id producer."),
+      "pandoraRazzled"
+    };
+
     Atom<string> RecoShowerSelectionLabel {
       Name("RecoShowerSelectionLabel"),
       Comment("Base label of shower selection vars producer."),

diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc
@@ -1686,6 +1686,10 @@ void CAFMaker::produce(art::Event& evt) noexcept {
       FindManyPStrict<sbn::MVAPID>(slcShowers, evt,
           fParams.ShowerRazzleLabel() + slice_tag_suff);
 
+    art::FindManyP<sbn::MVAPID> fmPFPRazzled =
+      FindManyPStrict<sbn::MVAPID>(fmPFPart, evt,
+          fParams.PFPRazzledLabel() + slice_tag_suff);
+
     art::FindOneP<sbn::PFPCNNScore> foCNNScores = 
       FindOnePStrict<sbn::PFPCNNScore>(fmPFPart, evt,
           fParams.CNNScoreLabel() + slice_tag_suff);
@@ -1916,6 +1920,10 @@ void CAFMaker::produce(art::Event& evt) noexcept {
           }
         }
 
+        if (fmPFPRazzled.isValid() && fmPFPRazzled.at(iPart).size()==1) {
+          FillPFPRazzled(fmPFPRazzled.at(iPart).front(), pfp);
+        }
+
         // fill all the stuff
         SRTrack& trk = pfp.trk;
         FillTrackVars(*thisTrack[0], producer, trk);

diff --git a/sbncode/CAFMaker/FillReco.cxx b/sbncode/CAFMaker/FillReco.cxx
@@ -996,6 +996,22 @@ namespace caf
 
   //......................................................................
 
+  void FillPFPRazzled(const art::Ptr<sbn::MVAPID> razzled,
+                      caf::SRPFP& srpfp,
+                      bool allowEmpty)
+  {
+    srpfp.razzled.electronScore = razzled->mvaScoreMap.at(11);
+    srpfp.razzled.muonScore = razzled->mvaScoreMap.at(13);
+    srpfp.razzled.photonScore = razzled->mvaScoreMap.at(22);
+    srpfp.razzled.pionScore = razzled->mvaScoreMap.at(211);
+    srpfp.razzled.protonScore = razzled->mvaScoreMap.at(2212);
+
+    srpfp.razzled.pdg = razzled->BestPDG();
+    srpfp.razzled.bestScore = razzled->BestScore();
+  }
+
+  //......................................................................
+
   void SetNuMuCCPrimary(std::vector<caf::StandardRecord> &recs,
                         std::vector<caf::SRTrueInteraction> &srneutrinos) {
   //   // set is_primary to true by default

diff --git a/sbncode/CAFMaker/FillReco.h b/sbncode/CAFMaker/FillReco.h
@@ -231,6 +231,10 @@ namespace caf
   void FillTPCPMTBarycenterMatch(const sbn::TPCPMTBarycenterMatch *matchInfo,
                            caf::SRSlice& slice);
 
+  void FillPFPRazzled(const art::Ptr<sbn::MVAPID> razzled,
+                      caf::SRPFP& srpfp,
+                      bool allowEmpty = false);
+
   template<class T, class U>
   void CopyPropertyIfSet( const std::map<std::string, T>& props, const std::string& search, U& value );
 }

diff --git a/sbncode/PID/CMakeLists.txt b/sbncode/PID/CMakeLists.txt
@@ -53,6 +53,19 @@ cet_build_plugin( Razzle art::module
   caf_RecoUtils
 )
 
+cet_build_plugin( Razzled art::module
+  LIBRARIES
+  ROOT::TMVA
+  canvas::canvas
+  art_root_io::TFileService_service
+  art::Persistency_Provenance
+  larsim::Utils
+  larsim::MCCheater_BackTrackerService_service
+  larsim::MCCheater_ParticleInventoryService_service
+  lardata::Utilities
+  sbnobj::Common_Reco
+)
+
 install_headers()
 install_fhicl()
 install_source()

diff --git a/sbncode/PID/Razzle_module.cc b/sbncode/PID/Razzle_module.cc
@@ -479,12 +479,12 @@ void Razzle::FillShowerMetrics(const recob::Shower& shower, const std::vector<ar
     showerPlanePitches[plane.ID().Plane] = plane.WirePitch() / cosgamma;
   }
 
-  // Fill only for the best plane, defined as the one with the most hits
-  // Prefer collection plane > 1st induction > 2nd induction
-  bestPlane = shower.best_plane();
-
-  if (bestPlane < 0 || bestPlane > 3)
-    throw cet::exception("Razzle") << "Best plane: " << bestPlane;
+  if(showerPlaneHits[2] >= showerPlaneHits[1] && showerPlaneHits[2] >= showerPlaneHits[0])
+    bestPlane = 2;
+  else if(showerPlaneHits[0] >= showerPlaneHits[1])
+    bestPlane = 0;
+  else
+    bestPlane = 1;
 
   bestdEdx = shower.dEdx()[bestPlane];
   bestdEdx = std::min(bestdEdx, 20.f);