From 31674fcfc03345d987d9d94769e68d628017939a Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 17 Dec 2025 14:58:53 -0500 Subject: [PATCH 1/8] custom trajectory class implementation --- README.md | 2 +- include/AnalysisManager.hh | 15 ++- include/FPFTrajectory.hh | 92 ++++++++++++++++++ src/AnalysisManager.cc | 46 +++++---- src/FPFTrajectory.cc | 188 +++++++++++++++++++++++++++++++++++++ src/TrackingAction.cc | 29 +++--- 6 files changed, 327 insertions(+), 45 deletions(-) create mode 100644 include/FPFTrajectory.hh create mode 100644 src/FPFTrajectory.cc diff --git a/README.md b/README.md index 1321443..875906f 100644 --- a/README.md +++ b/README.md @@ -186,7 +186,7 @@ Older versions of FORESEE output events in the HepMC2 format. To run over HepMC2 |Command |Description | |:--|:--| |/out/fileName | option for AnalysisManagerMessenger, set name of the file saving all analysis variables| -|/out/saveTrack | if `true` save all tracks, `false` by default, requires `\tracking\storeTrajectory 1`| +|/out/saveTrack | if `true` save full trajectories, `false` by default | |/out/flare/enableOutput | if `true` save FLArE output, `true` by default | |/out/flare/save3DEvd | if `true` save 3D spatial distribution of energy deposition, `false` by default | |/out/flare/save2DEvd | if `true` save 2D spatial distribution of energy deposition, `false` by default | diff --git a/include/AnalysisManager.hh b/include/AnalysisManager.hh index 672a3a1..86d5610 100644 --- a/include/AnalysisManager.hh +++ b/include/AnalysisManager.hh @@ -47,8 +47,8 @@ class AnalysisManager { void SetTrackPrimaryAncestor(G4int trackID, G4int ancestorID) { trackToPrimaryAncestor[trackID] = ancestorID; } G4int GetTrackPrimaryAncestor(G4int trackID) { return trackToPrimaryAncestor.at(trackID); } - // TODO: needed??? - void AddOnePrimaryTrack() { nTestNPrimaryTrack++; } + // return whether saving full tracks in trajectories + G4bool GetSaveTrack() { return fSaveTrack; } private: @@ -56,14 +56,14 @@ class AnalysisManager { // Book ROOT output TTrees // common + detector specific void bookEvtTree(); - void bookTrkTree(); + void bookParTree(); void bookPrimTree(); void bookFLArETrees(); void bookFASER2Trees(); void FillEventTree(const G4Event* event); void FillPrimariesTree(const G4Event* event); - void FillTrajectoriesTree(const G4Event* event); + void FillParticlesTree(const G4Event* event); void FillFLArEOutput(); void FillFLArEPseudoReco(); @@ -100,7 +100,7 @@ class AnalysisManager { hep_hpc::hdf5::File fH5file; TFile* fFile; TTree* fEvt; - TTree* fTrk; + TTree* fPar; TTree* fPrim; TDirectory* fFLArEDir; @@ -115,9 +115,6 @@ class AnalysisManager { // track to primary ancestor std::map trackToPrimaryAncestor; - // TODO: no longer needed? - G4int nTestNPrimaryTrack; - //--------------------------------------------------- // OUTPUT VARIABLES FOR COMMON TREES @@ -145,7 +142,7 @@ class AnalysisManager { double W; //--------------------------------------------------- - // Output variables for TRAJECTORIES tree + // Output variables for PARTICLES tree int trackTID; int trackPID; int trackPDG; diff --git a/include/FPFTrajectory.hh b/include/FPFTrajectory.hh new file mode 100644 index 0000000..ad0288b --- /dev/null +++ b/include/FPFTrajectory.hh @@ -0,0 +1,92 @@ +#ifndef FPFTRAJECTORY_HH +#define FPFTRAJECTORY_HH + +#include "G4VTrajectory.hh" +#include "G4VTrajectoryPoint.hh" +#include "G4Allocator.hh" +#include "G4ios.hh" +#include "G4ParticleDefinition.hh" +#include "G4TrajectoryPoint.hh" +#include "G4Track.hh" +#include "G4Step.hh" + +/// FPFTrajectory is a custom trajectory class that allows to disable +/// storing the full trajectory point collections (memory intensive) +/// while still bookkeeping the full hierarchial particle tree info. + +typedef std::vector TrajectoryPointContainer; +class FPFTrajectory : public G4VTrajectory +{ +public: + + FPFTrajectory(); + // if storePoints is enabled, full trajectories are saved + // if not (default), single steps are not saved + FPFTrajectory(const G4Track *aTrack, G4bool storePoint = false); + FPFTrajectory(FPFTrajectory &); + virtual ~FPFTrajectory(); + + inline void *operator new(size_t); + inline void operator delete(void *); + inline int operator==(const FPFTrajectory &right) const + { + return (this == &right); + } + + inline G4int GetTrackID() const { return fTrackID; } + inline G4int GetParentID() const { return fParentID; } + inline G4String GetParticleName() const { return fParticleName; } + inline G4double GetCharge() const { return fPDGCharge; } + inline G4int GetPDGEncoding() const { return fPDGEncoding; } + inline G4ThreeVector GetInitialMomentum() const { return fInitialMomentum; } + G4double GetInitialKineticEnergy() const; + inline G4String GetProcessName() const { return fProcessName; } + + virtual void ShowTrajectory(std::ostream& os=G4cout) const; + virtual void DrawTrajectory(G4int i_mode = 0) const; + + virtual void AppendStep(const G4Step *aStep); + + virtual int GetPointEntries() const { return fPositionRecord->size(); } + virtual G4VTrajectoryPoint *GetPoint(G4int i) const { return (*fPositionRecord)[i]; } + + virtual void MergeTrajectory(G4VTrajectory *secondTrajectory); + + G4ParticleDefinition *GetParticleDefinition() const; + + virtual const std::map *GetAttDefs() const; + virtual std::vector *CreateAttValues() const; + +private: + + // from standard G4Trajectory implementation + TrajectoryPointContainer *fPositionRecord; + G4int fTrackID; + G4int fParentID; + G4int fPDGEncoding; + G4double fPDGCharge; + G4String fParticleName; + G4ThreeVector fInitialMomentum; + + // custom additions + G4String fProcessName; // creator process + G4bool fStorePoints; // whether to save trajectory points + +}; + +#if defined G4TRACKING_ALLOC_EXPORT +extern G4DLLEXPORT G4Allocator aTrajAllocator; +#else +extern G4DLLIMPORT G4Allocator aTrajAllocator; +#endif + +inline void* FPFTrajectory::operator new(size_t) { + void* aTrajectory = (void*) aTrajAllocator.MallocSingle(); + return aTrajectory; +} + +inline void FPFTrajectory::operator delete(void* aTrajectory) { + aTrajAllocator.FreeSingle((FPFTrajectory*)aTrajectory); +} + +#endif \ No newline at end of file diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index ad9c7aa..2c8941f 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -11,7 +11,6 @@ #include #include #include -#include #include #include @@ -35,6 +34,8 @@ #include "reco/ShowerLID.hh" #include "reco/Barcode.hh" #include "FPFParticle.hh" +#include "FPFTrajectory.hh" + //--------------------------------------------------------------------- //--------------------------------------------------------------------- @@ -64,7 +65,7 @@ AnalysisManager::AnalysisManager() fMessenger = new AnalysisManagerMessenger(this); fEvt = nullptr; - fTrk = nullptr; + fPar = nullptr; fPrim = nullptr; fFLArEHits = nullptr; fFLArEHCALHits = nullptr; @@ -144,18 +145,23 @@ void AnalysisManager::bookPrimTree() fPrim->Branch("P", &primP, "P/F"); } -void AnalysisManager::bookTrkTree() +void AnalysisManager::bookParTree() { - fTrk = new TTree("trajectories", "trajectories info"); - fTrk->Branch("evtID", &evtID, "evtID/I"); - fTrk->Branch("trackTID", &trackTID, "trackTID/I"); - fTrk->Branch("trackPID", &trackPID, "trackPID/I"); - fTrk->Branch("trackPDG", &trackPDG, "trackPDG/I"); - fTrk->Branch("trackKinE", &trackKinE, "trackKinE/D"); - fTrk->Branch("trackNPoints", &trackNPoints, "trackNPoints/I"); - fTrk->Branch("trackPointX", &trackPointX); - fTrk->Branch("trackPointY", &trackPointY); - fTrk->Branch("trackPointZ", &trackPointZ); + fPar = new TTree("particles", "particle info"); + fPar->Branch("evtID", &evtID, "evtID/I"); + fPar->Branch("trackTID", &trackTID, "trackTID/I"); + fPar->Branch("trackPID", &trackPID, "trackPID/I"); + fPar->Branch("trackPDG", &trackPDG, "trackPDG/I"); + fPar->Branch("trackKinE", &trackKinE, "trackKinE/D"); + fPar->Branch("trackNPoints", &trackNPoints, "trackNPoints/I"); + + // if saving full trajectory, add vector of trajectory points + if (fSaveTrack) + { + fPar->Branch("trackPointX", &trackPointX); + fPar->Branch("trackPointY", &trackPointY); + fPar->Branch("trackPointZ", &trackPointZ); + } } //--------------------------------------------------------------------- @@ -327,7 +333,7 @@ void AnalysisManager::BeginOfRun() // Booking common output trees bookEvtTree(); bookPrimTree(); - if (fSaveTrack) bookTrkTree(); + bookParTree(); fSDNamelist = GeometricalParameters::Get()->GetSDNamelist(); G4cout << "Number of SDs : " << fSDNamelist.size() << G4endl; @@ -379,7 +385,7 @@ void AnalysisManager::EndOfRun() fFile->cd(); fEvt->Write(); fPrim->Write(); - if (fSaveTrack) fTrk->Write(); + fPar->Write(); // save detector-specif trees in their directories if (fFlareSDs.size()>0) @@ -495,7 +501,7 @@ void AnalysisManager::EndOfEvent(const G4Event *event) // FILL PRIMARIES/TRAJECTORIES TREE FillPrimariesTree(event); - if(fSaveTrack) FillTrajectoriesTree(event); + FillParticlesTree(event); //----------------------------------------------------------- @@ -635,7 +641,7 @@ void AnalysisManager::FillPrimariesTree(const G4Event *event) //--------------------------------------------------------------------- //--------------------------------------------------------------------- -void AnalysisManager::FillTrajectoriesTree(const G4Event* event) +void AnalysisManager::FillParticlesTree(const G4Event* event) { int count_tracks = 0; @@ -643,13 +649,13 @@ void AnalysisManager::FillTrajectoriesTree(const G4Event* event) auto trajectoryContainer = event->GetTrajectoryContainer(); if (!trajectoryContainer) { - G4cout << "No tracks found: did you enable their storage with '/tracking/storeTrajectory 1'?" << G4endl; + G4cout << "No tracks found in the event!" << G4endl; return; } for (size_t i = 0; i < trajectoryContainer->entries(); ++i) { - auto trajectory = static_cast((*trajectoryContainer)[i]); + auto trajectory = static_cast((*trajectoryContainer)[i]); trackTID = trajectory->GetTrackID(); trackPID = trajectory->GetParentID(); trackPDG = trajectory->GetPDGEncoding(); @@ -663,7 +669,7 @@ void AnalysisManager::FillTrajectoriesTree(const G4Event* event) trackPointY.push_back( pos.y() ); trackPointZ.push_back( pos.z() ); } - fTrk->Fill(); + fPar->Fill(); trackPointX.clear(); trackPointY.clear(); trackPointZ.clear(); diff --git a/src/FPFTrajectory.cc b/src/FPFTrajectory.cc new file mode 100644 index 0000000..6e5c268 --- /dev/null +++ b/src/FPFTrajectory.cc @@ -0,0 +1,188 @@ +#include "G4VProcess.hh" +#include "G4ParticleTable.hh" +#include "G4AttDefStore.hh" +#include "G4AttDef.hh" +#include "G4AttValue.hh" +#include "G4UnitsTable.hh" + +#include "FPFTrajectory.hh" + +G4Allocator aTrajAllocator; + +FPFTrajectory::FPFTrajectory() + : fPositionRecord(0), fTrackID(0), fParentID(0), + fPDGEncoding(0), fPDGCharge(0.0), fParticleName(""), + fProcessName(""), fInitialMomentum(G4ThreeVector()) +{} + +FPFTrajectory::FPFTrajectory(const G4Track *aTrack, G4bool storePoints) +{ + fTrackID = aTrack->GetTrackID(); + fParentID = aTrack->GetParentID(); + fInitialMomentum = aTrack->GetMomentum(); + + G4ParticleDefinition *fpParticleDefinition = aTrack->GetDefinition(); + fPDGEncoding = fpParticleDefinition->GetPDGEncoding(); + fPDGCharge = fpParticleDefinition->GetPDGCharge(); + fParticleName = fpParticleDefinition->GetParticleName(); + + // set creator process + const G4VProcess *proc = aTrack->GetCreatorProcess(); + if (proc) + fProcessName = proc->GetProcessName(); + else + fProcessName = "primary"; + + // set full trajectory storing + fStorePoints = storePoints; + + // anyways always initialize the container + fPositionRecord = new TrajectoryPointContainer(); + // and add the starting point (vertex) + fPositionRecord->push_back(new G4TrajectoryPoint(aTrack->GetPosition())); +} + +FPFTrajectory::FPFTrajectory(FPFTrajectory &right) : G4VTrajectory() +{ + fParticleName = right.fParticleName; + fPDGCharge = right.fPDGCharge; + fPDGEncoding = right.fPDGEncoding; + fTrackID = right.fTrackID; + fParentID = right.fParentID; + fInitialMomentum = right.fInitialMomentum; + fPositionRecord = new TrajectoryPointContainer(); + fStorePoints = right.fStorePoints; + fProcessName = right.fProcessName; + + for (size_t i = 0; i < right.fPositionRecord->size(); i++) + { + G4TrajectoryPoint *rightPoint = (G4TrajectoryPoint *)((*(right.fPositionRecord))[i]); + fPositionRecord->push_back(new G4TrajectoryPoint(*rightPoint)); + } +} + +FPFTrajectory::~FPFTrajectory() +{ + for (size_t i = 0; i < fPositionRecord->size(); i++) + { + delete (*fPositionRecord)[i]; + } + fPositionRecord->clear(); + delete fPositionRecord; +} + +void FPFTrajectory::ShowTrajectory(std::ostream &os) const +{ + G4VTrajectory::ShowTrajectory(os); +} + +void FPFTrajectory::DrawTrajectory(G4int i_mode) const +{ + //G4VTrajectory::DrawTrajectory(i_mode); + G4VTrajectory::DrawTrajectory(); +} + +const std::map *FPFTrajectory::GetAttDefs() const +{ + G4bool isNew; + std::map *store = G4AttDefStore::GetInstance("FPFTrajectory", isNew); + if (isNew) + { + G4String ID("ID"); + (*store)[ID] = G4AttDef(ID, "Track ID", "Bookkeeping", "", "G4int"); + + G4String PID("PID"); + (*store)[PID] = G4AttDef(PID, "Parent ID", "Bookkeeping", "", "G4int"); + + G4String PN("PN"); + (*store)[PN] = G4AttDef(PN, "Particle Name", "Physics", "", "G4String"); + + G4String Ch("Ch"); + (*store)[Ch] = G4AttDef(Ch, "Charge", "Physics", "", "G4double"); + + G4String PDG("PDG"); + (*store)[PDG] = G4AttDef(PDG, "PDG Encoding", "Physics", "", "G4int"); + + G4String IMom("IMom"); + (*store)[IMom] = G4AttDef(IMom, "Momentum of track at start of trajectory", "Physics", "", "G4ThreeVector"); + + G4String NTP("NTP"); + (*store)[NTP] = G4AttDef(NTP, "No. of points", "Physics", "", "G4int"); + } + return store; +} + +std::vector *FPFTrajectory::CreateAttValues() const +{ + std::string c; + std::ostringstream s(c); + std::vector *values = new std::vector; + + s.seekp(std::ios::beg); + s << fTrackID << std::ends; + values->push_back(G4AttValue("ID", c.c_str(), "")); + + s.seekp(std::ios::beg); + s << fParentID << std::ends; + values->push_back(G4AttValue("PID", c.c_str(), "")); + + values->push_back(G4AttValue("PN", fParticleName, "")); + + s.seekp(std::ios::beg); + s << fPDGCharge << std::ends; + values->push_back(G4AttValue("Ch", c.c_str(), "")); + + s.seekp(std::ios::beg); + s << fPDGEncoding << std::ends; + values->push_back(G4AttValue("PDG", c.c_str(), "")); + + s.seekp(std::ios::beg); + s << G4BestUnit(fInitialMomentum, "Energy") << std::ends; + values->push_back(G4AttValue("IMom", c.c_str(), "")); + + s.seekp(std::ios::beg); + s << GetPointEntries() << std::ends; + values->push_back(G4AttValue("NTP", c.c_str(), "")); + + return values; +} + +void FPFTrajectory::AppendStep(const G4Step *aStep) +{ + // If we are not interested in trajectory points, do nothing + if (!fStorePoints) + return; + + // otherwise, keep the standard behavior (store positions) + fPositionRecord->push_back(new G4TrajectoryPoint(aStep->GetPostStepPoint()->GetPosition())); +} + +G4ParticleDefinition *FPFTrajectory::GetParticleDefinition() const +{ + return (G4ParticleTable::GetParticleTable()->FindParticle(fParticleName)); +} + +void FPFTrajectory::MergeTrajectory(G4VTrajectory *secondTrajectory) +{ + if (!secondTrajectory) + return; + FPFTrajectory *second = (FPFTrajectory *)secondTrajectory; + G4int ent = second->GetPointEntries(); + // initial point of the second trajectory should not be merged + for (G4int i = 1; i < ent; ++i) + { + fPositionRecord->push_back((*(second->fPositionRecord))[i]); + } + delete (*second->fPositionRecord)[0]; + second->fPositionRecord->clear(); +} + +G4double FPFTrajectory::GetInitialKineticEnergy() const { + const G4ParticleDefinition* p = GetParticleDefinition(); + double mom = GetInitialMomentum().mag(); + if (!p) return mom; + double mass = p->GetPDGMass(); + double kin = std::sqrt(mom*mom + mass*mass) - mass; + if (kin<0.0) return 0.0; + return kin; +} \ No newline at end of file diff --git a/src/TrackingAction.cc b/src/TrackingAction.cc index 900d7fa..46e6db6 100644 --- a/src/TrackingAction.cc +++ b/src/TrackingAction.cc @@ -1,6 +1,7 @@ #include "TrackingAction.hh" #include "TrackInformation.hh" #include "AnalysisManager.hh" +#include "FPFTrajectory.hh" #include "G4TrackingManager.hh" #include "G4Track.hh" @@ -9,14 +10,23 @@ TrackingAction::TrackingAction() : G4UserTrackingAction() {;} void TrackingAction::PreUserTrackingAction(const G4Track* aTrack) { + // find out whether saving the full trajectory points or not + bool storeTrajectoryPoints = AnalysisManager::GetInstance()->GetSaveTrack(); + + // initialize our custom trajectory class + auto *traj = new FPFTrajectory(aTrack, storeTrajectoryPoints); + + // hand it over to the manager + fpTrackingManager->SetTrajectory(traj); + fpTrackingManager->SetStoreTrajectory(true); } void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) { - if (aTrack->GetParentID()==0) - { - AnalysisManager::GetInstance()->AddOnePrimaryTrack(); - } + /// the following blocks store additional information into the track + /// via the TrackInformation object (G4VTrackUserInformation) + /// it is currently used in some FLARE SD volume to change how to save things + /// TODO: consider for removal to simplify structure! if (aTrack->GetParentID()==0) { if (aTrack->GetParticleDefinition()->GetPDGEncoding()==111) @@ -36,7 +46,6 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) TrackInformation* info = new TrackInformation(); info->SetTrackIsFromPrimaryPizero(1); (*secondaries)[i]->SetUserInformation(info); - AnalysisManager::GetInstance()->AddOnePrimaryTrack(); } } } @@ -64,7 +73,6 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) TrackInformation* info = new TrackInformation(); info->SetTrackIsFromPrimaryLepton(1); (*secondaries)[i]->SetUserInformation(info); - AnalysisManager::GetInstance()->AddOnePrimaryTrack(); } } } @@ -90,19 +98,10 @@ void TrackingAction::PostUserTrackingAction(const G4Track* aTrack) TrackInformation* info = new TrackInformation(); info->SetTrackIsFromFSLPizero(1); (*secondaries)[i]->SetUserInformation(info); - AnalysisManager::GetInstance()->AddOnePrimaryTrack(); } } } } } } - - //TrackInformation* aTrackInfo = (TrackInformation*)(aTrack->GetUserInformation()); - //if (aTrackInfo) { - // if (aTrackInfo->IsTrackFromPrimaryTau() | aTrackInfo->IsTrackFromPrimaryPizero()) { - // std::cout<GetParentID()<<" "<GetParticleDefinition()->GetPDGEncoding()<Print(); - // } - //} } From 3d14e9981dfb45941954cf6f64be1afbda389a18 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Wed, 17 Dec 2025 16:25:20 -0500 Subject: [PATCH 2/8] remove primaries tree --- include/AnalysisManager.hh | 51 ++--- src/AnalysisManager.cc | 452 ++++++++++++++++++------------------- 2 files changed, 244 insertions(+), 259 deletions(-) diff --git a/include/AnalysisManager.hh b/include/AnalysisManager.hh index 86d5610..2841341 100644 --- a/include/AnalysisManager.hh +++ b/include/AnalysisManager.hh @@ -57,12 +57,10 @@ class AnalysisManager { // common + detector specific void bookEvtTree(); void bookParTree(); - void bookPrimTree(); void bookFLArETrees(); void bookFASER2Trees(); void FillEventTree(const G4Event* event); - void FillPrimariesTree(const G4Event* event); void FillParticlesTree(const G4Event* event); void FillFLArEOutput(); @@ -117,30 +115,45 @@ class AnalysisManager { //--------------------------------------------------- // OUTPUT VARIABLES FOR COMMON TREES + //--------------------------------------------------- + // Output variables for EVENT tree + // general vertex metadata G4int evtID; G4int vertexID; double weight; std::string genType; - std::string processName; + std::string processName; + // vertex position + double vtxX, vtxY, vtxZ, vtxT; + // initiator info (incoming particle) int initPDG; - double initX, initY, initZ, initT; double initPx, initPy, initPz, initE; double initM; - double initQ; - int intType; - int scatteringType; + double initQ; + // target info int fslPDG; int tgtPDG; int tgtA; int tgtZ; int hitnucPDG; + // interaction info + int intType; + int scatteringType; double xs; double Q2; double xBj; double y; double W; - + // primaries from vertex + int nPrimaries; + std::vector primTID; + std::vector primPDG; + std::vector primPx; + std::vector primPy; + std::vector primPz; + std::vector primE; + //--------------------------------------------------- // Output variables for PARTICLES tree int trackTID; @@ -152,28 +165,6 @@ class AnalysisManager { std::vector trackPointY; std::vector trackPointZ; - //--------------------------------------------------- - // Output variables for PRIMARIES tree - UInt_t primVtxID; - UInt_t primParticleID; - UInt_t primTrackID; - UInt_t primPDG; // why unsigned? - float_t primM; - float_t primQ; - float_t primEta; - float_t primPhi; - float_t primPt; - float_t primP; - float_t primVx; - float_t primVy; - float_t primVz; - float_t primVt; - float_t primPx; - float_t primPy; - float_t primPz; - float_t primE; - float_t primKE; - //--------------------------------------------------- // OUTPUT VARIABLES FOR FLArE TREES // TODO: merge hit variables? no need to use different names? diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index 2c8941f..2b1c2c4 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -2,6 +2,7 @@ #include #include #include +#include #include #include #include @@ -12,6 +13,7 @@ #include #include #include +#include #include #include @@ -36,7 +38,6 @@ #include "FPFParticle.hh" #include "FPFTrajectory.hh" - //--------------------------------------------------------------------- //--------------------------------------------------------------------- // AnalysisManager "singleton" instance @@ -66,13 +67,12 @@ AnalysisManager::AnalysisManager() fEvt = nullptr; fPar = nullptr; - fPrim = nullptr; fFLArEHits = nullptr; fFLArEHCALHits = nullptr; fFLArEPseudoReco = nullptr; fActsHitsTree = nullptr; fActsParticlesTree = nullptr; - + fSaveTrack = false; fSave3DEvd = false; fSave2DEvd = false; @@ -90,59 +90,47 @@ AnalysisManager::~AnalysisManager() {} void AnalysisManager::bookEvtTree() { fEvt = new TTree("event", "event info"); + // vertex generation metadata fEvt->Branch("evtID", &evtID, "evtID/I"); fEvt->Branch("vtxID", &vertexID, "vtxID/I"); fEvt->Branch("weight", &weight, "weight/D"); - fEvt->Branch("genType", &genType); + fEvt->Branch("generatorType", &genType); fEvt->Branch("processName", &processName); + // vertex position + fEvt->Branch("vtxX", &vtxX, "vtxX/D"); + fEvt->Branch("vtxY", &vtxY, "vtxY/D"); + fEvt->Branch("vtxZ", &vtxZ, "vtxZ/D"); + fEvt->Branch("vtxT", &vtxT, "vtxT/D"); + // initiator info (incoming particle) fEvt->Branch("initPDG", &initPDG, "initPDG/I"); - fEvt->Branch("initX", &initX, "initX/D"); - fEvt->Branch("initY", &initY, "initY/D"); - fEvt->Branch("initZ", &initZ, "initZ/D"); - fEvt->Branch("initT", &initT, "initT/D"); fEvt->Branch("initPx", &initPx, "initPx/D"); fEvt->Branch("initPy", &initPy, "initPy/D"); - fEvt->Branch("initPz", &initPz, "initPz/D"); + fEvt->Branch("initPz", &initPz, "initPz/D"); fEvt->Branch("initE", &initE, "initE/D"); fEvt->Branch("initM", &initM, "initM/D"); fEvt->Branch("initQ", &initQ, "initQ/D"); - fEvt->Branch("intType", &intType, "intType/I"); - fEvt->Branch("scatteringType", &scatteringType, "scatteringType/I"); - fEvt->Branch("fslPDG", &fslPDG, "fslPDG/I"); + // target info fEvt->Branch("tgtPDG", &tgtPDG, "tgtPDG/I"); fEvt->Branch("tgtA", &tgtA, "tgtA/I"); fEvt->Branch("tgtZ", &tgtZ, "tgtZ/I"); fEvt->Branch("hitnucPDG", &hitnucPDG, "hitnucPDG/I"); + fEvt->Branch("fslPDG", &fslPDG, "fslPDG/I"); + // interaction info + fEvt->Branch("intType", &intType, "intType/I"); + fEvt->Branch("scatteringType", &scatteringType, "scatteringType/I"); fEvt->Branch("xs", &xs, "xs/D"); fEvt->Branch("Q2", &Q2, "Q2/D"); fEvt->Branch("xBj", &xBj, "xBj/D"); fEvt->Branch("y", &y, "y/D"); fEvt->Branch("W", &W, "W/D"); -} - -void AnalysisManager::bookPrimTree() -{ - fPrim = new TTree("primaries", "primaries info"); - fPrim->Branch("evtID", &evtID, "evtID/I"); - fPrim->Branch("vtxID", &primVtxID, "vtxID/I"); - fPrim->Branch("PDG", &primPDG, "PDG/I"); - fPrim->Branch("trackID", &primTrackID, "trackID/I"); - fPrim->Branch("barcode", &primParticleID, "bardcode/I"); - fPrim->Branch("mass", &primM, "mass/F"); - fPrim->Branch("charge", &primQ, "charge/F"); - fPrim->Branch("Vx", &primVx, "Vx/F"); // position - fPrim->Branch("Vy", &primVy, "Vy/F"); - fPrim->Branch("Vz", &primVz, "Vz/F"); - fPrim->Branch("Vt", &primVt, "Vt/F"); - fPrim->Branch("Px", &primPx, "Px/F"); // momentum - fPrim->Branch("Py", &primPy, "Py/F"); - fPrim->Branch("Pz", &primPz, "Pz/F"); - fPrim->Branch("E", &primE, "E/F"); // initial total energy - fPrim->Branch("KE", &primKE, "KE/F"); // initial kinetic energy - fPrim->Branch("Eta", &primEta, "Eta/F"); - fPrim->Branch("Phi", &primPhi, "Phi/F"); - fPrim->Branch("Pt", &primPt, "Pt/F"); - fPrim->Branch("P", &primP, "P/F"); + // primaries from vertex + fEvt->Branch("nPrimaries", &nPrimaries, "nPrimaries/I"); + fEvt->Branch("primTID", &primTID); + fEvt->Branch("primPDG", &primPDG); + fEvt->Branch("primPx", &primPx); + fEvt->Branch("primPy", &primPy); + fEvt->Branch("primPz", &primPz); + fEvt->Branch("primE", &primE); } void AnalysisManager::bookParTree() @@ -170,7 +158,7 @@ void AnalysisManager::bookParTree() void AnalysisManager::bookFLArETrees() { // create subdirectory in file - fFLArEDir = fFile->mkdir("flare","FLArE output",kTRUE); + fFLArEDir = fFile->mkdir("flare", "FLArE output", kTRUE); fFile->cd(fFLArEDir->GetName()); // FLArE hits tree @@ -257,7 +245,7 @@ void AnalysisManager::bookFLArETrees() void AnalysisManager::bookFASER2Trees() { // create subdirectory in file - fFASER2Dir = fFile->mkdir("faser2","FASER2 output",kTRUE); + fFASER2Dir = fFile->mkdir("faser2", "FASER2 output", kTRUE); fFile->cd(fFASER2Dir->GetName()); //* Acts Hits Tree [i == unsigned int; F == float; l == Long unsigned 64 int] @@ -329,10 +317,9 @@ void AnalysisManager::BeginOfRun() // Preparing output file fFile = new TFile(fFilename.c_str(), "RECREATE"); - + // Booking common output trees bookEvtTree(); - bookPrimTree(); bookParTree(); fSDNamelist = GeometricalParameters::Get()->GetSDNamelist(); @@ -345,7 +332,7 @@ void AnalysisManager::BeginOfRun() G4cout << sdname.first << " " << sdname.second << G4endl; // if FLArE is enabled, book trees - if ( fEnableFLArE && sdname.second.find("FLArE") != std::string::npos ) + if (fEnableFLArE && sdname.second.find("FLArE") != std::string::npos) { fFlareSDs.push_back(sdname.first); bookFLArETrees(); @@ -353,7 +340,7 @@ void AnalysisManager::BeginOfRun() // if FASER2 is enabled, book ACTS trees // but give the option to disable output if required - else if (fSaveActs && sdname.second.find("FASER2") != std::string::npos ) + else if (fSaveActs && sdname.second.find("FASER2") != std::string::npos) { fFaser2SDs.push_back(sdname.first); bookFASER2Trees(); @@ -362,7 +349,7 @@ void AnalysisManager::BeginOfRun() // if FLArE is enabled & its geometry was found // prepare .h5 output file - if( fFlareSDs.size()>0 ) + if (fFlareSDs.size() > 0) { fH5Filename = fFilename; if (fH5Filename.find(".root") != std::string::npos) @@ -373,7 +360,6 @@ void AnalysisManager::BeginOfRun() fH5Filename += ".h5"; fH5file = hep_hpc::hdf5::File(fH5Filename, H5F_ACC_TRUNC); } - } //--------------------------------------------------------------------- @@ -384,20 +370,20 @@ void AnalysisManager::EndOfRun() // save common trees at the top of the output file fFile->cd(); fEvt->Write(); - fPrim->Write(); fPar->Write(); // save detector-specif trees in their directories - if (fFlareSDs.size()>0) + if (fFlareSDs.size() > 0) { fFile->cd(fFLArEDir->GetName()); fFLArEHits->Write(); fFLArEHCALHits->Write(); - if(fSavePseudoReco) fFLArEPseudoReco->Write(); + if (fSavePseudoReco) + fFLArEPseudoReco->Write(); fFile->cd(); // go back to top fH5file.close(); } - if (fFaser2SDs.size()>0) + if (fFaser2SDs.size() > 0) { fFile->cd(fFASER2Dir->GetName()); fActsHitsTree->Write(); @@ -419,9 +405,20 @@ void AnalysisManager::BeginOfEvent() primaries.clear(); primaryIDs.clear(); + // primaries in event tree + // (actually need reset at every vertex) + primTID.clear(); + primPDG.clear(); + primPx.clear(); + primPy.clear(); + primPz.clear(); + primE.clear(); + // track ID to primary ancestor association trackToPrimaryAncestor.clear(); + // trajectory points (if enabled) + // (actually need reset at every track) trackPointX.clear(); trackPointY.clear(); trackPointZ.clear(); @@ -461,12 +458,13 @@ void AnalysisManager::BeginOfEvent() TotalDedxLongitudinal.clear(); TrueTotalDedxLongitudinal.clear(); - TotalDedxLongitudinal.resize(3000, 0.0);; + TotalDedxLongitudinal.resize(3000, 0.0); + ; TrueTotalDedxLongitudinal.resize(3000, 0.0); primaryPDG.clear(); - primaryTrackLength.clear(); - primaryTrackLengthInTPC.clear(); + primaryTrackLength.clear(); + primaryTrackLengthInTPC.clear(); ProngEInLAr.clear(); ProngEInHadCal.clear(); ProngAngleToBeamDir.clear(); @@ -482,7 +480,6 @@ void AnalysisManager::BeginOfEvent() dir_coc_x.clear(); dir_coc_y.clear(); dir_coc_z.clear(); - } //--------------------------------------------------------------------- @@ -499,8 +496,7 @@ void AnalysisManager::EndOfEvent(const G4Event *event) //----------------------------------------------------------- - // FILL PRIMARIES/TRAJECTORIES TREE - FillPrimariesTree(event); + // FILL PARTICLES/TRAJECTORIES TREE FillParticlesTree(event); //----------------------------------------------------------- @@ -516,12 +512,13 @@ void AnalysisManager::EndOfEvent(const G4Event *event) //----------------------------------------------------------- - // FILL DETECTOR HITS - - if( fFlareSDs.size() > 0 ) FillFLArEOutput(); + // FILL DETECTOR HITS - if( fFaser2SDs.size() > 0 ) FillFASER2Output(); + if (fFlareSDs.size() > 0) + FillFLArEOutput(); + if (fFaser2SDs.size() > 0) + FillFASER2Output(); } //--------------------------------------------------------------------- @@ -529,148 +526,139 @@ void AnalysisManager::EndOfEvent(const G4Event *event) void AnalysisManager::FillEventTree(const G4Event *event) { - EventInformation* eventInfo = static_cast(event->GetUserInformation()); + EventInformation *eventInfo = static_cast(event->GetUserInformation()); eventInfo->Print(); auto metadata = eventInfo->GetEventMetadata(); - for(int i=0; iFill(); + int nVertexes = metadata.size(); + int nG4Vertexes = event->GetNumberOfPrimaryVertex(); + if (nVertexes != nG4Vertexes) + { + std::ostringstream oss; + oss << "Mismatch between vertexes in generator metadata vs G4Event (" << nVertexes << " vs " << nG4Vertexes << ")"; + G4Exception("AnalysisManager", "LogicError", JustWarning, oss.str().c_str()); } -} -//--------------------------------------------------------------------- -//--------------------------------------------------------------------- - -void AnalysisManager::FillPrimariesTree(const G4Event *event) -{ - nPrimaryVertex = event->GetNumberOfPrimaryVertex(); - G4cout << "\nNumber of primary vertices : " << nPrimaryVertex << G4endl; - /// loop over the vertices, and then over primary particles, - /// neutrino truth info from event generator. - for (G4int ivtx = 0; ivtx < event->GetNumberOfPrimaryVertex(); ++ivtx) + /// metadata info comes event generator + /// primaries come from G4event vertex + for (int ivtx = 0; ivtx < nVertexes; ivtx++) { - G4cout << "=== Vertex " << ivtx+1 << " of " << nPrimaryVertex << " -> " - << event->GetPrimaryVertex(ivtx)->GetNumberOfParticle() << " primaries ===" << G4endl; - for (G4int ipp = 0; ipp < event->GetPrimaryVertex(ivtx)->GetNumberOfParticle(); ++ipp) + vertexID = ivtx; + weight = metadata[ivtx].weight; + genType = metadata[ivtx].generatorType; + processName = metadata[ivtx].processName; + + vtxX = metadata[ivtx].x4.x(); + vtxY = metadata[ivtx].x4.y(); + vtxZ = metadata[ivtx].x4.z(); + vtxT = metadata[ivtx].x4.t(); + + initPDG = metadata[ivtx].pdg; + initPx = metadata[ivtx].p4.x(); + initPy = metadata[ivtx].p4.y(); + initPz = metadata[ivtx].p4.z(); + initE = metadata[ivtx].p4.e(); + initM = metadata[ivtx].mass; + initQ = metadata[ivtx].charge; + + fslPDG = metadata[ivtx].fsl_pdg; + tgtPDG = metadata[ivtx].tgt_pdg; + tgtZ = metadata[ivtx].tgt_Z; + tgtA = metadata[ivtx].tgt_A; + hitnucPDG = metadata[ivtx].hitnuc_pdg; + + intType = metadata[ivtx].intType; + scatteringType = metadata[ivtx].scatteringType; + xs = metadata[ivtx].xs; + Q2 = metadata[ivtx].Q2; + xBj = metadata[ivtx].xBj; + y = metadata[ivtx].y; + W = metadata[ivtx].W; + + nPrimaries = event->GetPrimaryVertex(ivtx)->GetNumberOfParticle(); + G4cout << "\n Vertex: " << ivtx << " -> number of primaries: " << nPrimaries << G4endl; + + primTID.clear(); + primPDG.clear(); + primPx.clear(); + primPy.clear(); + primPz.clear(); + primE.clear(); + + for (int ipp = 0; ipp < event->GetPrimaryVertex(ivtx)->GetNumberOfParticle(); ++ipp) { G4PrimaryParticle *primary_particle = event->GetPrimaryVertex(ivtx)->GetPrimary(ipp); if (primary_particle) { - - primVtxID = ivtx; - primTrackID = ipp + 1; // confirm matches track id? - - auto particleId = ActsFatras::Barcode(); - particleId.setVertexPrimary(ivtx); - particleId.setGeneration(0); - particleId.setSubParticle(0); - particleId.setParticle(primTrackID - 1); - - primParticleID = particleId.value(); - primPDG = primary_particle->GetPDGcode(); - primVx = event->GetPrimaryVertex(ivtx)->GetPosition().x(); - primVy = event->GetPrimaryVertex(ivtx)->GetPosition().y(); - primVz = event->GetPrimaryVertex(ivtx)->GetPosition().z(); - primVt = event->GetPrimaryVertex(ivtx)->GetT0(); - primPx = primary_particle->GetMomentum().x(); - primPy = primary_particle->GetMomentum().y(); - primPz = primary_particle->GetMomentum().z(); - primM = primary_particle->GetMass()/MeV; - primQ = primary_particle->GetCharge(); - - G4double energy = GetTotalEnergy(primPx, primPy, primPz, primM); - G4LorentzVector p4(primPx,primPy,primPz,energy); - primEta = p4.eta(); - primPhi = p4.phi(); - primPt = p4.perp(); - primP = p4.vect().mag(); - primE = energy; - primKE = energy - primM; + int tid = ipp + 1; // confirm matches track id? + int pdg = primary_particle->GetPDGcode(); + double px = primary_particle->GetMomentum().x(); + double py = primary_particle->GetMomentum().y(); + double pz = primary_particle->GetMomentum().z(); + double m = primary_particle->GetMass() / MeV; + double e = GetTotalEnergy(px, py, pz, m); + + primTID.push_back(tid); + primPDG.push_back(pdg); + primPx.push_back(px); + primPy.push_back(py); + primPz.push_back(pz); + primE.push_back(e); // store a copy as a FPFParticle for further processing - primaryIDs.push_back(primTrackID); //store to avoid duplicates - primaries.push_back(FPFParticle(primPDG, 0, - primTrackID, primaryIDs.size()-1, 1, - primM, - primVx, primVy, primVz, primVt, - primPx, primPy, primPz,energy)); - - G4cout << G4endl; - G4cout << "PrimaryParticleInfo: PDG code " << primPDG << G4endl - << "Particle unique ID : " << primTrackID << G4endl - << "Momentum : (" << primPx << ", " << primPy << ", " << primPz << ") MeV" << G4endl - << "Vertex : (" << primVx << ", " << primVy << ", " << primVz << ") mm" << G4endl; - - fPrim->Fill(); + // this accumulates across entire event (no vertex reset) + primaryIDs.push_back(tid); // store to avoid duplicates + primaries.push_back(FPFParticle(pdg, 0, + tid, primaryIDs.size() - 1, 1, + m, vtxX, vtxY, vtxZ, vtxT, + px, py, pz, e)); + + G4cout << "---" << G4endl; + G4cout << "PrimaryParticleInfo: PDG code " << pdg << G4endl + << "Particle TID : " << tid << G4endl + << "Momentum : (" << px << ", " << py << ", " << pz << ") MeV" << G4endl; } } + fEvt->Fill(); } - - G4cout << "\nNumber of primaries : " << primaryIDs.size() << G4endl; + G4cout << "\nTotal number of vertexes : " << nVertexes << G4endl; + G4cout << "Total number of primaries : " << primaryIDs.size() << G4endl; } //--------------------------------------------------------------------- //--------------------------------------------------------------------- -void AnalysisManager::FillParticlesTree(const G4Event* event) +void AnalysisManager::FillParticlesTree(const G4Event *event) { int count_tracks = 0; - G4cout << "==== Saving track information to tree ====" << G4endl; - auto trajectoryContainer = event->GetTrajectoryContainer(); + G4cout << "==== Saving track information to tree ====" << G4endl; + auto trajectoryContainer = event->GetTrajectoryContainer(); if (!trajectoryContainer) { G4cout << "No tracks found in the event!" << G4endl; return; } - for (size_t i = 0; i < trajectoryContainer->entries(); ++i) - { - auto trajectory = static_cast((*trajectoryContainer)[i]); + for (size_t i = 0; i < trajectoryContainer->entries(); ++i) + { + auto trajectory = static_cast((*trajectoryContainer)[i]); trackTID = trajectory->GetTrackID(); trackPID = trajectory->GetParentID(); - trackPDG = trajectory->GetPDGEncoding(); - trackKinE = trajectory->GetInitialKineticEnergy(); - trackNPoints = trajectory->GetPointEntries(); - count_tracks++; - for (size_t j = 0; j < trackNPoints; ++j) - { - G4ThreeVector pos = trajectory->GetPoint(j)->GetPosition(); - trackPointX.push_back( pos.x() ); - trackPointY.push_back( pos.y() ); - trackPointZ.push_back( pos.z() ); + trackPDG = trajectory->GetPDGEncoding(); + trackKinE = trajectory->GetInitialKineticEnergy(); + trackNPoints = trajectory->GetPointEntries(); + count_tracks++; + for (size_t j = 0; j < trackNPoints; ++j) + { + G4ThreeVector pos = trajectory->GetPoint(j)->GetPosition(); + trackPointX.push_back(pos.x()); + trackPointY.push_back(pos.y()); + trackPointZ.push_back(pos.z()); } fPar->Fill(); - trackPointX.clear(); + trackPointX.clear(); trackPointY.clear(); trackPointZ.clear(); } @@ -686,30 +674,30 @@ void AnalysisManager::FillFLArEOutput() // prepare the pixel map const double_t res_tpc[3] = {1, 5, 5}; // mm - + pm3D = new PixelMap3D(evtID, primaries.size(), initPDG, res_tpc); // boundaries in global coordinates - pm3D->SetPMBoundary(GeometricalParameters::Get()->GetFLArEPosition()/mm - - GeometricalParameters::Get()->GetTPCSizeXYZ()/mm/2, - GeometricalParameters::Get()->GetFLArEPosition()/mm + - GeometricalParameters::Get()->GetTPCSizeXYZ()/mm/2); + pm3D->SetPMBoundary(GeometricalParameters::Get()->GetFLArEPosition() / mm - + GeometricalParameters::Get()->GetTPCSizeXYZ() / mm / 2, + GeometricalParameters::Get()->GetFLArEPosition() / mm + + GeometricalParameters::Get()->GetTPCSizeXYZ() / mm / 2); pm3D->InitializePM(); // accumulate values per primary particle - ProngEInLAr.resize(primaries.size(),0.); - ProngEInHadCal.resize(primaries.size(),0.); - ShowerLength.resize(primaries.size(),0.); - ShowerLengthInLAr.resize(primaries.size(),0.); - ShowerWidth.resize(primaries.size(),0.); - ShowerWidthInLAr.resize(primaries.size(),0.); - primaryTrackLength.resize(primaries.size(),0.); - primaryTrackLengthInTPC.resize(primaries.size(),0.); + ProngEInLAr.resize(primaries.size(), 0.); + ProngEInHadCal.resize(primaries.size(), 0.); + ShowerLength.resize(primaries.size(), 0.); + ShowerLengthInLAr.resize(primaries.size(), 0.); + ShowerWidth.resize(primaries.size(), 0.); + ShowerWidthInLAr.resize(primaries.size(), 0.); + primaryTrackLength.resize(primaries.size(), 0.); + primaryTrackLengthInTPC.resize(primaries.size(), 0.); // loop over the detected FLArE sensitive volumes // for now, everything is lar_box int nLArHits = 0; int nHCALHits = 0; - for(const int sdId : fFlareSDs ) + for (const int sdId : fFlareSDs) { // Get and cast hit collection with LArBoxHits std::string sdName = fSDNamelist.at(sdId); @@ -719,7 +707,7 @@ void AnalysisManager::FillFLArEOutput() G4cout << "No hits recorded by " << sdName << G4endl; continue; } - + for (auto hit : *hitCollection->GetVector()) { flareTrackID = hit->GetTID(); @@ -745,8 +733,8 @@ void AnalysisManager::FillFLArEOutput() particleId.setParticle(flareTrackID); flareParticleID = particleId.value(); - double hit_position_xyz[3] = { flareX, flareY, flareZ }; - double vtx_xyz[3] = { primaries[0].Vx(), primaries[0].Vy(), primaries[0].Vz() }; + double hit_position_xyz[3] = {flareX, flareY, flareZ}; + double vtx_xyz[3] = {primaries[0].Vx(), primaries[0].Vy(), primaries[0].Vz()}; // which primary ancestor does this hit belong to? G4int whichPrim = GetTrackPrimaryAncestor(flareTrackID); @@ -755,24 +743,25 @@ void AnalysisManager::FillFLArEOutput() // pseudo reco: track/shower length and width primaryTrackLength[whichIndex] += hit->GetStepLength(); double ShowerP = primaries[whichIndex].P(); - double dsquare_hit_vtx = TMath::Power((flareX-primaries[whichIndex].Vx()),2)+ - TMath::Power((flareY-primaries[whichIndex].Vy()),2)+ - TMath::Power((flareZ-primaries[whichIndex].Vz()),2); - double product_hit_p = (flareX-primaries[whichIndex].Vx())*primaries[whichIndex].Px()+ - (flareY-primaries[whichIndex].Vy())*primaries[whichIndex].Py()+ - (flareZ-primaries[whichIndex].Vz())*primaries[whichIndex].Pz(); - double len_hit = TMath::Abs(product_hit_p)/ShowerP; - double width_hit = TMath::Sqrt((dsquare_hit_vtx - product_hit_p*product_hit_p/ShowerP/ShowerP)); - ShowerLength[whichIndex] = (ShowerLength[whichIndex]>len_hit) ? ShowerLength[whichIndex] : len_hit; - double weighted_width_hit = width_hit*flareEdep; - if (!std::isnan(weighted_width_hit)) ShowerWidth[whichIndex] += weighted_width_hit; + double dsquare_hit_vtx = TMath::Power((flareX - primaries[whichIndex].Vx()), 2) + + TMath::Power((flareY - primaries[whichIndex].Vy()), 2) + + TMath::Power((flareZ - primaries[whichIndex].Vz()), 2); + double product_hit_p = (flareX - primaries[whichIndex].Vx()) * primaries[whichIndex].Px() + + (flareY - primaries[whichIndex].Vy()) * primaries[whichIndex].Py() + + (flareZ - primaries[whichIndex].Vz()) * primaries[whichIndex].Pz(); + double len_hit = TMath::Abs(product_hit_p) / ShowerP; + double width_hit = TMath::Sqrt((dsquare_hit_vtx - product_hit_p * product_hit_p / ShowerP / ShowerP)); + ShowerLength[whichIndex] = (ShowerLength[whichIndex] > len_hit) ? ShowerLength[whichIndex] : len_hit; + double weighted_width_hit = width_hit * flareEdep; + if (!std::isnan(weighted_width_hit)) + ShowerWidth[whichIndex] += weighted_width_hit; if (sdName == "FLArEBoxSD/lar_box") { nLArHits++; fFLArEHits->Fill(); - if (fAddDiffusion == "toy") + if (fAddDiffusion == "toy") pm3D->FillEntryWithToyElectronTransportation(hit_position_xyz, vtx_xyz, flareEdep, whichIndex); else if (fAddDiffusion == "single") pm3D->FillEntryWithToySingleElectronTransportation(hit_position_xyz, vtx_xyz, flareEdep, whichIndex); @@ -783,19 +772,21 @@ void AnalysisManager::FillFLArEOutput() edepInLAr += flareEdep; ProngEInLAr[whichIndex] += flareEdep; primaryTrackLengthInTPC[whichIndex] += hit->GetStepLength(); - ShowerLengthInLAr[whichIndex] = (ShowerLengthInLAr[whichIndex]>len_hit) ? ShowerLengthInLAr[whichIndex] : len_hit; - if (!std::isnan(weighted_width_hit)) ShowerWidthInLAr[whichIndex] += weighted_width_hit; + ShowerLengthInLAr[whichIndex] = (ShowerLengthInLAr[whichIndex] > len_hit) ? ShowerLengthInLAr[whichIndex] : len_hit; + if (!std::isnan(weighted_width_hit)) + ShowerWidthInLAr[whichIndex] += weighted_width_hit; // accumulate dE/dx in LAr - float_t longitudinal_distance_to_vtx = ((flareX-vtx_xyz[0])*primaries[0].Px()+ - (flareY-vtx_xyz[1])*primaries[0].Py()+ - (flareZ-vtx_xyz[2])*primaries[0].Pz())/primaries[0].P(); - if (longitudinal_distance_to_vtx>=0 && longitudinal_distance_to_vtx<3000) // within 3000 mm + float_t longitudinal_distance_to_vtx = ((flareX - vtx_xyz[0]) * primaries[0].Px() + + (flareY - vtx_xyz[1]) * primaries[0].Py() + + (flareZ - vtx_xyz[2]) * primaries[0].Pz()) / + primaries[0].P(); + if (longitudinal_distance_to_vtx >= 0 && longitudinal_distance_to_vtx < 3000) // within 3000 mm TrueTotalDedxLongitudinal[Int_t(longitudinal_distance_to_vtx)] += hit->GetEdep(); } - else if (sdName == "FLArEHadCalXSD/lar_box" || + else if (sdName == "FLArEHadCalXSD/lar_box" || sdName == "FLArEMuonFinderXSD/lar_box" || - sdName == "FLArEBabyMINDHorBarSD/lar_box" ) + sdName == "FLArEBabyMINDHorBarSD/lar_box") { // specify this is an ZX hit nHCALHits++; @@ -806,11 +797,11 @@ void AnalysisManager::FillFLArEOutput() edepInHCALX += flareEdep; ProngEInHadCal[whichIndex] += flareEdep; } - else if (sdName == "FLArEHadCalYSD/lar_box" || + else if (sdName == "FLArEHadCalYSD/lar_box" || sdName == "FLArEMuonFinderYSD/lar_box" || sdName == "FLArEBabyMINDVerBarSD/lar_box" || sdName == "FLArEHadAbsorbSD/lar_box" || - sdName == "FLArEMuonFinderAbsorbSD/lar_box" ) + sdName == "FLArEMuonFinderAbsorbSD/lar_box") { nHCALHits++; flareIsZX = false; @@ -823,14 +814,16 @@ void AnalysisManager::FillFLArEOutput() } } - if (fSave2DEvd) pm3D->Write2DPMToFile(fFile,fFLArEDir); + if (fSave2DEvd) + pm3D->Write2DPMToFile(fFile, fFLArEDir); pm3D->Process3DPM(fH5file, initPDG, fslPDG, intType, scatteringType, initE, fSave3DEvd); sparseFractionMem = pm3D->GetSparseFractionMem(); sparseFractionBins = pm3D->GetSparseFractionBins(); - if( fSavePseudoReco ) FillFLArEPseudoReco(); + if (fSavePseudoReco) + FillFLArEPseudoReco(); G4cout << "Total FLArE recorded hits: " << nLArHits << G4endl; G4cout << "Total FLArE HCAL recorded hits: " << nHCALHits << G4endl; @@ -847,7 +840,7 @@ void AnalysisManager::FillFASER2Output() // loop over the detected FASER2 sensitive volumes int nHits = 0; - for(const int sdId : fFaser2SDs ) + for (const int sdId : fFaser2SDs) { // Get and cast hit collection with LArBoxHits std::string sdName = fSDNamelist.at(sdId); @@ -857,7 +850,7 @@ void AnalysisManager::FillFASER2Output() G4cout << "No hits recorded by " << sdName << G4endl; continue; } - + std::map sub_part_map{}; for (auto hit : *hitCollection->GetVector()) { @@ -923,7 +916,8 @@ void AnalysisManager::FillFASER2Output() isDuplicate = true; } } - if (isDuplicate) continue; // Skip this particle if it's already been added + if (isDuplicate) + continue; // Skip this particle if it's already been added ActsParticlesParticleId.push_back(particleId.value()); ActsParticlesParticleType.push_back(hit->GetPDGID()); @@ -979,26 +973,26 @@ void AnalysisManager::FillFLArEPseudoReco() << std::setw(12) << "Pz" << std::endl; nprimaries = primaries.size(); - primaryPDG.resize(nprimaries,0.); - ProngAngleToBeamDir.resize(nprimaries,0.); - ProngAvgdEdx.resize(nprimaries,0.); - ProngAvgdEdxInLAr.resize(nprimaries,0.); - dir_pol_x.resize(nprimaries,0.); - dir_pol_y.resize(nprimaries,0.); - dir_pol_z.resize(nprimaries,0.); - dir_coc_x.resize(nprimaries,0.); - dir_coc_y.resize(nprimaries,0.); - dir_coc_z.resize(nprimaries,0.); - - for (auto iPrim : primaryIDs ) + primaryPDG.resize(nprimaries, 0.); + ProngAngleToBeamDir.resize(nprimaries, 0.); + ProngAvgdEdx.resize(nprimaries, 0.); + ProngAvgdEdxInLAr.resize(nprimaries, 0.); + dir_pol_x.resize(nprimaries, 0.); + dir_pol_y.resize(nprimaries, 0.); + dir_pol_z.resize(nprimaries, 0.); + dir_coc_x.resize(nprimaries, 0.); + dir_coc_y.resize(nprimaries, 0.); + dir_coc_z.resize(nprimaries, 0.); + + for (auto iPrim : primaryIDs) { // trackIDs go from 1 to N, you need index: 0 to N-1 G4int iiPrim = iPrim - 1; primaryPDG[iiPrim] = primaries[iiPrim].PDG(); - float_t totProngE = ProngEInLAr[iiPrim]+ProngEInHadCal[iiPrim]; - if (totProngE>0) + float_t totProngE = ProngEInLAr[iiPrim] + ProngEInHadCal[iiPrim]; + if (totProngE > 0) { ShowerWidth[iiPrim] = ShowerWidth[iiPrim] / totProngE; } @@ -1012,8 +1006,8 @@ void AnalysisManager::FillFLArEPseudoReco() ProngAngleToBeamDir[iiPrim] = TMath::ACos(costheta); ProngAvgdEdx[iiPrim] = (ProngEInLAr[iiPrim] + - ProngEInHadCal[iiPrim]) / - ShowerLength[iiPrim]; + ProngEInHadCal[iiPrim]) / + ShowerLength[iiPrim]; ProngAvgdEdxInLAr[iiPrim] = ProngEInLAr[iiPrim] / ShowerLengthInLAr[iiPrim]; std::cout << std::setiosflags(std::ios::fixed) << std::setprecision(3); @@ -1030,7 +1024,7 @@ void AnalysisManager::FillFLArEPseudoReco() } slid::ShowerLID *shwlid = new slid::ShowerLID(pm3D->Get3DPixelMap(), - initX, initY, initZ, 0., 0., 1.); + vtxX, vtxY, vtxZ, 0., 0., 1.); Double_t *ptr_dedx = shwlid->GetTotalDedxLongitudinal(); std::copy(ptr_dedx, ptr_dedx + 3000, TotalDedxLongitudinal.begin()); @@ -1041,7 +1035,7 @@ void AnalysisManager::FillFLArEPseudoReco() pm3D->Get2DPixelMapZY(iPrim + 1), pm3D->Get2DVtxPixelMapZX(iPrim + 1), pm3D->Get2DVtxPixelMapZY(iPrim + 1), - initX, initY, initZ, + vtxX, vtxY, vtxZ, primaries[iPrim].Vx(), primaries[iPrim].Vy(), primaries[iPrim].Vz()); dir_pol_x[iPrim] = linFit->GetDir().X(); dir_pol_y[iPrim] = linFit->GetDir().Y(); From 8acb77e4814d387219b249f65dedeb5718c15a48 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Thu, 18 Dec 2025 00:13:33 -0500 Subject: [PATCH 3/8] set up unique particle/traj tree --- README.md | 1 + include/AnalysisManager.hh | 77 ++++---- include/AnalysisManagerMessenger.hh | 5 +- include/EventInformation.hh | 1 + include/FPFTrajectory.hh | 8 +- src/AnalysisManager.cc | 278 ++++++++++------------------ src/AnalysisManagerMessenger.cc | 11 +- src/EventInformation.cc | 15 +- src/FPFTrajectory.cc | 21 ++- 9 files changed, 173 insertions(+), 244 deletions(-) diff --git a/README.md b/README.md index 875906f..af95862 100644 --- a/README.md +++ b/README.md @@ -187,6 +187,7 @@ Older versions of FORESEE output events in the HepMC2 format. To run over HepMC2 |:--|:--| |/out/fileName | option for AnalysisManagerMessenger, set name of the file saving all analysis variables| |/out/saveTrack | if `true` save full trajectories, `false` by default | +|/out/parKinECut | set particle kinetic energy threshold for saving in output tree | |/out/flare/enableOutput | if `true` save FLArE output, `true` by default | |/out/flare/save3DEvd | if `true` save 3D spatial distribution of energy deposition, `false` by default | |/out/flare/save2DEvd | if `true` save 2D spatial distribution of energy deposition, `false` by default | diff --git a/include/AnalysisManager.hh b/include/AnalysisManager.hh index 2841341..c410939 100644 --- a/include/AnalysisManager.hh +++ b/include/AnalysisManager.hh @@ -35,6 +35,7 @@ class AnalysisManager { // functions for controlling from the configuration file void setFileName(std::string val) { fFilename = val; } void saveTrack(G4bool val) { fSaveTrack = val; } + void parKinECut(G4double val) { fParKinECut = val; } void saveActs(G4bool val) { fSaveActs = val; } void savePseudoReco(G4bool val) { fSavePseudoReco = val; } void addDiffusion(G4String val) { fAddDiffusion = val; } @@ -74,6 +75,7 @@ class AnalysisManager { AnalysisManagerMessenger* fMessenger; G4bool fSaveTrack; + G4double fParKinECut; G4bool fSave3DEvd; G4bool fSave2DEvd; G4bool fSavePseudoReco; @@ -108,36 +110,33 @@ class AnalysisManager { TDirectory* fFASER2Dir; TTree* fActsHitsTree; - TTree* fActsParticlesTree; - // track to primary ancestor + // track to primary ancestor (track id to track id) std::map trackToPrimaryAncestor; + // map trackID to its barcode + std::map trackIDtoParticleID; + //--------------------------------------------------- // OUTPUT VARIABLES FOR COMMON TREES //--------------------------------------------------- // Output variables for EVENT tree - // general vertex metadata G4int evtID; G4int vertexID; double weight; std::string genType; std::string processName; - // vertex position double vtxX, vtxY, vtxZ, vtxT; - // initiator info (incoming particle) int initPDG; double initPx, initPy, initPz, initE; double initM; double initQ; - // target info int fslPDG; int tgtPDG; int tgtA; int tgtZ; int hitnucPDG; - // interaction info int intType; int scatteringType; double xs; @@ -145,7 +144,6 @@ class AnalysisManager { double xBj; double y; double W; - // primaries from vertex int nPrimaries; std::vector primTID; std::vector primPDG; @@ -155,15 +153,31 @@ class AnalysisManager { std::vector primE; //--------------------------------------------------- - // Output variables for PARTICLES tree - int trackTID; - int trackPID; - int trackPDG; - double trackKinE; - int trackNPoints; - std::vector trackPointX; - std::vector trackPointY; - std::vector trackPointZ; + // Output variables for PARTICLES/TRAJECTORIES tree + + std::uint64_t particle_id; //barcode + int particle_TID; + int particle_PID; + int particle_PDG; + std::string particle_process; + float particle_vx; + float particle_vy; + float particle_vz; + float particle_vt; + float particle_px; + float particle_py; + float particle_pz; + float particle_m; + float particle_q; + float particle_eta; + float particle_phi; + float particle_pt; + float particle_p; + float particle_ke; + int traj_Npoints; + std::vector traj_pointX; + std::vector traj_pointY; + std::vector traj_pointZ; //--------------------------------------------------- // OUTPUT VARIABLES FOR FLArE TREES @@ -243,35 +257,6 @@ class AnalysisManager { UInt_t ActsHitsApproachID; UInt_t ActsHitsSensitiveID; - // Acts Particle Information - need the truth info on the particles in order to do the truth tracking - std::vector ActsParticlesParticleId; - std::vector ActsParticlesParticleType; - std::vector ActsParticlesProcess; - std::vector ActsParticlesVx; - std::vector ActsParticlesVy; - std::vector ActsParticlesVz; - std::vector ActsParticlesVt; - std::vector ActsParticlesPx; - std::vector ActsParticlesPy; - std::vector ActsParticlesPz; - std::vector ActsParticlesM; - std::vector ActsParticlesQ; - std::vector ActsParticlesEta; - std::vector ActsParticlesPhi; - std::vector ActsParticlesPt; - std::vector ActsParticlesP; - std::vector ActsParticlesVertexPrimary; - std::vector ActsParticlesVertexSecondary; - std::vector ActsParticlesParticle; - - std::vector ActsParticlesGeneration; - std::vector ActsParticlesSubParticle; - std::vector ActsParticlesELoss; - std::vector ActsParticlesPathInX0; - std::vector ActsParticlesPathInL0; - std::vector ActsParticlesNumberOfHits; - std::vector ActsParticlesOutcome; - }; #endif diff --git a/include/AnalysisManagerMessenger.hh b/include/AnalysisManagerMessenger.hh index ba305cb..47bc673 100644 --- a/include/AnalysisManagerMessenger.hh +++ b/include/AnalysisManagerMessenger.hh @@ -38,6 +38,8 @@ class G4UIcommand; class G4UIcmdWithAString; class G4UIcmdWithABool; class G4UIcmdWithAnInteger; +class G4UIcmdWithADoubleAndUnit; + //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -57,7 +59,8 @@ class AnalysisManagerMessenger: public G4UImessenger G4UIdirectory* fOutDir; G4UIcmdWithAString* fFileCmd; G4UIcmdWithABool* fSaveTrackCmd; - + G4UIcmdWithADoubleAndUnit* fParticleKinECutCmd; + G4UIdirectory* fFLArEDir; G4UIcmdWithABool* fEnableFLArEOutCmd; G4UIcmdWithABool* fSave3DEvdCmd; diff --git a/include/EventInformation.hh b/include/EventInformation.hh index acb08cc..9e35682 100644 --- a/include/EventInformation.hh +++ b/include/EventInformation.hh @@ -23,6 +23,7 @@ class EventInformation : public G4VUserEventInformation /// Prints the information about the event. virtual void Print() const; + virtual void Print(int i) const; private: /// Set of vertex metadata diff --git a/include/FPFTrajectory.hh b/include/FPFTrajectory.hh index ad0288b..4eded8f 100644 --- a/include/FPFTrajectory.hh +++ b/include/FPFTrajectory.hh @@ -9,6 +9,7 @@ #include "G4TrajectoryPoint.hh" #include "G4Track.hh" #include "G4Step.hh" +#include "G4LorentzVector.hh" /// FPFTrajectory is a custom trajectory class that allows to disable /// storing the full trajectory point collections (memory intensive) @@ -38,9 +39,10 @@ public: inline G4String GetParticleName() const { return fParticleName; } inline G4double GetCharge() const { return fPDGCharge; } inline G4int GetPDGEncoding() const { return fPDGEncoding; } - inline G4ThreeVector GetInitialMomentum() const { return fInitialMomentum; } - G4double GetInitialKineticEnergy() const; inline G4String GetProcessName() const { return fProcessName; } + inline G4LorentzVector GetInitialP4() const {return fInitialP4; } + G4ThreeVector GetInitialMomentum() const; + G4double GetInitialKineticEnergy() const; virtual void ShowTrajectory(std::ostream& os=G4cout) const; virtual void DrawTrajectory(G4int i_mode = 0) const; @@ -66,7 +68,7 @@ private: G4int fPDGEncoding; G4double fPDGCharge; G4String fParticleName; - G4ThreeVector fInitialMomentum; + G4LorentzVector fInitialP4; // custom additions G4String fProcessName; // creator process diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index 2b1c2c4..8736b8c 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -71,9 +71,9 @@ AnalysisManager::AnalysisManager() fFLArEHCALHits = nullptr; fFLArEPseudoReco = nullptr; fActsHitsTree = nullptr; - fActsParticlesTree = nullptr; fSaveTrack = false; + fParKinECut = 50.0*keV; fSave3DEvd = false; fSave2DEvd = false; fSavePseudoReco = false; @@ -90,18 +90,15 @@ AnalysisManager::~AnalysisManager() {} void AnalysisManager::bookEvtTree() { fEvt = new TTree("event", "event info"); - // vertex generation metadata - fEvt->Branch("evtID", &evtID, "evtID/I"); - fEvt->Branch("vtxID", &vertexID, "vtxID/I"); + fEvt->Branch("event_id", &evtID, "event_id/I"); + fEvt->Branch("vertex_id", &vertexID, "vertex_id/I"); fEvt->Branch("weight", &weight, "weight/D"); fEvt->Branch("generatorType", &genType); fEvt->Branch("processName", &processName); - // vertex position fEvt->Branch("vtxX", &vtxX, "vtxX/D"); fEvt->Branch("vtxY", &vtxY, "vtxY/D"); fEvt->Branch("vtxZ", &vtxZ, "vtxZ/D"); fEvt->Branch("vtxT", &vtxT, "vtxT/D"); - // initiator info (incoming particle) fEvt->Branch("initPDG", &initPDG, "initPDG/I"); fEvt->Branch("initPx", &initPx, "initPx/D"); fEvt->Branch("initPy", &initPy, "initPy/D"); @@ -109,13 +106,11 @@ void AnalysisManager::bookEvtTree() fEvt->Branch("initE", &initE, "initE/D"); fEvt->Branch("initM", &initM, "initM/D"); fEvt->Branch("initQ", &initQ, "initQ/D"); - // target info fEvt->Branch("tgtPDG", &tgtPDG, "tgtPDG/I"); fEvt->Branch("tgtA", &tgtA, "tgtA/I"); fEvt->Branch("tgtZ", &tgtZ, "tgtZ/I"); fEvt->Branch("hitnucPDG", &hitnucPDG, "hitnucPDG/I"); fEvt->Branch("fslPDG", &fslPDG, "fslPDG/I"); - // interaction info fEvt->Branch("intType", &intType, "intType/I"); fEvt->Branch("scatteringType", &scatteringType, "scatteringType/I"); fEvt->Branch("xs", &xs, "xs/D"); @@ -123,7 +118,6 @@ void AnalysisManager::bookEvtTree() fEvt->Branch("xBj", &xBj, "xBj/D"); fEvt->Branch("y", &y, "y/D"); fEvt->Branch("W", &W, "W/D"); - // primaries from vertex fEvt->Branch("nPrimaries", &nPrimaries, "nPrimaries/I"); fEvt->Branch("primTID", &primTID); fEvt->Branch("primPDG", &primPDG); @@ -136,19 +130,32 @@ void AnalysisManager::bookEvtTree() void AnalysisManager::bookParTree() { fPar = new TTree("particles", "particle info"); - fPar->Branch("evtID", &evtID, "evtID/I"); - fPar->Branch("trackTID", &trackTID, "trackTID/I"); - fPar->Branch("trackPID", &trackPID, "trackPID/I"); - fPar->Branch("trackPDG", &trackPDG, "trackPDG/I"); - fPar->Branch("trackKinE", &trackKinE, "trackKinE/D"); - fPar->Branch("trackNPoints", &trackNPoints, "trackNPoints/I"); - - // if saving full trajectory, add vector of trajectory points - if (fSaveTrack) + fPar->Branch("event_id", &evtID, "event_id/I"); + fPar->Branch("particle_id", &particle_id, "particle_id/I"); + fPar->Branch("particle_pdg", &particle_PDG, "particle_pdg/I"); + fPar->Branch("track_id", &particle_TID, "track_id/I"); + fPar->Branch("parent_id", &particle_PID, "parent_id/I"); + fPar->Branch("process", &particle_process); + fPar->Branch("vx", &particle_vx, "vx/F"); + fPar->Branch("vy", &particle_vy, "vy/F"); + fPar->Branch("vz", &particle_vz, "vz/F"); + fPar->Branch("vt", &particle_vt, "vt/F"); + fPar->Branch("px", &particle_px, "px/F"); + fPar->Branch("py", &particle_py, "py/F"); + fPar->Branch("pz", &particle_pz, "pz/F"); + fPar->Branch("m", &particle_m, "m/F"); + fPar->Branch("q", &particle_q, "q/F"); + fPar->Branch("eta", &particle_eta, "eta/F"); + fPar->Branch("phi", &particle_phi, "phi/F"); + fPar->Branch("pt", &particle_pt, "pt/F"); + fPar->Branch("p", &particle_p, "p/F"); + fPar->Branch("ke", &particle_ke, "ke/F"); + if (fSaveTrack) // if saving full trajectory, add vector of trajectory points { - fPar->Branch("trackPointX", &trackPointX); - fPar->Branch("trackPointY", &trackPointY); - fPar->Branch("trackPointZ", &trackPointZ); + fPar->Branch("traj_Npoints", &traj_Npoints, "traj_Npoints/I"); + fPar->Branch("traj_pointX", &traj_pointX); + fPar->Branch("traj_pointY", &traj_pointY); + fPar->Branch("traj_pointZ", &traj_pointZ); } } @@ -272,36 +279,6 @@ void AnalysisManager::bookFASER2Trees() fActsHitsTree->Branch("approach_id", &ActsHitsApproachID, "approach_id/i"); fActsHitsTree->Branch("sensitive_id", &ActsHitsSensitiveID, "sensitive_id/i"); - //* Acts truth particle tree - fActsParticlesTree = new TTree("particles", "ActsParticlesTree"); - fActsParticlesTree->Branch("event_id", &ActsHitsEventID, "event_id/i"); - fActsParticlesTree->Branch("particle_id", &ActsParticlesParticleId); - fActsParticlesTree->Branch("particle_type", &ActsParticlesParticleType); - fActsParticlesTree->Branch("process", &ActsParticlesProcess); - fActsParticlesTree->Branch("vx", &ActsParticlesVx); - fActsParticlesTree->Branch("vy", &ActsParticlesVy); - fActsParticlesTree->Branch("vz", &ActsParticlesVz); - fActsParticlesTree->Branch("vt", &ActsParticlesVt); - fActsParticlesTree->Branch("px", &ActsParticlesPx); - fActsParticlesTree->Branch("py", &ActsParticlesPy); - fActsParticlesTree->Branch("pz", &ActsParticlesPz); - fActsParticlesTree->Branch("m", &ActsParticlesM); - fActsParticlesTree->Branch("q", &ActsParticlesQ); - fActsParticlesTree->Branch("eta", &ActsParticlesEta); - fActsParticlesTree->Branch("phi", &ActsParticlesPhi); - fActsParticlesTree->Branch("pt", &ActsParticlesPt); - fActsParticlesTree->Branch("p", &ActsParticlesP); - fActsParticlesTree->Branch("vertex_primary", &ActsParticlesVertexPrimary); - fActsParticlesTree->Branch("vertex_secondary", &ActsParticlesVertexSecondary); - fActsParticlesTree->Branch("particle", &ActsParticlesParticle); - fActsParticlesTree->Branch("generation", &ActsParticlesGeneration); - fActsParticlesTree->Branch("sub_particle", &ActsParticlesSubParticle); - fActsParticlesTree->Branch("e_loss", &ActsParticlesELoss); - fActsParticlesTree->Branch("total_x0", &ActsParticlesPathInX0); - fActsParticlesTree->Branch("total_l0", &ActsParticlesPathInL0); - fActsParticlesTree->Branch("number_of_hits", &ActsParticlesNumberOfHits); - fActsParticlesTree->Branch("outcome", &ActsParticlesOutcome); - fFile->cd(); } @@ -387,7 +364,6 @@ void AnalysisManager::EndOfRun() { fFile->cd(fFASER2Dir->GetName()); fActsHitsTree->Write(); - fActsParticlesTree->Write(); fFile->cd(); // go back to top } @@ -402,11 +378,10 @@ void AnalysisManager::BeginOfEvent() // reset vectors that need to be cleared for a new event // only reset arrays or vectors, tipically no need for other defaults - primaries.clear(); - primaryIDs.clear(); - // primaries in event tree // (actually need reset at every vertex) + primaries.clear(); + primaryIDs.clear(); primTID.clear(); primPDG.clear(); primPx.clear(); @@ -417,38 +392,14 @@ void AnalysisManager::BeginOfEvent() // track ID to primary ancestor association trackToPrimaryAncestor.clear(); + // track ID to particle ID association + trackIDtoParticleID.clear(); + // trajectory points (if enabled) // (actually need reset at every track) - trackPointX.clear(); - trackPointY.clear(); - trackPointZ.clear(); - - ActsParticlesParticleId.clear(); - ActsParticlesParticleType.clear(); - ActsParticlesProcess.clear(); - ActsParticlesVx.clear(); - ActsParticlesVy.clear(); - ActsParticlesVz.clear(); - ActsParticlesVt.clear(); - ActsParticlesPx.clear(); - ActsParticlesPy.clear(); - ActsParticlesPz.clear(); - ActsParticlesM.clear(); - ActsParticlesQ.clear(); - ActsParticlesEta.clear(); - ActsParticlesPhi.clear(); - ActsParticlesPt.clear(); - ActsParticlesP.clear(); - ActsParticlesVertexPrimary.clear(); - ActsParticlesVertexSecondary.clear(); - ActsParticlesParticle.clear(); - ActsParticlesGeneration.clear(); - ActsParticlesSubParticle.clear(); - ActsParticlesELoss.clear(); - ActsParticlesPathInX0.clear(); - ActsParticlesPathInL0.clear(); - ActsParticlesNumberOfHits.clear(); - ActsParticlesOutcome.clear(); + traj_pointX.clear(); + traj_pointY.clear(); + traj_pointZ.clear(); // these are used to accumulate // so need to be reset @@ -491,15 +442,16 @@ void AnalysisManager::EndOfEvent(const G4Event *event) /// evtID evtID = event->GetEventID(); + //----------------------------------------------------------- // FILL EVENT TREE FillEventTree(event); //----------------------------------------------------------- - // FILL PARTICLES/TRAJECTORIES TREE FillParticlesTree(event); //----------------------------------------------------------- + // FILL DETECTOR HITS // Get the hit collections // If there is no hit collection, there is nothing to be done @@ -510,10 +462,6 @@ void AnalysisManager::EndOfEvent(const G4Event *event) return; } - //----------------------------------------------------------- - - // FILL DETECTOR HITS - if (fFlareSDs.size() > 0) FillFLArEOutput(); @@ -527,7 +475,6 @@ void AnalysisManager::EndOfEvent(const G4Event *event) void AnalysisManager::FillEventTree(const G4Event *event) { EventInformation *eventInfo = static_cast(event->GetUserInformation()); - eventInfo->Print(); auto metadata = eventInfo->GetEventMetadata(); int nVertexes = metadata.size(); @@ -544,16 +491,16 @@ void AnalysisManager::FillEventTree(const G4Event *event) /// primaries come from G4event vertex for (int ivtx = 0; ivtx < nVertexes; ivtx++) { + eventInfo->Print(ivtx); + vertexID = ivtx; weight = metadata[ivtx].weight; genType = metadata[ivtx].generatorType; processName = metadata[ivtx].processName; - vtxX = metadata[ivtx].x4.x(); vtxY = metadata[ivtx].x4.y(); vtxZ = metadata[ivtx].x4.z(); vtxT = metadata[ivtx].x4.t(); - initPDG = metadata[ivtx].pdg; initPx = metadata[ivtx].p4.x(); initPy = metadata[ivtx].p4.y(); @@ -561,13 +508,11 @@ void AnalysisManager::FillEventTree(const G4Event *event) initE = metadata[ivtx].p4.e(); initM = metadata[ivtx].mass; initQ = metadata[ivtx].charge; - fslPDG = metadata[ivtx].fsl_pdg; tgtPDG = metadata[ivtx].tgt_pdg; tgtZ = metadata[ivtx].tgt_Z; tgtA = metadata[ivtx].tgt_A; hitnucPDG = metadata[ivtx].hitnuc_pdg; - intType = metadata[ivtx].intType; scatteringType = metadata[ivtx].scatteringType; xs = metadata[ivtx].xs; @@ -577,8 +522,6 @@ void AnalysisManager::FillEventTree(const G4Event *event) W = metadata[ivtx].W; nPrimaries = event->GetPrimaryVertex(ivtx)->GetNumberOfParticle(); - G4cout << "\n Vertex: " << ivtx << " -> number of primaries: " << nPrimaries << G4endl; - primTID.clear(); primPDG.clear(); primPx.clear(); @@ -586,6 +529,8 @@ void AnalysisManager::FillEventTree(const G4Event *event) primPz.clear(); primE.clear(); + G4cout << "\n-> number of primaries: " << nPrimaries << G4endl; + for (int ipp = 0; ipp < event->GetPrimaryVertex(ivtx)->GetNumberOfParticle(); ++ipp) { G4PrimaryParticle *primary_particle = event->GetPrimaryVertex(ivtx)->GetPrimary(ipp); @@ -614,10 +559,7 @@ void AnalysisManager::FillEventTree(const G4Event *event) m, vtxX, vtxY, vtxZ, vtxT, px, py, pz, e)); - G4cout << "---" << G4endl; - G4cout << "PrimaryParticleInfo: PDG code " << pdg << G4endl - << "Particle TID : " << tid << G4endl - << "Momentum : (" << px << ", " << py << ", " << pz << ") MeV" << G4endl; + G4cout << "TID: " << tid << ", PDG: " << pdg << ", p=(" << px << ", " << py << ", " << pz << ") MeV" << G4endl; } } fEvt->Fill(); @@ -632,37 +574,73 @@ void AnalysisManager::FillEventTree(const G4Event *event) void AnalysisManager::FillParticlesTree(const G4Event *event) { int count_tracks = 0; - - G4cout << "==== Saving track information to tree ====" << G4endl; + G4cout << "==== Saving particle information to tree ====" << G4endl; auto trajectoryContainer = event->GetTrajectoryContainer(); if (!trajectoryContainer) { - G4cout << "No tracks found in the event!" << G4endl; + G4cout << "No trajectories found in the event!" << G4endl; return; } + std::map sub_part_map{}; for (size_t i = 0; i < trajectoryContainer->entries(); ++i) { auto trajectory = static_cast((*trajectoryContainer)[i]); - trackTID = trajectory->GetTrackID(); - trackPID = trajectory->GetParentID(); - trackPDG = trajectory->GetPDGEncoding(); - trackKinE = trajectory->GetInitialKineticEnergy(); - trackNPoints = trajectory->GetPointEntries(); - count_tracks++; - for (size_t j = 0; j < trackNPoints; ++j) + + // do not save if below cut threshold + particle_ke = trajectory->GetInitialKineticEnergy(); + if(particle_ke < fParKinECut) continue; + + particle_TID = trajectory->GetTrackID(); + particle_PID = trajectory->GetParentID(); + particle_PDG = trajectory->GetPDGEncoding(); + particle_process = trajectory->GetProcessName(); + + auto particleId = ActsFatras::Barcode(); + particleId.setVertexPrimary(1); + particleId.setVertexSecondary(0); + particleId.setParticle(particle_TID - 1); // The track ID is the primary particle index plus one + particleId.setGeneration(particle_PID); + sub_part_map.try_emplace(particle_TID- 1, sub_part_map.size()); + // This is a fudge - assumes that the secondary particles are always sub-particles of the primary particle + particleId.setSubParticle(particle_PID == 0 ? 0 : sub_part_map[particle_TID - 1]); + particle_id = particleId.value(); + + // store in map for use in other trees! + trackIDtoParticleID[particle_TID] = particle_id; + + // first point is vertex (always stored even if traj is disabled) + particle_vx = trajectory->GetPoint(0)->GetPosition().x(); + particle_vy = trajectory->GetPoint(0)->GetPosition().y(); + particle_vz = trajectory->GetPoint(0)->GetPosition().z(); + particle_vt = 0; + particle_px = trajectory->GetInitialP4().px(); + particle_py = trajectory->GetInitialP4().py(); + particle_pz = trajectory->GetInitialP4().pz(); + particle_m = trajectory->GetInitialP4().m(); + particle_q = trajectory->GetCharge(); + particle_eta = trajectory->GetInitialP4().eta(); + particle_phi = trajectory->GetInitialP4().phi(); + particle_pt = trajectory->GetInitialP4().vect().perp(); + particle_p = trajectory->GetInitialP4().vect().mag(); + + traj_Npoints = trajectory->GetPointEntries(); + traj_pointX.clear(); + traj_pointY.clear(); + traj_pointZ.clear(); + + for (size_t j = 0; j < traj_Npoints; ++j) { G4ThreeVector pos = trajectory->GetPoint(j)->GetPosition(); - trackPointX.push_back(pos.x()); - trackPointY.push_back(pos.y()); - trackPointZ.push_back(pos.z()); + traj_pointX.push_back(pos.x()); + traj_pointY.push_back(pos.y()); + traj_pointZ.push_back(pos.z()); } + count_tracks++; fPar->Fill(); - trackPointX.clear(); - trackPointY.clear(); - trackPointZ.clear(); } - G4cout << "Total number of recorded track: " << count_tracks << std::endl; + + G4cout << "Total number of recorded particles: " << count_tracks << std::endl; } //--------------------------------------------------------------------- @@ -725,13 +703,7 @@ void AnalysisManager::FillFLArEOutput() flareDeltaPy = hit->GetDeltaMom().y(); flareDeltaPz = hit->GetDeltaMom().z(); flareEdep = hit->GetEdep(); - - auto particleId = ActsFatras::Barcode(); - particleId.setVertexPrimary(1); // fix this value - particleId.setGeneration(flareParentID); - particleId.setSubParticle(0); - particleId.setParticle(flareTrackID); - flareParticleID = particleId.value(); + flareParticleID = trackIDtoParticleID.at(flareTrackID); double hit_position_xyz[3] = {flareX, flareY, flareZ}; double vtx_xyz[3] = {primaries[0].Vx(), primaries[0].Vy(), primaries[0].Vz()}; @@ -851,7 +823,6 @@ void AnalysisManager::FillFASER2Output() continue; } - std::map sub_part_map{}; for (auto hit : *hitCollection->GetVector()) { if (hit->GetCharge() == 0) @@ -870,19 +841,7 @@ void AnalysisManager::FillFASER2Output() ActsHitsGeometryID = 0; int hitID = hit->GetTrackID(); - int nPrimaries = ActsParticlesParticleId.size(); - - auto particleId = ActsFatras::Barcode(); - particleId.setVertexPrimary(1); - particleId.setVertexSecondary(0); - particleId.setParticle(hit->GetTrackID() - 1); // The track ID is the primary particle index plus one - particleId.setGeneration(hit->GetParentID()); - - sub_part_map.try_emplace(hit->GetTrackID() - 1, sub_part_map.size()); - - // This is a fudge - assumes that that the secondary particles are always sub-particles of the primary particle - particleId.setSubParticle(hit->GetParentID() == 0 ? 0 : sub_part_map[hit->GetTrackID() - 1]); - ActsHitsParticleID = particleId.value(); + ActsHitsParticleID = trackIDtoParticleID.at(hitID); ActsHitsX = hit->GetX(); ActsHitsY = hit->GetY(); @@ -907,49 +866,8 @@ void AnalysisManager::FillFASER2Output() ActsHitsSensitiveID = 1; fActsHitsTree->Fill(); - // Now fill the Acts particles tree - bool isDuplicate = false; - for (const auto &id : ActsParticlesParticleId) - { - if (id == particleId.value()) - { - isDuplicate = true; - } - } - if (isDuplicate) - continue; // Skip this particle if it's already been added - - ActsParticlesParticleId.push_back(particleId.value()); - ActsParticlesParticleType.push_back(hit->GetPDGID()); - ActsParticlesProcess.push_back(0); - ActsParticlesVx.push_back(hit->GetTrackVertex().x()); - ActsParticlesVy.push_back(hit->GetTrackVertex().y()); - ActsParticlesVz.push_back(hit->GetTrackVertex().z()); - ActsParticlesVt.push_back(0); - ActsParticlesPx.push_back(hit->GetTrackP4().px()); - ActsParticlesPy.push_back(hit->GetTrackP4().py()); - ActsParticlesPz.push_back(hit->GetTrackP4().pz()); - ActsParticlesM.push_back(hit->GetTrackP4().m()); - ActsParticlesQ.push_back(hit->GetCharge()); - - ActsParticlesEta.push_back(hit->GetTrackP4().eta()); - ActsParticlesPhi.push_back(hit->GetTrackP4().phi()); - ActsParticlesPt.push_back(pow(pow(hit->GetTrackP4().px(), 2) + pow(hit->GetTrackP4().py(), 2), 0.5)); - ActsParticlesP.push_back(pow(pow(hit->GetTrackP4().px(), 2) + pow(hit->GetTrackP4().py(), 2) + pow(hit->GetTrackP4().pz(), 2), 0.5)); - ActsParticlesVertexPrimary.push_back(hit->GetIsPrimaryTrack()); //? These variables need to be filled, but are unused by Acts - ActsParticlesVertexSecondary.push_back(hit->GetIsSecondaryTrack()); //? These variables need to be filled, but are unused by Acts - ActsParticlesParticle.push_back(1); //? These variables need to be filled, but are unused by Acts - ActsParticlesGeneration.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesSubParticle.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesELoss.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesPathInX0.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesPathInL0.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesNumberOfHits.push_back(0); //? These variables need to be filled, but are unused by Acts - ActsParticlesOutcome.push_back(0); //? These variables need to be filled, but are unused by Acts } // end of loop over hits - fActsParticlesTree->Fill(); } - G4cout << "Total FASER2 recorded hits: " << nHits << G4endl; } diff --git a/src/AnalysisManagerMessenger.cc b/src/AnalysisManagerMessenger.cc index 428eed5..db67641 100644 --- a/src/AnalysisManagerMessenger.cc +++ b/src/AnalysisManagerMessenger.cc @@ -36,6 +36,7 @@ #include "G4UIcmdWithAString.hh" #include "G4UIcmdWithAnInteger.hh" #include "G4UIcmdWithABool.hh" +#include "G4UIcmdWithADoubleAndUnit.hh" //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... @@ -50,10 +51,16 @@ fFileCmd->AvailableForStates(G4State_PreInit,G4State_Idle); fSaveTrackCmd = new G4UIcmdWithABool("/out/saveTrack", this); - fSaveTrackCmd->SetGuidance("whether save the information of all tracks"); + fSaveTrackCmd->SetGuidance("whether save the full trajectories"); fSaveTrackCmd->SetParameterName("saveTrack", true); fSaveTrackCmd->SetDefaultValue(false); + fParticleKinECutCmd = new G4UIcmdWithADoubleAndUnit("/out/parKinECut", this); + fParticleKinECutCmd->SetGuidance("set kinetic energy cut for particles tree"); + fParticleKinECutCmd->SetUnitCandidates("MeV GeV keV"); + fParticleKinECutCmd->SetUnitCategory("Energy"); + fParticleKinECutCmd->SetDefaultUnit("keV"); + fFLArEDir = new G4UIdirectory("/out/flare/"); fFLArEDir->SetGuidance("flare output control"); @@ -97,6 +104,7 @@ AnalysisManagerMessenger::~AnalysisManagerMessenger() { delete fFileCmd; delete fSaveTrackCmd; + delete fParticleKinECutCmd; delete fSave3DEvdCmd; delete fSave2DEvdCmd; delete fAddDiffusionCmd; @@ -114,6 +122,7 @@ void AnalysisManagerMessenger::SetNewValue(G4UIcommand* command,G4String newValu { if (command == fFileCmd) fAnalysisManager->setFileName(newValues); if (command == fSaveTrackCmd) fAnalysisManager->saveTrack(fSaveTrackCmd->GetNewBoolValue(newValues)); + if (command == fParticleKinECutCmd) fAnalysisManager->parKinECut(fParticleKinECutCmd->ConvertToDimensionedDouble(newValues)); if (command == fEnableFLArEOutCmd) fAnalysisManager->enableFLArE(fEnableFLArEOutCmd->GetNewBoolValue(newValues)); if (command == fSave3DEvdCmd) fAnalysisManager->save3DEvd(fSave3DEvdCmd->GetNewBoolValue(newValues)); diff --git a/src/EventInformation.cc b/src/EventInformation.cc index 5487745..7c2341a 100644 --- a/src/EventInformation.cc +++ b/src/EventInformation.cc @@ -13,13 +13,9 @@ EventInformation::EventInformation(std::vector genMetad EventInformation::~EventInformation() {} -void EventInformation::Print() const +void EventInformation::Print(int i) const { G4cout << G4endl; - G4cout << "EventInformation: " << fGenMetadata.size() << " vertex(es)" << G4endl; - for(int i=0; i aTrajAllocator; FPFTrajectory::FPFTrajectory() : fPositionRecord(0), fTrackID(0), fParentID(0), fPDGEncoding(0), fPDGCharge(0.0), fParticleName(""), - fProcessName(""), fInitialMomentum(G4ThreeVector()) + fProcessName(""), fInitialP4(G4LorentzVector()) {} FPFTrajectory::FPFTrajectory(const G4Track *aTrack, G4bool storePoints) { fTrackID = aTrack->GetTrackID(); fParentID = aTrack->GetParentID(); - fInitialMomentum = aTrack->GetMomentum(); + fInitialP4 = aTrack->GetDynamicParticle()->Get4Momentum(); G4ParticleDefinition *fpParticleDefinition = aTrack->GetDefinition(); fPDGEncoding = fpParticleDefinition->GetPDGEncoding(); @@ -49,7 +50,7 @@ FPFTrajectory::FPFTrajectory(FPFTrajectory &right) : G4VTrajectory() fPDGEncoding = right.fPDGEncoding; fTrackID = right.fTrackID; fParentID = right.fParentID; - fInitialMomentum = right.fInitialMomentum; + fInitialP4 = right.fInitialP4; fPositionRecord = new TrajectoryPointContainer(); fStorePoints = right.fStorePoints; fProcessName = right.fProcessName; @@ -104,7 +105,7 @@ const std::map *FPFTrajectory::GetAttDefs() const (*store)[PDG] = G4AttDef(PDG, "PDG Encoding", "Physics", "", "G4int"); G4String IMom("IMom"); - (*store)[IMom] = G4AttDef(IMom, "Momentum of track at start of trajectory", "Physics", "", "G4ThreeVector"); + (*store)[IMom] = G4AttDef(IMom, "4Momentum of track at start of trajectory", "Physics", "", "G4LorentzVector"); G4String NTP("NTP"); (*store)[NTP] = G4AttDef(NTP, "No. of points", "Physics", "", "G4int"); @@ -137,7 +138,7 @@ std::vector *FPFTrajectory::CreateAttValues() const values->push_back(G4AttValue("PDG", c.c_str(), "")); s.seekp(std::ios::beg); - s << G4BestUnit(fInitialMomentum, "Energy") << std::ends; + s << G4BestUnit(fInitialP4, "Energy") << std::ends; values->push_back(G4AttValue("IMom", c.c_str(), "")); s.seekp(std::ios::beg); @@ -164,7 +165,7 @@ G4ParticleDefinition *FPFTrajectory::GetParticleDefinition() const void FPFTrajectory::MergeTrajectory(G4VTrajectory *secondTrajectory) { - if (!secondTrajectory) + if (!secondTrajectory) return; FPFTrajectory *second = (FPFTrajectory *)secondTrajectory; G4int ent = second->GetPointEntries(); @@ -177,7 +178,13 @@ void FPFTrajectory::MergeTrajectory(G4VTrajectory *secondTrajectory) second->fPositionRecord->clear(); } -G4double FPFTrajectory::GetInitialKineticEnergy() const { +G4ThreeVector FPFTrajectory::GetInitialMomentum() const +{ + return G4ThreeVector(fInitialP4.px(), fInitialP4.py(), fInitialP4.pz()); +} + +G4double FPFTrajectory::GetInitialKineticEnergy() const +{ const G4ParticleDefinition* p = GetParticleDefinition(); double mom = GetInitialMomentum().mag(); if (!p) return mom; From f4e86eb0e548a8abadf6e4a83f726ffc636ecb59 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Thu, 18 Dec 2025 00:27:33 -0500 Subject: [PATCH 4/8] add ancestor id, fix map access --- include/AnalysisManager.hh | 1 + src/AnalysisManager.cc | 10 ++++++---- 2 files changed, 7 insertions(+), 4 deletions(-) diff --git a/include/AnalysisManager.hh b/include/AnalysisManager.hh index c410939..d554e3e 100644 --- a/include/AnalysisManager.hh +++ b/include/AnalysisManager.hh @@ -159,6 +159,7 @@ class AnalysisManager { int particle_TID; int particle_PID; int particle_PDG; + int particle_ancestor; std::string particle_process; float particle_vx; float particle_vy; diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index 8736b8c..7991725 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -135,6 +135,7 @@ void AnalysisManager::bookParTree() fPar->Branch("particle_pdg", &particle_PDG, "particle_pdg/I"); fPar->Branch("track_id", &particle_TID, "track_id/I"); fPar->Branch("parent_id", &particle_PID, "parent_id/I"); + fPar->Branch("ancestor_id", &particle_ancestor, "ancestor_id/I"); fPar->Branch("process", &particle_process); fPar->Branch("vx", &particle_vx, "vx/F"); fPar->Branch("vy", &particle_vy, "vy/F"); @@ -586,14 +587,11 @@ void AnalysisManager::FillParticlesTree(const G4Event *event) for (size_t i = 0; i < trajectoryContainer->entries(); ++i) { auto trajectory = static_cast((*trajectoryContainer)[i]); - - // do not save if below cut threshold - particle_ke = trajectory->GetInitialKineticEnergy(); - if(particle_ke < fParKinECut) continue; particle_TID = trajectory->GetTrackID(); particle_PID = trajectory->GetParentID(); particle_PDG = trajectory->GetPDGEncoding(); + particle_ancestor = trackToPrimaryAncestor.at(particle_TID); particle_process = trajectory->GetProcessName(); auto particleId = ActsFatras::Barcode(); @@ -623,6 +621,10 @@ void AnalysisManager::FillParticlesTree(const G4Event *event) particle_phi = trajectory->GetInitialP4().phi(); particle_pt = trajectory->GetInitialP4().vect().perp(); particle_p = trajectory->GetInitialP4().vect().mag(); + + // do not save if below cut threshold + particle_ke = trajectory->GetInitialKineticEnergy(); + if(particle_ke < fParKinECut) continue; traj_Npoints = trajectory->GetPointEntries(); traj_pointX.clear(); From 9e5140353105c8ceff10fe96ff97c02da345515f Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Thu, 18 Dec 2025 09:31:48 -0500 Subject: [PATCH 5/8] add ancestor fo flare tree --- include/AnalysisManager.hh | 1 + src/AnalysisManager.cc | 7 +++++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/include/AnalysisManager.hh b/include/AnalysisManager.hh index d554e3e..02a86b7 100644 --- a/include/AnalysisManager.hh +++ b/include/AnalysisManager.hh @@ -189,6 +189,7 @@ class AnalysisManager { UInt_t flareTrackID; UInt_t flareParticleID; UInt_t flareParentID; + UInt_t flareAncestorID; UInt_t flarePDG; UInt_t flareCopyNum; float_t flareT; diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index 7991725..b122760 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -175,6 +175,7 @@ void AnalysisManager::bookFLArETrees() fFLArEHits->Branch("flareTrackID", &flareTrackID, "flareTrackID/I"); fFLArEHits->Branch("flareBarcode", &flareParticleID, "flareParticleID/I"); fFLArEHits->Branch("flareParentID", &flareParentID, "flareParentID/I"); + fFLArEHits->Branch("flareAncestorID", &flareAncestorID, "flareAncestorID/I"); fFLArEHits->Branch("flarePDG", &flarePDG, "flarePDG/I"); fFLArEHits->Branch("flareCopyNum", &flareCopyNum, "flareCopyNum/I"); fFLArEHits->Branch("flareT", &flareT, "flareT/I"); @@ -196,6 +197,7 @@ void AnalysisManager::bookFLArETrees() fFLArEHCALHits->Branch("hadTrackID", &flareTrackID, "hadTrackID/I"); fFLArEHCALHits->Branch("hadBarcode", &flareParticleID, "hadParticleID/I"); fFLArEHCALHits->Branch("hadParentID", &flareParentID, "hadParentID/I"); + fFLArEHCALHits->Branch("hadAncestorID", &flareAncestorID, "hadAncestorID/I"); fFLArEHCALHits->Branch("hadPDG", &flarePDG, "hadPDG/I"); fFLArEHCALHits->Branch("hadCopyNum", &flareCopyNum, "hadCopyNum/I"); fFLArEHCALHits->Branch("hadT", &flareT, "hadT/I"); @@ -591,7 +593,7 @@ void AnalysisManager::FillParticlesTree(const G4Event *event) particle_TID = trajectory->GetTrackID(); particle_PID = trajectory->GetParentID(); particle_PDG = trajectory->GetPDGEncoding(); - particle_ancestor = trackToPrimaryAncestor.at(particle_TID); + particle_ancestor = GetTrackPrimaryAncestor(particle_TID); particle_process = trajectory->GetProcessName(); auto particleId = ActsFatras::Barcode(); @@ -692,6 +694,7 @@ void AnalysisManager::FillFLArEOutput() { flareTrackID = hit->GetTID(); flareParentID = hit->GetPID(); + flareAncestorID = GetTrackPrimaryAncestor(flareTrackID); flarePDG = hit->GetPDG(); flareCopyNum = hit->GetCopyNum(); flareT = hit->GetTime(); @@ -711,7 +714,7 @@ void AnalysisManager::FillFLArEOutput() double vtx_xyz[3] = {primaries[0].Vx(), primaries[0].Vy(), primaries[0].Vz()}; // which primary ancestor does this hit belong to? - G4int whichPrim = GetTrackPrimaryAncestor(flareTrackID); + G4int whichPrim = flareAncestorID; G4int whichIndex = whichPrim - 1; // need to start from zero // pseudo reco: track/shower length and width From d82bf297ffac9fec441acb8f53c6efa607d01bff Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Thu, 18 Dec 2025 14:13:12 -0500 Subject: [PATCH 6/8] change default cut --- README.md | 2 +- macros/backgrounds.mac | 4 ++-- src/AnalysisManager.cc | 2 +- src/AnalysisManagerMessenger.cc | 2 +- 4 files changed, 5 insertions(+), 5 deletions(-) diff --git a/README.md b/README.md index af95862..8aade8d 100644 --- a/README.md +++ b/README.md @@ -187,7 +187,7 @@ Older versions of FORESEE output events in the HepMC2 format. To run over HepMC2 |:--|:--| |/out/fileName | option for AnalysisManagerMessenger, set name of the file saving all analysis variables| |/out/saveTrack | if `true` save full trajectories, `false` by default | -|/out/parKinECut | set particle kinetic energy threshold for saving in output tree | +|/out/parKinECut | set particle kinetic energy threshold for saving in output tree (default: `100 keV`) | |/out/flare/enableOutput | if `true` save FLArE output, `true` by default | |/out/flare/save3DEvd | if `true` save 3D spatial distribution of energy deposition, `false` by default | |/out/flare/save2DEvd | if `true` save 2D spatial distribution of energy deposition, `false` by default | diff --git a/macros/backgrounds.mac b/macros/backgrounds.mac index 219576d..e7181d4 100644 --- a/macros/backgrounds.mac +++ b/macros/backgrounds.mac @@ -16,9 +16,9 @@ /out/flare/addDiffusion false /out/fileName test_backgrounds.root -# store and save trajectories for EVD -/tracking/storeTrajectory 0 +# store full trajectories for EVD /out/saveTrack false +/out/parKinECut 350 keV # shoot background equivalent to 1 "spill" /run/beamOn 1 diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index b122760..fd838f2 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -73,7 +73,7 @@ AnalysisManager::AnalysisManager() fActsHitsTree = nullptr; fSaveTrack = false; - fParKinECut = 50.0*keV; + fParKinECut = 100.0*keV; fSave3DEvd = false; fSave2DEvd = false; fSavePseudoReco = false; diff --git a/src/AnalysisManagerMessenger.cc b/src/AnalysisManagerMessenger.cc index db67641..2966384 100644 --- a/src/AnalysisManagerMessenger.cc +++ b/src/AnalysisManagerMessenger.cc @@ -56,7 +56,7 @@ fSaveTrackCmd->SetDefaultValue(false); fParticleKinECutCmd = new G4UIcmdWithADoubleAndUnit("/out/parKinECut", this); - fParticleKinECutCmd->SetGuidance("set kinetic energy cut for particles tree"); + fParticleKinECutCmd->SetGuidance("set kinetic energy cut for particle tree"); fParticleKinECutCmd->SetUnitCandidates("MeV GeV keV"); fParticleKinECutCmd->SetUnitCategory("Energy"); fParticleKinECutCmd->SetDefaultUnit("keV"); From c6da497a58d3da6966cbb37abeb326bce6cfe1ea Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 22 Dec 2025 13:20:44 -0500 Subject: [PATCH 7/8] save only particles that left hits --- README.md | 4 +- include/AnalysisManager.hh | 39 +++++-- include/AnalysisManagerMessenger.hh | 4 +- macros/backgrounds.mac | 4 +- src/AnalysisManager.cc | 165 +++++++++++++++++++++++----- src/AnalysisManagerMessenger.cc | 25 ++--- src/StackingAction.cc | 13 ++- src/TrackingAction.cc | 2 +- 8 files changed, 193 insertions(+), 63 deletions(-) diff --git a/README.md b/README.md index 8aade8d..0cddbc1 100644 --- a/README.md +++ b/README.md @@ -186,8 +186,8 @@ Older versions of FORESEE output events in the HepMC2 format. To run over HepMC2 |Command |Description | |:--|:--| |/out/fileName | option for AnalysisManagerMessenger, set name of the file saving all analysis variables| -|/out/saveTrack | if `true` save full trajectories, `false` by default | -|/out/parKinECut | set particle kinetic energy threshold for saving in output tree (default: `100 keV`) | +|/out/saveAllParticles | if `true` save all particles in the event, `false` by default | +|/out/saveTrajectories | if `true` save full trajectories, `false` by default | |/out/flare/enableOutput | if `true` save FLArE output, `true` by default | |/out/flare/save3DEvd | if `true` save 3D spatial distribution of energy deposition, `false` by default | |/out/flare/save2DEvd | if `true` save 2D spatial distribution of energy deposition, `false` by default | diff --git a/include/AnalysisManager.hh b/include/AnalysisManager.hh index 02a86b7..02afb67 100644 --- a/include/AnalysisManager.hh +++ b/include/AnalysisManager.hh @@ -4,6 +4,7 @@ #include #include #include +#include #include "globals.hh" #include "G4Event.hh" @@ -34,8 +35,8 @@ class AnalysisManager { //------------------------------------------------ // functions for controlling from the configuration file void setFileName(std::string val) { fFilename = val; } - void saveTrack(G4bool val) { fSaveTrack = val; } - void parKinECut(G4double val) { fParKinECut = val; } + void saveAllParticles(G4bool val) { fSaveAllParticles = val; } + void saveTrajectories(G4bool val) { fSaveTrajectories = val; } void saveActs(G4bool val) { fSaveActs = val; } void savePseudoReco(G4bool val) { fSavePseudoReco = val; } void addDiffusion(G4String val) { fAddDiffusion = val; } @@ -48,8 +49,17 @@ class AnalysisManager { void SetTrackPrimaryAncestor(G4int trackID, G4int ancestorID) { trackToPrimaryAncestor[trackID] = ancestorID; } G4int GetTrackPrimaryAncestor(G4int trackID) { return trackToPrimaryAncestor.at(trackID); } + // build TID to parentID association + // filled progressively from StackingAction + void SetTrackParentID(G4int trackID, G4int parentID) { trackIDtoParentID[trackID] = parentID; } + G4int GetTrackParentID(G4int trackID) { return trackIDtoParentID.at(trackID); } + // return whether saving full tracks in trajectories - G4bool GetSaveTrack() { return fSaveTrack; } + G4bool GetSaveTrajectories() { return fSaveTrajectories; } + + // register track and its ancestors for saving + void RegisterTrackAndAncestors(const G4int trackID); + std::uint64_t GetOrBuildParticleID(const G4int trackID); private: @@ -74,8 +84,8 @@ class AnalysisManager { static AnalysisManager* fInstance; AnalysisManagerMessenger* fMessenger; - G4bool fSaveTrack; - G4double fParKinECut; + G4bool fSaveAllParticles; + G4bool fSaveTrajectories; G4bool fSave3DEvd; G4bool fSave2DEvd; G4bool fSavePseudoReco; @@ -93,6 +103,19 @@ class AnalysisManager { std::vector primaries; std::vector primaryIDs; + // track to primary ancestor (track id to track id) + std::map trackToPrimaryAncestor; + + // map track id to its barcode + std::map trackIDtoParticleID; + std::map, unsigned int> nextSubIndex; + + // map track id to its parent id + std::map trackIDtoParentID; + + // list of track ids to be saved + std::unordered_set trackIDsToKeep; + //------------------------------------------------ // output files and trees std::string fFilename; @@ -111,12 +134,6 @@ class AnalysisManager { TDirectory* fFASER2Dir; TTree* fActsHitsTree; - // track to primary ancestor (track id to track id) - std::map trackToPrimaryAncestor; - - // map trackID to its barcode - std::map trackIDtoParticleID; - //--------------------------------------------------- // OUTPUT VARIABLES FOR COMMON TREES //--------------------------------------------------- diff --git a/include/AnalysisManagerMessenger.hh b/include/AnalysisManagerMessenger.hh index 47bc673..c85b5c2 100644 --- a/include/AnalysisManagerMessenger.hh +++ b/include/AnalysisManagerMessenger.hh @@ -58,8 +58,8 @@ class AnalysisManagerMessenger: public G4UImessenger G4UIdirectory* fOutDir; G4UIcmdWithAString* fFileCmd; - G4UIcmdWithABool* fSaveTrackCmd; - G4UIcmdWithADoubleAndUnit* fParticleKinECutCmd; + G4UIcmdWithABool* fSaveParticlesCmd; + G4UIcmdWithABool* fSaveTrajectoriesCmd; G4UIdirectory* fFLArEDir; G4UIcmdWithABool* fEnableFLArEOutCmd; diff --git a/macros/backgrounds.mac b/macros/backgrounds.mac index e7181d4..f07cbf2 100644 --- a/macros/backgrounds.mac +++ b/macros/backgrounds.mac @@ -17,8 +17,8 @@ /out/fileName test_backgrounds.root # store full trajectories for EVD -/out/saveTrack false -/out/parKinECut 350 keV +/out/saveAllParticles false +/out/saveTrajectories false # shoot background equivalent to 1 "spill" /run/beamOn 1 diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index fd838f2..5e1df88 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -6,6 +6,7 @@ #include #include #include +#include #include #include @@ -72,8 +73,8 @@ AnalysisManager::AnalysisManager() fFLArEPseudoReco = nullptr; fActsHitsTree = nullptr; - fSaveTrack = false; - fParKinECut = 100.0*keV; + fSaveAllParticles = false; + fSaveTrajectories = false; fSave3DEvd = false; fSave2DEvd = false; fSavePseudoReco = false; @@ -151,7 +152,7 @@ void AnalysisManager::bookParTree() fPar->Branch("pt", &particle_pt, "pt/F"); fPar->Branch("p", &particle_p, "p/F"); fPar->Branch("ke", &particle_ke, "ke/F"); - if (fSaveTrack) // if saving full trajectory, add vector of trajectory points + if (fSaveTrajectories) // if saving full trajectory, add vector of trajectory points { fPar->Branch("traj_Npoints", &traj_Npoints, "traj_Npoints/I"); fPar->Branch("traj_pointX", &traj_pointX); @@ -192,7 +193,6 @@ void AnalysisManager::bookFLArETrees() // FLArE HCAL hits fFLArEHCALHits = new TTree("hcal_hits", "FLArE HCAL hits info"); - fFLArEHCALHits->Branch("hadEvtID", &evtID, "hadEvtID/I"); fFLArEHCALHits->Branch("hadTrackID", &flareTrackID, "hadTrackID/I"); fFLArEHCALHits->Branch("hadBarcode", &flareParticleID, "hadParticleID/I"); @@ -311,8 +311,14 @@ void AnalysisManager::BeginOfRun() { G4cout << sdname.first << " " << sdname.second << G4endl; + // make it lowercase for case-insensitive check + std::string SDname = sdname.second; + std::transform(SDname.begin(), SDname.end(), SDname.begin(), + [](unsigned char c) { return std::tolower(c); }); + // if FLArE is enabled, book trees - if (fEnableFLArE && sdname.second.find("FLArE") != std::string::npos) + // but give option to disable output if requires + if (fEnableFLArE && SDname.find("flare") != std::string::npos) { fFlareSDs.push_back(sdname.first); bookFLArETrees(); @@ -320,7 +326,7 @@ void AnalysisManager::BeginOfRun() // if FASER2 is enabled, book ACTS trees // but give the option to disable output if required - else if (fSaveActs && sdname.second.find("FASER2") != std::string::npos) + else if (fSaveActs && SDname.find("faser") != std::string::npos) { fFaser2SDs.push_back(sdname.first); bookFASER2Trees(); @@ -395,8 +401,15 @@ void AnalysisManager::BeginOfEvent() // track ID to primary ancestor association trackToPrimaryAncestor.clear(); + // track ID to parent ID association (daughter-father) + trackIDtoParentID.clear(); + // track ID to particle ID association trackIDtoParticleID.clear(); + nextSubIndex.clear(); + + // track IDs to be stored in output tree + trackIDsToKeep.clear(); // trajectory points (if enabled) // (actually need reset at every track) @@ -449,10 +462,6 @@ void AnalysisManager::EndOfEvent(const G4Event *event) // FILL EVENT TREE FillEventTree(event); - //----------------------------------------------------------- - // FILL PARTICLES/TRAJECTORIES TREE - FillParticlesTree(event); - //----------------------------------------------------------- // FILL DETECTOR HITS @@ -470,6 +479,14 @@ void AnalysisManager::EndOfEvent(const G4Event *event) if (fFaser2SDs.size() > 0) FillFASER2Output(); + + //----------------------------------------------------------- + // FILL PARTICLES/TRAJECTORIES TREE + + // by default, only particles that left hits in SDs + their ancestors are saved + // so this needs the be filled last so we know what to save + FillParticlesTree(event); + } //--------------------------------------------------------------------- @@ -557,6 +574,7 @@ void AnalysisManager::FillEventTree(const G4Event *event) // store a copy as a FPFParticle for further processing // this accumulates across entire event (no vertex reset) primaryIDs.push_back(tid); // store to avoid duplicates + RegisterTrackAndAncestors(tid); // store for saving in particles tree primaries.push_back(FPFParticle(pdg, 0, tid, primaryIDs.size() - 1, 1, m, vtxX, vtxY, vtxZ, vtxT, @@ -574,6 +592,98 @@ void AnalysisManager::FillEventTree(const G4Event *event) //--------------------------------------------------------------------- //--------------------------------------------------------------------- +std::uint64_t AnalysisManager::GetOrBuildParticleID(const G4int trackID) +{ + // If we already have it, just return it + auto it = trackIDtoParticleID.find(trackID); + if (it != trackIDtoParticleID.end()) + return it->second; + + // Find the primary ancestor (your existing helper) + G4int primaryTID = GetTrackPrimaryAncestor(trackID); + + // Primary index (0-based) = Acts "particle" index + int primaryIndex = primaryTID - 1; + + // compute generation = distance from this track back to its primary + unsigned int generation = 0; + { + G4int tid = trackID; + // Walk up parents until we reach the primary + while (true) + { + auto itPar = trackIDtoParentID.find(tid); + if (itPar == trackIDtoParentID.end()) // should never happen + { + G4cout << "Something fishy going on in the 'trackIDtoParentID' map..." << G4endl; + break; + } + + G4int parentID = itPar->second; + if (parentID <= 0) // reached Geant4 primary + break; + + generation++; + if (parentID == primaryTID) // reached our primary ancestor + break; + + tid = parentID; + } + } + + // Assign sub-particle index + // According to the doc: primaries must have generation=0 & sub=0. + unsigned int subParticle = 0; + + int vertexPrimary = 1; // could map to true generator vertex idx? + int vertexSecondary = 0; // don't encode secondary vertices for now + + if (generation > 0) { + auto key = std::make_tuple(vertexPrimary, primaryIndex, (int)generation); + auto& counter = nextSubIndex[key]; // default-inits to 0 when first seen + subParticle = counter; + ++counter; + } + + // Build the ActsFatras::Barcode + ActsFatras::Barcode bc; + bc.setVertexPrimary(vertexPrimary); + bc.setVertexSecondary(vertexSecondary); + bc.setParticle(primaryIndex); // which primary + bc.setGeneration(generation); // depth from that primary + bc.setSubParticle(subParticle); // sibling index within generation + + std::uint64_t value = bc.value(); + + // Cache and return + trackIDtoParticleID[trackID] = value; + return value; +} + +//--------------------------------------------------------------------- + +void AnalysisManager::RegisterTrackAndAncestors(const G4int trackID) +{ + // is this track already registered? + // try to insert this track into the keep-set + auto [it, inserted] = trackIDsToKeep.insert(trackID); + if (!inserted) + { + // if already inserted, then its ancestors are as well + // just return + return; + } + + // now we go up checking its parents + G4int parentID = trackIDtoParentID.at(trackID); + if( parentID > 0 ) + { + RegisterTrackAndAncestors(parentID); + } +} + +//--------------------------------------------------------------------- + void AnalysisManager::FillParticlesTree(const G4Event *event) { int count_tracks = 0; @@ -581,7 +691,7 @@ void AnalysisManager::FillParticlesTree(const G4Event *event) auto trajectoryContainer = event->GetTrajectoryContainer(); if (!trajectoryContainer) { - G4cout << "No trajectories found in the event!" << G4endl; + G4cout << "No particle trajectories found in the event!" << G4endl; return; } @@ -595,19 +705,17 @@ void AnalysisManager::FillParticlesTree(const G4Event *event) particle_PDG = trajectory->GetPDGEncoding(); particle_ancestor = GetTrackPrimaryAncestor(particle_TID); particle_process = trajectory->GetProcessName(); + particle_id = GetOrBuildParticleID(particle_TID); - auto particleId = ActsFatras::Barcode(); - particleId.setVertexPrimary(1); - particleId.setVertexSecondary(0); - particleId.setParticle(particle_TID - 1); // The track ID is the primary particle index plus one - particleId.setGeneration(particle_PID); - sub_part_map.try_emplace(particle_TID- 1, sub_part_map.size()); - // This is a fudge - assumes that the secondary particles are always sub-particles of the primary particle - particleId.setSubParticle(particle_PID == 0 ? 0 : sub_part_map[particle_TID - 1]); - particle_id = particleId.value(); - - // store in map for use in other trees! - trackIDtoParticleID[particle_TID] = particle_id; + // NOTE: we don't save all particles (unless specifically requested) + // we save only what we registered for saving while scanning hits + if( !fSaveAllParticles ) + { + // does it exist in the registry? + if (trackIDsToKeep.find(particle_TID) == trackIDsToKeep.end()) { + continue; // if not, skip saving this particle entirely + } + } // first point is vertex (always stored even if traj is disabled) particle_vx = trajectory->GetPoint(0)->GetPosition().x(); @@ -623,10 +731,7 @@ void AnalysisManager::FillParticlesTree(const G4Event *event) particle_phi = trajectory->GetInitialP4().phi(); particle_pt = trajectory->GetInitialP4().vect().perp(); particle_p = trajectory->GetInitialP4().vect().mag(); - - // do not save if below cut threshold particle_ke = trajectory->GetInitialKineticEnergy(); - if(particle_ke < fParKinECut) continue; traj_Npoints = trajectory->GetPointEntries(); traj_pointX.clear(); @@ -693,6 +798,8 @@ void AnalysisManager::FillFLArEOutput() for (auto hit : *hitCollection->GetVector()) { flareTrackID = hit->GetTID(); + RegisterTrackAndAncestors(flareTrackID); // register for saving in particle tree + flareParentID = hit->GetPID(); flareAncestorID = GetTrackPrimaryAncestor(flareTrackID); flarePDG = hit->GetPDG(); @@ -708,7 +815,7 @@ void AnalysisManager::FillFLArEOutput() flareDeltaPy = hit->GetDeltaMom().y(); flareDeltaPz = hit->GetDeltaMom().z(); flareEdep = hit->GetEdep(); - flareParticleID = trackIDtoParticleID.at(flareTrackID); + flareParticleID = GetOrBuildParticleID(flareTrackID); double hit_position_xyz[3] = {flareX, flareY, flareZ}; double vtx_xyz[3] = {primaries[0].Vx(), primaries[0].Vy(), primaries[0].Vz()}; @@ -846,7 +953,9 @@ void AnalysisManager::FillFASER2Output() ActsHitsGeometryID = 0; int hitID = hit->GetTrackID(); - ActsHitsParticleID = trackIDtoParticleID.at(hitID); + RegisterTrackAndAncestors(hitID); // register for saving in particle tree + + ActsHitsParticleID = GetOrBuildParticleID(hitID); ActsHitsX = hit->GetX(); ActsHitsY = hit->GetY(); diff --git a/src/AnalysisManagerMessenger.cc b/src/AnalysisManagerMessenger.cc index 2966384..3bab6bc 100644 --- a/src/AnalysisManagerMessenger.cc +++ b/src/AnalysisManagerMessenger.cc @@ -50,16 +50,15 @@ fFileCmd->SetGuidance("set name for the histograms file"); fFileCmd->AvailableForStates(G4State_PreInit,G4State_Idle); - fSaveTrackCmd = new G4UIcmdWithABool("/out/saveTrack", this); - fSaveTrackCmd->SetGuidance("whether save the full trajectories"); - fSaveTrackCmd->SetParameterName("saveTrack", true); - fSaveTrackCmd->SetDefaultValue(false); + fSaveParticlesCmd = new G4UIcmdWithABool("/out/saveAllParticles", this); + fSaveParticlesCmd->SetGuidance("whether saving all particles in the event"); + fSaveParticlesCmd->SetParameterName("saveAllParticles", true); + fSaveParticlesCmd->SetDefaultValue(false); - fParticleKinECutCmd = new G4UIcmdWithADoubleAndUnit("/out/parKinECut", this); - fParticleKinECutCmd->SetGuidance("set kinetic energy cut for particle tree"); - fParticleKinECutCmd->SetUnitCandidates("MeV GeV keV"); - fParticleKinECutCmd->SetUnitCategory("Energy"); - fParticleKinECutCmd->SetDefaultUnit("keV"); + fSaveTrajectoriesCmd = new G4UIcmdWithABool("/out/saveTrajectories", this); + fSaveTrajectoriesCmd->SetGuidance("whether saving the full trajectories"); + fSaveTrajectoriesCmd->SetParameterName("saveTrajectories", true); + fSaveTrajectoriesCmd->SetDefaultValue(false); fFLArEDir = new G4UIdirectory("/out/flare/"); fFLArEDir->SetGuidance("flare output control"); @@ -103,8 +102,8 @@ AnalysisManagerMessenger::~AnalysisManagerMessenger() { delete fFileCmd; - delete fSaveTrackCmd; - delete fParticleKinECutCmd; + delete fSaveParticlesCmd; + delete fSaveTrajectoriesCmd; delete fSave3DEvdCmd; delete fSave2DEvdCmd; delete fAddDiffusionCmd; @@ -121,8 +120,8 @@ AnalysisManagerMessenger::~AnalysisManagerMessenger() void AnalysisManagerMessenger::SetNewValue(G4UIcommand* command,G4String newValues) { if (command == fFileCmd) fAnalysisManager->setFileName(newValues); - if (command == fSaveTrackCmd) fAnalysisManager->saveTrack(fSaveTrackCmd->GetNewBoolValue(newValues)); - if (command == fParticleKinECutCmd) fAnalysisManager->parKinECut(fParticleKinECutCmd->ConvertToDimensionedDouble(newValues)); + if (command == fSaveParticlesCmd) fAnalysisManager->saveAllParticles(fSaveParticlesCmd->GetNewBoolValue(newValues)); + if (command == fSaveTrajectoriesCmd) fAnalysisManager->saveTrajectories(fSaveTrajectoriesCmd->GetNewBoolValue(newValues)); if (command == fEnableFLArEOutCmd) fAnalysisManager->enableFLArE(fEnableFLArEOutCmd->GetNewBoolValue(newValues)); if (command == fSave3DEvdCmd) fAnalysisManager->save3DEvd(fSave3DEvdCmd->GetNewBoolValue(newValues)); diff --git a/src/StackingAction.cc b/src/StackingAction.cc index ba00128..5af353c 100644 --- a/src/StackingAction.cc +++ b/src/StackingAction.cc @@ -12,9 +12,12 @@ StackingAction::StackingAction(RunAction* aRunAction, EventAction* aEventAction) G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack (const G4Track* aTrack) { - // for each track, build track ID to primary ancestor - // primaries have themselves as ancestor - // go up the tree until original primary for everything else + // for each track: + // 1. build track ID to parent ID association + // 2. build track ID to primary ancestor association + // primaries have themselves as ancestor + // go up the tree until original primary for everything else + // 3. register it with EventAction (total primaries/secondaries in event) G4int trackID = aTrack->GetTrackID(); G4int parentID = aTrack->GetParentID(); G4int ancestorID = -1; @@ -39,8 +42,10 @@ G4ClassificationOfNewTrack StackingAction::ClassifyNewTrack (const G4Track* aTra ancestorID = AnalysisManager::GetInstance()->GetTrackPrimaryAncestor(parentID); } - // add track with its ancestor!!! + // add ancestor for this track in the map! AnalysisManager::GetInstance()->SetTrackPrimaryAncestor(trackID,ancestorID); + // add direct parent for this track in the map! + AnalysisManager::GetInstance()->SetTrackParentID(trackID, parentID); // Do not affect track classification. Just return what would have // been returned by the base class diff --git a/src/TrackingAction.cc b/src/TrackingAction.cc index 46e6db6..8b40d5d 100644 --- a/src/TrackingAction.cc +++ b/src/TrackingAction.cc @@ -11,7 +11,7 @@ TrackingAction::TrackingAction() : G4UserTrackingAction() {;} void TrackingAction::PreUserTrackingAction(const G4Track* aTrack) { // find out whether saving the full trajectory points or not - bool storeTrajectoryPoints = AnalysisManager::GetInstance()->GetSaveTrack(); + bool storeTrajectoryPoints = AnalysisManager::GetInstance()->GetSaveTrajectories(); // initialize our custom trajectory class auto *traj = new FPFTrajectory(aTrack, storeTrajectoryPoints); From e97b6fbc3064eee33665d0c18468d403ea086ccb Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 22 Dec 2025 13:21:00 -0500 Subject: [PATCH 8/8] fix branch type --- include/AnalysisManager.hh | 4 ++-- src/AnalysisManager.cc | 6 +++--- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/include/AnalysisManager.hh b/include/AnalysisManager.hh index 02afb67..9a163e6 100644 --- a/include/AnalysisManager.hh +++ b/include/AnalysisManager.hh @@ -172,7 +172,7 @@ class AnalysisManager { //--------------------------------------------------- // Output variables for PARTICLES/TRAJECTORIES tree - std::uint64_t particle_id; //barcode + ULong64_t particle_id; //barcode int particle_TID; int particle_PID; int particle_PDG; @@ -204,7 +204,7 @@ class AnalysisManager { PixelMap3D* pm3D; UInt_t flareTrackID; - UInt_t flareParticleID; + ULong64_t flareParticleID; UInt_t flareParentID; UInt_t flareAncestorID; UInt_t flarePDG; diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index 5e1df88..ee295d1 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -132,7 +132,7 @@ void AnalysisManager::bookParTree() { fPar = new TTree("particles", "particle info"); fPar->Branch("event_id", &evtID, "event_id/I"); - fPar->Branch("particle_id", &particle_id, "particle_id/I"); + fPar->Branch("particle_id", &particle_id, "particle_id/l"); fPar->Branch("particle_pdg", &particle_PDG, "particle_pdg/I"); fPar->Branch("track_id", &particle_TID, "track_id/I"); fPar->Branch("parent_id", &particle_PID, "parent_id/I"); @@ -174,7 +174,7 @@ void AnalysisManager::bookFLArETrees() fFLArEHits = new TTree("flare_hits", "FLArE hits info"); fFLArEHits->Branch("flareEvtID", &evtID, "flareEvtID/I"); fFLArEHits->Branch("flareTrackID", &flareTrackID, "flareTrackID/I"); - fFLArEHits->Branch("flareBarcode", &flareParticleID, "flareParticleID/I"); + fFLArEHits->Branch("flareBarcode", &flareParticleID, "flareParticleID/l"); fFLArEHits->Branch("flareParentID", &flareParentID, "flareParentID/I"); fFLArEHits->Branch("flareAncestorID", &flareAncestorID, "flareAncestorID/I"); fFLArEHits->Branch("flarePDG", &flarePDG, "flarePDG/I"); @@ -195,7 +195,7 @@ void AnalysisManager::bookFLArETrees() fFLArEHCALHits = new TTree("hcal_hits", "FLArE HCAL hits info"); fFLArEHCALHits->Branch("hadEvtID", &evtID, "hadEvtID/I"); fFLArEHCALHits->Branch("hadTrackID", &flareTrackID, "hadTrackID/I"); - fFLArEHCALHits->Branch("hadBarcode", &flareParticleID, "hadParticleID/I"); + fFLArEHCALHits->Branch("hadBarcode", &flareParticleID, "hadParticleID/l"); fFLArEHCALHits->Branch("hadParentID", &flareParentID, "hadParentID/I"); fFLArEHCALHits->Branch("hadAncestorID", &flareAncestorID, "hadAncestorID/I"); fFLArEHCALHits->Branch("hadPDG", &flarePDG, "hadPDG/I");