diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index b38ed1f6a..06ab89f3e 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" @@ -1510,14 +1512,136 @@ void CAFMaker::produce(art::Event& evt) noexcept { //####################################################### // Fill truths & fake reco //####################################################### - + caf::SRTruthBranch srtruthbranch; if (mc_particles.isValid()) { + 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; + 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; + } + for (const simb::MCParticle &part: *mc_particles) { + + std::optional missed_parent_id = std::nullopt; + // 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") { + + // 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(); + + for (int ipart = 0; ipart < inventoryTruth->NParticles(); ++ipart) { + const simb::MCParticle& genpart = inventoryTruth->GetParticle(ipart); + 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) { + + // 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) { + 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(), synthetic_missed_parent_id) == mother_ids.end()) { + mother_ids.push_back(synthetic_missed_parent_id); + 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); + + // 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() != 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(); + mother_ids.push_back(synthetic_next_parent_id); + 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); + + 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; + } + } + } + } + } + } // matched Genie Particle + } // found inventoryTruth + } // primary particle with no mother + true_particles.emplace_back(); FillTrueG4Particle(part, @@ -1528,7 +1652,8 @@ void CAFMaker::produce(art::Event& evt) noexcept { *bt_serv, *pi_serv, mctruths, - true_particles.back()); + true_particles.back(), missed_parent_id); + } } diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 4d1d794d0..29dd24b82 100644 --- a/sbncode/CAFMaker/FillTrue.cxx +++ b/sbncode/CAFMaker/FillTrue.cxx @@ -9,6 +9,7 @@ #include "CLHEP/Random/RandGauss.h" +#include #include #include @@ -480,7 +481,6 @@ namespace caf { // 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); - // total up the deposited energy for(int p = 0; p < 3; ++p) { for (int i_cryo = 0; i_cryo < 2; i_cryo++) { @@ -639,7 +639,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; @@ -806,7 +806,11 @@ 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(); + srparticle.parent = static_cast(requested_parent); + } // Set the initial cryostat srparticle.cryostat = -1; @@ -836,6 +840,223 @@ 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, 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; + + 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()); + srparticle.end_process = GetG4ProcessID(particle.EndProcess()); + + // 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; + 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)); + } + + // Set interaction_id to the matched MCTruth + srparticle.interaction_id = interaction_id; + } //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, diff --git a/sbncode/CAFMaker/FillTrue.h b/sbncode/CAFMaker/FillTrue.h index 419ed6ee9..7e6d581a6 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,23 @@ 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 + 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, + int id_offset); + + bool IsInitialStateParticle(const simb::MCParticle& particle, + const simb::MCTruth& truth); void FillMeVPrtlTruth(const evgen::ldm::MeVPrtlTruth &truth, const std::vector &active_volumes,