From 54e23b4e5fc776de0f2d3587277f738985a6bb40 Mon Sep 17 00:00:00 2001 From: John Back Date: Tue, 20 Aug 2024 17:45:56 +0100 Subject: [PATCH 01/10] Add PandoraLArRecoND reco branch filler & example pandora.fcl. Updated ndcaf_setup.sh and build.sh to use a consistent environment. --- build.sh | 20 +- cfg/pandora.fcl | 7 + ndcaf_setup.sh | 5 +- src/Params.h | 1 + src/makeCAF.C | 8 + src/reco/PandoraLArRecoNDBranchFiller.cxx | 306 ++++++++++++++++++++++ src/reco/PandoraLArRecoNDBranchFiller.h | 79 ++++++ 7 files changed, 406 insertions(+), 20 deletions(-) create mode 100644 cfg/pandora.fcl create mode 100644 src/reco/PandoraLArRecoNDBranchFiller.cxx create mode 100644 src/reco/PandoraLArRecoNDBranchFiller.h diff --git a/build.sh b/build.sh index 3f27fbc..c6a26e3 100755 --- a/build.sh +++ b/build.sh @@ -6,23 +6,7 @@ FORCE=no if [[ $# == 1 && x$1 == x-f ]]; then FORCE=yes; fi # set up software -source /cvmfs/dune.opensciencegrid.org/products/dune/setup_dune.sh -setup cmake v3_22_2 -setup gcc v9_3_0 -setup pycurl -setup ifdhc -setup geant4 v4_11_0_p01c -q e20:debug -setup dk2nugenie v01_10_01k -q debug:e20 -setup genie_xsec v3_04_00 -q AR2320i00000:e1000:k250 -setup genie_phyopt v3_04_00 -q dkcharmtau -setup jobsub_client -setup eigen v3_3_5 -setup duneanaobj v03_01_00 -q e20:prof -setup hdf5 v1_12_0b -q e20:prof - -# edep-sim needs to know where a certain GEANT .cmake file is... -G4_cmake_file=`find ${GEANT4_FQ_DIR}/lib -name 'Geant4Config.cmake'` -export Geant4_DIR=`dirname $G4_cmake_file` +source ndcaf_setup.sh # Just use the edep-sim UPS product, don't clone master branch off repos! # Get edep-sim and build it @@ -74,6 +58,6 @@ cd ${TOPDIR} export PYTHONPATH=${PYTHONPATH}:${TOPDIR}/DUNE_ND_GeoEff/lib/ # make tarballs of edep-sim and nusystematics for grid jobs -tar -zcf edep-sim.tar.gz edep-sim +#tar -zcf edep-sim.tar.gz edep-sim #tar -zcf nusystematics.tar.gz nusystematics tar -zcf DUNE_ND_GeoEff.tar.gz DUNE_ND_GeoEff diff --git a/cfg/pandora.fcl b/cfg/pandora.fcl new file mode 100644 index 0000000..f6945cf --- /dev/null +++ b/cfg/pandora.fcl @@ -0,0 +1,7 @@ +#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.PandoraLArRecoNDFile: "/exp/dune/data/users/jback/ND_CAFs/MiniRun5/LArRecoND_MR5_0000001.root" +nd_cafmaker.CAFMakerSettings.OutputFile: "CAF_MR5_0000001.root" +nd_cafmaker.CAFMakerSettings.Verbosity: VERBOSE +nd_cafmaker.CAFMakerSettings.NumEvts: -1 diff --git a/ndcaf_setup.sh b/ndcaf_setup.sh index aaf7397..c1eb823 100644 --- a/ndcaf_setup.sh +++ b/ndcaf_setup.sh @@ -12,10 +12,11 @@ setup eigen v3_3_5 setup duneanaobj v03_05_00 -q debug:e26 setup hdf5 v1_10_5a -q e20 setup fhiclcpp v4_18_04 -q debug:e26 - +setup root v6_28_12 -q e26:p3915:prof +setup genie v3_04_02 -q e26:prof # edep-sim needs to know where a certain GEANT .cmake file is... -G4_cmake_file=`find ${GEANT4_FQ_DIR}/lib -name 'Geant4Config.cmake'` +G4_cmake_file=`find ${GEANT4_FQ_DIR}/lib64 -name 'Geant4Config.cmake'` export Geant4_DIR=`dirname $G4_cmake_file` # edep-sim needs to have the GEANT bin directory in the path diff --git a/src/Params.h b/src/Params.h index 8ec4622..2d9e9a5 100644 --- a/src/Params.h +++ b/src/Params.h @@ -43,6 +43,7 @@ namespace cafmaker fhicl::OptionalAtom tmsRecoFile { fhicl::Name{"TMSRecoFile"}, fhicl::Comment("Input TMS reco .root file") }; fhicl::OptionalAtom sandRecoFile { fhicl::Name{"SANDRecoFile"}, fhicl::Comment("Input SAND reco .root file") }; fhicl::OptionalAtom minervaRecoFile { fhicl::Name{"MINERVARecoFile"}, fhicl::Comment("Input MINERVA reco .root file") }; + fhicl::OptionalAtom pandoraLArRecoNDFile { fhicl::Name{"PandoraLArRecoNDFile"}, fhicl::Comment("Input Pandora LArRecoND .root file") }; // this is optional by way of the default value. Will result in an extra output file if enabled fhicl::Atom makeFlatCAF { fhicl::Name{"MakeFlatCAF"}, fhicl::Comment("Make 'flat' CAF in addition to structured CAF?"), true }; diff --git a/src/makeCAF.C b/src/makeCAF.C index 9b11db6..027f3a0 100644 --- a/src/makeCAF.C +++ b/src/makeCAF.C @@ -26,6 +26,7 @@ #include "reco/NDLArTMSMatchRecoFiller.h" #include "reco/SANDRecoBranchFiller.h" #include "reco/MINERvARecoBranchFiller.h" +#include "reco/PandoraLArRecoNDBranchFiller.h" #include "truth/FillTruth.h" #include "util/GENIEQuiet.h" #include "util/Logger.h" @@ -134,6 +135,13 @@ std::vector> getRecoFillers(const c recoFillers.emplace_back(std::make_unique(sandFile)); std::cout << " SAND\n"; } + // Pandora LArRecoND + std::string pandoraFile; + if (par().cafmaker().pandoraLArRecoNDFile(pandoraFile)) + { + recoFillers.emplace_back(std::make_unique(pandoraFile)); + std::cout << " Pandora LArRecoND\n"; + } // next: did we do TMS reco? std::string tmsFile; diff --git a/src/reco/PandoraLArRecoNDBranchFiller.cxx b/src/reco/PandoraLArRecoNDBranchFiller.cxx new file mode 100644 index 0000000..2d3857d --- /dev/null +++ b/src/reco/PandoraLArRecoNDBranchFiller.cxx @@ -0,0 +1,306 @@ +#include "PandoraLArRecoNDBranchFiller.h" + +namespace cafmaker +{ + + PandoraLArRecoNDBranchFiller::~PandoraLArRecoNDBranchFiller() + { + if (m_LArRecoNDFile && m_LArRecoNDFile->IsOpen()) + { + delete m_LArRecoNDTree; m_LArRecoNDTree = nullptr; + } + delete m_LArRecoNDFile; m_LArRecoNDFile = nullptr; + } + + PandoraLArRecoNDBranchFiller::PandoraLArRecoNDBranchFiller(const std::string &pandoraLArRecoNDFilename) + : IRecoBranchFiller("PandoraLArRecoND"), + m_Triggers(), + m_LastTriggerReqd(m_Triggers.end()) + { + // Open Pandora LArRecoND hierarchy analysis ROOT file + m_LArRecoNDFile = TFile::Open(pandoraLArRecoNDFilename.c_str(), "read"); + if (m_LArRecoNDFile && m_LArRecoNDFile->IsOpen()) { + + LOG.VERBOSE() << " Using PandoraLArRecoND file " << pandoraLArRecoNDFilename << "\n"; + + // Input tree + m_LArRecoNDTree = dynamic_cast(m_LArRecoNDFile->Get("LArRecoND")); + if (!m_LArRecoNDTree) { + LOG.FATAL() << "Did not find LArRecoND tree in input file " << pandoraLArRecoNDFilename << "\n"; + throw; + } + + // Set branch addresses + m_LArRecoNDTree->SetBranchAddress("event", &m_eventId); + m_LArRecoNDTree->SetBranchAddress("run", &m_run); + m_LArRecoNDTree->SetBranchAddress("subRun", &m_subRun); + m_LArRecoNDTree->SetBranchAddress("startTime", &m_startTime); + m_LArRecoNDTree->SetBranchAddress("sliceId", &m_sliceIdVect); + m_LArRecoNDTree->SetBranchAddress("startX", &m_startXVect); + m_LArRecoNDTree->SetBranchAddress("startY", &m_startYVect); + m_LArRecoNDTree->SetBranchAddress("startZ", &m_startZVect); + m_LArRecoNDTree->SetBranchAddress("endX", &m_endXVect); + m_LArRecoNDTree->SetBranchAddress("endY", &m_endYVect); + m_LArRecoNDTree->SetBranchAddress("endZ", &m_endZVect); + m_LArRecoNDTree->SetBranchAddress("dirX", &m_dirXVect); + m_LArRecoNDTree->SetBranchAddress("dirY", &m_dirYVect); + m_LArRecoNDTree->SetBranchAddress("dirZ", &m_dirZVect); + m_LArRecoNDTree->SetBranchAddress("energy", &m_energyVect); + m_LArRecoNDTree->SetBranchAddress("n3DHits", &m_n3DHitsVect); + m_LArRecoNDTree->SetBranchAddress("length1", &m_length1Vect); + m_LArRecoNDTree->SetBranchAddress("mcNuId", &m_mcNuIdVect); + m_LArRecoNDTree->SetBranchAddress("isPrimary", &m_isPrimaryVect); + m_LArRecoNDTree->SetBranchAddress("mcId", &m_mcIdVect); + m_LArRecoNDTree->SetBranchAddress("completeness", &m_completenessVect); + + // We have setup the input tree + SetConfigured(true); + } + } + + // Copy all of the Pandora LArRecoND info to the PandoraLArRecoND branch of the StandardRecord object + void PandoraLArRecoNDBranchFiller::_FillRecoBranches(const Trigger &trigger, + caf::StandardRecord &sr, + const cafmaker::Params &par, + const TruthMatcher *truthMatch) const + { + // Figure out where in our list of triggers this event index is. + // We should always be looking forwards, since we expect to be traversing in that direction + auto it_start = (m_LastTriggerReqd == m_Triggers.end()) ? m_Triggers.cbegin() : m_LastTriggerReqd; + auto itTrig = std::find(it_start, m_Triggers.cend(), trigger); + if (itTrig == m_Triggers.end()) + { + LOG.FATAL() << " Reco branch filler '" << GetName() << "' could not find trigger with evtID == " + << trigger.evtID << "! Abort.\n"; + abort(); + } + std::size_t idx = std::distance(m_Triggers.cbegin(), itTrig); + + LOG.VERBOSE() << " Reco branch filler '" << GetName() << "', trigger.evtID == " << trigger.evtID + << ", internal evt idx = " << idx << ".\n"; + + // Get the event entry + m_LArRecoNDTree->GetEntry(idx); + + // Set the event and run numbers + sr.meta.nd_lar.enabled = true; + sr.meta.nd_lar.event = m_eventId; + sr.meta.nd_lar.run = m_run; + sr.meta.nd_lar.subrun = m_subRun; + + // Set the size for the Pandora NDLAr standard record for this trigger/event + const int nClusters = (m_sliceIdVect != nullptr) ? m_sliceIdVect->size() : 0; + sr.nd.lar.pandora.resize(nClusters); + sr.nd.lar.npandora = sr.nd.lar.pandora.size(); + + // 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); + } + + // ------------------------------------------------------------------------------ + void PandoraLArRecoNDBranchFiller::FillTracks(caf::StandardRecord& sr, const int nClusters) const + { + // Create tracks for each PFO (cluster) in the event + LOG.VERBOSE() << " Pandora LArRecoND FillTracks using " << nClusters <<" PFO clusters\n"; + + for (int i = 0; i < nClusters; i++) + { + // Slice id of the PFO cluster + const int sliceId = (*m_sliceIdVect)[i]; + + // Create standard record track + caf::SRTrack track; + // Starting position (vertex or first hit location) + const float startX = (m_startXVect != nullptr) ? (*m_startXVect)[i] : 0.0; + const float startY = (m_startYVect != nullptr) ? (*m_startYVect)[i] : 0.0; + const float startZ = (m_startZVect != nullptr) ? (*m_startZVect)[i] : 0.0; + track.start = caf::SRVector3D(startX, startY, startZ); + + // End position (length along principal direction or last hit location) + const float endX = (m_endXVect != nullptr) ? (*m_endXVect)[i] : 0.0; + const float endY = (m_endYVect != nullptr) ? (*m_endYVect)[i] : 0.0; + const float endZ = (m_endZVect != nullptr) ? (*m_endZVect)[i] : 0.0; + track.end = caf::SRVector3D(endX, endY, endZ); + + // Principal axis direction + const float dirX = (m_dirXVect != nullptr) ? (*m_dirXVect)[i] : 0.0; + const float dirY = (m_dirYVect != nullptr) ? (*m_dirYVect)[i] : 0.0; + const float dirZ = (m_dirZVect != nullptr) ? (*m_dirZVect)[i] : 0.0; + track.dir = caf::SRVector3D(dirX, dirY, dirZ); + track.enddir = caf::SRVector3D(dirX, dirY, dirZ); + + // Energy (GeV) + const float energy = (m_energyVect != nullptr) ? (*m_energyVect)[i] : 0.0; + track.Evis = energy; + track.E = energy; + + // Total number of 3D hits in the cluster + const int n3DHits = (m_n3DHitsVect != nullptr) ? (*m_n3DHitsVect)[i] : 0; + track.qual = n3DHits*1.0; + + // Cluster length along principal axis direction (cm) + const float length1 = (m_length1Vect != nullptr) ? (*m_length1Vect)[i] : 0.0; + track.len_cm = length1; + + // Cluster length multiplied by LAr density (g/cm2) + track.len_gcm2 = length1*m_LArRho; + + // Use truth matching info from Pandora's LArContent hierarchy tools. + // For LArRecoND MC SpacePoints, we offset the MCId's to make them all unique: + // nuId = orig_nuId + 10^8, so orig_nuId = nuId - 10^8 + // mcId = orig_mcId + nuIndex*10^6, where nuIndex = 0 to N-1 neutrinos + // nuIndex = int(mcId/10^6), so orig_mcId = mcId - nuIndex*10^6 + // Use the original MCId values + const int mcNuId = (m_mcNuIdVect != nullptr) ? (*m_mcNuIdVect)[i] : 0; + const int origMCNuId = (mcNuId > m_nuIdOffset) ? mcNuId - m_nuIdOffset : mcNuId; + 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 origMCId = mcId - nuIndex*m_maxMCId; + + caf::TrueParticleID trueID; + trueID.ixn = origMCNuId; + if (isPrimary == 1) { + trueID.type = caf::TrueParticleID::kPrimary; + } else if (isPrimary == -1) { + trueID.type = caf::TrueParticleID::kUnknown; + } else { + trueID.type = caf::TrueParticleID::kSecondary; + } + trueID.part = origMCId; + + // Just store the best MC match + std::vector trueIDVect; + trueIDVect.emplace_back(trueID); + track.truth = trueIDVect; + + // 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++; + } + } + + // ------------------------------------------------------------------------------ + void PandoraLArRecoNDBranchFiller::FillShowers(caf::StandardRecord& sr, const int nClusters) const + { + // Create showers for each PFO (cluster) in the event + LOG.VERBOSE() << " Pandora LArRecoND FillShowers using " << nClusters <<" PFO clusters\n"; + + for (int i = 0; i < nClusters; i++) + { + // Slice id of the PFO cluster + const int sliceId = (*m_sliceIdVect)[i]; + + // Create standard record shower + caf::SRShower shower; + // Starting position + const float startX = (m_startXVect != nullptr) ? (*m_startXVect)[i] : 0.0; + const float startY = (m_startYVect != nullptr) ? (*m_startYVect)[i] : 0.0; + const float startZ = (m_startZVect != nullptr) ? (*m_startZVect)[i] : 0.0; + shower.start = caf::SRVector3D(startX, startY, startZ); + + // Principal axis direction + const float dirX = (m_dirXVect != nullptr) ? (*m_dirXVect)[i] : 0.0; + const float dirY = (m_dirYVect != nullptr) ? (*m_dirYVect)[i] : 0.0; + const float dirZ = (m_dirZVect != nullptr) ? (*m_dirZVect)[i] : 0.0; + shower.direction = caf::SRVector3D(dirX, dirY, dirZ); + + // Energy (GeV) + const float energy = (m_energyVect != nullptr) ? (*m_energyVect)[i] : 0.0; + shower.Evis = energy; + + // Use truth matching info from Pandora's LArContent hierarchy tools. + // For LArRecoND MC SpacePoints, we offset the MCId's to make them all unique: + // nuId = orig_nuId + 10^8, so orig_nuId = nuId - 10^8 + // mcId = orig_mcId + nuIndex*10^6, where nuIndex = 0 to N-1 neutrinos + // nuIndex = int(mcId/10^6), so orig_mcId = mcId - nuIndex*10^6 + // Use the original MCId values + const int mcNuId = (m_mcNuIdVect != nullptr) ? (*m_mcNuIdVect)[i] : 0; + const int origMCNuId = (mcNuId > m_nuIdOffset) ? mcNuId - m_nuIdOffset : mcNuId; + 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 origMCId = mcId - nuIndex*m_maxMCId; + + caf::TrueParticleID trueID; + trueID.ixn = origMCNuId; + if (isPrimary == 1) { + trueID.type = caf::TrueParticleID::kPrimary; + } else if (isPrimary == -1) { + trueID.type = caf::TrueParticleID::kUnknown; + } else { + trueID.type = caf::TrueParticleID::kSecondary; + } + trueID.part = origMCId; + + // Just store the best MC match + std::vector trueIDVect; + trueIDVect.emplace_back(trueID); + shower.truth = trueIDVect; + + // 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); + shower.truthOverlap = truthOverlap; + + // Add shower to the record + sr.nd.lar.pandora[sliceId].showers.emplace_back(std::move(shower)); + sr.nd.lar.pandora[sliceId].nshowers++; + } + } + + // ------------------------------------------------------------------------------ + std::deque PandoraLArRecoNDBranchFiller::GetTriggers(int triggerType) const + { + if (m_Triggers.empty()) + { + const int nEvents = m_LArRecoNDTree->GetEntries(); + LOG.DEBUG() << "Loading triggers with type " << triggerType << " within branch filler '" << GetName() + << "' from " << nEvents << " Pandora LArRecoND tree entries:\n"; + + m_Triggers.reserve(nEvents); + for (int entry = 0; entry < nEvents; entry++) + { + m_LArRecoNDTree->GetEntry(entry); + + m_Triggers.emplace_back(); + Trigger &trig = m_Triggers.back(); + // Event number + trig.evtID = m_eventId; + // Pandora SpacePoint (SP) H5Flow-to-ROOT format doesn't store trigger type, so just select "all" + trig.triggerType = -1; + // Pandora SP format uses timestamp ticks from /charge/events/ts_start + trig.triggerTime_s = m_startTime; + // Use the same (integer) time for now + trig.triggerTime_ns = m_startTime; + + LOG.VERBOSE() << " added trigger: evtID = " << trig.evtID + << ", triggerType = " << trig.triggerType + << ", triggerTime_s = " << trig.triggerTime_s + << ", triggerTime_ns = " << trig.triggerTime_ns + << "\n"; + } + // Since we just modified the list, any iterators have been invalidated + m_LastTriggerReqd = m_Triggers.end(); + } + + std::deque triggers; + for (const Trigger &trigger : m_Triggers) + { + if (triggerType < 0 || triggerType == m_Triggers.back().triggerType) + triggers.push_back(trigger); + } + + return triggers; + } + +} // end namespace diff --git a/src/reco/PandoraLArRecoNDBranchFiller.h b/src/reco/PandoraLArRecoNDBranchFiller.h new file mode 100644 index 0000000..ccd19e7 --- /dev/null +++ b/src/reco/PandoraLArRecoNDBranchFiller.h @@ -0,0 +1,79 @@ +/// Fill Pandora LArRecoND branches +/// +/// \author John Back +/// \date Aug 2024 + +#ifndef ND_CAFMAKER_PandoraLArRecoNDBranchFiller_H +#define ND_CAFMAKER_PandoraLArRecoNDBranchFiller_H + +#include +#include +#include + +// The virtual base class +#include "reco/IRecoBranchFiller.h" +#include "truth/FillTruth.h" + +// ROOT headers +#include "TFile.h" +#include "TTree.h" + +// duneanaobj +#include "duneanaobj/StandardRecord/StandardRecord.h" + +namespace cafmaker +{ + class PandoraLArRecoNDBranchFiller : public cafmaker::IRecoBranchFiller + { + public: + PandoraLArRecoNDBranchFiller(const std::string & pandoraLArRecoNDFilename); + + std::deque GetTriggers(int triggerType) const override; + + ~PandoraLArRecoNDBranchFiller(); + + private: + void _FillRecoBranches(const Trigger &trigger, + caf::StandardRecord &sr, + 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; + + TFile *m_LArRecoNDFile; + TTree *m_LArRecoNDTree; + + int m_eventId; + int m_run; + int m_subRun; + int m_startTime; + std::vector *m_sliceIdVect = nullptr; + std::vector *m_startXVect = nullptr; + std::vector *m_startYVect = nullptr; + std::vector *m_startZVect = nullptr; + std::vector *m_endXVect = nullptr; + std::vector *m_endYVect = nullptr; + std::vector *m_endZVect = nullptr; + std::vector *m_dirXVect = nullptr; + std::vector *m_dirYVect = nullptr; + std::vector *m_dirZVect = nullptr; + std::vector *m_energyVect = nullptr; + std::vector *m_n3DHitsVect = nullptr; + std::vector *m_length1Vect = nullptr; + std::vector *m_mcNuIdVect = nullptr; + std::vector *m_isPrimaryVect = nullptr; + std::vector *m_mcIdVect = nullptr; + std::vector *m_completenessVect = nullptr; + + float m_LArRho{1.3973f}; // LAr density (g/cm3) + + int m_nuIdOffset{100000000}; + int m_maxMCId{1000000}; + + mutable std::vector m_Triggers; + mutable decltype(m_Triggers)::const_iterator m_LastTriggerReqd; ///< the last trigger requested using _FillRecoBranches + }; + +} +#endif //ND_CAFMAKER_PandoraLArRecoNDBranchFiller_H From 58cd6b9161d1c229899214a5e3a541ffd58511c8 Mon Sep 17 00:00:00 2001 From: John Back Date: Thu, 22 Aug 2024 15:44:19 +0100 Subject: [PATCH 02/10] Set Pandora track length using start and end points. --- src/reco/PandoraLArRecoNDBranchFiller.cxx | 13 +++++++------ src/reco/PandoraLArRecoNDBranchFiller.h | 1 - 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/src/reco/PandoraLArRecoNDBranchFiller.cxx b/src/reco/PandoraLArRecoNDBranchFiller.cxx index 2d3857d..7e6ea6e 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.cxx +++ b/src/reco/PandoraLArRecoNDBranchFiller.cxx @@ -47,7 +47,6 @@ namespace cafmaker m_LArRecoNDTree->SetBranchAddress("dirZ", &m_dirZVect); m_LArRecoNDTree->SetBranchAddress("energy", &m_energyVect); m_LArRecoNDTree->SetBranchAddress("n3DHits", &m_n3DHitsVect); - m_LArRecoNDTree->SetBranchAddress("length1", &m_length1Vect); m_LArRecoNDTree->SetBranchAddress("mcNuId", &m_mcNuIdVect); m_LArRecoNDTree->SetBranchAddress("isPrimary", &m_isPrimaryVect); m_LArRecoNDTree->SetBranchAddress("mcId", &m_mcIdVect); @@ -118,7 +117,7 @@ namespace cafmaker const float startZ = (m_startZVect != nullptr) ? (*m_startZVect)[i] : 0.0; track.start = caf::SRVector3D(startX, startY, startZ); - // End position (length along principal direction or last hit location) + // End position const float endX = (m_endXVect != nullptr) ? (*m_endXVect)[i] : 0.0; const float endY = (m_endYVect != nullptr) ? (*m_endYVect)[i] : 0.0; const float endZ = (m_endZVect != nullptr) ? (*m_endZVect)[i] : 0.0; @@ -140,12 +139,14 @@ namespace cafmaker const int n3DHits = (m_n3DHitsVect != nullptr) ? (*m_n3DHitsVect)[i] : 0; track.qual = n3DHits*1.0; - // Cluster length along principal axis direction (cm) - const float length1 = (m_length1Vect != nullptr) ? (*m_length1Vect)[i] : 0.0; - track.len_cm = length1; + // Cluster length from start and end points + const float dX = endX - startX; + const float dY = endY - startY; + const float dZ = endZ - startZ; + track.len_cm = sqrt(dX*dX + dY*dY + dZ*dZ); // Cluster length multiplied by LAr density (g/cm2) - track.len_gcm2 = length1*m_LArRho; + track.len_gcm2 = track.len_cm*m_LArRho; // Use truth matching info from Pandora's LArContent hierarchy tools. // For LArRecoND MC SpacePoints, we offset the MCId's to make them all unique: diff --git a/src/reco/PandoraLArRecoNDBranchFiller.h b/src/reco/PandoraLArRecoNDBranchFiller.h index ccd19e7..b41706c 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.h +++ b/src/reco/PandoraLArRecoNDBranchFiller.h @@ -60,7 +60,6 @@ namespace cafmaker std::vector *m_dirZVect = nullptr; std::vector *m_energyVect = nullptr; std::vector *m_n3DHitsVect = nullptr; - std::vector *m_length1Vect = nullptr; std::vector *m_mcNuIdVect = nullptr; std::vector *m_isPrimaryVect = nullptr; std::vector *m_mcIdVect = nullptr; From 4bd318d2cb14bf44559e5d36e9e82e80cda64a64 Mon Sep 17 00:00:00 2001 From: John Back Date: Tue, 27 Aug 2024 17:16:10 +0100 Subject: [PATCH 03/10] Use isShower flag to check if PFOs are showers or tracks. --- src/reco/PandoraLArRecoNDBranchFiller.cxx | 11 +++++++++++ src/reco/PandoraLArRecoNDBranchFiller.h | 1 + 2 files changed, 12 insertions(+) diff --git a/src/reco/PandoraLArRecoNDBranchFiller.cxx b/src/reco/PandoraLArRecoNDBranchFiller.cxx index 7e6ea6e..8a91054 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.cxx +++ b/src/reco/PandoraLArRecoNDBranchFiller.cxx @@ -35,6 +35,7 @@ namespace cafmaker m_LArRecoNDTree->SetBranchAddress("run", &m_run); m_LArRecoNDTree->SetBranchAddress("subRun", &m_subRun); m_LArRecoNDTree->SetBranchAddress("startTime", &m_startTime); + m_LArRecoNDTree->SetBranchAddress("isShower", &m_isShowerVect); m_LArRecoNDTree->SetBranchAddress("sliceId", &m_sliceIdVect); m_LArRecoNDTree->SetBranchAddress("startX", &m_startXVect); m_LArRecoNDTree->SetBranchAddress("startY", &m_startYVect); @@ -106,6 +107,11 @@ namespace cafmaker for (int i = 0; i < nClusters; i++) { + // Check that the PFO is a track and not a shower + const int isShower = (*m_isShowerVect)[i]; + if (isShower == 1) + continue; + // Slice id of the PFO cluster const int sliceId = (*m_sliceIdVect)[i]; @@ -197,6 +203,11 @@ namespace cafmaker for (int i = 0; i < nClusters; i++) { + // Check that the PFO is a shower + const int isShower = (*m_isShowerVect)[i]; + if (isShower == 0) + continue; + // Slice id of the PFO cluster const int sliceId = (*m_sliceIdVect)[i]; diff --git a/src/reco/PandoraLArRecoNDBranchFiller.h b/src/reco/PandoraLArRecoNDBranchFiller.h index b41706c..fe02c43 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.h +++ b/src/reco/PandoraLArRecoNDBranchFiller.h @@ -48,6 +48,7 @@ namespace cafmaker int m_run; int m_subRun; int m_startTime; + std::vector *m_isShowerVect = nullptr; std::vector *m_sliceIdVect = nullptr; std::vector *m_startXVect = nullptr; std::vector *m_startYVect = nullptr; From 8f9264fdae08b07c3749ed50bc9e4771e8b1f01a Mon Sep 17 00:00:00 2001 From: John Back Date: Thu, 3 Oct 2024 14:36:55 +0100 Subject: [PATCH 04/10] Add unix_ts trigger time and update MCId variable names --- cfg/pandora.fcl | 4 +-- src/reco/PandoraLArRecoNDBranchFiller.cxx | 39 +++++++++++------------ src/reco/PandoraLArRecoNDBranchFiller.h | 1 + 3 files changed, 22 insertions(+), 22 deletions(-) diff --git a/cfg/pandora.fcl b/cfg/pandora.fcl index f6945cf..ebb7f14 100644 --- a/cfg/pandora.fcl +++ b/cfg/pandora.fcl @@ -1,7 +1,7 @@ #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.PandoraLArRecoNDFile: "/exp/dune/data/users/jback/ND_CAFs/MiniRun5/LArRecoND_MR5_0000001.root" -nd_cafmaker.CAFMakerSettings.OutputFile: "CAF_MR5_0000001.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 nd_cafmaker.CAFMakerSettings.NumEvts: -1 diff --git a/src/reco/PandoraLArRecoNDBranchFiller.cxx b/src/reco/PandoraLArRecoNDBranchFiller.cxx index 8a91054..c179d3d 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.cxx +++ b/src/reco/PandoraLArRecoNDBranchFiller.cxx @@ -34,6 +34,7 @@ namespace cafmaker m_LArRecoNDTree->SetBranchAddress("event", &m_eventId); m_LArRecoNDTree->SetBranchAddress("run", &m_run); m_LArRecoNDTree->SetBranchAddress("subRun", &m_subRun); + m_LArRecoNDTree->SetBranchAddress("unixTime", &m_unixTime); m_LArRecoNDTree->SetBranchAddress("startTime", &m_startTime); m_LArRecoNDTree->SetBranchAddress("isShower", &m_isShowerVect); m_LArRecoNDTree->SetBranchAddress("sliceId", &m_sliceIdVect); @@ -156,19 +157,18 @@ 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: - // nuId = orig_nuId + 10^8, so orig_nuId = nuId - 10^8 - // mcId = orig_mcId + nuIndex*10^6, where nuIndex = 0 to N-1 neutrinos - // nuIndex = int(mcId/10^6), so orig_mcId = mcId - nuIndex*10^6 - // Use the original MCId values + // 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 origMCNuId = (mcNuId > m_nuIdOffset) ? mcNuId - m_nuIdOffset : mcNuId; + const int vertex_id = (mcNuId > m_nuIdOffset) ? mcNuId - m_nuIdOffset : mcNuId; 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 origMCId = mcId - nuIndex*m_maxMCId; + const int traj_id = mcId - nuIndex*m_maxMCId; caf::TrueParticleID trueID; - trueID.ixn = origMCNuId; + trueID.ixn = vertex_id; if (isPrimary == 1) { trueID.type = caf::TrueParticleID::kPrimary; } else if (isPrimary == -1) { @@ -176,7 +176,7 @@ namespace cafmaker } else { trueID.type = caf::TrueParticleID::kSecondary; } - trueID.part = origMCId; + trueID.part = traj_id; // Just store the best MC match std::vector trueIDVect; @@ -231,19 +231,18 @@ 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: - // nuId = orig_nuId + 10^8, so orig_nuId = nuId - 10^8 - // mcId = orig_mcId + nuIndex*10^6, where nuIndex = 0 to N-1 neutrinos - // nuIndex = int(mcId/10^6), so orig_mcId = mcId - nuIndex*10^6 - // Use the original MCId values + // 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 origMCNuId = (mcNuId > m_nuIdOffset) ? mcNuId - m_nuIdOffset : mcNuId; + const int vertex_id = (mcNuId > m_nuIdOffset) ? mcNuId - m_nuIdOffset : mcNuId; 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 origMCId = mcId - nuIndex*m_maxMCId; + const int traj_id = mcId - nuIndex*m_maxMCId; caf::TrueParticleID trueID; - trueID.ixn = origMCNuId; + trueID.ixn = vertex_id; if (isPrimary == 1) { trueID.type = caf::TrueParticleID::kPrimary; } else if (isPrimary == -1) { @@ -251,7 +250,7 @@ namespace cafmaker } else { trueID.type = caf::TrueParticleID::kSecondary; } - trueID.part = origMCId; + trueID.part = traj_id; // Just store the best MC match std::vector trueIDVect; @@ -290,10 +289,10 @@ namespace cafmaker trig.evtID = m_eventId; // Pandora SpacePoint (SP) H5Flow-to-ROOT format doesn't store trigger type, so just select "all" trig.triggerType = -1; - // Pandora SP format uses timestamp ticks from /charge/events/ts_start - trig.triggerTime_s = m_startTime; - // Use the same (integer) time for now - trig.triggerTime_ns = m_startTime; + // unix_ts trigger time (seconds) + trig.triggerTime_s = m_unixTime; + // ts_start ticks (0.1 microseconds) converted to nanoseconds + trig.triggerTime_ns = m_startTime * 100; LOG.VERBOSE() << " added trigger: evtID = " << trig.evtID << ", triggerType = " << trig.triggerType diff --git a/src/reco/PandoraLArRecoNDBranchFiller.h b/src/reco/PandoraLArRecoNDBranchFiller.h index fe02c43..b9e2bfa 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.h +++ b/src/reco/PandoraLArRecoNDBranchFiller.h @@ -47,6 +47,7 @@ namespace cafmaker int m_eventId; int m_run; int m_subRun; + int m_unixTime; int m_startTime; std::vector *m_isShowerVect = nullptr; std::vector *m_sliceIdVect = nullptr; From 611c0118051e8811294951663535b33e8c57b3ff Mon Sep 17 00:00:00 2001 From: John Back Date: Thu, 3 Oct 2024 15:53:15 +0100 Subject: [PATCH 05/10] Add RecoFillerType to PandoraLArRecoNDBranchFiller. --- src/reco/PandoraLArRecoNDBranchFiller.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/reco/PandoraLArRecoNDBranchFiller.h b/src/reco/PandoraLArRecoNDBranchFiller.h index b9e2bfa..661d664 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.h +++ b/src/reco/PandoraLArRecoNDBranchFiller.h @@ -30,6 +30,8 @@ namespace cafmaker std::deque GetTriggers(int triggerType) const override; + RecoFillerType FillerType() const override { return RecoFillerType::BaseReco; } + ~PandoraLArRecoNDBranchFiller(); private: From 27b5638bba71396fb0b429c77a1fee7afa3e0c53 Mon Sep 17 00:00:00 2001 From: John Back Date: Wed, 9 Oct 2024 18:00:38 +0100 Subject: [PATCH 06/10] 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; From 03b711843289597260ed6930cccf28ab68de7de0 Mon Sep 17 00:00:00 2001 From: John Back Date: Wed, 30 Oct 2024 16:08:23 +0000 Subject: [PATCH 07/10] Update src/makeCAF.C from main branch --- src/makeCAF.C | 135 +++++++------------------------------------------- 1 file changed, 18 insertions(+), 117 deletions(-) diff --git a/src/makeCAF.C b/src/makeCAF.C index a09603b..460ed83 100644 --- a/src/makeCAF.C +++ b/src/makeCAF.C @@ -1,11 +1,5 @@ #include #include -#include -#include -#include -#include -#include -#include #include "boost/program_options/options_description.hpp" #include "boost/program_options/parsers.hpp" @@ -33,9 +27,9 @@ #include "reco/NDLArTMSMatchRecoFiller.h" #include "reco/NDLArMINERvAMatchRecoFiller.h" -#include "reco/PandoraLArRecoNDBranchFiller.h" #include "reco/SANDRecoBranchFiller.h" #include "truth/FillTruth.h" +#include "beam/IFBeam.h" #include "util/GENIEQuiet.h" #include "util/Logger.h" #include "util/Progress.h" @@ -144,14 +138,6 @@ std::vector> getRecoFillers(const c std::cout << " SAND\n"; } - // Pandora LArRecoND - std::string pandoraFile; - if (par().cafmaker().pandoraLArRecoNDFile(pandoraFile)) - { - recoFillers.emplace_back(std::make_unique(pandoraFile)); - std::cout << " Pandora LArRecoND\n"; - } - // next: did we do TMS reco? std::string tmsFile; if (par().cafmaker().tmsRecoFile(tmsFile)) @@ -179,7 +165,7 @@ std::vector> getRecoFillers(const c { recoFillers.emplace_back(std::make_unique(par().cafmaker().trackMatchExtrapolatedZ(), par().cafmaker().trackMatchdX(), par().cafmaker().trackMatchdY(), par().cafmaker().trackMatchdThetaX(), par().cafmaker().trackMatchdThetaY())); std::cout << " ND-LAr + MINERvA matching\n"; - } + } // for now all the fillers get the same threshold. // if we decide we need to do it differently later // we can adjust the FCL params... @@ -347,104 +333,6 @@ buildTriggerList(std::map; - -bool loadBeamSpills(const std::string& filename, BeamSpills& beam_spills) { - std::ifstream file(filename); - if (!file.is_open()) { - std::cerr << "Could not open the file " << filename << "\n"; - return false; - } - - std::string line; - while (std::getline(file, line)) { - std::istringstream iss(line); - double time; - double pot; - if (iss >> time >> pot) - beam_spills[time] = pot; - } - file.close(); - return true; -} - -double GetTriggerTime(const cafmaker::Trigger& trigger) { - return trigger.triggerTime_s + 1e-9 * trigger.triggerTime_ns; -} - -//Return pot from text file if all of the trigger times match the beam times within a certain dt, else fill given pot -double getPOT(const cafmaker::Params& par, std::vector>& groupedTrigger, int ii) { - std::string potFile; - BeamSpills beam_spills; - - double pot = 0.0; - - if (par().cafmaker().POTFile(potFile) && loadBeamSpills(potFile, beam_spills)) { //If there is a POT file in config that is readable, check if all trigger times match any of the beam times - auto it = std::find_if(beam_spills.begin(), beam_spills.end(), - [par, &groupedTrigger](const auto& spill) { - return std::all_of(groupedTrigger.cbegin(), groupedTrigger.cend(), - [par, &spill](const auto& groupedTrigger) { - return std::abs(GetTriggerTime(groupedTrigger.second) - spill.first) < par().cafmaker().beamMatchDT(); //fixme: shouldn't be abs after 2x2 trigger times are fixed - }); - }); - - if (it != beam_spills.end()) { - pot = it->second; - } - else { //Check if any of the triggers match the beam time if there is no beam time match in the previous search - bool any_matched = false; - std::vector> matched_triggers; - std::vector> unmatched_triggers; - - for (auto trig : groupedTrigger) { - bool matched = std::any_of(beam_spills.cbegin(), beam_spills.cend(), - [par, &trig](const auto& spill) { - return std::abs(GetTriggerTime(trig.second) - spill.first) < par().cafmaker().beamMatchDT(); //fixme: shouldn't be abs after 2x2 trigger times are fixed - }); - - if (matched) { - any_matched = true; - matched_triggers.push_back(trig); - } - else { - unmatched_triggers.push_back(trig); - } - } - - auto LOG = [&]() -> const cafmaker::Logger & { return cafmaker::LOG_S("Beam spill matching"); }; - std::stringstream log_message; - - if (any_matched) { //If only some of the trigger times in a grouped trigger match the beam spill, abort! - log_message << "Only some triggers match beam spill for trigger group " << ii << ":\n" - << "Matched triggers: \n"; - for (auto trig : matched_triggers) - log_message << std::fixed << trig.first->GetName() << " " << GetTriggerTime(trig.second) << "\n"; - - log_message << "Unmatched triggers: \n"; - for (auto trig : unmatched_triggers) - log_message << std::fixed << trig.first->GetName() << " " << GetTriggerTime(trig.second) << "\n"; - - LOG().ERROR() << log_message.str() << "\n"; - std::abort(); - } - else { //If none of the trigger times match give a warning - log_message << "No matching spill found for trigger group " << ii << " with triggers: \n"; - for (auto trig : groupedTrigger) - log_message << std::fixed << trig.first->GetName() << " " << GetTriggerTime(trig.second) << "\n"; - LOG().WARNING() << log_message.str() << "\n"; - } - } - } - else { //else get POT from config - pot = par().runInfo().POTPerSpill() * 1e13; - } - - return pot; -} // ------------------------------------------------- // main loop function void loop(CAF &caf, @@ -485,6 +373,11 @@ void loop(CAF &caf, abort(); } + bool useIFBeam = false; + if (ghepFilenames.empty() && edepsimFilename.empty() && !par().cafmaker().ForceDisableIFBeam()) useIFBeam = true; + + cafmaker::IFBeam beamManager(groupedTriggers, useIFBeam); //initialize IFBeam manager if data and when IFBeam is not force disabled + // Main event loop cafmaker::Progress progBar("Processing " + std::to_string(N - start) + " triggers"); for( int ii = start; ii < start + N; ++ii ) @@ -514,14 +407,22 @@ void loop(CAF &caf, filler->FillRecoBranches(groupedTriggers[ii][0].second, caf.sr, par, &truthMatcher); } } - + //Fill POT - double pot = getPOT(par, groupedTriggers[ii], ii); + double pot = 0.0; + if (useIFBeam) + { + pot = beamManager.getPOT(par, groupedTriggers[ii], ii); + } + else + { + pot = par().runInfo().POTPerSpill() * 1e13; + caf.sr.beam.ismc = true; + } if (std::isnan(caf.pot)) caf.pot = 0; caf.pot += pot; caf.sr.beam.pulsepot = pot; - caf.sr.beam.ismc = par().cafmaker().POTFile.hasValue(); // fixme: when we have proper IFDB interface, should use the same mechanism as however we decide when to use that caf.fill(); } progBar.Done(); From 9d172eea7e0ef32f1a71a5c3cbecde88ee395e99 Mon Sep 17 00:00:00 2001 From: John Back Date: Wed, 30 Oct 2024 16:51:19 +0000 Subject: [PATCH 08/10] Update MC Ids used for Pandora CAF truth matching --- src/makeCAF.C | 19 +++++--- src/reco/PandoraLArRecoNDBranchFiller.cxx | 55 +++++++++++------------ src/reco/PandoraLArRecoNDBranchFiller.h | 4 +- 3 files changed, 42 insertions(+), 36 deletions(-) diff --git a/src/makeCAF.C b/src/makeCAF.C index 460ed83..68deaac 100644 --- a/src/makeCAF.C +++ b/src/makeCAF.C @@ -27,6 +27,7 @@ #include "reco/NDLArTMSMatchRecoFiller.h" #include "reco/NDLArMINERvAMatchRecoFiller.h" +#include "reco/PandoraLArRecoNDBranchFiller.h" #include "reco/SANDRecoBranchFiller.h" #include "truth/FillTruth.h" #include "beam/IFBeam.h" @@ -138,6 +139,14 @@ std::vector> getRecoFillers(const c std::cout << " SAND\n"; } + // Pandora LArRecoND + std::string pandoraFile; + if (par().cafmaker().pandoraLArRecoNDFile(pandoraFile)) + { + recoFillers.emplace_back(std::make_unique(pandoraFile)); + std::cout << " Pandora LArRecoND\n"; + } + // next: did we do TMS reco? std::string tmsFile; if (par().cafmaker().tmsRecoFile(tmsFile)) @@ -165,7 +174,7 @@ std::vector> getRecoFillers(const c { recoFillers.emplace_back(std::make_unique(par().cafmaker().trackMatchExtrapolatedZ(), par().cafmaker().trackMatchdX(), par().cafmaker().trackMatchdY(), par().cafmaker().trackMatchdThetaX(), par().cafmaker().trackMatchdThetaY())); std::cout << " ND-LAr + MINERvA matching\n"; - } + } // for now all the fillers get the same threshold. // if we decide we need to do it differently later // we can adjust the FCL params... @@ -375,7 +384,7 @@ void loop(CAF &caf, bool useIFBeam = false; if (ghepFilenames.empty() && edepsimFilename.empty() && !par().cafmaker().ForceDisableIFBeam()) useIFBeam = true; - + cafmaker::IFBeam beamManager(groupedTriggers, useIFBeam); //initialize IFBeam manager if data and when IFBeam is not force disabled // Main event loop @@ -407,17 +416,17 @@ void loop(CAF &caf, filler->FillRecoBranches(groupedTriggers[ii][0].second, caf.sr, par, &truthMatcher); } } - + //Fill POT double pot = 0.0; if (useIFBeam) { - pot = beamManager.getPOT(par, groupedTriggers[ii], ii); + pot = beamManager.getPOT(par, groupedTriggers[ii], ii); } else { pot = par().runInfo().POTPerSpill() * 1e13; - caf.sr.beam.ismc = true; + caf.sr.beam.ismc = true; } if (std::isnan(caf.pot)) caf.pot = 0; diff --git a/src/reco/PandoraLArRecoNDBranchFiller.cxx b/src/reco/PandoraLArRecoNDBranchFiller.cxx index 4e4a206..a53a5b4 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.cxx +++ b/src/reco/PandoraLArRecoNDBranchFiller.cxx @@ -49,9 +49,9 @@ namespace cafmaker m_LArRecoNDTree->SetBranchAddress("dirZ", &m_dirZVect); m_LArRecoNDTree->SetBranchAddress("energy", &m_energyVect); m_LArRecoNDTree->SetBranchAddress("n3DHits", &m_n3DHitsVect); - m_LArRecoNDTree->SetBranchAddress("mcVertexId", &m_mcVertexIdVect); + m_LArRecoNDTree->SetBranchAddress("mcNuId", &m_mcNuIdVect); + m_LArRecoNDTree->SetBranchAddress("mcLocalId", &m_mcLocalIdVect); m_LArRecoNDTree->SetBranchAddress("isPrimary", &m_isPrimaryVect); - m_LArRecoNDTree->SetBranchAddress("mcId", &m_mcIdVect); m_LArRecoNDTree->SetBranchAddress("completeness", &m_completenessVect); // We have setup the input tree @@ -89,7 +89,8 @@ namespace cafmaker sr.meta.nd_lar.run = m_run; sr.meta.nd_lar.subrun = m_subRun; - // Set the size for the Pandora NDLAr standard record for this trigger/event + // Set the size for the Pandora NDLAr standard record for this trigger/event. + // The number of clusters will be equal to the size of the sliceID vector const int nClusters = (m_sliceIdVect != nullptr) ? m_sliceIdVect->size() : 0; sr.nd.lar.pandora.resize(nClusters); sr.nd.lar.npandora = sr.nd.lar.pandora.size(); @@ -157,23 +158,21 @@ namespace cafmaker track.len_gcm2 = track.len_cm*m_LArRho; // Use truth matching info from Pandora's LArContent hierarchy tools. - // For LArRecoND MC SpacePoints, we offset the MCId's to make them all unique: - // 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 - // The vertex_id's are the same as those in the HDF5 files - - const long vertex_id = (m_mcVertexIdVect != nullptr) ? (*m_mcVertexIdVect)[i] : 0; + // For LArRecoND MC SpacePoints, we use the unique file_traj_id's + // for the MCParticles and the neutrino ID is the unique vertex_id. + // We also store the local traj_id's for the CAF truth matching tools + + // Neutrino ID = vertex_id + const long mcNuId = (m_mcNuIdVect != nullptr) ? (*m_mcNuIdVect)[i] : 0; + // MC particle ID = traj_id + const long mcId = (m_mcLocalIdVect != nullptr) ? (*m_mcLocalIdVect)[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 truePartID; if (isPrimary == 1) { truePartID.type = caf::TrueParticleID::kPrimary; - truePartID.part = traj_id; + truePartID.part = mcId; } else if (isPrimary == -1) { truePartID.type = caf::TrueParticleID::kUnknown; } else { @@ -183,7 +182,7 @@ namespace cafmaker try { // Get the true interaction in the stack - caf::SRTrueInteraction &srTrueInt = truthMatch->GetTrueInteraction(sr, vertex_id); + caf::SRTrueInteraction &srTrueInt = truthMatch->GetTrueInteraction(sr, mcNuId); 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)); @@ -191,7 +190,7 @@ namespace cafmaker // 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; }; + const auto pred = [&mcId](const caf::SRTrueParticle& part) { return part.G4ID == mcId; }; truePartID.part = std::distance(srTrueInt.sec.begin(), std::find_if(srTrueInt.sec.begin(), srTrueInt.sec.end(), pred)); } @@ -254,23 +253,21 @@ namespace cafmaker shower.Evis = energy; // Use truth matching info from Pandora's LArContent hierarchy tools. - // For LArRecoND MC SpacePoints, we offset the MCId's to make them all unique: - // 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 - // The vertex_id's are the same as those in the HDF5 files - - const long vertex_id = (m_mcVertexIdVect != nullptr) ? (*m_mcVertexIdVect)[i] : 0; + // For LArRecoND MC SpacePoints, we use the unique file_traj_id's + // for the MCParticles and the neutrino ID is the unique vertex_id. + // We also store the local traj_id's for the CAF truth matching tools + + // Neutrino ID = vertex_id + const long mcNuId = (m_mcNuIdVect != nullptr) ? (*m_mcNuIdVect)[i] : 0; + // MC particle ID = traj_id + const long mcId = (m_mcLocalIdVect != nullptr) ? (*m_mcLocalIdVect)[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 truePartID; if (isPrimary == 1) { truePartID.type = caf::TrueParticleID::kPrimary; - truePartID.part = traj_id; + truePartID.part = mcId; } else if (isPrimary == -1) { truePartID.type = caf::TrueParticleID::kUnknown; } else { @@ -280,7 +277,7 @@ namespace cafmaker try { // Get the true interaction in the stack - caf::SRTrueInteraction &srTrueInt = truthMatch->GetTrueInteraction(sr, vertex_id); + caf::SRTrueInteraction &srTrueInt = truthMatch->GetTrueInteraction(sr, mcNuId); 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)); @@ -288,7 +285,7 @@ namespace cafmaker // 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; }; + const auto pred = [&mcId](const caf::SRTrueParticle& part) { return part.G4ID == mcId; }; truePartID.part = std::distance(srTrueInt.sec.begin(), std::find_if(srTrueInt.sec.begin(), srTrueInt.sec.end(), pred)); } diff --git a/src/reco/PandoraLArRecoNDBranchFiller.h b/src/reco/PandoraLArRecoNDBranchFiller.h index 0f77b16..8d279a8 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.h +++ b/src/reco/PandoraLArRecoNDBranchFiller.h @@ -64,9 +64,9 @@ namespace cafmaker std::vector *m_dirZVect = nullptr; std::vector *m_energyVect = nullptr; std::vector *m_n3DHitsVect = nullptr; - std::vector *m_mcVertexIdVect = nullptr; + std::vector *m_mcNuIdVect = nullptr; + std::vector *m_mcLocalIdVect = nullptr; std::vector *m_isPrimaryVect = nullptr; - std::vector *m_mcIdVect = nullptr; std::vector *m_completenessVect = nullptr; float m_LArRho{1.3973f}; // LAr density (g/cm3) From 4048dcebd97c621fc6137140924ae561565ee122 Mon Sep 17 00:00:00 2001 From: John Back Date: Thu, 31 Oct 2024 11:42:48 +0000 Subject: [PATCH 09/10] Check if mcNuId != 0 for Pandora truth info. --- src/reco/PandoraLArRecoNDBranchFiller.cxx | 12 ++---------- 1 file changed, 2 insertions(+), 10 deletions(-) diff --git a/src/reco/PandoraLArRecoNDBranchFiller.cxx b/src/reco/PandoraLArRecoNDBranchFiller.cxx index a53a5b4..bfae72d 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.cxx +++ b/src/reco/PandoraLArRecoNDBranchFiller.cxx @@ -179,7 +179,7 @@ namespace cafmaker truePartID.type = caf::TrueParticleID::kSecondary; } - try + if (mcNuId != 0) { // Get the true interaction in the stack caf::SRTrueInteraction &srTrueInt = truthMatch->GetTrueInteraction(sr, mcNuId); @@ -195,10 +195,6 @@ namespace cafmaker std::find_if(srTrueInt.sec.begin(), srTrueInt.sec.end(), pred)); } } - catch (...) - { - LOG.VERBOSE() << "Could not fill some of the truth information" << "\n"; - } // Just store the best MC match std::vector truePartIDVect; @@ -274,7 +270,7 @@ namespace cafmaker truePartID.type = caf::TrueParticleID::kSecondary; } - try + if (mcNuId != 0) { // Get the true interaction in the stack caf::SRTrueInteraction &srTrueInt = truthMatch->GetTrueInteraction(sr, mcNuId); @@ -290,10 +286,6 @@ namespace cafmaker std::find_if(srTrueInt.sec.begin(), srTrueInt.sec.end(), pred)); } } - catch (...) - { - LOG.VERBOSE() << "Could not fill some of the truth information" << "\n"; - } // Just store the best MC match std::vector truePartIDVect; From 4eb842bdd86ffa4b2956593c2c11a3defc1c96c1 Mon Sep 17 00:00:00 2001 From: John Back Date: Fri, 8 Nov 2024 16:37:04 +0000 Subject: [PATCH 10/10] Address Jeremy's comments for Pandora's PR76: Unique pointers for the ROOT file/tree & LAr density pararmeter --- ndcaf_setup.sh | 1 - src/Params.h | 2 +- src/makeCAF.C | 3 ++- src/reco/PandoraLArRecoNDBranchFiller.cxx | 32 ++++++++++------------- src/reco/PandoraLArRecoNDBranchFiller.h | 21 +++++---------- 5 files changed, 24 insertions(+), 35 deletions(-) diff --git a/ndcaf_setup.sh b/ndcaf_setup.sh index 132acbf..eaee523 100644 --- a/ndcaf_setup.sh +++ b/ndcaf_setup.sh @@ -14,7 +14,6 @@ setup duneanaobj v03_06_01b -q e20:prof setup hdf5 v1_10_5a -q e20 setup fhiclcpp v4_15_03 -q debug:e20 setup edepsim v3_2_0c -q debug:e20 -setup root v6_26_06b -q e20:p3913:prof export LD_LIBRARY_PATH=$CURL_ROOT/lib:$LD_LIBRARY_PATH diff --git a/src/Params.h b/src/Params.h index 3788fca..9195e17 100644 --- a/src/Params.h +++ b/src/Params.h @@ -45,7 +45,6 @@ namespace cafmaker fhicl::OptionalAtom sandRecoFile { fhicl::Name{"SANDRecoFile"}, fhicl::Comment("Input SAND reco .root file") }; fhicl::OptionalAtom minervaRecoFile { fhicl::Name{"MINERVARecoFile"}, fhicl::Comment("Input MINERVA reco .root file") }; fhicl::OptionalAtom pandoraLArRecoNDFile { fhicl::Name{"PandoraLArRecoNDFile"}, fhicl::Comment("Input Pandora LArRecoND .root file") }; - // this is optional by way of the default value. Will result in an extra output file if enabled fhicl::Atom makeFlatCAF { fhicl::Name{"MakeFlatCAF"}, fhicl::Comment("Make 'flat' CAF in addition to structured CAF?"), true }; @@ -91,6 +90,7 @@ namespace cafmaker fhicl::Atom gastpc_B { fhicl::Name("gastpc_B"), fhicl::Comment("Gas TPC B field strength (Tesla)"), 0.4 }; fhicl::Atom gastpc_padPitch { fhicl::Name("gastpc_padPitch"), fhicl::Comment("(Fixed) pad pitch of gas TPC (cm)"), 0.1 }; // Actual pad pitch varies, which is going to be impossible to implement fhicl::Atom gastpc_X0 { fhicl::Name("gastpc_X0"), fhicl::Comment("Gas TPC radiation length (cm)"), 1300. }; + fhicl::Atom LArDensity { fhicl::Name("LArDensity"), fhicl::Comment("LAr density (g/cm3)"), 1.3973 }; }; /// FHICL table specifying which params are accepted diff --git a/src/makeCAF.C b/src/makeCAF.C index 68deaac..fef4dd8 100644 --- a/src/makeCAF.C +++ b/src/makeCAF.C @@ -143,7 +143,8 @@ std::vector> getRecoFillers(const c std::string pandoraFile; if (par().cafmaker().pandoraLArRecoNDFile(pandoraFile)) { - recoFillers.emplace_back(std::make_unique(pandoraFile)); + const float LArDensity = par().pseudoReco().LArDensity(); // For track lengths in g/cm2 + recoFillers.emplace_back(std::make_unique(pandoraFile, LArDensity)); std::cout << " Pandora LArRecoND\n"; } diff --git a/src/reco/PandoraLArRecoNDBranchFiller.cxx b/src/reco/PandoraLArRecoNDBranchFiller.cxx index bfae72d..adec501 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.cxx +++ b/src/reco/PandoraLArRecoNDBranchFiller.cxx @@ -1,30 +1,26 @@ #include "PandoraLArRecoNDBranchFiller.h" +#include "Params.h" + namespace cafmaker { - PandoraLArRecoNDBranchFiller::~PandoraLArRecoNDBranchFiller() - { - if (m_LArRecoNDFile && m_LArRecoNDFile->IsOpen()) - { - delete m_LArRecoNDTree; m_LArRecoNDTree = nullptr; - } - delete m_LArRecoNDFile; m_LArRecoNDFile = nullptr; - } - - PandoraLArRecoNDBranchFiller::PandoraLArRecoNDBranchFiller(const std::string &pandoraLArRecoNDFilename) + PandoraLArRecoNDBranchFiller::PandoraLArRecoNDBranchFiller(const std::string &pandoraLArRecoNDFilename, + const float LArDensity) : IRecoBranchFiller("PandoraLArRecoND"), m_Triggers(), - m_LastTriggerReqd(m_Triggers.end()) + m_LastTriggerReqd(m_Triggers.end()), + m_LArDensity(LArDensity) { // Open Pandora LArRecoND hierarchy analysis ROOT file - m_LArRecoNDFile = TFile::Open(pandoraLArRecoNDFilename.c_str(), "read"); + m_LArRecoNDFile.reset(TFile::Open(pandoraLArRecoNDFilename.c_str(), "read")); + if (m_LArRecoNDFile && m_LArRecoNDFile->IsOpen()) { LOG.VERBOSE() << " Using PandoraLArRecoND file " << pandoraLArRecoNDFilename << "\n"; // Input tree - m_LArRecoNDTree = dynamic_cast(m_LArRecoNDFile->Get("LArRecoND")); + m_LArRecoNDTree.reset(dynamic_cast(m_LArRecoNDFile->Get("LArRecoND"))); if (!m_LArRecoNDTree) { LOG.FATAL() << "Did not find LArRecoND tree in input file " << pandoraLArRecoNDFilename << "\n"; throw; @@ -102,8 +98,8 @@ namespace cafmaker } // ------------------------------------------------------------------------------ - void PandoraLArRecoNDBranchFiller::FillTracks(caf::StandardRecord &sr, const int nClusters, - const TruthMatcher *truthMatch) 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"; @@ -155,7 +151,7 @@ namespace cafmaker track.len_cm = sqrt(dX*dX + dY*dY + dZ*dZ); // Cluster length multiplied by LAr density (g/cm2) - track.len_gcm2 = track.len_cm*m_LArRho; + track.len_gcm2 = track.len_cm*m_LArDensity; // Use truth matching info from Pandora's LArContent hierarchy tools. // For LArRecoND MC SpacePoints, we use the unique file_traj_id's @@ -214,8 +210,8 @@ namespace cafmaker } // ------------------------------------------------------------------------------ - void PandoraLArRecoNDBranchFiller::FillShowers(caf::StandardRecord &sr, const int nClusters, - const TruthMatcher *truthMatch) 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"; diff --git a/src/reco/PandoraLArRecoNDBranchFiller.h b/src/reco/PandoraLArRecoNDBranchFiller.h index 8d279a8..e7d41dd 100644 --- a/src/reco/PandoraLArRecoNDBranchFiller.h +++ b/src/reco/PandoraLArRecoNDBranchFiller.h @@ -6,9 +6,7 @@ #ifndef ND_CAFMAKER_PandoraLArRecoNDBranchFiller_H #define ND_CAFMAKER_PandoraLArRecoNDBranchFiller_H -#include #include -#include // The virtual base class #include "reco/IRecoBranchFiller.h" @@ -26,25 +24,24 @@ namespace cafmaker class PandoraLArRecoNDBranchFiller : public cafmaker::IRecoBranchFiller { public: - PandoraLArRecoNDBranchFiller(const std::string & pandoraLArRecoNDFilename); + PandoraLArRecoNDBranchFiller(const std::string &pandoraLArRecoNDFilename, + const float LArDensity = 1.3973); - std::deque GetTriggers(int triggerType) const override; + std::deque GetTriggers(int triggerType) const override; RecoFillerType FillerType() const override { return RecoFillerType::BaseReco; } - ~PandoraLArRecoNDBranchFiller(); - private: void _FillRecoBranches(const Trigger &trigger, caf::StandardRecord &sr, const cafmaker::Params &par, - const TruthMatcher *truthMatch= nullptr) const override; + const TruthMatcher *truthMatch = nullptr) const override; 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; + std::unique_ptr m_LArRecoNDFile; + std::unique_ptr m_LArRecoNDTree; int m_eventId; int m_run; @@ -69,13 +66,9 @@ namespace cafmaker std::vector *m_isPrimaryVect = nullptr; std::vector *m_completenessVect = nullptr; - float m_LArRho{1.3973f}; // LAr density (g/cm3) - - int m_nuIdOffset{100000000}; - int m_maxMCId{1000000}; - mutable std::vector m_Triggers; mutable decltype(m_Triggers)::const_iterator m_LastTriggerReqd; ///< the last trigger requested using _FillRecoBranches + const float m_LArDensity; }; }