From 27b5638bba71396fb0b429c77a1fee7afa3e0c53 Mon Sep 17 00:00:00 2001 From: John Back Date: Wed, 9 Oct 2024 18:00:38 +0100 Subject: [PATCH] Truth matching updates for the Pandora LArRecoND CAFs --- cfg/pandora.fcl | 25 ++++- src/reco/PandoraLArRecoNDBranchFiller.cxx | 108 +++++++++++++++------- src/reco/PandoraLArRecoNDBranchFiller.h | 6 +- 3 files changed, 104 insertions(+), 35 deletions(-) diff --git a/cfg/pandora.fcl b/cfg/pandora.fcl index ebb7f14..11c386d 100644 --- a/cfg/pandora.fcl +++ b/cfg/pandora.fcl @@ -1,6 +1,29 @@ #include "NDCAFMaker.fcl" nd_cafmaker: @local::standard_nd_cafmaker -nd_cafmaker.CAFMakerSettings.GHEPFiles: ["/exp/dune/data/users/jback/ND_CAFs/MiniRun5/MiniRun5_1E19_RHC.genie.nu.0000001.GHEP.root"] + +nd_cafmaker.CAFMakerSettings.GHEPFiles: [ + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.nu.0000010.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.nu.0000011.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.nu.0000012.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.nu.0000013.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.nu.0000014.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.nu.0000015.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.nu.0000016.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.nu.0000017.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.nu.0000018.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.nu.0000019.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.rock.0000010.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.rock.0000011.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.rock.0000012.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.rock.0000013.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.rock.0000014.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.rock.0000015.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.rock.0000016.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.rock.0000017.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.rock.0000018.GHEP.root", + "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/MiniRun5_1E19_RHC.genie.rock.0000019.GHEP.root" +] + nd_cafmaker.CAFMakerSettings.PandoraLArRecoNDFile: "/exp/dune/data/users/jback/ND_CAFs/MiniRun6/LArRecoND_MR6_0000001.root" nd_cafmaker.CAFMakerSettings.OutputFile: "CAF_MR6_0000001.root" nd_cafmaker.CAFMakerSettings.Verbosity: VERBOSE diff --git a/src/reco/PandoraLArRecoNDBranchFiller.cxx b/src/reco/PandoraLArRecoNDBranchFiller.cxx index c179d3d..4e4a206 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.cxx +++ b/src/reco/PandoraLArRecoNDBranchFiller.cxx @@ -49,7 +49,7 @@ namespace cafmaker m_LArRecoNDTree->SetBranchAddress("dirZ", &m_dirZVect); m_LArRecoNDTree->SetBranchAddress("energy", &m_energyVect); m_LArRecoNDTree->SetBranchAddress("n3DHits", &m_n3DHitsVect); - m_LArRecoNDTree->SetBranchAddress("mcNuId", &m_mcNuIdVect); + m_LArRecoNDTree->SetBranchAddress("mcVertexId", &m_mcVertexIdVect); m_LArRecoNDTree->SetBranchAddress("isPrimary", &m_isPrimaryVect); m_LArRecoNDTree->SetBranchAddress("mcId", &m_mcIdVect); m_LArRecoNDTree->SetBranchAddress("completeness", &m_completenessVect); @@ -96,12 +96,13 @@ namespace cafmaker // Fill track and shower info. Both use the same clusters (PFOs), and no // distinction is made (yet) to identify which are tracks or showers - FillTracks(sr, nClusters); - FillShowers(sr, nClusters); + FillTracks(sr, nClusters, truthMatch); + FillShowers(sr, nClusters, truthMatch); } // ------------------------------------------------------------------------------ - void PandoraLArRecoNDBranchFiller::FillTracks(caf::StandardRecord& sr, const int nClusters) const + void PandoraLArRecoNDBranchFiller::FillTracks(caf::StandardRecord &sr, const int nClusters, + const TruthMatcher *truthMatch) const { // Create tracks for each PFO (cluster) in the event LOG.VERBOSE() << " Pandora LArRecoND FillTracks using " << nClusters <<" PFO clusters\n"; @@ -157,38 +158,60 @@ namespace cafmaker // Use truth matching info from Pandora's LArContent hierarchy tools. // For LArRecoND MC SpacePoints, we offset the MCId's to make them all unique: - // mcNuId = vertex_id + 10^8, so vertex_id = mcNuId - 10^8 // mcId = traj_id + nuIndex*10^6, where nuIndex = 0 to N-1 neutrinos // nuIndex = int(mcId/10^6), so traj_id = mcId - nuIndex*10^6 - const int mcNuId = (m_mcNuIdVect != nullptr) ? (*m_mcNuIdVect)[i] : 0; - const int vertex_id = (mcNuId > m_nuIdOffset) ? mcNuId - m_nuIdOffset : mcNuId; + // The vertex_id's are the same as those in the HDF5 files + + const long vertex_id = (m_mcVertexIdVect != nullptr) ? (*m_mcVertexIdVect)[i] : 0; const int isPrimary = (m_isPrimaryVect != nullptr) ? (*m_isPrimaryVect)[i] : -1; + const int mcId = (m_mcIdVect != nullptr) ? (*m_mcIdVect)[i] : 0; const int nuIndex = int(mcId/m_maxMCId); const int traj_id = mcId - nuIndex*m_maxMCId; - - caf::TrueParticleID trueID; - trueID.ixn = vertex_id; + + caf::TrueParticleID truePartID; + if (isPrimary == 1) { - trueID.type = caf::TrueParticleID::kPrimary; + truePartID.type = caf::TrueParticleID::kPrimary; + truePartID.part = traj_id; } else if (isPrimary == -1) { - trueID.type = caf::TrueParticleID::kUnknown; + truePartID.type = caf::TrueParticleID::kUnknown; } else { - trueID.type = caf::TrueParticleID::kSecondary; + truePartID.type = caf::TrueParticleID::kSecondary; + } + + try + { + // Get the true interaction in the stack + caf::SRTrueInteraction &srTrueInt = truthMatch->GetTrueInteraction(sr, vertex_id); + const auto predicate = [&srTrueInt](const caf::SRTrueInteraction& ixn) { return ixn.id == srTrueInt.id; }; + // Get the truth interaction index + const int srTrueIntIdx = std::distance(sr.mc.nu.begin(), std::find_if(sr.mc.nu.begin(), sr.mc.nu.end(), predicate)); + truePartID.ixn = srTrueIntIdx; + + // If the particle is not a primary, we might want to create a new particle if it wasn't created originally + if (isPrimary != 1) { + const auto pred = [&traj_id](const caf::SRTrueParticle& part) { return part.G4ID == traj_id; }; + truePartID.part = std::distance(srTrueInt.sec.begin(), + std::find_if(srTrueInt.sec.begin(), srTrueInt.sec.end(), pred)); + } + } + catch (...) + { + LOG.VERBOSE() << "Could not fill some of the truth information" << "\n"; } - trueID.part = traj_id; // Just store the best MC match - std::vector trueIDVect; - trueIDVect.emplace_back(trueID); - track.truth = trueIDVect; + std::vector truePartIDVect; + truePartIDVect.emplace_back(truePartID); + track.truth = truePartIDVect; // Fraction of true MC hits that are captured by the reconstructed cluster const float completeness = (m_completenessVect != nullptr) ? (*m_completenessVect)[i] : 0.0; std::vector truthOverlap; truthOverlap.emplace_back(completeness); track.truthOverlap = truthOverlap; - + // Add track to the record sr.nd.lar.pandora[sliceId].tracks.emplace_back(std::move(track)); sr.nd.lar.pandora[sliceId].ntracks++; @@ -196,7 +219,8 @@ namespace cafmaker } // ------------------------------------------------------------------------------ - void PandoraLArRecoNDBranchFiller::FillShowers(caf::StandardRecord& sr, const int nClusters) const + void PandoraLArRecoNDBranchFiller::FillShowers(caf::StandardRecord &sr, const int nClusters, + const TruthMatcher *truthMatch) const { // Create showers for each PFO (cluster) in the event LOG.VERBOSE() << " Pandora LArRecoND FillShowers using " << nClusters <<" PFO clusters\n"; @@ -231,31 +255,53 @@ namespace cafmaker // Use truth matching info from Pandora's LArContent hierarchy tools. // For LArRecoND MC SpacePoints, we offset the MCId's to make them all unique: - // mcNuId = vertex_id + 10^8, so vertex_id = mcNuId - 10^8 // mcId = traj_id + nuIndex*10^6, where nuIndex = 0 to N-1 neutrinos // nuIndex = int(mcId/10^6), so traj_id = mcId - nuIndex*10^6 - const int mcNuId = (m_mcNuIdVect != nullptr) ? (*m_mcNuIdVect)[i] : 0; - const int vertex_id = (mcNuId > m_nuIdOffset) ? mcNuId - m_nuIdOffset : mcNuId; + // The vertex_id's are the same as those in the HDF5 files + + const long vertex_id = (m_mcVertexIdVect != nullptr) ? (*m_mcVertexIdVect)[i] : 0; const int isPrimary = (m_isPrimaryVect != nullptr) ? (*m_isPrimaryVect)[i] : -1; + const int mcId = (m_mcIdVect != nullptr) ? (*m_mcIdVect)[i] : 0; const int nuIndex = int(mcId/m_maxMCId); const int traj_id = mcId - nuIndex*m_maxMCId; - caf::TrueParticleID trueID; - trueID.ixn = vertex_id; + caf::TrueParticleID truePartID; + if (isPrimary == 1) { - trueID.type = caf::TrueParticleID::kPrimary; + truePartID.type = caf::TrueParticleID::kPrimary; + truePartID.part = traj_id; } else if (isPrimary == -1) { - trueID.type = caf::TrueParticleID::kUnknown; + truePartID.type = caf::TrueParticleID::kUnknown; } else { - trueID.type = caf::TrueParticleID::kSecondary; + truePartID.type = caf::TrueParticleID::kSecondary; + } + + try + { + // Get the true interaction in the stack + caf::SRTrueInteraction &srTrueInt = truthMatch->GetTrueInteraction(sr, vertex_id); + const auto predicate = [&srTrueInt](const caf::SRTrueInteraction& ixn) { return ixn.id == srTrueInt.id; }; + // Get the truth interaction index + const int srTrueIntIdx = std::distance(sr.mc.nu.begin(), std::find_if(sr.mc.nu.begin(), sr.mc.nu.end(), predicate)); + truePartID.ixn = srTrueIntIdx; + + // If the particle is not a primary, we might want to create a new particle if it wasn't created originally + if (isPrimary != 1) { + const auto pred = [&traj_id](const caf::SRTrueParticle& part) { return part.G4ID == traj_id; }; + truePartID.part = std::distance(srTrueInt.sec.begin(), + std::find_if(srTrueInt.sec.begin(), srTrueInt.sec.end(), pred)); + } + } + catch (...) + { + LOG.VERBOSE() << "Could not fill some of the truth information" << "\n"; } - trueID.part = traj_id; // Just store the best MC match - std::vector trueIDVect; - trueIDVect.emplace_back(trueID); - shower.truth = trueIDVect; + std::vector truePartIDVect; + truePartIDVect.emplace_back(truePartID); + shower.truth = truePartIDVect; // Fraction of true MC hits that are captured by the reconstructed cluster const float completeness = (m_completenessVect != nullptr) ? (*m_completenessVect)[i] : 0.0; diff --git a/src/reco/PandoraLArRecoNDBranchFiller.h b/src/reco/PandoraLArRecoNDBranchFiller.h index 661d664..0f77b16 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.h +++ b/src/reco/PandoraLArRecoNDBranchFiller.h @@ -40,8 +40,8 @@ namespace cafmaker const cafmaker::Params &par, const TruthMatcher *truthMatch= nullptr) const override; - void FillTracks(caf::StandardRecord &sr, const int nClusters) const; - void FillShowers(caf::StandardRecord &sr, const int nClusters) const; + void FillTracks(caf::StandardRecord &sr, const int nClusters, const TruthMatcher *truthMatch) const; + void FillShowers(caf::StandardRecord &sr, const int nClusters, const TruthMatcher *truthMatch) const; TFile *m_LArRecoNDFile; TTree *m_LArRecoNDTree; @@ -64,7 +64,7 @@ namespace cafmaker std::vector *m_dirZVect = nullptr; std::vector *m_energyVect = nullptr; std::vector *m_n3DHitsVect = nullptr; - std::vector *m_mcNuIdVect = nullptr; + std::vector *m_mcVertexIdVect = nullptr; std::vector *m_isPrimaryVect = nullptr; std::vector *m_mcIdVect = nullptr; std::vector *m_completenessVect = nullptr;