From 9a259be36111ce3aa5cf380cbfe8ae9012c7b14a Mon Sep 17 00:00:00 2001 From: dhavvval Date: Thu, 9 Apr 2026 10:30:59 -0500 Subject: [PATCH 1/6] It includes ToolAnalysis related changes in one single commit for adding Neutron Acestor and DirectParentID related changes --- DataModel/Hit.h | 201 +++++++++++------- DataModel/LAPPDHit.h | 26 ++- DataModel/Particle.h | 16 +- Makefile | 7 +- .../ANNIEEventTreeMaker.cpp | 110 ++++++++++ .../ANNIEEventTreeMaker/ANNIEEventTreeMaker.h | 10 + UserTools/BackTracker/BackTracker.cpp | 143 ++++++++++++- UserTools/BackTracker/BackTracker.h | 20 ++ UserTools/ClusterFinder/ClusterFinder.cpp | 3 +- UserTools/LoadWCSim/LoadWCSim.cpp | 118 ++++++---- UserTools/LoadWCSim/LoadWCSim.h | 4 +- UserTools/LoadWCSimLAPPD/LoadWCSimLAPPD.cpp | 9 +- UserTools/PMTWaveformSim/PMTWaveformSim.cpp | 76 +++++-- UserTools/PMTWaveformSim/PMTWaveformSim.h | 7 +- .../ANNIEEventTreeMakerConfig | 17 +- .../BeamClusterAnalysisMC/BackTrackerConfig | 2 + .../LoadGenieEventConfig | 12 +- .../BeamClusterAnalysisMC/LoadWCSimConfig | 5 +- .../LoadWCSimLAPPDConfig | 4 +- .../PMTWaveformSimConfig | 8 + .../PhaseIIADCHitFinderConfig | 11 + configfiles/BeamClusterAnalysisMC/ToolsConfig | 13 +- configfiles/LoadWCSim/LoadWCSimConfig | 5 +- configfiles/LoadWCSim/LoadWCSimLAPPDConfig | 3 +- configfiles/LoadWCSim/ToolsConfig | 4 +- configfiles/PMTWaveformSim/BackTrackerConfig | 2 + .../PMTWaveformSim/ClusterFinderConfig | 2 +- configfiles/PMTWaveformSim/LoadWCSimConfig | 4 +- configfiles/PMTWaveformSim/ToolsConfig | 3 +- 29 files changed, 656 insertions(+), 189 deletions(-) create mode 100644 configfiles/BeamClusterAnalysisMC/BackTrackerConfig create mode 100644 configfiles/BeamClusterAnalysisMC/PMTWaveformSimConfig create mode 100644 configfiles/BeamClusterAnalysisMC/PhaseIIADCHitFinderConfig create mode 100644 configfiles/PMTWaveformSim/BackTrackerConfig diff --git a/DataModel/Hit.h b/DataModel/Hit.h index ba78f611c..9ef6059c7 100755 --- a/DataModel/Hit.h +++ b/DataModel/Hit.h @@ -6,105 +6,150 @@ #include -using namespace std; +class Hit : public SerialisableObject { + + friend class boost::serialization::access; + +public: + Hit() + : TubeId(0) + , Time(0) + , Charge(0) + { + serialise=true; + } -class Hit : public SerialisableObject{ - - friend class boost::serialization::access; - - public: - Hit() : TubeId(0), Time(0), Charge(0){serialise=true;} - Hit(int thetubeid, double thetime, double thecharge) : TubeId(thetubeid), Time(thetime), Charge(thecharge){serialise=true;} - virtual ~Hit(){}; + Hit(int thetubeid, double thetime, double thecharge) + : TubeId(thetubeid) + , Time(thetime) + , Charge(thecharge) + { + serialise=true; + } + + virtual ~Hit(){}; - inline int GetTubeId() const {return TubeId;} - inline double GetTime() const {return Time;} - inline double GetCharge() const {return Charge;} + inline int GetTubeId() const {return TubeId;} + inline double GetTime() const {return Time;} + inline double GetCharge() const {return Charge;} - inline void SetTubeId(int tubeid){TubeId=tubeid;} - inline void SetTime(double tc){Time=tc;} - inline void SetCharge(double chg){Charge=chg;} + inline void SetTubeId(int tubeid){TubeId=tubeid;} + inline void SetTime(double tc){Time=tc;} + inline void SetCharge(double chg){Charge=chg;} - bool Print() { - std::cout<<"TubeId : "< void serialize(Archive & ar, const unsigned int version){ - if(serialise){ - ar & TubeId; - ar & Time; - ar & Charge; - } + template void serialize(Archive & ar, const unsigned int version) + { + if (serialise) { + ar & TubeId; + ar & Time; + ar & Charge; } + } }; // Derived classes class MCHit : public Hit { - // XXX ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ XXX - // XXX ~~~~~~~~~~~~~~~~~~~~~~~~ UPDATING THIS CLASS? ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ XXX - // XXX ~~~~~ Everything added in this class must be duplicated in MCLAPPDHit!~~~~ XXX - // XXX ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ XXX + // XXX ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ XXX + // XXX ~~~~~~~~~~~~~~~~~~~~~~~~ UPDATING THIS CLASS? ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ XXX + // XXX ~~~~~ Everything added in this class must be duplicated in MCLAPPDHit!~~~~ XXX + // XXX ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ XXX - friend class boost::serialization::access; + + friend class boost::serialization::access; - public: - MCHit() : Hit(), Parents(std::vector{}) {serialise=true;} - MCHit(int tubeid, double thetime, double thecharge, std::vector theparents) : Hit(tubeid, thetime, thecharge), Parents(theparents) {serialise=true;} - virtual ~MCHit(){}; +public: + MCHit() + : Hit() + , Parents(std::vector{}) + , DirectParents(std::vector{}) + , StartTick(-5) + , EndTick(-5) + { + serialise=true; + } + + // Start and End ticks are only ever to be set after initialization + MCHit(int tubeid, double thetime, double thecharge, std::vector theparents, std::vector thedirectparents) + : Hit(tubeid, thetime, thecharge) + , Parents(theparents) + , DirectParents(thedirectparents) + , StartTick(-5) + , EndTick(-5) + { + serialise=true; + } + + virtual ~MCHit(){}; - const std::vector* GetParents() const { return &Parents; } - void SetParents(std::vector parentsin){ Parents = parentsin; } + const std::vector* GetParents() { return &Parents; } + const std::vector* GetDirectParents() { return &DirectParents; } + int GetStartTick() { return StartTick; } + int GetEndTick() { return EndTick; } + + void SetParents(std::vector parentsin) { Parents = parentsin; } + void SetDirectParents(std::vector directparentsin) { DirectParents = directparentsin; } + void SetStartTick(int tick) { StartTick = tick; } + void SetEndTick(int tick) { EndTick = tick; } + - bool Print(){ - std::cout<<"TubeId : "< Parents; + return true; + } - template void serialize(Archive & ar, const unsigned int version){ - if(serialise){ - ar & TubeId; - ar & Time; - ar & Charge; - // do not serialize parents; the indices by themselves are not meaningful - } +protected: + std::vector Parents; + std::vector DirectParents; + int StartTick; + int EndTick; + + template void serialize(Archive & ar, const unsigned int version) + { + if (serialise) { + ar & TubeId; + ar & Time; + ar & Charge; + + if (version > 0) { + ar & Parents; // Parents is now track IDs rather than index within vector + ar & DirectParents; + ar & StartTick; + ar & EndTick; + } } + } }; -/* -class RecoHit : public Hit { - public: - RecoHit(double thetime, double thecharge) : Time(thetime), Charge(thecharge){}; - - inline double GetCharge(){return Charge;} - inline void SetCharge(double chg){Charge=chg;} - - protected: - double Charge; -}; -*/ +BOOST_CLASS_VERSION(MCHit, 1) #endif diff --git a/DataModel/LAPPDHit.h b/DataModel/LAPPDHit.h index 251ac38b8..d84604f37 100755 --- a/DataModel/LAPPDHit.h +++ b/DataModel/LAPPDHit.h @@ -100,12 +100,15 @@ class MCLAPPDHit : public LAPPDHit friend class boost::serialization::access; public: - MCLAPPDHit() : LAPPDHit(), Parents(std::vector{}) { serialise = true; } - MCLAPPDHit(int thetubeid, double thetime, double thecharge, std::vector theposition, std::vector thelocalposition, std::vector theparents) : LAPPDHit(thetubeid, thetime, thecharge, theposition, thelocalposition), Parents(theparents) { serialise = true; } + MCLAPPDHit() : LAPPDHit(), Parents(std::vector{}), DirectParents(std::vector{}) { serialise = true; } + MCLAPPDHit(int thetubeid, double thetime, double thecharge, std::vector theposition, std::vector thelocalposition, std::vector theparents, std::vector thedirectparents) : LAPPDHit(thetubeid, thetime, thecharge, theposition, thelocalposition), Parents(theparents), DirectParents(thedirectparents) { serialise = true; } const std::vector *GetParents() const { return &Parents; } void SetParents(std::vector parentsin) { Parents = parentsin; } + const std::vector *GetDirectParents() const { return &DirectParents; } + void SetDirectParents(std::vector directparentsin) { DirectParents = directparentsin; } + bool Print() { cout << "TubeId : " << TubeId << endl; @@ -116,6 +119,7 @@ class MCLAPPDHit : public LAPPDHit cout << "Parallel Pos : " << LocalPosition.at(0) << endl; cout << "Transverse Pos : " << LocalPosition.at(1) << endl; cout << "Charge : " << Charge << endl; + if (Parents.size()) { cout << "Parent MCPartice indices: {"; @@ -131,6 +135,23 @@ class MCLAPPDHit : public LAPPDHit { cout << "No recorded parents" << endl; } + + if (DirectParents.size()) + { + cout << "Direct Parent MCPartice indices: {"; + for (int parenti = 0; parenti < (int)DirectParents.size(); ++parenti) + { + cout << DirectParents.at(parenti); + if ((parenti + 1) < (int)DirectParents.size()) + cout << ", "; + } + cout << "}" << endl; + } + else + { + cout << "##### No recorded Direct parents #####" << endl; + } + return true; } @@ -151,6 +172,7 @@ class MCLAPPDHit : public LAPPDHit protected: std::vector Parents; + std::vector DirectParents; }; /* diff --git a/DataModel/Particle.h b/DataModel/Particle.h index 2917e3e63..1df3d83bc 100644 --- a/DataModel/Particle.h +++ b/DataModel/Particle.h @@ -137,13 +137,13 @@ class MCParticle : public Particle { public: MCParticle() : Particle(0, 0., 0., Position(), Position(), 0., 0., Direction(), 0., - tracktype::UNCONTAINED), ParticleID(0), ParentPdg(0), StartsInFiducialVolume(false), TrackAngleX(0), TrackAngleY(0), TrackAngleFromBeam(0), EntersTank(false), TankEntryPoint(Position()), ExitsTank(false), TankExitPoint(Position()), TrackLengthInTank(0), EntersMrd(false), MrdEntryPoint(Position()), ExitsMrd(false), MrdExitPoint(Position()), PenetratesMrd(false), TrackLengthInMrd(0), MrdPenetration(0), MrdLayersPenetrated(0), MrdEnergyLoss(0), Flag(0), MCTriggerNum(0) {serialise=true;} + tracktype::UNCONTAINED), ParticleID(0), ParentPdg(0), PrimaryParentID(0), DirectParentID(0), StartsInFiducialVolume(false), TrackAngleX(0), TrackAngleY(0), TrackAngleFromBeam(0), EntersTank(false), TankEntryPoint(Position()), ExitsTank(false), TankExitPoint(Position()), TrackLengthInTank(0), EntersMrd(false), MrdEntryPoint(Position()), ExitsMrd(false), MrdExitPoint(Position()), PenetratesMrd(false), TrackLengthInMrd(0), MrdPenetration(0), MrdLayersPenetrated(0), MrdEnergyLoss(0), Flag(0), MCTriggerNum(0) {serialise=true;} MCParticle(int pdg, double sttE, double stpE, Position sttpos, Position stppos, double sttt, double stpt, Direction startdir, double len, tracktype tracktypein, - int partid, int parentpdg, int flagid, int triggernum) + int partid, int parentpdg, int primaryparentid, int directparentid, int flagid, int triggernum) : Particle(pdg, sttE, stpE, sttpos, stppos, sttt, stpt, startdir, len, tracktypein), - ParticleID(partid), ParentPdg(parentpdg), StartsInFiducialVolume(false), TrackAngleX(0), TrackAngleY(0), TrackAngleFromBeam(0), EntersTank(false), TankEntryPoint(Position()), ExitsTank(false), TankExitPoint(Position()), TrackLengthInTank(0), EntersMrd(false), MrdEntryPoint(Position()), ExitsMrd(false), MrdExitPoint(Position()), PenetratesMrd(false), TrackLengthInMrd(0), MrdPenetration(0), MrdLayersPenetrated(0), MrdEnergyLoss(0), Flag(flagid), MCTriggerNum(triggernum) + ParticleID(partid), ParentPdg(parentpdg), PrimaryParentID(primaryparentid), DirectParentID(directparentid), StartsInFiducialVolume(false), TrackAngleX(0), TrackAngleY(0), TrackAngleFromBeam(0), EntersTank(false), TankEntryPoint(Position()), ExitsTank(false), TankExitPoint(Position()), TrackLengthInTank(0), EntersMrd(false), MrdEntryPoint(Position()), ExitsMrd(false), MrdExitPoint(Position()), PenetratesMrd(false), TrackLengthInMrd(0), MrdPenetration(0), MrdLayersPenetrated(0), MrdEnergyLoss(0), Flag(flagid), MCTriggerNum(triggernum) { serialise=true; // override Hit tracktype @@ -160,6 +160,8 @@ class MCParticle : public Particle { inline int GetParticleID(){return ParticleID;} inline int GetParentPdg(){return ParentPdg;} + inline int GetPrimaryParentID(){return PrimaryParentID;} + inline int GetDirectParentID(){return DirectParentID;} inline int GetFlag(){return Flag;} inline int GetMCTriggerNum(){return MCTriggerNum;} @@ -188,6 +190,8 @@ class MCParticle : public Particle { inline void SetParticleID(int partidin){ParticleID=partidin;} inline void SetParentPdg(int parentpdgin){ParentPdg=parentpdgin;} + inline void SetPrimaryParentID(int primaryparentidin){PrimaryParentID=primaryparentidin;} + inline void SetDirectParentID(int directparentidin){DirectParentID=directparentidin;} inline void SetFlag(int flagidin){Flag=flagidin;} inline void SetMCTriggerNum(int triggernumin){MCTriggerNum=triggernumin;} @@ -228,6 +232,7 @@ class MCParticle : public Particle { bool Print() { std::cout<<"ParticlePDG : "<Branch("hitPMTType", &fHitPMTType); } + if (DirectParent_MCHit_fill){ + fANNIETree->Branch("DirectParent_PMTID", &fDirectParent_PMTID); + fANNIETree->Branch("DirectParent_HitTime", &fDirectParent_HitTime); + fANNIETree->Branch("DirectParent_TrackIDs", &fDirectParent_TrackIDs); + fANNIETree->Branch("DirectParent_PDGs", &fDirectParent_PDGs); + fANNIETree->Branch("DirectParent_NeutronAncestorTrackID", &fDirectParent_NeutronAncestorTrackID); + fANNIETree->Branch("DirectParent_NeutronAncestorPDG", &fDirectParent_NeutronAncestorPDG); + } + if (SiPMPulseInfo_fill) { fANNIETree->Branch("SiPMhitQ", &fSiPMHitQ); @@ -596,6 +606,12 @@ bool ANNIEEventTreeMaker::Execute() // this will fill all hits in this event LoadAllTankHits(); } + + //****************************** Fill MCHit DirectParent TrackIDs Info *************************************// + if (DirectParent_MCHit_fill) + { + LoadDirectParentIDsMCHits(); + } if (SiPMPulseInfo_fill) { LoadSiPMHits(); @@ -764,6 +780,14 @@ void ANNIEEventTreeMaker::ResetVariables() fHitChankeyMC.clear(); fHitPMTType.clear(); + // MCHit DirectParent TrackIDs info + fDirectParent_PMTID.clear(); + fDirectParent_HitTime.clear(); + fDirectParent_TrackIDs.clear(); + fDirectParent_PDGs.clear(); + fDirectParent_NeutronAncestorTrackID.clear(); + fDirectParent_NeutronAncestorPDG.clear(); + // SiPMPulse Info fSiPM1NPulses = 0; fSiPM2NPulses = 0; @@ -1384,6 +1408,92 @@ void ANNIEEventTreeMaker::LoadAllTankHits() return; } +// **************MCHit Directparent TrackIDs Info ************************** // + +void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ + //I will make changes here// + //It will help me to store the information about the direct parent track IDs for each MCHit in the tree + Log("ANNIEEventTreeMaker Tool: LoadDirectParentIDsMCHits", v_debug, ANNIEEventTreeMakerVerbosity); + std::map>> *fMCHitToDirectParents = nullptr; + std::vector *fMCParticles = nullptr; + std::map *fTrackIdToIndex = nullptr; + std::map>> *fMCHitToNeutronAncestor = nullptr; + + bool got_MCHitToDirectParents = m_data->Stores["ANNIEEvent"]->Get("MCHitToDirectParents", fMCHitToDirectParents); + if (!got_MCHitToDirectParents) { + std::cout << "No MCHitToDirectParents store in ANNIEEvent. Continuing to build tree " << std::endl; + return; + } + + bool got_MCParticles = m_data->Stores["ANNIEEvent"]->Get("MCParticles", fMCParticles); + if (!got_MCParticles) { + std::cout << "No MCParticles store in ANNIEEvent. Continuing to build tree " << std::endl; + return; + } + + bool got_TrackIdToIndex = m_data->Stores["ANNIEEvent"]->Get("TrackId_to_MCParticleIndex", fTrackIdToIndex); + if (!got_TrackIdToIndex) { + std::cout << "No TrackId_to_MCParticleIndex store in ANNIEEvent. Continuing to build tree " << std::endl; + return; + } + + bool got_neutronAncestor = m_data->Stores["ANNIEEvent"]->Get("MCHitToNeutronAncestor", fMCHitToNeutronAncestor); + if (!got_neutronAncestor) { + std::cout << "No MCHitToNeutronAncestor store in ANNIEEvent. Continuing to build tree " << std::endl; + return; + } + + for (auto const& apair : *fMCHitToDirectParents) { + unsigned long pmtID = apair.first; + for (auto const& hit_directparent_pair : apair.second){ + double hitTime = hit_directparent_pair.first; + std::vector const& directparentids = hit_directparent_pair.second; + + fDirectParent_PMTID.push_back(pmtID); + fDirectParent_HitTime.push_back(hitTime); + fDirectParent_TrackIDs.push_back(directparentids); + + std::vector pdgcodes; + if (got_MCParticles && got_TrackIdToIndex){ + for (int directparentid : directparentids){ + auto it = fTrackIdToIndex->find(directparentid); + if (it != fTrackIdToIndex->end()) { + int MCParticleIndex = it->second; + int pdg = fMCParticles->at(MCParticleIndex).GetPdgCode();; + + std::cout << "DEBUG DirectParent | " + << "TrackID(from hit)=" << directparentid + << ", MCParticle.GetPdgCode()=" << pdg + << std::endl; + + pdgcodes.push_back(pdg); + } + else { + std::cout << "DEBUG DirectParent | TrackID=" << directparentid + << " NOT FOUND in MCParticles (will use -999)" << std::endl; + pdgcodes.push_back(-999); + } + } + } + fDirectParent_PDGs.push_back(pdgcodes); + + int neutronAncestorTrackID = -5; + int neutronAncestorPDG = -5; + if (got_neutronAncestor && fMCHitToNeutronAncestor->find(pmtID) != fMCHitToNeutronAncestor->end()){ + auto const& pmtAncestors = fMCHitToNeutronAncestor->at(pmtID); + if (pmtAncestors.find(hitTime) != pmtAncestors.end()){ + auto const& ancestorPair = pmtAncestors.at(hitTime); + neutronAncestorTrackID = ancestorPair.first; + neutronAncestorPDG = ancestorPair.second; + } + } + fDirectParent_NeutronAncestorTrackID.push_back(neutronAncestorTrackID); + fDirectParent_NeutronAncestorPDG.push_back(neutronAncestorPDG); + } + } + return; +} + void ANNIEEventTreeMaker::LoadSiPMHits() { Log("ANNIEEventTreeMaker Tool: LoadSiPMHits", v_debug, ANNIEEventTreeMakerVerbosity); diff --git a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h index 864ac83dc..daae6708a 100644 --- a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h +++ b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h @@ -56,6 +56,7 @@ class ANNIEEventTreeMaker : public Tool void LoadRWMBRFInfo(); void LoadAllTankHits(); + void LoadDirectParentIDsMCHits(); void LoadSiPMHits(); void LoadLAPPDInfo(); @@ -218,6 +219,15 @@ class ANNIEEventTreeMaker : public Tool std::vector fHitChankeyMC; std::vector fHitPMTType; + // DirectParent_MCHit_fill + bool DirectParent_MCHit_fill = 0; + std::vector fDirectParent_PMTID; + std::vector fDirectParent_HitTime; + std::vector> fDirectParent_TrackIDs; + std::vector> fDirectParent_PDGs; + std::vector fDirectParent_NeutronAncestorTrackID; //Stored unique neutron ancestor track ID for each MCHit, if it exists. -5 if no neutron ancestor found. + std::vector fDirectParent_NeutronAncestorPDG; //For now, we are storing Neutron's truth information. But, it could be exapand to other particle types if needed. -5 if no neutron ancestor found. + // SiPMPulseInfo_fill int fSiPM1NPulses; int fSiPM2NPulses; diff --git a/UserTools/BackTracker/BackTracker.cpp b/UserTools/BackTracker/BackTracker.cpp index ce80ada9b..d0d1b316e 100644 --- a/UserTools/BackTracker/BackTracker.cpp +++ b/UserTools/BackTracker/BackTracker.cpp @@ -1,4 +1,5 @@ #include "BackTracker.h" +#include "ANNIEconstants.h" BackTracker::BackTracker():Tool(){} @@ -27,6 +28,9 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){ Log(logmessage, v_error, verbosity); } + bool gotUsePulseWindowMatching = m_variables.Get("UseDirectParentClockTickMatching", fDirectParentClockTickMatching); + if (!gotUsePulseWindowMatching) fDirectParentClockTickMatching = true; + // Set up the pointers we're going to save. No need to // delete them at Finalize, the store will handle it @@ -35,11 +39,12 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){ fClusterEfficiency = new std::map; fClusterPurity = new std::map; fClusterTotalCharge = new std::map; + fMCHitToDirectParents = new std::map>>; + fMCHitToNeutronAncestor = new std::map>>; return true; } -//------------------------------------------------------------------------------ bool BackTracker::Execute() { if (!LoadFromStores()) @@ -50,11 +55,19 @@ bool BackTracker::Execute() fClusterEfficiency ->clear(); fClusterPurity ->clear(); fClusterTotalCharge ->clear(); + fMCHitToDirectParents ->clear(); + // fMCHitToNeutronAncestor ->clear(); fParticleToTankTotalCharge.clear(); + SumParticleTankCharge(); - + if (fDirectParentClockTickMatching) { + // Required tool order: PMTWaveformSim -> PhaseIIADCHitFinder -> BackTracker + DirectParentsFromClockTickWindows(); + FindNeutronAncestors(); + } + // Loop over the clusters and do the things for (std::pair>&& apair : *fClusterMapMC) { int prtId = -5; @@ -70,6 +83,7 @@ bool BackTracker::Execute() fClusterEfficiency ->emplace(apair.first, eff); fClusterPurity ->emplace(apair.first, pur); fClusterTotalCharge ->emplace(apair.first, totalCharge); + } m_data->Stores.at("ANNIEEvent")->Set("ClusterToBestParticleID", fClusterToBestParticleID ); @@ -77,6 +91,8 @@ bool BackTracker::Execute() m_data->Stores.at("ANNIEEvent")->Set("ClusterEfficiency", fClusterEfficiency ); m_data->Stores.at("ANNIEEvent")->Set("ClusterPurity", fClusterPurity ); m_data->Stores.at("ANNIEEvent")->Set("ClusterTotalCharge", fClusterTotalCharge ); + m_data->Stores.at("ANNIEEvent")->Set("MCHitToDirectParents", fMCHitToDirectParents ); + m_data->Stores.at("ANNIEEvent")->Set("MCHitToNeutronAncestor", fMCHitToNeutronAncestor ); return true; } @@ -174,7 +190,101 @@ void BackTracker::MatchMCParticle(std::vector const &mchits, int &prtId, } -//------------------------------------------------------------------------------ +void BackTracker::DirectParentsFromClockTickWindows() +{ + if (fPMTToDirectParentMap.empty() || fRecoADCHits.empty()) return; + + const double prewindow_ns = static_cast(fPMTSimPrewindowTicks) * NS_PER_ADC_SAMPLE; + const double readout_ns = static_cast(fPMTSimReadoutWindowTicks) * NS_PER_ADC_SAMPLE; + + for (auto const& recoIt : fRecoADCHits) { + unsigned long pmtID = recoIt.first; + + auto parentMapIt = fPMTToDirectParentMap.find(pmtID); + if (parentMapIt == fPMTToDirectParentMap.end()) continue; + + std::map> const& hits_to_directparents_map = parentMapIt->second; + + for (std::vector const& minibufPulses : recoIt.second) { + for (ADCPulse const& pulse : minibufPulses) { + double pulseStart = pulse.start_time(); + double hitTime = pulse.peak_time(); + + double tmin = pulseStart - prewindow_ns; + double tmax = pulseStart + readout_ns; + + for (auto const& apair : hits_to_directparents_map) { + double mchitTime = static_cast(apair.first) * NS_PER_ADC_SAMPLE; + + if (mchitTime > tmin && mchitTime < tmax) { + (*fMCHitToDirectParents)[pmtID][hitTime].insert( + (*fMCHitToDirectParents)[pmtID][hitTime].end(), + apair.second.begin(), apair.second.end()); + } + } + } + } + } +} + +void BackTracker::FindNeutronAncestors() { + + std::map> trackMap; // trackId -> (ParentID, pdg) + for (auto& particle : *fMCParticles) { + int trackId = particle.GetParticleID(); + int parentId = particle.GetDirectParentID(); + int pdg = particle.GetPdgCode(); + + trackMap[trackId] = std::make_pair(parentId, pdg); + } + + for (const auto& pmtPair : *fMCHitToDirectParents) { + unsigned long pmtID = pmtPair.first; + for (const auto& hitPair : pmtPair.second) { + double hitTime = hitPair.first; + const std::vector& directParents = hitPair.second; + + if (directParents.empty()) continue; + + int neutronAncestorId = -5; + int neutronAncestorPdg = -5; + int startParentID = directParents[0]; // take the first direct parent as the starting point + int currentID = startParentID; + + + std::set visited; //Ensures that Ancestory finding doesn't get stuck in a loop if there are any circular references in the MCParticles + + while (trackMap.find(currentID) != trackMap.end()) { + if (visited.count(currentID) > 0) { + std::cerr << "WARNING: Circular reference detected at trackID " << currentID << std::endl; + break; + } + visited.insert(currentID); + + int parentID = trackMap[currentID].first; + int pdg = trackMap[currentID].second; + + if (pdg == 2112 && neutronAncestorId == -5) { // if it's a neutron and we haven't already found an ancestor, save it + neutronAncestorId = currentID; // This is the neutron's TRACK ID + neutronAncestorPdg = pdg; // Store the PDG code (2112 for neutron) + break; //We care about the immidiate neutron ancestor. + //Do we want to find the potential primary neutron ancestor, if there is one? (DJA) + } + + if (parentID == -1) break; // reached the end of the ancestry + currentID = parentID; + } + + (*fMCHitToNeutronAncestor)[pmtID][hitTime] = std::make_pair(neutronAncestorId, neutronAncestorPdg); + + } + } + + std::cout << "BackTracker::FindNeutronAncestors: found " + << fMCHitToNeutronAncestor->size() + << " PMTs with neutron ancestors." << std::endl; +} + bool BackTracker::LoadFromStores() { // Grab the stuff we need from the stores @@ -208,6 +318,31 @@ bool BackTracker::LoadFromStores() return false; } + if (fDirectParentClockTickMatching) { + fPMTToDirectParentMap.clear(); + fRecoADCHits.clear(); + + bool gotDirectParentMap = m_data->Stores.at("ANNIEEvent")->Get("PMTToDirectParentMap", fPMTToDirectParentMap); + if (!gotDirectParentMap) { + logmessage = "BackTracker: PMTToDirectParentMap missing, disabling pulse-window matching for this event."; + Log(logmessage, v_warning, verbosity); + } + + bool gotRecoADCHits = m_data->Stores.at("ANNIEEvent")->Get("RecoADCHits", fRecoADCHits); + if (!gotRecoADCHits) { + logmessage = "BackTracker: RecoADCHits missing, disabling pulse-window matching for this event."; + Log(logmessage, v_warning, verbosity); + } + + uint16_t prewindowTicks = fPMTSimPrewindowTicks; + uint16_t readoutTicks = fPMTSimReadoutWindowTicks; + if (m_data->Stores.at("ANNIEEvent")->Get("PMTSimPrewindowTicks", prewindowTicks)) { + fPMTSimPrewindowTicks = prewindowTicks; + } + if (m_data->Stores.at("ANNIEEvent")->Get("PMTSimReadoutWindowTicks", readoutTicks)) { + fPMTSimReadoutWindowTicks = readoutTicks; + } + } + return true; } - diff --git a/UserTools/BackTracker/BackTracker.h b/UserTools/BackTracker/BackTracker.h index 1cf5c7240..266ade505 100644 --- a/UserTools/BackTracker/BackTracker.h +++ b/UserTools/BackTracker/BackTracker.h @@ -3,8 +3,10 @@ #include #include +#include #include "Tool.h" +#include "ADCPulse.h" #include "Hit.h" #include "Particle.h" @@ -31,6 +33,8 @@ class BackTracker: public Tool { bool LoadFromStores(); ///< Does all the loading so I can move it away from the Execute function void SumParticleTankCharge(); void MatchMCParticle(std::vector const &mchits, int &prtId, int &prtPdg, double &eff, double &pur, double &totalCharge); ///< The meat and potatoes + void DirectParentsFromClockTickWindows(); + void FindNeutronAncestors(); private: @@ -39,6 +43,8 @@ class BackTracker: public Tool { std::map> *fClusterMapMC = nullptr; ///< Clusters that we will be linking MCParticles to std::vector *fMCParticles = nullptr; ///< The true particles from the event std::map *fMCParticleIndexMap = nullptr; ///< Map between the particle Id and it's position in MCParticles vector + std::map>> fPMTToDirectParentMap; ///< PMT -> t0 tick -> direct parent IDs from PMTWaveformSim + std::map>> fRecoADCHits; ///< Reconstructed ADCPulses from PhaseIIADCHitFinder // We'll calculate this map from MCHit parent particle to the total charge deposited throughout the tank // technically a MCHit could have multiple parents, but they don't appear to in practice @@ -57,6 +63,20 @@ class BackTracker: public Tool { std::map *fClusterPurity = nullptr; std::map *fClusterTotalCharge = nullptr; + // Cluster Time -> MCHit DirectParentIDs + std::map>> *fClusterHitToDirectParentTrackIDs = nullptr; + // PMT ID -> reco hit time -> direct parent track IDs + std::map>> *fMCHitToDirectParents = nullptr; + + // PMT ID -> reco hit time -> (neutron trackID, neutron PDG) + std::map>> *fMCHitToNeutronAncestor = nullptr; + + bool fDirectParentClockTickMatching = true; + uint16_t fPMTSimPrewindowTicks = 10; + uint16_t fPMTSimReadoutWindowTicks = 35; + + + /// \brief verbosity levels: if 'verbosity' < this level, the message type will be logged. int verbosity; int v_error=0; diff --git a/UserTools/ClusterFinder/ClusterFinder.cpp b/UserTools/ClusterFinder/ClusterFinder.cpp index 6d41ff085..1b05243a2 100644 --- a/UserTools/ClusterFinder/ClusterFinder.cpp +++ b/UserTools/ClusterFinder/ClusterFinder.cpp @@ -302,10 +302,11 @@ bool ClusterFinder::Execute(){ v_hittimes.push_back(datalike_hits.at(i_hit)); } std::vector parents = *(ThisPMTHits.at(0).GetParents()); + std::vector directparents = *(ThisPMTHits.at(0).GetDirectParents()); ThisPMTHits.clear(); std::vector newMCHits; for (int i_hit=0; i_hit < (int) datalike_hits.size(); i_hit++){ - newMCHits.push_back(MCHit(chankey,datalike_hits.at(i_hit),datalike_hits_charge.at(i_hit),parents)); + newMCHits.push_back(MCHit(chankey,datalike_hits.at(i_hit),datalike_hits_charge.at(i_hit),parents, directparents)); } MCHits->at(chankey) = newMCHits; } diff --git a/UserTools/LoadWCSim/LoadWCSim.cpp b/UserTools/LoadWCSim/LoadWCSim.cpp index 16ff3636a..c34df9311 100644 --- a/UserTools/LoadWCSim/LoadWCSim.cpp +++ b/UserTools/LoadWCSim/LoadWCSim.cpp @@ -8,16 +8,20 @@ bool LoadWCSim::Initialise(std::string configfile, DataModel &data) if (configfile!="") m_variables.Initialise(configfile); //loading config file //m_variables.Print(); + /////////////////// Useful header /////////////////////// + + if (verbosity) cout << "Initializing Tool LoadWCSim" << endl; + + if (configfile!="") m_variables.Initialise(configfile); //loading config file + //m_variables.Print(); + m_data = &data; //assigning transient data pointer // Get the Tool configuration variables and set defaults // ====================================================== if (!m_variables.Get("verbose", verbosity)) verbosity = 1; - logmessage = "LoadWCSim::Initialise: Initialising LoadWCSim!"; - Log(logmessage, v_warning, verbosity); - if (!m_variables.Get("MaxEntries", MaxEntries)) MaxEntries = -1; - + if (!m_variables.Get("InputFile", MCFile)) { logmessage = "LoadWCSim::Initialise: NO InputFile set in the config!"; Log(logmessage, v_error, verbosity); @@ -198,22 +202,23 @@ bool LoadWCSim::Initialise(std::string configfile, DataModel &data) // Short Stores README // ====================================================== - // n.b. m_data->vars is a Store (of ben's Store type) that is not saved to disk. + // n.b. m_data->vars is a Store (of ben's Store type) that is not saved to disk? // m_data->CStore is a single entry binary BoostStore that is not saved to disk. // m_data->Stores["StoreName"] is a map of binary BoostStores that are saved to disk. - // With BoostStore::Set("MyVariable", myvar), if myvar is not a pointer it will always - // be saved to disk, but if myvar is a pointer the persist flag (default true) controls - // whether it will be saved to disk. In either case the BoostStore becomes the owner - // of the object and will handle its deletion. + // If using Stores->BoostStore->Set("MyVariable") it will always be saved to disk + // Using Stores->BoostStore.Set("MyVariable",&myvar) if myvar is a pointer (to an object on + // the heap) puts myvar in the Store and it's deletion will be handled by the Store. + // (provided your class has a suitable destructor.) // Is 'BoostStore::Save' needed for single-entry stores? // ---------------- + // create a new BoostStore with key "ANNIEEvent" in the Stores std::map // BoostStore constructor args: typechecking (bool), m_format (0=binary, 1=ASCII, 2=multievent) // A BoostStore has a header where useful constants may be saved. The header is a BoostStore itself, // and can be accessed via: Store.Header->Get() and Store.Header->Set(). - // The method 'BoostStore::Save()' writes everything 'Set' since the last 'Save' to the current entry. - // BoostStore::Clear clears the map of the current entry, to start building a new one. - // Use 'BoostStore::GetEntry(int entrynum)' to load an entry to then be able to 'Get' it's contents. - // 'BoostStore->Header->Get("TotalEntries",NumEvents)' will load the num entries into NumEvents + // The method 'Store::Save()' writes everything 'Set' since the last 'Save' to the current entry. + // Store::Clear clears the map of the current entry, to start building a new one. + // Use 'Store::GetEntry(int entrynum)' to load an entry to then be able to 'Get' it's contents. + // 'Store->Header->Get("TotalEntries",NumEvents)' will load the num entries into NumEvents // ------------------ // When adding a BoostStore (or class object in general) to a BoostStore (such as ANNIEEvent) // you call BoostStore::Set("key",ObjectPointer) - BUT be aware that serialization happens when the @@ -628,7 +633,7 @@ bool LoadWCSim::Execute() // If we merged the subtriggers, report an MCTriggerNum of 0 // so downstream tools know it's a new event int reported_triggernum = splitSubtriggers ? MCTriggerNum-1 : 0; - m_data->Stores.at("ANNIEEvent")->Set("MCTriggerNum", reported_triggernum); + m_data->Stores.at("ANNIEEvent")->Set("MCTriggernum", reported_triggernum); m_data->Stores.at("ANNIEEvent")->Set("MCFile", MCFile); m_data->Stores.at("ANNIEEvent")->Set("MCFlag", true); m_data->Stores.at("ANNIEEvent")->Set("BeamStatus", beamstat); @@ -1101,12 +1106,12 @@ void LoadWCSim::LoadMCParticles(WCSimRootTrigger* firstTrig) logmessage += " tracks from trigger # " + std::to_string(trigIdx); Log(logmessage, v_message, verbosity); + std::cout<< "Event Number: " << aTrigTank->GetHeader()->GetEvtNum()<< std::endl; for (int trackIdx = 0; trackIdx < aTrigTank->GetNtrack(); trackIdx++) { logmessage = "LoadWCSim::LoadMCParticles: Getting WCSim track # " + std::to_string(trackIdx); Log(logmessage, v_message, verbosity); auto* nextTrack = (WCSimRootTrack*)aTrigTank->GetTracks()->At(trackIdx); - tracktype startStopType = tracktype::UNDEFINED; // Extract the neutrino information @@ -1120,14 +1125,16 @@ void LoadWCSim::LoadMCParticles(WCSimRootTrigger* firstTrig) double length = (stopPos-startPos).Mag(); MCParticle neutrino(nextTrack->GetIpnu(), nextTrack->GetE(), nextTrack->GetEndE(), - startPos, stopPos, startTime, stopTime, - Direction(nextTrack->GetDir(0), nextTrack->GetDir(1), nextTrack->GetDir(2)), - length, startStopType, - nextTrack->GetId(), - nextTrack->GetParenttype(), - nextTrack->GetFlag(), - trigIdx); - + startPos, stopPos, startTime, stopTime, + Direction(nextTrack->GetDir(0), nextTrack->GetDir(1), nextTrack->GetDir(2)), + length, startStopType, + nextTrack->GetId(), + nextTrack->GetParenttype(), + nextTrack->GetPrimaryParentID(), + nextTrack->GetDirectParentID(), + nextTrack->GetFlag(), + trigIdx); + // Save the neutrino own particle in the store m_data->Stores["ANNIEEvent"]->Set("NeutrinoParticle", neutrino); @@ -1148,6 +1155,7 @@ void LoadWCSim::LoadMCParticles(WCSimRootTrigger* firstTrig) logmessage = "LoadWCSim::LoadMCParticles: Loaded particle with PDG: " + std::to_string(nextTrack->GetIpnu()); logmessage += ", stop time: " + std::to_string(stopTime); logmessage += ", end process: " + nextTrack->GetEndProcess(); + logmessage += ", DirectParentID: " + std::to_string(nextTrack->GetDirectParentID()); Log(logmessage, v_debug, verbosity); // Record neutron primary/secondary @@ -1155,14 +1163,16 @@ void LoadWCSim::LoadMCParticles(WCSimRootTrigger* firstTrig) mapNeutronIsPrim->emplace(nextTrack->GetId(), (nextTrack->GetParenttype() == 0)); MCParticle thisparticle(nextTrack->GetIpnu(), nextTrack->GetE(), nextTrack->GetEndE(), - startPos, stopPos, startTime, stopTime, - Direction(nextTrack->GetDir(0), nextTrack->GetDir(1), nextTrack->GetDir(2)), - length, startStopType, - nextTrack->GetId(), - nextTrack->GetParenttype(), - nextTrack->GetFlag(), - trigIdx); - + startPos, stopPos, startTime, stopTime, + Direction(nextTrack->GetDir(0), nextTrack->GetDir(1), nextTrack->GetDir(2)), + length, startStopType, + nextTrack->GetId(), + nextTrack->GetParenttype(), + nextTrack->GetPrimaryParentID(), + nextTrack->GetDirectParentID(), + nextTrack->GetFlag(), + trigIdx); + // Exit point is not currently in constructor call so set it separately // Older WCSim files do not recor this info. This breaks backward compatibility Position exitPoint(nextTrack->GetTankExitPoint(0), @@ -1171,14 +1181,15 @@ void LoadWCSim::LoadMCParticles(WCSimRootTrigger* firstTrig) exitPoint.UnitToMeter(); thisparticle.SetTankExitPoint(exitPoint); - // Check if this is a primary muon. Only record the first one - if (nextTrack->GetIpnu() == 13 && nextTrack->GetParenttype() == 0 && + // Check if this is a primary muon. Only record the first one + if (nextTrack->GetIpnu() == 13 && nextTrack->GetParenttype() == 0 && //I think muon would always have directparent as muon and not any other particle? (DJA) nextTrack->GetFlag() == 0 && primaryMuonIndex < 0 ) primaryMuonIndex = MCParticles->size(); // Some print outs for "interesting" particles if (abs(nextTrack->GetIpnu()) == 13 || abs(nextTrack->GetIpnu()) == 211 || nextTrack->GetIpnu() == 111){ logmessage = "LoadWCSim::LoadMCParticles: Found " + std::to_string(nextTrack->GetIpnu()); + logmessage += ", with DirectParentID: " + std::to_string(nextTrack->GetDirectParentID()); logmessage += " with flag: " + std::to_string(nextTrack->GetFlag()); logmessage += ", parent type " + std::to_string(nextTrack->GetParenttype()); logmessage += ", Id " + std::to_string(nextTrack->GetId()); @@ -1198,6 +1209,13 @@ void LoadWCSim::LoadMCParticles(WCSimRootTrigger* firstTrig) logmessage = "LoadWCSim::LoadMCParticles: Loaded " + std::to_string(MCParticles->size()) + " MCParticles"; Log(logmessage, v_debug, verbosity); } // end loop over events + + std::cout << "DEBUG: All saved Track IDs in trackid_to_mcparticleindex: "; + for (auto const& pair : *trackid_to_mcparticleindex) { + std::cout << pair.first << " "; + } + +std::cout << std::endl; }// endif MCTriggerNum == 0 else { // if MCTrigger > 0 we need to update all the particle times @@ -1381,7 +1399,9 @@ bool LoadWCSim::LoadHits(WCSimRootTrigger* thisTrig, WCSimRootTrigger* firstTrig Log(logmessage, v_debug, verbosity); // Create the hit and put it in the correct map - MCHit nextHit(key, digiTime, digiQ, GetHitParentIdxs(digiHit, firstTrig)); + std::pair, std::vector> hitParentIDs = GetHitParentIDs(digiHit, firstTrig); + + MCHit nextHit(key, digiTime, digiQ, hitParentIDs.first, hitParentIDs.second); if (system == "Tank") { if (MCHits->count(key) == 0) MCHits->emplace(key, std::vector{nextHit}); @@ -1455,7 +1475,9 @@ void LoadWCSim::MakeParticleToPmtMap(WCSimRootTrigger* thistrig, auto* thehittimeobject = (WCSimRootCherenkovHitTime*)(firstTrig->GetCherenkovHitTimes()->At(thephotonsid)); // get the parent ID from the CherenkovHitTime - Int_t parentID = (thehittimeobject) ? thehittimeobject->GetParentID() : -1; + //I have changed GetParentID to GetDirectParentID in WCSimRootCherenkovHitTime, so this may need to be updated if we want direct parent IDs instead of primary parent IDs (DJA) + Int_t parentID = (thehittimeobject) ? thehittimeobject->GetPrimaryParentID() : -1; + Int_t directparentID = (thehittimeobject) ? thehittimeobject->GetDirectParentID() : -1; // We'll want a map of particle ID to channel keys, so convert WCSim TubeID to channelkey int chankey = tubeid_to_channelkey.at(tubeID); @@ -1484,9 +1506,11 @@ void LoadWCSim::MakeParticleToPmtMap(WCSimRootTrigger* thistrig, //////////////////////////////////////////////////////////////////////////////// // Get the ID of the primary MCParticle(s) that produced this digi. hit -std::vector LoadWCSim::GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig) +//Instead now it stores the direct parent IDs. What do we need primary MCParticles infor too? (DJA) +std::pair, std::vector> LoadWCSim::GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig) { std::vector parentIDs; // a hit could technically have more than one contrbuting particle + std::vector directParentIDs; // loop over the photons in this digit std::vector photonIdxs = digiHit->GetPhotonIds(); @@ -1504,27 +1528,37 @@ std::vector LoadWCSim::GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, logmessage = "LoadWCSim::GetHitParentIDs: HitTime object is NULL!!"; Log(logmessage, v_error, verbosity); } - else - parentIDs.push_back(theHitTimeObject->GetParentID()); + else { + parentIDs.push_back(theHitTimeObject->GetPrimaryParentID()); + directParentIDs.push_back(theHitTimeObject->GetDirectParentID()); + } + }// end loop over photons - return parentIDs; + return std::make_pair(parentIDs, directParentIDs); } //////////////////////////////////////////////////////////////////////////////// // Get the index within the the MCParticle vector of the primaries that produced this digi. hit -std::vector LoadWCSim::GetHitParentIdxs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig) +std::pair, std::vector> LoadWCSim::GetHitParentIdxs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig) { - std::vector parentIDs = GetHitParentIDs(digiHit, firstTrig); + std::pair, std::vector> bothIDs = GetHitParentIDs(digiHit, firstTrig); + std::vector parentIDs = bothIDs.first; + std::vector directParentIDs = bothIDs.second; std::vector parentIdxs; + std::vector directParentIdxs; // Check if the parent was recorded, and if so then translate ID to index for (int parentID : parentIDs) { if (trackid_to_mcparticleindex->count(parentID)) parentIdxs.push_back(trackid_to_mcparticleindex->at(parentID)); } + for (int directParentID : directParentIDs) { + if (trackid_to_mcparticleindex->count(directParentID)) + directParentIdxs.push_back(trackid_to_mcparticleindex->at(directParentID)); + } - return parentIdxs; + return std::make_pair(parentIdxs, directParentIdxs); } //////////////////////////////////////////////////////////////////////////////// diff --git a/UserTools/LoadWCSim/LoadWCSim.h b/UserTools/LoadWCSim/LoadWCSim.h index ad2069ec5..944fd516b 100644 --- a/UserTools/LoadWCSim/LoadWCSim.h +++ b/UserTools/LoadWCSim/LoadWCSim.h @@ -122,8 +122,8 @@ class LoadWCSim: public Tool { // Each MCHit will contain the idx of it's parent MCParticle's // position within the MCParticles vector std::map* trackid_to_mcparticleindex = nullptr; - std::vector GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig); - std::vector GetHitParentIdxs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig); + std::pair, std::vector> GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig); + std::pair, std::vector> GetHitParentIdxs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig); std::map timeArrayOffsetMap; void BuildTimeArrayOffsetMap(WCSimRootTrigger* firstTrig); diff --git a/UserTools/LoadWCSimLAPPD/LoadWCSimLAPPD.cpp b/UserTools/LoadWCSimLAPPD/LoadWCSimLAPPD.cpp index 0be6bcbbb..48f06fe42 100755 --- a/UserTools/LoadWCSimLAPPD/LoadWCSimLAPPD.cpp +++ b/UserTools/LoadWCSimLAPPD/LoadWCSimLAPPD.cpp @@ -291,6 +291,9 @@ bool LoadWCSimLAPPD::Execute(){ digits->Fill(relativedigitst); } std::vector parents; // info about which particle generated the photon for this hit + std::vector directparents; + //#### Do I need to add the DirectParentId for the TrackId_to_MCParticleIndex lookup? + // #### We may be interested in using TrackID_to_MCParticleIndex map for grabbing the pdg code of the given direct parent hit (DJA) if(TrackId_to_MCParticleIndex->count(LAPPDEntry->lappdhit_primaryParentID2->at(runningcount))){ parents.push_back(TrackId_to_MCParticleIndex->at(LAPPDEntry->lappdhit_primaryParentID2->at(runningcount))); } @@ -306,15 +309,15 @@ bool LoadWCSimLAPPD::Execute(){ relativedigitst<(posttriggerwindow) ){ //cout<<"LAPPD hit at absolute time "<count(key)==0){ - MCLAPPDHit nexthit(key, relativedigitst, digiq, globalpos, localpos, parents); + MCLAPPDHit nexthit(key, relativedigitst, digiq, globalpos, localpos, parents, directparents); MCLAPPDHits->emplace(key, std::vector{nexthit}); } else { MCLAPPDHits->at(key).emplace_back(key, relativedigitst, digiq, - globalpos, localpos, parents); + globalpos, localpos, parents, directparents); } if(verbosity>3) cout<<"new lappd digit added"< #include #include +#include // ANNIE includes #include "ANNIEconstants.h" @@ -102,6 +103,9 @@ bool PMTWaveformSim::Execute() // The container for the data that we'll put into the ANNIEEvent std::map> > RawADCDataMC; std::map> > CalADCDataMC; + std::map>> PMTToDirectParentMap; + std::map>> PMTToPrimaryParentMap; + // If MCHits is empty (load_status == 2), create one minimal baseline waveform so that the hit finder doesn't freak out @@ -132,6 +136,7 @@ bool PMTWaveformSim::Execute() std::vector> rawWaveforms; std::vector> calWaveforms; + std::vector directParentIDs; rawWaveforms.emplace_back(0, rawSamples); calWaveforms.emplace_back(0, calSamples, baseline, noiseSigma); @@ -145,14 +150,16 @@ bool PMTWaveformSim::Execute() for (auto mcHitsIt : *fMCHits) { // Loop over the hit PMTs int PMTID = mcHitsIt.first; - std::vector mcHits = mcHitsIt.second; + std::vector &mcHits = mcHitsIt.second; // Generate waveform samples from the MC hits // samples from hits that are close in time will be added together // key is hit time in clock ticks, value is amplitude std::map sample_map; - for (const auto& mcHit : mcHits) {// Loop through each MCHit in the vector + std::map> hits_to_directparents_map; + std::map> hits_to_primaryparents_map; + for (MCHit& mcHit : mcHits) {// Loop through each MCHit in the vector // skip negative hit times, what does that even mean if we're not using the smeared digit time? // skip hit times past 70 us since that's our longest readout if (mcHit.GetTime() < 0) continue; @@ -161,10 +168,12 @@ bool PMTWaveformSim::Execute() // Grab the hit time (also converted to clock ticks) and the charge double hit_t0 = mcHit.GetTime() + fTimeShift; double hit_charge = mcHit.GetCharge(); + const std::vector* directParentIDs = mcHit.GetDirectParents(); + const std::vector* primaryParentIDs = mcHit.GetParents(); - logmessage = "PMTWaveformSim:\n hit charge = " + std::to_string(hit_charge) + " p.e., hit time = " + std::to_string(hit_t0) + " for PMTID " + std::to_string(PMTID); + logmessage = "PMTWaveformSim:\n hit charge = " + std::to_string(hit_charge) + " p.e., hit time = " + std::to_string(hit_t0) + " for PMTID " + std::to_string(PMTID)+ "Direct parent track IDs: " + std::to_string(directParentIDs->size()); Log(logmessage, v_message, verbosity); - + // before "digitizing", add smearing based on the uncertainty extracted in the laser analysis if (fuseTimeSmearing) { double timesmear = TimeSmearing(PMTID); @@ -178,6 +187,10 @@ bool PMTWaveformSim::Execute() uint16_t start_clocktick = (t0_ticks > fPrewindow)? t0_ticks - fPrewindow : 0; uint16_t end_clocktick = start_clocktick + fReadoutWindow; + // Put these ticks into the actual MCHit + mcHit.SetStartTick(start_clocktick); + mcHit.SetEndTick(end_clocktick); + // Randomly Sample the PMT parameters for each MCHit SampleFitParameters(PMTID); @@ -188,19 +201,45 @@ bool PMTWaveformSim::Execute() std::stringstream logmessage; logmessage << " --> clocktick = " << clocktick << ", sample = " << sample; Log(logmessage.str(), v_message, verbosity); - + // check if this hit time has been recorded // either set it or add to it if (sample_map.find(clocktick) == sample_map.end()) sample_map[clocktick] = sample; else - sample_map[clocktick] += sample; - }// end loop over clock ticks + sample_map[clocktick] += sample; + + }// end loop over clock ticks + + // Store parent IDs once per MCHit (at t0 tick) to avoid repeating the same parent info for every sample tick + if (directParentIDs->size() > 0) { + std::vector unique_direct_parent_ids = *directParentIDs; + std::sort(unique_direct_parent_ids.begin(), unique_direct_parent_ids.end()); + unique_direct_parent_ids.erase(std::unique(unique_direct_parent_ids.begin(), + unique_direct_parent_ids.end()), + unique_direct_parent_ids.end()); + hits_to_directparents_map[t0_ticks].insert(hits_to_directparents_map[t0_ticks].end(), + unique_direct_parent_ids.begin(), + unique_direct_parent_ids.end()); + } + + if (primaryParentIDs->size() > 0) { + std::vector unique_primary_parent_ids = *primaryParentIDs; + std::sort(unique_primary_parent_ids.begin(), unique_primary_parent_ids.end()); + unique_primary_parent_ids.erase(std::unique(unique_primary_parent_ids.begin(), + unique_primary_parent_ids.end()), + unique_primary_parent_ids.end()); + hits_to_primaryparents_map[t0_ticks].insert(hits_to_primaryparents_map[t0_ticks].end(), + unique_primary_parent_ids.begin(), + unique_primary_parent_ids.end()); + } + }// end loop over mcHits // If there are no samples for this PMT then no need to do the rest if (sample_map.empty()) continue; - + + // Set the noise envelope and baseline for this PMT // The noise std dev appears to be normally distributed around 1 with sigma 0.25 @@ -211,23 +250,31 @@ bool PMTWaveformSim::Execute() // convert the sample map into a vector of Waveforms and put them into the container std::vector> rawWaveforms; std::vector> calWaveforms; - ConvertMapToWaveforms(sample_map, rawWaveforms, calWaveforms, noiseSigma, basline); + ConvertMapToWaveforms(sample_map, hits_to_directparents_map, rawWaveforms, calWaveforms, noiseSigma, basline); RawADCDataMC.emplace(PMTID, rawWaveforms); CalADCDataMC.emplace(PMTID, calWaveforms); - }// end loop over PMTs + PMTToDirectParentMap[PMTID] = hits_to_directparents_map; + PMTToPrimaryParentMap[PMTID] = hits_to_primaryparents_map; + } // end loop over PMTs + + std::cout << "PMTWaveformSim: Finished looping over MCHits, now publishing waveforms to ANNIEEvent..." << std::endl; // Publish the waveforms to the ANNIEEvent store if we have them m_data->Stores.at("ANNIEEvent")->Set("RawADCDataMC", RawADCDataMC); m_data->Stores.at("ANNIEEvent")->Set("CalibratedADCData", CalADCDataMC); + m_data->Stores.at("ANNIEEvent")->Set("PMTToDirectParentMap", PMTToDirectParentMap); + m_data->Stores.at("ANNIEEvent")->Set("PMTToPrimaryParentMap", PMTToPrimaryParentMap); + m_data->Stores.at("ANNIEEvent")->Set("PMTSimPrewindowTicks", fPrewindow); + m_data->Stores.at("ANNIEEvent")->Set("PMTSimReadoutWindowTicks", fReadoutWindow); + if (fDebug) FillDebugGraphs(RawADCDataMC); return true; } - //------------------------------------------------------------------------------ bool PMTWaveformSim::Finalise() { @@ -441,6 +488,7 @@ uint16_t PMTWaveformSim::CustomLogNormalPulse(double hit_t0, uint16_t clocktick, //------------------------------------------------------------------------------ void PMTWaveformSim::ConvertMapToWaveforms(const std::map &sample_map, + const std::map> &hits_to_directparents_map, std::vector> &rawWaveforms, std::vector> &calWaveforms, double noiseSigma, int baseline) @@ -497,7 +545,6 @@ int PMTWaveformSim::LoadFromStores() return 2; } - return 1; } @@ -505,7 +552,7 @@ int PMTWaveformSim::LoadFromStores() void PMTWaveformSim::FillDebugGraphs(const std::map> > &RawADCDataMC) { for (auto itpair : RawADCDataMC) { - std::string chanString = std::to_string(itpair.first); + std::string chanString = "RawADCData" + std::to_string(itpair.first); // Get/make the directory for this PMT TDirectory* dir = fOutFile->GetDirectory(chanString.c_str()); @@ -563,6 +610,3 @@ double PMTWaveformSim::TimeSmearing(int pmtid) } - - - diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.h b/UserTools/PMTWaveformSim/PMTWaveformSim.h index 8ca705ae1..b3f62be0d 100644 --- a/UserTools/PMTWaveformSim/PMTWaveformSim.h +++ b/UserTools/PMTWaveformSim/PMTWaveformSim.h @@ -44,9 +44,10 @@ class PMTWaveformSim: public Tool { bool SampleFitParameters(int pmtid); ///< sample fit parameters in a way that preserves covariance uint16_t CustomLogNormalPulse(double hit_t0, uint16_t t0_clocktick, double hit_charge); //< construct simulated ADC pulses using sampled fit values void ConvertMapToWaveforms(const std::map &sample_map, ///< construct a waveform with the simulated pulse + baseline modulation, to feed into the PhaseIIADCHitFinder - std::vector> &rawWaveforms, - std::vector> &calWaveforms, - double noiseSigma, int baseline); + const std::map> &hits_to_directparents_map, + std::vector> &rawWaveforms, + std::vector> &calWaveforms, + double noiseSigma, int baseline); void FillDebugGraphs(const std::map> > &RawADCDataMC); ///< debugging double TimeSmearing(int pmtid); ///< prior to sampling the fits, we can add realistic time smearing (instead of relying on WCSim's time smearing) to the MCHit (true) time diff --git a/configfiles/BeamClusterAnalysisMC/ANNIEEventTreeMakerConfig b/configfiles/BeamClusterAnalysisMC/ANNIEEventTreeMakerConfig index 3917232a3..7694fbcf2 100644 --- a/configfiles/BeamClusterAnalysisMC/ANNIEEventTreeMakerConfig +++ b/configfiles/BeamClusterAnalysisMC/ANNIEEventTreeMakerConfig @@ -3,20 +3,21 @@ ANNIEEventTreeMakerVerbosity 0 OutputFile ANNIETree_MC.root -fillAllTriggers 1 +ifillAllTriggers 1 fill_singleTrigger 0 fillLAPPDEventsOnly 0 TankCluster_fill 1 cluster_TankHitInfo_fill 1 TankReco_fill 0 -RingCounting_fill 1 +RingCounting_fill 0 -TankClusterProcessing 1 -MRDClusterProcessing 1 -TriggerProcessing 1 +TankClusterProcessing 0 +MRDClusterProcessing 0 +TriggerProcessing 0 TankHitInfo_fill 1 -MRDHitInfo_fill 1 +DirectParent_MCHit_fill 1 +MRDHitInfo_fill 0 MRDReco_fill 1 SiPMPulseInfo_fill 0 fillCleanEventsOnly 0 @@ -25,11 +26,11 @@ Reco_fill 0 RecoDebug_fill 0 muonTruthRecoDiff_fill 0 isData 0 -HasGenie 1 +HasGenie 0 LAPPDData_fill 0 LAPPDReco_fill 0 RWMBRF_fill 0 LAPPD_PPS_fill 0 LAPPD_Waveform_fill 0 -LAPPD_MC_fill 1 +LAPPD_MC_fill 0 diff --git a/configfiles/BeamClusterAnalysisMC/BackTrackerConfig b/configfiles/BeamClusterAnalysisMC/BackTrackerConfig new file mode 100644 index 000000000..ed1c651c5 --- /dev/null +++ b/configfiles/BeamClusterAnalysisMC/BackTrackerConfig @@ -0,0 +1,2 @@ +verbosity 1 +UseDirectParentClockTickMatching 1 \ No newline at end of file diff --git a/configfiles/BeamClusterAnalysisMC/LoadGenieEventConfig b/configfiles/BeamClusterAnalysisMC/LoadGenieEventConfig index 4da19c900..1a4203631 100644 --- a/configfiles/BeamClusterAnalysisMC/LoadGenieEventConfig +++ b/configfiles/BeamClusterAnalysisMC/LoadGenieEventConfig @@ -1,4 +1,4 @@ -verbosity 0 +verbosity 2 # Refers to format of flux file used to generate GENIE sample (bnb or gsimple) FluxVersion 1 # use 0 to load genie files based on bnb_annie_0000.root etc files (bnb/redecay/dk2nu flux format) @@ -6,19 +6,19 @@ FluxVersion 1 # use 0 to load genie files based on bnb_annie_00 # - All GENIE (and WCSim) files generated by James uses gsimple (1) format # Path to directory of GENIE files -#FileDir NA # special option to load path from associated WCSim (Not recommended) +FileDir NA # special option to load path from associated WCSim (Not recommended) #FileDir /pnfs/annie/persistent/simulations/genie3/G1810a0211a/standard/ -FileDir /pnfs/annie/persistent/simulations/genie3/G1810a0211a/standard/tank - +#FileDir /pnfs/annie/persistent/simulations/genie3/G1810a0211a/standard/tank +FileDir /pnfs/annie/persistent/simulations/genie3/G1810a0211a/standardv1.0/tank/ # Name of GENIE file to open -FilePattern gntp.99.ghep.root # gntp.[run_number].ghep.root +FilePattern gntp.0.ghep.root # gntp.[run_number].ghep.root #FilePattern LoadWCSimTool # special option to load GENIE events corresponding to WCSim events # Option to match GENIE events to WCSim events without using file path saved in WCSim file ManualFileMatching 1 # 0 (false) - no manual matching, 1 (true) - uses FileDir and info from WCSim file name to find GENIE files/events - strongly recommended to use (1)! (as of Aug 2023) # Number of events in the WCSim file (used for offsetting the GENIE file event when used in conjunction with ManualFileMatching) -FileEvents 1000 # 1000 for James' WCSim files, 500 for Marcus' files +FileEvents 2000 # 1000 for James' WCSim files, 500 for Marcus' files # Number to offset the first loaded GENIE event (for usage without ManualFileMatching) - 0 by default EventOffset 0 diff --git a/configfiles/BeamClusterAnalysisMC/LoadWCSimConfig b/configfiles/BeamClusterAnalysisMC/LoadWCSimConfig index 717035a0c..a2acea139 100644 --- a/configfiles/BeamClusterAnalysisMC/LoadWCSimConfig +++ b/configfiles/BeamClusterAnalysisMC/LoadWCSimConfig @@ -7,7 +7,10 @@ verbose 0 # where X = run number for the corresponding GENIE file (gntp.X.ghep.root) # Y = offset multiple for GENIE event number -InputFile /pnfs/annie/persistent/simulations/wcsim/G1810a0211a/standard/tank/pmt/wcsim_0.99.1.root +#InputFile /pnfs/annie/persistent/simulations/wcsim/G1810a0211a/standard/tank/pmt/wcsim_0.99.1.root +#InputFile /exp/annie/app/users/dajana/WCSim/build/wcsim_0.root +InputFile /exp/annie/app/users/dajana/WCSim/build/wcsimwithsectrueformuon_0.root +#InputFile /pnfs/annie/persistent/simulations/wcsim/wcsim_QE_retuning_AmBe_neutrons/QE_1.25/port5_z0_QE_1.25/wcsim_0_999.root WCSimVersion 3 ## should reflect the WCSim version of the files being loaded HistoricTriggeroffset 0 ## time offset of digits relative to the trigger diff --git a/configfiles/BeamClusterAnalysisMC/LoadWCSimLAPPDConfig b/configfiles/BeamClusterAnalysisMC/LoadWCSimLAPPDConfig index 4655cb1f0..0ac7d030e 100644 --- a/configfiles/BeamClusterAnalysisMC/LoadWCSimLAPPDConfig +++ b/configfiles/BeamClusterAnalysisMC/LoadWCSimLAPPDConfig @@ -10,8 +10,8 @@ verbose 0 #/pnfs/annie/persistent/users/mnieslon/wcsim/output/tankonly/wcsim_ANNIEp2v7_throughgoing/wcsim_throughgoing_muon_R2614_lappd_0.0.root #InputFile /pnfs/annie/persistent/simulations/wcsim/wcsim_ANNIEp2v7_beam/lappd-files/LAPPD_wcsim_beam_gst_1_99_0.199.root -InputFile /pnfs/annie/persistent/simulations/wcsim/G1810a0211a/standard/tank/lappd/wcsim_lappd_0.99.1.root - +#InputFile /pnfs/annie/persistent/simulations/wcsim/G1810a0212a/standard/tank/lappd/wcsim_lappd_0.99.1.root +InputFile /exp/annie/app/users/dajana/WCSim/build/wcsimwithsectrueformuon_lappd_0.root WCSimVersion 3 ## should reflect the WCSim version of the files being loaded InnerStructureRadius 1.3545 ## octagonal inner structure radius in m (from drawings 106.64") DrawDebugGraphs 0 ## whether to draw TPolyMarker3D's of hits diff --git a/configfiles/BeamClusterAnalysisMC/PMTWaveformSimConfig b/configfiles/BeamClusterAnalysisMC/PMTWaveformSimConfig new file mode 100644 index 000000000..4233a2a9c --- /dev/null +++ b/configfiles/BeamClusterAnalysisMC/PMTWaveformSimConfig @@ -0,0 +1,8 @@ +verbosity 0 +PMTParameterFile configfiles/PMTWaveformSim/PMTWaveformLognormFit.csv +useTimeSmearing 0 +Prewindow 80 +ReadoutWindow 200 +T0Offset 25 +TimeShift 0 +MakeDebugFile 0 diff --git a/configfiles/BeamClusterAnalysisMC/PhaseIIADCHitFinderConfig b/configfiles/BeamClusterAnalysisMC/PhaseIIADCHitFinderConfig new file mode 100644 index 000000000..d7df3353e --- /dev/null +++ b/configfiles/BeamClusterAnalysisMC/PhaseIIADCHitFinderConfig @@ -0,0 +1,11 @@ +verbosity 0 + +UseLEDWaveforms 0 + +PulseFindingApproach threshold +PulseWindowType Fixed_2023_Gains +DefaultADCThreshold 7 +DefaultThresholdType relative + +EventBuilding 0 +MCWaveforms 1 diff --git a/configfiles/BeamClusterAnalysisMC/ToolsConfig b/configfiles/BeamClusterAnalysisMC/ToolsConfig index 1d60f1abb..682b2eb8e 100644 --- a/configfiles/BeamClusterAnalysisMC/ToolsConfig +++ b/configfiles/BeamClusterAnalysisMC/ToolsConfig @@ -1,7 +1,7 @@ myLoadGeometry LoadGeometry ./configfiles/LoadGeometry/LoadGeometryConfig myLoadWCSim LoadWCSim ./configfiles/BeamClusterAnalysisMC/LoadWCSimConfig myLoadWCSimLAPPD LoadWCSimLAPPD ./configfiles/BeamClusterAnalysisMC/LoadWCSimLAPPDConfig -myLoadGenieEvent LoadGenieEvent ./configfiles/BeamClusterAnalysisMC/LoadGenieEventConfig +#myLoadGenieEvent LoadGenieEvent ./configfiles/BeamClusterAnalysisMC/LoadGenieEventConfig myMCParticleProperties MCParticleProperties ./configfiles/BeamClusterAnalysisMC/MCParticlePropertiesConfig myMCRecoEventLoader MCRecoEventLoader ./configfiles/BeamClusterAnalysisMC/MCRecoEventLoaderConfig @@ -12,8 +12,11 @@ myFindMrdTracks FindMrdTracks configfiles/BeamClusterAnalysisMC/FindMrdTracksCon myClusterFinder ClusterFinder ./configfiles/BeamClusterAnalysisMC/ClusterFinderConfig myClusterClassifiers ClusterClassifiers ./configfiles/BeamClusterAnalysisMC/ClusterClassifiersConfig myEventSelector EventSelector ./configfiles/BeamClusterAnalysisMC/EventSelectorConfig -myCNNImage CNNImage ./configfiles/BeamClusterAnalysisMC/CNNImageConfig -myRingCounting PythonScript configfiles/BeamClusterAnalysisMC/RingCountingConfig +#myCNNImage CNNImage ./configfiles/BeamClusterAnalysisMC/CNNImageConfig +#myRingCounting PythonScript configfiles/BeamClusterAnalysisMC/RingCountingConfig +myPMTWaveformSim PMTWaveformSim configfiles/BeamClusterAnalysisMC/PMTWaveformSimConfig +myPhaseIIADCHitFinder PhaseIIADCHitFinder configfiles/BeamClusterAnalysisMC/PhaseIIADCHitFinderConfig +myBackTracker BackTracker configfiles/BeamClusterAnalysisMC/BackTrackerConfig myANNIEEventTreeMaker ANNIEEventTreeMaker ./configfiles/BeamClusterAnalysisMC/ANNIEEventTreeMakerConfig -myAssignBunchTimingMC AssignBunchTimingMC ./configfiles/BeamClusterAnalysisMC/AssignBunchTimingMCConfig -myPhaseIITreeMaker PhaseIITreeMaker ./configfiles/BeamClusterAnalysisMC/PhaseIITreeMakerConfig +#myAssignBunchTimingMC AssignBunchTimingMC ./configfiles/BeamClusterAnalysisMC/AssignBunchTimingMCConfig +#myPhaseIITreeMaker PhaseIITreeMaker ./configfiles/BeamClusterAnalysisMC/PhaseIITreeMakerConfig diff --git a/configfiles/LoadWCSim/LoadWCSimConfig b/configfiles/LoadWCSim/LoadWCSimConfig index a25d47f59..d8fbd4192 100644 --- a/configfiles/LoadWCSim/LoadWCSimConfig +++ b/configfiles/LoadWCSim/LoadWCSimConfig @@ -7,8 +7,9 @@ verbose 1 # where X = run number for the corresponding GENIE file (gntp.X.ghep.root) # Y = offset multiple for GENIE event number -InputFile /pnfs/annie/persistent/users/moflaher/wcsim/multipmt/tankonly/wcsim_25_04_19_ANNIEp2v6_nodigit_BNB_Water_10k_22-05-17/wcsim_0.1.9.root - +#InputFile /pnfs/annie/persistent/users/moflaher/wcsim/multipmt/tankonly/wcsim_25_04_19_ANNIEp2v6_nodigit_BNB_Water_10k_22-05-17/wcsim_0.1.9.root +#InputFile /pnfs/annie/persistent/simulations/wcsim/G1810a0211a/standardv2.0.0/tank/pmt/wcsim_0.352.1.root +InputFile /exp/annie/app/users/dajana/WCSim/build/wcsim_0.root WCSimVersion 3 ## should reflect the WCSim version of the files being loaded HistoricTriggeroffset 0 ## time offset of digits relative to the trigger UseDigitSmearedTime 1 ## whether to use smeared digit time (T), or true time of first photon (F) diff --git a/configfiles/LoadWCSim/LoadWCSimLAPPDConfig b/configfiles/LoadWCSim/LoadWCSimLAPPDConfig index fec5bd80b..b7d1bd4ed 100644 --- a/configfiles/LoadWCSim/LoadWCSimLAPPDConfig +++ b/configfiles/LoadWCSim/LoadWCSimLAPPDConfig @@ -4,8 +4,7 @@ verbose 1 #InputFile /pnfs/annie/persistent/users/moflaher/wcsim/lappd/tankonly/wcsim_lappd_tankonly_24-09-17_BNB_Water_10k_22-05-17/wcsim_lappd_0.0.0.root #InputFile /pnfs/annie/persistent/users/moflaher/wcsim/lappd/tankonly/wcsim_lappd_tankonly_03-05-17_rhatcher/wcsim_lappd_0.1000.root ## first of the DOE proposal files -InputFile /pnfs/annie/persistent/users/moflaher/wcsim/multipmt/tankonly/wcsim_3-12-18_ANNIEp2v6_BNB_Water_10k_22-05-17/wcsim_lappd_0.0.0.root - +InputFile /pnfs/annie/persistent/simulations/wcsim/G1810a0211a/standardv2.0.0/tank/lappd/wcsim_lappd_0.352.0.root WCSimVersion 3 ## should reflect the WCSim version of the files being loaded InnerStructureRadius 1.3545 ## octagonal inner structure radius in m (from drawings 106.64") DrawDebugGraphs 0 ## whether to draw TPolyMarker3D's of hits diff --git a/configfiles/LoadWCSim/ToolsConfig b/configfiles/LoadWCSim/ToolsConfig index a03c1c489..6c8326d62 100644 --- a/configfiles/LoadWCSim/ToolsConfig +++ b/configfiles/LoadWCSim/ToolsConfig @@ -1,7 +1,7 @@ myLoadWCSim LoadWCSim ./configfiles/LoadWCSim/LoadWCSimConfig -myLoadWCSimLAPPD LoadWCSimLAPPD ./configfiles/LoadWCSim/LoadWCSimLAPPDConfig +#myLoadWCSimLAPPD LoadWCSimLAPPD ./configfiles/LoadWCSim/LoadWCSimLAPPDConfig #myPrintANNIEEvent PrintANNIEEvent ./configfiles/LoadWCSim/PrintANNIEEventConfig -myHitTimeResiduals HitResiduals ./configfiles/LoadWCSim/HitResidualsConfig +#myHitTimeResiduals HitResiduals ./configfiles/LoadWCSim/HitResidualsConfig #myPlotLAPPDTimesFromStore PlotLAPPDTimesFromStore ./configfiles/LoadWCSim/PlotLAPPDTimesFromStoreConfig myCheckDetectorCounts CheckDetectorCounts ./configfiles/LoadWCSim/CheckDetectorCountsConfig diff --git a/configfiles/PMTWaveformSim/BackTrackerConfig b/configfiles/PMTWaveformSim/BackTrackerConfig new file mode 100644 index 000000000..ed1c651c5 --- /dev/null +++ b/configfiles/PMTWaveformSim/BackTrackerConfig @@ -0,0 +1,2 @@ +verbosity 1 +UseDirectParentClockTickMatching 1 \ No newline at end of file diff --git a/configfiles/PMTWaveformSim/ClusterFinderConfig b/configfiles/PMTWaveformSim/ClusterFinderConfig index ded568144..e0d3c4120 100644 --- a/configfiles/PMTWaveformSim/ClusterFinderConfig +++ b/configfiles/PMTWaveformSim/ClusterFinderConfig @@ -1,7 +1,7 @@ # ClusterFinder Config File verbosity 0 -HitStore Hits #Either MCHits or Hits (accessed in ANNIEEvent store) +HitStore MCHits #Either MCHits or Hits (accessed in ANNIEEvent store) OutputFile BeamRun_ClusterFinder_DefaultOutput #Output root prefix name for the current run ClusterFindingWindow 40 # in ns, size of the window used to "clusterize" AcqTimeWindow 70000 # in ns, size of the acquisition window diff --git a/configfiles/PMTWaveformSim/LoadWCSimConfig b/configfiles/PMTWaveformSim/LoadWCSimConfig index 810e60949..c456453e7 100644 --- a/configfiles/PMTWaveformSim/LoadWCSimConfig +++ b/configfiles/PMTWaveformSim/LoadWCSimConfig @@ -1,6 +1,6 @@ verbose 1 -InputFile wcsim_0.root - +#InputFile wcsim_0.root +InputFile /exp/annie/app/users/dajana/WCSim/build/wcsim_0.root WCSimVersion 3 HistoricTriggeroffset 0 UseDigitSmearedTime 0 diff --git a/configfiles/PMTWaveformSim/ToolsConfig b/configfiles/PMTWaveformSim/ToolsConfig index 782b19507..2e18d8e2f 100644 --- a/configfiles/PMTWaveformSim/ToolsConfig +++ b/configfiles/PMTWaveformSim/ToolsConfig @@ -2,4 +2,5 @@ LoadGeometry LoadGeometry configfiles/LoadGeometry/LoadGeometryConfig LoadWCSim LoadWCSim configfiles/PMTWaveformSim/LoadWCSimConfig PMTWaveformSim PMTWaveformSim configfiles/PMTWaveformSim/PMTWaveformSimConfig PhaseIIADCHitFinder PhaseIIADCHitFinder configfiles/PMTWaveformSim/PhaseIIADCHitFinderConfig -ClusterFinder ClusterFinder configfiles/PMTWaveformSim/ClusterFinderConfig \ No newline at end of file +ClusterFinder ClusterFinder configfiles/PMTWaveformSim/ClusterFinderConfig +BackTracker BackTracker configfiles/PMTWaveformSim/BackTrackerConfig From 7fd59e8d4343395c4d7b5018aff928136a203b69 Mon Sep 17 00:00:00 2001 From: dhavvval Date: Thu, 9 Apr 2026 14:21:07 -0500 Subject: [PATCH 2/6] Ideation of tagging primary and secondary neutrons aka Neutron Classification --- .../ANNIEEventTreeMaker.cpp | 31 ++++++++++++- .../ANNIEEventTreeMaker/ANNIEEventTreeMaker.h | 3 ++ UserTools/BackTracker/BackTracker.cpp | 45 +++++++++++++++---- UserTools/BackTracker/BackTracker.h | 4 ++ 4 files changed, 74 insertions(+), 9 deletions(-) diff --git a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp index 5a1d6532a..d6419e48c 100644 --- a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp +++ b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp @@ -189,6 +189,9 @@ bool ANNIEEventTreeMaker::Initialise(std::string configfile, DataModel &data) fANNIETree->Branch("DirectParent_PDGs", &fDirectParent_PDGs); fANNIETree->Branch("DirectParent_NeutronAncestorTrackID", &fDirectParent_NeutronAncestorTrackID); fANNIETree->Branch("DirectParent_NeutronAncestorPDG", &fDirectParent_NeutronAncestorPDG); + fANNIETree->Branch("DirectParent_NeutronAncestorClass", &fDirectParent_NeutronAncestorClass); + fANNIETree->Branch("DirectParent_NeutronParentTrackID", &fDirectParent_NeutronParentTrackID); + fANNIETree->Branch("DirectParent_NeutronParentPDG", &fDirectParent_NeutronParentPDG); } if (SiPMPulseInfo_fill) @@ -787,6 +790,9 @@ void ANNIEEventTreeMaker::ResetVariables() fDirectParent_PDGs.clear(); fDirectParent_NeutronAncestorTrackID.clear(); fDirectParent_NeutronAncestorPDG.clear(); + fDirectParent_NeutronAncestorClass.clear(); + fDirectParent_NeutronParentTrackID.clear(); + fDirectParent_NeutronParentPDG.clear(); // SiPMPulse Info fSiPM1NPulses = 0; @@ -1418,6 +1424,8 @@ void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ std::vector *fMCParticles = nullptr; std::map *fTrackIdToIndex = nullptr; std::map>> *fMCHitToNeutronAncestor = nullptr; + std::map> *fMCHitToNeutronAncestorClass = nullptr; + std::map>> *fMCHitToNeutronParent = nullptr; bool got_MCHitToDirectParents = m_data->Stores["ANNIEEvent"]->Get("MCHitToDirectParents", fMCHitToDirectParents); if (!got_MCHitToDirectParents) { @@ -1442,6 +1450,8 @@ void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ std::cout << "No MCHitToNeutronAncestor store in ANNIEEvent. Continuing to build tree " << std::endl; return; } + bool got_neutronAncestorClass = m_data->Stores["ANNIEEvent"]->Get("MCHitToNeutronAncestorClass", fMCHitToNeutronAncestorClass); + bool got_neutronParent = m_data->Stores["ANNIEEvent"]->Get("MCHitToNeutronParent", fMCHitToNeutronParent); for (auto const& apair : *fMCHitToDirectParents) { unsigned long pmtID = apair.first; @@ -1479,6 +1489,9 @@ void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ int neutronAncestorTrackID = -5; int neutronAncestorPDG = -5; + int neutronAncestorClass = -5; + int neutronParentTrackID = -5; + int neutronParentPDG = -5; if (got_neutronAncestor && fMCHitToNeutronAncestor->find(pmtID) != fMCHitToNeutronAncestor->end()){ auto const& pmtAncestors = fMCHitToNeutronAncestor->at(pmtID); if (pmtAncestors.find(hitTime) != pmtAncestors.end()){ @@ -1487,8 +1500,25 @@ void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ neutronAncestorPDG = ancestorPair.second; } } + if (got_neutronAncestorClass && fMCHitToNeutronAncestorClass->find(pmtID) != fMCHitToNeutronAncestorClass->end()){ + auto const& pmtClasses = fMCHitToNeutronAncestorClass->at(pmtID); + if (pmtClasses.find(hitTime) != pmtClasses.end()){ + neutronAncestorClass = pmtClasses.at(hitTime); + } + } + if (got_neutronParent && fMCHitToNeutronParent->find(pmtID) != fMCHitToNeutronParent->end()){ + auto const& pmtParents = fMCHitToNeutronParent->at(pmtID); + if (pmtParents.find(hitTime) != pmtParents.end()){ + auto const& parentPair = pmtParents.at(hitTime); + neutronParentTrackID = parentPair.first; + neutronParentPDG = parentPair.second; + } + } fDirectParent_NeutronAncestorTrackID.push_back(neutronAncestorTrackID); fDirectParent_NeutronAncestorPDG.push_back(neutronAncestorPDG); + fDirectParent_NeutronAncestorClass.push_back(neutronAncestorClass); + fDirectParent_NeutronParentTrackID.push_back(neutronParentTrackID); + fDirectParent_NeutronParentPDG.push_back(neutronParentPDG); } } return; @@ -2913,4 +2943,3 @@ tuple ANNIEEventTreeMaker::queryNearestACCID(const vector> fDirectParent_PDGs; std::vector fDirectParent_NeutronAncestorTrackID; //Stored unique neutron ancestor track ID for each MCHit, if it exists. -5 if no neutron ancestor found. std::vector fDirectParent_NeutronAncestorPDG; //For now, we are storing Neutron's truth information. But, it could be exapand to other particle types if needed. -5 if no neutron ancestor found. + std::vector fDirectParent_NeutronAncestorClass; // -5 none, 1 primary, 2 secondary-from-proton, 3 secondary-from-neutron, 4 secondary-other + std::vector fDirectParent_NeutronParentTrackID; // direct parent track ID of the stored neutron ancestor + std::vector fDirectParent_NeutronParentPDG; // direct parent PDG of the stored neutron ancestor // SiPMPulseInfo_fill int fSiPM1NPulses; diff --git a/UserTools/BackTracker/BackTracker.cpp b/UserTools/BackTracker/BackTracker.cpp index d0d1b316e..680d40c44 100644 --- a/UserTools/BackTracker/BackTracker.cpp +++ b/UserTools/BackTracker/BackTracker.cpp @@ -39,8 +39,10 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){ fClusterEfficiency = new std::map; fClusterPurity = new std::map; fClusterTotalCharge = new std::map; - fMCHitToDirectParents = new std::map>>; - fMCHitToNeutronAncestor = new std::map>>; + fMCHitToDirectParents = new std::map>>; + fMCHitToNeutronAncestor = new std::map>>; + fMCHitToNeutronAncestorClass = new std::map>; + fMCHitToNeutronParent = new std::map>>; return true; } @@ -56,7 +58,9 @@ bool BackTracker::Execute() fClusterPurity ->clear(); fClusterTotalCharge ->clear(); fMCHitToDirectParents ->clear(); - // fMCHitToNeutronAncestor ->clear(); + fMCHitToNeutronAncestor ->clear(); + fMCHitToNeutronAncestorClass->clear(); + fMCHitToNeutronParent->clear(); fParticleToTankTotalCharge.clear(); @@ -93,6 +97,8 @@ bool BackTracker::Execute() m_data->Stores.at("ANNIEEvent")->Set("ClusterTotalCharge", fClusterTotalCharge ); m_data->Stores.at("ANNIEEvent")->Set("MCHitToDirectParents", fMCHitToDirectParents ); m_data->Stores.at("ANNIEEvent")->Set("MCHitToNeutronAncestor", fMCHitToNeutronAncestor ); + m_data->Stores.at("ANNIEEvent")->Set("MCHitToNeutronAncestorClass", fMCHitToNeutronAncestorClass); + m_data->Stores.at("ANNIEEvent")->Set("MCHitToNeutronParent", fMCHitToNeutronParent); //It stores return true; } @@ -229,13 +235,13 @@ void BackTracker::DirectParentsFromClockTickWindows() void BackTracker::FindNeutronAncestors() { - std::map> trackMap; // trackId -> (ParentID, pdg) + std::map> trackMap; // trackId -> (DirectParentID, pdg) for (auto& particle : *fMCParticles) { int trackId = particle.GetParticleID(); - int parentId = particle.GetDirectParentID(); + int directParentId = particle.GetDirectParentID(); int pdg = particle.GetPdgCode(); - trackMap[trackId] = std::make_pair(parentId, pdg); + trackMap[trackId] = std::make_pair(directParentId, pdg); } for (const auto& pmtPair : *fMCHitToDirectParents) { @@ -248,7 +254,15 @@ void BackTracker::FindNeutronAncestors() { int neutronAncestorId = -5; int neutronAncestorPdg = -5; - int startParentID = directParents[0]; // take the first direct parent as the starting point + int neutronAncestorClass = -5; + //Neutron Classification: + // 1: primary neutron from initial interaction boundary + // 2: secondary neutron from proton + // 3: secondary neutron from neutron + // 4: secondary neutron from other parent type + int neutronParentTrackId = -5; + int neutronParentPdg = -5; + int startParentID = directParents[0]; // take the first direct parent of neutron as the starting point int currentID = startParentID; @@ -267,8 +281,21 @@ void BackTracker::FindNeutronAncestors() { if (pdg == 2112 && neutronAncestorId == -5) { // if it's a neutron and we haven't already found an ancestor, save it neutronAncestorId = currentID; // This is the neutron's TRACK ID neutronAncestorPdg = pdg; // Store the PDG code (2112 for neutron) + neutronParentTrackId = parentID; + neutronParentPdg = (trackMap.count(parentID) ? trackMap[parentID].second : -5); + + // Neutron classification based on parentage + if (neutronParentTrackId == 0) { + neutronAncestorClass = 1; // primary neutron from initial interaction boundary + } else if (neutronParentPdg == 2212) { + neutronAncestorClass = 2; // secondary neutron from proton + } else if (neutronParentPdg == 2112) { + neutronAncestorClass = 3; // secondary neutron from neutron + } else { + neutronAncestorClass = 4; // secondary neutron from other parent type + } break; //We care about the immidiate neutron ancestor. - //Do we want to find the potential primary neutron ancestor, if there is one? (DJA) + //Do we want to find the potential primary neutron ancestor, if there is one? (DJA) - If the classification portion works that this question is answered. } if (parentID == -1) break; // reached the end of the ancestry @@ -276,6 +303,8 @@ void BackTracker::FindNeutronAncestors() { } (*fMCHitToNeutronAncestor)[pmtID][hitTime] = std::make_pair(neutronAncestorId, neutronAncestorPdg); + (*fMCHitToNeutronAncestorClass)[pmtID][hitTime] = neutronAncestorClass; + (*fMCHitToNeutronParent)[pmtID][hitTime] = std::make_pair(neutronParentTrackId, neutronParentPdg); } } diff --git a/UserTools/BackTracker/BackTracker.h b/UserTools/BackTracker/BackTracker.h index 266ade505..43b6455fb 100644 --- a/UserTools/BackTracker/BackTracker.h +++ b/UserTools/BackTracker/BackTracker.h @@ -70,6 +70,10 @@ class BackTracker: public Tool { // PMT ID -> reco hit time -> (neutron trackID, neutron PDG) std::map>> *fMCHitToNeutronAncestor = nullptr; + // PMT ID -> reco hit time -> neutron class (-5 none, 1 primary, 2 secondary-from-proton, 3 secondary-from-neutron, 4 secondary-other) + std::map> *fMCHitToNeutronAncestorClass = nullptr; + // PMT ID -> reco hit time -> (neutron direct parent trackID, neutron direct parent PDG) + std::map>> *fMCHitToNeutronParent = nullptr; bool fDirectParentClockTickMatching = true; uint16_t fPMTSimPrewindowTicks = 10; From 94a05ff78040f09a8e39b3a5367aa06d9cbfe2d9 Mon Sep 17 00:00:00 2001 From: Dhavalkumar Ajana Date: Wed, 22 Apr 2026 11:33:31 -0500 Subject: [PATCH 3/6] Added a feature about storing the darknoise MCHits to train OPTICS against it for Neutron clustering --- DataModel/Hit.h | 23 ++++++++--- DataModel/LAPPDHit.h | 18 +++++++-- .../ANNIEEventTreeMaker.cpp | 13 +++++++ .../ANNIEEventTreeMaker/ANNIEEventTreeMaker.h | 3 +- UserTools/BackTracker/BackTracker.cpp | 26 ++++++++++++- UserTools/BackTracker/BackTracker.h | 16 ++++++-- UserTools/LoadWCSim/LoadWCSim.cpp | 38 ++++++++++++------- UserTools/LoadWCSim/LoadWCSim.h | 4 +- 8 files changed, 112 insertions(+), 29 deletions(-) diff --git a/DataModel/Hit.h b/DataModel/Hit.h index 9ef6059c7..b654e7bf9 100755 --- a/DataModel/Hit.h +++ b/DataModel/Hit.h @@ -78,6 +78,7 @@ class MCHit : public Hit { , DirectParents(std::vector{}) , StartTick(-5) , EndTick(-5) + , IsDarknoise(false) { serialise=true; } @@ -89,6 +90,7 @@ class MCHit : public Hit { , DirectParents(thedirectparents) , StartTick(-5) , EndTick(-5) + , IsDarknoise(false) { serialise=true; } @@ -99,11 +101,13 @@ class MCHit : public Hit { const std::vector* GetDirectParents() { return &DirectParents; } int GetStartTick() { return StartTick; } int GetEndTick() { return EndTick; } - + bool GetIsDarknoise() const { return IsDarknoise; } + void SetParents(std::vector parentsin) { Parents = parentsin; } void SetDirectParents(std::vector directparentsin) { DirectParents = directparentsin; } void SetStartTick(int tick) { StartTick = tick; } void SetEndTick(int tick) { EndTick = tick; } + void SetIsDarknoise(bool v) { IsDarknoise = v; } bool Print() @@ -123,16 +127,19 @@ class MCHit : public Hit { } else { std::cout << "No recorded parents" << std::endl; } - + + std::cout << "IsDarknoise : " << (IsDarknoise ? "true" : "false") << std::endl; + return true; } - + protected: std::vector Parents; std::vector DirectParents; int StartTick; int EndTick; - + bool IsDarknoise; + template void serialize(Archive & ar, const unsigned int version) { if (serialise) { @@ -142,14 +149,18 @@ class MCHit : public Hit { if (version > 0) { ar & Parents; // Parents is now track IDs rather than index within vector - ar & DirectParents; + ar & DirectParents; ar & StartTick; ar & EndTick; } + + if (version > 1) { + ar & IsDarknoise; + } } } }; -BOOST_CLASS_VERSION(MCHit, 1) +BOOST_CLASS_VERSION(MCHit, 2) #endif diff --git a/DataModel/LAPPDHit.h b/DataModel/LAPPDHit.h index d84604f37..565c84545 100755 --- a/DataModel/LAPPDHit.h +++ b/DataModel/LAPPDHit.h @@ -100,8 +100,8 @@ class MCLAPPDHit : public LAPPDHit friend class boost::serialization::access; public: - MCLAPPDHit() : LAPPDHit(), Parents(std::vector{}), DirectParents(std::vector{}) { serialise = true; } - MCLAPPDHit(int thetubeid, double thetime, double thecharge, std::vector theposition, std::vector thelocalposition, std::vector theparents, std::vector thedirectparents) : LAPPDHit(thetubeid, thetime, thecharge, theposition, thelocalposition), Parents(theparents), DirectParents(thedirectparents) { serialise = true; } + MCLAPPDHit() : LAPPDHit(), Parents(std::vector{}), DirectParents(std::vector{}), IsDarknoise(false) { serialise = true; } + MCLAPPDHit(int thetubeid, double thetime, double thecharge, std::vector theposition, std::vector thelocalposition, std::vector theparents, std::vector thedirectparents) : LAPPDHit(thetubeid, thetime, thecharge, theposition, thelocalposition), Parents(theparents), DirectParents(thedirectparents), IsDarknoise(false) { serialise = true; } const std::vector *GetParents() const { return &Parents; } void SetParents(std::vector parentsin) { Parents = parentsin; } @@ -109,6 +109,9 @@ class MCLAPPDHit : public LAPPDHit const std::vector *GetDirectParents() const { return &DirectParents; } void SetDirectParents(std::vector directparentsin) { DirectParents = directparentsin; } + bool GetIsDarknoise() const { return IsDarknoise; } + void SetIsDarknoise(bool v) { IsDarknoise = v; } + bool Print() { cout << "TubeId : " << TubeId << endl; @@ -151,7 +154,9 @@ class MCLAPPDHit : public LAPPDHit { cout << "##### No recorded Direct parents #####" << endl; } - + + cout << "IsDarknoise : " << (IsDarknoise ? "true" : "false") << endl; + return true; } @@ -167,14 +172,21 @@ class MCLAPPDHit : public LAPPDHit ar & Charge; // n.b. at time of writing MCHit stores no additional persistent members // - it only adds parent MCParticle indices, and these aren't saved... + + if (version > 0) { + ar & IsDarknoise; + } } } protected: std::vector Parents; std::vector DirectParents; + bool IsDarknoise; }; +BOOST_CLASS_VERSION(MCLAPPDHit, 1) + /* class TDCHit : public Hit { public: diff --git a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp index d6419e48c..f3f26cd61 100644 --- a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp +++ b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp @@ -192,6 +192,7 @@ bool ANNIEEventTreeMaker::Initialise(std::string configfile, DataModel &data) fANNIETree->Branch("DirectParent_NeutronAncestorClass", &fDirectParent_NeutronAncestorClass); fANNIETree->Branch("DirectParent_NeutronParentTrackID", &fDirectParent_NeutronParentTrackID); fANNIETree->Branch("DirectParent_NeutronParentPDG", &fDirectParent_NeutronParentPDG); + fANNIETree->Branch("DirectParent_IsDarknoise", &fDirectParent_IsDarknoise); } if (SiPMPulseInfo_fill) @@ -793,6 +794,7 @@ void ANNIEEventTreeMaker::ResetVariables() fDirectParent_NeutronAncestorClass.clear(); fDirectParent_NeutronParentTrackID.clear(); fDirectParent_NeutronParentPDG.clear(); + fDirectParent_IsDarknoise.clear(); // SiPMPulse Info fSiPM1NPulses = 0; @@ -1426,6 +1428,7 @@ void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ std::map>> *fMCHitToNeutronAncestor = nullptr; std::map> *fMCHitToNeutronAncestorClass = nullptr; std::map>> *fMCHitToNeutronParent = nullptr; + std::map> *fMCHitToIsDarknoise = nullptr; bool got_MCHitToDirectParents = m_data->Stores["ANNIEEvent"]->Get("MCHitToDirectParents", fMCHitToDirectParents); if (!got_MCHitToDirectParents) { @@ -1452,6 +1455,7 @@ void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ } bool got_neutronAncestorClass = m_data->Stores["ANNIEEvent"]->Get("MCHitToNeutronAncestorClass", fMCHitToNeutronAncestorClass); bool got_neutronParent = m_data->Stores["ANNIEEvent"]->Get("MCHitToNeutronParent", fMCHitToNeutronParent); + bool got_isDarknoise = m_data->Stores["ANNIEEvent"]->Get("MCHitToIsDarknoise", fMCHitToIsDarknoise); for (auto const& apair : *fMCHitToDirectParents) { unsigned long pmtID = apair.first; @@ -1519,6 +1523,15 @@ void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ fDirectParent_NeutronAncestorClass.push_back(neutronAncestorClass); fDirectParent_NeutronParentTrackID.push_back(neutronParentTrackID); fDirectParent_NeutronParentPDG.push_back(neutronParentPDG); + + int isDarknoise = 0; + if (got_isDarknoise && fMCHitToIsDarknoise->find(pmtID) != fMCHitToIsDarknoise->end()){ + auto const& pmtNoise = fMCHitToIsDarknoise->at(pmtID); + if (pmtNoise.find(hitTime) != pmtNoise.end()){ + isDarknoise = pmtNoise.at(hitTime) ? 1 : 0; + } + } + fDirectParent_IsDarknoise.push_back(isDarknoise); } } return; diff --git a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h index ca9222d6b..056294b35 100644 --- a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h +++ b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h @@ -227,9 +227,10 @@ class ANNIEEventTreeMaker : public Tool std::vector> fDirectParent_PDGs; std::vector fDirectParent_NeutronAncestorTrackID; //Stored unique neutron ancestor track ID for each MCHit, if it exists. -5 if no neutron ancestor found. std::vector fDirectParent_NeutronAncestorPDG; //For now, we are storing Neutron's truth information. But, it could be exapand to other particle types if needed. -5 if no neutron ancestor found. - std::vector fDirectParent_NeutronAncestorClass; // -5 none, 1 primary, 2 secondary-from-proton, 3 secondary-from-neutron, 4 secondary-other + std::vector fDirectParent_NeutronAncestorClass; // 0 dark-noise, -5 non-neutron physics, 1 primary, 2 secondary-from-proton, 3 secondary-from-neutron, 4 secondary-other std::vector fDirectParent_NeutronParentTrackID; // direct parent track ID of the stored neutron ancestor std::vector fDirectParent_NeutronParentPDG; // direct parent PDG of the stored neutron ancestor + std::vector fDirectParent_IsDarknoise; // 1 if the pulse is pure dark noise, 0 otherwise // SiPMPulseInfo_fill int fSiPM1NPulses; diff --git a/UserTools/BackTracker/BackTracker.cpp b/UserTools/BackTracker/BackTracker.cpp index 680d40c44..ec2da4a51 100644 --- a/UserTools/BackTracker/BackTracker.cpp +++ b/UserTools/BackTracker/BackTracker.cpp @@ -1,5 +1,6 @@ #include "BackTracker.h" #include "ANNIEconstants.h" +#include BackTracker::BackTracker():Tool(){} @@ -43,6 +44,7 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){ fMCHitToNeutronAncestor = new std::map>>; fMCHitToNeutronAncestorClass = new std::map>; fMCHitToNeutronParent = new std::map>>; + fMCHitToIsDarknoise = new std::map>; return true; } @@ -61,6 +63,7 @@ bool BackTracker::Execute() fMCHitToNeutronAncestor ->clear(); fMCHitToNeutronAncestorClass->clear(); fMCHitToNeutronParent->clear(); + fMCHitToIsDarknoise->clear(); fParticleToTankTotalCharge.clear(); @@ -98,7 +101,8 @@ bool BackTracker::Execute() m_data->Stores.at("ANNIEEvent")->Set("MCHitToDirectParents", fMCHitToDirectParents ); m_data->Stores.at("ANNIEEvent")->Set("MCHitToNeutronAncestor", fMCHitToNeutronAncestor ); m_data->Stores.at("ANNIEEvent")->Set("MCHitToNeutronAncestorClass", fMCHitToNeutronAncestorClass); - m_data->Stores.at("ANNIEEvent")->Set("MCHitToNeutronParent", fMCHitToNeutronParent); //It stores + m_data->Stores.at("ANNIEEvent")->Set("MCHitToNeutronParent", fMCHitToNeutronParent); //It stores + m_data->Stores.at("ANNIEEvent")->Set("MCHitToIsDarknoise", fMCHitToIsDarknoise); return true; } @@ -252,16 +256,36 @@ void BackTracker::FindNeutronAncestors() { if (directParents.empty()) continue; + // Dark-noise check: a pulse is pure dark noise iff every contributing + // MCHit was pure noise. Noise photons carry direct parent -1 (WCSim + // convention), so the pulse's directParents vector is all -1 only when + // no real physics MCHit contributed within the pulse window. This is + // equivalent to MCHit::GetIsDarknoise() set in LoadWCSim, but keyed on + // the pulse peak_time which is what fMCHitToDirectParents uses. + bool isDarknoise = std::all_of(directParents.begin(), directParents.end(), + [](int id) { return id == -1; }); + (*fMCHitToIsDarknoise)[pmtID][hitTime] = isDarknoise; + int neutronAncestorId = -5; int neutronAncestorPdg = -5; int neutronAncestorClass = -5; //Neutron Classification: + // 0: dark noise (pure-noise pulse) // 1: primary neutron from initial interaction boundary // 2: secondary neutron from proton // 3: secondary neutron from neutron // 4: secondary neutron from other parent type int neutronParentTrackId = -5; int neutronParentPdg = -5; + + if (isDarknoise) { + neutronAncestorClass = 0; + (*fMCHitToNeutronAncestor)[pmtID][hitTime] = std::make_pair(neutronAncestorId, neutronAncestorPdg); + (*fMCHitToNeutronAncestorClass)[pmtID][hitTime] = neutronAncestorClass; + (*fMCHitToNeutronParent)[pmtID][hitTime] = std::make_pair(neutronParentTrackId, neutronParentPdg); + continue; + } + int startParentID = directParents[0]; // take the first direct parent of neutron as the starting point int currentID = startParentID; diff --git a/UserTools/BackTracker/BackTracker.h b/UserTools/BackTracker/BackTracker.h index 43b6455fb..3473e5b05 100644 --- a/UserTools/BackTracker/BackTracker.h +++ b/UserTools/BackTracker/BackTracker.h @@ -16,9 +16,9 @@ * * A tool to link reco info to the paticle(s) that generated the light * -* $Author: A.Sutton $ -* $Date: 2024/06/16 $ -* Contact: atcsutton@gmail.com +* $Author: A.Sutton $ -> D. Ajana +* $Date: 2024/06/16 $ -> 2026/01/26 +* Contact: atcsutton@gmail.com -> dja23@fsu.edu */ class BackTracker: public Tool { @@ -70,10 +70,18 @@ class BackTracker: public Tool { // PMT ID -> reco hit time -> (neutron trackID, neutron PDG) std::map>> *fMCHitToNeutronAncestor = nullptr; - // PMT ID -> reco hit time -> neutron class (-5 none, 1 primary, 2 secondary-from-proton, 3 secondary-from-neutron, 4 secondary-other) + // PMT ID -> reco hit time -> neutron class + // 0 dark noise (pure-noise pulse; all contributing MCHits have primary parent -1) + // 1 primary neutron from initial interaction boundary + // 2 secondary neutron from proton + // 3 secondary neutron from neutron + // 4 secondary neutron from other parent type + // -5 non-neutron physics background aka MCHits from other particles (no neutron in ancestry) std::map> *fMCHitToNeutronAncestorClass = nullptr; // PMT ID -> reco hit time -> (neutron direct parent trackID, neutron direct parent PDG) std::map>> *fMCHitToNeutronParent = nullptr; + // PMT ID -> reco hit time -> true if pulse is pure dark noise + std::map> *fMCHitToIsDarknoise = nullptr; bool fDirectParentClockTickMatching = true; uint16_t fPMTSimPrewindowTicks = 10; diff --git a/UserTools/LoadWCSim/LoadWCSim.cpp b/UserTools/LoadWCSim/LoadWCSim.cpp index c34df9311..b7a998b42 100644 --- a/UserTools/LoadWCSim/LoadWCSim.cpp +++ b/UserTools/LoadWCSim/LoadWCSim.cpp @@ -1399,9 +1399,13 @@ bool LoadWCSim::LoadHits(WCSimRootTrigger* thisTrig, WCSimRootTrigger* firstTrig Log(logmessage, v_debug, verbosity); // Create the hit and put it in the correct map - std::pair, std::vector> hitParentIDs = GetHitParentIDs(digiHit, firstTrig); + std::vector primaryParents; + std::vector directParents; + bool hitIsDarknoise = false; + std::tie(primaryParents, directParents, hitIsDarknoise) = GetHitParentIDs(digiHit, firstTrig); - MCHit nextHit(key, digiTime, digiQ, hitParentIDs.first, hitParentIDs.second); + MCHit nextHit(key, digiTime, digiQ, primaryParents, directParents); + nextHit.SetIsDarknoise(hitIsDarknoise); if (system == "Tank") { if (MCHits->count(key) == 0) MCHits->emplace(key, std::vector{nextHit}); @@ -1507,44 +1511,52 @@ void LoadWCSim::MakeParticleToPmtMap(WCSimRootTrigger* thistrig, //////////////////////////////////////////////////////////////////////////////// // Get the ID of the primary MCParticle(s) that produced this digi. hit //Instead now it stores the direct parent IDs. What do we need primary MCParticles infor too? (DJA) -std::pair, std::vector> LoadWCSim::GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig) +std::tuple, std::vector, bool> LoadWCSim::GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig) { std::vector parentIDs; // a hit could technically have more than one contrbuting particle - std::vector directParentIDs; - + std::vector directParentIDs; + // loop over the photons in this digit std::vector photonIdxs = digiHit->GetPhotonIds(); + + // Dark-noise photons in WCSim carry GetPrimaryParentID() == -1; signal photons >= 1. + // Flag the MCHit as dark-noise only if every contributing photon is a noise photon. + bool isDarknoise = !photonIdxs.empty(); + for (int photonIdx : photonIdxs) { // Special offset for older WCSim if (WCSimVersion < 2) { if (timeArrayOffsetMap.size() == 0) BuildTimeArrayOffsetMap(firstTrig); photonIdx += timeArrayOffsetMap.at(digiHit->GetTubeId()); } - + // Get the CherenkovHitTime objects themselves, which contain the primary parent IDs auto* theHitTimeObject = (WCSimRootCherenkovHitTime*)(firstTrig->GetCherenkovHitTimes()->At(photonIdx)); if (theHitTimeObject == nullptr) { logmessage = "LoadWCSim::GetHitParentIDs: HitTime object is NULL!!"; Log(logmessage, v_error, verbosity); + isDarknoise = false; } else { - parentIDs.push_back(theHitTimeObject->GetPrimaryParentID()); + int primaryParentID = theHitTimeObject->GetPrimaryParentID(); + parentIDs.push_back(primaryParentID); directParentIDs.push_back(theHitTimeObject->GetDirectParentID()); + if (primaryParentID != -1) isDarknoise = false; } - - }// end loop over photons - return std::make_pair(parentIDs, directParentIDs); + + }// end loop over photons + return std::make_tuple(parentIDs, directParentIDs, isDarknoise); } //////////////////////////////////////////////////////////////////////////////// // Get the index within the the MCParticle vector of the primaries that produced this digi. hit std::pair, std::vector> LoadWCSim::GetHitParentIdxs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig) { - std::pair, std::vector> bothIDs = GetHitParentIDs(digiHit, firstTrig); - std::vector parentIDs = bothIDs.first; - std::vector directParentIDs = bothIDs.second; + auto bothIDs = GetHitParentIDs(digiHit, firstTrig); + std::vector parentIDs = std::get<0>(bothIDs); + std::vector directParentIDs = std::get<1>(bothIDs); std::vector parentIdxs; std::vector directParentIdxs; diff --git a/UserTools/LoadWCSim/LoadWCSim.h b/UserTools/LoadWCSim/LoadWCSim.h index 944fd516b..752411092 100644 --- a/UserTools/LoadWCSim/LoadWCSim.h +++ b/UserTools/LoadWCSim/LoadWCSim.h @@ -4,6 +4,7 @@ #include #include +#include //#include // for boost::algorithm::to_lower(string) #include "Tool.h" @@ -121,8 +122,9 @@ class LoadWCSim: public Tool { // Each MCHit will contain the idx of it's parent MCParticle's // position within the MCParticles vector + //********** Note:- ParentID serves the purpose of the PrimaryParentID in the LoadWCSim tool, I did't want to break other things by just changing the variable name. ********** // std::map* trackid_to_mcparticleindex = nullptr; - std::pair, std::vector> GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig); + std::tuple, std::vector, bool> GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig); std::pair, std::vector> GetHitParentIdxs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig); std::map timeArrayOffsetMap; From 176f80af19d5597c207ae5d361b7fc392be64479 Mon Sep 17 00:00:00 2001 From: Dhavalkumar Ajana Date: Thu, 17 Sep 2026 18:09:17 -0500 Subject: [PATCH 4/6] Propagate the WCSim interaction mode through to each PMT hit LoadWCSim now records WCSimRootTrigger::GetMode() for every trigger in the entry and publishes the result to ANNIEEvent as "WCSimInteractionModes". BackTracker picks that vector up in LoadFromStores and, for each MC PMT hit, maps the hit's direct parent back to its MCTriggerNum to look up the interaction mode of the trigger that produced it. The result is published as "MCHitToInteractionMode" (-999 for hits taken from the pre-clock-tick path, -9999 when the mode cannot be resolved). ANNIEEventTreeMaker exposes this as the trueWCSimMode branch alongside the existing DirectParent hit-level branches, so the true GENIE/WCSim interaction channel behind every hit is available in the ntuple. --- .../ANNIEEventTreeMaker.cpp | 24 +++++++++++++++++++ .../ANNIEEventTreeMaker/ANNIEEventTreeMaker.h | 2 ++ UserTools/BackTracker/BackTracker.cpp | 24 ++++++++++++++++++- UserTools/BackTracker/BackTracker.h | 5 ++++ UserTools/LoadWCSim/LoadWCSim.cpp | 7 +++++- 5 files changed, 60 insertions(+), 2 deletions(-) diff --git a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp index f3f26cd61..a839b7d17 100644 --- a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp +++ b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp @@ -193,6 +193,7 @@ bool ANNIEEventTreeMaker::Initialise(std::string configfile, DataModel &data) fANNIETree->Branch("DirectParent_NeutronParentTrackID", &fDirectParent_NeutronParentTrackID); fANNIETree->Branch("DirectParent_NeutronParentPDG", &fDirectParent_NeutronParentPDG); fANNIETree->Branch("DirectParent_IsDarknoise", &fDirectParent_IsDarknoise); + fANNIETree->Branch("DirectParent_InteractionMode", &fDirectParent_InteractionMode); } if (SiPMPulseInfo_fill) @@ -484,6 +485,7 @@ bool ANNIEEventTreeMaker::Initialise(std::string configfile, DataModel &data) fANNIETree->Branch("trueKPlusCher", &fTrueKPlusCher, "trueKPlusCher/I"); fANNIETree->Branch("trueKMinus", &fTrueKMinus, "trueKMinus/I"); fANNIETree->Branch("trueKMinusCher", &fTrueKMinusCher, "trueKMinusCher/I"); + fANNIETree->Branch("trueWCSimMode", &fTrueWCSimMode, "trueWCSimMode/I"); } // Reconstructed variables after full Muon Reco Analysis @@ -795,6 +797,7 @@ void ANNIEEventTreeMaker::ResetVariables() fDirectParent_NeutronParentTrackID.clear(); fDirectParent_NeutronParentPDG.clear(); fDirectParent_IsDarknoise.clear(); + fDirectParent_InteractionMode.clear(); // SiPMPulse Info fSiPM1NPulses = 0; @@ -1038,6 +1041,7 @@ void ANNIEEventTreeMaker::ResetVariables() fTrueKPlusCher = -9999; fTrueKMinus = -9999; fTrueKMinusCher = -9999; + fTrueWCSimMode = -9999; // TankReco_fill fRecoVtxX = -9999; @@ -1456,6 +1460,8 @@ void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ bool got_neutronAncestorClass = m_data->Stores["ANNIEEvent"]->Get("MCHitToNeutronAncestorClass", fMCHitToNeutronAncestorClass); bool got_neutronParent = m_data->Stores["ANNIEEvent"]->Get("MCHitToNeutronParent", fMCHitToNeutronParent); bool got_isDarknoise = m_data->Stores["ANNIEEvent"]->Get("MCHitToIsDarknoise", fMCHitToIsDarknoise); + std::map>* fMCHitToInteractionMode = nullptr; + bool got_interactionMode = m_data->Stores["ANNIEEvent"]->Get("MCHitToInteractionMode", fMCHitToInteractionMode); for (auto const& apair : *fMCHitToDirectParents) { unsigned long pmtID = apair.first; @@ -1532,6 +1538,16 @@ void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ } } fDirectParent_IsDarknoise.push_back(isDarknoise); + + int interactionMode = -9999; + if (got_interactionMode && fMCHitToInteractionMode) { + auto pmtIt = fMCHitToInteractionMode->find(pmtID); + if (pmtIt != fMCHitToInteractionMode->end()) { + auto modeIt = pmtIt->second.find(hitTime); + if (modeIt != pmtIt->second.end()) interactionMode = modeIt->second; + } + } + fDirectParent_InteractionMode.push_back(interactionMode); } } return; @@ -2460,6 +2476,14 @@ bool ANNIEEventTreeMaker::FillMCTruthInfo() fiMCTriggerNum = (int)fMCTriggerNum; + { + std::vector wcSimModes; + if (m_data->Stores.at("ANNIEEvent")->Get("WCSimInteractionModes", wcSimModes) && !wcSimModes.empty()) { + int trigIdx = fiMCTriggerNum; + fTrueWCSimMode = (trigIdx >= 0 && trigIdx < (int)wcSimModes.size()) ? wcSimModes[trigIdx] : wcSimModes[0]; + } + } + std::map> MCNeutCap; bool get_neutcap = m_data->Stores.at("ANNIEEvent")->Get("MCNeutCap", MCNeutCap); if (!get_neutcap) diff --git a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h index 056294b35..9fa6a420a 100644 --- a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h +++ b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h @@ -231,6 +231,7 @@ class ANNIEEventTreeMaker : public Tool std::vector fDirectParent_NeutronParentTrackID; // direct parent track ID of the stored neutron ancestor std::vector fDirectParent_NeutronParentPDG; // direct parent PDG of the stored neutron ancestor std::vector fDirectParent_IsDarknoise; // 1 if the pulse is pure dark noise, 0 otherwise + std::vector fDirectParent_InteractionMode; // WCSim Nuance mode per hit; -999 darknoise; -9999 unavailable // SiPMPulseInfo_fill int fSiPM1NPulses; @@ -481,6 +482,7 @@ class ANNIEEventTreeMaker : public Tool int fTrueKPlusCher; int fTrueKMinus; int fTrueKMinusCher; + int fTrueWCSimMode = -9999; // raw WCSim Nuance interaction mode from GetMode() // TankReco_fill double fRecoVtxX; diff --git a/UserTools/BackTracker/BackTracker.cpp b/UserTools/BackTracker/BackTracker.cpp index ec2da4a51..799dc472c 100644 --- a/UserTools/BackTracker/BackTracker.cpp +++ b/UserTools/BackTracker/BackTracker.cpp @@ -45,7 +45,8 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){ fMCHitToNeutronAncestorClass = new std::map>; fMCHitToNeutronParent = new std::map>>; fMCHitToIsDarknoise = new std::map>; - + fMCHitToInteractionMode = new std::map>; + return true; } @@ -64,6 +65,7 @@ bool BackTracker::Execute() fMCHitToNeutronAncestorClass->clear(); fMCHitToNeutronParent->clear(); fMCHitToIsDarknoise->clear(); + fMCHitToInteractionMode->clear(); fParticleToTankTotalCharge.clear(); @@ -103,6 +105,7 @@ bool BackTracker::Execute() m_data->Stores.at("ANNIEEvent")->Set("MCHitToNeutronAncestorClass", fMCHitToNeutronAncestorClass); m_data->Stores.at("ANNIEEvent")->Set("MCHitToNeutronParent", fMCHitToNeutronParent); //It stores m_data->Stores.at("ANNIEEvent")->Set("MCHitToIsDarknoise", fMCHitToIsDarknoise); + m_data->Stores.at("ANNIEEvent")->Set("MCHitToInteractionMode", fMCHitToInteractionMode); return true; } @@ -283,6 +286,7 @@ void BackTracker::FindNeutronAncestors() { (*fMCHitToNeutronAncestor)[pmtID][hitTime] = std::make_pair(neutronAncestorId, neutronAncestorPdg); (*fMCHitToNeutronAncestorClass)[pmtID][hitTime] = neutronAncestorClass; (*fMCHitToNeutronParent)[pmtID][hitTime] = std::make_pair(neutronParentTrackId, neutronParentPdg); + (*fMCHitToInteractionMode)[pmtID][hitTime] = -999; continue; } @@ -330,6 +334,17 @@ void BackTracker::FindNeutronAncestors() { (*fMCHitToNeutronAncestorClass)[pmtID][hitTime] = neutronAncestorClass; (*fMCHitToNeutronParent)[pmtID][hitTime] = std::make_pair(neutronParentTrackId, neutronParentPdg); + int interactionMode = -9999; + if (!directParents.empty() && !fWCSimInteractionModes.empty()) { + auto particleIt = fMCParticleIndexMap->find(directParents[0]); + if (particleIt != fMCParticleIndexMap->end()) { + int trigNum = fMCParticles->at(particleIt->second).GetMCTriggerNum(); + if (trigNum >= 0 && trigNum < (int)fWCSimInteractionModes.size()) + interactionMode = fWCSimInteractionModes[trigNum]; + } + } + (*fMCHitToInteractionMode)[pmtID][hitTime] = interactionMode; + } } @@ -397,5 +412,12 @@ bool BackTracker::LoadFromStores() } } + fWCSimInteractionModes.clear(); + bool gotModes = m_data->Stores.at("ANNIEEvent")->Get("WCSimInteractionModes", fWCSimInteractionModes); + if (!gotModes) { + logmessage = "BackTracker: WCSimInteractionModes not in ANNIEEvent; interaction mode will be -9999 for all hits."; + Log(logmessage, v_warning, verbosity); + } + return true; } diff --git a/UserTools/BackTracker/BackTracker.h b/UserTools/BackTracker/BackTracker.h index 3473e5b05..1ae2aae7a 100644 --- a/UserTools/BackTracker/BackTracker.h +++ b/UserTools/BackTracker/BackTracker.h @@ -82,6 +82,11 @@ class BackTracker: public Tool { std::map>> *fMCHitToNeutronParent = nullptr; // PMT ID -> reco hit time -> true if pulse is pure dark noise std::map> *fMCHitToIsDarknoise = nullptr; + // PMT ID -> reco hit time -> WCSim Nuance interaction mode + // -999 = dark noise (no physics parent) + // -9999 = mode unavailable (LoadWCSim not run or BackTracker skipped) + std::map> *fMCHitToInteractionMode = nullptr; + std::vector fWCSimInteractionModes; // one entry per trigger, from LoadWCSim bool fDirectParentClockTickMatching = true; uint16_t fPMTSimPrewindowTicks = 10; diff --git a/UserTools/LoadWCSim/LoadWCSim.cpp b/UserTools/LoadWCSim/LoadWCSim.cpp index b7a998b42..2a07a5f89 100644 --- a/UserTools/LoadWCSim/LoadWCSim.cpp +++ b/UserTools/LoadWCSim/LoadWCSim.cpp @@ -495,7 +495,10 @@ bool LoadWCSim::Execute() int nMRDTriggers = WCSimEntry->wcsimrootevent_mrd->GetNumberOfEvents(); int nVetoTriggers = WCSimEntry->wcsimrootevent_facc->GetNumberOfEvents(); - + + std::vector WCSimInteractionModes; + WCSimInteractionModes.reserve(trigsInEntry); + // Loop over over the triggers // ============================================= while (MCTriggerNum < MaxEventNr) { @@ -504,6 +507,7 @@ bool LoadWCSim::Execute() Log(logmessage, v_message, verbosity); WCSimRootTrigger* aTrigTank = WCSimEntry->wcsimrootevent->GetTrigger(MCTriggerNum); + WCSimInteractionModes.push_back(aTrigTank->GetMode()); WCSimRootTrigger* aTrigMRD = ( (MCTriggerNum < nMRDTriggers) ? WCSimEntry->wcsimrootevent_mrd->GetTrigger(MCTriggerNum) : nullptr ); @@ -639,6 +643,7 @@ bool LoadWCSim::Execute() m_data->Stores.at("ANNIEEvent")->Set("BeamStatus", beamstat); m_data->Stores.at("ANNIEEvent")->Set("MCNeutCap", MCNeutCap); m_data->Stores.at("ANNIEEvent")->Set("MCNeutCapGammas", MCNeutCapGammas); + m_data->Stores.at("ANNIEEvent")->Set("WCSimInteractionModes", WCSimInteractionModes); m_data->CStore.Set("NumTriggersThisMCEvt", trigsInEntry); // auxilliary information about MC Truth particles From 405c878fd20ec5151f615c5b7217ff0abb53e837 Mon Sep 17 00:00:00 2001 From: Dhavalkumar Ajana Date: Thu, 17 Sep 2026 18:09:31 -0500 Subject: [PATCH 5/6] Trace the ancestry chain behind each tank PMT hit Extends BackTracker's hit-level backtracking in three ways and surfaces all of it through ANNIEEventTreeMaker: * Dark-noise tagging: a hit whose direct parents are all -1 was treated as dark noise unconditionally. In older WCSim files (e.g. the AmBe wcsim_0_999 samples) -1 is a legitimate track ID, so if -1 is present in the track map the hit is kept as real physics instead of being discarded as noise. * Immediate ancestor: for every MC PMT hit the background particle directly behind the detected photon is resolved and published as "MCHitToImmediateAncestor" / "MCHitToImmediateAncestorClass". This is what is needed to separate genuine neutron-capture light from the background particles landing in the same time window. * Full lineage: rather than stopping at the immediate ancestor, the parent links are walked all the way back, and the resulting chain is published as "MCHitToLineage" and "MCHitToLineageStatus", with the end of the chain in "MCHitToRootAncestor". This makes it possible to attribute a tank optical photon to the primary interaction product it ultimately came from. ANNIEEventTreeMaker gains the matching per-hit branches -- DirectParent_ImmediateAncestor{TrackID,PDG,Class}, DirectParent_RootAncestor{TrackID,PDG} and DirectParent_Lineage{PDG,TrackID,Depth,Status} -- so the chain is written out with the rest of the DirectParent hit information. --- .../ANNIEEventTreeMaker.cpp | 111 ++++++++++++++ .../ANNIEEventTreeMaker/ANNIEEventTreeMaker.h | 33 +++++ UserTools/BackTracker/BackTracker.cpp | 135 ++++++++++++++++++ UserTools/BackTracker/BackTracker.h | 68 ++++++++- 4 files changed, 346 insertions(+), 1 deletion(-) diff --git a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp index a839b7d17..a12598993 100644 --- a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp +++ b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.cpp @@ -194,6 +194,18 @@ bool ANNIEEventTreeMaker::Initialise(std::string configfile, DataModel &data) fANNIETree->Branch("DirectParent_NeutronParentPDG", &fDirectParent_NeutronParentPDG); fANNIETree->Branch("DirectParent_IsDarknoise", &fDirectParent_IsDarknoise); fANNIETree->Branch("DirectParent_InteractionMode", &fDirectParent_InteractionMode); + fANNIETree->Branch("DirectParent_ImmediateAncestorTrackID", &fDirectParent_ImmediateAncestorTrackID); + fANNIETree->Branch("DirectParent_ImmediateAncestorPDG", &fDirectParent_ImmediateAncestorPDG); + fANNIETree->Branch("DirectParent_ImmediateAncestorClass", &fDirectParent_ImmediateAncestorClass); + // DISABLED: PrimaryAncestor + // fANNIETree->Branch("DirectParent_PrimaryAncestorTrackID", &fDirectParent_PrimaryAncestorTrackID); + // fANNIETree->Branch("DirectParent_PrimaryAncestorPDG", &fDirectParent_PrimaryAncestorPDG); + fANNIETree->Branch("DirectParent_RootAncestorTrackID", &fDirectParent_RootAncestorTrackID); + fANNIETree->Branch("DirectParent_RootAncestorPDG", &fDirectParent_RootAncestorPDG); + fANNIETree->Branch("DirectParent_LineagePDG", &fDirectParent_LineagePDG); + fANNIETree->Branch("DirectParent_LineageTrackID", &fDirectParent_LineageTrackID); + fANNIETree->Branch("DirectParent_LineageDepth", &fDirectParent_LineageDepth); + fANNIETree->Branch("DirectParent_LineageStatus", &fDirectParent_LineageStatus); } if (SiPMPulseInfo_fill) @@ -798,6 +810,18 @@ void ANNIEEventTreeMaker::ResetVariables() fDirectParent_NeutronParentPDG.clear(); fDirectParent_IsDarknoise.clear(); fDirectParent_InteractionMode.clear(); + fDirectParent_ImmediateAncestorTrackID.clear(); + fDirectParent_ImmediateAncestorPDG.clear(); + fDirectParent_ImmediateAncestorClass.clear(); + // DISABLED: PrimaryAncestor + // fDirectParent_PrimaryAncestorTrackID.clear(); + // fDirectParent_PrimaryAncestorPDG.clear(); + fDirectParent_RootAncestorTrackID.clear(); + fDirectParent_RootAncestorPDG.clear(); + fDirectParent_LineagePDG.clear(); + fDirectParent_LineageTrackID.clear(); + fDirectParent_LineageDepth.clear(); + fDirectParent_LineageStatus.clear(); // SiPMPulse Info fSiPM1NPulses = 0; @@ -1462,6 +1486,19 @@ void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ bool got_isDarknoise = m_data->Stores["ANNIEEvent"]->Get("MCHitToIsDarknoise", fMCHitToIsDarknoise); std::map>* fMCHitToInteractionMode = nullptr; bool got_interactionMode = m_data->Stores["ANNIEEvent"]->Get("MCHitToInteractionMode", fMCHitToInteractionMode); + std::map>>* fMCHitToImmediateAncestor = nullptr; + bool got_immediateAncestor = m_data->Stores["ANNIEEvent"]->Get("MCHitToImmediateAncestor", fMCHitToImmediateAncestor); + std::map>* fMCHitToImmediateAncestorClass = nullptr; + bool got_immediateAncestorClass = m_data->Stores["ANNIEEvent"]->Get("MCHitToImmediateAncestorClass", fMCHitToImmediateAncestorClass); + // DISABLED: PrimaryAncestor + // std::map>>* fMCHitToPrimaryAncestor = nullptr; + // bool got_primaryAncestor = m_data->Stores["ANNIEEvent"]->Get("MCHitToPrimaryAncestor", fMCHitToPrimaryAncestor); + std::map>>* fMCHitToRootAncestor = nullptr; + bool got_rootAncestor = m_data->Stores["ANNIEEvent"]->Get("MCHitToRootAncestor", fMCHitToRootAncestor); + std::map>>>* fMCHitToLineage = nullptr; + bool got_lineage = m_data->Stores["ANNIEEvent"]->Get("MCHitToLineage", fMCHitToLineage); + std::map>* fMCHitToLineageStatus = nullptr; + bool got_lineageStatus = m_data->Stores["ANNIEEvent"]->Get("MCHitToLineageStatus", fMCHitToLineageStatus); for (auto const& apair : *fMCHitToDirectParents) { unsigned long pmtID = apair.first; @@ -1548,6 +1585,80 @@ void ANNIEEventTreeMaker::LoadDirectParentIDsMCHits(){ } } fDirectParent_InteractionMode.push_back(interactionMode); + + int immediateAncestorTrackID = -5; + int immediateAncestorPDG = -5; + if (got_immediateAncestor && fMCHitToImmediateAncestor->find(pmtID) != fMCHitToImmediateAncestor->end()){ + auto const& pmtAncestors = fMCHitToImmediateAncestor->at(pmtID); + if (pmtAncestors.find(hitTime) != pmtAncestors.end()){ + auto const& ancestorPair = pmtAncestors.at(hitTime); + immediateAncestorTrackID = ancestorPair.first; + immediateAncestorPDG = ancestorPair.second; + } + } + fDirectParent_ImmediateAncestorTrackID.push_back(immediateAncestorTrackID); + fDirectParent_ImmediateAncestorPDG.push_back(immediateAncestorPDG); + + int immediateAncestorClass = -5; + if (got_immediateAncestorClass && fMCHitToImmediateAncestorClass->find(pmtID) != fMCHitToImmediateAncestorClass->end()){ + auto const& pmtClasses = fMCHitToImmediateAncestorClass->at(pmtID); + if (pmtClasses.find(hitTime) != pmtClasses.end()){ + immediateAncestorClass = pmtClasses.at(hitTime); + } + } + fDirectParent_ImmediateAncestorClass.push_back(immediateAncestorClass); + + // --------------------------------------------------------------------------------- + // DISABLED: PrimaryAncestor. Original fill, kept verbatim for future retrieval. + // int primaryAncestorTrackID = -5; + // int primaryAncestorPDG = -5; + // if (got_primaryAncestor && fMCHitToPrimaryAncestor->find(pmtID) != fMCHitToPrimaryAncestor->end()){ + // auto const& pmtPrimaries = fMCHitToPrimaryAncestor->at(pmtID); + // if (pmtPrimaries.find(hitTime) != pmtPrimaries.end()){ + // auto const& primaryPair = pmtPrimaries.at(hitTime); + // primaryAncestorTrackID = primaryPair.first; + // primaryAncestorPDG = primaryPair.second; + // } + // } + // fDirectParent_PrimaryAncestorTrackID.push_back(primaryAncestorTrackID); + // fDirectParent_PrimaryAncestorPDG.push_back(primaryAncestorPDG); + // --------------------------------------------------------------------------------- + + int rootAncestorTrackID = -5; + int rootAncestorPDG = -5; + if (got_rootAncestor && fMCHitToRootAncestor->find(pmtID) != fMCHitToRootAncestor->end()){ + auto const& pmtRoots = fMCHitToRootAncestor->at(pmtID); + if (pmtRoots.find(hitTime) != pmtRoots.end()){ + auto const& rootPair = pmtRoots.at(hitTime); + rootAncestorTrackID = rootPair.first; + rootAncestorPDG = rootPair.second; + } + } + fDirectParent_RootAncestorTrackID.push_back(rootAncestorTrackID); + fDirectParent_RootAncestorPDG.push_back(rootAncestorPDG); + + // Flatten this hit's lineage onto the concatenated arrays and record its depth, so the + // per-hit chain stays recoverable without a nested-vector branch. + int lineageDepth = 0; + if (got_lineage && fMCHitToLineage->find(pmtID) != fMCHitToLineage->end()){ + auto const& pmtLineages = fMCHitToLineage->at(pmtID); + if (pmtLineages.find(hitTime) != pmtLineages.end()){ + auto const& chain = pmtLineages.at(hitTime); + for (auto const& step : chain){ + fDirectParent_LineageTrackID.push_back(step.first); + fDirectParent_LineagePDG.push_back(step.second); + } + lineageDepth = (int)chain.size(); + } + } + fDirectParent_LineageDepth.push_back(lineageDepth); + + int lineageStatus = -5; + if (got_lineageStatus && fMCHitToLineageStatus->find(pmtID) != fMCHitToLineageStatus->end()){ + auto const& pmtStatus = fMCHitToLineageStatus->at(pmtID); + if (pmtStatus.find(hitTime) != pmtStatus.end()) lineageStatus = pmtStatus.at(hitTime); + } + fDirectParent_LineageStatus.push_back(lineageStatus); } } return; diff --git a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h index 9fa6a420a..7c9466730 100644 --- a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h +++ b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h @@ -232,6 +232,39 @@ class ANNIEEventTreeMaker : public Tool std::vector fDirectParent_NeutronParentPDG; // direct parent PDG of the stored neutron ancestor std::vector fDirectParent_IsDarknoise; // 1 if the pulse is pure dark noise, 0 otherwise std::vector fDirectParent_InteractionMode; // WCSim Nuance mode per hit; -999 darknoise; -9999 unavailable + // Immediate background particle: generalizes the neutron-only ancestry walk above to + // any species via a shallow (at most one-step) walk that skips a leading e-/e+ direct + // parent, since electrons/positrons are the ubiquitous last-step Cherenkov/ionization + // carriers in a water Cherenkov detector and aren't informative as "the background particle." + std::vector fDirectParent_ImmediateAncestorTrackID; // -5 if dark noise or untraced + std::vector fDirectParent_ImmediateAncestorPDG; // -5 if dark noise or untraced + std::vector fDirectParent_ImmediateAncestorClass; // 0 dark-noise,1 neutron,2 muon,3 charged pion,4 proton,5 photon,6 kaon,7 electron/positron,8 other,-5 untraced + // DISABLED: PrimaryAncestor (PrimaryParentID-based, never committed - kept for retrieval). + // Replaced by the RootAncestor fields below, which reach the same top-of-tree particle by + // walking DirectParentID, the same mechanism as the neutron scheme. + // std::vector fDirectParent_PrimaryAncestorTrackID; // -5 if dark noise or unresolved + // std::vector fDirectParent_PrimaryAncestorPDG; // -5 if dark noise or unresolved + + // Root ancestor: the particle at the TOP of the DirectParentID chain, i.e. the + // generator-level particle out of the neutrino interaction that this deposit descends + // from. Derived by walking DirectParentID to its end -- the same mechanism as the + // neutron scheme -- NOT from PrimaryParentID. Trust only where LineageStatus == 1. + std::vector fDirectParent_RootAncestorTrackID; // -5 if dark noise or unresolved + std::vector fDirectParent_RootAncestorPDG; // -5 if dark noise or unresolved + + // FULL lineage per hit, the species-general analogue of the neutron scheme's + // {NeutronAncestor, NeutronParent} pair but as a whole chain. Stored flattened + // (CSR-style) rather than as vector>, which would need a custom ROOT + // dictionary. To read hit i, take Depth[i] entries starting at sum(Depth[0..i-1]): + // chain_pdgs_of_hit_i = LineagePDG[offset : offset + LineageDepth[i]] + // Order is nearest-first: [0] the hit's direct parent (nearly always e-/e+), + // [1] its parent (e.g. a capture gamma), ... up to the generator primary. + std::vector fDirectParent_LineagePDG; // concatenated over all hits in the event + std::vector fDirectParent_LineageTrackID; // concatenated, parallel to LineagePDG + std::vector fDirectParent_LineageDepth; // per hit: number of entries contributed + // per hit: 1 chain reached a generator primary (complete); 0 truncated because WCSim + // did not save the next ancestor; -1 circular reference or depth cap; -5 dark noise + std::vector fDirectParent_LineageStatus; // SiPMPulseInfo_fill int fSiPM1NPulses; diff --git a/UserTools/BackTracker/BackTracker.cpp b/UserTools/BackTracker/BackTracker.cpp index 799dc472c..9edc747eb 100644 --- a/UserTools/BackTracker/BackTracker.cpp +++ b/UserTools/BackTracker/BackTracker.cpp @@ -46,6 +46,13 @@ bool BackTracker::Initialise(std::string configfile, DataModel &data){ fMCHitToNeutronParent = new std::map>>; fMCHitToIsDarknoise = new std::map>; fMCHitToInteractionMode = new std::map>; + fMCHitToImmediateAncestor = new std::map>>; + fMCHitToImmediateAncestorClass = new std::map>; + // DISABLED: PrimaryAncestor (PrimaryParentID-based; see BackTracker.h) + // fMCHitToPrimaryAncestor = new std::map>>; + fMCHitToRootAncestor = new std::map>>; + fMCHitToLineage = new std::map>>>; + fMCHitToLineageStatus = new std::map>; return true; } @@ -66,6 +73,13 @@ bool BackTracker::Execute() fMCHitToNeutronParent->clear(); fMCHitToIsDarknoise->clear(); fMCHitToInteractionMode->clear(); + fMCHitToImmediateAncestor->clear(); + fMCHitToImmediateAncestorClass->clear(); + // DISABLED: PrimaryAncestor + // fMCHitToPrimaryAncestor->clear(); + fMCHitToRootAncestor->clear(); + fMCHitToLineage->clear(); + fMCHitToLineageStatus->clear(); fParticleToTankTotalCharge.clear(); @@ -106,6 +120,13 @@ bool BackTracker::Execute() m_data->Stores.at("ANNIEEvent")->Set("MCHitToNeutronParent", fMCHitToNeutronParent); //It stores m_data->Stores.at("ANNIEEvent")->Set("MCHitToIsDarknoise", fMCHitToIsDarknoise); m_data->Stores.at("ANNIEEvent")->Set("MCHitToInteractionMode", fMCHitToInteractionMode); + m_data->Stores.at("ANNIEEvent")->Set("MCHitToImmediateAncestor", fMCHitToImmediateAncestor); + m_data->Stores.at("ANNIEEvent")->Set("MCHitToImmediateAncestorClass", fMCHitToImmediateAncestorClass); + // DISABLED: PrimaryAncestor + // m_data->Stores.at("ANNIEEvent")->Set("MCHitToPrimaryAncestor", fMCHitToPrimaryAncestor); + m_data->Stores.at("ANNIEEvent")->Set("MCHitToRootAncestor", fMCHitToRootAncestor); + m_data->Stores.at("ANNIEEvent")->Set("MCHitToLineage", fMCHitToLineage); + m_data->Stores.at("ANNIEEvent")->Set("MCHitToLineageStatus", fMCHitToLineageStatus); return true; } @@ -267,6 +288,11 @@ void BackTracker::FindNeutronAncestors() { // the pulse peak_time which is what fMCHitToDirectParents uses. bool isDarknoise = std::all_of(directParents.begin(), directParents.end(), [](int id) { return id == -1; }); + // Temporary fix for older WCSim files (e.g. AmBe wcsim_0_999) where -1 is + // stored as the neutron's actual track ID rather than a dark-noise sentinel. + // If -1 exists in trackMap it is a real physics particle, not noise. + if (isDarknoise && trackMap.find(-1) != trackMap.end()) + isDarknoise = false; (*fMCHitToIsDarknoise)[pmtID][hitTime] = isDarknoise; int neutronAncestorId = -5; @@ -287,6 +313,13 @@ void BackTracker::FindNeutronAncestors() { (*fMCHitToNeutronAncestorClass)[pmtID][hitTime] = neutronAncestorClass; (*fMCHitToNeutronParent)[pmtID][hitTime] = std::make_pair(neutronParentTrackId, neutronParentPdg); (*fMCHitToInteractionMode)[pmtID][hitTime] = -999; + (*fMCHitToImmediateAncestor)[pmtID][hitTime] = std::make_pair(-5, -5); + (*fMCHitToImmediateAncestorClass)[pmtID][hitTime] = 0; + // DISABLED: PrimaryAncestor + // (*fMCHitToPrimaryAncestor)[pmtID][hitTime] = std::make_pair(-5, -5); + (*fMCHitToRootAncestor)[pmtID][hitTime] = std::make_pair(-5, -5); + (*fMCHitToLineage)[pmtID][hitTime] = std::vector>(); + (*fMCHitToLineageStatus)[pmtID][hitTime] = -5; continue; } @@ -334,6 +367,94 @@ void BackTracker::FindNeutronAncestors() { (*fMCHitToNeutronAncestorClass)[pmtID][hitTime] = neutronAncestorClass; (*fMCHitToNeutronParent)[pmtID][hitTime] = std::make_pair(neutronParentTrackId, neutronParentPdg); + // Immediate background particle: walks the same DirectParentID chain as the + // neutron search above, but stops at the nearest ancestor that actually names a + // species instead of at PDG 2112. Electrons/positrons are skipped over, because + // in a water Cherenkov detector essentially every hit is mediated by an e-/e+ at + // the last step (Cherenkov radiation, Compton scattering, pair production, the + // photoelectric effect), so reporting "electron" carries no information. + // + // The skip is repeated, not single-step: an EM shower can put several e-/e+ + // generations in a row on the chain, and one step back would stop short of the + // species that seeded the shower. If the whole recorded lineage is e-/e+ we + // report the last one reached (class 7), and if nothing on the chain resolves at + // all the hit stays untraced (-5). + int immediateAncestorId = -5; + int immediateAncestorPdg = -5; + { + std::set skipVisited; // guards against circular DirectParentID references + int walkId = startParentID; + auto walkIt = trackMap.find(walkId); + while (walkIt != trackMap.end() && skipVisited.count(walkId) == 0) { + skipVisited.insert(walkId); + int walkPdg = walkIt->second.second; + immediateAncestorId = walkId; + immediateAncestorPdg = walkPdg; + if (walkPdg != 11 && walkPdg != -11) break; // found a species-naming ancestor + walkId = walkIt->second.first; // still a lepton: keep going up + walkIt = trackMap.find(walkId); + } + } + (*fMCHitToImmediateAncestor)[pmtID][hitTime] = std::make_pair(immediateAncestorId, immediateAncestorPdg); + (*fMCHitToImmediateAncestorClass)[pmtID][hitTime] = ClassifyBackgroundPDG(immediateAncestorPdg); + + // --------------------------------------------------------------------------------- + // DISABLED: PrimaryAncestor. Original PrimaryParentID-based lookup, kept verbatim for + // future retrieval (it was never committed, so git cannot recover it). Superseded by + // the DirectParentID walk below, which reaches the same top-of-tree particle using the + // same mechanism as the neutron scheme. + // + // int primaryAncestorId = -5; + // int primaryAncestorPdg = -5; + // { + // auto startIt = fMCParticleIndexMap->find(startParentID); + // if (startIt != fMCParticleIndexMap->end()) { + // primaryAncestorId = fMCParticles->at(startIt->second).GetPrimaryParentID(); + // auto primaryIt = fMCParticleIndexMap->find(primaryAncestorId); + // if (primaryIt != fMCParticleIndexMap->end()) + // primaryAncestorPdg = fMCParticles->at(primaryIt->second).GetPdgCode(); + // } + // } + // (*fMCHitToPrimaryAncestor)[pmtID][hitTime] = std::make_pair(primaryAncestorId, primaryAncestorPdg); + // --------------------------------------------------------------------------------- + + // Full lineage of this optical photon, walking the SAME DirectParentID chain the + // neutron search above uses -- no PrimaryParentID, so every branch here is derived + // by one consistent mechanism. Ordered nearest-first: [0] is the hit's direct parent + // (almost always an e-/e+), [1] its parent (e.g. a capture gamma), and so on up to + // the generator primary. This is what lets a class-(-5) hit be read as a chain, + // "e- <- gamma <- neutron <- proton", rather than as a single label. + std::vector> lineage; + int lineageStatus = 0; // 0 = truncated at an unsaved track, until proven otherwise + { + std::set lineageVisited; + int walkId = startParentID; + while (true) { + if ((int)lineage.size() >= kMaxLineageDepth) { lineageStatus = -1; break; } + if (lineageVisited.count(walkId)) { + std::cerr << "WARNING: circular DirectParentID in lineage at trackID " + << walkId << std::endl; + lineageStatus = -1; + break; + } + auto walkIt = trackMap.find(walkId); + if (walkIt == trackMap.end()) break; // WCSim never saved this track + lineageVisited.insert(walkId); + lineage.push_back(std::make_pair(walkId, walkIt->second.second)); + int nextId = walkIt->second.first; + if (nextId <= 0) { lineageStatus = 1; break; } // reached a generator primary + walkId = nextId; + } + } + (*fMCHitToLineage)[pmtID][hitTime] = lineage; + (*fMCHitToLineageStatus)[pmtID][hitTime] = lineage.empty() ? 0 : lineageStatus; + + // Root of that chain: the generator-level particle the deposit descends from. + // Trustworthy only when lineageStatus == 1. + int rootAncestorId = lineage.empty() ? -5 : lineage.back().first; + int rootAncestorPdg = lineage.empty() ? -5 : lineage.back().second; + (*fMCHitToRootAncestor)[pmtID][hitTime] = std::make_pair(rootAncestorId, rootAncestorPdg); + int interactionMode = -9999; if (!directParents.empty() && !fWCSimInteractionModes.empty()) { auto particleIt = fMCParticleIndexMap->find(directParents[0]); @@ -353,6 +474,20 @@ void BackTracker::FindNeutronAncestors() { << " PMTs with neutron ancestors." << std::endl; } +int BackTracker::ClassifyBackgroundPDG(int pdg) const { + switch (pdg) { + case 2112: return 1; // neutron + case 13: case -13: return 2; // muon + case 211: case -211: return 3; // charged pion + case 2212: return 4; // proton + case 22: return 5; // photon + case 321: case -321: case 311: case -311: return 6; // kaon + case 11: case -11: return 7; // electron/positron (skip landed on another lepton) + case -5: return -5; // untraced + default: return 8; // other identified species + } +} + bool BackTracker::LoadFromStores() { // Grab the stuff we need from the stores diff --git a/UserTools/BackTracker/BackTracker.h b/UserTools/BackTracker/BackTracker.h index 1ae2aae7a..340cc0aec 100644 --- a/UserTools/BackTracker/BackTracker.h +++ b/UserTools/BackTracker/BackTracker.h @@ -35,9 +35,13 @@ class BackTracker: public Tool { void MatchMCParticle(std::vector const &mchits, int &prtId, int &prtPdg, double &eff, double &pur, double &totalCharge); ///< The meat and potatoes void DirectParentsFromClockTickWindows(); void FindNeutronAncestors(); - + private: + // Classifies the PDG of an immediate-background-particle lookup (see + // fMCHitToImmediateAncestorClass below for the class codes). + int ClassifyBackgroundPDG(int pdg) const; + // Things we need to pull out of the store std::map> *fMCHitsMap = nullptr; ///< All of the MCHits keyed by channel number std::map> *fClusterMapMC = nullptr; ///< Clusters that we will be linking MCParticles to @@ -88,6 +92,68 @@ class BackTracker: public Tool { std::map> *fMCHitToInteractionMode = nullptr; std::vector fWCSimInteractionModes; // one entry per trigger, from LoadWCSim + // PMT ID -> reco hit time -> (immediate background ancestor trackID, PDG) + // Generalization of the neutron-only walk above to any species: walks up the + // DirectParentID chain from the hit's direct parent. If that direct parent is an + // e-/e+ (PDG +-11) -- in a water Cherenkov detector essentially every hit is + // mediated by an electron/positron at the last step (Cherenkov radiation, + // Compton/pair-production, photoelectric effect), so naming it as "the + // background particle" carries no information -- look exactly one step further + // up to that electron's own direct parent and report that instead. + // -5 = dark noise, or ancestor track ID not present in MCParticles (untraced) + std::map>> *fMCHitToImmediateAncestor = nullptr; + // PMT ID -> reco hit time -> classification of fMCHitToImmediateAncestor's PDG + // 0 dark noise + // 1 neutron (2112) + // 2 muon (+-13) + // 3 charged pion (+-211) + // 4 proton (2212) + // 5 photon (22) -- e.g. non-capture gamma, such as from a pi0 decay or bremsstrahlung + // 6 kaon (+-321, 311, -311) + // 7 electron/positron (+-11) -- the one-step skip landed on another lepton + // 8 other identified species (not covered above) + // -5 untraced -- ancestor track ID not found in MCParticles (not saved by WCSim) + std::map> *fMCHitToImmediateAncestorClass = nullptr; + + // --------------------------------------------------------------------------------- + // DISABLED (kept for future retrieval, never committed anywhere): PrimaryParentID-based + // ancestor. This read MCParticle::GetPrimaryParentID() instead of walking DirectParentID. + // Turned off deliberately because PrimaryParentID is a separate WCSim mechanism (inherited + // top-down through WCSimTrackInformation), so it answers a question derived differently + // from the neutron scheme and the rest of these branches. fMCHitToRootAncestor below is + // the walk-derived replacement. To re-enable, uncomment this plus the four blocks marked + // "DISABLED: PrimaryAncestor" in BackTracker.cpp and ANNIEEventTreeMaker.{h,cpp}. + // + // std::map>> *fMCHitToPrimaryAncestor = nullptr; + // --------------------------------------------------------------------------------- + + // PMT ID -> reco hit time -> (root ancestor trackID, PDG) + // The particle at the TOP of the DirectParentID chain below -- i.e. the generator-level + // particle out of the neutrino interaction that this charge deposit descends from. + // Deliberately derived by walking DirectParentID to its end, exactly like the neutron + // scheme, and NOT from MCParticle::GetPrimaryParentID(): PrimaryParentID is a separate + // WCSim mechanism (inherited top-down through WCSimTrackInformation) and would answer a + // question derived differently from the rest of these branches. + // Only meaningful when the corresponding LineageStatus == 1 (chain reached a primary). + // -5 = dark noise, or nothing on the chain resolved + std::map>> *fMCHitToRootAncestor = nullptr; + + // PMT ID -> reco hit time -> FULL lineage of the optical photon, ordered nearest-first: + // [0] = the hit's direct parent (in a water Cherenkov detector nearly always an e-/e+), + // [1] = that particle's parent (e.g. the capture gamma), ... up to the generator primary. + // Each entry is (trackID, PDG). This is the species-general analogue of the neutron + // scheme's {NeutronAncestor, NeutronParent} pair, but as a whole chain rather than a + // single hop, so a class-(-5) hit can be read as "e- <- gamma <- neutron <- proton". + std::map>>> *fMCHitToLineage = nullptr; + // PMT ID -> reco hit time -> how the walk terminated: + // 1 reached a generator-level primary (DirectParentID == 0) -- lineage is complete + // 0 truncated: the next ancestor is not in MCParticles (WCSim did not save it) + // -1 aborted on a circular DirectParentID reference or the depth cap + // -5 dark noise: no lineage + std::map> *fMCHitToLineageStatus = nullptr; + // Safety cap on the ancestry walk; real chains are only a few generations deep. + static const int kMaxLineageDepth = 64; + bool fDirectParentClockTickMatching = true; uint16_t fPMTSimPrewindowTicks = 10; uint16_t fPMTSimReadoutWindowTicks = 35; From fec3e2242fe78b351a9e6b53bab97286cc61e517 Mon Sep 17 00:00:00 2001 From: Dhavalkumar Ajana Date: Thu, 17 Sep 2026 18:10:11 -0500 Subject: [PATCH 6/6] Reduce debug verbosity across the MC tools These tools printed on every event regardless of the configured verbosity, which made long MC productions unreadable and dominated the log files. * BackTracker: gate the FindNeutronAncestors summary print on verbosity > 1. * LoadWCSim: drop the unconditional per-event "Event Number" print and the dump of every saved track ID in trackid_to_mcparticleindex. * PMTWaveformSim: drop the unconditional "finished looping over MCHits" print. * DigitBuilder: an absent MCLAPPDHits store entry is expected when running without LAPPD MC, so log it at v_debug rather than v_error. * MCParticleProperties: track-length and tank entry/exit intercept warnings fire routinely for particles that clip the tank, so log them at v_debug rather than v_error. No behaviour changes -- only logging levels and removed prints. --- UserTools/BackTracker/BackTracker.cpp | 8 +++++--- UserTools/DigitBuilder/DigitBuilder.cpp | 2 +- UserTools/LoadWCSim/LoadWCSim.cpp | 14 +++++++------- .../MCParticleProperties/MCParticleProperties.cpp | 6 +++--- UserTools/PMTWaveformSim/PMTWaveformSim.cpp | 2 +- 5 files changed, 17 insertions(+), 15 deletions(-) diff --git a/UserTools/BackTracker/BackTracker.cpp b/UserTools/BackTracker/BackTracker.cpp index 9edc747eb..f164caa5c 100644 --- a/UserTools/BackTracker/BackTracker.cpp +++ b/UserTools/BackTracker/BackTracker.cpp @@ -469,9 +469,11 @@ void BackTracker::FindNeutronAncestors() { } } - std::cout << "BackTracker::FindNeutronAncestors: found " - << fMCHitToNeutronAncestor->size() - << " PMTs with neutron ancestors." << std::endl; + if (verbosity > 1) { + std::cout << "BackTracker::FindNeutronAncestors: found " + << fMCHitToNeutronAncestor->size() + << " PMTs with neutron ancestors." << std::endl; + } } int BackTracker::ClassifyBackgroundPDG(int pdg) const { diff --git a/UserTools/DigitBuilder/DigitBuilder.cpp b/UserTools/DigitBuilder/DigitBuilder.cpp index 62f591f66..7f570d637 100644 --- a/UserTools/DigitBuilder/DigitBuilder.cpp +++ b/UserTools/DigitBuilder/DigitBuilder.cpp @@ -161,7 +161,7 @@ bool DigitBuilder::Execute(){ } auto get_mclappdhits = m_data->Stores.at("ANNIEEvent")->Get("MCLAPPDHits",fMCLAPPDHits); if(!get_mclappdhits){ - Log("DigitBuilder Tool: Error retrieving MCLAPPDHits from ANNIEEvent!",v_error,verbosity); + Log("DigitBuilder Tool: Error retrieving MCLAPPDHits from ANNIEEvent!",v_debug,verbosity); return false; } } else { diff --git a/UserTools/LoadWCSim/LoadWCSim.cpp b/UserTools/LoadWCSim/LoadWCSim.cpp index 2a07a5f89..cb69c7e5e 100644 --- a/UserTools/LoadWCSim/LoadWCSim.cpp +++ b/UserTools/LoadWCSim/LoadWCSim.cpp @@ -1111,7 +1111,7 @@ void LoadWCSim::LoadMCParticles(WCSimRootTrigger* firstTrig) logmessage += " tracks from trigger # " + std::to_string(trigIdx); Log(logmessage, v_message, verbosity); - std::cout<< "Event Number: " << aTrigTank->GetHeader()->GetEvtNum()<< std::endl; + // std::cout<< "Event Number: " << aTrigTank->GetHeader()->GetEvtNum()<< std::endl; for (int trackIdx = 0; trackIdx < aTrigTank->GetNtrack(); trackIdx++) { logmessage = "LoadWCSim::LoadMCParticles: Getting WCSim track # " + std::to_string(trackIdx); Log(logmessage, v_message, verbosity); @@ -1215,12 +1215,12 @@ void LoadWCSim::LoadMCParticles(WCSimRootTrigger* firstTrig) Log(logmessage, v_debug, verbosity); } // end loop over events - std::cout << "DEBUG: All saved Track IDs in trackid_to_mcparticleindex: "; - for (auto const& pair : *trackid_to_mcparticleindex) { - std::cout << pair.first << " "; - } - -std::cout << std::endl; + // std::cout << "DEBUG: All saved Track IDs in trackid_to_mcparticleindex: "; + // for (auto const& pair : *trackid_to_mcparticleindex) { + // std::cout << pair.first << " "; + // } + // + // std::cout << std::endl; }// endif MCTriggerNum == 0 else { // if MCTrigger > 0 we need to update all the particle times diff --git a/UserTools/MCParticleProperties/MCParticleProperties.cpp b/UserTools/MCParticleProperties/MCParticleProperties.cpp index ccbb4d9d5..5105505be 100644 --- a/UserTools/MCParticleProperties/MCParticleProperties.cpp +++ b/UserTools/MCParticleProperties/MCParticleProperties.cpp @@ -373,7 +373,7 @@ bool MCParticleProperties::Execute(){ cout<<"c.f. max possible tank track length is "< maxtanktracklength){ - Log("MCParticleProperties Tool: Track length is impossibly long!",v_error,verbosity); + Log("MCParticleProperties Tool: Track length is impossibly long!",v_debug,verbosity); //return false; } if(atracklengthintank > differencevector.Mag()){ @@ -762,7 +762,7 @@ bool MCParticleProperties::CheckTankIntercepts( Position startvertex, Position s } } // else track did not start outside tank x bounds: no wall entry - if(!entryfound) Log("MCParticleProperties tool: Could not find track entry point!? ("+std::to_string(startvertex.X())+","+std::to_string(startvertex.Y())+","+std::to_string(startvertex.Z())+") --> ("+std::to_string(stopvertex.X())+","+std::to_string(stopvertex.Y())+","+std::to_string(stopvertex.Z())+")",v_error,verbosity); + if(!entryfound) Log("MCParticleProperties tool: Could not find track entry point!? ("+std::to_string(startvertex.X())+","+std::to_string(startvertex.Y())+","+std::to_string(startvertex.Z())+") --> ("+std::to_string(stopvertex.X())+","+std::to_string(stopvertex.Y())+","+std::to_string(stopvertex.Z())+")",v_debug,verbosity); else Hit2.SetZ(startvertex.Z()); if(verbose){ if(entryfound) cout<<"setting entry Z to "<<(Hit2.Z()-tank_start-tank_radius)< ("+std::to_string(stopvertex.X())+","+std::to_string(stopvertex.Y())+","+std::to_string(stopvertex.Z())+")",v_error,verbosity); + if(!exitfound) Log("MCParticleProperties tool: Could not find track exit point!? ("+std::to_string(startvertex.X())+","+std::to_string(startvertex.Y())+","+std::to_string(startvertex.Z())+") --> ("+std::to_string(stopvertex.X())+","+std::to_string(stopvertex.Y())+","+std::to_string(stopvertex.Z())+")",v_debug,verbosity); else Hit.SetZ(startvertex.Z()); if(verbose){ if(exitfound) cout<<"setting exit Z to "<<(Hit.Z()-tank_start-tank_radius)<