diff --git a/sbncode/CAFMaker/CAFMaker_module.cc b/sbncode/CAFMaker/CAFMaker_module.cc index b38ed1f6a..be18ac04a 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,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 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(); 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() == 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 + 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; + + // 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(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(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); + } } diff --git a/sbncode/CAFMaker/FillTrue.cxx b/sbncode/CAFMaker/FillTrue.cxx index 4d1d794d0..cb4fddbfd 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; @@ -803,10 +806,18 @@ 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(); + 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; @@ -836,6 +847,231 @@ 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()); + if (particle.PdgCode() == 221) std::cout << "Eta srparticle start process: " << srparticle.start_process << std::endl; + + 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, @@ -1660,3 +1896,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..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,24 @@ 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,