diff --git a/CMakeLists.txt b/CMakeLists.txt index 3335c6e..df44485 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -38,17 +38,33 @@ find_package (Boost REQUIRED COMPONENTS log) find_package (fire REQUIRED) find_package (Eigen3 REQUIRED NO_MODULE) +set (biasing_sources + ${g4fire_SOURCE_DIR}/src/g4fire/biasing/DarkBrem.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/biasing/ElectroNuclear.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/biasing/GammaToMuPair.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/biasing/K0LongInelastic.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/biasing/NeutronInelastic.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/biasing/PhotoNuclear.cxx +) + set (dark_brem_sources - ${g4fire_SOURCE_DIR}/src/g4fire/DarkBrem/APrimePhysics.cxx - ${g4fire_SOURCE_DIR}/src/g4fire/DarkBrem/G4APrime.cxx - ${g4fire_SOURCE_DIR}/src/g4fire/DarkBrem/DarkBremVertexLibraryModel.cxx - ${g4fire_SOURCE_DIR}/src/g4fire/DarkBrem/G4eDarkBremsstrahlung.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/darkbrem/APrimePhysics.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/darkbrem/G4APrime.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/darkbrem/DarkBremVertexLibraryModel.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/darkbrem/G4eDarkBremsstrahlung.cxx +) + +set (event_sources + ${g4fire_SOURCE_DIR}/src/g4fire/event/EventBuilder.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/event/SimParticle.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/event/SimCalorimeterHit.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/event/SimTrackerHit.cxx ) set (geo_sources - ${g4fire_SOURCE_DIR}/src/g4fire/Geo/AuxInfoReader.cxx - ${g4fire_SOURCE_DIR}/src/g4fire/Geo/GDMLParser.cxx - ${g4fire_SOURCE_DIR}/src/g4fire/Geo/ParserFactory.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/geo/AuxInfoReader.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/geo/GDMLParser.cxx + ${g4fire_SOURCE_DIR}/src/g4fire/geo/ParserFactory.cxx ) set (sim_sources @@ -82,23 +98,57 @@ set (sim_sources ) add_library(g4fire SHARED + ${biasing_sources} ${dark_brem_sources} - ${sim_sources} + ${event_sources} ${geo_sources} + ${sim_sources} ) target_link_libraries(g4fire PUBLIC - Geant4::Interface fire::framework Eigen3::Eigen) -target_include_directories(g4fire PUBLIC ${g4fire_SOURCE_DIR}/include/) + Geant4::Interface fire::io fire::framework Eigen3::Eigen) + +target_include_directories(g4fire PUBLIC + "$" + "$" +) install(TARGETS g4fire + EXPORT g4fire-targets LIBRARY DESTINATION lib + RUNTIME DESTINATION bin INCLUDES DESTINATION include ) +# need pattern to exclude version header configured by CMake above +install(DIRECTORY include/ DESTINATION include/ FILES_MATCHING PATTERN "*.h") + +configure_file(${g4fire_SOURCE_DIR}/python/g4fire/_bias_operators.py.in + ${CMAKE_CURRENT_BINARY_DIR}/python/g4fire/_bias_operators.py) +install(FILES ${CMAKE_CURRENT_BINARY_DIR}/python/g4fire/_bias_operators.py + DESTINATION ${CMAKE_INSTALL_PREFIX}/python/g4fire) install(DIRECTORY python/ DESTINATION python FILES_MATCHING PATTERN "*.py" ) +install(EXPORT g4fire-targets + FILE g4fire-targets.cmake + NAMESPACE g4fire:: + DESTINATION lib/cmake/g4fire +) + +include(CMakePackageConfigHelpers) +configure_package_config_file( + "${CMAKE_CURRENT_SOURCE_DIR}/g4fireConfig.cmake.in" + "${CMAKE_CURRENT_BINARY_DIR}/g4fireConfig.cmake" + INSTALL_DESTINATION "lib/cmake/g4fire" +) + +install(FILES + "${CMAKE_CURRENT_BINARY_DIR}/g4fireConfig.cmake" + DESTINATION "lib/cmake/g4fire" +) + + # Unpack the example dark brem vertex library (or libraries) #file(GLOB vertex_libraries data/*.tar.gz) diff --git a/g4fireConfig.cmake.in b/g4fireConfig.cmake.in new file mode 100644 index 0000000..345b8d5 --- /dev/null +++ b/g4fireConfig.cmake.in @@ -0,0 +1,3 @@ +@PACKAGE_INIT@ + +include("${CMAKE_CURRENT_LIST_DIR}/g4fire-targets.cmake") diff --git a/include/g4fire/CalorimeterSD.h b/include/g4fire/CalorimeterSD.h deleted file mode 100644 index d82b4be..0000000 --- a/include/g4fire/CalorimeterSD.h +++ /dev/null @@ -1,79 +0,0 @@ -#ifndef SIMCORE_CALORIMETERSD_H -#define SIMCORE_CALORIMETERSD_H - -// Geant4 -#include "G4VSensitiveDetector.hh" - -// LDMX -#include "DetDescr/DetectorID.h" -#include "Framework/Event.h" -#include "g4fire/G4CalorimeterHit.h" - -using ldmx::DetectorID; - -namespace g4fire { - -/** - * @class CalorimeterSD - * @brief Basic calorimeter sensitive detector - */ -class CalorimeterSD : public G4VSensitiveDetector { - public: - /** - * Class constructor. - * @param name The name of the calorimeter. - * @param theCollectionName The name of the hits collection. - */ - CalorimeterSD(G4String name, G4String theCollectionName); - - /** - * Class destructor. - */ - virtual ~CalorimeterSD(); - - /** - * Set the depth into the geometric hierarchy to the layer volume. - * @param layerDepth The depth into the geometric hierarchy to the layer - * volume. - */ - void setLayerDepth(int layerDepth) { this->layerDepth_ = layerDepth; } - - /** - * Process a step by creating a hit. - * @param aStep The step information - * @param ROhist The readout history. - */ - virtual G4bool ProcessHits(G4Step* aStep, G4TouchableHistory* ROhist) = 0; - - /** - * Initialize the sensitive detector. - * @param hcEvent The hits collections of the event. - */ - void Initialize(G4HCofThisEvent* hcEvent); - - /** - * End of event hook. - * @param hcEvent The hits collections of the event. - */ - void EndOfEvent(G4HCofThisEvent* hcEvent); - - protected: - /** - * The hits collections of the sensitive detector. - */ - G4CalorimeterHitsCollection* hitsCollection_; - - /** - * The subdetector ID. - */ - int subdet_; - - /** - * The depth to the layer volume. - */ - int layerDepth_{2}; -}; - -} // namespace g4fire - -#endif diff --git a/include/g4fire/DetectorConstruction.h b/include/g4fire/DetectorConstruction.h index 5d6aebb..917fb48 100644 --- a/include/g4fire/DetectorConstruction.h +++ b/include/g4fire/DetectorConstruction.h @@ -7,7 +7,7 @@ #include "fire/config/Parameters.h" -#include "g4fire/Geo/Parser.h" +#include "g4fire/geo/Parser.h" namespace g4fire::geo { class Parser; diff --git a/include/g4fire/EcalHitIO.h b/include/g4fire/EcalHitIO.h deleted file mode 100644 index 29e0c34..0000000 --- a/include/g4fire/EcalHitIO.h +++ /dev/null @@ -1,102 +0,0 @@ - -#ifndef SIMCORE_ECALHITIO_H_ -#define SIMCORE_ECALHITIO_H_ - -// LDMX -#include "DetDescr/EcalHexReadout.h" -#include "g4fire/ConditionsInterface.h" -#include "g4fire/Event/SimCalorimeterHit.h" -#include "g4fire/G4CalorimeterHit.h" - -// STL -#include -#include - -namespace g4fire { - -/** - * @class EcalHitIO - * @brief Provides hit readout for simulated ECal detector - * - * @note - * This class uses the EcalHexReadout to transform the G4CalorimeterHit hits - * collection into a collection of SimCalorimeterHit objects assigned to - * hexagonal cells by their ID and positions. Energy depositions in the same - * cells are combined together into single hits. - * - * @par - * It can be run in three modes: - *
    - *
  • Compressed hit contributions combined by matching SimParticle and PDG - * code (default mode)
  • Full hit contributions with one record per step (when - * compressHitContribs_ is false)
  • No hit contribution information where - * energy is combined but vectors are not filled (when enableHitContribs_ is - * false) - *
