Feature/aantonak caf truthfix

sbncode/CAFMaker/CAFMaker_module.cc

@@ -35,6 +35,7 @@
 #include <string>
 #include <vector>
 #include <array>
+#include <cstdlib>
 
 #ifdef DARWINBUILD
 #include <libgen.h>
@@ -79,6 +80,7 @@
 #include "canvas/Persistency/Common/FindOneP.h"
 #include "canvas/Persistency/Common/Ptr.h"
 #include "canvas/Persistency/Common/PtrVector.h"
+#include "canvas/Utilities/Exception.h"
 
 #include "cetlib_except/exception.h"
 #include "cetlib_except/demangle.h"
@@ -1510,14 +1512,173 @@ void CAFMaker::produce(art::Event& evt) noexcept {
   //#######################################################
   // Fill truths & fake reco
   //#######################################################
-
+  
   caf::SRTruthBranch                  srtruthbranch;
 
   if (mc_particles.isValid()) {
+
+
+    std::cout << std::endl;
+    std::cout << std::endl;
+    std::cout << std::endl;
+    std::cout << std::endl;
+    std::cout << std::endl;
+
+    std::cout << "//------------------------------------------------------------------------------- //" << std::endl;
+    std::cout << "//------------------------------------------------------------------------------- //" << std::endl;
+    std::cout << "//------------------------------------------------------------------------------- //" << std::endl;
+    std::cout << "//------------------------------------------------------------------------------- //" << std::endl;
+
+    std::cout << "DEBUG DEBUG DEBUG" << std::endl;
+    std::cout << "DEBUG DEBUG DEBUG" << std::endl;
+    std::cout << "DEBUG DEBUG DEBUG" << std::endl;
+    std::cout << "DEBUG DEBUG DEBUG" << std::endl;
+    std::cout << "DEBUG DEBUG DEBUG" << std::endl;
+    std::cout << "DEBUG DEBUG DEBUG" << std::endl;
+    std::cout << "DEBUG DEBUG DEBUG" << std::endl;
+    std::cout << "DEBUG DEBUG DEBUG" << std::endl;
+    std::cout << "DEBUG DEBUG DEBUG" << std::endl;
+    std::cout << "DEBUG DEBUG DEBUG" << std::endl;
+
+    std::cout << "//------------------------------------------------------------------------------- //" << std::endl;
+    std::cout << "//------------------------------------------------------------------------------- //" << std::endl;
+    std::cout << "//------------------------------------------------------------------------------- //" << std::endl;
+    std::cout << "//------------------------------------------------------------------------------- //" << std::endl;
+
+
+    std::cout << std::endl;
+    std::cout << std::endl;
+    std::cout << std::endl;
+    std::cout << std::endl;
+    std::cout << std::endl;
+
     art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
     art::ServiceHandle<cheat::BackTrackerService> bt_serv;
 
+    // We need to add a vector of particle mothers IDS that we can check against to make sure we are not double counting! 
+    std::vector<int> mother_ids; // these are the missed mothers
+    int max_g4_track_id = 0;
+    int min_g4_track_id = std::numeric_limits<int>::max();
     for (const simb::MCParticle &part: *mc_particles) {
+      if (part.TrackId() > max_g4_track_id) max_g4_track_id = part.TrackId();
+      if (part.TrackId() < min_g4_track_id) min_g4_track_id = part.TrackId();
+    }
+    std::cout << "Minimum G4 Track ID " << min_g4_track_id << std::endl;
+    for (const simb::MCParticle &part: *mc_particles) {
+
+      std::optional<int> missed_parent_id = std::nullopt;
+      // Now we need to check if the Mother is zero and the parent is not the neutrino/initial state particle. 
+      // If a particle passed to G4 is primary and it's parent is zero, then the parent was not propagated to G4. This is a missed particle of interest that we need to fill in the CAF.
+      //if (part.Mother() == 0 && part.Process() == "primary") {
+      if (part.Mother() == min_g4_track_id-1 && part.Process() == "primary") {
+        std::cout << "Found a primary particle with no mother in FillTrueG4Particle !!!!" << std::endl;
+
+        // Grab the MCTruth associated to this particle
+        const art::Ptr<simb::MCTruth> inventoryTruth = pi_serv->TrackIdToMCTruth_P(part.TrackId());
+        // Loop over the particles in the MCTruth to first find this particle and then check it's Mother again to find it's missed parent
+        if (inventoryTruth) {
+
+          // Loop over the particles in the MCTruth to find this particle and then check it's Mother again to find it's missed parent
+          const simb::MCParticle* matchedGenie = nullptr;
+          double bestScore = std::numeric_limits<double>::infinity();
+
+          int max_genie_track_id = max_g4_track_id;
+          for (int ipart = 0; ipart < inventoryTruth->NParticles(); ++ipart) {
+            const simb::MCParticle& genpart = inventoryTruth->GetParticle(ipart);
+            if (genpart.TrackId() > max_genie_track_id) max_genie_track_id = genpart.TrackId();
+            if (genpart.PdgCode() != part.PdgCode())
+              continue;
+            const auto& gp = part.Momentum(0);
+            const auto& tp = genpart.Momentum(0);
+            const double score =
+              std::pow(gp.Px() - tp.Px(), 2) +
+              std::pow(gp.Py() - tp.Py(), 2) +
+              std::pow(gp.Pz() - tp.Pz(), 2) +
+              std::pow(gp.E()  - tp.E(),  2);
+
+            if (score < bestScore) {
+              bestScore = score;
+              matchedGenie = &genpart;
+            }
+          }
+          if (matchedGenie) {
+            std::cout << "Found the matched GENIE particle in the MCTruth!" << std::endl;
+            std::cout << "Best score: " << bestScore << std::endl;
+
+            // Now we can fill this missed parent in the CAF using our custom FillTrueGENIEParticle function
+            //if (matchedGenie->Mother() != 0) {
+            if (matchedGenie->Mother() != min_g4_track_id-1) {
+              const simb::MCParticle& missedParent = inventoryTruth->GetParticle(matchedGenie->Mother());
+              bool isInitialStateParticle = IsInitialStateParticle(missedParent, *inventoryTruth);
+              if (isInitialStateParticle) {
+                std::cout << "Missed parent is an initial state particle --> reject!" << std::endl;
+              } else {
+                //need_missed_parent = true;
+                const int special_id_offset = max_genie_track_id + max_g4_track_id + 1;
+                const int special_parent_id = special_id_offset + missedParent.TrackId();
+                missed_parent_id = special_parent_id;
+                // need to check if the Mother ID is already in the list of missed mothers to avoid double counting
+                if (std::find(mother_ids.begin(), mother_ids.end(), missedParent.Mother()) != mother_ids.end()) {
+                  std::cout << "This missed parent has already been filled in the CAF, skipping ..." << std::endl;
+                } else {
+                  std::cout << "About to use my custom fill function for the missed parent ..." << std::endl;
+                  int interaction_id = -1;
+                  for (unsigned iTruth = 0; iTruth < mctruths.size(); iTruth++) {
+                    if (inventoryTruth.get() == mctruths[iTruth].get()) { // TODO this gave an error ‘const class simb::MCTruth’ has no member named ‘get
+                      interaction_id = iTruth;
+                      break;
+                    }
+                  }
+                  mother_ids.push_back(missedParent.TrackId());
+                  true_particles.emplace_back();
+                  FillTrueGENIEParticle(missedParent,
+                            fActiveVolumes,
+                            fTPCVolumes,
+                            id_to_ide_map,
+                            id_to_truehit_map,
+                            *bt_serv,
+                            *pi_serv,
+                            mctruths,
+                            true_particles.back(), static_cast<int>(interaction_id), special_id_offset); // TODO --> need to ge the itruth index 
+                  std::cout << "Made it through my fill function for the missed parent" << std::endl;
+
+                  // Now, In principle the parent's parent could also be missed and so forth. Let's while loop until we find a parent that is either an initial state particle or has a mother of zero.
+                  // First check if the next parent is already in the Mother list
+                  if (std::find(mother_ids.begin(), mother_ids.end(), missedParent.Mother()) != mother_ids.end()) {
+                    std::cout << "first grandparent has already been filled in the CAF, skipping ..." << std::endl;
+                  } else {
+                    const simb::MCParticle* currentParent = &missedParent;
+                    //while (currentParent->Mother() != 0 && !IsInitialStateParticle(inventoryTruth->GetParticle(currentParent->Mother()), *inventoryTruth)) {
+                    while (currentParent->Mother() != min_g4_track_id-1 && !IsInitialStateParticle(inventoryTruth->GetParticle(currentParent->Mother()), *inventoryTruth)) {
+                      const simb::MCParticle& nextParent = inventoryTruth->GetParticle(currentParent->Mother());         
+                      mother_ids.push_back(nextParent.TrackId());
+                      true_particles.emplace_back();
+                      FillTrueGENIEParticle(nextParent,
+                                fActiveVolumes,
+                                fTPCVolumes,
+                                id_to_ide_map,
+                                id_to_truehit_map,
+                                *bt_serv,
+                                *pi_serv,
+                                mctruths,
+                                true_particles.back(), static_cast<int>(interaction_id), special_id_offset);
+                      std::cout << "Made it through my fill function for the missed grandparent" << std::endl;
+
+                      if (std::find(mother_ids.begin(), mother_ids.end(), nextParent.Mother()) != mother_ids.end()) {
+                        std::cout << "next grandparent has already been filled in the CAF, exiting while loop ..." << std::endl;
+                        break;
+                      }
+
+                      currentParent = &nextParent;
+                    }
+                  }
+                }
+              }
+            }
+          } // matched Genie Particle
+        } // found inventoryTruth
+      } // primary particle with no motherx
+
       true_particles.emplace_back();
 
       FillTrueG4Particle(part,
@@ -1528,7 +1689,8 @@ void CAFMaker::produce(art::Event& evt) noexcept {
                          *bt_serv,
                          *pi_serv,
                          mctruths,
-                         true_particles.back());
+                         true_particles.back(), missed_parent_id);
+
     }
   }