diff --git a/DataModel/Hit.h b/DataModel/Hit.h index ba78f611c..b654e7bf9 100755 --- a/DataModel/Hit.h +++ b/DataModel/Hit.h @@ -6,105 +6,161 @@ #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{ + Hit(int thetubeid, double thetime, double thecharge) + : TubeId(thetubeid) + , Time(thetime) + , Charge(thecharge) + { + serialise=true; + } + + virtual ~Hit(){}; - friend class boost::serialization::access; + inline int GetTubeId() const {return TubeId;} + inline double GetTime() const {return Time;} + inline double GetCharge() const {return Charge;} - 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(){}; + inline void SetTubeId(int tubeid){TubeId=tubeid;} + inline void SetTime(double tc){Time=tc;} + inline void SetCharge(double chg){Charge=chg;} - inline int GetTubeId() const {return TubeId;} - inline double GetTime() const {return Time;} - inline double GetCharge() const {return Charge;} + 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 - - friend class boost::serialization::access; + // XXX ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ XXX + // XXX ~~~~~~~~~~~~~~~~~~~~~~~~ UPDATING THIS CLASS? ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ XXX + // XXX ~~~~~ Everything added in this class must be duplicated in MCLAPPDHit!~~~~ XXX + // XXX ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ XXX - 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(){}; - - const std::vector* GetParents() const { return &Parents; } - void SetParents(std::vector parentsin){ Parents = parentsin; } + + friend class boost::serialization::access; - bool Print(){ - std::cout<<"TubeId : "<{}) + , DirectParents(std::vector{}) + , StartTick(-5) + , EndTick(-5) + , IsDarknoise(false) + { + 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) + , IsDarknoise(false) + { + serialise=true; + } + + virtual ~MCHit(){}; - protected: - std::vector Parents; + const std::vector* GetParents() { return &Parents; } + 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; } + - 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 - } + bool Print() + { + std::cout << "TubeId : " << TubeId << std::endl; + std::cout << "Time : " << Time << std::endl; + std::cout << "Charge : " << Charge << std::endl; + + if (Parents.size()) { + std::cout << "Parent MCPartice TrackIDs: {"; + for (uint idx = 0; idx < Parents.size(); ++idx) { + std::cout << Parents.at(idx); + + if ((idx+1) < Parents.size()) std::cout << ", "; + } + std::cout << "}" << std::endl; + } else { + std::cout << "No recorded parents" << std::endl; } -}; -/* -class RecoHit : public Hit { - public: - RecoHit(double thetime, double thecharge) : Time(thetime), Charge(thecharge){}; + std::cout << "IsDarknoise : " << (IsDarknoise ? "true" : "false") << std::endl; - inline double GetCharge(){return Charge;} - inline void SetCharge(double chg){Charge=chg;} + return true; + } - protected: - double Charge; +protected: + std::vector Parents; + std::vector DirectParents; + int StartTick; + int EndTick; + bool IsDarknoise; + + 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; + } + + if (version > 1) { + ar & IsDarknoise; + } + } + } }; -*/ + +BOOST_CLASS_VERSION(MCHit, 2) #endif diff --git a/DataModel/LAPPDHit.h b/DataModel/LAPPDHit.h index 251ac38b8..565c84545 100755 --- a/DataModel/LAPPDHit.h +++ b/DataModel/LAPPDHit.h @@ -100,12 +100,18 @@ 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{}), 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; } + 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; @@ -116,6 +122,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 +138,25 @@ 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; + } + + cout << "IsDarknoise : " << (IsDarknoise ? "true" : "false") << endl; + return true; } @@ -146,13 +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/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); + fANNIETree->Branch("DirectParent_NeutronAncestorClass", &fDirectParent_NeutronAncestorClass); + fANNIETree->Branch("DirectParent_NeutronParentTrackID", &fDirectParent_NeutronParentTrackID); + 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) { fANNIETree->Branch("SiPMhitQ", &fSiPMHitQ); @@ -463,6 +490,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 @@ -589,6 +617,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(); @@ -757,6 +791,31 @@ 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(); + fDirectParent_NeutronAncestorClass.clear(); + fDirectParent_NeutronParentTrackID.clear(); + 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; fSiPM2NPulses = 0; @@ -999,6 +1058,7 @@ void ANNIEEventTreeMaker::ResetVariables() fTrueKPlusCher = -9999; fTrueKMinus = -9999; fTrueKMinusCher = -9999; + fTrueWCSimMode = -9999; // TankReco_fill fRecoVtxX = -9999; @@ -1377,6 +1437,226 @@ 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; + 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) { + 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; + } + 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); + 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; + 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; + 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()){ + auto const& ancestorPair = pmtAncestors.at(hitTime); + neutronAncestorTrackID = ancestorPair.first; + 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); + + 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); + + 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); + + 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; +} + void ANNIEEventTreeMaker::LoadSiPMHits() { Log("ANNIEEventTreeMaker Tool: LoadSiPMHits", v_debug, ANNIEEventTreeMakerVerbosity); @@ -2328,6 +2608,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 44a321dba..887b89a59 100644 --- a/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h +++ b/UserTools/ANNIEEventTreeMaker/ANNIEEventTreeMaker.h @@ -50,6 +50,7 @@ class ANNIEEventTreeMaker : public Tool void LoadRWMBRFInfo(); void LoadAllTankHits(); + void LoadDirectParentIDsMCHits(); void LoadSiPMHits(); void LoadLAPPDInfo(); @@ -212,6 +213,53 @@ 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. + 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 + 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; int fSiPM2NPulses; @@ -461,6 +509,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 ce80ada9b..f164caa5c 100644 --- a/UserTools/BackTracker/BackTracker.cpp +++ b/UserTools/BackTracker/BackTracker.cpp @@ -1,4 +1,6 @@ #include "BackTracker.h" +#include "ANNIEconstants.h" +#include BackTracker::BackTracker():Tool(){} @@ -27,6 +29,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 +40,23 @@ 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>>; + fMCHitToNeutronAncestorClass = new std::map>; + 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; } -//------------------------------------------------------------------------------ bool BackTracker::Execute() { if (!LoadFromStores()) @@ -50,11 +67,30 @@ bool BackTracker::Execute() fClusterEfficiency ->clear(); fClusterPurity ->clear(); fClusterTotalCharge ->clear(); + fMCHitToDirectParents ->clear(); + fMCHitToNeutronAncestor ->clear(); + fMCHitToNeutronAncestorClass->clear(); + fMCHitToNeutronParent->clear(); + fMCHitToIsDarknoise->clear(); + fMCHitToInteractionMode->clear(); + fMCHitToImmediateAncestor->clear(); + fMCHitToImmediateAncestorClass->clear(); + // DISABLED: PrimaryAncestor + // fMCHitToPrimaryAncestor->clear(); + fMCHitToRootAncestor->clear(); + fMCHitToLineage->clear(); + fMCHitToLineageStatus->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 +106,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 +114,19 @@ 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 ); + 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); + 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; } @@ -174,7 +224,272 @@ 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 -> (DirectParentID, pdg) + for (auto& particle : *fMCParticles) { + int trackId = particle.GetParticleID(); + int directParentId = particle.GetDirectParentID(); + int pdg = particle.GetPdgCode(); + + trackMap[trackId] = std::make_pair(directParentId, 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; + + // 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; }); + // 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; + 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); + (*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; + } + + int startParentID = directParents[0]; // take the first direct parent of neutron 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) + 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) - If the classification portion works that this question is answered. + } + + if (parentID == -1) break; // reached the end of the ancestry + currentID = parentID; + } + + (*fMCHitToNeutronAncestor)[pmtID][hitTime] = std::make_pair(neutronAncestorId, neutronAncestorPdg); + (*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]); + 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; + + } + } + + if (verbosity > 1) { + std::cout << "BackTracker::FindNeutronAncestors: found " + << fMCHitToNeutronAncestor->size() + << " 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 @@ -208,6 +523,38 @@ 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; + } + } + + 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 1cf5c7240..340cc0aec 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" @@ -14,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 { @@ -31,14 +33,22 @@ 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: + // 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 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 +67,99 @@ 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; + // 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; + // 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 + + // 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; + + + /// \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/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 da9d07573..cb69c7e5e 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 @@ -490,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) { @@ -499,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 ); @@ -634,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 @@ -1101,12 +1111,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 +1130,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 +1160,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 +1168,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 +1186,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 +1214,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 +1404,13 @@ 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::vector primaryParents; + std::vector directParents; + bool hitIsDarknoise = false; + std::tie(primaryParents, directParents, hitIsDarknoise) = GetHitParentIDs(digiHit, firstTrig); + + MCHit nextHit(key, digiTime, digiQ, primaryParents, directParents); + nextHit.SetIsDarknoise(hitIsDarknoise); if (system == "Tank") { if (MCHits->count(key) == 0) MCHits->emplace(key, std::vector{nextHit}); @@ -1455,7 +1484,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,47 +1515,67 @@ 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::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; + // 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->GetParentID()); - - }// end loop over photons - return parentIDs; + else { + 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_tuple(parentIDs, directParentIDs, isDarknoise); } //////////////////////////////////////////////////////////////////////////////// // 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); + 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; // 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..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,9 +122,10 @@ 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::vector GetHitParentIDs(WCSimRootCherenkovDigiHit* digiHit, WCSimRootTrigger* firstTrig); - std::vector GetHitParentIdxs(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; 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"< 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)< #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: 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) { @@ -199,6 +208,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); @@ -215,13 +228,39 @@ bool PMTWaveformSim::Execute() 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 @@ -232,23 +271,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() { @@ -457,6 +504,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) @@ -513,7 +561,6 @@ int PMTWaveformSim::LoadFromStores() return 2; } - return 1; } @@ -521,7 +568,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()); @@ -569,5 +616,3 @@ bool PMTWaveformSim::TimeSmearing(int pmtid) } - - diff --git a/UserTools/PMTWaveformSim/PMTWaveformSim.h b/UserTools/PMTWaveformSim/PMTWaveformSim.h index 9b7b9e1ba..bd529f07c 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 bool 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