From e89c34ca133653321a71355710ea3a8df3221c54 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Fri, 22 Aug 2025 17:22:54 -0400 Subject: [PATCH 1/8] preliminary metadata struct --- include/generators/GeneratorBase.hh | 14 ++++++-- include/generators/GeneratorVertexMetadata.hh | 32 +++++++++++++++++++ src/PrimaryGeneratorAction.cc | 3 ++ 3 files changed, 47 insertions(+), 2 deletions(-) create mode 100644 include/generators/GeneratorVertexMetadata.hh diff --git a/include/generators/GeneratorBase.hh b/include/generators/GeneratorBase.hh index 48b9e85..f282b74 100644 --- a/include/generators/GeneratorBase.hh +++ b/include/generators/GeneratorBase.hh @@ -1,8 +1,10 @@ #ifndef GENERATOR_BASE_HH #define GENERATOR_BASE_HH +#include #include "G4Event.hh" #include "G4UImessenger.hh" +#include "generators/GeneratorVertexMetadata.hh" class GeneratorBase { @@ -19,13 +21,21 @@ class GeneratorBase virtual void GeneratePrimaries(G4Event *event) = 0; // return name of current generator - G4String GetGeneratorName() { return fGeneratorName; } + const G4String GetGeneratorName() { return fGeneratorName; } + + // reset the list of metadata + void ResetEventMetadata() { fVertexMetadata.clear(); } + + // return full event vertex metadata + const std::vector GetEventMetadata() { return fVertexMetadata; } + // return single vertex metadata + const GeneratorVertexMetadata GetEventMetadataPerVertex(G4int i) { return fVertexMetadata.at(i); } protected : G4String fGeneratorName; G4UImessenger* fMessenger; - + std::vector fVertexMetadata; }; #endif diff --git a/include/generators/GeneratorVertexMetadata.hh b/include/generators/GeneratorVertexMetadata.hh new file mode 100644 index 0000000..ae3bc80 --- /dev/null +++ b/include/generators/GeneratorVertexMetadata.hh @@ -0,0 +1,32 @@ +#ifndef GENERATOR_VERTEX_METDATA_HH +#define GENERATOR_VERTEX_METDATA_HH + +#include +#include "TLorentzVector.h" + +// Metadata information on each generated vertex +// Store input interaction/decay process information +struct GeneratorVertexMetadata +{ + std::string generatorType; ///< which generator class + std::string processName; ///< which process + double weight; ///< process weight + + int pdg; ///< initiator pdg + TLorentzVector x4; ///< initiator 4-position (vertex) + TLorentzVector p4; ///< initiator 4-momentum + double mass; ///< initiator mass + double charge; ///< initiator charge + + int tgt_pdg; ///< target pdg + int tgt_A; ///< nuclear target A + int tgt_Z; ///< nuclear target Z + + double xs = 0.0; ///< cross-section + double Q2 = -1.0; ///< momentum transfer + double xBj = 0.0; ///< Bjorken x + double y = -1.0; ///< inelasticity + double W = -1.0; ///< hadronic invariant mass +}; + +#endif \ No newline at end of file diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc index cd659b0..0a11d38 100644 --- a/src/PrimaryGeneratorAction.cc +++ b/src/PrimaryGeneratorAction.cc @@ -66,6 +66,9 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) G4cout << G4endl; G4cout << "===oooOOOooo=== Event Generator (# " << anEvent->GetEventID(); + // reset event metadata + fGenerator->ResetEventMetadata(); + // produce an event with current generator fGenerator->GeneratePrimaries(anEvent); From 8c8e8b24c63686f2645e14d7c5618dd4e6aef3cd Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 25 Aug 2025 17:00:40 -0400 Subject: [PATCH 2/8] propagate metadata to ouptut for genie --- include/AnalysisManager.hh | 28 ++++-- include/EventInformation.hh | 32 ++++++ include/PixelMap3D.hh | 4 +- include/PrimaryParticleInformation.hh | 56 +---------- include/generators/GENIEGenerator.hh | 25 +---- include/generators/GeneratorBase.hh | 6 +- include/generators/GeneratorVertexMetadata.hh | 20 ++-- src/AnalysisManager.cc | 97 +++++++++++++++---- src/EventInformation.cc | 32 ++++++ src/PixelMap3D.cc | 5 +- src/PrimaryGeneratorAction.cc | 33 +++---- src/PrimaryParticleInformation.cc | 25 +---- src/generators/GENIEGenerator.cc | 93 +++++++++++++----- 13 files changed, 274 insertions(+), 182 deletions(-) create mode 100644 include/EventInformation.hh create mode 100644 src/EventInformation.cc diff --git a/include/AnalysisManager.hh b/include/AnalysisManager.hh index cbe7829..5693994 100644 --- a/include/AnalysisManager.hh +++ b/include/AnalysisManager.hh @@ -62,6 +62,7 @@ class AnalysisManager { void bookFLArETrees(); void bookFASER2Trees(); + void FillEventTree(const G4Event* event); void FillPrimariesTree(const G4Event* event); void FillTrajectoriesTree(const G4Event* event); @@ -120,13 +121,29 @@ class AnalysisManager { //--------------------------------------------------- // OUTPUT VARIABLES FOR COMMON TREES - // TODO: review carefully - // TODO: need to make evt tree less FLARE-centric - // TODO: remove pseudo-reco or add it as FLARE-only tree - // TODO: turn arrays in std::vector! G4int evtID; - FPFNeutrino neutrino; //TODO: remove?? + G4int vertexID; + double weight; + std::string genType; + std::string processName; + int initPDG; + double initX, initY, initZ, initT; + double initPx, initPy, initPz, initE; + double initM; + double initQ; + int intType; + int scatteringType; + int fslPDG; + int tgtPDG; + int tgtA; + int tgtZ; + int hitnucPDG; + double xs; + double Q2; + double xBj; + double y; + double W; //--------------------------------------------------- // Output variables for TRAJECTORIES tree @@ -164,7 +181,6 @@ class AnalysisManager { //--------------------------------------------------- // OUTPUT VARIABLES FOR FLArE TREES // TODO: merge hit variables? no need to use different names? - // TODO: somehow need to add back here info from evt tree PixelMap3D* pm3D; diff --git a/include/EventInformation.hh b/include/EventInformation.hh new file mode 100644 index 0000000..acb08cc --- /dev/null +++ b/include/EventInformation.hh @@ -0,0 +1,32 @@ +#ifndef EventInformation_HH +#define EventInformation_HH + +#include "G4VUserEventInformation.hh" +#include "generators/GeneratorVertexMetadata.hh" +#include "globals.hh" +#include + +// Event information +class EventInformation : public G4VUserEventInformation +{ + public: + + EventInformation(); + EventInformation(std::vector genMetadata); + virtual ~EventInformation(); + + /// Gets vertex metadata full vector + inline std::vector GetEventMetadata() const { return fGenMetadata; } + + /// Gets metadata per vertex + inline GeneratorVertexMetadata GetMetadataPerVertex(int i) const { return fGenMetadata.at(i); } + + /// Prints the information about the event. + virtual void Print() const; + + private: + /// Set of vertex metadata + std::vector fGenMetadata; +}; + +#endif diff --git a/include/PixelMap3D.hh b/include/PixelMap3D.hh index 9c92578..0e5d849 100644 --- a/include/PixelMap3D.hh +++ b/include/PixelMap3D.hh @@ -10,8 +10,6 @@ #include "G4ThreeVector.hh" -#include "FPFNeutrino.hh" - #include "hep_hpc/hdf5/File.hpp" //#include "hep_hpc/hdf5/Ntuple.hpp" @@ -27,7 +25,7 @@ class PixelMap3D { void FillEntryWithToyElectronTransportation(const Double_t* pos_xyz, const Double_t* vtx_xyz, Double_t edep, const Int_t idxPrim); void FillEntryWithToySingleElectronTransportation(const Double_t* pos_xyz, const Double_t* vtx_xyz, Double_t edep, const Int_t idxPrim); void Write2DPMToFile(TFile* thefile, TDirectory *thedir); - void Process3DPM(hep_hpc::hdf5::File &h5file, FPFNeutrino neutrino, G4bool save3D); + void Process3DPM(hep_hpc::hdf5::File &h5file, G4int initPDG, G4int fslPDG, G4int intType, G4int scatType, G4double initE, G4bool save3D); // this should go to a Geometry Service class G4double DistanceToAnode(G4double x); diff --git a/include/PrimaryParticleInformation.hh b/include/PrimaryParticleInformation.hh index 767dc05..4bb5592 100644 --- a/include/PrimaryParticleInformation.hh +++ b/include/PrimaryParticleInformation.hh @@ -16,19 +16,8 @@ class PrimaryParticleInformation : public G4VUserPrimaryParticleInformation { /// @param aMass A mass of the particle. /// @param aMomentum An initial particle momentum (at the primary vertex). /// @param aVertex An initial particle vertex. - /// @param neutrino index of the neutrino interaction if it is from GENIE - /// @param neutrino PDG of the neutrino interaction if it is from GENIE - /// @param neutrino P4 of the neutrino interaction if it is from GENIE - /// @param neutrino X4 of the neutrino interaction input to the detector - /// @param interaction type of the neutrino interaction if it is from GENIE - /// @param scattering type of the neutrino interaction if it is from GENIE - /// @param fsl PDG of the neutrino interaction if it is from GENIE - /// @param fsl P4 of the neutrino interaction if it is from GENIE PrimaryParticleInformation(G4int aID, G4int aPDG, G4double aMass, G4double aCharge, - G4ThreeVector aMomentum, G4ThreeVector aVertex, - G4int aneuIdx, G4int aneuPDG, TLorentzVector aneuP4, TLorentzVector aneuX4, - G4int aInttype, G4int aScatteringtype, G4double aW, - G4int afslPDG, TLorentzVector afslP4); + G4ThreeVector aMomentum, G4ThreeVector aVertex); virtual ~PrimaryParticleInformation(); @@ -38,7 +27,6 @@ class PrimaryParticleInformation : public G4VUserPrimaryParticleInformation { /// Get the particle PDG code inline G4int GetPDG() const { return fPDG; }; - /// Get the particle charge inline G4double GetCharge() const { return fCharge; }; @@ -51,17 +39,6 @@ class PrimaryParticleInformation : public G4VUserPrimaryParticleInformation { /// Get the particle initial vertex. inline G4ThreeVector GetVertexMC() const { return fVertexMC; }; - /// Get the truth information of the neutrino interaction if it's from GENIE - inline G4int GetNeuIdx() const { return fNeuIdx; }; - inline G4int GetNeuPDG() const { return fNeuPDG; }; - inline TLorentzVector GetNeuP4() const { return fNeuP4; }; - inline TLorentzVector GetNeuX4() const { return fNeuX4; }; - inline G4int GetInteractionTypeId() const { return fInteractionType; }; - inline G4int GetScatteringTypeId() const { return fScatteringType; }; - inline G4double GetW() const { return fW; }; - inline G4int GetFSLPDG() const { return fFSLPDG; }; - inline TLorentzVector GetFSLP4() const { return fFSLP4; }; - /// Prints the information about the particle. virtual void Print() const; @@ -82,37 +59,6 @@ class PrimaryParticleInformation : public G4VUserPrimaryParticleInformation { /// A particle initial vertex (from particle generator) G4ThreeVector fVertexMC; - /// neutrino index of the neutrino interaction if it is from GENIE - G4int fNeuIdx; - - /// neutrino PDG of the neutrino interaction if it is from GENIE - G4int fNeuPDG; - - /// neutrino P4 of the neutrino interaction if it is from GENIE - TLorentzVector fNeuP4; - - /// neutrino X4 of the neutrino interaction input to the detector - TLorentzVector fNeuX4; - - /// Interaction type of the neutrino interaction if it is from GENIE - G4int fInteractionType; - - /// Scattering type of the neutrino interaction if it is from GENIE - G4int fScatteringType; - - /// invariant hadronic mass - G4double fW; - - /// fsl PDG of the neutrino interaction if it is from GENIE - G4int fFSLPDG; - - /// fsl P4 of the neutrino interaction if it is from GENIE - TLorentzVector fFSLP4; - - /// fsl X4 of the neutrino interaction if it is from GENIE - TLorentzVector fFSLX4; - - /// particle's charge G4double fCharge; diff --git a/include/generators/GENIEGenerator.hh b/include/generators/GENIEGenerator.hh index 81a063b..bfbe0a8 100644 --- a/include/generators/GENIEGenerator.hh +++ b/include/generators/GENIEGenerator.hh @@ -6,7 +6,6 @@ #include "TFile.h" #include "TTree.h" -#include "TLorentzVector.h" #include "globals.hh" class G4Event; @@ -26,17 +25,6 @@ class GENIEGenerator : public GeneratorBase void SetEvtStartIdx(G4int val) { fEvtStartIdx = val; } void SetRandomVertex(G4bool val) { fRandomVtx = val; } - // get methods - G4int GetNeuIdx() { return m_neuIdx; }; - G4int GetNeuPDG() { return m_neuPDG; }; - TLorentzVector GetNeuP4() { return m_neuP4; }; - TLorentzVector GetNeuX4() { return m_neuX4; }; - G4int GetInteractionTypeId() { return m_int_type; }; - G4int GetScatteringTypeId() { return m_scattering_type; }; - G4int GetFSLPDG() { return m_fslPDG; }; - TLorentzVector GetFSLP4() { return m_fslP4; }; - G4double GetW() { return m_W; }; - private: G4String fGSTFilename; G4int fEventCounter; @@ -49,6 +37,7 @@ class GENIEGenerator : public GeneratorBase // define the branches we are interested in from the GST tree const G4int fkNPmax = 250; G4int fNEntries; + G4double m_wght; G4double m_Ev, m_pxv, m_pyv, m_pzv; G4double m_El, m_pxl, m_pyl, m_pzl; G4int m_nf; @@ -58,20 +47,14 @@ class GENIEGenerator : public GeneratorBase G4bool m_imdanh, m_nuel, m_em, m_cc, m_nc, m_charm; G4bool m_singlek, m_amnugamma; G4int m_neuPDG, m_fslPDG; - - // for output - G4int m_neuIdx; - TLorentzVector m_neuP4; - TLorentzVector m_neuX4; - G4int m_int_type; - G4int m_scattering_type; - TLorentzVector m_fslP4; - G4double m_W; + G4int m_tgt, m_Z, m_A, m_hitnuc; + G4double m_x, m_y, m_Q2, m_W; // specific internal functions G4bool FindParticleDefinition(G4int pdg, G4ParticleDefinition* &particleDefinition) const; G4int DecodeInteractionType() const; G4int DecodeScatteringType() const; + G4String EncodeProcessName() const; }; #endif diff --git a/include/generators/GeneratorBase.hh b/include/generators/GeneratorBase.hh index f282b74..de93f6e 100644 --- a/include/generators/GeneratorBase.hh +++ b/include/generators/GeneratorBase.hh @@ -21,15 +21,15 @@ class GeneratorBase virtual void GeneratePrimaries(G4Event *event) = 0; // return name of current generator - const G4String GetGeneratorName() { return fGeneratorName; } + G4String GetGeneratorName() const { return fGeneratorName; } // reset the list of metadata void ResetEventMetadata() { fVertexMetadata.clear(); } // return full event vertex metadata - const std::vector GetEventMetadata() { return fVertexMetadata; } + std::vector GetEventMetadata() const { return fVertexMetadata; } // return single vertex metadata - const GeneratorVertexMetadata GetEventMetadataPerVertex(G4int i) { return fVertexMetadata.at(i); } + GeneratorVertexMetadata GetEventMetadataPerVertex(G4int i) const { return fVertexMetadata.at(i); } protected : diff --git a/include/generators/GeneratorVertexMetadata.hh b/include/generators/GeneratorVertexMetadata.hh index ae3bc80..985c18b 100644 --- a/include/generators/GeneratorVertexMetadata.hh +++ b/include/generators/GeneratorVertexMetadata.hh @@ -15,18 +15,24 @@ struct GeneratorVertexMetadata int pdg; ///< initiator pdg TLorentzVector x4; ///< initiator 4-position (vertex) TLorentzVector p4; ///< initiator 4-momentum - double mass; ///< initiator mass - double charge; ///< initiator charge + double mass; ///< initiator mass + double charge; ///< initiator charge - int tgt_pdg; ///< target pdg - int tgt_A; ///< nuclear target A - int tgt_Z; ///< nuclear target Z + int intType = -1; ///< interaction type + int scatteringType = -1; ///< scattering type + int fsl_pdg = -1; ///< final state lepton pdg + + int tgt_pdg = -1; ///< target pdg + int tgt_A = -1; ///< nuclear target A + int tgt_Z = -1; ///< nuclear target Z + int hitnuc_pdg = -1; ///< hit nucleon pdg - double xs = 0.0; ///< cross-section + double xs = -1.0; ///< cross-section double Q2 = -1.0; ///< momentum transfer - double xBj = 0.0; ///< Bjorken x + double xBj = 0.0; ///< Bjorken x double y = -1.0; ///< inelasticity double W = -1.0; ///< hadronic invariant mass + }; #endif \ No newline at end of file diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index 7c1e96b..3a472f2 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -27,13 +27,14 @@ #include "FASER2TrackerSD.hh" #include "LArBoxHit.hh" #include "PrimaryParticleInformation.hh" +#include "EventInformation.hh" +#include "generators/GeneratorVertexMetadata.hh" #include "reco/PCAAnalysis3D.hh" #include "reco/Cluster3D.hh" #include "reco/LinearFit.hh" #include "reco/ShowerLID.hh" #include "reco/Barcode.hh" #include "FPFParticle.hh" -#include "FPFNeutrino.hh" //--------------------------------------------------------------------- //--------------------------------------------------------------------- @@ -89,9 +90,33 @@ void AnalysisManager::bookEvtTree() { fEvt = new TTree("event", "event info"); fEvt->Branch("evtID", &evtID, "evtID/I"); - - //TODO EXPAND!!!! - + fEvt->Branch("vtxID", &vertexID, "vtxID/I"); + fEvt->Branch("weight", &weight, "weight/D"); + fEvt->Branch("genType", &genType); + fEvt->Branch("processName", &processName); + 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("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"); + 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("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() @@ -382,12 +407,9 @@ void AnalysisManager::EndOfRun() void AnalysisManager::BeginOfEvent() { - // TODO: REVIEW, CLEAN-UP UNUSED VARIABLES - // TODO: change arrays into std::vectors // reset vectors that need to be cleared for a new event - // only reset arrays or vectors, no need for other defaults?? + // only reset arrays or vectors, tipically no need for other defaults - neutrino = FPFNeutrino(); // remove??? primaries.clear(); primaryIDs.clear(); @@ -467,9 +489,7 @@ void AnalysisManager::EndOfEvent(const G4Event *event) evtID = event->GetEventID(); // FILL EVENT TREE - // nothing else except evtID for now... - // TODO: need to find good way to propagate generator-level info - fEvt->Fill(); + FillEventTree(event); //----------------------------------------------------------- @@ -501,6 +521,48 @@ void AnalysisManager::EndOfEvent(const G4Event *event) //--------------------------------------------------------------------- //--------------------------------------------------------------------- +void AnalysisManager::FillEventTree(const G4Event *event) +{ + EventInformation* eventInfo = static_cast(event->GetUserInformation()); + eventInfo->Print(); + auto metadata = eventInfo->GetEventMetadata(); + for(int i=0; iFill(); + } +} + +//--------------------------------------------------------------------- +//--------------------------------------------------------------------- + void AnalysisManager::FillPrimariesTree(const G4Event *event) { nPrimaryVertex = event->GetNumberOfPrimaryVertex(); @@ -614,12 +676,7 @@ void AnalysisManager::FillFLArEOutput() // prepare the pixel map const double_t res_tpc[3] = {1, 5, 5}; // mm - // TODO: rethink once generator info is included... - /*if (neutrino.NuPDG()!=0) { - pm3D = new PixelMap3D(evtID, nPrimaryParticle, neutrino.NuPDG(), res_tpc); - } else { */ - - pm3D = new PixelMap3D(evtID, primaries.size(), primaries[0].PDG(), res_tpc); + pm3D = new PixelMap3D(evtID, primaries.size(), initPDG, res_tpc); // boundaries in global coordinates pm3D->SetPMBoundary(GeometricalParameters::Get()->GetFLArEPosition()/mm - GeometricalParameters::Get()->GetTPCSizeXYZ()/mm/2, @@ -756,7 +813,7 @@ void AnalysisManager::FillFLArEOutput() if (fSave2DEvd) pm3D->Write2DPMToFile(fFile,fFLArEDir); - pm3D->Process3DPM(fH5file, neutrino, fSave3DEvd); + pm3D->Process3DPM(fH5file, initPDG, fslPDG, intType, scatteringType, initE, fSave3DEvd); sparseFractionMem = pm3D->GetSparseFractionMem(); sparseFractionBins = pm3D->GetSparseFractionBins(); @@ -958,7 +1015,7 @@ void AnalysisManager::FillFLArEPseudoReco() } slid::ShowerLID *shwlid = new slid::ShowerLID(pm3D->Get3DPixelMap(), - neutrino.NuVx(), neutrino.NuVy(), neutrino.NuVz(), 0., 0., 1.); + initX, initY, initZ, 0., 0., 1.); Double_t *ptr_dedx = shwlid->GetTotalDedxLongitudinal(); std::copy(ptr_dedx, ptr_dedx + 3000, TotalDedxLongitudinal.begin()); @@ -969,7 +1026,7 @@ void AnalysisManager::FillFLArEPseudoReco() pm3D->Get2DPixelMapZY(iPrim + 1), pm3D->Get2DVtxPixelMapZX(iPrim + 1), pm3D->Get2DVtxPixelMapZY(iPrim + 1), - neutrino.NuVx(), neutrino.NuVy(), neutrino.NuVz(), + initX, initY, initZ, primaries[iPrim].Vx(), primaries[iPrim].Vy(), primaries[iPrim].Vz()); dir_pol_x[iPrim] = linFit->GetDir().X(); dir_pol_y[iPrim] = linFit->GetDir().Y(); diff --git a/src/EventInformation.cc b/src/EventInformation.cc new file mode 100644 index 0000000..6aed599 --- /dev/null +++ b/src/EventInformation.cc @@ -0,0 +1,32 @@ +#include "EventInformation.hh" +#include "generators/GeneratorVertexMetadata.hh" +#include + +EventInformation::EventInformation() +{} + +EventInformation::EventInformation(std::vector genMetadata) +{ + fGenMetadata = genMetadata; +} + +EventInformation::~EventInformation() +{} + +void EventInformation::Print() const +{ + G4cout << G4endl; + G4cout << "EventInformation: " << fGenMetadata.size() << " vertex(es)" << G4endl; + for(int i=0; i("mode"), make_scalar_column("nuE")); - evt_data.insert(fEvtID, nu.NuPDG(), nu.NuFSLPDG(), - nu.NuIntType(), nu.NuScatteringType(), nu.NuE()); + evt_data.insert(fEvtID, initPDG, fslPDG, intType, scatType, initE); static auto pm_data = make_ntuple({h5file, "pm_data"}, diff --git a/src/PrimaryGeneratorAction.cc b/src/PrimaryGeneratorAction.cc index 0a11d38..83ad525 100644 --- a/src/PrimaryGeneratorAction.cc +++ b/src/PrimaryGeneratorAction.cc @@ -10,6 +10,7 @@ #include "geometry/GeometricalParameters.hh" #include "PrimaryParticleInformation.hh" +#include "EventInformation.hh" #include "G4Event.hh" #include "G4Exception.hh" @@ -72,23 +73,10 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) // produce an event with current generator fGenerator->GeneratePrimaries(anEvent); - // save additional truth information alongside primary particles - // this makes it available for the output tree (e.g: neutrino info for genie) - - // TODO: make it more generic for other generators? - bool isGenie = (fGenerator->GetGeneratorName() == "genie"); - // downcast: if it fails, we won't use it anyway... - GENIEGenerator *genieGen = dynamic_cast(fGenerator); - G4int neuidx = (isGenie) ? genieGen->GetNeuIdx() : -1; - G4int neupdg = (isGenie) ? genieGen->GetNeuPDG() : -1; - TLorentzVector neup4 = (isGenie) ? genieGen->GetNeuP4() : TLorentzVector(); - TLorentzVector neux4 = (isGenie) ? genieGen->GetNeuX4() : TLorentzVector(); - G4int int_type = (isGenie) ? genieGen->GetInteractionTypeId() : -1; - G4int scattering_type = (isGenie) ? genieGen->GetScatteringTypeId() : -1; - G4double w = (isGenie) ? genieGen->GetW() : -1; - G4int fslpdg = (isGenie) ? genieGen->GetFSLPDG() : -1; - TLorentzVector fslp4 = (isGenie) ? genieGen->GetFSLP4() : TLorentzVector(); + // save vertex metadata information into the event + anEvent->SetUserInformation(new EventInformation(fGenerator->GetEventMetadata())); + // save additional truth information alongside primary particles // loop over the vertices, and then over primary particles, // and for each primary particle create an info object G4int count_particles = 0; @@ -99,14 +87,17 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) if (primary_particle) { primary_particle->SetUserInformation(new PrimaryParticleInformation( - count_particles, primary_particle->GetPDGcode(), primary_particle->GetMass(), - primary_particle->GetCharge(), - primary_particle->GetMomentum(), anEvent->GetPrimaryVertex(ivtx)->GetPosition(), - neuidx, neupdg, neup4, neux4, int_type, scattering_type, w, - fslpdg, fslp4)); + count_particles, + primary_particle->GetPDGcode(), + primary_particle->GetMass(), + primary_particle->GetCharge(), + primary_particle->GetMomentum(), + anEvent->GetPrimaryVertex(ivtx)->GetPosition() + )); count_particles++; } } } + } diff --git a/src/PrimaryParticleInformation.cc b/src/PrimaryParticleInformation.cc index 2413429..cbf6441 100644 --- a/src/PrimaryParticleInformation.cc +++ b/src/PrimaryParticleInformation.cc @@ -1,21 +1,13 @@ #include "PrimaryParticleInformation.hh" PrimaryParticleInformation::PrimaryParticleInformation(G4int aID, G4int aPDG, G4double aMass, G4double aCharge, - G4ThreeVector aMomentum, G4ThreeVector aVertex, - G4int aneuIdx, G4int aneuPDG, TLorentzVector aneuP4, TLorentzVector aneuX4, - G4int aInttype, G4int aScatteringtype, G4double aW, - G4int afslPDG, TLorentzVector afslP4 ) : + G4ThreeVector aMomentum, G4ThreeVector aVertex) : fPartID(aID), fPDG(aPDG), fMass(aMass), fCharge{aCharge}, - fMomentumMC(aMomentum), fVertexMC(aVertex), - fNeuIdx(aneuIdx), fNeuPDG(aneuPDG), fNeuP4(aneuP4), fNeuX4(aneuX4), - fInteractionType(aInttype), fScatteringType(aScatteringtype), fW(aW), - fFSLPDG(afslPDG), fFSLP4(afslP4) -{ -} + fMomentumMC(aMomentum), fVertexMC(aVertex) +{} PrimaryParticleInformation::~PrimaryParticleInformation() -{ -} +{} void PrimaryParticleInformation::Print() const { G4cout<SetBranchAddress("pzf",&m_pzf); // hadrons pz (GeV) fGSTTree->SetBranchAddress("W",&m_W); // invariant hadronic mass (GeV) + fGSTTree->SetBranchAddress("Q2",&m_Q2); // momentum transfer (GeV^2) + fGSTTree->SetBranchAddress("x",&m_x); // Bjorken x + fGSTTree->SetBranchAddress("y",&m_y); // inelasticity + + fGSTTree->SetBranchAddress("wght",&m_wght); // event weigth + + fGSTTree->SetBranchAddress("tgt",&m_tgt); // nuclear target pdg + fGSTTree->SetBranchAddress("Z",&m_Z); // nuclear target Z + fGSTTree->SetBranchAddress("A",&m_A); // nuclear target A + fGSTTree->SetBranchAddress("hitnuc",&m_hitnuc); // hit nucleon pfg } @@ -145,41 +157,37 @@ void GENIEGenerator::GeneratePrimaries(G4Event* anEvent) // compute/repackage what is not directly available from the tree // position is randomly extracted in the detector fiducial volume - // or set to the center according to config parameter - - m_neuIdx = currentIdx; - m_neuP4.SetPxPyPzE(m_pxv,m_pyv,m_pzv,m_Ev); - m_fslP4.SetPxPyPzE(m_pxl,m_pyl,m_pzl,m_El); - - m_int_type = DecodeInteractionType(); - m_scattering_type = DecodeScatteringType(); + // or set to the center according to config parameter + TLorentzVector neuP4(m_pxv,m_pyv,m_pzv,m_Ev); + TLorentzVector fslP4(m_pxl,m_pyl,m_pzl,m_El); + TLorentzVector neuX4; G4Random::setTheSeed(currentIdx+1); if(fRandomVtx){ - m_neuX4.SetX(GeometricalParameters::Get()->GetFLArEPosition().x() + - (G4UniformRand()-0.5) * GeometricalParameters::Get()->GetFLArEFidVolSize().x()); - m_neuX4.SetY(GeometricalParameters::Get()->GetFLArEPosition().y() + - (G4UniformRand()-0.5) * GeometricalParameters::Get()->GetFLArEFidVolSize().y()); - m_neuX4.SetZ(GeometricalParameters::Get()->GetFLArEPosition().z() + - (G4UniformRand()-0.5) * GeometricalParameters::Get()->GetFLArEFidVolSize().z()); - m_neuX4.SetT(0.); + neuX4.SetX(GeometricalParameters::Get()->GetFLArEPosition().x() + + (G4UniformRand()-0.5) * GeometricalParameters::Get()->GetFLArEFidVolSize().x()); + neuX4.SetY(GeometricalParameters::Get()->GetFLArEPosition().y() + + (G4UniformRand()-0.5) * GeometricalParameters::Get()->GetFLArEFidVolSize().y()); + neuX4.SetZ(GeometricalParameters::Get()->GetFLArEPosition().z() + + (G4UniformRand()-0.5) * GeometricalParameters::Get()->GetFLArEFidVolSize().z()); + neuX4.SetT(0.); } else { - m_neuX4.SetX(0.*m); - m_neuX4.SetY(0.*m); - m_neuX4.SetZ(GeometricalParameters::Get()->GetFLArEPosition().z() - - GeometricalParameters::Get()->GetFLArEFidVolSize().z()/2); - m_neuX4.SetT(0.); + neuX4.SetX(0.*m); + neuX4.SetY(0.*m); + neuX4.SetZ(GeometricalParameters::Get()->GetFLArEPosition().z() - + GeometricalParameters::Get()->GetFLArEFidVolSize().z()/2); + neuX4.SetT(0.); } // create primary vertex (neutrino) - G4PrimaryVertex* vtx = new G4PrimaryVertex(m_neuX4.X(), m_neuX4.Y(), m_neuX4.Z(), m_neuX4.T()); // in mm + G4PrimaryVertex* vtx = new G4PrimaryVertex(neuX4.X(), neuX4.Y(), neuX4.Z(), neuX4.T()); // in mm // now add all the final state particles into the vertex // - final state lepton (if NC, it's the neutrino!) G4ParticleDefinition* particleDefinition; if ( FindParticleDefinition(m_fslPDG, particleDefinition) ){ - G4PrimaryParticle* plepton = new G4PrimaryParticle(particleDefinition, m_fslP4.X()*GeV, m_fslP4.Y()*GeV, m_fslP4.Z()*GeV, m_fslP4.E()*GeV); //in GeV + G4PrimaryParticle* plepton = new G4PrimaryParticle(particleDefinition, fslP4.X()*GeV, fslP4.Y()*GeV, fslP4.Z()*GeV, fslP4.E()*GeV); //in GeV /* G4cout << "Lepton PDG " << m_fslPDG << " mass " << particleDefinition->GetPDGMass()*MeV << G4endl; G4cout << "p4 " << m_fslP4.X() << " " << m_fslP4.Y() << " " << m_fslP4.Z() << " " << m_fslP4.E() << G4endl; G4cout << "kinE " << ( m_fslP4.E()*GeV - particleDefinition->GetPDGMass()*MeV) << G4endl; */ @@ -200,6 +208,29 @@ void GENIEGenerator::GeneratePrimaries(G4Event* anEvent) } + // package and ship metadata + GeneratorVertexMetadata metadata; + metadata.generatorType = fGeneratorName; + metadata.processName = EncodeProcessName(); + metadata.weight = m_wght; + metadata.pdg = m_neuPDG; + metadata.x4 = neuX4; + metadata.p4 = neuP4; + metadata.mass = 0.; // neutrinos always massless + metadata.charge = 0.; // neutrinos always no charge + metadata.intType = DecodeInteractionType(); + metadata.scatteringType = DecodeScatteringType(); + metadata.fsl_pdg = m_fslPDG; + metadata.tgt_pdg = m_tgt; + metadata.tgt_A = m_A; + metadata.tgt_Z = m_Z; + metadata.hitnuc_pdg = m_hitnuc; + metadata.Q2 = m_Q2; + metadata.xBj = m_x; + metadata.y = m_y; + metadata.W = m_W; + fVertexMetadata.push_back(metadata); + anEvent->AddPrimaryVertex(vtx); fEventCounter++; @@ -278,3 +309,21 @@ G4int GENIEGenerator::DecodeScatteringType() const return kScUnknown; } + +G4String GENIEGenerator::EncodeProcessName() const +{ + G4String process = ""; + if (m_cc) process += "WeakCC"; + else if (m_nc) process += "WeakNC"; + else if (m_em) process += "EM"; + + if(m_qel) process += " QE"; + else if(m_res) process += " RES"; + else if(m_dis) process += " DIS"; + else if(m_coh) process += " COH"; + else if(m_imd) process += " IMD"; + else if(m_mec) process += " MEC"; + else if(m_nuel) process += " nuELASTIC"; + + return process; +} \ No newline at end of file From d2688ece5c59520ff9e50fbe26f6916910c89da6 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 25 Aug 2025 17:19:32 -0400 Subject: [PATCH 3/8] removing PrimaryParticleInformation --- include/AnalysisManager.hh | 3 +- include/PrimaryParticleInformation.hh | 67 --------------------------- src/AnalysisManager.cc | 34 ++++++++------ src/EventInformation.cc | 5 +- src/PrimaryGeneratorAction.cc | 25 ---------- src/PrimaryParticleInformation.cc | 19 -------- 6 files changed, 23 insertions(+), 130 deletions(-) delete mode 100644 include/PrimaryParticleInformation.hh delete mode 100644 src/PrimaryParticleInformation.cc diff --git a/include/AnalysisManager.hh b/include/AnalysisManager.hh index 5693994..672a3a1 100644 --- a/include/AnalysisManager.hh +++ b/include/AnalysisManager.hh @@ -14,7 +14,6 @@ #include "AnalysisManagerMessenger.hh" #include "PixelMap3D.hh" #include "FPFParticle.hh" -#include "FPFNeutrino.hh" #include "hep_hpc/hdf5/File.hpp" @@ -122,7 +121,7 @@ class AnalysisManager { //--------------------------------------------------- // OUTPUT VARIABLES FOR COMMON TREES - G4int evtID; + G4int evtID; G4int vertexID; double weight; std::string genType; diff --git a/include/PrimaryParticleInformation.hh b/include/PrimaryParticleInformation.hh deleted file mode 100644 index 4bb5592..0000000 --- a/include/PrimaryParticleInformation.hh +++ /dev/null @@ -1,67 +0,0 @@ -#ifndef PrimaryParticleInformation_HH -#define PrimaryParticleInformation_HH - -#include -#include "G4VUserPrimaryParticleInformation.hh" -#include "G4ThreeVector.hh" -#include "globals.hh" - -// Primary particle information -class PrimaryParticleInformation : public G4VUserPrimaryParticleInformation { - - public: - /// a constructor - /// @param aID A unique particle ID within event. - /// @param aPDG A PDG code of the particle. - /// @param aMass A mass of the particle. - /// @param aMomentum An initial particle momentum (at the primary vertex). - /// @param aVertex An initial particle vertex. - PrimaryParticleInformation(G4int aID, G4int aPDG, G4double aMass, G4double aCharge, - G4ThreeVector aMomentum, G4ThreeVector aVertex); - - virtual ~PrimaryParticleInformation(); - - /// Get the particle unique ID (within event). Can be set only in the constructor. - inline G4int GetPartID() const { return fPartID; }; - - /// Get the particle PDG code - inline G4int GetPDG() const { return fPDG; }; - - /// Get the particle charge - inline G4double GetCharge() const { return fCharge; }; - - /// Get the particle mass in MeV - inline G4double GetMass() const { return fMass; }; - - /// Get the particle initial momentum. - inline G4ThreeVector GetMomentumMC() const { return fMomentumMC; }; - - /// Get the particle initial vertex. - inline G4ThreeVector GetVertexMC() const { return fVertexMC; }; - - /// Prints the information about the particle. - virtual void Print() const; - - private: - - /// A particle unique ID. - G4int fPartID; - - /// A particle type (PDG code). - G4int fPDG; - - /// A particle mass - G4double fMass; - - /// A particle initial momentum (from particle generator) - G4ThreeVector fMomentumMC; - - /// A particle initial vertex (from particle generator) - G4ThreeVector fVertexMC; - - /// particle's charge - G4double fCharge; - -}; - -#endif diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index 3a472f2..886ee6a 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -26,7 +26,6 @@ #include "LArBoxSD.hh" #include "FASER2TrackerSD.hh" #include "LArBoxHit.hh" -#include "PrimaryParticleInformation.hh" #include "EventInformation.hh" #include "generators/GeneratorVertexMetadata.hh" #include "reco/PCAAnalysis3D.hh" @@ -579,11 +578,9 @@ void AnalysisManager::FillPrimariesTree(const G4Event *event) G4PrimaryParticle *primary_particle = event->GetPrimaryVertex(ivtx)->GetPrimary(ipp); if (primary_particle) { - PrimaryParticleInformation *primary_particle_info = dynamic_cast(primary_particle->GetUserInformation()); - primary_particle_info->Print(); - + primVtxID = ivtx; - primTrackID = primary_particle_info->GetPartID(); // confirm matches track's id later? + primTrackID = ipp; // confirm matches track id? auto particleId = ActsFatras::Barcode(); particleId.setVertexPrimary(ivtx); @@ -592,16 +589,16 @@ void AnalysisManager::FillPrimariesTree(const G4Event *event) particleId.setParticle(primTrackID - 1); primParticleID = particleId.value(); - primPDG = primary_particle_info->GetPDG(); - primVx = primary_particle_info->GetVertexMC().x(); - primVy = primary_particle_info->GetVertexMC().y(); - primVz = primary_particle_info->GetVertexMC().z(); - primVt = 0; - primPx = primary_particle_info->GetMomentumMC().x(); - primPy = primary_particle_info->GetMomentumMC().y(); - primPz = primary_particle_info->GetMomentumMC().z(); - primM = primary_particle_info->GetMass()/MeV; - primQ = primary_particle_info->GetCharge(); + 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); TLorentzVector p4(primPx,primPy,primPz,energy); @@ -619,6 +616,13 @@ void AnalysisManager::FillPrimariesTree(const G4Event *event) 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(); } } diff --git a/src/EventInformation.cc b/src/EventInformation.cc index 6aed599..54d9286 100644 --- a/src/EventInformation.cc +++ b/src/EventInformation.cc @@ -19,8 +19,9 @@ void EventInformation::Print() const G4cout << "EventInformation: " << fGenMetadata.size() << " vertex(es)" << G4endl; for(int i=0; iSetUserInformation(new EventInformation(fGenerator->GetEventMetadata())); - // save additional truth information alongside primary particles - // loop over the vertices, and then over primary particles, - // and for each primary particle create an info object - G4int count_particles = 0; - for (G4int ivtx = 0; ivtx < anEvent->GetNumberOfPrimaryVertex(); ++ivtx) { - for (G4int ipp = 0; ipp < anEvent->GetPrimaryVertex(ivtx)->GetNumberOfParticle(); ++ipp) { - - G4PrimaryParticle* primary_particle = anEvent->GetPrimaryVertex(ivtx)->GetPrimary(ipp); - - if (primary_particle) { - primary_particle->SetUserInformation(new PrimaryParticleInformation( - count_particles, - primary_particle->GetPDGcode(), - primary_particle->GetMass(), - primary_particle->GetCharge(), - primary_particle->GetMomentum(), - anEvent->GetPrimaryVertex(ivtx)->GetPosition() - )); - - count_particles++; - } - } - } - } diff --git a/src/PrimaryParticleInformation.cc b/src/PrimaryParticleInformation.cc deleted file mode 100644 index cbf6441..0000000 --- a/src/PrimaryParticleInformation.cc +++ /dev/null @@ -1,19 +0,0 @@ -#include "PrimaryParticleInformation.hh" - -PrimaryParticleInformation::PrimaryParticleInformation(G4int aID, G4int aPDG, G4double aMass, G4double aCharge, - G4ThreeVector aMomentum, G4ThreeVector aVertex) : - fPartID(aID), fPDG(aPDG), fMass(aMass), fCharge{aCharge}, - fMomentumMC(aMomentum), fVertexMC(aVertex) -{} - -PrimaryParticleInformation::~PrimaryParticleInformation() -{} - -void PrimaryParticleInformation::Print() const { - G4cout< Date: Mon, 25 Aug 2025 17:22:45 -0400 Subject: [PATCH 4/8] remove FPFNeutrino --- CMakeLists.txt | 6 ++-- include/FPFNeutrino.hh | 73 ------------------------------------------ include/LinkDef.h | 1 - src/FPFNeutrino.cc | 71 ---------------------------------------- 4 files changed, 2 insertions(+), 149 deletions(-) delete mode 100644 include/FPFNeutrino.hh delete mode 100644 src/FPFNeutrino.cc diff --git a/CMakeLists.txt b/CMakeLists.txt index b7e5dd8..f983040 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -90,14 +90,12 @@ endforeach() #---------------------------------------------------------------------------- # generate dictionaries for the custom classes in the ROOT tree -ROOT_GENERATE_DICTIONARY(FPFClassesDict ${PROJECT_SOURCE_DIR}/include/FPFNeutrino.hh - ${PROJECT_SOURCE_DIR}/include/FPFParticle.hh +ROOT_GENERATE_DICTIONARY(FPFClassesDict ${PROJECT_SOURCE_DIR}/include/FPFParticle.hh LINKDEF ${PROJECT_SOURCE_DIR}/include/LinkDef.h MODULE FPFClasses) # create a shared library that includes the dictionary -add_library(FPFClasses SHARED ${PROJECT_SOURCE_DIR}/src/FPFNeutrino.cc - ${PROJECT_SOURCE_DIR}/src/FPFParticle.cc +add_library(FPFClasses SHARED ${PROJECT_SOURCE_DIR}/src/FPFParticle.cc FPFClassesDict.cxx) # include directories for headers diff --git a/include/FPFNeutrino.hh b/include/FPFNeutrino.hh deleted file mode 100644 index 571d711..0000000 --- a/include/FPFNeutrino.hh +++ /dev/null @@ -1,73 +0,0 @@ -#ifndef FPFNEUTRINO_HH -#define FPFNEUTRINO_HH - -class FPFNeutrino -{ - public: - - FPFNeutrino(); - - FPFNeutrino(const int pdg, - const double Vx, - const double Vy, - const double Vz, - const double T, - const double Px, - const double Py, - const double Pz, - const double E, - const int nuIdx, - const int nuIntType, - const int nuScatteringType, - const double nuW, - const int nuFSLPDG, - const double nuFSLPx, - const double nuFSLPy, - const double nuFSLPz, - const double nuFSLE); - - ~FPFNeutrino(); - - public: - // methods to access data members and other information - int NuPDG() const { return fNuPDG; } - double NuVx() const { return fNuVx; } - double NuVy() const { return fNuVy; } - double NuVz() const { return fNuVz; } - double NuT() const { return fNuT; } - double NuPx() const { return fNuPx; } - double NuPy() const { return fNuPy; } - double NuPz() const { return fNuPz; } - double NuE() const { return fNuE; } - int NuIdx() const { return fNuIdx; } - int NuIntType() const { return fNuIntType; } - int NuScatteringType() const { return fNuScatteringType;} - double NuW() const { return fNuW; } - int NuFSLPDG() const { return fNuFSLPDG; } - double NuFSLPx() const { return fNuFSLPx; } - double NuFSLPy() const { return fNuFSLPy; } - double NuFSLPz() const { return fNuFSLPz; } - double NuFSLE() const { return fNuFSLE; } - - private: - int fNuPDG; - double fNuVx; - double fNuVy; - double fNuVz; - double fNuT; - double fNuPx; - double fNuPy; - double fNuPz; - double fNuE; - int fNuIdx; - int fNuIntType; - int fNuScatteringType; - double fNuW; - int fNuFSLPDG; - double fNuFSLPx; - double fNuFSLPy; - double fNuFSLPz; - double fNuFSLE; -}; - -#endif diff --git a/include/LinkDef.h b/include/LinkDef.h index a5916d9..147adf5 100644 --- a/include/LinkDef.h +++ b/include/LinkDef.h @@ -1,5 +1,4 @@ #ifdef __CLING__ -#pragma link C++ class FPFNeutrino+; #pragma link C++ class FPFParticle+; #pragma link C++ class vector+; #endif diff --git a/src/FPFNeutrino.cc b/src/FPFNeutrino.cc deleted file mode 100644 index 2b2abec..0000000 --- a/src/FPFNeutrino.cc +++ /dev/null @@ -1,71 +0,0 @@ -#include -#include -#include - -#include "FPFNeutrino.hh" - -FPFNeutrino::FPFNeutrino() - : fNuPDG(0) - , fNuVx(-999) - , fNuVy(-999) - , fNuVz(-999) - , fNuT(-999) - , fNuPx(-999) - , fNuPy(-999) - , fNuPz(-999) - , fNuE(-999) - , fNuIdx(0) - , fNuIntType(0) - , fNuScatteringType(0) - , fNuW(-999) - , fNuFSLPDG(0) - , fNuFSLPx(-999) - , fNuFSLPy(-999) - , fNuFSLPz(-999) - , fNuFSLE(-999) -{ -} - -FPFNeutrino::FPFNeutrino( int pdg, - double Vx, - double Vy, - double Vz, - double T, - double Px, - double Py, - double Pz, - double E, - int nuIdx, - int nuIntType, - int nuScatteringType, - double nuW, - int nuFSLPDG, - double nuFSLPx, - double nuFSLPy, - double nuFSLPz, - double nuFSLE) - : fNuPDG(pdg) - , fNuVx(Vx) - , fNuVy(Vy) - , fNuVz(Vz) - , fNuT(T) - , fNuPx(Px) - , fNuPy(Py) - , fNuPz(Pz) - , fNuE(E) - , fNuIdx(nuIdx) - , fNuIntType(nuIntType) - , fNuScatteringType(nuScatteringType) - , fNuW(nuW) - , fNuFSLPDG(nuFSLPDG) - , fNuFSLPx(nuFSLPx) - , fNuFSLPy(nuFSLPy) - , fNuFSLPz(nuFSLPz) - , fNuFSLE(nuFSLE) -{ -} - -FPFNeutrino::~FPFNeutrino() -{ -} - From 28d6ce1f7e8d1dbb3e98be422211baae28972c70 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 25 Aug 2025 17:30:50 -0400 Subject: [PATCH 5/8] drop FPFClasses dictionary generation --- CMakeLists.txt | 19 +------------------ include/LinkDef.h | 4 ---- 2 files changed, 1 insertion(+), 22 deletions(-) delete mode 100644 include/LinkDef.h diff --git a/CMakeLists.txt b/CMakeLists.txt index f983040..3637c15 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -58,7 +58,6 @@ link_directories("$ENV{hep_hpc_path}/lib") set(HEP_HPC_LIBRARIES "$ENV{hep_hpc_path}/lib/libhep_hpc_hdf5.so;$ENV{hep_hpc_path}/lib/libhep_hpc_Utilities.so;$ENV{hep_hpc_path}/lib/libhep_hpc_concat_hdf5.so") message(STATUS "HEP_HPC lib: ${HEP_HPC_LIBRARIES}") - # Locate sources and headers for this project # We presume the existence of three directories @@ -88,19 +87,6 @@ foreach(_file ${macros} ${genie} ${analysis} ${gridutils} ${setup}) #message(STATUS "Copying ${_file} in ${PROJECT_BINARY_DIR}/${_file}") endforeach() -#---------------------------------------------------------------------------- -# generate dictionaries for the custom classes in the ROOT tree -ROOT_GENERATE_DICTIONARY(FPFClassesDict ${PROJECT_SOURCE_DIR}/include/FPFParticle.hh - LINKDEF ${PROJECT_SOURCE_DIR}/include/LinkDef.h - MODULE FPFClasses) - -# create a shared library that includes the dictionary -add_library(FPFClasses SHARED ${PROJECT_SOURCE_DIR}/src/FPFParticle.cc - FPFClassesDict.cxx) - -# include directories for headers -target_include_directories(FPFClasses PUBLIC ${PROJECT_SOURCE_DIR}/include) - #---------------------------------------------------------------------------- # Add the executable, and link it to the Geant4 libraries # @@ -112,12 +98,9 @@ target_link_libraries(FPFSim ${HDF5_LIBRARIES} ${HEP_HPC_LIBRARIES} ${HEPMC3_LIBRARIES} ${HEPMC3_FIO_LIBRARIES} ${HEPMC3_LIB} - FPFClasses ) #---------------------------------------------------------------------------- # Install the executable to 'bin' directory under CMAKE_INSTALL_PREFIX # -install(TARGETS FPFSim DESTINATION ${PROJECT_SOURCE_DIR}/bin) -install(TARGETS FPFClasses DESTINATION ${PROJECT_SOURCE_DIR}/lib) -install(FILES ${PROJECT_SOURCE_DIR}/build/libFPFClasses_rdict.pcm DESTINATION ${PROJECT_SOURCE_DIR}/lib) +install(TARGETS FPFSim DESTINATION ${PROJECT_SOURCE_DIR}/bin) \ No newline at end of file diff --git a/include/LinkDef.h b/include/LinkDef.h deleted file mode 100644 index 147adf5..0000000 --- a/include/LinkDef.h +++ /dev/null @@ -1,4 +0,0 @@ -#ifdef __CLING__ -#pragma link C++ class FPFParticle+; -#pragma link C++ class vector+; -#endif From fec6f3b22071690d50b1f88e93c75600e3f171a0 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 25 Aug 2025 22:45:25 -0400 Subject: [PATCH 6/8] add support for other generators --- include/generators/BackgroundGenerator.hh | 5 ++- include/generators/GeneratorVertexMetadata.hh | 14 +++---- macros/backgrounds.mac | 6 +-- src/AnalysisManager.cc | 27 +++++++------- src/EventInformation.cc | 8 ++-- src/generators/BackgroundGenerator.cc | 37 +++++++++++++------ src/generators/GENIEGenerator.cc | 31 ++++++++-------- src/generators/GPSGenerator.cc | 18 ++++++++- src/generators/HepMCGenerator.cc | 12 ++++++ 9 files changed, 100 insertions(+), 58 deletions(-) diff --git a/include/generators/BackgroundGenerator.hh b/include/generators/BackgroundGenerator.hh index 8e5c259..0769a5e 100644 --- a/include/generators/BackgroundGenerator.hh +++ b/include/generators/BackgroundGenerator.hh @@ -3,9 +3,10 @@ #include "generators/GeneratorBase.hh" +#include "G4LorentzVector.hh" + #include "TFile.h" #include "TH3D.h" -#include "TLorentzVector.h" #include "globals.hh" class G4Event; @@ -45,7 +46,7 @@ class BackgroundGenerator : public GeneratorBase G4double TPC_driftTime_s = 187.5e-6; // 187.5 us for 30cm drift // specific internal functions - void ShootParticle(G4Event* anEvent, G4int pdg, TLorentzVector x4, TLorentzVector p4) const; + void ShootParticle(G4Event* anEvent, G4int pdg, G4LorentzVector x4, G4LorentzVector p4); void SampleDirectionCosines(double& xdircos, double& ydircos, double E) const; int ExtractBackgroundParticles() const; diff --git a/include/generators/GeneratorVertexMetadata.hh b/include/generators/GeneratorVertexMetadata.hh index 985c18b..c6f022f 100644 --- a/include/generators/GeneratorVertexMetadata.hh +++ b/include/generators/GeneratorVertexMetadata.hh @@ -2,7 +2,7 @@ #define GENERATOR_VERTEX_METDATA_HH #include -#include "TLorentzVector.h" +#include "G4LorentzVector.hh" // Metadata information on each generated vertex // Store input interaction/decay process information @@ -10,13 +10,13 @@ struct GeneratorVertexMetadata { std::string generatorType; ///< which generator class std::string processName; ///< which process - double weight; ///< process weight + double weight = 1.0; ///< process weight - int pdg; ///< initiator pdg - TLorentzVector x4; ///< initiator 4-position (vertex) - TLorentzVector p4; ///< initiator 4-momentum - double mass; ///< initiator mass - double charge; ///< initiator charge + int pdg = -1; ///< initiator pdg + G4LorentzVector x4 = G4LorentzVector(); ///< initiator 4-position (vertex) + G4LorentzVector p4 = G4LorentzVector(); ///< initiator 4-momentum + double mass = -1; ///< initiator mass + double charge = 0; ///< initiator charge int intType = -1; ///< interaction type int scatteringType = -1; ///< scattering type diff --git a/macros/backgrounds.mac b/macros/backgrounds.mac index a65515a..219576d 100644 --- a/macros/backgrounds.mac +++ b/macros/backgrounds.mac @@ -13,12 +13,12 @@ /out/flare/save3DEvd false /out/flare/save2DEvd false -/out/addDiffusion false +/out/flare/addDiffusion false /out/fileName test_backgrounds.root # store and save trajectories for EVD -/tracking/storeTrajectory 1 -/out/saveTrack true +/tracking/storeTrajectory 0 +/out/saveTrack false # shoot background equivalent to 1 "spill" /run/beamOn 1 diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index 886ee6a..45b32fe 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -12,6 +12,7 @@ #include #include #include +#include #include #include @@ -532,14 +533,14 @@ void AnalysisManager::FillEventTree(const G4Event *event) genType = metadata[i].generatorType; processName = metadata[i].processName; initPDG = metadata[i].pdg; - initX = metadata[i].x4.X(); - initY = metadata[i].x4.Y(); - initZ = metadata[i].x4.Z(); - initT = metadata[i].x4.T(); - initPx = metadata[i].p4.X(); - initPy = metadata[i].p4.Y(); - initPz = metadata[i].p4.Z(); - initE = metadata[i].p4.E(); + initX = metadata[i].x4.x(); + initY = metadata[i].x4.y(); + initZ = metadata[i].x4.z(); + initT = metadata[i].x4.t(); + initPx = metadata[i].p4.x(); + initPy = metadata[i].p4.y(); + initPz = metadata[i].p4.z(); + initE = metadata[i].p4.e(); initM = metadata[i].mass; initQ = metadata[i].charge; intType = metadata[i].intType; @@ -601,11 +602,11 @@ void AnalysisManager::FillPrimariesTree(const G4Event *event) primQ = primary_particle->GetCharge(); G4double energy = GetTotalEnergy(primPx, primPy, primPz, primM); - TLorentzVector p4(primPx,primPy,primPz,energy); - primEta = p4.Eta(); - primPhi = p4.Phi(); - primPt = p4.Pt(); - primP = p4.P(); + G4LorentzVector p4(primPx,primPy,primPz,energy); + primEta = p4.eta(); + primPhi = p4.phi(); + primPt = p4.perp(); + primP = p4.vect().mag(); primE = energy; primKE = energy - primM; diff --git a/src/EventInformation.cc b/src/EventInformation.cc index 54d9286..5487745 100644 --- a/src/EventInformation.cc +++ b/src/EventInformation.cc @@ -24,10 +24,10 @@ void EventInformation::Print() const << "Generator : " << fGenMetadata[i].generatorType << G4endl << "Process name : " << fGenMetadata[i].processName << G4endl << "Initiator PDG : " << fGenMetadata[i].pdg << G4endl - << "Initiator p4 : (" << fGenMetadata[i].p4.X() << ", " << fGenMetadata[i].p4.Y() << ", " - << fGenMetadata[i].p4.Z() << ", " << fGenMetadata[i].p4.E() << ")" << G4endl - << "Initiator x4 : (" << fGenMetadata[i].x4.X() << ", " << fGenMetadata[i].x4.Y() << ", " - << fGenMetadata[i].x4.Z() << ", " << fGenMetadata[i].x4.T() << ")" << G4endl; + << "Initiator p4 : (" << fGenMetadata[i].p4.x() << ", " << fGenMetadata[i].p4.y() << ", " + << fGenMetadata[i].p4.z() << ", " << fGenMetadata[i].p4.e() << ")" << G4endl + << "Initiator x4 : (" << fGenMetadata[i].x4.x() << ", " << fGenMetadata[i].x4.y() << ", " + << fGenMetadata[i].x4.z() << ", " << fGenMetadata[i].x4.t() << ")" << G4endl; } } diff --git a/src/generators/BackgroundGenerator.cc b/src/generators/BackgroundGenerator.cc index 1ab1741..5d4b404 100644 --- a/src/generators/BackgroundGenerator.cc +++ b/src/generators/BackgroundGenerator.cc @@ -1,6 +1,7 @@ #include "generators/GeneratorBase.hh" #include "generators/BackgroundGenerator.hh" #include "generators/BackgroundGeneratorMessenger.hh" +#include "generators/GeneratorVertexMetadata.hh" #include "geometry/GeometricalParameters.hh" @@ -9,9 +10,9 @@ #include "G4IonTable.hh" #include "G4SystemOfUnits.hh" #include "G4Poisson.hh" -#include "Randomize.hh" +#include "G4LorentzVector.hh" -#include "TLorentzVector.h" +#include "Randomize.hh" #include "TMath.h" #include "TFile.h" #include "TTree.h" @@ -68,25 +69,37 @@ void BackgroundGenerator::LoadData() } -void BackgroundGenerator::ShootParticle(G4Event* anEvent, G4int pdg, TLorentzVector x4, TLorentzVector p4) const +void BackgroundGenerator::ShootParticle(G4Event* anEvent, G4int pdg, G4LorentzVector x4, G4LorentzVector p4) { // prepare a particle with the extracted starting position and momentum // once ready, shoot it with the gun! G4ParticleDefinition* particleDefinition = G4ParticleTable::GetParticleTable()->FindParticle(pdg); - /* uncomment for debugging.. + /* uncomment for debugging.. G4cout << "Particle PDG " << pdg << " mass " << particleDefinition->GetPDGMass()*MeV << G4endl; - G4cout << " x4 " << x4.X() << " " << x4.Y() << " " << x4.Z() << " " << x4.T() << G4endl; - G4cout << " p4 " << p4.X() << " " << p4.Y() << " " << p4.Z() << " " << p4.E() << G4endl; - G4cout << " kinE " << (p4.E() - particleDefinition->GetPDGMass()*MeV) << G4endl; + G4cout << " x4 " << x4.x() << " " << x4.y() << " " << x4.z() << " " << x4.t() << G4endl; + G4cout << " p4 " << p4.x() << " " << p4.y() << " " << p4.z() << " " << p4.e() << G4endl; + G4cout << " kinE " << (p4.e() - particleDefinition->GetPDGMass()*MeV) << G4endl; */ + // save metadata info... + GeneratorVertexMetadata metadata; + metadata.generatorType = fGeneratorName; + metadata.processName = "Background"; + metadata.weight = fBkgTimeWindow/s; + metadata.pdg = pdg; + metadata.x4 = x4; + metadata.p4 = p4; + metadata.mass = particleDefinition->GetPDGMass()*MeV; + metadata.charge = particleDefinition->GetPDGCharge(); + fVertexMetadata.push_back(metadata); + // load the gun... fGPS->SetParticleDefinition(particleDefinition); - fGPS->GetCurrentSource()->GetEneDist()->SetMonoEnergy(( p4.E() - particleDefinition->GetPDGMass()*MeV)); // kinetic energy - fGPS->GetCurrentSource()->GetAngDist()->SetParticleMomentumDirection(G4ThreeVector(p4.X(), p4.Y(), p4.Z())); + fGPS->GetCurrentSource()->GetEneDist()->SetMonoEnergy(( p4.e() - particleDefinition->GetPDGMass()*MeV)); // kinetic energy + fGPS->GetCurrentSource()->GetAngDist()->SetParticleMomentumDirection(G4ThreeVector(p4.x(), p4.y(), p4.z())); fGPS->GetCurrentSource()->GetPosDist()->SetPosDisType("Point"); - fGPS->GetCurrentSource()->GetPosDist()->SetCentreCoords( G4ThreeVector(x4.X(), x4.Y(), x4.Z()) ); + fGPS->GetCurrentSource()->GetPosDist()->SetCentreCoords( G4ThreeVector(x4.x(), x4.y(), x4.z()) ); fGPS->GeneratePrimaryVertex(anEvent); // ...and shoot! } @@ -181,12 +194,12 @@ void BackgroundGenerator::GeneratePrimaries(G4Event* anEvent) double zdircos = TMath::Sqrt(zdircos2); // always positive // now prepare inputs to load the gun - TLorentzVector x4(x*cm, y*cm, z, t); // tell G4 that zy are cm + G4LorentzVector x4(x*cm, y*cm, z, t); // tell G4 that zy are cm double mass = G4ParticleTable::GetParticleTable()->FindParticle(pdg)->GetPDGMass()*MeV; // in MeV double totE = E*GeV + mass; //total energy double p = TMath::Sqrt( totE*totE - mass*mass ); - TLorentzVector p4( xdircos*p, ydircos*p, zdircos*p, totE); + G4LorentzVector p4( xdircos*p, ydircos*p, zdircos*p, totE); // you may fire when ready now ShootParticle( anEvent, pdg, x4, p4); diff --git a/src/generators/GENIEGenerator.cc b/src/generators/GENIEGenerator.cc index 079564f..2f8198b 100644 --- a/src/generators/GENIEGenerator.cc +++ b/src/generators/GENIEGenerator.cc @@ -158,36 +158,35 @@ void GENIEGenerator::GeneratePrimaries(G4Event* anEvent) // compute/repackage what is not directly available from the tree // position is randomly extracted in the detector fiducial volume // or set to the center according to config parameter - TLorentzVector neuP4(m_pxv,m_pyv,m_pzv,m_Ev); - TLorentzVector fslP4(m_pxl,m_pyl,m_pzl,m_El); - TLorentzVector neuX4; + G4LorentzVector neuP4(m_pxv*GeV,m_pyv*GeV,m_pzv*GeV,m_Ev*GeV); + G4LorentzVector fslP4(m_pxl*GeV,m_pyl*GeV,m_pzl*GeV,m_El*GeV); + G4LorentzVector neuX4; G4Random::setTheSeed(currentIdx+1); if(fRandomVtx){ - neuX4.SetX(GeometricalParameters::Get()->GetFLArEPosition().x() + + neuX4.setX(GeometricalParameters::Get()->GetFLArEPosition().x() + (G4UniformRand()-0.5) * GeometricalParameters::Get()->GetFLArEFidVolSize().x()); - neuX4.SetY(GeometricalParameters::Get()->GetFLArEPosition().y() + + neuX4.setY(GeometricalParameters::Get()->GetFLArEPosition().y() + (G4UniformRand()-0.5) * GeometricalParameters::Get()->GetFLArEFidVolSize().y()); - neuX4.SetZ(GeometricalParameters::Get()->GetFLArEPosition().z() + + neuX4.setZ(GeometricalParameters::Get()->GetFLArEPosition().z() + (G4UniformRand()-0.5) * GeometricalParameters::Get()->GetFLArEFidVolSize().z()); - neuX4.SetT(0.); + neuX4.setT(0.); } else { - neuX4.SetX(0.*m); - neuX4.SetY(0.*m); - neuX4.SetZ(GeometricalParameters::Get()->GetFLArEPosition().z() - + neuX4.setX(0.*m); + neuX4.setY(0.*m); + neuX4.setZ(GeometricalParameters::Get()->GetFLArEPosition().z() - GeometricalParameters::Get()->GetFLArEFidVolSize().z()/2); - neuX4.SetT(0.); + neuX4.setT(0.); } // create primary vertex (neutrino) - G4PrimaryVertex* vtx = new G4PrimaryVertex(neuX4.X(), neuX4.Y(), neuX4.Z(), neuX4.T()); // in mm - + G4PrimaryVertex* vtx = new G4PrimaryVertex(neuX4.x(), neuX4.y(), neuX4.z(), neuX4.t()); // now add all the final state particles into the vertex // - final state lepton (if NC, it's the neutrino!) G4ParticleDefinition* particleDefinition; if ( FindParticleDefinition(m_fslPDG, particleDefinition) ){ - G4PrimaryParticle* plepton = new G4PrimaryParticle(particleDefinition, fslP4.X()*GeV, fslP4.Y()*GeV, fslP4.Z()*GeV, fslP4.E()*GeV); //in GeV + G4PrimaryParticle* plepton = new G4PrimaryParticle(particleDefinition, fslP4.x(), fslP4.y(), fslP4.z(), fslP4.e()); /* G4cout << "Lepton PDG " << m_fslPDG << " mass " << particleDefinition->GetPDGMass()*MeV << G4endl; G4cout << "p4 " << m_fslP4.X() << " " << m_fslP4.Y() << " " << m_fslP4.Z() << " " << m_fslP4.E() << G4endl; G4cout << "kinE " << ( m_fslP4.E()*GeV - particleDefinition->GetPDGMass()*MeV) << G4endl; */ @@ -198,9 +197,9 @@ void GENIEGenerator::GeneratePrimaries(G4Event* anEvent) // - all final state hadrons for (int ipar=0; iparGetPDGMass()*MeV << G4endl; G4cout << "p4 " << p.X() << " " << p.Y() << " " << p.Z() << " " << p.E() << G4endl; G4cout << "kinE " << (p.E()*GeV - particleDefinition->GetPDGMass()*MeV) << G4endl; */ diff --git a/src/generators/GPSGenerator.cc b/src/generators/GPSGenerator.cc index 43074ea..4a8b87d 100644 --- a/src/generators/GPSGenerator.cc +++ b/src/generators/GPSGenerator.cc @@ -1,11 +1,13 @@ #include "generators/GeneratorBase.hh" #include "generators/GPSGenerator.hh" +#include "generators/GeneratorVertexMetadata.hh" #include "geometry/GeometricalParameters.hh" #include "G4GeneralParticleSource.hh" #include "G4SystemOfUnits.hh" #include "G4ParticleTable.hh" +#include "G4LorentzVector.hh" #include "Randomize.hh" GPSGenerator::GPSGenerator() @@ -56,8 +58,22 @@ void GPSGenerator::GeneratePrimaries(G4Event* anEvent) { // complete line from PrimaryGeneratorAction.. G4cout << "): General Particle Source ===oooOOOooo===" << G4endl; - fGPS->GeneratePrimaryVertex(anEvent); + // preparing to ship metadata + GeneratorVertexMetadata metadata; + metadata.generatorType = fGeneratorName; + metadata.processName = "Gun"; + metadata.weight = 1.0; + metadata.pdg = fGPS->GetParticleDefinition()->GetPDGEncoding(); + metadata.mass = fGPS->GetParticleDefinition()->GetPDGMass(); + metadata.charge = fGPS->GetParticleDefinition()->GetPDGCharge(); + G4PrimaryVertex* vtx = anEvent->GetPrimaryVertex(); + G4PrimaryParticle* pp = vtx->GetPrimary(0); + metadata.x4 = G4LorentzVector(vtx->GetX0(),vtx->GetY0(),vtx->GetZ0(),vtx->GetT0()); + metadata.p4 = G4LorentzVector(pp->GetPx(),pp->GetPy(),pp->GetPz(),pp->GetTotalEnergy()); + fVertexMetadata.push_back(metadata); + + fGPS->GeneratePrimaryVertex(anEvent); } diff --git a/src/generators/HepMCGenerator.cc b/src/generators/HepMCGenerator.cc index 4bda094..3eeb111 100644 --- a/src/generators/HepMCGenerator.cc +++ b/src/generators/HepMCGenerator.cc @@ -1,6 +1,8 @@ #include "generators/GeneratorBase.hh" #include "generators/HepMCGenerator.hh" #include "generators/HepMCGeneratorMessenger.hh" +#include "generators/GeneratorVertexMetadata.hh" + #include "geometry/GeometricalParameters.hh" #include "HepMC3/ReaderAscii.h" @@ -168,6 +170,16 @@ void HepMCGenerator::HepMC2G4(const std::shared_ptr hepmcevt, } } + // FIXME: event header in original file is very simple + // need to save more useful metadata in event header! + GeneratorVertexMetadata metadata; + metadata.generatorType = fGeneratorName; + metadata.processName = "Decay"; + metadata.weight = hepmcevt->weights()[0]; + metadata.x4 = xvtx; + metadata.xs = hepmcevt->cross_section()->xsec(); + fVertexMetadata.push_back(metadata); + g4event->AddPrimaryVertex(g4vtx); } } From 23c3b7b930914e8c5aa3d6e89c027ba4c078f023 Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 25 Aug 2025 22:58:31 -0400 Subject: [PATCH 7/8] fix include --- src/generators/GENIEGenerator.cc | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/generators/GENIEGenerator.cc b/src/generators/GENIEGenerator.cc index 2f8198b..b75af37 100644 --- a/src/generators/GENIEGenerator.cc +++ b/src/generators/GENIEGenerator.cc @@ -11,12 +11,12 @@ #include "G4IonTable.hh" #include "G4SystemOfUnits.hh" #include "G4Exception.hh" +#include "G4LorentzVector.hh" #include "Randomize.hh" #include "TMath.h" #include "TFile.h" #include "TTree.h" -#include "TLorentzVector.h" GENIEGenerator::GENIEGenerator() { From 50ce1b2b1ccc8e615058809f19d9a9c45dac7fba Mon Sep 17 00:00:00 2001 From: Matteo Vicenzi Date: Mon, 15 Sep 2025 15:12:21 -0400 Subject: [PATCH 8/8] fix PixelMap3D issues --- src/AnalysisManager.cc | 90 ++++++++++++++++++++++-------------------- 1 file changed, 47 insertions(+), 43 deletions(-) diff --git a/src/AnalysisManager.cc b/src/AnalysisManager.cc index 45b32fe..ad9c7aa 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -581,7 +581,7 @@ void AnalysisManager::FillPrimariesTree(const G4Event *event) { primVtxID = ivtx; - primTrackID = ipp; // confirm matches track id? + primTrackID = ipp + 1; // confirm matches track id? auto particleId = ActsFatras::Barcode(); particleId.setVertexPrimary(ivtx); @@ -744,21 +744,22 @@ void AnalysisManager::FillFLArEOutput() // which primary ancestor does this hit belong to? G4int whichPrim = GetTrackPrimaryAncestor(flareTrackID); + G4int whichIndex = whichPrim - 1; // need to start from zero // pseudo reco: track/shower length and width - primaryTrackLength[whichPrim] += hit->GetStepLength(); - double ShowerP = primaries[whichPrim].P(); - double dsquare_hit_vtx = TMath::Power((flareX-primaries[whichPrim].Vx()),2)+ - TMath::Power((flareY-primaries[whichPrim].Vy()),2)+ - TMath::Power((flareZ-primaries[whichPrim].Vz()),2); - double product_hit_p = (flareX-primaries[whichPrim].Vx())*primaries[whichPrim].Px()+ - (flareY-primaries[whichPrim].Vy())*primaries[whichPrim].Py()+ - (flareZ-primaries[whichPrim].Vz())*primaries[whichPrim].Pz(); + 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[whichPrim] = (ShowerLength[whichPrim]>len_hit) ? ShowerLength[whichPrim] : len_hit; + ShowerLength[whichIndex] = (ShowerLength[whichIndex]>len_hit) ? ShowerLength[whichIndex] : len_hit; double weighted_width_hit = width_hit*flareEdep; - if (!std::isnan(weighted_width_hit)) ShowerWidth[whichPrim] += weighted_width_hit; + if (!std::isnan(weighted_width_hit)) ShowerWidth[whichIndex] += weighted_width_hit; if (sdName == "FLArEBoxSD/lar_box") { @@ -766,18 +767,18 @@ void AnalysisManager::FillFLArEOutput() fFLArEHits->Fill(); if (fAddDiffusion == "toy") - pm3D->FillEntryWithToyElectronTransportation(hit_position_xyz, vtx_xyz, flareEdep, whichPrim); + pm3D->FillEntryWithToyElectronTransportation(hit_position_xyz, vtx_xyz, flareEdep, whichIndex); else if (fAddDiffusion == "single") - pm3D->FillEntryWithToySingleElectronTransportation(hit_position_xyz, vtx_xyz, flareEdep, whichPrim); - else if (sdName == "lArBoxSD/lar_box") - pm3D->FillEntry(hit_position_xyz, vtx_xyz, flareEdep, whichPrim); + pm3D->FillEntryWithToySingleElectronTransportation(hit_position_xyz, vtx_xyz, flareEdep, whichIndex); + else + pm3D->FillEntry(hit_position_xyz, vtx_xyz, flareEdep, whichIndex); // accumulate Edep in LAr edepInLAr += flareEdep; - ProngEInLAr[whichPrim] += flareEdep; - primaryTrackLengthInTPC[whichPrim] += hit->GetStepLength(); - ShowerLengthInLAr[whichPrim] = (ShowerLengthInLAr[whichPrim]>len_hit) ? ShowerLengthInLAr[whichPrim] : len_hit; - if (!std::isnan(weighted_width_hit)) ShowerWidthInLAr[whichPrim] += weighted_width_hit; + 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; // accumulate dE/dx in LAr float_t longitudinal_distance_to_vtx = ((flareX-vtx_xyz[0])*primaries[0].Px()+ @@ -797,7 +798,7 @@ void AnalysisManager::FillFLArEOutput() // accumulate Edep in HCAL edepInHCALX += flareEdep; - ProngEInHadCal[whichPrim] += flareEdep; + ProngEInHadCal[whichIndex] += flareEdep; } else if (sdName == "FLArEHadCalYSD/lar_box" || sdName == "FLArEMuonFinderYSD/lar_box" || @@ -811,7 +812,7 @@ void AnalysisManager::FillFLArEOutput() // accumulate Edep in HCAL edepInHCALY += flareEdep; - ProngEInHadCal[whichPrim] += flareEdep; + ProngEInHadCal[whichIndex] += flareEdep; } } } @@ -985,38 +986,41 @@ void AnalysisManager::FillFLArEPseudoReco() for (auto iPrim : primaryIDs ) { - primaryPDG[iPrim] = primaries[iPrim].PDG(); + // trackIDs go from 1 to N, you need index: 0 to N-1 + G4int iiPrim = iPrim - 1; - float_t totProngE = ProngEInLAr[iPrim]+ProngEInHadCal[iPrim]; + primaryPDG[iiPrim] = primaries[iiPrim].PDG(); + + float_t totProngE = ProngEInLAr[iiPrim]+ProngEInHadCal[iiPrim]; if (totProngE>0) { - ShowerWidth[iPrim] = ShowerWidth[iPrim] / totProngE; + ShowerWidth[iiPrim] = ShowerWidth[iiPrim] / totProngE; } - if (ProngEInLAr[iPrim] > 0) + if (ProngEInLAr[iiPrim] > 0) { - ShowerWidthInLAr[iPrim] = ShowerWidthInLAr[iPrim] / ProngEInLAr[iPrim]; + ShowerWidthInLAr[iiPrim] = ShowerWidthInLAr[iiPrim] / ProngEInLAr[iiPrim]; } - double ShowerP = primaries[iPrim].P(); - double costheta = primaries[iPrim].Pz() / ShowerP; - ProngAngleToBeamDir[iPrim] = TMath::ACos(costheta); + double ShowerP = primaries[iiPrim].P(); + double costheta = primaries[iiPrim].Pz() / ShowerP; + ProngAngleToBeamDir[iiPrim] = TMath::ACos(costheta); - ProngAvgdEdx[iPrim] = (ProngEInLAr[iPrim] + - ProngEInHadCal[iPrim]) / - ShowerLength[iPrim]; - ProngAvgdEdxInLAr[iPrim] = ProngEInLAr[iPrim] / ShowerLengthInLAr[iPrim]; + ProngAvgdEdx[iiPrim] = (ProngEInLAr[iiPrim] + + ProngEInHadCal[iiPrim]) / + ShowerLength[iiPrim]; + ProngAvgdEdxInLAr[iiPrim] = ProngEInLAr[iiPrim] / ShowerLengthInLAr[iiPrim]; std::cout << std::setiosflags(std::ios::fixed) << std::setprecision(3); - std::cout << std::setw(10) << primaries[iPrim].PDG(); - std::cout << std::setw(12) << ProngAngleToBeamDir[iPrim]; - std::cout << std::setw(13) << primaryTrackLength[iPrim]; - std::cout << std::setw(13) << ShowerLength[iPrim]; - std::cout << std::setw(18) << ShowerWidthInLAr[iPrim]; - std::cout << std::setw(12) << ProngEInLAr[iPrim]; - std::cout << std::setw(12) << ProngEInHadCal[iPrim]; - std::cout << std::setw(12) << ProngAvgdEdxInLAr[iPrim]; - std::cout << std::setw(10) << primaries[iPrim].ProngType(); - std::cout << std::setw(12) << primaries[iPrim].Pz() << std::endl; + std::cout << std::setw(10) << primaries[iiPrim].PDG(); + std::cout << std::setw(12) << ProngAngleToBeamDir[iiPrim]; + std::cout << std::setw(13) << primaryTrackLength[iiPrim]; + std::cout << std::setw(13) << ShowerLength[iiPrim]; + std::cout << std::setw(18) << ShowerWidthInLAr[iiPrim]; + std::cout << std::setw(12) << ProngEInLAr[iiPrim]; + std::cout << std::setw(12) << ProngEInHadCal[iiPrim]; + std::cout << std::setw(12) << ProngAvgdEdxInLAr[iiPrim]; + std::cout << std::setw(10) << primaries[iiPrim].ProngType(); + std::cout << std::setw(12) << primaries[iiPrim].Pz() << std::endl; } slid::ShowerLID *shwlid = new slid::ShowerLID(pm3D->Get3DPixelMap(),