diff --git a/CMakeLists.txt b/CMakeLists.txt index b7e5dd8..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,21 +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/FPFNeutrino.hh - ${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 - 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 # @@ -114,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/AnalysisManager.hh b/include/AnalysisManager.hh index cbe7829..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" @@ -62,6 +61,7 @@ class AnalysisManager { void bookFLArETrees(); void bookFASER2Trees(); + void FillEventTree(const G4Event* event); void FillPrimariesTree(const G4Event* event); void FillTrajectoriesTree(const G4Event* event); @@ -120,13 +120,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 evtID; + 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 +180,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/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 deleted file mode 100644 index a5916d9..0000000 --- a/include/LinkDef.h +++ /dev/null @@ -1,5 +0,0 @@ -#ifdef __CLING__ -#pragma link C++ class FPFNeutrino+; -#pragma link C++ class FPFParticle+; -#pragma link C++ class vector+; -#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 deleted file mode 100644 index 767dc05..0000000 --- a/include/PrimaryParticleInformation.hh +++ /dev/null @@ -1,121 +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. - /// @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); - - 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; }; - - /// 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; - - 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; - - /// 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; - -}; - -#endif 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/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 48b9e85..de93f6e 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; } + G4String GetGeneratorName() const { return fGeneratorName; } + + // reset the list of metadata + void ResetEventMetadata() { fVertexMetadata.clear(); } + + // return full event vertex metadata + std::vector GetEventMetadata() const { return fVertexMetadata; } + // return single vertex metadata + GeneratorVertexMetadata GetEventMetadataPerVertex(G4int i) const { 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..c6f022f --- /dev/null +++ b/include/generators/GeneratorVertexMetadata.hh @@ -0,0 +1,38 @@ +#ifndef GENERATOR_VERTEX_METDATA_HH +#define GENERATOR_VERTEX_METDATA_HH + +#include +#include "G4LorentzVector.hh" + +// 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 = 1.0; ///< process weight + + 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 + 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 = -1.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/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 7c1e96b..ad9c7aa 100644 --- a/src/AnalysisManager.cc +++ b/src/AnalysisManager.cc @@ -12,6 +12,7 @@ #include #include #include +#include #include #include @@ -26,14 +27,14 @@ #include "LArBoxSD.hh" #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(); @@ -517,11 +579,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 + 1; // confirm matches track id? auto particleId = ActsFatras::Barcode(); particleId.setVertexPrimary(ivtx); @@ -530,23 +590,23 @@ 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); - 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; @@ -557,6 +617,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(); } } @@ -614,12 +681,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, @@ -682,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") { @@ -704,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()+ @@ -735,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" || @@ -749,14 +812,14 @@ void AnalysisManager::FillFLArEOutput() // accumulate Edep in HCAL edepInHCALY += flareEdep; - ProngEInHadCal[whichPrim] += flareEdep; + ProngEInHadCal[whichIndex] += flareEdep; } } } 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(); @@ -923,42 +986,45 @@ 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; + + primaryPDG[iiPrim] = primaries[iiPrim].PDG(); - float_t totProngE = ProngEInLAr[iPrim]+ProngEInHadCal[iPrim]; + 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(), - 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 +1035,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..5487745 --- /dev/null +++ b/src/EventInformation.cc @@ -0,0 +1,33 @@ +#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 -#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() -{ -} - diff --git a/src/PixelMap3D.cc b/src/PixelMap3D.cc index 5e6feb6..3e20e6f 100644 --- a/src/PixelMap3D.cc +++ b/src/PixelMap3D.cc @@ -259,7 +259,7 @@ void PixelMap3D::FillEntryWithToySingleElectronTransportation(const Double_t* po } } -void PixelMap3D::Process3DPM(File &h5file, FPFNeutrino nu, G4bool save3D) +void PixelMap3D::Process3DPM(File &h5file, G4int initPDG, G4int fslPDG, G4int intType, G4int scatType, G4double initE, G4bool save3D) { static auto evt_data = make_ntuple({h5file, "evt_data"}, @@ -270,8 +270,7 @@ void PixelMap3D::Process3DPM(File &h5file, FPFNeutrino nu, G4bool save3D) make_scalar_column("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 cd659b0..476d353 100644 --- a/src/PrimaryGeneratorAction.cc +++ b/src/PrimaryGeneratorAction.cc @@ -9,7 +9,7 @@ #include "geometry/GeometricalParameters.hh" -#include "PrimaryParticleInformation.hh" +#include "EventInformation.hh" #include "G4Event.hh" #include "G4Exception.hh" @@ -66,44 +66,13 @@ 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); - // 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(); - - // 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(), - neuidx, neupdg, neup4, neux4, int_type, scattering_type, w, - fslpdg, fslp4)); - - count_particles++; - } - } - } + // save vertex metadata information into the event + anEvent->SetUserInformation(new EventInformation(fGenerator->GetEventMetadata())); + } diff --git a/src/PrimaryParticleInformation.cc b/src/PrimaryParticleInformation.cc deleted file mode 100644 index 2413429..0000000 --- a/src/PrimaryParticleInformation.cc +++ /dev/null @@ -1,36 +0,0 @@ -#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 ) : - 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) -{ -} - -PrimaryParticleInformation::~PrimaryParticleInformation() -{ -} - -void PrimaryParticleInformation::Print() const { - G4cout<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 4e66234..b75af37 100644 --- a/src/generators/GENIEGenerator.cc +++ b/src/generators/GENIEGenerator.cc @@ -1,6 +1,7 @@ #include "generators/GeneratorBase.hh" #include "generators/GENIEGenerator.hh" #include "generators/GENIEGeneratorMessenger.hh" +#include "generators/GeneratorVertexMetadata.hh" #include "geometry/GeometricalParameters.hh" @@ -10,6 +11,7 @@ #include "G4IonTable.hh" #include "G4SystemOfUnits.hh" #include "G4Exception.hh" +#include "G4LorentzVector.hh" #include "Randomize.hh" #include "TMath.h" @@ -88,6 +90,16 @@ void GENIEGenerator::LoadData() fGSTTree->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,36 @@ 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 + 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){ - 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()); // 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(), 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; */ @@ -190,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; */ @@ -200,6 +207,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 +308,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 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); } }