- */ -class EcalHitIO { - public: - /** - * Layer-cell pair. - */ - typedef std::pair LayerCellPair; - - public: - /** - * Class constructor. - */ - EcalHitIO(ConditionsInterface& ci) : conditionsIntf_(ci) {} - - /** - * Class destructor. - */ - ~EcalHitIO() {} - - /** - * Configure the EcalHitIO using the passed parameters - */ - void configure(const framework::config::Parameters& ps); - - /** - * Write out a Geant4 hits collection to the provided ROOT array. - * @param hc The input hits collection. - * @param outputColl The output collection in ROOT. - */ - void writeHitsCollection(G4CalorimeterHitsCollection* hc, - std::vector& outputColl); - - /** - * Set whether hit contributions should be enabled for the output hits. - * @param enableHitContribs True to enable hit contributions. - */ - void setEnableHitContribs(bool enableHitContribs) { - enableHitContribs_ = enableHitContribs; - } - - /** - * Set whether hit contributions should be compressed by particle and PDG - * code. - * @param compressHitContribs True to compress hit contribution information. - */ - void setCompressHitContribs(bool compressHitContribs) { - compressHitContribs_ = compressHitContribs; - } - - private: - /// ConditionsInterface - ConditionsInterface& conditionsIntf_; - - /** - * Enable hit contribution output. - */ - bool enableHitContribs_{true}; - - /** - * Enable compression of hit contributions by SimParticle and PDG code. - */ - bool compressHitContribs_{true}; -}; - -} // namespace g4fire - -#endif diff --git a/include/g4fire/EcalSD.h b/include/g4fire/EcalSD.h deleted file mode 100644 index 93cf833..0000000 --- a/include/g4fire/EcalSD.h +++ /dev/null @@ -1,79 +0,0 @@ -/** - * @file EcalSD.h - * @brief Class defining an ECal sensitive detector using an EcalHexReadout to - * create the hits - * @author Jeremy McCormick, SLAC National Accelerator Laboratory - */ - -#ifndef SIMCORE_ECALSD_H_ -#define SIMCORE_ECALSD_H_ - -// LDMX -#include "DetDescr/EcalHexReadout.h" -#include "g4fire/CalorimeterSD.h" -#include "g4fire/ConditionsInterface.h" - -// ROOT -#include "TMath.h" - -// Geant4 -#include "G4Polyhedra.hh" - -namespace g4fire { - -/** - * @class EcalSD - * @brief ECal sensitive detector that uses an EcalHexReadout to create the hits - */ -class EcalSD : public CalorimeterSD { - public: - /** - * Class constructor. - * @param name The name of the sensitive detector. - * @param theCollectionName The name of the hits collection. - * @param subDetID The subdetector ID. - */ - EcalSD(G4String name, G4String theCollectionName, int subDetID, - ConditionsInterface& ci); - - /** - * Class destructor. - */ - virtual ~EcalSD(); - - /** - * Process steps to create hits. - * @param aStep The step information. - * @param ROhist The readout history. - */ - G4bool ProcessHits(G4Step* aStep, G4TouchableHistory* ROhist); - - private: - /** - * Return the hit position of a step. - * X and Y are computed from the midpoint of the step. - * Z corresponds to the volume's center. - * @return The hit position from the step. - * @todo This function is probably slow due to it inspecting the - * geometry to get the Z position so this should be sped up somehow. - */ - G4ThreeVector getHitPosition(G4Step* aStep); - - private: - /** - * The hex readout defining the cell grid. - */ - std::unique_ptr hitMap_; - - /** - * Map of polygonal layers for getting Z positions. - */ - std::map polyMap_; - - /// ConditionsInterface - ConditionsInterface& conditionsIntf_; -}; - -} // namespace g4fire - -#endif diff --git a/include/g4fire/Event/SimCalorimeterHit.h b/include/g4fire/Event/SimCalorimeterHit.h deleted file mode 100644 index c87b691..0000000 --- a/include/g4fire/Event/SimCalorimeterHit.h +++ /dev/null @@ -1,255 +0,0 @@ -/** - * @file SimCalorimeterHit.h - * @brief Class which stores simulated calorimeter hit information - * @author Jeremy McCormick, SLAC National Accelerator Laboratory - */ - -#ifndef SIMCORE_EVENT_SIMCALORIMETERHIT_H_ -#define SIMCORE_EVENT_SIMCALORIMETERHIT_H_ - -// ROOT -#include "TObject.h" //For ClassDef - -// LDMX -#include "g4fire/Event/SimParticle.h" - -namespace ldmx { - -/** - * @class SimCalorimeterHit - * @brief Stores simulated calorimeter hit information - * - * @note - * This class represents simulated hit information from a calorimeter detector. - * It provides access to the cell ID, energy deposition, cell position and time. - * Additionally, individual depositions or steps from MC particles are tabulated - * as contributions stored in vectors. Contribution information includes a - * reference to the relevant SimParticle, the PDG code of the actual particle - * which deposited energy (may be different from the actual SimParticle), the - * time of the contribution and the energy deposition. - */ -class SimCalorimeterHit { - public: - /** - * @class Contrib - * @brief Information about a contribution to the hit in the associated cell - */ - struct Contrib { - /** - * trackID of incident particle that is an ancestor of the contributor - * - * The incident ancestor is found in TrackMap::findIncident where the - * ancestry is looped upwards until a particle is found that matches - * the criteria. - * (1) particle will be saved to output file AND - * (2) particle originates in a region outside the CalorimeterRegion - * If no particle is found matching these criteria, the primary particle - * that is this trackID's ancestor is chosen. - */ - int incidentID{-1}; - - /// track ID of this contributor - int trackID{-1}; - - /// PDG ID of this contributor - int pdgCode{0}; - - /// Energy depostied by this contributor - float edep{0}; - - /// Time this contributor made the hit (global Geant4 time) - float time{0}; - }; - - /** - * Class constructor. - */ - SimCalorimeterHit(); - - /** - * Class destructor. - */ - virtual ~SimCalorimeterHit(); - - /** - * Clear the data in the object. - */ - void Clear(); - - /** - * Print out the object. - */ - void Print() const; - - /** - * Get the detector ID. - * @return The detector ID. - */ - int getID() const { return id_; } - - /** - * Set the detector ID. - * @id The detector ID. - */ - void setID(const int id) { this->id_ = id; } - - /** - * Get the energy deposition of the hit [MeV]. - * @return The energy deposition of the hit. - */ - float getEdep() const { return edep_; } - - /** - * Set the energy deposition of the hit [MeV]. - * @param edep The energy deposition of the hit. - */ - void setEdep(const float edep) { this->edep_ = edep; } - - /** - * Get the XYZ position of the hit [mm]. - * @return The XYZ position of the hit. - */ - std::vector getPosition() const { return {x_, y_, z_}; } - - /** - * Set the XYZ position of the hit [mm]. - * @param x The X position. - * @param y The Y position. - * @param z The Z position. - */ - void setPosition(const float x, const float y, const float z) { - this->x_ = x; - this->y_ = y; - this->z_ = z; - } - - /** - * Get the global time of the hit [ns]. - * @return The global time of the hit. - */ - float getTime() const { return time_; } - - /** - * Set the time of the hit [ns]. - * @param time The time of the hit. - */ - void setTime(const float time) { this->time_ = time; } - - /** - * Get the number of hit contributions. - * @return The number of hit contributions. - */ - unsigned getNumberOfContribs() const { return nContribs_; } - - /** - * Add a hit contribution from a SimParticle. - * @param incidentID the Geant4 track ID for the particle's parent incident on - * the Calorimeter region - * @param trackID the Geant4 track ID for the particle - * @param pdgCode The PDG code of the actual track. - * @param edep The energy deposition of the hit [MeV]. - * @param time The time of the hit [ns]. - */ - void addContrib(int incidentID, int trackID, int pdgCode, float edep, - float time); - - /** - * Get a hit contribution by index. - * @param i The index of the hit contribution. - * @return The hit contribution at the index. - */ - Contrib getContrib(int i) const; - - /** - * Find the index of a hit contribution from a SimParticle and PDG code. - * @param trackID the track ID of the particle causing the hit - * @param pdgCode The PDG code of the contribution. - * @return The index of the contribution or -1 if none exists. - */ - int findContribIndex(int trackID, int pdgCode) const; - - /** - * Update an existing hit contribution by incrementing its edep and setting - * the time if the new time is less than the old one. - * @param i The index of the contribution. - * @param edep The additional energy contribution [MeV]. - * @param time The time of the contribution [ns]. - */ - void updateContrib(int i, float edep, float time); - - /** - * Sort by time of hit - */ - bool operator<(const SimCalorimeterHit &rhs) const { - return this->getTime() < rhs.getTime(); - } - - private: - /** - * The detector ID. - */ - int id_{0}; - - /** - * The energy deposition. - */ - float edep_{0}; - - /** - * The X position. - */ - float x_{0}; - - /** - * The Y position. - */ - float y_{0}; - - /** - * The Z position. - */ - float z_{0}; - - /** - * The global time of the hit. - */ - float time_{0}; - - /** - * The list of track IDs contributing to the hit. - */ - std::vector trackIDContribs_; - - /** - * The list of incident IDs contributing to the hit - */ - std::vector incidentIDContribs_; - - /** - * The list of PDG codes contributing to the hit. - */ - std::vector pdgCodeContribs_; - - /** - * The list of energy depositions contributing to the hit. - */ - std::vector edepContribs_; - - /** - * The list of times contributing to the hit. - */ - std::vector timeContribs_; - - /** - * The number of hit contributions. - */ - unsigned nContribs_{0}; - - /** - * ROOT class definition. - */ - ClassDef(SimCalorimeterHit, 3) -}; -} // namespace ldmx - -#endif diff --git a/include/g4fire/Event/SimTrackerHit.h b/include/g4fire/Event/SimTrackerHit.h deleted file mode 100644 index b9f0da0..0000000 --- a/include/g4fire/Event/SimTrackerHit.h +++ /dev/null @@ -1,278 +0,0 @@ -/** - * @file SimTrackerHit.h - * @brief Class which encapsulates information from a hit in a simulated - * tracking detector - * @author Omar Moreno, SLAC National Accelerator Laboratory - * @author Jeremy McCormick, SLAC National Accelerator Laboratory - */ - -#ifndef SIMCORE_EVENT_SIMTRACKERHIT_H_ -#define SIMCORE_EVENT_SIMTRACKERHIT_H_ - -// ROOT -#include "TObject.h" //For ClassDef - -// STL -#include - -namespace ldmx { - -/** - * @class SimTrackerHit - * @brief Represents a simulated tracker hit in the simulation - */ -class SimTrackerHit { - public: - /** - * Class constructor. - */ - SimTrackerHit(); - - /** - * Class destructor. - */ - virtual ~SimTrackerHit(); - - /** - * Print a description of this object. - */ - void Print() const; - - /** - * Reset the SimTrackerHit object. - */ - void Clear(); - - /** - * Get the detector ID of the hit. - * @return The detector ID of the hit. - */ - int getID() const { return id_; }; - - /** - * Get the geometric layer ID of the hit. - * @return The layer ID of the hit. - */ - int getLayerID() const { return layerID_; }; - - /** - * Get the module ID associated with a hit. This is used to - * uniquely identify a sensor within a layer. - * @return The module ID associated with a hit. - */ - int getModuleID() const { return moduleID_; }; - - /** - * Get the XYZ position of the hit [mm]. - * @return The position of the hit. - */ - std::vector getPosition() const { return {x_, y_, z_}; }; - - /** - * Get the energy deposited on the hit [MeV]. - * @return The energy deposited on the hit. - */ - float getEdep() const { return edep_; }; - - /** - * Get the energy - * @return The energy of the hit. - */ - float getEnergy() const { return energy_; }; - - /** - * Get the global time of the hit [ns]. - * @return The global time of the hit. - */ - float getTime() const { return time_; }; - - /** - * Get the path length between the start and end points of the - * hit [mm]. - * @return The path length of the hit. - */ - float getPathLength() const { return pathLength_; }; - - /** - * Get the XYZ momentum of the particle at the position at which - * the hit took place [MeV]. - * @return The momentum of the particle. - */ - std::vector getMomentum() const { return {px_, py_, pz_}; }; - - /** - * Get the Sim particle track ID of the hit. - * @return The Sim particle track ID of the hit. - */ - int getTrackID() const { return trackID_; }; - - /** - * Get the Sim particle track ID of the hit. - * @return The Sim particle track ID of the hit. - */ - int getPdgID() const { return pdgID_; }; - - /** - * Set the detector ID of the hit. - * @param id The detector ID of the hit. - */ - void setID(const long id) { this->id_ = id; }; - - /** - * Set the geometric layer ID of the hit. - * @param layerID The layer ID of the hit. - */ - void setLayerID(const int layerID) { this->layerID_ = layerID; }; - - /** - * Set the module ID associated with a hit. This is used to - * uniquely identify a sensor within a layer. - * @return moduleID The module ID associated with a hit. - */ - void setModuleID(const int moduleID) { this->moduleID_ = moduleID; }; - - /** - * Set the position of the hit [mm]. - * @param x The X position. - * @param y The Y position. - * @param z The Z position. - */ - void setPosition(const float x, const float y, const float z); - - /** - * Set the energy deposited on the hit [MeV]. - * @param edep The energy deposited on the hit. - */ - void setEdep(const float edep) { this->edep_ = edep; }; - - /** - * Set the energy of the hit. - * @param e The energy of the hit. - */ - void setEnergy(const float energy) { energy_ = energy; }; - - /** - * Set the global time of the hit [ns]. - * @param time The global time of the hit. - */ - void setTime(const float time) { this->time_ = time; }; - - /** - * Set the path length of the hit [mm]. - * @param pathLength The path length of the hit. - */ - void setPathLength(const float pathLength) { - this->pathLength_ = pathLength; - }; - - /** - * Set the momentum of the particle at the position at which - * the hit took place [GeV]. - * @param px The X momentum. - * @param py The Y momentum. - * @param pz The Z momentum. - */ - void setMomentum(const float px, const float py, const float pz); - - /** - * Set the Sim particle track ID of the hit. - * @return The Sim particle track ID of the hit. - */ - void setTrackID(const int simTrackID) { this->trackID_ = simTrackID; }; - - /** - * Set the Sim particle track ID of the hit. - * @return The Sim particle track ID of the hit. - */ - void setPdgID(const int simPdgID) { this->pdgID_ = simPdgID; }; - - /** - * Sort by time of hit - */ - bool operator<(const ldmx::SimTrackerHit &rhs) const { - return this->getTime() < rhs.getTime(); - } - - private: - /** - * The detector ID. - */ - int id_{0}; - - /** - * The layer ID. - */ - int layerID_{0}; - - /** The module ID. */ - int moduleID_{0}; - - /** - * The energy deposited on the hit. - */ - float edep_{0}; - - /** - * The global time of the hit. - */ - float time_{0}; - - /** - * The X momentum. - */ - float px_{0}; - - /** - * The Y momentum. - */ - float py_{0}; - - /** - * The Z momentum. - */ - float pz_{0}; - - /** - * The total energy. - */ - float energy_{0}; - - /** - * The X position. - */ - float x_{0}; - - /** - * The Y position. - */ - float y_{0}; - - /** - * The Z position. - */ - float z_{0}; - - /** - * The path length of the hit. - */ - float pathLength_{0}; - - /** - * The Sim Track ID. - */ - int trackID_{0}; - - /** - * The Sim PDG ID. - */ - int pdgID_{0}; - - /** - * The ROOT class definition. - */ - ClassDef(SimTrackerHit, 3); - -}; // SimTrackerHit -} // namespace ldmx - -#endif // EVENT_SIMTRACKERHIT_H_ diff --git a/include/g4fire/HcalSD.h b/include/g4fire/HcalSD.h deleted file mode 100644 index 0f093fc..0000000 --- a/include/g4fire/HcalSD.h +++ /dev/null @@ -1,52 +0,0 @@ -#ifndef SIMCORE_HCALSD_H -#define SIMCORE_HCALSD_H - -/*~~~~~~~~~~~~~~*/ -/* DetDescr */ -/*~~~~~~~~~~~~~~*/ -#include "g4fire/CalorimeterSD.h" - -// Forward declarations -class G4Step; -class G4TouchableHistory; - -namespace g4fire { - -/** - * Class defining a sensitive detector of type HCal. - */ -class HcalSD : public CalorimeterSD { - public: - /** - * Class constructor. - * - * @param[in] name The namme of the sensitive detector. - * @param[in] collectionName Name of the colleciton of hits. - * #param[in] subDetID The subdetectorID - */ - HcalSD(G4String name, G4String collectionName, int subDetID); - - /// Destructor - ~HcalSD(); - - /** - * Create a hit out of the energy deposition deposited during a - * step. - * - * @param[in] step The current step. - * @param[in] history The readout history. - */ - G4bool ProcessHits(G4Step* aStep, G4TouchableHistory* ROhist); - - private: - // TODO: document! - double birksc1_; - - // TODO: document! - double birksc2_; - -}; // HcalSD - -} // namespace g4fire - -#endif // SIMCORE_HCALSD_H diff --git a/include/g4fire/ParallelWorld.h b/include/g4fire/ParallelWorld.h index e5cb330..dbb427d 100644 --- a/include/g4fire/ParallelWorld.h +++ b/include/g4fire/ParallelWorld.h @@ -7,7 +7,7 @@ #include "G4String.hh" #include "G4VUserParallelWorld.hh" -#include "g4fire/Geo/AuxInfoReader.h" +#include "g4fire/geo/AuxInfoReader.h" namespace g4fire { diff --git a/include/g4fire/Parameters.h b/include/g4fire/Parameters.h new file mode 100644 index 0000000..66d5b0c --- /dev/null +++ b/include/g4fire/Parameters.h @@ -0,0 +1,155 @@ +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include + +#include "fire/exception/Exception.h" + +namespace g4fire { + +/** + * Class encapsulating parameters for configuring a processor. + */ +class Parameters { +public: + /// Constructor + Parameters() = default; + + /// Destructor + ~Parameters() = default; + + /** + * Set the mapping of parameter names to value. + * + * @param [in, out] parameters mapping between parameter names and the + * corresponding value. + */ + void set(std::map parameters) { + parameters_ = parameters; + } + + /** + * Set the value of an existing parameter in the list. If the parameter + * does not exists, throw an exception. + * + * @param name Name of the parameter + * @param value The value of the parameter + * @throw Exception if the parameter is not in the list of parameters. + */ + template + void set(const std::string &name, const T &value) { + if (!exists(name)) { + throw fire::Exception("NonExistParam", + "Parameter '" + name + + "' does not exist in list of parameters.", + false); + } + parameters_[name] = value; + } + + /** + * Add a parameter to the parameter list. If the parameter already + * exists in the list, throw an exception. + * + * @param[in] name Name of the parameter. + * @param[in] value The value of the parameter. + * @throw Exception if a parameter by that name already exist in + * the list. + */ + template + void add(const std::string &name, const T &value) { + if (exists(name)) { + throw fire::Exception("ParameterExist", + "The parameter " + name + + " already exists in the list of parameters.", + false); + } + parameters_[name] = value; + } + + /** + * Check to see if a parameter exists + * + * @param[in] name name of parameter to check + * @return true if parameter exists in configuration set + */ + bool exists(const std::string &name) const { + return parameters_.find(name) != parameters_.end(); + } + + /** + * Retrieve the parameter of the given name. + * + * @throw Exception if parameter of the given name isn't found + * + * @throw Exception if parameter is found but not of the input type + * + * @param T the data type to cast the parameter to. + * + * @param[in] name the name of the parameter value to retrieve. + * + * @return The user specified parameter of type T. + */ + template T get(const std::string &name) const { + // Check if the variable exists in the map. If it doesn't, + // raise an exception. + if (not exists(name)) { + throw fire::Exception("NonExistParam", + "Parameter '" + name + + "' does not exist in list of parameters.", + false); + } + + try { + auto parameter = std::any_cast(parameters_.at(name)); + return parameter; + } catch (const std::bad_any_cast &e) { + throw fire::Exception("BadTypeParam", + "Parameter '" + name + "' of type '" + + parameters_.at(name).type().name() + + "' is being cast to incorrect type '" + + typeid(T).name() + "'.", + false); + } + } + + /** + * Retrieve a parameter with a default specified. + * + * Return the input default if a parameter is not found in map. + * + * @return the user parameter of type T + */ + template + T get(const std::string &name, const T &def) const { + if (not exists(name)) + return def; + + // get here knowing that name exists in parameters_ + return get(name); + } + + /** + * Get a list of the keys available. + * This may be helpful in debugging to make sure the parameters are spelled + * correctly. + */ + std::vector keys() const { + std::vector key; + for (auto i : parameters_) + key.push_back(i.first); + return key; + } + +private: + /// Parameters + std::map parameters_; + +}; // Parameters +} // namespace g4fire diff --git a/include/g4fire/Persist/RootPersistencyManager.h b/include/g4fire/Persist/RootPersistencyManager.h deleted file mode 100644 index 10668f0..0000000 --- a/include/g4fire/Persist/RootPersistencyManager.h +++ /dev/null @@ -1,210 +0,0 @@ -#ifndef SIMCORE_PERSIST_ROOTPERSISTENCYMANAGER_H_ -#define SIMCORE_PERSIST_ROOTPERSISTENCYMANAGER_H - -/*~~~~~~~~~~~~~~~~*/ -/* C++ StdLib */ -/*~~~~~~~~~~~~~~~~*/ -#include -#include - -/*~~~~~~~~~~~~~~~*/ -/* Framework */ -/*~~~~~~~~~~~~~~~*/ -#include "Framework/Configure/Parameters.h" -#include "Framework/EventFile.h" - -/*~~~~~~~~~~~~~*/ -/* g4fire */ -/*~~~~~~~~~~~~~*/ -#include "g4fire/EcalHitIO.h" -#include "g4fire/G4CalorimeterHit.h" -#include "g4fire/G4TrackerHit.h" - -/*~~~~~~~~~~~~*/ -/* Geant4 */ -/*~~~~~~~~~~~~*/ -#include "G4PersistencyCenter.hh" -#include "G4PersistencyManager.hh" - -// Forward declarations -class G4Run; - -// Forward declarations within the ldmx namespace -namespace framework { -class Event; -class RunHeader; -} // namespace framework - -namespace g4fire { -namespace persist { - -/** - * @class RootPersistencyManager - * - * @note - * Output is written at the end of each event. An EventFile is used to - * write from an Event buffer object into an output branch within a tree. - * The event buffer is cleared after the event is written. A - * SimParticleBuilder is used to build a set of SimParticle objects from - * the Trajectory objects which were created during event processing. An - * EcalHitIO instance provides translation of G4CalorimeterHit objects in - * the ECal to an output SimCalorimeterHit collection, transforming the - * individual steps into cell energy depositions. The tracker hit - * collections of G4TrackerHit objects are translated directly into - * output SimTrackerHit collections. - */ -class RootPersistencyManager : public G4PersistencyManager { - public: - /** - * Class constructor. - * - * @param eventFile file to put output events into - * @param parameters configuration parameters from Simulator - * @param runNumber current run identifer from Process - */ - RootPersistencyManager(framework::EventFile &file, - framework::config::Parameters ¶meters, - const int &runNumber, ConditionsInterface &ci); - - /// Destructor - virtual ~RootPersistencyManager() {} - - /** - * Get the current ROOT persistency manager or nullptr if not - * registered. - * - * @return The ROOT persistency manager. - */ - static RootPersistencyManager *getInstance() { - return static_cast( - G4PersistencyCenter::GetPersistencyCenter() - ->CurrentPersistencyManager()); - } - - /** - * Builds the output ROOT event. - * - * @param anEvent The Geant4 event data. - */ - G4bool Store(const G4Event *anEvent); - - /** - * This gets called automatically at the end of the run and is used to write - * out the run header and close the writer. - * - * @param aRun The Geant4 run data (not used right now) - * - * @return True if event is stored (function is hard-coded to return true). - */ - G4bool Store(const G4Run *aRun); - - /** - * Implementing this makes an "overloaded-virtual" compiler warning go away. - */ - G4bool Store(const G4VPhysicalVolume *) { return false; } - - /** - * This is called "manually" in UserRunAction to open the ROOT writer for the - * run. - */ - void Initialize(); - - /** - * Set the current ldmx-sw event. This is used by the persistency - * manager to retrieve and fill the containers that will be - * persisted. - * - * @param event Event buffer for the current event. - */ - void setCurrentEvent(framework::Event *event) { event_ = event; } - - /** - * Set the number of events began and completed. - * - * These two numbers may or may not be equal - * depending on if the simulation ran with any filters - * that would abort events early. - * - * These numbers are helpful for evaluating filtering - * performance, so we put them both in the RunHeader. - * - * @param[in] began number of events that were started - * @param[in] completed number of events that were completed without an abort - * signal - */ - void setNumEvents(int began, int completed) { - eventsBegan_ = began; - eventsCompleted_ = completed; - } - - public: - /** - * Build an output event from the current Geant4 event. - * - * @param anEvent The Geant4 event. - * @param outputEvent The output event. - */ - void buildEvent(const G4Event *anEvent); - - /** - * Write header info into the output event from Geant4. - * - * @param anEvent The Geant4 event. - */ - void writeHeader(const G4Event *anEvent); - - /** - * Write hits collections from Geant4 into a ROOT event. - * - * @param anEvent The Geant4 event. - * @param outputEvent The output event. - */ - void writeHitsCollections(const G4Event *anEvent, - framework::Event *outputEvent); - - /** - * Write a collection of tracker hits to an output collection. - * - * @param hc The collection of G4TrackerHits. - * @param outputColl The output collection of SimTrackerHits. - */ - void writeTrackerHitsCollection(G4TrackerHitsCollection *hc, - std::vector &outputColl); - - /** - * Write a collection of tracker hits to an output collection. - * - * @param hc The collection of G4CalorimeterHits. - * @param outputColl The output collection of SimCalorimeterHits. - */ - void writeCalorimeterHitsCollection( - G4CalorimeterHitsCollection *hc, - std::vector &outputColl); - - private: - /// Configuration parameters passed to Simulator - framework::config::Parameters parameters_; - - /// Run Number, given to us by Simulator from Process - int run_; - - /// Number of events started on this production run - int eventsBegan_{-1}; - - /// Number of events completed without being aborted (due to filters) - int eventsCompleted_{-1}; - - /// The output file. - framework::EventFile &file_; - - /// The event container used to manage the tree/branches/collections. - framework::Event *event_{nullptr}; - - /// Handles ECal hit readout and IO. - EcalHitIO ecalHitIO_; -}; - -} // namespace persist -} // namespace g4fire - -#endif diff --git a/include/g4fire/RootCompleteReSim.h b/include/g4fire/RootCompleteReSim.h deleted file mode 100644 index de73026..0000000 --- a/include/g4fire/RootCompleteReSim.h +++ /dev/null @@ -1,95 +0,0 @@ -/** - * @file RootCompleteReSim.h - * @brief Primary generator used to generate primaries from SimParticles. - * @author Nhan Tran, Fermilab - * @author Omar Moreno, SLAC National Accelerator Laboratory - * @author Tom Eichlersmith, University of Minnesota - */ - -#ifndef SIMCORE_ROOTCOMPLETERESIM_H -#define SIMCORE_ROOTCOMPLETERESIM_H - -//----------------// -// C++ StdLib // -//----------------// -#include -#include -#include - -//----------// -// ROOT // -//----------// -#include "TFile.h" -#include "TTree.h" -#include "TVector3.h" - -//-------------// -// ldmx-sw // -//-------------// -#include "Framework/Event.h" -#include "Framework/EventFile.h" -#include "g4fire/PrimaryGenerator.h" - -class G4Event; - -namespace g4fire { - -class Parameters; - -/** - * @class RootCompleteReSim - * - * PrimaryGenerator that gets primaries and event seeds and - * inputs them into current event as primaries with the exact same kinematics. - */ -class RootCompleteReSim : public PrimaryGenerator { - public: - /** - * Class constructor. - * @param name the name of the generator - * @param parameters configuration parameters - * - * Parameters: - * filePath : path to root file to re-sim - * simParticleCollName : name of collection of SimParticles - * simParticlePassName : name of pass of SimParticles - */ - RootCompleteReSim(const std::string& name, - framework::config::Parameters& parameters); - - /** - * Class destructor. - */ - virtual ~RootCompleteReSim(); - - /** - * Generate vertices in the Geant4 event. - * @param anEvent The Geant4 event. - */ - void GeneratePrimaryVertex(G4Event* anEvent); - - private: - /** - * Name of SimParticles collection - */ - std::string simParticleCollName_; - - /** - * Name of SimParticles pass - */ - std::string simParticlePassName_; - - /** - * The input root file - */ - std::unique_ptr ifile_; - - /** - * The input ldmx event bus - */ - framework::Event ievent_; -}; - -} // namespace g4fire - -#endif // SIMCORE_ROOTCOMPLETERESIM_H diff --git a/include/g4fire/RootSimFromEcalSP.h b/include/g4fire/RootSimFromEcalSP.h deleted file mode 100644 index acadedd..0000000 --- a/include/g4fire/RootSimFromEcalSP.h +++ /dev/null @@ -1,106 +0,0 @@ -/** - * @file RootSimFromEcalSP.h - * @brief Primary generator used to generate primaries from SimParticles. - * @author Nhan Tran, Fermilab - * @author Omar Moreno, SLAC National Accelerator Laboratory - * @author Tom Eichlersmith, University of Minnesota - */ - -#ifndef SIMCORE_ROOTSIMFROMECALSP_H -#define SIMCORE_ROOTSIMFROMECALSP_H - -//----------------// -// C++ StdLib // -//----------------// -#include -#include -#include - -//----------// -// ROOT // -//----------// -#include "TFile.h" -#include "TTree.h" -#include "TVector3.h" - -//-------------// -// ldmx-sw // -//-------------// -#include "Framework/Event.h" -#include "Framework/EventFile.h" -#include "g4fire/PrimaryGenerator.h" - -class G4Event; - -namespace g4fire { - -class Parameters; - -/** - * @class RootSimFromEcalSP - * - * Generate primaries that correspond to particles leaving - * the ECal (passing through the Ecal Scoring Planes) before - * a time cutoff. - */ -class RootSimFromEcalSP : public PrimaryGenerator { - public: - /** - * Class constructor. - * @param name name of this generator - * @param parameters configuration parameters - * - * Parameters: - * filePath : path to root file containing events to re-sim - * timeCutoff : maximum time that a particle can pass through Ecal Scoring - * Planes and be included in re-sim (ns) ecalSPHitsCollName : name of - * collection for Ecal Scoring Plane hits ecalSPHitsPassName : name of pass - * for Ecal Scoring Plane hits - */ - RootSimFromEcalSP(const std::string& name, - framework::config::Parameters& parameters); - - /** - * Class destructor. - */ - virtual ~RootSimFromEcalSP(); - - /** - * Generate vertices in the Geant4 event. - * @param anEvent The Geant4 event. - */ - void GeneratePrimaryVertex(G4Event* anEvent); - - private: - /** - * The cutoff time - * - * Any particle passing through the ECal scoring planes at a time - * greater than this time is NOT re-used as a primary. - */ - double timeCutoff_; - - /** - * The Ecal Scoring Planes Hits collection name - */ - std::string ecalSPHitsCollName_; - - /** - * The Ecal Scoring Planes Hits pass name - */ - std::string ecalSPHitsPassName_; - - /** - * The input root file - */ - std::unique_ptr ifile_; - - /** - * The input ldmx event bus - */ - framework::Event ievent_; -}; - -} // namespace g4fire - -#endif // SIMCORE_ROOTSIMFROMECALSP_H diff --git a/include/g4fire/ScoringPlaneSD.h b/include/g4fire/ScoringPlaneSD.h deleted file mode 100644 index 87d9309..0000000 --- a/include/g4fire/ScoringPlaneSD.h +++ /dev/null @@ -1,74 +0,0 @@ -#ifndef SIMCORE_SCORINGPLANESD_H -#define SIMCORE_SCORINGPLANESD_H - -/*~~~~~~~~~~~~*/ -/* Geant4 */ -/*~~~~~~~~~~~~*/ -#include "G4VSensitiveDetector.hh" - -/*~~~~~~~~~~~~~~*/ -/* DetDescr */ -/*~~~~~~~~~~~~~~*/ -#include "DetDescr/DetectorID.h" - -/*~~~~~~~~~~~~~~~~~~~~*/ -/* g4fire */ -/*~~~~~~~~~~~~~~~~~~~~*/ -#include "g4fire/G4TrackerHit.h" - -// Forward declaration -class G4Step; - -namespace g4fire { - -/** - * Class defining a basic sensitive detector for scoring planes. - */ -class ScoringPlaneSD : public G4VSensitiveDetector { - public: - /** - * Constructor - * - * @param name The name of the sensitive detector. - * @param collectionID The name of the hits collection. - * @param subDetID The subdetector ID. - */ - ScoringPlaneSD(G4String name, G4String colName, int subDetID); - - /** Destructor */ - ~ScoringPlaneSD(); - - /** - * Process a step and create a hit out of it. - * - * @param step The current step. - * @param history The readout history. - */ - G4bool ProcessHits(G4Step* step, G4TouchableHistory* history); - - /** - * Initialize the sensitive detector. - * - * @param hcEvent The hits collections associated with this event. - */ - void Initialize(G4HCofThisEvent* hcEvent); - - /** - * End of event hook. - * - * @param hcEvent The hits collections associated with this event. - */ - void EndOfEvent(G4HCofThisEvent* hcEvent); - - private: - /** Output hits collection */ - G4TrackerHitsCollection* hitsCollection_{nullptr}; - - /** The detector ID. */ - // DetectorID* detID_{new DefaultDetectorID()}; - -}; // ScoringPlaneSD - -} // namespace g4fire - -#endif // SIMCORE_SCORINGPLANESD_H_ diff --git a/include/g4fire/Simulator.h b/include/g4fire/Simulator.h index f2a5e30..d3298b3 100644 --- a/include/g4fire/Simulator.h +++ b/include/g4fire/Simulator.h @@ -1,5 +1,4 @@ -#ifndef G4FIRE_SIMULATOR_H -#define G4FIRE_SIMULATOR_H +#pragma once #include #include @@ -8,69 +7,70 @@ #include "fire/Processor.h" #include "fire/config/Parameters.h" -#include "fire/config/Parameters.h" -//#include "Framework/EventDef.h" -#include "fire/Processor.h" #include "g4fire/ConditionsInterface.h" +#include "g4fire/event/EventBuilder.h" class G4UImanager; class G4UIsession; -class G4GDMLParser; -class G4GDMLMessenger; -class G4CascadeParameters; namespace g4fire { class RunManager; -class EventFile; -class ParameterSet; -class DetectorConstruction; /** * Geant4 simulation wrapped within a fire producer. */ class Simulator : public fire::Processor { -public: + public: + /** + * Constructor. + * + * @param params The parameters used to configure the simulation. + */ Simulator(const fire::config::Parameters ¶ms); /// Default destructor ~Simulator() = default; + /** + * Generate an event by firing the particles of interest through the detector + * using the generator specified in the configuration. + * + * This method is called once per event and initializes the propagation of + * particles through the detector. All Geant4 tracks (particles) generated in + * the propagation are stored along with basic information (e.g. momentum, + * start position). Once all tracks have been fully propagated, the event is + * passed to the EventBuilder which creates objects to persist and adds them + * to the event bus. + * + * @param[in] event The current event being processed. + */ void process(fire::Event &event) final override; /** - * Given a non-const reference to the new RunHeader, - * we can add parameters from the simulation here - * before the run starts. + * Given a non-const reference to the new RunHeader, we can add parameters + * from the simulation here before the run starts. * - * @param header of new run + * @param[in] header The header of the new run. */ void beforeNewRun(fire::RunHeader &header) final override; /** - * Before the run starts (but after the conditions are configured) - * set up the random seeds for this run. + * Before the run starts (but after the conditions are configured) set up the + * random seeds for this run. * * @param[in] header RunHeader for this run, unused */ void onNewRun(const fire::RunHeader &header) final override; /** - * Callback for the EventProcessor to take any necessary action - * when a new file is opened. + * Callback for the Processor to take any necessary action when a event input + * file is closed. * - * @param eventFile The input/output file. + * @param[in] file_name The input event file name. */ - //void onFileOpen(fire::EventFile &eventFile) final override; - - /** - * Callback for the EventProcessor to take any necessary action - * when a file is closed. - * - * @param eventFile The intput/output file. - */ - // void onFileClose(fire::EventFile &eventFile) final override; + void onFileClose(const std::string &file_name) final override; /** * Initialization of simulation @@ -85,7 +85,7 @@ class Simulator : public fire::Processor { /// Callback called once processing is complete. void onProcessEnd() final override; -private: + private: /** * Configure the simulation given the set of parameters passed by the user * at runtime. @@ -124,15 +124,14 @@ class Simulator : public fire::Processor { /// User interface handle G4UImanager *ui_manager_{nullptr}; - /// PersistencyManager - // std::unique_ptr - // persistencyManager_; + /// The event builder used to create the objects to persist. + g4fire::event::EventBuilder eb_; /// Handle to the G4Session -> how to deal with G4cout and G4cerr std::unique_ptr session_handle_; - /// Commands not allowed to be passed from python config file - /// This is because Simulator already runs them. + /// Commands not allowed to be passed from python config file. This is + /// because Simulator already runs them. static const std::vector invalid_cmds; /// Number of events started @@ -144,10 +143,6 @@ class Simulator : public fire::Processor { /// Conditions interface ConditionsInterface conditions_intf_; - /********************************************************* - * Python Configuration Parameters - *********************************************************/ - /// The parameters used to configure the simulation fire::config::Parameters params_; @@ -156,4 +151,3 @@ class Simulator : public fire::Processor { }; // Simulator } // namespace g4fire -#endif // G4FIRE_SIMULATOR_H diff --git a/include/g4fire/TrackMap.h b/include/g4fire/TrackMap.h index 8252c33..364dd42 100644 --- a/include/g4fire/TrackMap.h +++ b/include/g4fire/TrackMap.h @@ -2,12 +2,11 @@ #include -#include "G4Event.hh" -#include "G4Track.hh" - -//#include "g4fire/Event/SimParticle.h" -#include "g4fire/UserTrackInformation.h" #include "g4fire/UserPrimaryParticleInformation.h" +#include "g4fire/UserTrackInformation.h" +#include "g4fire/event/SimParticle.h" + +class G4Track; namespace g4fire { @@ -21,19 +20,21 @@ namespace g4fire { * parent and children faithfully recorded in the output file. */ class TrackMap { - public: +public: /** * Add a record in the map for the input track. * - * @param track G4Track to insert + * @param track G4Track to insert. */ - void insert(const G4Track* track); + void insert(const G4Track *track); /** * Check if the passed track has already been inserted * into the track map. + * + * @param track G4Track to check for. */ - inline bool contains(const G4Track* track) const { + inline bool contains(const G4Track *track) const { return ancestry_.find(track->GetTrackID()) != ancestry_.end(); } @@ -49,15 +50,14 @@ class TrackMap { int findIncident(int track_id) const; /** - * Return true if the given track ID is saved - * i.e. will be stored in output file + * Return true if the given track ID will be persisted. * * @param track_id The track ID. * @return True if the track ID has been inserted in output particle map */ - //inline bool isSaved(int track_id) const { - // return particle_map_.find(track_id) != particle_map_.end(); - //} + inline bool isSaved(int track_id) const { + return particle_map_.find(track_id) != particle_map_.end(); + } /** * Add a track to be stored into output map @@ -66,7 +66,7 @@ class TrackMap { * kinematics. * @param track G4Track to store into output */ - void save(const G4Track* track); + void save(const G4Track *track); /** * Trace the ancestry for the particles that will be stored. @@ -86,23 +86,23 @@ class TrackMap { void clear(); /** - * Get the map of particles to be stored in output event. + * @return The map of particles to be persisted. */ - //std::map &getParticleMap() { - // return particle_map_; - //} + std::map &particleMap() { + return particle_map_; + } - private: +private: /** * Was the input track generated inside the calorimeter region? * * We rely on the fact that the calorimeter region is named - * 'CalorimeterRegion' - * and no other region names contain the string 'Calorimeter' + * 'CalorimeterRegion' and no other region names contain the string + * 'Calorimeter'. */ - bool isInCalorimeterRegion(const G4Track* track) const; + bool isInCalorimeterRegion(const G4Track *track) const; - private: +private: /** * ancestry map of particles in event (child -> parent) * @@ -121,12 +121,12 @@ class TrackMap { * @see isInCalorimeterRegion for how we check if a track * originated in the calorimeter region. */ - std::unordered_map> ancestry_; + std::unordered_map> ancestry_; /// descendents map of particles in event (parent -> children) - std::unordered_map> descendents_; + std::unordered_map> descendents_; /// map of SimParticles that will be stored - //std::map particle_map_; + std::map particle_map_; }; -} // namespace g4fire +} // namespace g4fire diff --git a/include/g4fire/TrackerSD.h b/include/g4fire/TrackerSD.h deleted file mode 100644 index 1ea1d8a..0000000 --- a/include/g4fire/TrackerSD.h +++ /dev/null @@ -1,72 +0,0 @@ -#ifndef SIMCORE_TRACKERSD_H -#define SIMCORE_TRACKERSD_H - -/*~~~~~~~~~~~~*/ -/* Geant4 */ -/*~~~~~~~~~~~~*/ -#include "G4VSensitiveDetector.hh" - -/*~~~~~~~~~~~~~~*/ -/* DetDescr */ -/*~~~~~~~~~~~~~~*/ -#include "DetDescr/TrackerID.h" - -/*~~~~~~~~~~~~~~~~~~~~*/ -/* g4fire */ -/*~~~~~~~~~~~~~~~~~~~~*/ -#include "g4fire/G4TrackerHit.h" - -namespace g4fire { - -/** - * @class TrackerSD - * @brief Basic sensitive detector for trackers - * - * @note - * This class creates a G4TrackerHit for each step within the subdetector. - */ -class TrackerSD : public G4VSensitiveDetector { - public: - /** - * Class constructor. - * @param[in] name The name of the sensitive detector. - * @param[in] collectionName The name of the hits collection. - * @param[in] subDetID The subdetector ID. - */ - TrackerSD(G4String name, G4String collectionName, int subDetID); - - /// Destructor - ~TrackerSD(); - - /** - * Process a step by creating a hit. - * - * @param step The step information - * @param history The readout history. - */ - G4bool ProcessHits(G4Step* step, G4TouchableHistory* history); - - /** - * Initialize the sensitive detector. - * @param hcEvent The hits collections of the event. - */ - void Initialize(G4HCofThisEvent* hcEvent); - - /** - * End of event hook. - * @param hcEvent The hits collections of the event. - */ - void EndOfEvent(G4HCofThisEvent* hcEvent); - - private: - /// The output hits collection of G4TrackerHits. - G4TrackerHitsCollection* hitsCollection_; - - /// The detector ID - ldmx::SubdetectorIDType subDetID_; - -}; // TrackerID - -} // namespace g4fire - -#endif // SIMCORE_TRACKERSD_H diff --git a/include/g4fire/TrigScintSD.h b/include/g4fire/TrigScintSD.h deleted file mode 100644 index 40550be..0000000 --- a/include/g4fire/TrigScintSD.h +++ /dev/null @@ -1,46 +0,0 @@ -#ifndef SIMCORE_TRIGSD_H -#define SIMCORE_TRIGSD_H - -/*~~~~~~~~~~~~~~~~~~~~*/ -/* g4fire */ -/*~~~~~~~~~~~~~~~~~~~~*/ -#include "g4fire/CalorimeterSD.h" - -// Forward declarations -class G4Step; -class G4TouchableHistory; - -namespace g4fire { - -/** - * Class defining a sensitive detector of type trigger scintillator. - */ -class TrigScintSD : public CalorimeterSD { - public: - /** - * Class constructor. - * - * @param[in] name The name of the sensitive detector. - * @param[in] theCollectionName The name of the hits collection. - * @param[in] subdet The subdetector ID. - */ - TrigScintSD(G4String name, G4String collectionName, int subDetID); - - /// Destructor - ~TrigScintSD(); - - /** - * Process steps to create hits. - * - * @param[in] step The step information. - * @param[in] history The readout history. - */ - G4bool ProcessHits(G4Step* step, G4TouchableHistory* history); - - private: - int moduleId_; -}; - -} // namespace g4fire - -#endif diff --git a/include/g4fire/UserEventInformation.h b/include/g4fire/UserEventInformation.h index 2eadd8d..ca82577 100644 --- a/include/g4fire/UserEventInformation.h +++ b/include/g4fire/UserEventInformation.h @@ -1,14 +1,15 @@ -#ifndef G4FIRE_USEREVENTINFORMATION_H -#define G4FIRE_USEREVENTINFORMATION_H +#pragma once #include "G4VUserEventInformation.hh" +#include "g4fire/Parameters.h" + namespace g4fire { /** * Encapsulates user defined information associated with a Geant4 event. */ -class UserEventInformation : public G4VUserEventInformation { +class UserEventInformation : public G4VUserEventInformation, Parameters { public: /// Constructor UserEventInformation() = default; @@ -17,7 +18,6 @@ class UserEventInformation : public G4VUserEventInformation { ~UserEventInformation() = default; /// Print the information associated with the track - // TODO(OM): Use stream operator instead void Print() const final override; /// Increment the number of brem candidates in an event. @@ -36,7 +36,7 @@ class UserEventInformation : public G4VUserEventInformation { /** * @return The event weight */ - double getWeight() { return weight_; } + double getWeight() { return weight_; } /** * Increment the event weight by the input weight @@ -144,5 +144,3 @@ class UserEventInformation : public G4VUserEventInformation { bool last_step_en_{false}; }; } // namespace g4fire - -#endif // G4FIRE_USEREVENTINFORMATION_H diff --git a/include/g4fire/UserRegionInformation.h b/include/g4fire/UserRegionInformation.h index 17928ce..d0f1e22 100644 --- a/include/g4fire/UserRegionInformation.h +++ b/include/g4fire/UserRegionInformation.h @@ -1,5 +1,4 @@ -#ifndef G4FIRE_USERREGIONINFORMATION_H_ -#define G4FIRE_USERREGIONINFORMATION_H_ +#pragma once #include "G4VUserRegionInformation.hh" @@ -18,7 +17,7 @@ class UserRegionInformation : public G4VUserRegionInformation { public: UserRegionInformation(bool store_secondaries); - virtual ~UserRegionInformation(); + ~UserRegionInformation() = default; void Print() const; @@ -29,5 +28,3 @@ class UserRegionInformation : public G4VUserRegionInformation { }; } // namespace g4fire - -#endif diff --git a/include/g4fire/XsecBiasingOperator.h b/include/g4fire/XsecBiasingOperator.h index d8a8db9..086017f 100644 --- a/include/g4fire/XsecBiasingOperator.h +++ b/include/g4fire/XsecBiasingOperator.h @@ -1,8 +1,7 @@ -#ifndef G4FIRE_XSECBIASINGOPERATOR_H_ -#define G4FIRE_XSECBIASINGOPERATOR_H_ +#pragma once #include "fire/config/Parameters.h" -//#include "Framework/RunHeader.h" +#include "fire/RunHeader.h" #include "G4BOptnChangeCrossSection.hh" #include "G4BiasingProcessInterface.hh" @@ -126,7 +125,7 @@ class XsecBiasingOperator : public G4VBiasingOperator { * * @param[in,out] header RunHeader to write configuration to */ - //virtual void RecordConfig(ldmx::RunHeader& header) const = 0; + virtual void RecordConfig(fire::RunHeader& header) const = 0; protected: /** @@ -193,5 +192,3 @@ class XsecBiasingOperator : public G4VBiasingOperator { g4fire::XsecBiasingOperator::declare( \ std::string(#NS) + "::" + std::string(#CLASS), &CLASS##Builder); \ } - -#endif // G4FIRE_XSECBIASINGOPERATOR_H_ diff --git a/include/g4fire/BiasOperators/DarkBrem.h b/include/g4fire/biasing/DarkBrem.h similarity index 79% rename from include/g4fire/BiasOperators/DarkBrem.h rename to include/g4fire/biasing/DarkBrem.h index be4eaa0..8b21837 100644 --- a/include/g4fire/BiasOperators/DarkBrem.h +++ b/include/g4fire/biasing/DarkBrem.h @@ -1,34 +1,29 @@ -#ifndef SIMCORE_BIASOPERATORS_DARKBREM_H -#define SIMCORE_BIASOPERATORS_DARKBREM_H +#pragma once -//----------// -// LDMX // -//----------// -#include "g4fire/DarkBrem/G4eDarkBremsstrahlung.h" #include "g4fire/XsecBiasingOperator.h" +#include "g4fire/darkbrem/G4eDarkBremsstrahlung.h" -namespace g4fire { -namespace biasoperators { +namespace g4fire::biasing { /** * Bias operator for the dark brem process */ class DarkBrem : public XsecBiasingOperator { - public: +public: /** * Constructor * * Calls base class constructor and allows * access to configuration parameters. */ - DarkBrem(std::string name, const framework::config::Parameters& p); + DarkBrem(std::string name, const fire::config::Parameters &p); /** * Destructor * * Blank right now */ - ~DarkBrem() {} + ~DarkBrem() = default; /** * Calculate the biased cross section given the @@ -43,9 +38,9 @@ class DarkBrem : public XsecBiasingOperator { * @return Method that returns the biasing operation that will be used * to bias the occurence of events. */ - G4VBiasingOperation* ProposeOccurenceBiasingOperation( - const G4Track* track, - const G4BiasingProcessInterface* callingProcess) final override; + G4VBiasingOperation *ProposeOccurenceBiasingOperation( + const G4Track *track, + const G4BiasingProcessInterface *callingProcess) final override; /// Return the name of the process this operator biases virtual std::string getProcessToBias() const { @@ -63,9 +58,9 @@ class DarkBrem : public XsecBiasingOperator { * * @param[in,out] header RunHeader to record configuration to */ - virtual void RecordConfig(ldmx::RunHeader& header) const; + virtual void RecordConfig(fire::RunHeader &header) const; - protected: +protected: /** * DEBUG FUNCTION * This function is called by the biasing interface class during PostStepDoIt. @@ -92,7 +87,7 @@ class DarkBrem : public XsecBiasingOperator { } */ - private: +private: /// volume we want to bias in std::string volume_; @@ -102,8 +97,5 @@ class DarkBrem : public XsecBiasingOperator { /// should we bias all electrons? (or only the primary) bool bias_all_; -}; // DarkBrem -} // namespace biasoperators -} // namespace g4fire - -#endif // SIMCORE_DARKBREM_DARKBREMXSECBIASINGOPERATOR_H +}; // DarkBrem +} // namespace g4fire::biasing diff --git a/include/g4fire/BiasOperators/ElectroNuclear.h b/include/g4fire/biasing/ElectroNuclear.h similarity index 67% rename from include/g4fire/BiasOperators/ElectroNuclear.h rename to include/g4fire/biasing/ElectroNuclear.h index 1759d09..4ff2c0a 100644 --- a/include/g4fire/BiasOperators/ElectroNuclear.h +++ b/include/g4fire/biasing/ElectroNuclear.h @@ -1,10 +1,8 @@ -#ifndef SIMCORE_BIASOPERATORS_ELECTRONUCLEAR_H_ -#define SIMCORE_BIASOPERATORS_ELECTRONUCLEAR_H_ +#pragma once #include "g4fire/XsecBiasingOperator.h" -namespace g4fire { -namespace biasoperators { +namespace g4fire::biasing { /** * Bias the Electron-Nuclear process @@ -17,10 +15,10 @@ class ElectroNuclear : public XsecBiasingOperator { * Calls parent constructor and allows * accesss to configuration parameters. */ - ElectroNuclear(std::string name, const framework::config::Parameters& p); + ElectroNuclear(std::string name, const fire::config::Parameters& p); /** Destructor */ - ~ElectroNuclear() {} + ~ElectroNuclear() = default; /** * @return Method that returns the biasing operation that will be used @@ -44,10 +42,10 @@ class ElectroNuclear : public XsecBiasingOperator { * * @param[in,out] header RunHeader to record to */ - virtual void RecordConfig(ldmx::RunHeader& header) const { - header.setStringParameter("BiasOperator::ElectroNuclear::Volume", volume_); - header.setFloatParameter("BiasOperator::ElectroNuclear::Factor", factor_); - header.setFloatParameter("BiasOperator::ElectroNuclear::Threshold", + virtual void RecordConfig(fire::RunHeader& header) const { + header.set("biasing::ElectroNuclear::Volume", volume_); + header.set("biasing::ElectroNuclear::Factor", factor_); + header.set("biasing::ElectroNuclear::Threshold", threshold_); } @@ -62,8 +60,4 @@ class ElectroNuclear : public XsecBiasingOperator { double threshold_; }; // ElectroNuclear - -} // namespace biasoperators } // namespace g4fire - -#endif // SIMCORE_BIASOPERATORS_ELECTRONUCLEAR_H_ diff --git a/include/g4fire/BiasOperators/GammaToMuPair.h b/include/g4fire/biasing/GammaToMuPair.h similarity index 56% rename from include/g4fire/BiasOperators/GammaToMuPair.h rename to include/g4fire/biasing/GammaToMuPair.h index 25bacec..eb6ac58 100644 --- a/include/g4fire/BiasOperators/GammaToMuPair.h +++ b/include/g4fire/biasing/GammaToMuPair.h @@ -1,23 +1,19 @@ -#ifndef SIMCORE_BIASOPERATORS_GAMMATOMUPAIR_H_ -#define SIMCORE_BIASOPERATORS_GAMMATOMUPAIR_H_ - #include "g4fire/XsecBiasingOperator.h" -namespace g4fire { -namespace biasoperators { +namespace g4fire::biasing { /** * Bias the Gamma to Mu Pair process */ class GammaToMuPair : public XsecBiasingOperator { - public: +public: /** * Constructor * * Calls parent constructor and allows * accesss to configuration parameters. */ - GammaToMuPair(std::string name, const framework::config::Parameters& p); + GammaToMuPair(std::string name, const fire::config::Parameters &p); /** Destructor */ ~GammaToMuPair() = default; @@ -26,9 +22,9 @@ class GammaToMuPair : public XsecBiasingOperator { * @return Method that returns the biasing operation that will be used * to bias the conversion of gammas to muon pairs. */ - G4VBiasingOperation* ProposeOccurenceBiasingOperation( - const G4Track* track, - const G4BiasingProcessInterface* callingProcess) final override; + G4VBiasingOperation *ProposeOccurenceBiasingOperation( + const G4Track *track, + const G4BiasingProcessInterface *callingProcess) final override; /// Return the process to bias virtual std::string getProcessToBias() const { return "GammaToMuPair"; } @@ -44,14 +40,13 @@ class GammaToMuPair : public XsecBiasingOperator { * * @param[in,out] header RunHeader to record to */ - virtual void RecordConfig(ldmx::RunHeader& header) const { - header.setStringParameter("BiasOperator::GammaToMuPair::Volume", volume_); - header.setFloatParameter("BiasOperator::GammaToMuPair::Factor", factor_); - header.setFloatParameter("BiasOperator::GammaToMuPair::Threshold", - threshold_); + virtual void RecordConfig(fire::RunHeader &header) const { + header.set("biasing::GammaToMuPair::Volume", volume_); + header.set("biasing::GammaToMuPair::Factor", factor_); + header.set("biasing::GammaToMuPair::Threshold", threshold_); } - private: +private: /// The volume to bias in std::string volume_; @@ -61,8 +56,4 @@ class GammaToMuPair : public XsecBiasingOperator { /// Minimum kinetic energy [MeV] to allow a track to be biased double threshold_; }; - -} // namespace biasoperators -} // namespace g4fire - -#endif // SIMCORE_BIASOPERATORS_GAMMATOMUPAIR_H_ +} // namespace g4fire::biasing diff --git a/include/g4fire/BiasOperators/K0LongInelastic.h b/include/g4fire/biasing/K0LongInelastic.h similarity index 55% rename from include/g4fire/BiasOperators/K0LongInelastic.h rename to include/g4fire/biasing/K0LongInelastic.h index c964e37..68aa1b3 100644 --- a/include/g4fire/BiasOperators/K0LongInelastic.h +++ b/include/g4fire/biasing/K0LongInelastic.h @@ -1,23 +1,19 @@ -#ifndef SIMCORE_BIASOPERATORS_K0LONGINELASTIC_H_ -#define SIMCORE_BIASOPERATORS_K0LONGINELASTIC_H_ - #include "g4fire/XsecBiasingOperator.h" -namespace g4fire { -namespace biasoperators { +namespace g4fire::biasing { /** * Bias the k0 long inelastic collisions */ class K0LongInelastic : public XsecBiasingOperator { - public: +public: /** * Constructor * * Calls parent constructor and allows * accesss to configuration parameters. */ - K0LongInelastic(std::string name, const framework::config::Parameters& p); + K0LongInelastic(std::string name, const fire::config::Parameters &p); /** Destructor */ ~K0LongInelastic() = default; @@ -26,9 +22,9 @@ class K0LongInelastic : public XsecBiasingOperator { * @return Method that returns the biasing operation that will be used * to bias the K0 Long inelactic hadronic interactions. */ - G4VBiasingOperation* ProposeOccurenceBiasingOperation( - const G4Track* track, - const G4BiasingProcessInterface* callingProcess) final override; + G4VBiasingOperation *ProposeOccurenceBiasingOperation( + const G4Track *track, + const G4BiasingProcessInterface *callingProcess) final override; /// Return the process to bias virtual std::string getProcessToBias() const { return "kaon0LInelastic"; } @@ -44,26 +40,20 @@ class K0LongInelastic : public XsecBiasingOperator { * * @param[in,out] header RunHeader to record to */ - virtual void RecordConfig(ldmx::RunHeader& header) const { - header.setStringParameter("BiasOperator::K0LongInelastic::Volume", volume_); - header.setFloatParameter("BiasOperator::K0LongInelastic::Factor", factor_); - header.setFloatParameter("BiasOperator::K0LongInelastic::Threshold", - threshold_); + virtual void RecordConfig(fire::RunHeader &header) const { + header.set("biasing::K0LongInelastic::Volume", volume_); + header.set("biasing::K0LongInelastic::Factor", factor_); + header.set("biasing::K0LongInelastic::Threshold", threshold_); } - private: +private: /// The volume to bias in std::string volume_; - /// The biasing factor + /// The bias factor double factor_; /// Minimum kinetic energy [MeV] to allow a track to be biased double threshold_; - }; - -} // namespace biasoperators -} // namespace g4fire - -#endif // SIMCORE_BIASOPERATORS_K0LONGINELASTIC_H_ +} // namespace g4fire::biasing diff --git a/include/g4fire/BiasOperators/NeutronInelastic.h b/include/g4fire/biasing/NeutronInelastic.h similarity index 56% rename from include/g4fire/BiasOperators/NeutronInelastic.h rename to include/g4fire/biasing/NeutronInelastic.h index 7309c32..c76ad48 100644 --- a/include/g4fire/BiasOperators/NeutronInelastic.h +++ b/include/g4fire/biasing/NeutronInelastic.h @@ -1,23 +1,19 @@ -#ifndef SIMCORE_BIASOPERATORS_NEUTRONINELASTIC_H_ -#define SIMCORE_BIASOPERATORS_NEUTRONINELASTIC_H_ - #include "g4fire/XsecBiasingOperator.h" -namespace g4fire { -namespace biasoperators { +namespace g4fire::biasing { /** * Bias the neutron inelastic collsions */ class NeutronInelastic : public XsecBiasingOperator { - public: +public: /** * Constructor * * Calls parent constructor and allows * accesss to configuration parameters. */ - NeutronInelastic(std::string name, const framework::config::Parameters& p); + NeutronInelastic(std::string name, const fire::config::Parameters &p); /** Destructor */ ~NeutronInelastic() = default; @@ -26,9 +22,9 @@ class NeutronInelastic : public XsecBiasingOperator { * @return Method that returns the biasing operation that will be used * to bias the neutron inelactic hadronic interactions. */ - G4VBiasingOperation* ProposeOccurenceBiasingOperation( - const G4Track* track, - const G4BiasingProcessInterface* callingProcess) final override; + G4VBiasingOperation *ProposeOccurenceBiasingOperation( + const G4Track *track, + const G4BiasingProcessInterface *callingProcess) final override; /// Return the process to bias virtual std::string getProcessToBias() const { return "neutronInelastic"; } @@ -44,14 +40,13 @@ class NeutronInelastic : public XsecBiasingOperator { * * @param[in,out] header RunHeader to record to */ - virtual void RecordConfig(ldmx::RunHeader& header) const { - header.setStringParameter("BiasOperator::NeutronInelastic::Volume", volume_); - header.setFloatParameter("BiasOperator::NeutronInelastic::Factor", factor_); - header.setFloatParameter("BiasOperator::NeutronInelastic::Threshold", - threshold_); + virtual void RecordConfig(fire::RunHeader &header) const { + header.set("biasing::NeutronInelastic::Volume", volume_); + header.set("biasing::NeutronInelastic::Factor", factor_); + header.set("biasing::NeutronInelastic::Threshold", threshold_); } - private: +private: /// The volume to bias in std::string volume_; @@ -60,10 +55,5 @@ class NeutronInelastic : public XsecBiasingOperator { /// Minimum kinetic energy [MeV] to allow a track to be biased double threshold_; - }; - -} // namespace biasoperators -} // namespace g4fire - -#endif // SIMCORE_BIASOPERATORS_NEUTRONINELASTIC_H_ +} // namespace g4fire::biasing diff --git a/include/g4fire/BiasOperators/PhotoNuclear.h b/include/g4fire/biasing/PhotoNuclear.h similarity index 57% rename from include/g4fire/BiasOperators/PhotoNuclear.h rename to include/g4fire/biasing/PhotoNuclear.h index 030d50a..ee97424 100644 --- a/include/g4fire/BiasOperators/PhotoNuclear.h +++ b/include/g4fire/biasing/PhotoNuclear.h @@ -1,19 +1,15 @@ -#ifndef SIMCORE_BIASOPERATORS_PHOTONUCLEAR_H_ -#define SIMCORE_BIASOPERATORS_PHOTONUCLEAR_H_ - #include "g4fire/XsecBiasingOperator.h" -namespace g4fire { -namespace biasoperators { +namespace g4fire::biasing { /** * Bias the Photon-Nuclear process */ class PhotoNuclear : public XsecBiasingOperator { - public: +public: /** Constructor */ - PhotoNuclear(std::string name, const framework::config::Parameters& p); + PhotoNuclear(std::string name, const fire::config::Parameters &p); /** Method called at the beginning of a run. */ void StartRun(); @@ -22,9 +18,9 @@ class PhotoNuclear : public XsecBiasingOperator { * @return Method that returns the biasing operation that will be used * to bias the occurence of photonuclear events. */ - G4VBiasingOperation* ProposeOccurenceBiasingOperation( - const G4Track* track, - const G4BiasingProcessInterface* callingProcess) final override; + G4VBiasingOperation *ProposeOccurenceBiasingOperation( + const G4Track *track, + const G4BiasingProcessInterface *callingProcess) final override; /// return the process we want to bias virtual std::string getProcessToBias() const { return "photonNuclear"; } @@ -36,22 +32,21 @@ class PhotoNuclear : public XsecBiasingOperator { virtual std::string getVolumeToBias() const { return volume_; } /// record the configuration into the run header - virtual void RecordConfig(ldmx::RunHeader& h) const { - h.setStringParameter("BiasOperators::PhotoNuclear::Volume", volume_); - h.setFloatParameter("BiasOperators::PhotoNuclear::Threshold", threshold_); - h.setFloatParameter("BiasOperators::PhotoNuclear::Factor", factor_); - h.setIntParameter("BiasOperators::PhotoNuclear::Bias Conv Down", - down_bias_conv_); - h.setIntParameter("BiasOperators::PhotoNuclear::Only Children Of Primary", - only_children_of_primary_); + virtual void RecordConfig(fire::RunHeader &h) const { + h.set("biasings::PhotoNuclear::Volume", volume_); + h.set("biasings::PhotoNuclear::Threshold", threshold_); + h.set("biasings::PhotoNuclear::Factor", factor_); + h.set("biasings::PhotoNuclear::Bias Conv Down", down_bias_conv_); + h.set("biasings::PhotoNuclear::Only Children Of Primary", + only_children_of_primary_); } - private: +private: /** Geant4 gamma conversion process name. */ static const std::string CONVERSION_PROCESS; /** Cross-section biasing operation for conversion process */ - G4BOptnChangeCrossSection* emXsecOperation{nullptr}; + G4BOptnChangeCrossSection *emXsecOperation{nullptr}; /** Unbiased photonuclear xsec. */ double pnXsecUnbiased_{0}; @@ -74,8 +69,5 @@ class PhotoNuclear : public XsecBiasingOperator { /// Should we restrict biasing to only children of primary? bool only_children_of_primary_; -}; // PhotoNuclear -} // namespace biasoperators -} // namespace g4fire - -#endif // SIMCORE_BIASOPERATORS_PHOTONUCLEAR_H_ +}; // PhotoNuclear +} // namespace g4fire::biasing diff --git a/include/g4fire/DarkBrem/APrimePhysics.h b/include/g4fire/darkbrem/APrimePhysics.h similarity index 97% rename from include/g4fire/DarkBrem/APrimePhysics.h rename to include/g4fire/darkbrem/APrimePhysics.h index fe72a35..1344995 100644 --- a/include/g4fire/DarkBrem/APrimePhysics.h +++ b/include/g4fire/darkbrem/APrimePhysics.h @@ -3,7 +3,7 @@ #include "G4VPhysicsConstructor.hh" #include "fire/config/Parameters.h" -#include "g4fire/DarkBrem/G4eDarkBremsstrahlung.h" +#include "g4fire/darkbrem/G4eDarkBremsstrahlung.h" namespace g4fire { diff --git a/include/g4fire/DarkBrem/DarkBremVertexLibraryModel.h b/include/g4fire/darkbrem/DarkBremVertexLibraryModel.h similarity index 99% rename from include/g4fire/DarkBrem/DarkBremVertexLibraryModel.h rename to include/g4fire/darkbrem/DarkBremVertexLibraryModel.h index 553ba50..aefd2ab 100644 --- a/include/g4fire/DarkBrem/DarkBremVertexLibraryModel.h +++ b/include/g4fire/darkbrem/DarkBremVertexLibraryModel.h @@ -3,7 +3,7 @@ #include #include "fire/config/Parameters.h" -#include "g4fire/DarkBrem/G4eDarkBremsstrahlung.h" +#include "g4fire/darkbrem/G4eDarkBremsstrahlung.h" namespace g4fire { diff --git a/include/g4fire/DarkBrem/G4APrime.h b/include/g4fire/darkbrem/G4APrime.h similarity index 100% rename from include/g4fire/DarkBrem/G4APrime.h rename to include/g4fire/darkbrem/G4APrime.h diff --git a/include/g4fire/DarkBrem/G4eDarkBremsstrahlung.h b/include/g4fire/darkbrem/G4eDarkBremsstrahlung.h similarity index 100% rename from include/g4fire/DarkBrem/G4eDarkBremsstrahlung.h rename to include/g4fire/darkbrem/G4eDarkBremsstrahlung.h diff --git a/include/g4fire/event/EventBuilder.h b/include/g4fire/event/EventBuilder.h new file mode 100644 index 0000000..351e964 --- /dev/null +++ b/include/g4fire/event/EventBuilder.h @@ -0,0 +1,28 @@ +#pragma once + +#include "fire/Event.h" + +class G4Event; + +namespace g4fire::event { + +/** + */ +class EventBuilder { + public: + /** + */ + EventBuilder() = default; + + /// Destructor + ~EventBuilder() = default; + + /** + */ + void writeEvent(const G4Event *g4event, fire::Event &event); + + /** + */ + void writeHeader(const G4Event *g4event, fire::Event &event); +}; +} // namespace g4fire::event diff --git a/include/g4fire/event/SimCalorimeterHit.h b/include/g4fire/event/SimCalorimeterHit.h new file mode 100644 index 0000000..79faec1 --- /dev/null +++ b/include/g4fire/event/SimCalorimeterHit.h @@ -0,0 +1,213 @@ +#pragma once + +#include +#include + +#include "g4fire/event/SimParticle.h" + +#include "fire/io/Data.h" + +namespace g4fire::event { + +/** + * @brief Stores simulated calorimeter hit information + * + * @note + * This class represents simulated hit information from a calorimeter detector. + * It provides access to the cell ID, energy deposition, cell position and time. + * Additionally, individual depositions or steps from MC particles are tabulated + * as contributions stored in vectors. Contribution information includes a + * reference to the relevant SimParticle, the PDG code of the actual particle + * which deposited energy (may be different from the actual SimParticle), the + * time of the contribution and the energy deposition. + */ +class SimCalorimeterHit { + public: + /** + * @brief Information about a contribution to the hit in the associated cell + */ + struct Contrib { + /** + * track_id of incident particle that is an ancestor of the contributor + * + * The incident ancestor is found in TrackMap::findIncident where the + * ancestry is looped upwards until a particle is found that matches + * the criteria. + * (1) particle will be saved to output file AND + * (2) particle originates in a region outside the CalorimeterRegion + * If no particle is found matching these criteria, the primary particle + * that is this track_id's ancestor is chosen. + */ + int incident_id{-1}; + + /// track ID of this contributor + int track_id{-1}; + + /// PDG ID of this contributor + int pdg_id{0}; + + /// Energy depostied by this contributor + float edep{0}; + + /// Time this contributor made the hit (global Geant4 time) + float time{0}; + }; + + /// Constructor + SimCalorimeterHit() = default; + + /// Destructor + ~SimCalorimeterHit(); + + /// Reset an instance of this class by clearing all of its data. + void clear(); + + /// @return The detector ID. + int id() const { return id_; } + + /** + * Set the detector ID. + * + * @param[in] id The detector ID. + */ + void setID(const int id) { id_ = id; } + + /// @return The energy deposition of the hit in MeV. + float edep() const { return edep_; } + + /** + * Set the energy deposition of the hit [MeV]. + * + * @param[in] edep The energy deposition of the hit. + */ + void setEdep(const float edep) { edep_ = edep; } + + /// @return A vector containing the x, y, and z positions of the hit in mm. + std::vector position() const { return {x_, y_, z_}; } + + /** + * Set the XYZ position of the hit in mm. + * + * @param[in] x The x position. + * @param[in] y The y position. + * @param[in] z The z position. + */ + void setPosition(const float x, const float y, const float z) { + x_ = x; + y_ = y; + z_ = z; + } + + /// @return The global time of the hit in ns. + float time() const { return time_; } + + /** + * Set the time of the hit [ns]. + * + * @param[in] time The time of the hit. + */ + void setTime(const float time) { time_ = time; } + + /// @return The number of hit contributions. + unsigned contribCount() const { return contrib_count_; } + + /** + * Add a hit contribution from a SimParticle. + * + * @param[in] incident_id The Geant4 track ID for the particle's parent + * incident on the calorimeter region + * @param[in] track_id The Geant4 track ID for the particle + * @param[in] pdg_id The PDG code of the actual track. + * @param[in] edep The energy deposition of the hit [MeV]. + * @param[in] time The time of the hit [ns]. + */ + void addContrib(const int &incident_id, const int &track_id, + const int &pdg_id, const float &edep, const float &time); + + /** + * Get a hit contribution by index. + * + * @param[i] i The index of the hit contribution. + * @return The hit contribution at the index. + */ + Contrib getContrib(const int &i) const; + + /** + * Find the index of a hit contribution from a SimParticle and PDG code. + * + * @param[in] track_id The track ID of the particle causing the hit + * @param[in] pdg_id The PDG code of the contribution. + * @return The index of the contribution or -1 if none exists. + */ + int findContribIndex(const int &track_id, const int &pdg_id) const; + + /** + * Update an existing hit contribution by incrementing its edep and setting + * the time if the new time is less than the old one. + * + * @param[in] i The index of the contribution. + * @param[in] edep The additional energy contribution [MeV]. + * @param[in] time The time of the contribution [ns]. + */ + void updateContrib(const int &i, const float &edep, const float &time); + + /// Sort by time of hit + bool operator<(const SimCalorimeterHit &rhs) const { + return time() < rhs.time(); + } + + /** + * Overload the stream insertion operator to output a string representation + * of this SimCalorimeterHit. + * + * @param[in] output The output stream where the string representation will + * be inserted. + * @param[in] particle The SimCalorimeterHit to stringify. + * + * @return[out] An ostream object with the string representation of + * SimCalorimeterHit inserted. + */ + friend std::ostream &operator<<(std::ostream &output, + const SimCalorimeterHit &hit); + + private: + /// The detector ID. + int id_{0}; + + /// The energy deposition in MeV. + float edep_{0}; + + /// The X position in mm. + float x_{0}; + + /// The Y position in mm. + float y_{0}; + + /// The Z position in mm. + float z_{0}; + + /// The global time of the hit in ns. + float time_{0}; + + /// The list of track IDs contributing to the hit. + std::vector contribs_track_id_; + + /// The list of incident IDs contributing to the hit + std::vector contribs_incident_id_; + + /// The list of PDG codes contributing to the hit. + std::vector contribs_pdg_id_; + + /// The list of energy depositions contributing to the hit. + std::vector contribs_edeps_; + + /// The list of times contributing to the hit. + std::vector contribs_time_; + + /// The number of hit contributions. + unsigned contrib_count_{0}; + + friend class fire::io::Data; + void attach(fire::io::Data &d); +}; // SimCalorimeterHit +} // namespace g4fire::event diff --git a/include/g4fire/Event/SimParticle.h b/include/g4fire/event/SimParticle.h similarity index 65% rename from include/g4fire/Event/SimParticle.h rename to include/g4fire/event/SimParticle.h index 14b4158..23eab89 100644 --- a/include/g4fire/Event/SimParticle.h +++ b/include/g4fire/event/SimParticle.h @@ -1,19 +1,12 @@ -#ifndef SIMCORE_EVENT_SIMPARTICLE_H -#define SIMCORE_EVENT_SIMPARTICLE_H +#pragma one -/*~~~~~~~~~~*/ -/* ROOT */ -/*~~~~~~~~~~*/ -#include "TObject.h" - -/*~~~~~~~~~~~~~~~~*/ -/* C++ StdLib */ -/*~~~~~~~~~~~~~~~~*/ #include #include #include -namespace ldmx { +#include "fire/io/Data.h" + +namespace g4fire::event { /** * Class representing a simulated particle. @@ -23,7 +16,8 @@ namespace ldmx { class SimParticle { public: /** - * Enum for interesting process types. + * Enum for interesting process types. The names of these process + * types are mapped directly to their Geant4 counterpart. */ enum ProcessType { unknown = 0, @@ -45,30 +39,27 @@ class SimParticle { typedef std::map ProcessTypeMap; /// Constructor - SimParticle(); + SimParticle() = default; /// Destructor - ~SimParticle(); + ~SimParticle() = default; /// Reset an instance of this class by clearing all of its data. - void Clear(); - - /// Print a string representation of this object. - void Print() const; + void clear(); /** * Get the energy of this particle [MeV]. * * @return The energy of this particle. */ - double getEnergy() const { return energy_; } + double energy() const { return energy_; } /** * Get the PDG ID of this particle. * * @return The PDG ID of this particle. */ - int getPdgID() const { return pdgID_; } + int pdgID() const { return pdg_id_; } /** * Get the generator status of this particle. A non-zero status @@ -77,14 +68,14 @@ class SimParticle { * * @return The generator status. */ - int getGenStatus() const { return genStatus_; } + short genStatus() const { return gen_status_; } /** * Get the global time of this particle's creation [ns]. * * @return The global time of this particle's creation. */ - double getTime() const { return time_; } + float time() const { return time_; } /** * Get a vector containing the vertex of this particle in mm. @@ -96,7 +87,7 @@ class SimParticle { * * @return The vertex of this particle. */ - std::vector getVertex() const { return {x_, y_, z_}; } + std::vector vertex() const { return {x_, y_, z_}; } /** * Get the volume name in which this particle was created in. @@ -105,7 +96,7 @@ class SimParticle { * * @return The volume name in which this particle was created in. */ - std::string getVertexVolume() const { return vertexVolume_; } + std::string vertexVolume() const { return vertex_vol_; } /** * Get the endpoint of this particle where it was destroyed @@ -113,7 +104,7 @@ class SimParticle { * * @return The endpoint of this particle */ - std::vector getEndPoint() const { return {endX_, endY_, endZ_}; } + std::vector endPoint() const { return {endx_, endy_, endz_}; } /** * Get a vector containing the momentum of this particle [MeV]. @@ -122,21 +113,21 @@ class SimParticle { * * @return The momentum of this particle. */ - std::vector getMomentum() const { return {px_, py_, pz_}; } + std::vector momentum() const { return {px_, py_, pz_}; } /** * Get the mass of this particle [GeV]. * * @return The mass of this particle in GeV. */ - double getMass() const { return mass_; } + float mass() const { return mass_; } /** * Get the charge of this particle. * * @return The charge of this particle. */ - double getCharge() const { return charge_; } + short charge() const { return charge_; } /** * Get a vector containing the track IDs of all daughter particles. @@ -144,42 +135,42 @@ class SimParticle { * @return A vector containing the track IDs of all daughter * particles. */ - std::vector getDaughters() const { return daughters_; } + std::vector daughters() const { return daughters_; } /** * Get a vector containing the track IDs of the parent particles. * * @return A vector containing the track IDs the parent particles. */ - std::vector getParents() const { return parents_; } + std::vector parents() const { return parents_; } /** * Set the energy of this particle [MeV]. * * @param[in] energy the energy of this particle. */ - void setEnergy(const double& energy) { energy_ = energy; } + void setEnergy(const double &energy) { energy_ = energy; } /** * Set the PDG ID of this particle. * - * @param[in] pdgID the PDG ID of the hit. + * @param[in] pdg_id the PDG ID of the hit. */ - void setPdgID(const int& pdgID) { pdgID_ = pdgID; } + void setPdgID(const int &pdg_id) { pdg_id_ = pdg_id; } /** * Set the generator status of this particle. * - * @param[in] genStatus the generator status of the hit. + * @param[in] gen_status the generator status of the hit. */ - void setGenStatus(const int& genStatus) { genStatus_ = genStatus; } + void setGenStatus(const short &gen_status) { gen_status_ = gen_status; } /** * Set the global time of this particle's creation [ns]. * * @param[in] time The global time of this particle's creation. */ - void setTime(const double& time) { time_ = time; } + void setTime(const float &time) { time_ = time; } /** * Set the vertex of this particle [mm]. @@ -188,7 +179,7 @@ class SimParticle { * @param[in] y The vertex y position. * @param[in] z The vertex z position. */ - void setVertex(const double& x, const double& y, const double& z) { + void setVertex(const float &x, const float &y, const float &z) { x_ = x; y_ = y; z_ = z; @@ -197,24 +188,24 @@ class SimParticle { /** * Set the name of the volume that this particle was created in. * - * @param[in] vertexVolume volume name that this particle was + * @param[in] vertex_vol volume name that this particle was * created in. */ - void setVertexVolume(const std::string& vertexVolume) { - vertexVolume_ = vertexVolume; + void setVertexVolume(const std::string &vertex_vol) { + vertex_vol_ = vertex_vol; } /** * Set the end point position of this particle [mm]. * - * @param[in] endX The x position of the end point. - * @param[in] endY The y position of the end point. - * @param[in] endZ The z position of the end point. + * @param[in] endx The x position of the end point. + * @param[in] endy The y position of the end point. + * @param[in] endz The z position of the end point. */ - void setEndPoint(const double& endX, const double& endY, const double& endZ) { - endX_ = endX; - endY_ = endY; - endZ_ = endZ; + void setEndPoint(const float &endx, const float &endy, const float &endz) { + endx_ = endx; + endy_ = endy; + endz_ = endz; } /** @@ -224,7 +215,7 @@ class SimParticle { * @param[in] py The y momentum component. * @param[in] pz The z momentum component. */ - void setMomentum(const double& px, const double& py, const double& pz) { + void setMomentum(const double &px, const double &py, const double &pz) { px_ = px; py_ = py; pz_ = pz; @@ -235,14 +226,14 @@ class SimParticle { * * @param[in] mass The mass of this particle. */ - void setMass(const double& mass) { mass_ = mass; } + void setMass(const float &mass) { mass_ = mass; } /** * Set the charge of this particle. * * @param[in] charge The charge of this particle. */ - void setCharge(const double& charge) { charge_ = charge; } + void setCharge(const short &charge) { charge_ = charge; } /** * Add a reference to a daughter particle by its track ID. @@ -252,7 +243,7 @@ class SimParticle { * * @param[in] daughterTrackID The daughter particle track ID. */ - void addDaughter(const int& daughterTrackID) { + void addDaughter(const short &daughterTrackID) { daughters_.push_back(daughterTrackID); } @@ -261,7 +252,7 @@ class SimParticle { * * @param[in] parentTrackID The track ID of the parent particle. */ - void addParent(const int& parentTrackID) { + void addParent(const short &parentTrackID) { parents_.push_back(parentTrackID); } @@ -270,14 +261,14 @@ class SimParticle { * * @return The creator process type of this particle. */ - int getProcessType() const { return processType_; } + int processType() const { return process_type_; } /** * Set the creator process type of this particle. * - * @param[in] processType the creator process type of this particle. + * @param[in] process_type the creator process type of this particle. */ - void setProcessType(const int& processType) { processType_ = processType; } + void setProcessType(const short &process_type) { process_type_ = process_type; } /** * Set the momentum at this particle's end point. @@ -286,8 +277,8 @@ class SimParticle { * @param[in] endpy The y component of the momentum. * @param[in] endpz The z component of the momentum. */ - void setEndPointMomentum(const double& endpx, const double& endpy, - const double& endpz) { + void setEndPointMomentum(const double &endpx, const double &endpy, + const double &endpz) { endpx_ = endpx; endpy_ = endpy; endpz_ = endpz; @@ -298,7 +289,7 @@ class SimParticle { * * @return The momentum at this particle's end point as a vector. */ - std::vector getEndPointMomentum() const { + std::vector endPointMomentum() const { return {endpx_, endpy_, endpz_}; } @@ -309,39 +300,52 @@ class SimParticle { */ static ProcessType findProcessType(std::string processName); + /** + * Overload the stream insertion operator to output a string representation + * of this SimParticle. + * + * @param[in] output The output stream where the string representation will + * be inserted. + * @param[in] particle The SimParticle to stringify. + * + * @return[out] An ostream object with the string representation of + * SimParticle inserted. + */ + friend std::ostream &operator<<(std::ostream &output, + const SimParticle &particle); + private: static ProcessTypeMap createProcessTypeMap(); - private: /// The energy of this particle. double energy_{0}; /// The PDG ID of this particle. - int pdgID_{0}; + int pdg_id_{0}; /// The generator status. - int genStatus_{-1}; + short gen_status_{-1}; /// The global creation time. - double time_{0}; + float time_{0}; /// The x component of the vertex. - double x_{0}; + float x_{0}; /// The y component of the vertex. - double y_{0}; + float y_{0}; /// The z component of the vertex. - double z_{0}; + float z_{0}; /// The x component of the end point. - double endX_{0}; + float endx_{0}; /// The y component of the end point. - double endY_{0}; + float endy_{0}; /// The z component of the end point. - double endZ_{0}; + float endz_{0}; /// The x component of the momentum. double px_{0}; @@ -362,29 +366,28 @@ class SimParticle { double endpz_{0}; /// The particle's mass. - double mass_{0}; + float mass_{0}; /// The particle's charge. - double charge_{0}; + short charge_{0}; /// The list of daughter particle track IDs. - std::vector daughters_; + std::vector daughters_; /// The list of parent particles track IDs. - std::vector parents_; + std::vector parents_; /// Encoding of Geant4 process type. - int processType_{-1}; + short process_type_{-1}; /// Volume the track was created in. - std::string vertexVolume_{""}; + std::string vertex_vol_{""}; /// Map containing the process types. static ProcessTypeMap PROCESS_MAP; - ClassDef(SimParticle, 7); - -}; // SimParticle -} // namespace ldmx + friend class fire::io::Data; + void attach(fire::io::Data &d); -#endif // SIMCORE_EVENT_SIMPARTICLE_H +}; // SimParticle +} // namespace g4fire::event diff --git a/include/g4fire/event/SimTrackerHit.h b/include/g4fire/event/SimTrackerHit.h new file mode 100644 index 0000000..9fa5c41 --- /dev/null +++ b/include/g4fire/event/SimTrackerHit.h @@ -0,0 +1,235 @@ +#pragma once + +#include + +#include "fire/io/Data.h" + +namespace g4fire::event { + +/** + * @brief Represents a simulated tracker hit in the simulation + */ +class SimTrackerHit { +public: + /// Constructor + SimTrackerHit() = default; + + /// Destructor + ~SimTrackerHit() = default; + + /// Reset an instance of this class by clearing all of its data. + void clear(); + + /// @return The detector ID of the hit. + int id() const { return id_; }; + + /// @return The geometric layer ID of the hit. + short layerID() const { return layer_id_; }; + + /** + * Get the module ID associated with a hit. This is used to + * uniquely identify a sensor within a layer. + * + * @return The module ID associated with a hit. + */ + short moduleID() const { return module_id_; }; + + /** + * Get the (x, y, z) position of the hit in mm. + * + * @return A vector representation of the x, y, z position of a hit. + */ + std::vector position() const { return {x_, y_, z_}; }; + + /// @return The energy deposited on the hit in MeV. + float edep() const { return edep_; }; + + /// @return The energy of the hit in MeV. + float energy() const { return energy_; }; + + /// @return The global time of the hit in ns. + float time() const { return time_; }; + + /** + * Get the path length between the start and end points of the + * hit [mm]. + * + * @return The path length of the hit. + */ + float pathLength() const { return path_length_; }; + + /** + * Get the momentum of the particle at the position at which + * the hit took place [MeV]. + * + * @return A vector representation of the momentum of the particle. + */ + std::vector momentum() const { return {px_, py_, pz_}; }; + + /// @return The Geant4 track ID of the particle associated with this hit. + short trackID() const { return track_id_; }; + + /// @return The PDG ID of the SimParticle associated with this hit. + int pdgID() const { return pdg_id_; }; + + /** + * Set the detector ID of the hit. + * + * @param[in] id The detector ID of the hit. + */ + void setID(const long id) { id_ = id; }; + + /** + * Set the geometric layer ID of the hit. + * + * @param[in] layer_id The layer ID of the hit. + */ + void setLayerID(const int layer_id) { layer_id_ = layer_id; }; + + /** + * Set the module ID associated with a hit. This is used to + * uniquely identify a sensor within a layer. + * + * @param[in] module_id The module ID associated with a hit. + */ + void setModuleID(const int module_id) { module_id_ = module_id; }; + + /** + * Set the position of the hit in mm. + * + * @param[in] x The x position. + * @param[in] y The y position. + * @param[in] z The z position. + */ + inline void setPosition(const float x, const float y, const float z) { + x_ = x; + y_ = y; + z_ = z; + } + + /** + * Set the energy deposited on the hit in MeV. + * + * @param[in] edep The energy deposited on the hit. + */ + void setEdep(const float edep) { edep_ = edep; }; + + /** + * Set the energy of the hit. + * + * @param[in] energy The energy of the hit. + */ + void setEnergy(const float energy) { energy_ = energy; }; + + /** + * Set the global time of the hit in ns. + * + * @param[in] time The global time of the hit. + */ + void setTime(const float time) { time_ = time; }; + + /** + * Set the path length of the hit in mm. + * + * @param[in] path_length The path length of the hit. + */ + void setPathLength(const float path_length) { path_length_ = path_length; }; + + /** + * Set the momentum of the particle at the position at which + * the hit took place in MeV. + * + * @param[in] px The X momentum. + * @param[in] py The Y momentum. + * @param[in] pz The Z momentum. + */ + inline void setMomentum(const float px, const float py, const float pz) { + px_ = px; + py_ = py; + pz_ = pz; + }; + + /** + * Set the Geant4 track ID of the particle associted with this hit. + * + * @param[in] track_id The Sim particle track ID of the hit. + */ + void setTrackID(const int track_id) { track_id_ = track_id; }; + + /** + * Set the PDG ID of the SimParticle associated with this hit. + * + * @param[in] pdg_id The PDG ID of the SimParticle associated with this hit. + */ + void setPdgID(const int pdg_id) { pdg_id_ = pdg_id; }; + + /// Sort by time of hit + bool operator<(const g4fire::event::SimTrackerHit &rhs) const { + return time() < rhs.time(); + } + + /** + * Overload the stream insertion operator to output a string representation + * of this SimTrackerHit. + * + * @param[in] output The output stream where the string representation will + * be inserted. + * @param[in] particle The SimTrackerHit to stringify. + * + * @return[out] An ostream object with the string representation of + * SimTrackerHit inserted. + */ + friend std::ostream &operator<<(std::ostream &output, + const SimTrackerHit &hit); + +private: + /// The detector ID. + int id_{0}; + + /// The layer ID. + int layer_id_{0}; + + /// The module ID. + int module_id_{0}; + + /// The energy deposited on the hit. + float edep_{0}; + + /// The global time of the hit. + float time_{0}; + + /// The x component of the momentum. + float px_{0}; + + /// The y component of the momentum. + float py_{0}; + + /// The z component of the momentum. + float pz_{0}; + + /// The hit energy. + float energy_{0}; + + /// The x position. + float x_{0}; + + /// The y position. + float y_{0}; + + /// The z position. + float z_{0}; + + /// The path length of the hit. + float path_length_{0}; + + /// The Geant4 track ID. + int track_id_{0}; + + /// The PDG ID of the particle associated with this track. + int pdg_id_{0}; + + friend class fire::io::Data; + void attach(fire::io::Data &d); + +}; // SimTrackerHit +} // namespace g4fire::event diff --git a/include/g4fire/Geo/AuxInfoReader.h b/include/g4fire/geo/AuxInfoReader.h similarity index 100% rename from include/g4fire/Geo/AuxInfoReader.h rename to include/g4fire/geo/AuxInfoReader.h diff --git a/include/g4fire/Geo/GDMLParser.h b/include/g4fire/geo/GDMLParser.h similarity index 96% rename from include/g4fire/Geo/GDMLParser.h rename to include/g4fire/geo/GDMLParser.h index dcaf0a1..0add80d 100644 --- a/include/g4fire/Geo/GDMLParser.h +++ b/include/g4fire/geo/GDMLParser.h @@ -5,8 +5,8 @@ #include "fire/config/Parameters.h" -#include "g4fire/Geo/AuxInfoReader.h" -#include "g4fire/Geo/Parser.h" +#include "g4fire/geo/AuxInfoReader.h" +#include "g4fire/geo/Parser.h" class G4VPhysicalVolume; diff --git a/include/g4fire/Geo/Parser.h b/include/g4fire/geo/Parser.h similarity index 100% rename from include/g4fire/Geo/Parser.h rename to include/g4fire/geo/Parser.h diff --git a/include/g4fire/Geo/ParserFactory.h b/include/g4fire/geo/ParserFactory.h similarity index 97% rename from include/g4fire/Geo/ParserFactory.h rename to include/g4fire/geo/ParserFactory.h index bb7709a..2f59a12 100644 --- a/include/g4fire/Geo/ParserFactory.h +++ b/include/g4fire/geo/ParserFactory.h @@ -1,7 +1,7 @@ #ifndef G4FIRE_GEO_PARSERFACTORY_H_ #define G4FIRE_GEO_PARSERFACTORY_H_ -#include "g4fire/Geo/Parser.h" +#include "g4fire/geo/Parser.h" #include diff --git a/python/g4fire/_bias_operators.py.in b/python/g4fire/_bias_operators.py.in new file mode 100644 index 0000000..c03f101 --- /dev/null +++ b/python/g4fire/_bias_operators.py.in @@ -0,0 +1,190 @@ +"""Biasing operator templates for use throughout g4fire + +Mainly focused on reducing the number of places that certain parameter and class names +are hardcoded into the python configuration. +""" + + + +class XsecBiasingOperator: + """Object that stores parameters for a XsecBiasingOperator + + Parameters + ---------- + instance_name : str + Unique name for this particular instance of a XsecBiasingOperator + class_name : str + Name of C++ class that this XsecBiasingOperator should be + module_name : str, optional + Name of module that the class was compiled into + """ + def __init__(self, instance_name, class_name, module_name='g4fire'): + self.class_name = class_name + self.instance_name = instance_name + + from fire.cfg._process import Process + Process.addLibrary('@CMAKE_INSTALL_PREFIX@/lib/lib%s.so' % module_name) + + def __str__(self): + """Stringify this XsecBiasingOperator + + Returns + ------- + str + A human-readable version of this XsecBiasingOperator printing all its attributes + """ + + string = "XsecBiasingOperator (" + self.__repr__() + ") {" + for k, v in self.__dict__.items(): + string += " %s : %s" % (k, v) + string += " }" + + return string + + def __repr__(self): + """A shorter string representation of this XsecBiasingOperator + + Returns + ------- + str + Just printing its instance and class names + """ + + return '%s of class %s' % (self.instance_name, self.class_name) + + +class PhotoNuclear(XsecBiasingOperator): + """Bias photo nuclear process + + Parameters + ---------- + vol : str + name of volume to bias within + factor : float + biasing factor to mutliply by + threshold : float, optional + minimum kinetic energy [MeV] to bias tracks + down_bias_conv : bool, optional + Should we down bias the gamma conversion process? + only_children_of_primary : bool, optional + Should we only bias photons that are children of the primary particle? + """ + def __init__(self, + vol, + factor, + thresh=0., + down_bias_conv=True, + only_children_of_primary=False): + super().__init__('%s_bias_pn' % vol, + 'g4fire::biasing::PhotoNuclear') + + self.volume = vol + self.factor = factor + self.threshold = thresh + self.down_bias_conv = down_bias_conv + self.only_children_of_primary = only_children_of_primary + + +class ElectroNuclear(XsecBiasingOperator): + """Bias electro nuclear process + + Parameters + ---------- + vol : str + name of volume to bias within + factor : float + biasing factor to mutliply by + threshold : float, optional + minimum kinetic energy [MeV] to bias tracks + """ + def __init__(self, vol, factor, thresh=0.): + super().__init__('%s_bias_en' % vol, + 'g4fire::biasing::ElectroNuclear') + + self.volume = vol + self.factor = factor + self.threshold = thresh + + +class GammaToMuPair(XsecBiasingOperator): + """Bias gamma -> mu+ mu- process + + Parameters + ---------- + vol : str + name of volume to bias within + factor : float + biasing factor to multiply by + threshold : float, optional + minimum kinetic energy [MeV] to bias tracks + """ + def __init__(self, vol, factor, thresh=0.): + super().__init__('%s_bias_mumu' % vol, + 'g4fire::biasing::GammaToMuPair') + + self.volume = vol + self.factor = factor + self.threshold = thresh + + +class NeutronInelastic(XsecBiasingOperator): + """Bias neutron inelastic collisions + + Parameters + ---------- + vol : str + name of volume to bias within + factor : float + biasing factor to multiply by + threshold : float, optional + minimum kinetic energy [MeV] to bias tracks + """ + def __init__(self, vol, factor, thresh=0.): + super().__init__('%s_bias_neutroninelastic' % vol, + 'g4fire::biasing::NeutronInelastic') + + self.volume = vol + self.factor = factor + self.threshold = thresh + + +class K0LongInelastic(XsecBiasingOperator): + """Bias K0 Long inelastic collisions + + Parameters + ---------- + vol : str + name of volume to bias within + factor : float + biasing factor to multiply by + threshold : float, optional + minimum kinetic energy [MeV] to bias tracks + """ + def __init__(self, vol, factor, thresh=0.): + super().__init__('%s_bias_k0longinelastic' % vol, + 'g4fire::biasing::K0LongInelastic') + + self.volume = vol + self.factor = factor + self.threshold = thresh + + +class DarkBrem(XsecBiasingOperator): + """Bias dark brem process + + Parameters + ---------- + vol : str + name of volume to bias within + bial_all : bool + Should we bias all electrons or just the primary? + factor : float + biasing factor to mutliply by + """ + def __init__(self, vol, bias_all, factor): + super().__init__('%s_bias_darkbrem' % vol, + 'g4fire::biasing::DarkBrem') + + self.volume = vol + self.bias_all = bias_all + self.factor = factor diff --git a/python/g4fire/_dark_brem.py b/python/g4fire/_dark_brem.py new file mode 100644 index 0000000..ca82199 --- /dev/null +++ b/python/g4fire/_dark_brem.py @@ -0,0 +1,110 @@ +"""Configuration module for dark brem simulation""" + +class DarkBremModel() : + """Storage for parameters of a dark brem model + + All other models should inherit from this class + in order to keep the correct internal parameters. + + Parameters + ---------- + name : str + Name of this dark brem model + """ + + def __init__(self,name) : + self.name = name + + def __str__(self) : + string = '%s {'%self.name + for key in self.__dict__ : + if key is not self.name : + string += ' %s=%s'%( key , self.__dict__[key] ) + string += ' }' + return string + +class VertexLibraryModel(DarkBremModel) : + """Configuration for the vertex library dark brem model + + Parameters + ---------- + library_path : str + Path to directory of LHE files containing dark brem vertices + + Attributes + ---------- + method : str + Interpretation method for LHE files + threshold : float + Minimum energy [GeV] that electron should have for dark brem to have nonzero xsec + epsilon : float + Epsilon for dark brem xsec calculation + """ + + def __init__(self, library_path) : + super().__init__('vertex_library') + self.library_path = library_path + self.method = 'forward_only' + self.threshold = 2.0 #GeV + self.epsilon = 0.01 + +class DarkBrem: + """Storage for parameters of dark brem process + + Attributes + ---------- + ap_mass : float + Mass of A' in MeV + enable : bool + Should we use the custom Geant4 dark brem process? (Default: No) + only_one_per_event : bool + Should we deactivate the process after one dark brem or allow for more than one? (Default: No) + cache_xsec : bool + Should we cache the xsec's computed from the model? (Default: yes) + model : DarkBremModel + The model that should be use for dark bremming + """ + + def __init__(self) : + self.ap_mass = 0. + self.only_one_per_event = False + self.enable = False #off by default + self.cache_xsec = True + self.model = DarkBremModel('UNDEFINED') + + def activate(self, ap_mass, model = None) : + """Activate the dark brem process with the input A' mass [MeV] and dark brem model + + If no dark brem model is given, we do not activate the process + and only define the A' mass. This allows for some backwards + compatibility by allowing users to use the LHEPrimaryGenerator + with A' particles. + """ + + self.ap_mass = ap_mass + + if model is not None : + if not isinstance(model,DarkBremModel) : + raise Exception('Dark brem process needs to be configured with an associated DarkBremModel.') + + self.enable = True + self.model = model + + def __str__(self): + """Stringify the DarkBrem configuration + + Returns + ------- + str + A human-readable version of all its attributes + """ + + string = "{ Enabled: %r"%self.enable + if self.enable : + string += ", Mass: %.1f MeV"%self.ap_mass + string += ", Only One Per Event: %r"%self.only_one_per_event + string += ", Model: %s"%self.model + + string += " }" + + return string diff --git a/python/g4fire/_simulator.py b/python/g4fire/_simulator.py index 2b5fc80..d37202d 100644 --- a/python/g4fire/_simulator.py +++ b/python/g4fire/_simulator.py @@ -5,7 +5,7 @@ """ from fire.cfg import Processor - +from g4fire import _dark_brem class Simulator(Processor): """A instance of the simulation configuration @@ -73,7 +73,9 @@ def __init__(self, instance_name, detector, description, generators, biasing_operators=[], logging_prefix='', validate_detector=False, - verbosity = 0): + verbosity = 0, + dark_brem = _dark_brem.DarkBrem() + ): super().__init__(instance_name, "g4fire::Simulator", detector=detector, @@ -91,11 +93,9 @@ def __init__(self, instance_name, detector, description, generators, biasing_operators=biasing_operators, logging_prefix=logging_prefix, validate_detector=validate_detector, - verbosity=verbosity) - - #Dark Brem stuff - #from LDMX.g4fire import dark_brem - #self.dark_brem = dark_brem.DarkBrem() + verbosity=verbosity, + dark_brem=dark_brem + ) #def setDetector(self, det_name , include_scoring_planes = False ) : """Set the detector description with the option to include the scoring planes diff --git a/src/g4fire/BiasOperators/NeutronInelastic.cxx b/src/g4fire/BiasOperators/NeutronInelastic.cxx deleted file mode 100644 index 24c4e19..0000000 --- a/src/g4fire/BiasOperators/NeutronInelastic.cxx +++ /dev/null @@ -1,37 +0,0 @@ - -#include "g4fire/BiasOperators/NeutronInelastic.h" - -namespace g4fire { -namespace biasoperators { - -NeutronInelastic::NeutronInelastic(std::string name, const framework::config::Parameters& p) - : XsecBiasingOperator(name, p) { - volume_ = p.getParameter("volume"); - factor_ = p.getParameter("factor"); - threshold_ = p.getParameter("threshold"); -} - -G4VBiasingOperation* NeutronInelastic::ProposeOccurenceBiasingOperation( - const G4Track* track, const G4BiasingProcessInterface* callingProcess) { - - if (track->GetKineticEnergy() < threshold_) return 0; - - std::string currentProcess = - callingProcess->GetWrappedProcess()->GetProcessName(); - if (currentProcess.compare(this->getProcessToBias()) == 0) { - G4double interactionLength = - callingProcess->GetWrappedProcess()->GetCurrentInteractionLength(); - - double neutInXsecUnbiased = 1. / interactionLength; - - double neutInXsecBiased = neutInXsecUnbiased * factor_; - - return BiasedXsec(neutInXsecBiased); - } else - return 0; -} - -} // namespace biasoperators -} // namespace g4fire - -DECLARE_XSECBIASINGOPERATOR(g4fire::biasoperators, NeutronInelastic) diff --git a/src/g4fire/CalorimeterSD.cxx b/src/g4fire/CalorimeterSD.cxx deleted file mode 100644 index f62c224..0000000 --- a/src/g4fire/CalorimeterSD.cxx +++ /dev/null @@ -1,49 +0,0 @@ -#include "g4fire/CalorimeterSD.h" - -// STL -#include - -// Geant4 -#include "G4ChargedGeantino.hh" -#include "G4Geantino.hh" -#include "G4SDManager.hh" -#include "G4Step.hh" -#include "G4StepPoint.hh" - -namespace g4fire { - -CalorimeterSD::CalorimeterSD(G4String name, G4String theCollectionName) - : G4VSensitiveDetector(name), hitsCollection_(0) { - // Add the collection name to vector of names. - this->collectionName.push_back(theCollectionName); - - // Register this SD with the manager. - G4SDManager::GetSDMpointer()->AddNewDetector(this); -} - -CalorimeterSD::~CalorimeterSD() {} - -void CalorimeterSD::Initialize(G4HCofThisEvent* hce) { - // Setup hits collection and the HC ID. - hitsCollection_ = - new G4CalorimeterHitsCollection(SensitiveDetectorName, collectionName[0]); - G4int hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); - hce->AddHitsCollection(hcID, hitsCollection_); -} - -void CalorimeterSD::EndOfEvent(G4HCofThisEvent*) { - // Print number of hits. - if (this->verboseLevel > 0) { - std::cout << GetName() << " had " << hitsCollection_->entries() - << " hits in event" << std::endl; - } - - // Print each hit in hits collection. - if (this->verboseLevel > 1) { - for (unsigned iHit = 0; iHit < hitsCollection_->GetSize(); iHit++) { - (*hitsCollection_)[iHit]->Print(); - } - } -} - -} // namespace g4fire diff --git a/src/g4fire/EcalHitIO.cxx b/src/g4fire/EcalHitIO.cxx deleted file mode 100644 index dd1eead..0000000 --- a/src/g4fire/EcalHitIO.cxx +++ /dev/null @@ -1,103 +0,0 @@ -#include "g4fire/EcalHitIO.h" - -#include "g4fire/UserTrackingAction.h" //to get handle on track map - -// STL -#include - -// LDMX -#include "g4fire/Event/SimCalorimeterHit.h" -#include "g4fire/Event/SimParticle.h" - -namespace g4fire { - -void EcalHitIO::configure(const framework::config::Parameters& ps) { - enableHitContribs_ = ps.getParameter("enableHitContribs"); - compressHitContribs_ = ps.getParameter("compressHitContribs"); -} - -void EcalHitIO::writeHitsCollection( - G4CalorimeterHitsCollection* hc, - std::vector& outputColl) { - const ldmx::EcalHexReadout& hexReadout = - conditionsIntf_.getCondition( - ldmx::EcalHexReadout::CONDITIONS_OBJECT_NAME); - - // get ancestral mapping of tracks - auto trackMap{UserTrackingAction::getUserTrackingAction()->getTrackMap()}; - - int nHits = hc->GetSize(); - std::map hitMap; - - // Loop over input hits from Geant4. - for (int iHit = 0; iHit < nHits; iHit++) { - // Get the hit and its ID. - G4CalorimeterHit* g4hit = (G4CalorimeterHit*)hc->GetHit(iHit); - int hitID = g4hit->getID(); - - // See if hit exists in map already. - std::map::iterator it = hitMap.find(hitID); - - // Is it a new hit? - if (it == hitMap.end()) { - // Create sim hit and assign the ID. - hitMap[hitID] = ldmx::SimCalorimeterHit(); - hitMap[hitID].setID(hitID); - - /** - * Assign XY position to the hit using the ECal hex readout. - * Z position is set from the original hit, which should be the middle of - * the sensor. - */ - double x, y, z; - hexReadout.getCellAbsolutePosition(hitID, x, y, z); - hitMap[hitID].setPosition(x, y, z); - } - - // Get info from the G4 hit. - float edep = g4hit->getEdep(); - float time = g4hit->getTime(); - int pdgCode = g4hit->getPdgCode(); - - // Is hit contrib output enabled? - if (enableHitContribs_) { - // Find the SimParticle associated with this hit. - int trackID = g4hit->getTrackID(); - - // Find if there is an existing hit contrib. - int contribIndex = hitMap[hitID].findContribIndex(trackID, pdgCode); - - // Is contrib output being compressed and a record exists for this - // SimParticle and PDG code? - if (compressHitContribs_ && contribIndex != -1) { - // Update an existing hit contrib. - hitMap[hitID].updateContrib(contribIndex, edep, time); - - } else { - // Add a hit contrib because all steps are being saved or there is not - // an existing record. - hitMap[hitID].addContrib(trackMap->findIncident(trackID), trackID, - pdgCode, edep, time); - } - - } else { - // Hit contributions are not being saved so manually increment the edep - // and set time. - hitMap[hitID].setEdep(hitMap[hitID].getEdep() + edep); - if (time < hitMap[hitID].getTime() || hitMap[hitID].getTime() == 0) { - hitMap[hitID].setTime(time); - } - - } // contrib output enabled or not - } // loop through geant4 hits - - // copy aggregated hits into output collection - outputColl.clear(); - for (auto& mapHit : hitMap) { - outputColl.push_back(mapHit.second); - } - - return; -} - -} // namespace g4fire diff --git a/src/g4fire/EcalSD.cxx b/src/g4fire/EcalSD.cxx deleted file mode 100644 index 79232ef..0000000 --- a/src/g4fire/EcalSD.cxx +++ /dev/null @@ -1,144 +0,0 @@ -#include "g4fire/EcalSD.h" - -// Geant4 -#include "G4ChargedGeantino.hh" -#include "G4Geantino.hh" -#include "G4Polyhedron.hh" -#include "G4Step.hh" -#include "G4StepPoint.hh" -#include "G4VSolid.hh" - -/*~~~~~~~~~~~~~~*/ -/* DetDescr */ -/*~~~~~~~~~~~~~~*/ -#include "DetDescr/EcalID.h" - -namespace g4fire { - -EcalSD::EcalSD(G4String name, G4String theCollectionName, int subDetID, - ConditionsInterface& ci) - : CalorimeterSD(name, theCollectionName), conditionsIntf_(ci) {} - -EcalSD::~EcalSD() {} - -G4bool EcalSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) { - const ldmx::EcalHexReadout& hitMap = - conditionsIntf_.getCondition( - ldmx::EcalHexReadout::CONDITIONS_OBJECT_NAME); - - // Determine if current particle of this step is a Geantino. - G4ParticleDefinition* pdef = aStep->GetTrack()->GetDefinition(); - bool isGeantino = false; - if (pdef == G4Geantino::Definition() || - pdef == G4ChargedGeantino::Definition()) { - isGeantino = true; - } - - // Get the edep from the step. - G4double edep = aStep->GetTotalEnergyDeposit(); - - // Skip steps with no energy dep which come from non-Geantino particles. - if (edep == 0.0 && !isGeantino) { - if (verboseLevel > 2) { - G4cout << "CalorimeterSD skipping step with zero edep." << G4endl - << G4endl; - } - return false; - } - - // Create a new cal hit. - G4CalorimeterHit* hit = new G4CalorimeterHit(); - - // Set the edep. - hit->setEdep(edep); - - // Compute the hit position using the utility function. - G4ThreeVector hitPosition = getHitPosition(aStep); - hit->setPosition(hitPosition.x(), hitPosition.y(), hitPosition.z()); - - // Set the global time. - hit->setTime(aStep->GetTrack()->GetGlobalTime()); - - // Create the ID for the hit. - int cpynum = aStep->GetPreStepPoint() - ->GetTouchableHandle() - ->GetHistory() - ->GetVolume(layerDepth_) - ->GetCopyNo(); - int layerNumber; - layerNumber = int(cpynum / 7); - int module_position = cpynum % 7; - - ldmx::EcalID partialId = - hitMap.getCellModuleID(hitPosition[0], hitPosition[1]); - ldmx::EcalID id(layerNumber, module_position, partialId.cell()); - hit->setID(id.raw()); - - // Set the track ID on the hit. - hit->setTrackID(aStep->GetTrack()->GetTrackID()); - - // Set the PDG code from the track. - hit->setPdgCode(aStep->GetTrack()->GetParticleDefinition()->GetPDGEncoding()); - - if (this->verboseLevel > 2) { - G4cout << "Created new SimCalorimeterHit in detector " << this->GetName() - << " with subdet ID " << id << " ..."; - hit->Print(); - G4cout << G4endl; - } - - // Insert the hit into the hits collection. - hitsCollection_->insert(hit); - - return true; -} - -G4ThreeVector EcalSD::getHitPosition(G4Step* aStep) { - /** - * Set initial hit position from midpoint of the step. - */ - G4StepPoint* prePoint = aStep->GetPreStepPoint(); - G4StepPoint* postPoint = aStep->GetPostStepPoint(); - G4ThreeVector position = - 0.5 * (prePoint->GetPosition() + postPoint->GetPosition()); - - /* - * Get the volume position in global coordinates, which for the ECal is the - * center of the front face of the sensor. - */ - G4ThreeVector volumePosition = aStep->GetPreStepPoint() - ->GetTouchableHandle() - ->GetHistory() - ->GetTopTransform() - .Inverse() - .TransformPoint(G4ThreeVector()); - - // Get the solid from this step. - G4VSolid* solid = prePoint->GetTouchableHandle() - ->GetVolume() - ->GetLogicalVolume() - ->GetSolid(); - auto it = polyMap_.find(solid); - G4Polyhedron* poly; - if (it == polyMap_.end()) { - poly = solid->CreatePolyhedron(); - polyMap_[solid] = poly; - } else { - poly = polyMap_[solid]; - } - - /** - * Use facet info of the solid for setting Z to the sensor midpoint. - */ - G4Point3D iNodes[4]; - G4int n; - poly->GetFacet(1, n, iNodes); - G4double zstart = iNodes[1][2]; - G4double zend = iNodes[0][2]; - G4double zmid = (zstart - zend) / 2; - position.setZ(volumePosition.z() + zmid); - - return position; -} - -} // namespace g4fire diff --git a/src/g4fire/Event/SimCalorimeterHit.cxx b/src/g4fire/Event/SimCalorimeterHit.cxx deleted file mode 100644 index c8b66e9..0000000 --- a/src/g4fire/Event/SimCalorimeterHit.cxx +++ /dev/null @@ -1,81 +0,0 @@ -#include "g4fire/Event/SimCalorimeterHit.h" - -// STL -#include - -ClassImp(ldmx::SimCalorimeterHit) - - namespace ldmx { - SimCalorimeterHit::SimCalorimeterHit() {} - - SimCalorimeterHit::~SimCalorimeterHit() {} - - void SimCalorimeterHit::Clear() { - incidentIDContribs_.clear(); - trackIDContribs_.clear(); - pdgCodeContribs_.clear(); - edepContribs_.clear(); - timeContribs_.clear(); - - nContribs_ = 0; - id_ = 0; - edep_ = 0; - x_ = 0; - y_ = 0; - z_ = 0; - time_ = 0; - } - - void SimCalorimeterHit::Print() const { - std::cout << "SimCalorimeterHit { " - << "id: " << id_ << ", edep: " << edep_ - << ", " - "position: ( " - << x_ << ", " << y_ << ", " << z_ - << " ), num contribs: " << nContribs_ << " }" << std::endl; - } - - void SimCalorimeterHit::addContrib(int incidentID, int trackID, int pdgCode, - float edep, float time) { - incidentIDContribs_.push_back(incidentID); - trackIDContribs_.push_back(trackID); - pdgCodeContribs_.push_back(pdgCode); - edepContribs_.push_back(edep); - timeContribs_.push_back(time); - edep_ += edep; - if (time < time_ || time_ == 0) { - time_ = time; - } - ++nContribs_; - } - - SimCalorimeterHit::Contrib SimCalorimeterHit::getContrib(int i) const { - Contrib contrib; - contrib.incidentID = incidentIDContribs_.at(i); - contrib.trackID = trackIDContribs_.at(i); - contrib.edep = edepContribs_.at(i); - contrib.time = timeContribs_.at(i); - contrib.pdgCode = pdgCodeContribs_.at(i); - return contrib; - } - - int SimCalorimeterHit::findContribIndex(int trackID, int pdgCode) const { - int contribIndex = -1; - for (int iContrib = 0; iContrib < nContribs_; iContrib++) { - Contrib contrib = getContrib(iContrib); - if (contrib.trackID == trackID && contrib.pdgCode == pdgCode) { - contribIndex = iContrib; - break; - } - } - return contribIndex; - } - - void SimCalorimeterHit::updateContrib(int i, float edep, float time) { - this->edepContribs_[i] += edep; - if (time < this->timeContribs_.at(i)) { - this->timeContribs_[i] = time; - } - edep_ += edep; - } -} // namespace ldmx diff --git a/src/g4fire/Event/SimParticle.cxx b/src/g4fire/Event/SimParticle.cxx deleted file mode 100644 index a0cc224..0000000 --- a/src/g4fire/Event/SimParticle.cxx +++ /dev/null @@ -1,103 +0,0 @@ -#include "g4fire/Event/SimParticle.h" - -/*~~~~~~~~~~~~~~~~*/ -/* C++ StdLib */ -/*~~~~~~~~~~~~~~~~*/ -#include - -ClassImp(ldmx::SimParticle) - - namespace ldmx { - SimParticle::ProcessTypeMap SimParticle::createProcessTypeMap() { - ProcessTypeMap procMap; - /// e Z --> e Z gamma - procMap["eBrem"] = ProcessType::eBrem; - /// gamma --> e+ e- - procMap["conv"] = ProcessType::conv; - /// e+ e- --> gamma gamma - procMap["annihil"] = ProcessType::annihil; - /// gamma e --> gamma e - procMap["compt"] = ProcessType::compt; - /// gamma Z --> e- Z - procMap["phot"] = ProcessType::phot; - /// Electron ionization - procMap["eIoni"] = ProcessType::eIoni; - /// Multiple scattering - procMap["msc"] = ProcessType::msc; - /// gamma Z --> Z + X - procMap["photonNuclear"] = ProcessType::photonNuclear; - /// e Z --> e Z + X - procMap["electronNuclear"] = ProcessType::electronNuclear; - /// gamma --> mu+ mu- - procMap["GammaToMuPair"] = ProcessType::GammaToMuPair; - /// e- Z --> e- Z A' - procMap["eDarkBrem"] = ProcessType::eDarkBrem; - return procMap; - } - - SimParticle::ProcessTypeMap SimParticle::PROCESS_MAP = - SimParticle::createProcessTypeMap(); - - SimParticle::SimParticle() {} - - SimParticle::~SimParticle() {} - - void SimParticle::Clear() { - daughters_.clear(); - parents_.clear(); - - energy_ = 0; - pdgID_ = 0; - genStatus_ = -1; - time_ = 0; - x_ = 0; - y_ = 0; - z_ = 0; - endX_ = 0; - endY_ = 0; - endZ_ = 0; - px_ = 0; - py_ = 0; - pz_ = 0; - endpx_ = 0; - endpy_ = 0; - endpz_ = 0; - mass_ = 0; - charge_ = 0; - processType_ = ProcessType::unknown; - vertexVolume_ = ""; - } - - void SimParticle::Print() const { - std::cout << "SimParticle { " - << "energy: " << energy_ << ", " - << "PDG ID: " << pdgID_ << ", " - << "genStatus: " << genStatus_ << ", " - << "time: " << time_ << ", " - << "vertex: ( " << x_ << ", " << y_ << ", " << z_ << " ), " - << "endPoint: ( " << endX_ << ", " << endY_ << ", " << endZ_ - << " ), " - << "momentum: ( " << px_ << ", " << py_ << ", " << pz_ << " ), " - << "endPointMomentum: ( " << endpx_ << ", " << endpy_ << ", " - << endpz_ << " ), " - << "mass: " << mass_ << ", " - << "nDaughters: " << daughters_.size() << ", " - << "nParents: " << parents_.size() << ", " - << "processType: " << processType_ << ", " - << "vertex volume: " << vertexVolume_ << " }" << std::endl; - } - - SimParticle::ProcessType SimParticle::findProcessType( - std::string processName) { - if (processName.find("biasWrapper") != std::string::npos) { - std::size_t pos = processName.find_first_of("(") + 1; - processName = processName.substr(pos, processName.size() - pos - 1); - } - - if (PROCESS_MAP.find(processName) != PROCESS_MAP.end()) { - return PROCESS_MAP[processName]; - } else { - return ProcessType::unknown; - } - } -} // namespace ldmx diff --git a/src/g4fire/Event/SimTrackerHit.cxx b/src/g4fire/Event/SimTrackerHit.cxx deleted file mode 100644 index fd4ac66..0000000 --- a/src/g4fire/Event/SimTrackerHit.cxx +++ /dev/null @@ -1,52 +0,0 @@ -#include "g4fire/Event/SimTrackerHit.h" - -ClassImp(ldmx::SimTrackerHit) - - namespace ldmx { - SimTrackerHit::SimTrackerHit() {} - - SimTrackerHit::~SimTrackerHit() { Clear(); } - - void SimTrackerHit::Print() const { - std::cout << "SimTrackerHit { " - << "id: " << id_ << ", " - << "layerID: " << layerID_ << ", " - << "moduleID: " << moduleID_ << ", " - << "position: ( " << x_ << ", " << y_ << ", " << z_ << " ), " - << "edep: " << edep_ << ", " - << "time: " << time_ << ", " - << "momentum: ( " << px_ << ", " << py_ << ", " << pz_ << " )" - << " }" << std::endl; - } - - void SimTrackerHit::Clear() { - id_ = 0; - layerID_ = 0; - moduleID_ = 0; - edep_ = 0; - time_ = 0; - px_ = 0; - py_ = 0; - pz_ = 0; - x_ = 0; - y_ = 0; - z_ = 0; - energy_ = 0; - pathLength_ = 0; - trackID_ = -1; - pdgID_ = 0; - } - - void SimTrackerHit::setPosition(const float x, const float y, const float z) { - this->x_ = x; - this->y_ = y; - this->z_ = z; - } - - void SimTrackerHit::setMomentum(const float px, const float py, - const float pz) { - this->px_ = px; - this->py_ = py; - this->pz_ = pz; - } -} // namespace ldmx diff --git a/src/g4fire/HcalSD.cxx b/src/g4fire/HcalSD.cxx deleted file mode 100644 index 617b967..0000000 --- a/src/g4fire/HcalSD.cxx +++ /dev/null @@ -1,190 +0,0 @@ -#include "g4fire/HcalSD.h" - -/*~~~~~~~~~~~~~~*/ -/* DetDescr */ -/*~~~~~~~~~~~~~~*/ -#include "DetDescr/HcalID.h" - -// STL -#include - -// Geant4 -#include "G4Box.hh" -#include "G4ChargedGeantino.hh" -#include "G4Geantino.hh" -#include "G4ParticleDefinition.hh" -#include "G4ParticleTypes.hh" -#include "G4SDManager.hh" -#include "G4Step.hh" -#include "G4StepPoint.hh" - -namespace g4fire { - -HcalSD::HcalSD(G4String name, G4String collectionName, int subDetID) - : CalorimeterSD(name, collectionName), - birksc1_(1.29e-2), - birksc2_(9.59e-6) {} - -HcalSD::~HcalSD() {} - -G4bool HcalSD::ProcessHits(G4Step* aStep, G4TouchableHistory* ROhist) { - // Determine if current particle of this step is a Geantino. - G4ParticleDefinition* pdef = aStep->GetTrack()->GetDefinition(); - bool isGeantino = false; - if (pdef == G4Geantino::Definition() || - pdef == G4ChargedGeantino::Definition()) { - isGeantino = true; - } - - // Get the edep from the step. - G4double edep = aStep->GetTotalEnergyDeposit(); - - // Skip steps with no energy dep which come from non-Geantino particles. - if (edep == 0.0 && !isGeantino) { - if (verboseLevel > 2) { - std::cout << "CalorimeterSD skipping step with zero edep." << std::endl - << std::endl; - } - return false; - } - - //--------------------------------------------------------------------------------------------------- - // Birks' Law - // =========== - // - // In the case of Scintillator as active medium, we can - // describe the quenching effects with the Birks' law, - // using the expression and the coefficients taken from - // the paper NIM 80 (1970) 239-244 for the organic - // scintillator NE-102: - // S*dE/dr - // dL/dr = ----------------------------------- - // 1 + C1*(dE/dr) - // with: - // S=1 - // C1 = 1.29 x 10^-2 g*cm^-2*MeV^-1 - // C2 = 9.59 x 10^-6 g^2*cm^-4*MeV^-2 - // These are the same values used by ATLAS TileCal - // and CMS HCAL (and also the default in Geant3). - // You can try different values for these parameters, - // to have an idea on the uncertainties due to them, - // by uncommenting one of the lines below. - // To get the "dE/dr" that appears in the formula, - // which has the dimensions - // [ dE/dr ] = MeV * cm^2 / g - // we have to divide the energy deposit in MeV by the - // product of the step length (in cm) and the density - // of the scintillator: - - G4double birksFactor(1.0); - G4double stepLength = aStep->GetStepLength() / CLHEP::cm; - // Do not apply Birks for gamma deposits! - if (stepLength > 1.0e-6) // Check, cut if necessary. - { - G4double rho = aStep->GetPreStepPoint()->GetMaterial()->GetDensity() / - (CLHEP::g / CLHEP::cm3); - G4double dedx = edep / (rho * stepLength); //[MeV*cm^2/g] - birksFactor = 1.0 / (1.0 + birksc1_ * dedx + birksc2_ * dedx * dedx); - if (aStep->GetTrack()->GetDefinition() == G4Gamma::GammaDefinition()) - birksFactor = 1.0; - if (aStep->GetTrack()->GetDefinition() == G4Neutron::NeutronDefinition()) - birksFactor = 1.0; - } - - // Create a new cal hit. - G4CalorimeterHit* hit = new G4CalorimeterHit(); - - // Set the edep. - hit->setEdep(edep * birksFactor); - - // Get the scintillator solid box - G4Box* scint = static_cast(aStep->GetPreStepPoint() - ->GetTouchableHandle() - ->GetVolume() - ->GetLogicalVolume() - ->GetSolid()); - - // Set the step mid-point as the hit position. - G4StepPoint* prePoint = aStep->GetPreStepPoint(); - G4StepPoint* postPoint = aStep->GetPostStepPoint(); - G4ThreeVector position = - 0.5 * (prePoint->GetPosition() + postPoint->GetPosition()); - G4ThreeVector localPosition = aStep->GetPreStepPoint() - ->GetTouchableHandle() - ->GetHistory() - ->GetTopTransform() - .TransformPoint(position); - hit->setPosition(position[0], position[1], position[2]); - - // Set the global time. - hit->setTime(aStep->GetTrack()->GetGlobalTime()); - - // Create the ID for the hit. - int copyNum = aStep->GetPreStepPoint() - ->GetTouchableHandle() - ->GetHistory() - ->GetVolume(layerDepth_) - ->GetCopyNo(); - int section = copyNum / 1000; - int layer = copyNum % 1000; - - // stripID: back Hcal, segmented along y direction for now every 10 cm -- - // alternate x-y in the future? - // left/right side hcal: segmented along x direction every 10 cm - // top/bottom side hcal: segmented along y direction every 10 cm - - int stripID = -1; - // Odd layers have bars horizontal - // Even layers have bars vertical - // 5cm wide bars are HARD-CODED - if (section == ldmx::HcalID::BACK && layer % 2 == 1) - stripID = int((localPosition.y() + scint->GetYHalfLength()) / 50.0); - else if (section == ldmx::HcalID::BACK && layer % 2 == 0) - stripID = int((localPosition.x() + scint->GetXHalfLength()) / 50.0); - else - stripID = int((localPosition.z() + scint->GetZHalfLength()) / 50.0); - - // std::cout << "---" << std::endl; - // std::cout << "GetXHalfLength = " << scint->GetXHalfLength() << "\t - // GetYHalfLength = " << scint->GetYHalfLength() << "\t GetZHalfLength = " << - // scint->GetZHalfLength() << std::endl; std::cout << "xpos = " << - // localPosition.x() << "\t ypos = " << localPosition.y() << "\t zpos = " << - // localPosition.z() << std::endl; std::cout << "xpos_g = " << position.x() << - // "\t ypos_g = " << position.y() << "\t zpos_g = " << position.z() << - // std::endl; std::cout << "Layer = " << layer << "\t section = " << section - // << "\t strip = " << stripID << std::endl; - - ldmx::HcalID id(section, layer, stripID); - hit->setID(id.raw()); - - // Set the track ID on the hit. - hit->setTrackID(aStep->GetTrack()->GetTrackID()); - - // Set the PDG code from the track. - hit->setPdgCode(aStep->GetTrack()->GetParticleDefinition()->GetPDGEncoding()); - - // do we want to set the hit coordinate in the middle of the absorber? - // G4ThreeVector volumePosition = - // aStep->GetPreStepPoint()->GetTouchableHandle()->GetHistory()->GetTopTransform().Inverse().TransformPoint(G4ThreeVector()); - // if (section==ldmx::HcalID::BACK) hit->setPosition(position[0], position[1], - // volumePosition.z()); elseif (section==ldmx::HcalID::TOP || - // section==ldmx::HcalID::BOTTOM) hit->setPosition(position[0], - // volumePosition.y(), position[2]); elseif (section==ldmx::HcalID::LEFT || - // section==ldmx::HcalID::RIGHT) - // hit->setPosition(volumePosition.x(),position[1] , position[2]); - - if (this->verboseLevel > 2) { - std::cout << "Created new SimCalorimeterHit in detector " << this->GetName() - << " subdet ID <" << subdet_ << ">, layer <" << layer - << "> and section <" << section << ">, copynum <" << copyNum - << ">" << std::endl; - hit->Print(); - std::cout << std::endl; - } - - // Insert the hit into the hits collection. - hitsCollection_->insert(hit); - - return true; -} -} // namespace g4fire diff --git a/src/g4fire/Persist/RootPersistencyManager.cxx b/src/g4fire/Persist/RootPersistencyManager.cxx deleted file mode 100644 index f6fe095..0000000 --- a/src/g4fire/Persist/RootPersistencyManager.cxx +++ /dev/null @@ -1,229 +0,0 @@ - -#include "g4fire/Persist/RootPersistencyManager.h" - -/*~~~~~~~~~~~~~~~~*/ -/* C++ StdLib */ -/*~~~~~~~~~~~~~~~~*/ -#include -#include - -/*~~~~~~~~~~~*/ -/* Event */ -/*~~~~~~~~~~~*/ -#include "Recon/Event/EventConstants.h" - -/*~~~~~~~~~~~~~~~*/ -/* Framework */ -/*~~~~~~~~~~~~~~~*/ -#include "Framework/EventHeader.h" -#include "Framework/RunHeader.h" -#include "Framework/Version.h" - -/*~~~~~~~~~~~~~*/ -/* g4fire */ -/*~~~~~~~~~~~~~*/ -#include "g4fire/DetectorConstruction.h" -#include "g4fire/Event/SimTrackerHit.h" -#include "g4fire/RunManager.h" -#include "g4fire/UserEventInformation.h" -#include "g4fire/UserTrackingAction.h" - -/*~~~~~~~~~~~~*/ -/* Geant4 */ -/*~~~~~~~~~~~~*/ -#include "G4Run.hh" -#include "G4RunManager.hh" -#include "G4RunManagerKernel.hh" - -namespace g4fire { -namespace persist { - -RootPersistencyManager::RootPersistencyManager( - framework::EventFile &file, framework::config::Parameters ¶meters, - const int &runNumber, ConditionsInterface &ci) - : G4PersistencyManager(G4PersistencyCenter::GetPersistencyCenter(), - "RootPersistencyManager"), - file_(file), - ecalHitIO_(ci) { - // Let Geant4 know what to use this persistency manager - G4PersistencyCenter::GetPersistencyCenter()->RegisterPersistencyManager(this); - G4PersistencyCenter::GetPersistencyCenter()->SetPersistencyManager( - this, "RootPersistencyManager"); - - // Set the parameters, used laster when printing run header - parameters_ = parameters; - - ecalHitIO_.configure(parameters_); - - run_ = runNumber; -} - -G4bool RootPersistencyManager::Store(const G4Event *anEvent) { - // Check if the event has been aborted. If so, skip storage of the - // event. - if (G4RunManager::GetRunManager()->GetCurrentEvent()->IsAborted()) - return false; - - // Build the output collections. - buildEvent(anEvent); - - return true; -} - -G4bool RootPersistencyManager::Store(const G4Run *) { - // NOTE: This method is called once the run is terminated through - // the run manager. - - // throws an exception if not correct run number - ldmx::RunHeader &runHeader = file_.getRunHeader(run_); - - // Set parameter value with number of events processed. - runHeader.setIntParameter("Event Count", eventsCompleted_); - runHeader.setIntParameter("Events Began", eventsBegan_); - - return true; -} - -void RootPersistencyManager::Initialize() {} - -void RootPersistencyManager::buildEvent(const G4Event *anEvent) { - // Set basic event information. - writeHeader(anEvent); - - TrackMap* tracks{UserTrackingAction::getUserTrackingAction()->getTrackMap()}; - - tracks->traceAncestry(); - event_->add("SimParticles", tracks->getParticleMap()); - - // Copy hit objects from SD hit collections into the output event. - writeHitsCollections(anEvent, event_); -} - -void RootPersistencyManager::writeHeader(const G4Event *anEvent) { - // Retrieve a mutable version of the event header - ldmx::EventHeader &eventHeader = event_->getEventHeader(); - - auto event_info{ - static_cast(anEvent->GetUserInformation())}; - - eventHeader.setWeight(event_info->getWeight()); - eventHeader.setFloatParameter("total_photonuclear_energy", - event_info->getPNEnergy()); - eventHeader.setFloatParameter("total_electronuclear_energy", - event_info->getENEnergy()); - - // Save the state of the random engine to an output stream. A string - // is then extracted and saved to the event header. - std::ostringstream stream; - G4Random::saveFullState(stream); - // std::cout << stream.str() << std::endl; - eventHeader.setStringParameter("eventSeed", stream.str()); -} - -void RootPersistencyManager::writeHitsCollections( - const G4Event *anEvent, framework::Event *outputEvent) { - // Get the HC of this event. - G4HCofThisEvent *hce = anEvent->GetHCofThisEvent(); - // check for no hit collections for this event - if (!hce) return; - - // Loop over all hits collections. - int nColl = hce->GetNumberOfCollections(); - for (int iColl = 0; iColl < nColl; iColl++) { - // Get a hits collection and its name. - G4VHitsCollection *hc = hce->GetHC(iColl); - if (!hc) { - EXCEPTION_RAISE("G4HitColl", "G4VHitsCollection indexed " + - std::to_string(iColl) + - " returned a nullptr."); - } - - std::string collName = hc->GetName(); - - if (dynamic_cast(hc) != nullptr) { - // Write G4TrackerHit collection to output SimTrackerHit collection. - G4TrackerHitsCollection *trackerHitsColl = - dynamic_cast(hc); - std::vector outputColl; - writeTrackerHitsCollection(trackerHitsColl, outputColl); - - // Add hits collection to output event. - outputEvent->add(collName, outputColl); - - } else if (dynamic_cast(hc) != nullptr) { - G4CalorimeterHitsCollection *calHitsColl = - dynamic_cast(hc); - std::vector outputColl; - if (collName == ldmx::EventConstants::ECAL_SIM_HITS) { - // Write ECal G4CalorimeterHit collection to output SimCalorimeterHit - // collection using helper class. - ecalHitIO_.writeHitsCollection(calHitsColl, outputColl); - } else { - // Write generic G4CalorimeterHit collection to output SimCalorimeterHit - // collection. - writeCalorimeterHitsCollection(calHitsColl, outputColl); - } - - // Add hits collection to output event. - outputEvent->add(collName, outputColl); - } // switch on type of hit collection - - } // loop through geant4 hit collections - - return; -} - -void RootPersistencyManager::writeTrackerHitsCollection( - G4TrackerHitsCollection *hc, std::vector &outputColl) { - outputColl.clear(); - int nHits = hc->GetSize(); - for (int iHit = 0; iHit < nHits; iHit++) { - G4TrackerHit *g4hit = (G4TrackerHit *)hc->GetHit(iHit); - const G4ThreeVector &momentum = g4hit->getMomentum(); - const G4ThreeVector &position = g4hit->getPosition(); - - ldmx::SimTrackerHit simTrackerHit; - simTrackerHit.setID(g4hit->getID()); - simTrackerHit.setTime(g4hit->getTime()); - simTrackerHit.setLayerID(g4hit->getLayerID()); - simTrackerHit.setModuleID(g4hit->getModuleID()); - simTrackerHit.setEdep(g4hit->getEdep()); - simTrackerHit.setEnergy(g4hit->getEnergy()); - simTrackerHit.setPathLength(g4hit->getPathLength()); - simTrackerHit.setTrackID(g4hit->getTrackID()); - simTrackerHit.setPdgID(g4hit->getPdgID()); - simTrackerHit.setPosition(position.x(), position.y(), position.z()); - simTrackerHit.setMomentum(momentum.x(), momentum.y(), momentum.z()); - - outputColl.push_back(simTrackerHit); - } - - return; -} - -void RootPersistencyManager::writeCalorimeterHitsCollection( - G4CalorimeterHitsCollection *hc, - std::vector &outputColl) { - // get ancestral tracking information - auto trackMap{UserTrackingAction::getUserTrackingAction()->getTrackMap()}; - - int nHits = hc->GetSize(); - for (int iHit = 0; iHit < nHits; iHit++) { - G4CalorimeterHit *g4hit = (G4CalorimeterHit *)hc->GetHit(iHit); - const G4ThreeVector &pos = g4hit->getPosition(); - - ldmx::SimCalorimeterHit simHit; - simHit.setID(g4hit->getID()); - simHit.addContrib(trackMap->findIncident(g4hit->getTrackID()), - g4hit->getTrackID(), g4hit->getPdgCode(), - g4hit->getEdep(), g4hit->getTime()); - simHit.setPosition(pos.x(), pos.y(), pos.z()); - - outputColl.push_back(simHit); - } - - return; -} - -} // namespace persist -} // namespace g4fire diff --git a/src/g4fire/RootCompleteReSim.cxx b/src/g4fire/RootCompleteReSim.cxx deleted file mode 100644 index f736b1f..0000000 --- a/src/g4fire/RootCompleteReSim.cxx +++ /dev/null @@ -1,102 +0,0 @@ -/** - * @file RootCompleteReSim.cxx - * @brief Primary generator used to generate primaries from SimParticles. - * @author Nhan Tran, Fermilab - * @author Omar Moreno, SLAC National Accelerator Laboratory - * @author Tom Eichlersmith, University of Minnesota - */ - -#include "g4fire/RootCompleteReSim.h" - -//------------// -// Geant4 // -//------------// -#include "G4Event.hh" -#include "G4IonTable.hh" -#include "G4PhysicalConstants.hh" -#include "G4RunManager.hh" -#include "G4SystemOfUnits.hh" - -//-------------// -// ldmx-sw // -//-------------// -#include "Framework/Configure/Parameters.h" -#include "g4fire/Event/SimParticle.h" - -namespace g4fire { - -RootCompleteReSim::RootCompleteReSim(const std::string& name, - framework::config::Parameters& parameters) - : PrimaryGenerator(name, parameters), ievent_("InputReSim") { - std::string filename = parameters_.getParameter("filePath"); - ifile_ = std::make_unique(parameters, filename); - ifile_->setupEvent(&ievent_); - - simParticleCollName_ = - parameters.getParameter("collection_name"); - simParticlePassName_ = parameters.getParameter("pass_name"); -} - -RootCompleteReSim::~RootCompleteReSim() { - ifile_->close(); - ifile_.reset(nullptr); -} - -void RootCompleteReSim::GeneratePrimaryVertex(G4Event* anEvent) { - // go to next event - // if there isn't another event ==> EventFile::nextEvent returns false - if (not ifile_->nextEvent()) { - std::cout << "[ RootSimFromEcalSP ]: End of file reached." << std::endl; - G4RunManager::GetRunManager()->AbortRun(true); - anEvent->SetEventAborted(); - } - - auto simParticles{ievent_.getMap( - simParticleCollName_, simParticlePassName_)}; - std::vector vertices; // vertices already put into Geant4 - for (const auto& [trackID, sp] : simParticles) { - // check if particle has status 1 - // skip if not primary gen status - if (sp.getGenStatus() != 1) continue; - - bool vertexExists = false; - G4PrimaryVertex* curvertex = new G4PrimaryVertex(); - for (unsigned int iV = 0; iV < vertices.size(); ++iV) { - double cur_vx = sp.getVertex()[0]; - double cur_vy = sp.getVertex()[1]; - double cur_vz = sp.getVertex()[2]; - double i_vx = vertices.at(iV)->GetX0(); - double i_vy = vertices.at(iV)->GetY0(); - double i_vz = vertices.at(iV)->GetZ0(); - if ((cur_vx == i_vx) && (cur_vy == i_vy) && (cur_vz == i_vz)) { - vertexExists = true; - curvertex = vertices.at(iV); - } - } - if (vertexExists == false) { - curvertex->SetPosition(sp.getVertex()[0], sp.getVertex()[1], - sp.getVertex()[2]); - curvertex->SetWeight(1.); - anEvent->AddPrimaryVertex(curvertex); - } - - G4PrimaryParticle* primary = new G4PrimaryParticle(); - primary->SetPDGcode(sp.getPdgID()); - primary->SetMomentum(sp.getMomentum()[0] * MeV, sp.getMomentum()[1] * MeV, - sp.getMomentum()[2] * MeV); - primary->SetMass(sp.getMass() * MeV); - - curvertex->SetPrimary(primary); - } - - // Create an input stream with a copy of the event seed as content. - // The input stream is then passed to the random engine to restore - // the state. - std::istringstream iss( - ievent_.getEventHeader().getStringParameter("eventSeed")); - G4Random::restoreFullState(iss); -} - -} // namespace g4fire - -DECLARE_GENERATOR(g4fire, RootCompleteReSim) diff --git a/src/g4fire/RootSimFromEcalSP.cxx b/src/g4fire/RootSimFromEcalSP.cxx deleted file mode 100644 index dfeaa35..0000000 --- a/src/g4fire/RootSimFromEcalSP.cxx +++ /dev/null @@ -1,126 +0,0 @@ -/** - * @file RootSimFromEcalSP.cxx - * @brief Primary generator used to generate primaries from SimParticles. - * @author Nhan Tran, Fermilab - * @author Omar Moreno, SLAC National Accelerator Laboratory - * @author Tom Eichlersmith, University of Minnesota - */ - -#include "g4fire/RootSimFromEcalSP.h" - -//----------------// -// C++ StdLib // -//----------------// -#include - -//------------// -// Geant4 // -//------------// -#include "G4Event.hh" -#include "G4IonTable.hh" -#include "G4PhysicalConstants.hh" -#include "G4RunManager.hh" -#include "G4SystemOfUnits.hh" - -//-------------// -// ldmx-sw // -//-------------// -#include "Framework/Configure/Parameters.h" -#include "Recon/Event/EventConstants.h" -#include "g4fire/Event/SimTrackerHit.h" - -namespace g4fire { - -RootSimFromEcalSP::RootSimFromEcalSP(const std::string& name, - framework::config::Parameters& parameters) - : PrimaryGenerator(name, parameters), ievent_("InputReSim") { - std::string filename = parameters_.getParameter("filePath"); - ifile_ = std::make_unique(parameters, filename); - ifile_->setupEvent(&ievent_); - - timeCutoff_ = parameters_.getParameter("time_cutoff"); - - ecalSPHitsCollName_ = - parameters_.getParameter("collection_name"); - ecalSPHitsPassName_ = parameters_.getParameter("pass_name"); -} - -RootSimFromEcalSP::~RootSimFromEcalSP() { - ifile_->close(); - ifile_.reset(nullptr); -} - -void RootSimFromEcalSP::GeneratePrimaryVertex(G4Event* anEvent) { - // go to next event - // if there isn't another event ==> EventFile::nextEvent returns false - if (not ifile_->nextEvent()) { - std::cout << "[ RootSimFromEcalSP ]: End of file reached." << std::endl; - G4RunManager::GetRunManager()->AbortRun(true); - anEvent->SetEventAborted(); - } - - // In this mode, we need to loop through all ECal scoring plane hits - // and find the subset of unique hits created by particles exiting - // the ECal. These particles will be stored in a container and - // re-fired into the HCal. - std::unordered_map spHits; - - auto ecalSPParticles{ievent_.getCollection( - ecalSPHitsCollName_, ecalSPHitsPassName_)}; - // Loop through all of the ECal scoring plane hits. - for (const ldmx::SimTrackerHit& spHit : ecalSPParticles) { - // First, start by skipping all hits that were created by - // particles entering the ECal volume. - if (spHit.getLayerID() == 1 and spHit.getMomentum()[2] > 0) continue; - if (spHit.getLayerID() == 2 and spHit.getMomentum()[2] < 0) continue; - if (spHit.getLayerID() == 3 and spHit.getMomentum()[1] < 0) continue; - if (spHit.getLayerID() == 4 and spHit.getMomentum()[1] > 0) continue; - if (spHit.getLayerID() == 5 and spHit.getMomentum()[0] > 0) continue; - if (spHit.getLayerID() == 6 and spHit.getMomentum()[0] < 0) continue; - - // Don't consider particles created outside of the HCal readout - // window. Currently, this is estimated to be 50 ns. - if (spHit.getTime() > timeCutoff_) continue; - - if (spHits.find(spHit.getTrackID()) == spHits.end()) { - spHits[spHit.getTrackID()] = &spHit; - } else { - float currentPMag = - sqrt(pow(spHit.getMomentum()[0], 2) + pow(spHit.getMomentum()[1], 2) + - pow(spHit.getMomentum()[2], 2)); - float pMag = sqrt(pow(spHits[spHit.getTrackID()]->getMomentum()[0], 2) + - pow(spHits[spHit.getTrackID()]->getMomentum()[1], 2) + - pow(spHits[spHit.getTrackID()]->getMomentum()[2], 2)); - - if (pMag < currentPMag) spHits[spHit.getTrackID()] = &spHit; - } - } - - for (auto const& spHit : spHits) { - auto cVertex{new G4PrimaryVertex()}; - cVertex->SetPosition(spHit.second->getPosition()[0] * mm, - spHit.second->getPosition()[1] * mm, - spHit.second->getPosition()[2] * mm); - cVertex->SetWeight(1.); - - auto primary{new G4PrimaryParticle()}; - primary->SetPDGcode(spHit.second->getPdgID()); - primary->SetMomentum(spHit.second->getMomentum()[0] * MeV, - spHit.second->getMomentum()[1] * MeV, - spHit.second->getMomentum()[2] * MeV); - - cVertex->SetPrimary(primary); - anEvent->AddPrimaryVertex(cVertex); - } - - // Create an input stream with a copy of the event seed as content. - // The input stream is then passed to the random engine to restore - // the state. - std::istringstream iss( - ievent_.getEventHeader().getStringParameter("eventSeed")); - G4Random::restoreFullState(iss); -} - -} // namespace g4fire - -DECLARE_GENERATOR(g4fire, RootSimFromEcalSP) diff --git a/src/g4fire/RunManager.cxx b/src/g4fire/RunManager.cxx index e8ab6c7..269f995 100644 --- a/src/g4fire/RunManager.cxx +++ b/src/g4fire/RunManager.cxx @@ -8,18 +8,18 @@ #include "G4VModularPhysicsList.hh" #include "g4fire/ConditionsInterface.h" -#include "g4fire/DarkBrem/APrimePhysics.h" -#include "g4fire/DarkBrem/G4eDarkBremsstrahlung.h" //for process name #include "g4fire/DetectorConstruction.h" #include "g4fire/GammaPhysics.h" #include "g4fire/ParallelWorld.h" #include "g4fire/PluginFactory.h" +#include "g4fire/PrimaryGeneratorAction.h" #include "g4fire/USteppingAction.h" +#include "g4fire/UserEventAction.h" #include "g4fire/UserRunAction.h" #include "g4fire/UserStackingAction.h" #include "g4fire/UserTrackingAction.h" -#include "g4fire/PrimaryGeneratorAction.h" -#include "g4fire/UserEventAction.h" +#include "g4fire/darkbrem/APrimePhysics.h" +#include "g4fire/darkbrem/G4eDarkBremsstrahlung.h" //for process name namespace g4fire { @@ -41,8 +41,8 @@ void RunManager::setupPhysics() { std::cout << "setting up physics." << std::endl; auto physics_list{physics_list_factory_.GetReferencePhysList("FTFP_BERT")}; physics_list->RegisterPhysics(new GammaPhysics); - /*physics_list->RegisterPhysics(new darkbrem::APrimePhysics( - params_.get("dark_brem")));*/ + physics_list->RegisterPhysics(new darkbrem::APrimePhysics( + params_.get("dark_brem"))); parallel_world_path_ = params_.get("parallel_world", {}); pw_enabled_ = !parallel_world_path_.empty(); @@ -51,8 +51,7 @@ void RunManager::setupPhysics() { std::cout << "[ RunManager ]: Parallel worlds physics list has been registered." << std::endl; - physics_list->RegisterPhysics( - new G4ParallelWorldPhysics("parallel_world")); + physics_list->RegisterPhysics(new G4ParallelWorldPhysics("parallel_world")); } auto biasing_operators{params_.get>( @@ -104,7 +103,7 @@ void RunManager::Initialize() { } // This is where the physics lists are told to construct their particles and - // their processes. They are constructed in order, so it is important to + // their processes. They are constructed in order, so it is important to // register the biasing physics *after* any other processes that need to be // able to be biased G4RunManager::Initialize(); @@ -149,7 +148,7 @@ void RunManager::TerminateOneEvent() { // while the table is silenced. If the table isn't silenced, // the process that isn't in the table will cause the table // to throw a "not found" warning. - /*std::vector dark_brem_processes = { + std::vector dark_brem_processes = { darkbrem::G4eDarkBremsstrahlung::PROCESS_NAME, "biasWrapper(" + darkbrem::G4eDarkBremsstrahlung::PROCESS_NAME + ")"}; ptable->SetVerboseLevel( @@ -157,13 +156,12 @@ void RunManager::TerminateOneEvent() { for (auto const &name : dark_brem_processes) ptable->SetProcessActivation(name, true); ptable->SetVerboseLevel(verbosity); - -- Up-to-date: /home/omoreno/projects/ldmx/softwa if (this->GetVerboseLevel() > 1) { std::cout << "[ RunManager ] : " << "Reset the dark brem process (if it was activated)." << std::endl; } - ptable->SetVerboseLevel(verbosity); */ + ptable->SetVerboseLevel(verbosity); } DetectorConstruction *RunManager::getDetectorConstruction() { diff --git a/src/g4fire/ScoringPlaneSD.cxx b/src/g4fire/ScoringPlaneSD.cxx deleted file mode 100644 index 5684a26..0000000 --- a/src/g4fire/ScoringPlaneSD.cxx +++ /dev/null @@ -1,118 +0,0 @@ - -#include "g4fire/ScoringPlaneSD.h" -#include "DetDescr/SimSpecialID.h" - -/*----------------*/ -/* C++ StdLib */ -/*----------------*/ -#include - -/*~~~~~~~~~~~~*/ -/* Geant4 */ -/*~~~~~~~~~~~~*/ -#include "G4ChargedGeantino.hh" -#include "G4Geantino.hh" -#include "G4SDManager.hh" -#include "G4Step.hh" -#include "G4StepPoint.hh" - -namespace g4fire { - -ScoringPlaneSD::ScoringPlaneSD(G4String name, G4String colName, int subDetID) - : G4VSensitiveDetector(name) { - // Add the collection name to vector of names. - this->collectionName.push_back(colName); - - // Register this SD with the manager. - G4SDManager::GetSDMpointer()->AddNewDetector(this); - - // at some point, confirm that the subDetID is as expected... -} - -ScoringPlaneSD::~ScoringPlaneSD() {} - -G4bool ScoringPlaneSD::ProcessHits(G4Step* step, G4TouchableHistory* history) { - // Get the edep from the step. - G4double edep = step->GetTotalEnergyDeposit(); - - // Create a new hit object. - G4TrackerHit* hit = new G4TrackerHit(); - - // Assign track ID for finding the SimParticle in post event processing. - hit->setTrackID(step->GetTrack()->GetTrackID()); - hit->setPdgID(step->GetTrack()->GetDynamicParticle()->GetPDGcode()); - - // Set the edep. - hit->setEdep(edep); - - // Set the start position. - G4StepPoint* prePoint = step->GetPreStepPoint(); - - // Set the end position. - G4StepPoint* postPoint = step->GetPostStepPoint(); - - G4ThreeVector start = prePoint->GetPosition(); - G4ThreeVector end = postPoint->GetPosition(); - - // Set the mid position. - G4ThreeVector mid = 0.5 * (start + end); - hit->setPosition(mid.x(), mid.y(), mid.z()); - - // Compute path length. - G4double pathLength = - sqrt(pow(start.x() - end.x(), 2) + pow(start.y() - end.y(), 2) + - pow(start.z() - end.z(), 2)); - hit->setPathLength(pathLength); - - // Set the global time. - hit->setTime(step->GetTrack()->GetGlobalTime()); - - // Set the momentum - G4ThreeVector p = postPoint->GetMomentum(); - hit->setMomentum(p.x(), p.y(), p.z()); - hit->setEnergy(postPoint->GetTotalEnergy()); - - /* - * Set the 32-bit ID on the hit. - */ - int cpNumber = prePoint->GetTouchableHandle()->GetCopyNumber(); - ldmx::SimSpecialID id = ldmx::SimSpecialID::ScoringPlaneID(cpNumber); - hit->setID(id.raw()); - - /* - * Debug print. - */ - if (this->verboseLevel > 2) { - hit->Print(); - std::cout << std::endl; - } - - // Insert hit into current hits collection. - hitsCollection_->insert(hit); - - return true; -} - -void ScoringPlaneSD::Initialize(G4HCofThisEvent* hce) { - // Setup hits collection and the HC ID. - hitsCollection_ = - new G4TrackerHitsCollection(SensitiveDetectorName, collectionName[0]); - int hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); - hce->AddHitsCollection(hcID, hitsCollection_); -} - -void ScoringPlaneSD::EndOfEvent(G4HCofThisEvent*) { - // Print number of hits. - if (this->verboseLevel > 0) { - std::cout << GetName() << " had " << hitsCollection_->entries() - << " hits in event" << std::endl; - } - - // Print each hit in hits collection. - if (this->verboseLevel > 1) { - for (unsigned iHit = 0; iHit < hitsCollection_->GetSize(); iHit++) { - (*hitsCollection_)[iHit]->Print(); - } - } -} -} // namespace g4fire diff --git a/src/g4fire/Simulator.cxx b/src/g4fire/Simulator.cxx index c7c7bb9..dc5157e 100644 --- a/src/g4fire/Simulator.cxx +++ b/src/g4fire/Simulator.cxx @@ -5,12 +5,12 @@ #include "fire/exception/Exception.h" #include "fire/version/Version.h" -#include "g4fire/DarkBrem/G4eDarkBremsstrahlung.h" +#include "g4fire/darkbrem/G4eDarkBremsstrahlung.h" #include "g4fire/DetectorConstruction.h" #include "g4fire/G4Session.h" -#include "g4fire/Geo/ParserFactory.h" #include "g4fire/PluginFactory.h" #include "g4fire/RunManager.h" +#include "g4fire/geo/ParserFactory.h" #include "G4CascadeParameters.hh" #include "G4Electron.hh" @@ -42,8 +42,6 @@ Simulator::Simulator(const fire::config::Parameters ¶ms) } void Simulator::configure(const fire::config::Parameters ¶ms) { - std::cout << "Configuring ..." << std::endl; - // parameters used to configure the simulation params_ = params; @@ -86,12 +84,11 @@ void Simulator::configure(const fire::config::Parameters ¶ms) { parser->read(); run_manager_->DefineWorldVolume(parser->GetWorldVolume()); - auto pre_init_cmds = - params_.get>("pre_init_cmds", {}); - for (const std::string &cmd : pre_init_cmds) { + auto pre_init_cmds{ + params_.get>("pre_init_cmds", {})}; + for (const auto &cmd : pre_init_cmds) { if (allowed(cmd)) { - int g4ret = ui_manager_->ApplyCommand(cmd); - if (g4ret > 0) { + if (auto g4ret{ui_manager_->ApplyCommand(cmd)}; g4ret > 0) { throw fire::Exception("PreInitCmd", "Pre Initialization command '" + cmd + "' returned a failue status from Geant4: " + @@ -108,19 +105,10 @@ void Simulator::configure(const fire::config::Parameters ¶ms) { } } -/*void Simulator::onFileOpen(fire::EventFile& file) { - // Initialize persistency manager and connect it to the current EventFile - persistencyManager_ = - std::make_unique( - file, params_, this->getRunNumber(), conditions_intf_); - persistencyManager_->Initialize(); -}*/ - void Simulator::beforeNewRun(fire::RunHeader &header) { // Get the detector header from the user detector construction - DetectorConstruction *detector = - static_cast(RunManager::GetRunManager()) - ->getDetectorConstruction(); + auto detector{static_cast(RunManager::GetRunManager()) + ->getDetectorConstruction()}; if (!detector) throw fire::Exception("SimSetup", @@ -135,7 +123,7 @@ void Simulator::beforeNewRun(fire::RunHeader &header) { params_.get("enable_hit_contribs")); header.set("Compress calorimeter hit contribs", params_.get("compress_hit_contribs")); - header.set("Included Scoring Planes", + header.set("Included scoring planes", !params_.get("scoring_planes").empty()); // header.set("Use Random Seed from Event Header", // params_.get("rootPrimaryGenUseSeed")); @@ -150,7 +138,7 @@ void Simulator::beforeNewRun(fire::RunHeader &header) { auto beam_spot_delta{params_.get>("beam_spot_delta", {})}; if (!beam_spot_delta.empty()) - threeVectorDump("Smear Beam Spot [mm]", beam_spot_delta); + threeVectorDump("Beam spot delta [mm]", beam_spot_delta); // lambda function for dumping vectors of strings to the run header auto stringVectorDump = [&header](const std::string &name, @@ -161,17 +149,17 @@ void Simulator::beforeNewRun(fire::RunHeader &header) { } }; - stringVectorDump("Pre Init Command", + stringVectorDump("Pre init Command", params_.get>("pre_init_cmds", {})); - stringVectorDump("Post Init Command", + stringVectorDump("Post init Command", params_.get>("post_init_cmds", {})); auto bops{PluginFactory::getInstance().getBiasingOperators()}; for (const XsecBiasingOperator *bop : bops) { - // bop->RecordConfig(header); + bop->RecordConfig(header); } - /*auto dark_brem{params_.get("dark_brem")}; + auto dark_brem{params_.get("dark_brem")}; if (dark_brem.get("enable")) { // the dark brem process is enabled, find it and then record its // configuration @@ -192,7 +180,7 @@ void Simulator::beforeNewRun(fire::RunHeader &header) { break; } // this process is the dark brem process } // loop through electron processes - } */ // dark brem has been enabled + } // dark brem has been enabled auto generators{ params_.get>("generators")}; @@ -271,10 +259,6 @@ void Simulator::onNewRun(const fire::RunHeader &) { } void Simulator::process(fire::Event &event) { - std::cout << "Processing." << std::endl; - // Pass the current LDMX event object to the persistency manager. This - // is needed by the persistency manager to fill the current event. - // persistencyManager_->setCurrentEvent(&event); // Generate and process a Geant4 event. n_events_began_++; @@ -284,10 +268,14 @@ void Simulator::process(fire::Event &event) { // sequence. This will immediately force the simulation to move on to // the next event. if (run_manager_->GetCurrentEvent()->IsAborted()) { - run_manager_->TerminateOneEvent(); // clean up event objects + run_manager_->TerminateOneEvent(); // clean up event objects this->abortEvent(); // get out of processors loop + } else { + eb_.writeEvent(run_manager_->GetCurrentEvent(), event); } + // Write all hit objects to the event + /*if (this->getLogFrequency() > 0 and event.getEventHeader().getEventNumber() % this->getLogFrequency() == 0) { // print according to log frequency and verbosity @@ -302,22 +290,19 @@ void Simulator::process(fire::Event &event) { n_events_completed_++; run_manager_->TerminateOneEvent(); - return; + return; } void Simulator::onProcessStart() { - std::cout << "on process start" << std::endl; - // initialize run run_manager_->Initialize(); // Get the extra simulation configuring commands auto post_init_cmds{ params_.get>("post_init_cmds", {})}; - for (const std::string &cmd : post_init_cmds) { + for (const auto &cmd : post_init_cmds) { if (allowed(cmd)) { - int g4ret{ui_manager_->ApplyCommand(cmd)}; - if (g4ret > 0) { + if (auto g4ret{ui_manager_->ApplyCommand(cmd)}; g4ret > 0) { throw fire::Exception("PostInitCmd", "Post Initialization command '" + cmd + "' returned a failue status from Geant4: " + @@ -345,27 +330,19 @@ void Simulator::onProcessStart() { return; } -/* -void Simulator::onFileClose(fire::EventFile&) { +void Simulator::onFileClose(const std::string &file_name) { // End the current run and print out some basic statistics if verbose // level > 0. run_manager_->TerminateEventLoop(); // Pass the **real** number of events to the persistency manager - persistencyManager_->setNumEvents(n_events_began_, n_events_completed_); + // persistencyManager_->setNumEvents(n_events_began_, n_events_completed_); // Persist any remaining events, call the end of run action and // terminate the Geant4 kernel. run_manager_->RunTermination(); - - // Cleanup persistency manager - // Geant4 expects us to handle the persistency manager - // In order to avoid segfaulting nonsense, I delete it here - // so that it is deleted before the EventFile it references - // is deleted - persistencyManager_.reset(nullptr); } -*/ + void Simulator::onProcessEnd() { std::cout << "[ Simulator ] : " << "Started " << n_events_began_ << " events to produce " diff --git a/src/g4fire/TrackMap.cxx b/src/g4fire/TrackMap.cxx index 9f021c4..d52b1e8 100644 --- a/src/g4fire/TrackMap.cxx +++ b/src/g4fire/TrackMap.cxx @@ -2,6 +2,7 @@ #include "G4Event.hh" #include "G4EventManager.hh" +#include "G4Track.hh" namespace g4fire { @@ -30,62 +31,62 @@ int TrackMap::findIncident(G4int track_id) const { } void TrackMap::save(const G4Track *track) { - // create sim particle in map, keep reference to the newly created particle - // ldmx::SimParticle& particle{particle_map_[track->GetTrackID()]}; + auto &particle{particle_map_[track->GetTrackID()]}; // TODO: default gen status? - // particle.setGenStatus(0); + particle.setGenStatus(0); // Update the gen status from the primary particle. if (track->GetDynamicParticle()->GetPrimaryParticle() != NULL) { - G4VUserPrimaryParticleInformation *primary_info{track->GetDynamicParticle() - ->GetPrimaryParticle() - ->GetUserInformation()}; + auto primary_info{track->GetDynamicParticle() + ->GetPrimaryParticle() + ->GetUserInformation()}; + if (primary_info != NULL) { - //particle.setGenStatus( - // ((UserPrimaryParticleInformation *)primary_info)->getHepEvtStatus()); + particle.setGenStatus( + static_cast(primary_info) + ->getHepEvtStatus()); } } - /*auto particle_def{track->GetDefinition()}; + auto particle_def{track->GetDefinition()}; particle.setPdgID(particle_def->GetPDGEncoding()); particle.setCharge(particle_def->GetPDGCharge()); particle.setMass(track->GetDynamicParticle()->GetMass()); particle.setEnergy(track->GetVertexKineticEnergy() + track->GetDynamicParticle()->GetMass()); - */ auto track_info{UserTrackInformation::get(track)}; - //particle.setVertexVolume(track_info->getVertexVolume()); + particle.setVertexVolume(track_info->getVertexVolume()); auto vert{track->GetVertexPosition()}; - //particle.setVertex(vert.x(), vert.y(), vert.z()); - //particle.setTime(track_info->getVertexTime()); + particle.setVertex(vert.x(), vert.y(), vert.z()); + particle.setTime(track_info->getVertexTime()); auto init_momentum{track_info->getInitialMomentum()}; - //particle.setMomentum(init_momentum.x(), init_momentum.y(), init_momentum.z()); + particle.setMomentum(init_momentum.x(), init_momentum.y(), init_momentum.z()); const G4VProcess *process{track->GetCreatorProcess()}; if (process) { const G4String &name{process->GetProcessName()}; - //particle.setProcessType(ldmx::SimParticle::findProcessType(name)); + particle.setProcessType(g4fire::event::SimParticle::findProcessType(name)); } else { - //particle.setProcessType(ldmx::SimParticle::ProcessType::unknown); + particle.setProcessType(g4fire::event::SimParticle::ProcessType::unknown); } // track's current kinematics is its end point kinematics // because we are assuming this track is being stopped/killed auto momentum{track->GetMomentum()}; - //particle.setEndPointMomentum(momentum.x(), momentum.y(), momentum.z()); + particle.setEndPointMomentum(momentum.x(), momentum.y(), momentum.z()); auto end_pt{track->GetPosition()}; - //particle.setEndPoint(end_pt.x(), end_pt.y(), end_pt.z()); + particle.setEndPoint(end_pt.x(), end_pt.y(), end_pt.z()); } void TrackMap::traceAncestry() { - //for (auto &[id, particle] : particle_map_) { - // particle.addParent(ancestry_.at(id).first); + for (auto &[id, particle] : particle_map_) { + particle.addParent(ancestry_.at(id).first); /** * Use [] instead of at() for descendents_ @@ -93,16 +94,16 @@ void TrackMap::traceAncestry() { * we will just silently create an empty * vector and move on. */ - // for (auto &child : descendents_[id]) { - // particle.addDaughter(child); - // } - //} + for (auto &child : descendents_[id]) { + particle.addDaughter(child); + } + } } void TrackMap::clear() { ancestry_.clear(); descendents_.clear(); - //particle_map_.clear(); + particle_map_.clear(); } bool TrackMap::isInCalorimeterRegion(const G4Track *track) const { diff --git a/src/g4fire/TrackerSD.cxx b/src/g4fire/TrackerSD.cxx deleted file mode 100644 index 1b1583c..0000000 --- a/src/g4fire/TrackerSD.cxx +++ /dev/null @@ -1,151 +0,0 @@ -#include "g4fire/TrackerSD.h" - -// STL -#include - -// Geant4 -#include "G4ChargedGeantino.hh" -#include "G4Geantino.hh" -#include "G4SDManager.hh" -#include "G4Step.hh" -#include "G4StepPoint.hh" - -// LDMX -#include "DetDescr/IDField.h" - -namespace g4fire { - -TrackerSD::TrackerSD(G4String name, G4String theCollectionName, int subDetID) - : G4VSensitiveDetector(name), hitsCollection_(0) { - // Add the collection name to vector of names. - this->collectionName.push_back(theCollectionName); - - // Register this SD with the manager. - G4SDManager::GetSDMpointer()->AddNewDetector(this); - - // Set the subdet ID as it will always be the same for every hit. - subDetID_ = ldmx::SubdetectorIDType(subDetID); -} - -TrackerSD::~TrackerSD() {} - -G4bool TrackerSD::ProcessHits(G4Step* aStep, G4TouchableHistory*) { - // Determine if current particle of this step is a Geantino. - G4ParticleDefinition* pdef = aStep->GetTrack()->GetDefinition(); - bool isGeantino = false; - if (pdef == G4Geantino::Definition() || - pdef == G4ChargedGeantino::Definition()) { - isGeantino = true; - } - - // Get the edep from the step. - G4double edep = aStep->GetTotalEnergyDeposit(); - - // Skip steps with no energy dep which come from non-Geantino particles. - if (edep == 0.0 && !isGeantino) { - if (verboseLevel > 2) { - std::cout << "TrackerSD skipping step with zero edep" << std::endl - << std::endl; - } - return false; - } - - // Create a new hit object. - G4TrackerHit* hit = new G4TrackerHit(); - - // Assign track ID for finding the SimParticle in post event processing. - hit->setTrackID(aStep->GetTrack()->GetTrackID()); - - // Set the edep. - hit->setEdep(edep); - - // Set the start position. - G4StepPoint* prePoint = aStep->GetPreStepPoint(); - // hit->setStartPosition(prePoint->GetPosition()); - - // Set the end position. - G4StepPoint* postPoint = aStep->GetPostStepPoint(); - // hit->setEndPosition(postPoint->GetPosition()); - - G4ThreeVector start = prePoint->GetPosition(); - G4ThreeVector end = postPoint->GetPosition(); - - // Set the mid position. - G4ThreeVector mid = 0.5 * (start + end); - hit->setPosition(mid.x(), mid.y(), mid.z()); - - // Compute path length. - G4double pathLength = - sqrt(pow(start.x() - end.x(), 2) + pow(start.y() - end.y(), 2) + - pow(start.z() - end.z(), 2)); - hit->setPathLength(pathLength); - - // Set the global time. - hit->setTime(aStep->GetTrack()->GetGlobalTime()); - - /* - * Compute and set the momentum. - */ - /* - double mag = (prePoint->GetMomentum().mag() + postPoint->GetMomentum().mag()) - / 2; G4ThreeVector p = (postPoint->GetPosition() - prePoint->GetPosition()); - if (mag > 0) { - p.setMag(mag); - } - */ - G4ThreeVector p = postPoint->GetMomentum(); - hit->setMomentum(p.x(), p.y(), p.z()); - - /* - * Set the 32-bit ID on the hit. - */ - int copyNum = - prePoint->GetTouchableHandle()->GetHistory()->GetVolume(2)->GetCopyNo(); - int layer = copyNum / 10; - int module = copyNum % 10; - ldmx::TrackerID id(subDetID_, layer, module); - hit->setID(id.raw()); - hit->setLayerID(layer); - hit->setModuleID(module); - - // Set energy and pdg code of SimParticle (common things requested) - hit->setEnergy(postPoint->GetTotalEnergy()); - hit->setPdgID(aStep->GetTrack()->GetDynamicParticle()->GetPDGcode()); - - /* - * Debug print. - */ - if (this->verboseLevel > 2) { - hit->Print(); - } - - // Insert hit into current hits collection. - hitsCollection_->insert(hit); - - return true; -} - -void TrackerSD::Initialize(G4HCofThisEvent* hce) { - // Setup hits collection and the HC ID. - hitsCollection_ = - new G4TrackerHitsCollection(SensitiveDetectorName, collectionName[0]); - int hcID = G4SDManager::GetSDMpointer()->GetCollectionID(collectionName[0]); - hce->AddHitsCollection(hcID, hitsCollection_); -} - -void TrackerSD::EndOfEvent(G4HCofThisEvent*) { - // Print number of hits. - if (this->verboseLevel > 0) { - std::cout << GetName() << " had " << hitsCollection_->entries() - << " hits in event" << std::endl; - } - - // Print each hit in hits collection. - if (this->verboseLevel > 1) { - for (unsigned iHit = 0; iHit < hitsCollection_->GetSize(); iHit++) { - (*hitsCollection_)[iHit]->Print(); - } - } -} - -} // namespace g4fire diff --git a/src/g4fire/TrigScintSD.cxx b/src/g4fire/TrigScintSD.cxx deleted file mode 100644 index 0d324b9..0000000 --- a/src/g4fire/TrigScintSD.cxx +++ /dev/null @@ -1,74 +0,0 @@ -#include "g4fire/TrigScintSD.h" - -/*~~~~~~~~~~~~~~*/ -/* DetDescr */ -/*~~~~~~~~~~~~~~*/ -#include "DetDescr/TrigScintID.h" - -/*~~~~~~~~~~~~*/ -/* Geant4 */ -/*~~~~~~~~~~~~*/ -#include "G4ChargedGeantino.hh" -#include "G4Geantino.hh" -#include "G4SDManager.hh" -#include "G4Step.hh" -#include "G4StepPoint.hh" - -namespace g4fire { - -TrigScintSD::TrigScintSD(G4String name, G4String theCollectionName, - int subDetID) - : CalorimeterSD(name, theCollectionName), moduleId_{subDetID} {} - -TrigScintSD::~TrigScintSD() {} - -G4bool TrigScintSD::ProcessHits(G4Step* step, G4TouchableHistory* history) { - // Get the energy deposited by the particle during the step - auto energy{step->GetTotalEnergyDeposit()}; - - // If a non-Geantino particle doesn't deposit energy during the step, - // skip processing it. - if (auto particleDef{step->GetTrack()->GetDefinition()}; - (energy == 0) && ((particleDef != G4Geantino::Definition()) || - (particleDef != G4ChargedGeantino::Definition()))) - return false; - - // Create a new instance of a calorimeter hit - auto hit{new G4CalorimeterHit()}; - - // Set the energy deposition - hit->setEdep(energy); - - // Set the hit position - auto position{0.5 * (step->GetPreStepPoint()->GetPosition() + - step->GetPostStepPoint()->GetPosition())}; - auto volumePosition{step->GetPreStepPoint() - ->GetTouchableHandle() - ->GetHistory() - ->GetTopTransform() - .Inverse() - .TransformPoint(G4ThreeVector())}; - hit->setPosition(position[0], position[1], volumePosition.z()); - - // Get the track associated with this step - auto track{step->GetTrack()}; - - // Set the global time. - hit->setTime(track->GetGlobalTime()); - - // Set the ID on the hit. - auto bar{track->GetVolume()->GetCopyNo()}; - ldmx::TrigScintID id(moduleId_, bar); - hit->setID(id.raw()); - - // Set the track ID on the hit. - hit->setTrackID(track->GetTrackID()); - - // Set the PDG code from the track. - hit->setPdgCode(track->GetParticleDefinition()->GetPDGEncoding()); - - hitsCollection_->insert(hit); - - return true; -} -} // namespace g4fire diff --git a/src/g4fire/USteppingAction.cxx b/src/g4fire/USteppingAction.cxx index f002afc..84ecc0f 100644 --- a/src/g4fire/USteppingAction.cxx +++ b/src/g4fire/USteppingAction.cxx @@ -8,15 +8,18 @@ void USteppingAction::UserSteppingAction(const G4Step *step) { // get the track weights before this step and after this step // ** these weights include the factors of all upstream step weights ** - double track_weight_pre_step = step->GetPreStepPoint()->GetWeight(); - double track_weight_post_step = step->GetPostStepPoint()->GetWeight(); + auto track_weight_pre_step{step->GetPreStepPoint()->GetWeight()}; + auto track_weight_post_step{step->GetPostStepPoint()->GetWeight()}; // so, to get _this_ step's weight, we divide post_weight by pre_weight - double weight_of_this_step_alone = - track_weight_post_step / track_weight_pre_step; + auto weight_of_this_step_alone{track_weight_post_step / + track_weight_pre_step}; event_info->incWeight(weight_of_this_step_alone); + // + // TODO(OM) This should just move to it's own processor. + // const std::vector *secondaries{ step->GetSecondaryInCurrentStep()}; diff --git a/src/g4fire/UserRegionInformation.cxx b/src/g4fire/UserRegionInformation.cxx index 8f8d5c4..1a3402f 100644 --- a/src/g4fire/UserRegionInformation.cxx +++ b/src/g4fire/UserRegionInformation.cxx @@ -5,8 +5,6 @@ namespace g4fire { UserRegionInformation::UserRegionInformation(bool store_secondaries) : store_secondaries_(store_secondaries) {} -UserRegionInformation::~UserRegionInformation() {} - bool UserRegionInformation::getStoreSecondaries() const { return store_secondaries_; } diff --git a/src/g4fire/UserTrackingAction.cxx b/src/g4fire/UserTrackingAction.cxx index 9f22e7f..d61aabf 100644 --- a/src/g4fire/UserTrackingAction.cxx +++ b/src/g4fire/UserTrackingAction.cxx @@ -12,24 +12,23 @@ namespace g4fire { -void UserTrackingAction::PreUserTrackingAction(const G4Track* track) { +void UserTrackingAction::PreUserTrackingAction(const G4Track *track) { if (!track_map_.contains(track)) { - // New Track - - // get track information and initialize our new track - // this will create a new track info object if it doesn't exist + // If a track doesn't exists within the track map, create it. + // This involves retrieving the track information and instantiating + // a new track object to ecanpsulate that information. + auto track_info{UserTrackInformation::get(track)}; track_info->initialize(track); // Get the region info for where the track was created (could be NULL) - auto region_info = (UserRegionInformation*)track->GetLogicalVolumeAtVertex() - ->GetRegion() - ->GetUserInformation(); + auto region_info{static_cast( + track->GetLogicalVolumeAtVertex()->GetRegion()->GetUserInformation())}; // Get the gen status if track was primary int cur_gen_status = -1; if (track->GetDynamicParticle()->GetPrimaryParticle()) { - auto primaryInfo = dynamic_cast( + auto primaryInfo = dynamic_cast( track->GetDynamicParticle() ->GetPrimaryParticle() ->GetUserInformation()); @@ -44,12 +43,13 @@ void UserTrackingAction::PreUserTrackingAction(const G4Track* track) { * DON'T change the save-status even if these are false * The track's save-status is false by default when the track-info * is constructed and the track's save-status could have been modified by a - * user action **prior** to the track being processed for the first time. + * user action **prior** to the track being processed for the first time. * For example, this happens if the user wants to save the * secondaries of a particular track. */ - if (cur_gen_status == 1 or !region_info or region_info->getStoreSecondaries()) { - track_info->setSaveFlag(true); + if (cur_gen_status == 1 || !region_info || + region_info->getStoreSecondaries()) { + track_info->setSaveFlag(true); } // insert this track into the event's track map @@ -57,13 +57,13 @@ void UserTrackingAction::PreUserTrackingAction(const G4Track* track) { } // Activate user tracking actions - for (auto& tracking_action : tracking_actions_) + for (auto &tracking_action : tracking_actions_) tracking_action->PreUserTrackingAction(track); } -void UserTrackingAction::PostUserTrackingAction(const G4Track* track) { +void UserTrackingAction::PostUserTrackingAction(const G4Track *track) { // Activate user tracking actions - for (auto& tracking_action : tracking_actions_) + for (auto &tracking_action : tracking_actions_) tracking_action->PostUserTrackingAction(track); /** @@ -80,4 +80,4 @@ void UserTrackingAction::PostUserTrackingAction(const G4Track* track) { } } -} // namespace g4fire +} // namespace g4fire diff --git a/src/g4fire/BiasOperators/DarkBrem.cxx b/src/g4fire/biasing/DarkBrem.cxx similarity index 61% rename from src/g4fire/BiasOperators/DarkBrem.cxx rename to src/g4fire/biasing/DarkBrem.cxx index 3136ad8..6910b65 100644 --- a/src/g4fire/BiasOperators/DarkBrem.cxx +++ b/src/g4fire/biasing/DarkBrem.cxx @@ -1,19 +1,16 @@ -#include "g4fire/BiasOperators/DarkBrem.h" +#include "g4fire/biasing/DarkBrem.h" -//------------// -// Geant4 // -//------------// #include "G4RunManager.hh" namespace g4fire { -namespace biasoperators { +namespace biasing { -DarkBrem::DarkBrem(std::string name, const framework::config::Parameters& p) +DarkBrem::DarkBrem(std::string name, const fire::config::Parameters& p) : XsecBiasingOperator(name, p) { - volume_ = p.getParameter("volume"); - factor_ = p.getParameter("factor"); - bias_all_ = p.getParameter("bias_all"); + volume_ = p.get("volume"); + factor_ = p.get("factor"); + bias_all_ = p.get("bias_all"); } G4VBiasingOperation* DarkBrem::ProposeOccurenceBiasingOperation( @@ -41,13 +38,13 @@ G4VBiasingOperation* DarkBrem::ProposeOccurenceBiasingOperation( return 0; } -void DarkBrem::RecordConfig(ldmx::RunHeader& h) const { - h.setIntParameter("BiasOperator::DarkBrem::Bias All Electrons", bias_all_); - h.setFloatParameter("BiasOperator::DarkBrem::Factor", factor_); - h.setStringParameter("BiasOperator::DarkBrem::Volume", volume_); +void DarkBrem::RecordConfig(fire::RunHeader& h) const { + h.set("biasing::DarkBrem::Bias All Electrons", bias_all_); + h.set("biasing::DarkBrem::Factor", factor_); + h.set("biasing::DarkBrem::Volume", volume_); } -} // namespace biasoperators +} // namespace biasing } // namespace g4fire -DECLARE_XSECBIASINGOPERATOR(g4fire::biasoperators, DarkBrem) +DECLARE_XSECBIASINGOPERATOR(g4fire::biasing, DarkBrem) diff --git a/src/g4fire/BiasOperators/ElectroNuclear.cxx b/src/g4fire/biasing/ElectroNuclear.cxx similarity index 71% rename from src/g4fire/BiasOperators/ElectroNuclear.cxx rename to src/g4fire/biasing/ElectroNuclear.cxx index 434ab0d..71affc0 100644 --- a/src/g4fire/BiasOperators/ElectroNuclear.cxx +++ b/src/g4fire/biasing/ElectroNuclear.cxx @@ -1,19 +1,18 @@ -#include "g4fire/BiasOperators/ElectroNuclear.h" +#include "g4fire/biasing/ElectroNuclear.h" -namespace g4fire { -namespace biasoperators { +namespace g4fire::biasing { ElectroNuclear::ElectroNuclear(std::string name, - const framework::config::Parameters& p) + const fire::config::Parameters &p) : XsecBiasingOperator(name, p) { - volume_ = p.getParameter("volume"); - factor_ = p.getParameter("factor"); - threshold_ = p.getParameter("threshold"); + volume_ = p.get("volume"); + factor_ = p.get("factor"); + threshold_ = p.get("threshold"); } -G4VBiasingOperation* ElectroNuclear::ProposeOccurenceBiasingOperation( - const G4Track* track, const G4BiasingProcessInterface* callingProcess) { +G4VBiasingOperation *ElectroNuclear::ProposeOccurenceBiasingOperation( + const G4Track *track, const G4BiasingProcessInterface *callingProcess) { /*std::cout << "[ ElectroNuclearXsecBiasingOperator ]: " << "Track ID: " << track->GetTrackID() << ", " << "Parent ID: " << track->GetParentID() << ", " @@ -22,7 +21,8 @@ G4VBiasingOperation* ElectroNuclear::ProposeOccurenceBiasingOperation( << "Currently in volume " << track->GetVolume()->GetName() << std::endl;*/ - if (track->GetKineticEnergy() < threshold_) return 0; + if (track->GetKineticEnergy() < threshold_) + return 0; /*std::cout << "[ ElectroNuclearXsecBiasingOperator ]: " << "Calling process: " @@ -51,7 +51,6 @@ G4VBiasingOperation* ElectroNuclear::ProposeOccurenceBiasingOperation( return 0; } -} // namespace biasoperators -} // namespace g4fire +} // namespace g4fire::biasing -DECLARE_XSECBIASINGOPERATOR(g4fire::biasoperators, ElectroNuclear) +DECLARE_XSECBIASINGOPERATOR(g4fire::biasing, ElectroNuclear) diff --git a/src/g4fire/BiasOperators/GammaToMuPair.cxx b/src/g4fire/biasing/GammaToMuPair.cxx similarity index 64% rename from src/g4fire/BiasOperators/GammaToMuPair.cxx rename to src/g4fire/biasing/GammaToMuPair.cxx index 17b3bcd..df45f8f 100644 --- a/src/g4fire/BiasOperators/GammaToMuPair.cxx +++ b/src/g4fire/biasing/GammaToMuPair.cxx @@ -1,15 +1,14 @@ -#include "g4fire/BiasOperators/GammaToMuPair.h" +#include "g4fire/biasing/GammaToMuPair.h" -namespace g4fire { -namespace biasoperators { +namespace g4fire::biasing { GammaToMuPair::GammaToMuPair(std::string name, - const framework::config::Parameters& p) + const fire::config::Parameters& p) : XsecBiasingOperator(name, p) { - volume_ = p.getParameter("volume"); - factor_ = p.getParameter("factor"); - threshold_ = p.getParameter("threshold"); + volume_ = p.get("volume"); + factor_ = p.get("factor"); + threshold_ = p.get("threshold"); } G4VBiasingOperation* GammaToMuPair::ProposeOccurenceBiasingOperation( @@ -31,7 +30,6 @@ G4VBiasingOperation* GammaToMuPair::ProposeOccurenceBiasingOperation( return 0; } -} // namespace biasoperators } // namespace g4fire -DECLARE_XSECBIASINGOPERATOR(g4fire::biasoperators, GammaToMuPair) +DECLARE_XSECBIASINGOPERATOR(g4fire::biasing, GammaToMuPair) diff --git a/src/g4fire/BiasOperators/K0LongInelastic.cxx b/src/g4fire/biasing/K0LongInelastic.cxx similarity index 62% rename from src/g4fire/BiasOperators/K0LongInelastic.cxx rename to src/g4fire/biasing/K0LongInelastic.cxx index a5f2538..02afc14 100644 --- a/src/g4fire/BiasOperators/K0LongInelastic.cxx +++ b/src/g4fire/biasing/K0LongInelastic.cxx @@ -1,14 +1,13 @@ -#include "g4fire/BiasOperators/K0LongInelastic.h" +#include "g4fire/biasing/K0LongInelastic.h" -namespace g4fire { -namespace biasoperators { +namespace g4fire::biasing { -K0LongInelastic::K0LongInelastic(std::string name, const framework::config::Parameters& p) +K0LongInelastic::K0LongInelastic(std::string name, const fire::config::Parameters& p) : XsecBiasingOperator(name, p) { - volume_ = p.getParameter("volume"); - factor_ = p.getParameter("factor"); - threshold_ = p.getParameter("threshold"); + volume_ = p.get("volume"); + factor_ = p.get("factor"); + threshold_ = p.get("threshold"); } G4VBiasingOperation* K0LongInelastic::ProposeOccurenceBiasingOperation( @@ -31,7 +30,6 @@ G4VBiasingOperation* K0LongInelastic::ProposeOccurenceBiasingOperation( return 0; } -} // namespace biasoperators } // namespace g4fire -DECLARE_XSECBIASINGOPERATOR(g4fire::biasoperators, K0LongInelastic) +DECLARE_XSECBIASINGOPERATOR(g4fire::biasing, K0LongInelastic) diff --git a/src/g4fire/biasing/NeutronInelastic.cxx b/src/g4fire/biasing/NeutronInelastic.cxx new file mode 100644 index 0000000..6f9dde0 --- /dev/null +++ b/src/g4fire/biasing/NeutronInelastic.cxx @@ -0,0 +1,37 @@ + +#include "g4fire/biasing/NeutronInelastic.h" + +namespace g4fire::biasing { + +NeutronInelastic::NeutronInelastic(std::string name, + const fire::config::Parameters &p) + : XsecBiasingOperator(name, p) { + volume_ = p.get("volume"); + factor_ = p.get("factor"); + threshold_ = p.get("threshold"); +} + +G4VBiasingOperation *NeutronInelastic::ProposeOccurenceBiasingOperation( + const G4Track *track, const G4BiasingProcessInterface *callingProcess) { + + if (track->GetKineticEnergy() < threshold_) + return 0; + + std::string currentProcess = + callingProcess->GetWrappedProcess()->GetProcessName(); + if (currentProcess.compare(this->getProcessToBias()) == 0) { + G4double interactionLength = + callingProcess->GetWrappedProcess()->GetCurrentInteractionLength(); + + double neutInXsecUnbiased = 1. / interactionLength; + + double neutInXsecBiased = neutInXsecUnbiased * factor_; + + return BiasedXsec(neutInXsecBiased); + } else + return 0; +} + +} // namespace g4fire::biasing + +DECLARE_XSECBIASINGOPERATOR(g4fire::biasing, NeutronInelastic) diff --git a/src/g4fire/BiasOperators/PhotoNuclear.cxx b/src/g4fire/biasing/PhotoNuclear.cxx similarity index 73% rename from src/g4fire/BiasOperators/PhotoNuclear.cxx rename to src/g4fire/biasing/PhotoNuclear.cxx index 87f38c7..8172ca1 100644 --- a/src/g4fire/BiasOperators/PhotoNuclear.cxx +++ b/src/g4fire/biasing/PhotoNuclear.cxx @@ -1,18 +1,19 @@ -#include "g4fire/BiasOperators/PhotoNuclear.h" -namespace g4fire { -namespace biasoperators { +#include "g4fire/biasing/PhotoNuclear.h" + +#include "fire/exception/Exception.h" + +namespace g4fire::biasing { const std::string PhotoNuclear::CONVERSION_PROCESS = "conv"; -PhotoNuclear::PhotoNuclear(std::string name, - const framework::config::Parameters& p) +PhotoNuclear::PhotoNuclear(std::string name, const fire::config::Parameters &p) : XsecBiasingOperator(name, p) { - volume_ = p.getParameter("volume"); - threshold_ = p.getParameter("threshold"); - factor_ = p.getParameter("factor"); - down_bias_conv_ = p.getParameter("down_bias_conv"); - only_children_of_primary_ = p.getParameter("only_children_of_primary"); + volume_ = p.get("volume"); + threshold_ = p.get("threshold"); + factor_ = p.get("factor"); + down_bias_conv_ = p.get("down_bias_conv"); + only_children_of_primary_ = p.get("only_children_of_primary"); } void PhotoNuclear::StartRun() { @@ -21,24 +22,27 @@ void PhotoNuclear::StartRun() { if (processIsBiased(CONVERSION_PROCESS)) { emXsecOperation = new G4BOptnChangeCrossSection("changeXsec-conv"); } else if (down_bias_conv_) { - EXCEPTION_RAISE( - "PhotoNuclearBiasing", - "Gamma Conversion process '" + CONVERSION_PROCESS + "' is not biased!"); + throw fire::Exception("PhotoNuclearBiasing", + "Gamma Conversion process '" + CONVERSION_PROCESS + + "' is not biased!", + false); } } -G4VBiasingOperation* PhotoNuclear::ProposeOccurenceBiasingOperation( - const G4Track* track, const G4BiasingProcessInterface* callingProcess) { +G4VBiasingOperation *PhotoNuclear::ProposeOccurenceBiasingOperation( + const G4Track *track, const G4BiasingProcessInterface *callingProcess) { /*std::cout << "[ PhotoNuclearXsecBiasingOperator ]: " << "Kinetic energy: " << track->GetKineticEnergy() << " MeV" << std::endl;*/ // if we want to only bias children of primary, leave if this track is NOT a // child of the primary - if (only_children_of_primary_ and track->GetParentID() != 1) return 0; + if (only_children_of_primary_ and track->GetParentID() != 1) + return 0; // is this track too low energy to be biased? - if (track->GetKineticEnergy() < threshold_) return 0; + if (track->GetKineticEnergy() < threshold_) + return 0; /*std::cout << "[ PhotoNuclearXsecBiasingOperator ]: " << "Calling process: " @@ -93,8 +97,6 @@ G4VBiasingOperation* PhotoNuclear::ProposeOccurenceBiasingOperation( } else return 0; } +} // namespace g4fire::biasing -} // namespace biasoperators -} // namespace g4fire - -DECLARE_XSECBIASINGOPERATOR(g4fire::biasoperators, PhotoNuclear) +DECLARE_XSECBIASINGOPERATOR(g4fire::biasing, PhotoNuclear) diff --git a/src/g4fire/DarkBrem/APrimePhysics.cxx b/src/g4fire/darkbrem/APrimePhysics.cxx similarity index 95% rename from src/g4fire/DarkBrem/APrimePhysics.cxx rename to src/g4fire/darkbrem/APrimePhysics.cxx index bd6f199..7f1ef8c 100644 --- a/src/g4fire/DarkBrem/APrimePhysics.cxx +++ b/src/g4fire/darkbrem/APrimePhysics.cxx @@ -1,6 +1,6 @@ -#include "g4fire/DarkBrem/APrimePhysics.h" +#include "g4fire/darkbrem/APrimePhysics.h" -#include "g4fire/DarkBrem/G4APrime.h" +#include "g4fire/darkbrem/G4APrime.h" #include "G4Electron.hh" #include "G4ProcessManager.hh" diff --git a/src/g4fire/DarkBrem/DarkBremVertexLibraryModel.cxx b/src/g4fire/darkbrem/DarkBremVertexLibraryModel.cxx similarity index 99% rename from src/g4fire/DarkBrem/DarkBremVertexLibraryModel.cxx rename to src/g4fire/darkbrem/DarkBremVertexLibraryModel.cxx index d3cd8cd..68ef660 100644 --- a/src/g4fire/DarkBrem/DarkBremVertexLibraryModel.cxx +++ b/src/g4fire/darkbrem/DarkBremVertexLibraryModel.cxx @@ -1,4 +1,4 @@ -#include "g4fire/DarkBrem/DarkBremVertexLibraryModel.h" +#include "g4fire/darkbrem/DarkBremVertexLibraryModel.h" #include #include @@ -11,7 +11,7 @@ #include "fire/exception/Exception.h" //#include "Framework/Logger.h" -#include "g4fire/DarkBrem/G4APrime.h" +#include "g4fire/darkbrem/G4APrime.h" #include "G4Electron.hh" #include "G4EventManager.hh" diff --git a/src/g4fire/DarkBrem/G4APrime.cxx b/src/g4fire/darkbrem/G4APrime.cxx similarity index 96% rename from src/g4fire/DarkBrem/G4APrime.cxx rename to src/g4fire/darkbrem/G4APrime.cxx index a1326d0..aeb99b2 100644 --- a/src/g4fire/DarkBrem/G4APrime.cxx +++ b/src/g4fire/darkbrem/G4APrime.cxx @@ -1,4 +1,4 @@ -#include "g4fire/DarkBrem/G4APrime.h" +#include "g4fire/darkbrem/G4APrime.h" //#include "Framework/Exception/Exception.h" #include "G4ParticleTable.hh" diff --git a/src/g4fire/DarkBrem/G4eDarkBremsstrahlung.cxx b/src/g4fire/darkbrem/G4eDarkBremsstrahlung.cxx similarity index 97% rename from src/g4fire/DarkBrem/G4eDarkBremsstrahlung.cxx rename to src/g4fire/darkbrem/G4eDarkBremsstrahlung.cxx index 9f1d571..da35087 100644 --- a/src/g4fire/DarkBrem/G4eDarkBremsstrahlung.cxx +++ b/src/g4fire/darkbrem/G4eDarkBremsstrahlung.cxx @@ -1,4 +1,4 @@ -#include "g4fire/DarkBrem/G4eDarkBremsstrahlung.h" +#include "g4fire/darkbrem/G4eDarkBremsstrahlung.h" #include "fire/RunHeader.h" @@ -8,8 +8,8 @@ #include "G4ProcessType.hh" //for type of process #include "G4RunManager.hh" //for VerboseLevel -#include "g4fire/DarkBrem/DarkBremVertexLibraryModel.h" -#include "g4fire/DarkBrem/G4APrime.h" +#include "g4fire/darkbrem/DarkBremVertexLibraryModel.h" +#include "g4fire/darkbrem/G4APrime.h" namespace g4fire { namespace darkbrem { diff --git a/src/g4fire/DarkBrem/print_dark_brem_xsec_table.cxx.in b/src/g4fire/darkbrem/print_dark_brem_xsec_table.cxx.in similarity index 96% rename from src/g4fire/DarkBrem/print_dark_brem_xsec_table.cxx.in rename to src/g4fire/darkbrem/print_dark_brem_xsec_table.cxx.in index acfbb67..8012f3b 100644 --- a/src/g4fire/DarkBrem/print_dark_brem_xsec_table.cxx.in +++ b/src/g4fire/darkbrem/print_dark_brem_xsec_table.cxx.in @@ -9,7 +9,7 @@ // LDMX // //----------// #include "Framework/Configure/Parameters.h" -#include "g4fire/DarkBrem/G4eDarkBremsstrahlung.h" +#include "g4fire/darkbrem/G4eDarkBremsstrahlung.h" /** * @func printUsage diff --git a/src/g4fire/event/EventBuilder.cxx b/src/g4fire/event/EventBuilder.cxx new file mode 100644 index 0000000..683f903 --- /dev/null +++ b/src/g4fire/event/EventBuilder.cxx @@ -0,0 +1,39 @@ +#include "g4fire/event/EventBuilder.h" + +#include + +#include "G4Event.hh" +#include "G4RunManager.hh" + +#include "fire/EventHeader.h" + +#include "g4fire/UserTrackingAction.h" + +namespace g4fire::event { + +void EventBuilder::writeEvent(const G4Event *g4event, fire::Event &event) { + writeHeader(g4event, event); + + auto tracks{UserTrackingAction::getUserTrackingAction()->getTrackMap()}; + + tracks->traceAncestry(); + event.add("SimParticles", tracks->particleMap()); +} + +void EventBuilder::writeHeader(const G4Event *g4event, fire::Event &event) { + + // Retrieve a mutable version of the event header + fire::EventHeader &header{event.header()}; + + // auto event_info{ + // static_cast(g4event->GetUserInformation())}; + + // Save the state of the random engine to an output stream. A string + // is then extracted and saved to the event header. + std::ostringstream stream; + G4Random::saveFullState(stream); + // std::cout << stream.str() << std::endl; + header.set("event_seed", stream.str()); +} + +} // namespace g4fire::event diff --git a/src/g4fire/event/SimCalorimeterHit.cxx b/src/g4fire/event/SimCalorimeterHit.cxx new file mode 100644 index 0000000..04c42b4 --- /dev/null +++ b/src/g4fire/event/SimCalorimeterHit.cxx @@ -0,0 +1,96 @@ +#include "g4fire/event/SimCalorimeterHit.h" + +namespace g4fire::event { + +void SimCalorimeterHit::clear() { + + contribs_incident_id_.clear(); + contribs_track_id_.clear(); + contribs_pdg_id_.clear(); + contribs_edeps_.clear(); + contribs_time_.clear(); + + contrib_count_ = 0; + id_ = 0; + edep_ = 0; + x_ = 0; + y_ = 0; + z_ = 0; + time_ = 0; +} + +void SimCalorimeterHit::attach(fire::io::Data &d) { + d.attach("contrib_count", contrib_count_); + d.attach("id", id_); + d.attach("edep", edep_); + d.attach("x", x_); + d.attach("y", y_); + d.attach("z", z_); + d.attach("time", time_); + d.attach("contribs_incident_id", contribs_incident_id_); + d.attach("contribs_track_id", contribs_track_id_); + d.attach("contribs_pdg_id", contribs_pdg_id_); + d.attach("contribs_edeps", contribs_edeps_); + d.attach("contribs_time", contribs_time_); +} + +void SimCalorimeterHit::addContrib(const int &incident_id, const int &track_id, + const int &pdg_id, const float &edep, + const float &time) { + contribs_incident_id_.push_back(incident_id); + contribs_track_id_.push_back(track_id); + contribs_pdg_id_.push_back(pdg_id); + contribs_edeps_.push_back(edep); + contribs_time_.push_back(time); + + edep_ += edep; + + if (time < time_ || time_ == 0) { + time_ = time; + } + ++contrib_count_; +} + +SimCalorimeterHit::Contrib SimCalorimeterHit::getContrib(const int &i) const { + // TODO(OM): Should these be in a vector so they don't need to be created + // everytime getContrib is called? This will also simplify searching. + Contrib contrib; + contrib.incident_id = contribs_incident_id_.at(i); + contrib.track_id = contribs_track_id_.at(i); + contrib.edep = contribs_edeps_.at(i); + contrib.time = contribs_time_.at(i); + contrib.pdg_id = contribs_pdg_id_.at(i); + + return contrib; +} + +int SimCalorimeterHit::findContribIndex(const int& track_id, const int& pdg_id) const { + for (int jcontrib{0}; jcontrib < contrib_count_; ++jcontrib) { + Contrib contrib{getContrib(jcontrib)}; + if (contrib.track_id == track_id && contrib.pdg_id == pdg_id) { + return jcontrib; + } + } + return -1; +} + +void SimCalorimeterHit::updateContrib(const int& i, const float& edep, const float& time) { + contribs_edeps_[i] += edep; + if (time < contribs_time_.at(i)) { + contribs_time_[i] = time; + } + edep_ += edep; +} + +std::ostream &operator<<(std::ostream &output, const SimCalorimeterHit &hit) { + output << "---[ SimCalorimeterHit ] { \n" + << "\tID: " << hit.id() << "\n" + << "\tEnergy Deposition (MeV): " << hit.edep() << "\n" + << "\tPosition (mm): ( " << hit.position()[0] << ", " + << hit.position()[1] << ", " << hit.position()[2] << " )\n" + << "\tContribution count: " << hit.contribCount() << "\n" + << "}" << std::endl; + + return output; +} +} // namespace g4fire::event diff --git a/src/g4fire/event/SimParticle.cxx b/src/g4fire/event/SimParticle.cxx new file mode 100644 index 0000000..d258c10 --- /dev/null +++ b/src/g4fire/event/SimParticle.cxx @@ -0,0 +1,126 @@ +#include "g4fire/event/SimParticle.h" + +#include + +namespace g4fire::event { + +SimParticle::ProcessTypeMap SimParticle::createProcessTypeMap() { + + ProcessTypeMap procMap; + /// e Z --> e Z gamma + procMap["eBrem"] = ProcessType::eBrem; + /// gamma --> e+ e- + procMap["conv"] = ProcessType::conv; + /// e+ e- --> gamma gamma + procMap["annihil"] = ProcessType::annihil; + /// gamma e --> gamma e + procMap["compt"] = ProcessType::compt; + /// gamma Z --> e- Z + procMap["phot"] = ProcessType::phot; + /// Electron ionization + procMap["eIoni"] = ProcessType::eIoni; + /// Multiple scattering + procMap["msc"] = ProcessType::msc; + /// gamma Z --> Z + X + procMap["photonNuclear"] = ProcessType::photonNuclear; + /// e Z --> e Z + X + procMap["electronNuclear"] = ProcessType::electronNuclear; + /// gamma --> mu+ mu- + procMap["GammaToMuPair"] = ProcessType::GammaToMuPair; + /// e- Z --> e- Z A' + procMap["eDarkBrem"] = ProcessType::eDarkBrem; + return procMap; +} + +SimParticle::ProcessTypeMap SimParticle::PROCESS_MAP = + SimParticle::createProcessTypeMap(); + +void SimParticle::clear() { + daughters_.clear(); + parents_.clear(); + + energy_ = 0; + pdg_id_ = 0; + gen_status_ = -1; + time_ = 0; + x_ = 0; + y_ = 0; + z_ = 0; + endx_ = 0; + endy_ = 0; + endz_ = 0; + px_ = 0; + py_ = 0; + pz_ = 0; + endpx_ = 0; + endpy_ = 0; + endpz_ = 0; + mass_ = 0; + charge_ = 0; + process_type_ = ProcessType::unknown; + vertex_vol_ = ""; +} + +void SimParticle::attach(fire::io::Data &d) { + d.attach("energy", energy_); + d.attach("pdg_id", pdg_id_); + d.attach("gen_status", gen_status_); + d.attach("time", time_); + d.attach("x", x_); + d.attach("y", y_); + d.attach("z", z_); + d.attach("endx", endx_); + d.attach("endy", endy_); + d.attach("endz", endz_); + d.attach("px", px_); + d.attach("py", py_); + d.attach("pz", pz_); + d.attach("endpx", endpx_); + d.attach("endpy", endpy_); + d.attach("endpz", endpz_); + d.attach("mass", mass_); + d.attach("charge", charge_); + d.attach("process_type", process_type_); + d.attach("vertex_vol", vertex_vol_); + d.attach("daughters", daughters_); + d.attach("parents", parents_); +} + +std::ostream &operator<<(std::ostream &output, const SimParticle &particle) { + output << "---[ SimParticle ] { \n" + << "\tEnergy (MeV): " << particle.energy() << "\n" + << "\tPDG ID: " << particle.pdgID() << "\n" + << "\tGenerator Status: " << particle.genStatus() << "\n" + << "\tTime (ns): " << particle.time() << "\n" + << "\tVertex (mm): ( " << particle.vertex()[0] << ", " + << particle.vertex()[1] << ", " << particle.vertex()[2] << " )\n" + << "\tEndpoint (mm): ( " << particle.endPoint()[0] << ", " + << particle.endPoint()[1] << ", " << particle.endPoint()[2] << " )\n" + << "\tMomentum (MeV): ( " << particle.momentum()[0] << ", " + << particle.momentum()[1] << ", " << particle.momentum()[2] << " )\n" + << "\tEndpoint Momentum (MeV): ( " << particle.endPointMomentum()[0] + << ", " << particle.endPointMomentum()[1] << ", " + << particle.endPointMomentum()[2] << " )\n" + << "\tMass (MeV): " << particle.mass() << "\n" + << "\tDaughter count: " << particle.daughters().size() << "\n" + << "\tParent count: " << particle.parents().size() << "\n" + << "\tProcess type: " << particle.processType() << "\n" + << "\tVertex volume: " << particle.vertexVolume() << "\n}" + << std::endl; + + return output; +} + +SimParticle::ProcessType SimParticle::findProcessType(std::string processName) { + if (processName.find("biasWrapper") != std::string::npos) { + std::size_t pos = processName.find_first_of("(") + 1; + processName = processName.substr(pos, processName.size() - pos - 1); + } + + if (PROCESS_MAP.find(processName) != PROCESS_MAP.end()) { + return PROCESS_MAP[processName]; + } else { + return ProcessType::unknown; + } +} +} // namespace g4fire::event diff --git a/src/g4fire/event/SimTrackerHit.cxx b/src/g4fire/event/SimTrackerHit.cxx new file mode 100644 index 0000000..380868b --- /dev/null +++ b/src/g4fire/event/SimTrackerHit.cxx @@ -0,0 +1,60 @@ +#include "g4fire/event/SimTrackerHit.h" + +namespace g4fire::event { + +void SimTrackerHit::clear() { + id_ = 0; + layer_id_ = 0; + module_id_ = 0; + edep_ = 0; + time_ = 0; + px_ = 0; + py_ = 0; + pz_ = 0; + x_ = 0; + y_ = 0; + z_ = 0; + energy_ = 0; + path_length_ = 0; + track_id_ = -1; + pdg_id_ = 0; +} + +void SimTrackerHit::attach(fire::io::Data &d) { + d.attach("id", id_); + d.attach("layer_id", layer_id_); + d.attach("module_id", module_id_); + d.attach("edep", edep_); + d.attach("time", time_); + d.attach("x", x_); + d.attach("y", y_); + d.attach("z", z_); + d.attach("px", px_); + d.attach("py", py_); + d.attach("pz", pz_); + d.attach("energy", energy_); + d.attach("path_length", path_length_); + d.attach("track_id", track_id_); + d.attach("pdg_id", pdg_id_); +} + +std::ostream &operator<<(std::ostream &output, const SimTrackerHit &hit) { + output << "---[ SimTrackerHit ] {\n" + << "\tID: " << hit.id() << "\n" + << "\tLayer ID: " << hit.layerID() << "\n" + << "\tModule ID: " << hit.moduleID() << "\n" + << "\tPosition (mm): ( " << hit.position()[0] << ", " + << hit.position()[1] << ", " << hit.position()[2] << ")\n" + << "\tEnergy Deposited (MeV): " << hit.edep() << "\n" + << "\tEnergy (MeV): " << hit.energy() << "\n" + << "\tTime (ns): " << hit.time() << "\n" + << "\tMomentum (MeV): ( " << hit.momentum()[0] << ", " + << hit.momentum()[1] << ", " << hit.momentum()[2] << " )\n" + << "\tPath length: " << hit.pathLength() << "\n" + << "\tTrack ID: " << hit.trackID() << "\n" + << "\tPDG ID: " << hit.pdgID() << "\n" + << " }" << std::endl; + + return output; +} +} // namespace g4fire::event diff --git a/src/g4fire/Geo/AuxInfoReader.cxx b/src/g4fire/geo/AuxInfoReader.cxx similarity index 99% rename from src/g4fire/Geo/AuxInfoReader.cxx rename to src/g4fire/geo/AuxInfoReader.cxx index f311d96..40d65f8 100644 --- a/src/g4fire/Geo/AuxInfoReader.cxx +++ b/src/g4fire/geo/AuxInfoReader.cxx @@ -1,4 +1,4 @@ -#include "g4fire/Geo/AuxInfoReader.h" +#include "g4fire/geo/AuxInfoReader.h" #include #include diff --git a/src/g4fire/Geo/GDMLParser.cxx b/src/g4fire/geo/GDMLParser.cxx similarity index 95% rename from src/g4fire/Geo/GDMLParser.cxx rename to src/g4fire/geo/GDMLParser.cxx index 52b6efd..463ee00 100644 --- a/src/g4fire/Geo/GDMLParser.cxx +++ b/src/g4fire/geo/GDMLParser.cxx @@ -1,4 +1,4 @@ -#include "g4fire/Geo/GDMLParser.h" +#include "g4fire/geo/GDMLParser.h" namespace g4fire { namespace geo { diff --git a/src/g4fire/Geo/ParserFactory.cxx b/src/g4fire/geo/ParserFactory.cxx similarity index 88% rename from src/g4fire/Geo/ParserFactory.cxx rename to src/g4fire/geo/ParserFactory.cxx index f39e328..8fb5077 100644 --- a/src/g4fire/Geo/ParserFactory.cxx +++ b/src/g4fire/geo/ParserFactory.cxx @@ -1,10 +1,10 @@ -#include "g4fire/Geo/ParserFactory.h" +#include "g4fire/geo/ParserFactory.h" #include "fire/exception/Exception.h" -#include "g4fire/Geo/GDMLParser.h" -#include "g4fire/Geo/Parser.h" +#include "g4fire/geo/GDMLParser.h" +#include "g4fire/geo/Parser.h" namespace g4fire { namespace geo {