From 471e77d95d3c052c9cd3d33ce0ac6fcbbdbe5397 Mon Sep 17 00:00:00 2001 From: Alexander Antonakis Date: Wed, 29 Jul 2026 21:14:18 -0500 Subject: [PATCH 1/8] Added truth information for unstable primary particles that are missed by CAFMaker's FillTrueG4Particle function. Debug statements are currently included. --- sbncode/CAFMaker/CAFMaker_module.cc | 111 +++++++++++ sbncode/CAFMaker/FillTrue.cxx | 273 ++++++++++++++++++++++++++++ sbncode/CAFMaker/FillTrue.h | 16 ++ 3 files changed, 400 insertions(+) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index b38ed1f6a..4f6b16039 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -35,6 +35,7 @@ #include #include #include +#include #ifdef DARWINBUILD #include @@ -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" @@ -1404,6 +1406,43 @@ void CAFMaker::produce(art::Event& evt) noexcept { art::fill_ptr_vector(mctruths, mctruth_handle); } + + // Eta Debug Loop + for (size_t itruth = 0; itruth < mctruths.size(); ++itruth) { + + const simb::MCTruth& truth = *mctruths[itruth]; + + std::cout << "\nTruth " << itruth + << " has " << truth.NParticles() + << " particles\n"; + + + for (int i = 0; i < truth.NParticles(); ++i) { + + const simb::MCParticle& p = truth.GetParticle(i); + + if (std::abs(p.PdgCode()) == 221 || std::abs(p.PdgCode()) == 111 || std::abs(p.PdgCode()) == 22) { + + std::cout << "\n====================\n"; + std::cout << "Found GENIE eta or pi0 or gamma\n"; + std::cout << "pdg = " << p.PdgCode() << '\n'; + std::cout << "track = " << p.TrackId() << '\n'; + std::cout << "mother = " << p.Mother() << '\n'; + std::cout << "status = " << p.StatusCode() << '\n'; + std::cout << "process = " << p.Process() << '\n'; + std::cout << "E = " << p.E() << '\n'; + + std::cout << "vertex = " + << p.Vx() << " " + << p.Vy() << " " + << p.Vz() << '\n'; + } + } + } + + // End of Eta Debug + + // And associated GTruth objects art::FindManyP fmp_gtruth = FindManyPStrict(mctruths, evt, fParams.GenLabel()); @@ -1510,6 +1549,21 @@ void CAFMaker::produce(art::Event& evt) noexcept { //####################################################### // Fill truths & fake reco //####################################################### + + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << "Start of True G4 Loop: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; + std::cout << "Start of True G4 Loop: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; + std::cout << "Start of True G4 Loop: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; + std::cout << "Start of True G4 Loop: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; caf::SRTruthBranch srtruthbranch; @@ -1518,6 +1572,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { art::ServiceHandle bt_serv; for (const simb::MCParticle &part: *mc_particles) { + true_particles.emplace_back(); FillTrueG4Particle(part, @@ -1530,6 +1585,62 @@ void CAFMaker::produce(art::Event& evt) noexcept { mctruths, true_particles.back()); } + + // Add loop for unstable Genie particles that don't get propogated to G4 + for (std::size_t itruth = 0; itruth < mctruths.size(); ++itruth) { + const simb::MCTruth& truth = *mctruths.at(itruth); + + if (truth.Origin() != simb::kBeamNeutrino) + continue; + + if (!truth.NeutrinoSet()) + continue; + + for (int ipart = 0; ipart < truth.NParticles(); ++ipart) { + const simb::MCParticle& genpart = truth.GetParticle(ipart); + try { + const art::Ptr inventoryTruth = pi_serv->TrackIdToMCTruth_P(genpart.TrackId()); + if (inventoryTruth) { + // Should be contained in the FillTrueG4 loop + continue; + } + } + catch (const std::exception& e) { + //if (!inventoryTruth) { + // This may be a missed particle of interest! + std::cout << "Exception: " << e.what() << std::endl; + if (genpart.Process() == "primary") { + if (genpart.StatusCode() == 1) continue; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << "Check if it's an initial state particle ..." << std::endl; + bool isInitialStateParticle = IsInitialStateParticle(genpart, truth); + if (isInitialStateParticle) std::cout << "Is an initial state particle --> reject!" << std::endl; + if (isInitialStateParticle) continue; + std::cout << "Found a primary particle that failed track match !!!" << std::endl; + std::cout << "PDG Code: " << genpart.PdgCode() << std::endl; + std::cout << "genpart.Mother(): " << genpart.Mother() << std::endl; + std::cout << "genpart.TrackId(): " << genpart.TrackId() << std::endl; + std::cout << "genpart.StatusCode(): " << genpart.StatusCode() << std::endl; + std::cout << "genpart.Process(): " << genpart.Process() << std::endl; + + std::cout << "About to use my custom fill function ..." << std::endl; + true_particles.emplace_back(); + FillTrueGENIEParticle(genpart, + fActiveVolumes, + fTPCVolumes, + id_to_ide_map, + id_to_truehit_map, + *bt_serv, + *pi_serv, + mctruths, + true_particles.back(), static_cast(itruth)); + std::cout << "Made it through my fill function" << std::endl; + } // primary particles + } // invalid track match + }// loop over particles + } } std::vector> fmpewm; diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 4d1d794d0..b42c69f95 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -803,6 +803,8 @@ namespace caf { srparticle.endE = (exit_point >= 0) ? particle.Momentum(exit_point).E() : -9999.; srparticle.start_process = GetG4ProcessID(particle.Process()); + if (particle.PdgCode() == 221) std::cout << "Eta srparticle start process: " << srparticle.start_process << std::endl; + srparticle.end_process = GetG4ProcessID(particle.EndProcess()); srparticle.G4ID = particle.TrackId(); @@ -824,6 +826,23 @@ namespace caf { srparticle.daughters.push_back(particle.Daughter(i_d)); } + // Eta Debug + if (particle.PdgCode() == 221) { + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << "Found eta in FillTrueG4 Function !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + + } + + // See if this MCParticle matches a genie truth srparticle.interaction_id = -1; @@ -836,6 +855,255 @@ namespace caf { } } //FillTrueG4Particle + + + void FillTrueGENIEParticle(const simb::MCParticle &particle, + const std::vector &active_volumes, + const std::vector> &tpc_volumes, + const std::map>> &id_to_ide_map, + const std::map>> &id_to_truehit_map, + const cheat::BackTrackerService &backtracker, + const cheat::ParticleInventoryService &inventory_service, + const std::vector> &neutrinos, + caf::SRTrueParticle &srparticle, int interaction_id) { + + std::vector> empty; + const std::vector> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; + + std::vector> emptyHits; + const std::vector> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits; + + srparticle.length = 0.; + srparticle.crosses_tpc = false; + srparticle.wallin = caf::kWallNone; + srparticle.wallout = caf::kWallNone; + + for (unsigned c = 0; c < 2; c++) { + SRTrueParticlePlaneInfo init; + init.visE = 0.; + init.nhit = 0; + + for (int p = 0; p < 3; p++) { + srparticle.plane[c][p] = init; + } + } + + for (auto const &ide_pair: particle_ides) { + const geo::WireID &w = ide_pair.first; + const sim::IDE *ide = ide_pair.second; + + if(w.Plane >= 0 && w.Plane < 3 && w.Cryostat < 2){ + srparticle.plane[w.Cryostat][w.Plane].visE += ide->energy / 1000. /* MeV -> GeV*/; + } + } + + for (const art::Ptr h: particle_hits) { + const geo::WireID &w = h->WireID(); + + if(w.Plane >= 0 && w.Plane < 3 && w.Cryostat < 2) { + srparticle.plane[w.Cryostat][w.Plane].nhit ++; + } + } + + // if no trajectory points, then assume outside AV + srparticle.cont_tpc = particle.NumberTrajectoryPoints() > 0; + srparticle.contained = particle.NumberTrajectoryPoints() > 0; + + // Get the entry and exit points + int entry_point = -1; + + int cryostat_index = -1; + int tpc_index = -1; + + for (unsigned j = 0; j < particle.NumberTrajectoryPoints(); j++) { + for (unsigned i = 0; i < active_volumes.size(); i++) { + if (active_volumes.at(i).ContainsPosition(particle.Position(j).Vect())) { + entry_point = j; + cryostat_index = i; + break; + } + } + if (entry_point != -1) break; + } + // get the wall + if (entry_point > 0) { + srparticle.wallin = GetWallCross(active_volumes.at(cryostat_index), particle.Position(entry_point).Vect(), particle.Position(entry_point-1).Vect()); + } + + int exit_point = -1; + + // now setup the cryostat the particle is in + std::vector volumes; + if (entry_point >= 0) { + volumes = tpc_volumes.at(cryostat_index); + for (unsigned i = 0; i < volumes.size(); i++) { + if (volumes[i].ContainsPosition(particle.Position(entry_point).Vect())) { + tpc_index = i; + srparticle.cont_tpc = entry_point == 0; + break; + } + } + srparticle.contained = entry_point == 0; + } + // if we couldn't find the initial point, set not contained + else { + srparticle.contained = false; + } + if (tpc_index < 0) { + srparticle.cont_tpc = false; + } + + // setup aa volumes too for length calc + // Define the volume used for length calculation to be the cryostat volume in question + std::vector aa_volumes; + if (entry_point >= 0) { + const geo::BoxBoundedGeo &v = active_volumes.at(cryostat_index); + aa_volumes.emplace_back(v.MinX(), v.MinY(), v.MinZ(), v.MaxX(), v.MaxY(), v.MaxZ()); + } + + // Get the length and determine if any point leaves the active volume + // + // Use every trajectory point if possible + if (entry_point >= 0) { + // particle trajectory + const simb::MCTrajectory &trajectory = particle.Trajectory(); + TVector3 pos = trajectory.Position(entry_point).Vect(); + for (unsigned i = entry_point+1; i < particle.NumberTrajectoryPoints(); i++) { + TVector3 this_point = trajectory.Position(i).Vect(); + // get the exit point + // update if particle is contained + // check if particle has crossed TPC + if (!srparticle.crosses_tpc) { + for (unsigned j = 0; j < volumes.size(); j++) { + if (volumes[j].ContainsPosition(this_point) && tpc_index >= 0 && j != ((unsigned)tpc_index)) { + srparticle.crosses_tpc = true; + break; + } + } + } + // check if particle has left tpc + if (srparticle.cont_tpc) { + srparticle.cont_tpc = volumes[tpc_index].ContainsPosition(this_point); + } + + if (srparticle.contained) { + srparticle.contained = active_volumes.at(cryostat_index).ContainsPosition(this_point); + } + + // update length + srparticle.length += ContainedLength(this_point, pos, aa_volumes); + + if (!active_volumes.at(cryostat_index).ContainsPosition(this_point) && active_volumes.at(cryostat_index).ContainsPosition(pos)) { + exit_point = i-1; + } + + pos = trajectory.Position(i).Vect(); + } + } + if (exit_point < 0 && entry_point >= 0) { + exit_point = particle.NumberTrajectoryPoints() - 1; + } + if(exit_point >= 0 && entry_point >=0 && exit_point == entry_point && exit_point < static_cast(particle.NumberTrajectoryPoints()) - 1){ + exit_point++; // to avoid exactly the same start and end positions when single index is inside the active volumne + } + if (exit_point >= 0 && ((unsigned)exit_point) < particle.NumberTrajectoryPoints() - 1) { + srparticle.wallout = GetWallCross(active_volumes.at(cryostat_index), particle.Position(exit_point).Vect(), particle.Position(exit_point+1).Vect()); + } + + // other truth information + srparticle.pdg = particle.PdgCode(); + + srparticle.gen = particle.NumberTrajectoryPoints() ? particle.Position().Vect() : TVector3(-9999, -9999, -9999); + srparticle.genT = particle.NumberTrajectoryPoints() ? particle.Position().T() / 1000. /* ns -> us*/: -9999; + srparticle.genp = particle.NumberTrajectoryPoints() ? particle.Momentum().Vect(): TVector3(-9999, -9999, -9999); + srparticle.genE = particle.NumberTrajectoryPoints() ? particle.Momentum().E(): -9999; + + srparticle.start = (entry_point >= 0) ? particle.Position(entry_point).Vect(): TVector3(-9999, -9999, -9999); + srparticle.startT = (entry_point >= 0) ? particle.Position(entry_point).T() / 1000. /* ns-> us*/: -9999; + srparticle.end = (exit_point >= 0) ? particle.Position(exit_point).Vect(): TVector3(-9999, -9999, -9999); + srparticle.endT = (exit_point >= 0) ? particle.Position(exit_point).T() / 1000. /* ns -> us */ : -9999; + + srparticle.startp = (entry_point >= 0) ? particle.Momentum(entry_point).Vect() : TVector3(-9999, -9999, -9999); + srparticle.startE = (entry_point >= 0) ? particle.Momentum(entry_point).E() : -9999.; + srparticle.endp = (exit_point >= 0) ? particle.Momentum(exit_point).Vect() : TVector3(-9999, -9999, -9999); + srparticle.endE = (exit_point >= 0) ? particle.Momentum(exit_point).E() : -9999.; + + srparticle.start_process = GetG4ProcessID(particle.Process()); + if (particle.PdgCode() == 221) std::cout << "Eta srparticle start process: " << srparticle.start_process << std::endl; + + srparticle.end_process = GetG4ProcessID(particle.EndProcess()); + + srparticle.G4ID = particle.TrackId(); + srparticle.parent = particle.Mother(); + + // Set the initial cryostat + srparticle.cryostat = -1; + if (entry_point >= 0) { + for (unsigned c = 0; c < active_volumes.size(); c++) { + if (active_volumes[c].ContainsPosition(particle.Position(entry_point).Vect())) { + srparticle.cryostat = c; + break; + } + } + } + + // Save the daughter particles + for (int i_d = 0; i_d < particle.NumberDaughters(); i_d++) { + srparticle.daughters.push_back(particle.Daughter(i_d)); + } + + // Eta Debug + if (particle.PdgCode() == 221) { + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << "Found eta in FillTrueGENIEParticle Function !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + std::cout << std::endl; + + } + + // See if this MCParticle matches a genie truth + //srparticle.interaction_id = -1; + srparticle.interaction_id = interaction_id; + /* + art::Ptr truth = inventory_service.TrackIdToMCTruth_P(particle.TrackId()); + for (unsigned i = 0; i < neutrinos.size(); i++) { + if (truth.get() == neutrinos[i].get()) { + srparticle.interaction_id = i; + break; + } + } + */ + } //FillTrueGENIEParticle + + + bool IsInitialStateParticle(const simb::MCParticle& particle, + const simb::MCTruth& truth) + { + const int pdg = particle.PdgCode(); + + // Incoming neutrino. + + if (particle.TrackId() == + truth.GetNeutrino().Nu().TrackId()) + return true; + + // Nuclear PDG codes have the form ±10LZZZAAAI. + + if (std::abs(pdg) >= 1000000000) + return true; + + return false; + } + + + void FillFakeReco(const std::vector> &mctruths, const std::vector &srparticles, const std::vector> &mctracks, @@ -1660,3 +1928,8 @@ caf::SRTruthMatch MatchSlice2Truth(const std::vector> &hits } return ret; }//Slc2Truth + + + + + diff --git a/sbncode/CAFMaker/FillTrue.h b/sbncode/CAFMaker/FillTrue.h index 419ed6ee9..4863250a8 100644 --- a/sbncode/CAFMaker/FillTrue.h +++ b/sbncode/CAFMaker/FillTrue.h @@ -80,6 +80,22 @@ namespace caf const std::vector> &neutrinos, caf::SRTrueParticle &srparticle); + + // Added for unstable particles that don't propogate to G4 + void FillTrueGENIEParticle(const simb::MCParticle &particle, + const std::vector &active_volumes, + const std::vector> &tpc_volumes, + const std::map>> &id_to_ide_map, + const std::map>> &id_to_truehit_map, + const cheat::BackTrackerService &backtracker, + const cheat::ParticleInventoryService &inventory_service, + const std::vector> &neutrinos, + caf::SRTrueParticle &srparticle, + int interaction_id); + + bool IsInitialStateParticle(const simb::MCParticle& particle, + const simb::MCTruth& truth); + void FillMeVPrtlTruth(const evgen::ldm::MeVPrtlTruth &truth, const std::vector &active_volumes, caf::SRMeVPrtl &srtruth); From 1b6fee6c992a8e8dec7dc741ea82610f58906119 Mon Sep 17 00:00:00 2001 From: Alexander Antonakis Date: Mon, 3 Aug 2026 12:46:43 -0500 Subject: [PATCH 2/8] restored particle truth heirarchy for missing GENIE primary particles of interest --- sbncode/CAFMaker/CAFMaker_module.cc | 229 +++++++++++++++++++++++++--- sbncode/CAFMaker/FillTrue.cxx | 55 ++----- sbncode/CAFMaker/FillTrue.h | 6 +- 3 files changed, 232 insertions(+), 58 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 4f6b16039..475a854dd 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -1421,7 +1421,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { const simb::MCParticle& p = truth.GetParticle(i); - if (std::abs(p.PdgCode()) == 221 || std::abs(p.PdgCode()) == 111 || std::abs(p.PdgCode()) == 22) { + if (std::abs(p.PdgCode()) == 221 || std::abs(p.PdgCode()) == 111 || std::abs(p.PdgCode()) == 22 || p.TrackId() == 6) { std::cout << "\n====================\n"; std::cout << "Found GENIE eta or pi0 or gamma\n"; @@ -1550,20 +1550,12 @@ void CAFMaker::produce(art::Event& evt) noexcept { // Fill truths & fake reco //####################################################### + std::cout << std::endl; std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << "Start of True G4 Loop: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; - std::cout << "Start of True G4 Loop: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; std::cout << "Start of True G4 Loop: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; - std::cout << "Start of True G4 Loop: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; std::cout << std::endl; + caf::SRTruthBranch srtruthbranch; @@ -1571,8 +1563,18 @@ void CAFMaker::produce(art::Event& evt) noexcept { art::ServiceHandle pi_serv; art::ServiceHandle 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 mother_ids; // these are the missed mothers + int max_g4_track_id = 0; for (const simb::MCParticle &part: *mc_particles) { - + if (part.TrackId() > max_g4_track_id) max_g4_track_id = part.TrackId(); + } + std::cout << "Start of particle loop for FillTrueG4Particle ..." << std::endl; + int particle_counter = 0; + for (const simb::MCParticle &part: *mc_particles) { + std::cout << "Check particle " << particle_counter << " pdg " << part.PdgCode() << " trkID " << part.TrackId() << std::endl; + //if (part.PdgCode() == 212212) std::cout << "Weird parent is filled in FillTrueG4Particle !!!!" << std::endl; + /* true_particles.emplace_back(); FillTrueG4Particle(part, @@ -1584,41 +1586,231 @@ void CAFMaker::produce(art::Event& evt) noexcept { *pi_serv, mctruths, true_particles.back()); + */ + + //bool needs_missed_parent = false; + std::optional 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") { + std::cout << "Found a primary particle with no mother in FillTrueG4Particle !!!!" << std::endl; + std::cout << "PDG Code: " << part.PdgCode() << std::endl; + std::cout << "Track ID: " << part.TrackId() << std::endl; + std::cout << "Status Code: " << part.StatusCode() << std::endl; + std::cout << "Process: " << part.Process() << std::endl; + + + // Grab the MCTruth associated to this particle + const art::Ptr 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::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; + std::cout << "PDG Code: " << matchedGenie->PdgCode() << std::endl; + std::cout << "Track ID: " << matchedGenie->TrackId() << std::endl; + std::cout << "Mother: " << matchedGenie->Mother() << std::endl; + std::cout << "Status Code: " << matchedGenie->StatusCode() << std::endl; + std::cout << "Process: " << matchedGenie->Process() << std::endl; + + if (matchedGenie->Mother() != 0) { + std::cout << "Found a missed parent in the MCTruth for this particle !!!!" << std::endl; + std::cout << "Missed Parent PDG Code: " << inventoryTruth->GetParticle(matchedGenie->Mother()).PdgCode() << std::endl; + std::cout << "Missed Parent Track ID: " << inventoryTruth->GetParticle(matchedGenie->Mother()).TrackId() << std::endl; + std::cout << "Missed Parent Status Code: " << inventoryTruth->GetParticle(matchedGenie->Mother()).StatusCode() << std::endl; + std::cout << "Missed Parent Process: " << inventoryTruth->GetParticle(matchedGenie->Mother()).Process() << std::endl; + } + + // Now we can fill this missed parent in the CAF using our custom FillTrueGENIEParticle function + if (matchedGenie->Mother() != 0) { + 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(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)) { + const simb::MCParticle& nextParent = inventoryTruth->GetParticle(currentParent->Mother()); + std::cout << "Found a missed grandparent in the MCTruth for this particle !!!!" << std::endl; + std::cout << "Missed Grandparent PDG Code: " << nextParent.PdgCode() << std::endl; + std::cout << "Missed Grandparent Track ID: " << nextParent.TrackId() << std::endl; + std::cout << "Missed Grandparent Status Code: " << nextParent.StatusCode() << std::endl; + std::cout << "Missed Grandparent Process: " << nextParent.Process() << std::endl; + + 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(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 + std::cout << std::endl; + std::cout << std::endl; + } // primary particle with no motherx + std::cout << "About to Fill particle " << particle_counter << " in FillTrueG4Particle ..." << std::endl; + true_particles.emplace_back(); + + FillTrueG4Particle(part, + fActiveVolumes, + fTPCVolumes, + id_to_ide_map, + id_to_truehit_map, + *bt_serv, + *pi_serv, + mctruths, + true_particles.back(), missed_parent_id); + particle_counter++; } + std::cout << "End of loop for FillTrueG4Particle /////////////////////////////////////////" << std::endl; + std::cout << std::endl; + std::cout << std::endl; + + /* + std::cout << "Start of my loop over MC truths to find missed particles of interest ..." << std::endl; + std::cout << std::endl; // Add loop for unstable Genie particles that don't get propogated to G4 for (std::size_t itruth = 0; itruth < mctruths.size(); ++itruth) { - const simb::MCTruth& truth = *mctruths.at(itruth); + std::cout << "Analyzing mctruth " << itruth << std::endl; + const simb::MCTruth& truth = *mctruths.at(itruth); + if (truth.Origin() != simb::kBeamNeutrino) continue; if (!truth.NeutrinoSet()) continue; + std::cout << "mctruth is a beam neutrino interaction" << std::endl; + std::cout << std::endl; + std::cout << "Starting loop over true particles ..." << std::endl; for (int ipart = 0; ipart < truth.NParticles(); ++ipart) { const simb::MCParticle& genpart = truth.GetParticle(ipart); + std::cout << "Analyzing particle " << ipart << std::endl; + if (genpart.PdgCode() == 212212) std::cout << "Found missed parent! for event 138" << std::endl; + std::cout << "Check if particle has it's track ID matched by pi_serv->TrackIdToMCTruth_P" << std::endl; try { const art::Ptr inventoryTruth = pi_serv->TrackIdToMCTruth_P(genpart.TrackId()); if (inventoryTruth) { // Should be contained in the FillTrueG4 loop + // However, some seem to get missed somehow? + std::cout << "ipart " << ipart << " passed pi_serv->TrackIdToMCTruth_P" < reject!" << std::endl; if (isInitialStateParticle) continue; - std::cout << "Found a primary particle that failed track match !!!" << std::endl; + std::cout << "Found a primary particle candidate !!!" << std::endl; std::cout << "PDG Code: " << genpart.PdgCode() << std::endl; std::cout << "genpart.Mother(): " << genpart.Mother() << std::endl; std::cout << "genpart.TrackId(): " << genpart.TrackId() << std::endl; @@ -1639,8 +1831,11 @@ void CAFMaker::produce(art::Event& evt) noexcept { std::cout << "Made it through my fill function" << std::endl; } // primary particles } // invalid track match + std::cout << std::endl; + std::cout << std::endl; }// loop over particles } + */ } std::vector> fmpewm; diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index b42c69f95..6635cf943 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -9,6 +9,7 @@ #include "CLHEP/Random/RandGauss.h" +#include #include #include @@ -477,10 +478,12 @@ namespace caf { } for(const caf::SRTrueParticle& part: srparticles){ + if (part.pdg == 212212) std::cout << "Weird parent in FillTrueNeutrino" << std::endl; + // save the G4 particles that came from this interaction if(part.interaction_id == (int)i) { if(part.start_process == caf::kG4primary) srneutrino.prim.push_back(part); - + //std::cout << "succeeded FillNeutrino check" << std::endl; // total up the deposited energy for(int p = 0; p < 3; ++p) { for (int i_cryo = 0; i_cryo < 2; i_cryo++) { @@ -639,7 +642,7 @@ namespace caf { const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, const std::vector> &neutrinos, - caf::SRTrueParticle &srparticle) { + caf::SRTrueParticle &srparticle, std::optional new_mother) { std::vector> empty; const std::vector> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; @@ -808,7 +811,13 @@ namespace caf { srparticle.end_process = GetG4ProcessID(particle.EndProcess()); srparticle.G4ID = particle.TrackId(); - srparticle.parent = particle.Mother(); + srparticle.parent = particle.Mother() > 0 ? static_cast(particle.Mother()) : 0u; + if (new_mother.has_value()) { + const int requested_parent = new_mother.value(); + std::cout << "Changing parent of particle " << srparticle.G4ID << " from " << srparticle.parent << " to " << requested_parent << std::endl; + srparticle.parent = static_cast(requested_parent); + std::cout << "srparticle.parent is now " << srparticle.parent << std::endl; + } // Set the initial cryostat srparticle.cryostat = -1; @@ -826,23 +835,6 @@ namespace caf { srparticle.daughters.push_back(particle.Daughter(i_d)); } - // Eta Debug - if (particle.PdgCode() == 221) { - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << "Found eta in FillTrueG4 Function !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - - } - - // See if this MCParticle matches a genie truth srparticle.interaction_id = -1; @@ -865,7 +857,7 @@ namespace caf { const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, const std::vector> &neutrinos, - caf::SRTrueParticle &srparticle, int interaction_id) { + caf::SRTrueParticle &srparticle, int interaction_id, int id_offset) { std::vector> empty; const std::vector> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty; @@ -1033,8 +1025,9 @@ namespace caf { srparticle.end_process = GetG4ProcessID(particle.EndProcess()); - srparticle.G4ID = particle.TrackId(); - srparticle.parent = particle.Mother(); + // Special GENIE particles get a shifted positive ID so they do not overlap with G4 track IDs. + srparticle.G4ID = particle.TrackId() + id_offset; + srparticle.parent = particle.Mother() > 0 ? static_cast(particle.Mother() + id_offset) : 0u; // Set the initial cryostat srparticle.cryostat = -1; @@ -1052,22 +1045,6 @@ namespace caf { srparticle.daughters.push_back(particle.Daughter(i_d)); } - // Eta Debug - if (particle.PdgCode() == 221) { - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << "Found eta in FillTrueGENIEParticle Function !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - std::cout << std::endl; - - } - // See if this MCParticle matches a genie truth //srparticle.interaction_id = -1; srparticle.interaction_id = interaction_id; diff --git a/sbncode/CAFMaker/FillTrue.h b/sbncode/CAFMaker/FillTrue.h index 4863250a8..084152739 100644 --- a/sbncode/CAFMaker/FillTrue.h +++ b/sbncode/CAFMaker/FillTrue.h @@ -1,6 +1,7 @@ #ifndef CAF_FILLTRUE_H #define CAF_FILLTRUE_H +#include #include "TRandom.h" #include "TDatabasePDG.h" #include "CLHEP/Random/RandEngine.h" // CLHEP::HepRandomEngine @@ -78,7 +79,7 @@ namespace caf const cheat::BackTrackerService &backtracker, const cheat::ParticleInventoryService &inventory_service, const std::vector> &neutrinos, - caf::SRTrueParticle &srparticle); + caf::SRTrueParticle &srparticle, std::optional new_mother); // Added for unstable particles that don't propogate to G4 @@ -91,7 +92,8 @@ namespace caf const cheat::ParticleInventoryService &inventory_service, const std::vector> &neutrinos, caf::SRTrueParticle &srparticle, - int interaction_id); + int interaction_id, + int id_offset); bool IsInitialStateParticle(const simb::MCParticle& particle, const simb::MCTruth& truth); From c7ca06d900fbdda8ec75cbd5a3daa83c852a9dcf Mon Sep 17 00:00:00 2001 From: Alexander Antonakis Date: Fri, 21 Aug 2026 14:37:35 -0500 Subject: [PATCH 3/8] revised debug statements --- sbncode/CAFMaker/CAFMaker_module.cc | 233 +++++----------------------- sbncode/CAFMaker/FillTrue.cxx | 13 +- 2 files changed, 44 insertions(+), 202 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 475a854dd..10b83bdbd 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -1406,43 +1406,6 @@ void CAFMaker::produce(art::Event& evt) noexcept { art::fill_ptr_vector(mctruths, mctruth_handle); } - - // Eta Debug Loop - for (size_t itruth = 0; itruth < mctruths.size(); ++itruth) { - - const simb::MCTruth& truth = *mctruths[itruth]; - - std::cout << "\nTruth " << itruth - << " has " << truth.NParticles() - << " particles\n"; - - - for (int i = 0; i < truth.NParticles(); ++i) { - - const simb::MCParticle& p = truth.GetParticle(i); - - if (std::abs(p.PdgCode()) == 221 || std::abs(p.PdgCode()) == 111 || std::abs(p.PdgCode()) == 22 || p.TrackId() == 6) { - - std::cout << "\n====================\n"; - std::cout << "Found GENIE eta or pi0 or gamma\n"; - std::cout << "pdg = " << p.PdgCode() << '\n'; - std::cout << "track = " << p.TrackId() << '\n'; - std::cout << "mother = " << p.Mother() << '\n'; - std::cout << "status = " << p.StatusCode() << '\n'; - std::cout << "process = " << p.Process() << '\n'; - std::cout << "E = " << p.E() << '\n'; - - std::cout << "vertex = " - << p.Vx() << " " - << p.Vy() << " " - << p.Vz() << '\n'; - } - } - } - - // End of Eta Debug - - // And associated GTruth objects art::FindManyP fmp_gtruth = FindManyPStrict(mctruths, evt, fParams.GenLabel()); @@ -1550,16 +1513,45 @@ void CAFMaker::produce(art::Event& evt) noexcept { // Fill truths & fake reco //####################################################### - - std::cout << std::endl; - std::cout << std::endl; - std::cout << "Start of True G4 Loop: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; - std::cout << std::endl; - - 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 pi_serv; art::ServiceHandle bt_serv; @@ -1569,37 +1561,15 @@ void CAFMaker::produce(art::Event& evt) noexcept { for (const simb::MCParticle &part: *mc_particles) { if (part.TrackId() > max_g4_track_id) max_g4_track_id = part.TrackId(); } - std::cout << "Start of particle loop for FillTrueG4Particle ..." << std::endl; - int particle_counter = 0; + for (const simb::MCParticle &part: *mc_particles) { - std::cout << "Check particle " << particle_counter << " pdg " << part.PdgCode() << " trkID " << part.TrackId() << std::endl; - //if (part.PdgCode() == 212212) std::cout << "Weird parent is filled in FillTrueG4Particle !!!!" << std::endl; - /* - true_particles.emplace_back(); - - FillTrueG4Particle(part, - fActiveVolumes, - fTPCVolumes, - id_to_ide_map, - id_to_truehit_map, - *bt_serv, - *pi_serv, - mctruths, - true_particles.back()); - */ - //bool needs_missed_parent = false; std::optional 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") { std::cout << "Found a primary particle with no mother in FillTrueG4Particle !!!!" << std::endl; - std::cout << "PDG Code: " << part.PdgCode() << std::endl; - std::cout << "Track ID: " << part.TrackId() << std::endl; - std::cout << "Status Code: " << part.StatusCode() << std::endl; - std::cout << "Process: " << part.Process() << std::endl; - - + // Grab the MCTruth associated to this particle const art::Ptr 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 @@ -1631,20 +1601,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { if (matchedGenie) { std::cout << "Found the matched GENIE particle in the MCTruth!" << std::endl; std::cout << "Best score: " << bestScore << std::endl; - std::cout << "PDG Code: " << matchedGenie->PdgCode() << std::endl; - std::cout << "Track ID: " << matchedGenie->TrackId() << std::endl; - std::cout << "Mother: " << matchedGenie->Mother() << std::endl; - std::cout << "Status Code: " << matchedGenie->StatusCode() << std::endl; - std::cout << "Process: " << matchedGenie->Process() << std::endl; - - if (matchedGenie->Mother() != 0) { - std::cout << "Found a missed parent in the MCTruth for this particle !!!!" << std::endl; - std::cout << "Missed Parent PDG Code: " << inventoryTruth->GetParticle(matchedGenie->Mother()).PdgCode() << std::endl; - std::cout << "Missed Parent Track ID: " << inventoryTruth->GetParticle(matchedGenie->Mother()).TrackId() << std::endl; - std::cout << "Missed Parent Status Code: " << inventoryTruth->GetParticle(matchedGenie->Mother()).StatusCode() << std::endl; - std::cout << "Missed Parent Process: " << inventoryTruth->GetParticle(matchedGenie->Mother()).Process() << std::endl; - } - + // Now we can fill this missed parent in the CAF using our custom FillTrueGENIEParticle function if (matchedGenie->Mother() != 0) { const simb::MCParticle& missedParent = inventoryTruth->GetParticle(matchedGenie->Mother()); @@ -1688,13 +1645,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { } else { const simb::MCParticle* currentParent = &missedParent; while (currentParent->Mother() != 0 && !IsInitialStateParticle(inventoryTruth->GetParticle(currentParent->Mother()), *inventoryTruth)) { - const simb::MCParticle& nextParent = inventoryTruth->GetParticle(currentParent->Mother()); - std::cout << "Found a missed grandparent in the MCTruth for this particle !!!!" << std::endl; - std::cout << "Missed Grandparent PDG Code: " << nextParent.PdgCode() << std::endl; - std::cout << "Missed Grandparent Track ID: " << nextParent.TrackId() << std::endl; - std::cout << "Missed Grandparent Status Code: " << nextParent.StatusCode() << std::endl; - std::cout << "Missed Grandparent Process: " << nextParent.Process() << std::endl; - + const simb::MCParticle& nextParent = inventoryTruth->GetParticle(currentParent->Mother()); mother_ids.push_back(nextParent.TrackId()); true_particles.emplace_back(); FillTrueGENIEParticle(nextParent, @@ -1721,10 +1672,8 @@ void CAFMaker::produce(art::Event& evt) noexcept { } } // matched Genie Particle } // found inventoryTruth - std::cout << std::endl; - std::cout << std::endl; } // primary particle with no motherx - std::cout << "About to Fill particle " << particle_counter << " in FillTrueG4Particle ..." << std::endl; + true_particles.emplace_back(); FillTrueG4Particle(part, @@ -1736,106 +1685,8 @@ void CAFMaker::produce(art::Event& evt) noexcept { *pi_serv, mctruths, true_particles.back(), missed_parent_id); - particle_counter++; + } - std::cout << "End of loop for FillTrueG4Particle /////////////////////////////////////////" << std::endl; - std::cout << std::endl; - std::cout << std::endl; - - /* - std::cout << "Start of my loop over MC truths to find missed particles of interest ..." << std::endl; - std::cout << std::endl; - - // Add loop for unstable Genie particles that don't get propogated to G4 - for (std::size_t itruth = 0; itruth < mctruths.size(); ++itruth) { - - std::cout << "Analyzing mctruth " << itruth << std::endl; - const simb::MCTruth& truth = *mctruths.at(itruth); - - if (truth.Origin() != simb::kBeamNeutrino) - continue; - - if (!truth.NeutrinoSet()) - continue; - std::cout << "mctruth is a beam neutrino interaction" << std::endl; - - std::cout << std::endl; - std::cout << "Starting loop over true particles ..." << std::endl; - for (int ipart = 0; ipart < truth.NParticles(); ++ipart) { - const simb::MCParticle& genpart = truth.GetParticle(ipart); - std::cout << "Analyzing particle " << ipart << std::endl; - if (genpart.PdgCode() == 212212) std::cout << "Found missed parent! for event 138" << std::endl; - std::cout << "Check if particle has it's track ID matched by pi_serv->TrackIdToMCTruth_P" << std::endl; - try { - const art::Ptr inventoryTruth = pi_serv->TrackIdToMCTruth_P(genpart.TrackId()); - if (inventoryTruth) { - // Should be contained in the FillTrueG4 loop - // However, some seem to get missed somehow? - std::cout << "ipart " << ipart << " passed pi_serv->TrackIdToMCTruth_P" < reject!" << std::endl; - if (isInitialStateParticle) continue; - std::cout << "Found a primary particle candidate !!!" << std::endl; - std::cout << "PDG Code: " << genpart.PdgCode() << std::endl; - std::cout << "genpart.Mother(): " << genpart.Mother() << std::endl; - std::cout << "genpart.TrackId(): " << genpart.TrackId() << std::endl; - std::cout << "genpart.StatusCode(): " << genpart.StatusCode() << std::endl; - std::cout << "genpart.Process(): " << genpart.Process() << std::endl; - - std::cout << "About to use my custom fill function ..." << std::endl; - true_particles.emplace_back(); - FillTrueGENIEParticle(genpart, - fActiveVolumes, - fTPCVolumes, - id_to_ide_map, - id_to_truehit_map, - *bt_serv, - *pi_serv, - mctruths, - true_particles.back(), static_cast(itruth)); - std::cout << "Made it through my fill function" << std::endl; - } // primary particles - } // invalid track match - std::cout << std::endl; - std::cout << std::endl; - }// loop over particles - } - */ } std::vector> fmpewm; diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 6635cf943..cb4fddbfd 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -1045,18 +1045,9 @@ namespace caf { srparticle.daughters.push_back(particle.Daughter(i_d)); } - // See if this MCParticle matches a genie truth - //srparticle.interaction_id = -1; + // Set interaction_id to the matched MCTruth srparticle.interaction_id = interaction_id; - /* - art::Ptr truth = inventory_service.TrackIdToMCTruth_P(particle.TrackId()); - for (unsigned i = 0; i < neutrinos.size(); i++) { - if (truth.get() == neutrinos[i].get()) { - srparticle.interaction_id = i; - break; - } - } - */ + } //FillTrueGENIEParticle From c883a118d78cd1ad8ae40abdb1aa1574708935af Mon Sep 17 00:00:00 2001 From: Alexander Antonakis Date: Wed, 9 Sep 2026 22:18:09 -0500 Subject: [PATCH 4/8] fixed the parent track ID identifier associated with primary particles wtih missing parents --- sbncode/CAFMaker/CAFMaker_module.cc | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index 10b83bdbd..be18ac04a 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -1558,16 +1558,19 @@ void CAFMaker::produce(art::Event& evt) noexcept { // 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 mother_ids; // these are the missed mothers int max_g4_track_id = 0; + int min_g4_track_id = std::numeric_limits::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 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() == 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 @@ -1603,7 +1606,8 @@ void CAFMaker::produce(art::Event& evt) noexcept { 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() != 0) { + if (matchedGenie->Mother() != min_g4_track_id-1) { const simb::MCParticle& missedParent = inventoryTruth->GetParticle(matchedGenie->Mother()); bool isInitialStateParticle = IsInitialStateParticle(missedParent, *inventoryTruth); if (isInitialStateParticle) { @@ -1644,7 +1648,8 @@ void CAFMaker::produce(art::Event& evt) noexcept { 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() != 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(); From 8f858b49ab1c82a92a04593dfcdbf2b8cf65920c Mon Sep 17 00:00:00 2001 From: Alexander Antonakis Date: Wed, 30 Sep 2026 14:20:03 -0500 Subject: [PATCH 5/8] Fixed failure rate due to a logical error with the newly assigned missed parent IDs. The IDs were not guaranteed to be uniquie when multiple MCTruths are present. --- sbncode/CAFMaker/CAFMaker_module.cc | 67 +++++++++++++++++------------ sbncode/CAFMaker/FillTrue.cxx | 7 ++- 2 files changed, 46 insertions(+), 28 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index be18ac04a..c4e4f3709 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -1559,11 +1559,29 @@ void CAFMaker::produce(art::Event& evt) noexcept { std::vector mother_ids; // these are the missed mothers int max_g4_track_id = 0; int min_g4_track_id = std::numeric_limits::max(); + std::vector genie_track_id_offsets(mctruths.size(), 0); + int cumulative_genie_track_id_max = 0; 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(); } + for (unsigned iTruth = 0; iTruth < mctruths.size(); ++iTruth) { + const art::Ptr& truth = mctruths[iTruth]; + if (!truth) continue; + int max_truth_genie_track_id = 0; + for (int ipart = 0; ipart < truth->NParticles(); ++ipart) { + const simb::MCParticle& genpart = truth->GetParticle(ipart); + if (genpart.TrackId() > max_truth_genie_track_id) max_truth_genie_track_id = genpart.TrackId(); + } + genie_track_id_offsets[iTruth] = cumulative_genie_track_id_max; + cumulative_genie_track_id_max += max_truth_genie_track_id + 1; + } std::cout << "Minimum G4 Track ID " << min_g4_track_id << std::endl; + std::cout << "GENIE track ID offsets:"; + for (unsigned iTruth = 0; iTruth < genie_track_id_offsets.size(); ++iTruth) { + std::cout << " [" << iTruth << "]: " << genie_track_id_offsets[iTruth]; + } + std::cout << std::endl; for (const simb::MCParticle &part: *mc_particles) { std::optional missed_parent_id = std::nullopt; @@ -1582,10 +1600,8 @@ void CAFMaker::produce(art::Event& evt) noexcept { const simb::MCParticle* matchedGenie = nullptr; double bestScore = std::numeric_limits::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); @@ -1613,23 +1629,20 @@ void CAFMaker::produce(art::Event& evt) noexcept { 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; + int interaction_id = -1; + for (unsigned iTruth = 0; iTruth < mctruths.size(); iTruth++) { + if (inventoryTruth.get() == mctruths[iTruth].get()) { + interaction_id = iTruth; + break; + } + } + const int special_id_offset = genie_track_id_offsets[interaction_id] + max_g4_track_id + 1; const int special_parent_id = special_id_offset + missedParent.TrackId(); missed_parent_id = special_parent_id; + const int synthetic_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()); + if (std::find(mother_ids.begin(), mother_ids.end(), synthetic_missed_parent_id) == mother_ids.end()) { + mother_ids.push_back(synthetic_missed_parent_id); true_particles.emplace_back(); FillTrueGENIEParticle(missedParent, fActiveVolumes, @@ -1639,19 +1652,19 @@ void CAFMaker::produce(art::Event& evt) noexcept { *bt_serv, *pi_serv, mctruths, - true_particles.back(), static_cast(interaction_id), special_id_offset); // TODO --> need to ge the itruth index + true_particles.back(), static_cast(interaction_id), special_id_offset); 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 int synthetic_missed_grandparent_id = special_id_offset + missedParent.Mother(); + if (std::find(mother_ids.begin(), mother_ids.end(), synthetic_missed_grandparent_id) == mother_ids.end()) { 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()); + const simb::MCParticle& nextParent = inventoryTruth->GetParticle(currentParent->Mother()); + const int synthetic_next_parent_id = special_id_offset + nextParent.TrackId(); + mother_ids.push_back(synthetic_next_parent_id); true_particles.emplace_back(); FillTrueGENIEParticle(nextParent, fActiveVolumes, @@ -1662,13 +1675,11 @@ void CAFMaker::produce(art::Event& evt) noexcept { *pi_serv, mctruths, true_particles.back(), static_cast(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; + + const int synthetic_next_grandparent_id = special_id_offset + nextParent.Mother(); + if (std::find(mother_ids.begin(), mother_ids.end(), synthetic_next_grandparent_id) != mother_ids.end()) { break; } - currentParent = &nextParent; } } @@ -1677,6 +1688,8 @@ void CAFMaker::produce(art::Event& evt) noexcept { } } // matched Genie Particle } // found inventoryTruth + std::cout << std::endl; + std::cout << std::endl; } // primary particle with no motherx true_particles.emplace_back(); diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index cb4fddbfd..7ea8a955b 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -1020,11 +1020,14 @@ namespace caf { srparticle.endp = (exit_point >= 0) ? particle.Momentum(exit_point).Vect() : TVector3(-9999, -9999, -9999); srparticle.endE = (exit_point >= 0) ? particle.Momentum(exit_point).E() : -9999.; + std::cout << "Filling GENIE particle start_process for pdg " << particle.PdgCode() << std::endl; srparticle.start_process = GetG4ProcessID(particle.Process()); if (particle.PdgCode() == 221) std::cout << "Eta srparticle start process: " << srparticle.start_process << std::endl; - + std::cout << "Filling GENIE particle end_process for pdg " << particle.PdgCode() << std::endl; + srparticle.end_process = GetG4ProcessID(particle.EndProcess()); + std::cout << "Filling GENIE particle G4ID for pdg " << particle.PdgCode() << std::endl; // Special GENIE particles get a shifted positive ID so they do not overlap with G4 track IDs. srparticle.G4ID = particle.TrackId() + id_offset; srparticle.parent = particle.Mother() > 0 ? static_cast(particle.Mother() + id_offset) : 0u; @@ -1048,6 +1051,8 @@ namespace caf { // Set interaction_id to the matched MCTruth srparticle.interaction_id = interaction_id; + std::cout << std::endl; + } //FillTrueGENIEParticle From a14db9dc2a46f974aade6de247b4f0dc918ba134 Mon Sep 17 00:00:00 2001 From: Alexander Antonakis Date: Thu, 1 Oct 2026 02:12:39 -0500 Subject: [PATCH 6/8] Cleaned up debug statements --- sbncode/CAFMaker/CAFMaker_module.cc | 70 +++++------------------------ sbncode/CAFMaker/FillTrue.cxx | 12 ----- 2 files changed, 10 insertions(+), 72 deletions(-) diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index c4e4f3709..06ab89f3e 100644 --- a/sbncode/CAFMaker/CAFMaker_module.cc +++ b/sbncode/CAFMaker/CAFMaker_module.cc @@ -1517,41 +1517,6 @@ void CAFMaker::produce(art::Event& evt) noexcept { 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 pi_serv; art::ServiceHandle bt_serv; @@ -1576,21 +1541,15 @@ void CAFMaker::produce(art::Event& evt) noexcept { genie_track_id_offsets[iTruth] = cumulative_genie_track_id_max; cumulative_genie_track_id_max += max_truth_genie_track_id + 1; } - std::cout << "Minimum G4 Track ID " << min_g4_track_id << std::endl; - std::cout << "GENIE track ID offsets:"; - for (unsigned iTruth = 0; iTruth < genie_track_id_offsets.size(); ++iTruth) { - std::cout << " [" << iTruth << "]: " << genie_track_id_offsets[iTruth]; - } - std::cout << std::endl; + for (const simb::MCParticle &part: *mc_particles) { std::optional 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") { + // Now we need to check if the Mother is missing (== min_g4_track_id-1) and the parent is not the neutrino/initial state particle. + // If a particle passed to G4 is primary and it's parent is the smallest G4 track ID -1, 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() == 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 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 @@ -1618,17 +1577,12 @@ void CAFMaker::produce(art::Event& evt) noexcept { } } 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) { + // Now we can check if we need to fill this missed parent in the CAF 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 { + if (!isInitialStateParticle) { int interaction_id = -1; for (unsigned iTruth = 0; iTruth < mctruths.size(); iTruth++) { if (inventoryTruth.get() == mctruths[iTruth].get()) { @@ -1653,14 +1607,12 @@ void CAFMaker::produce(art::Event& evt) noexcept { *pi_serv, mctruths, true_particles.back(), static_cast(interaction_id), special_id_offset); - 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. + + // Now, In principle the parent's parent could also be missed and so forth. // First check if the next parent is already in the Mother list const int synthetic_missed_grandparent_id = special_id_offset + missedParent.Mother(); if (std::find(mother_ids.begin(), mother_ids.end(), synthetic_missed_grandparent_id) == mother_ids.end()) { 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()); const int synthetic_next_parent_id = special_id_offset + nextParent.TrackId(); @@ -1688,9 +1640,7 @@ void CAFMaker::produce(art::Event& evt) noexcept { } } // matched Genie Particle } // found inventoryTruth - std::cout << std::endl; - std::cout << std::endl; - } // primary particle with no motherx + } // primary particle with no mother true_particles.emplace_back(); diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 7ea8a955b..01793795d 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -806,17 +806,13 @@ namespace caf { srparticle.endE = (exit_point >= 0) ? particle.Momentum(exit_point).E() : -9999.; srparticle.start_process = GetG4ProcessID(particle.Process()); - if (particle.PdgCode() == 221) std::cout << "Eta srparticle start process: " << srparticle.start_process << std::endl; - srparticle.end_process = GetG4ProcessID(particle.EndProcess()); srparticle.G4ID = particle.TrackId(); srparticle.parent = particle.Mother() > 0 ? static_cast(particle.Mother()) : 0u; if (new_mother.has_value()) { const int requested_parent = new_mother.value(); - std::cout << "Changing parent of particle " << srparticle.G4ID << " from " << srparticle.parent << " to " << requested_parent << std::endl; srparticle.parent = static_cast(requested_parent); - std::cout << "srparticle.parent is now " << srparticle.parent << std::endl; } // Set the initial cryostat @@ -1020,14 +1016,9 @@ namespace caf { srparticle.endp = (exit_point >= 0) ? particle.Momentum(exit_point).Vect() : TVector3(-9999, -9999, -9999); srparticle.endE = (exit_point >= 0) ? particle.Momentum(exit_point).E() : -9999.; - std::cout << "Filling GENIE particle start_process for pdg " << particle.PdgCode() << std::endl; srparticle.start_process = GetG4ProcessID(particle.Process()); - if (particle.PdgCode() == 221) std::cout << "Eta srparticle start process: " << srparticle.start_process << std::endl; - std::cout << "Filling GENIE particle end_process for pdg " << particle.PdgCode() << std::endl; - srparticle.end_process = GetG4ProcessID(particle.EndProcess()); - std::cout << "Filling GENIE particle G4ID for pdg " << particle.PdgCode() << std::endl; // Special GENIE particles get a shifted positive ID so they do not overlap with G4 track IDs. srparticle.G4ID = particle.TrackId() + id_offset; srparticle.parent = particle.Mother() > 0 ? static_cast(particle.Mother() + id_offset) : 0u; @@ -1050,9 +1041,6 @@ namespace caf { // Set interaction_id to the matched MCTruth srparticle.interaction_id = interaction_id; - - std::cout << std::endl; - } //FillTrueGENIEParticle From afb9bb931f0899b67eb697dbb93ca78a16bde4b3 Mon Sep 17 00:00:00 2001 From: Alexander Antonakis Date: Tue, 6 Oct 2026 03:22:33 -0500 Subject: [PATCH 7/8] Removed debug statements --- sbncode/CAFMaker/FillTrue.cxx | 3 --- 1 file changed, 3 deletions(-) diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 01793795d..f32781c89 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -478,12 +478,9 @@ namespace caf { } for(const caf::SRTrueParticle& part: srparticles){ - if (part.pdg == 212212) std::cout << "Weird parent in FillTrueNeutrino" << std::endl; - // save the G4 particles that came from this interaction if(part.interaction_id == (int)i) { if(part.start_process == caf::kG4primary) srneutrino.prim.push_back(part); - //std::cout << "succeeded FillNeutrino check" << std::endl; // total up the deposited energy for(int p = 0; p < 3; ++p) { for (int i_cryo = 0; i_cryo < 2; i_cryo++) { From 7e2d19eb4e6a50427804b4929fd0e27a0ccc098e Mon Sep 17 00:00:00 2001 From: Alexander Antonakis Date: Tue, 6 Oct 2026 03:34:23 -0500 Subject: [PATCH 8/8] Fixed white space --- sbncode/CAFMaker/FillTrue.cxx | 10 ---------- sbncode/CAFMaker/FillTrue.h | 1 - 2 files changed, 11 deletions(-) diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index f32781c89..29dd24b82 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -840,8 +840,6 @@ namespace caf { } } //FillTrueG4Particle - - void FillTrueGENIEParticle(const simb::MCParticle &particle, const std::vector &active_volumes, const std::vector> &tpc_volumes, @@ -1040,7 +1038,6 @@ namespace caf { srparticle.interaction_id = interaction_id; } //FillTrueGENIEParticle - bool IsInitialStateParticle(const simb::MCParticle& particle, const simb::MCTruth& truth) { @@ -1060,8 +1057,6 @@ namespace caf { return false; } - - void FillFakeReco(const std::vector> &mctruths, const std::vector &srparticles, const std::vector> &mctracks, @@ -1886,8 +1881,3 @@ caf::SRTruthMatch MatchSlice2Truth(const std::vector> &hits } return ret; }//Slc2Truth - - - - - diff --git a/sbncode/CAFMaker/FillTrue.h b/sbncode/CAFMaker/FillTrue.h index 084152739..7e6d581a6 100644 --- a/sbncode/CAFMaker/FillTrue.h +++ b/sbncode/CAFMaker/FillTrue.h @@ -81,7 +81,6 @@ namespace caf const std::vector> &neutrinos, caf::SRTrueParticle &srparticle, std::optional new_mother); - // Added for unstable particles that don't propogate to G4 void FillTrueGENIEParticle(const simb::MCParticle &particle, const std::vector &active_volumes,