diff --git a/inc/TRestGeant4Event.h b/inc/TRestGeant4Event.h index 1e5b450..1da6782 100644 --- a/inc/TRestGeant4Event.h +++ b/inc/TRestGeant4Event.h @@ -29,15 +29,18 @@ #include #include #include -#include #include -#include #include #include #include #include +#include "TRestGeant4Track.h" + +class G4Event; +class G4Track; +class G4Step; class TRestGeant4Metadata; /// An event class to store geant4 generated event information @@ -103,54 +106,58 @@ class TRestGeant4Event : public TRestEvent { TH1D* GetZHistogram(Int_t gridElement, std::vector optList); #endif - TVector3 fPrimaryEventOrigin; + TVector3 fPrimaryPosition; + std::vector fPrimaryParticleNames; + std::vector fPrimaryEnergies; + std::vector fPrimaryDirections; - std::vector fPrimaryParticleName; - std::vector fPrimaryEventDirection; - std::vector fPrimaryEventEnergy; + TString fSubEventPrimaryParticleName; + Double_t fSubEventPrimaryEnergy; + TVector3 fSubEventPrimaryPosition; + TVector3 fSubEventPrimaryDirection; - Double_t fTotalDepositedEnergy; - Double_t fSensitiveVolumeEnergy; + Double_t fTotalDepositedEnergy = 0; + Double_t fSensitiveVolumeEnergy = 0; Int_t fNVolumes; std::vector fVolumeStored; std::vector fVolumeStoredNames; std::vector fVolumeDepositedEnergy; - - Int_t fNTracks; - std::vector fTrack; - - Int_t fMaxSubEventID; + std::map>> + fEnergyInVolumePerParticlePerProcess; + std::vector fTracks; public: void SetBoundaries(); void SetBoundaries(Double_t xMin, Double_t xMax, Double_t yMin, Double_t yMax, Double_t zMin, Double_t zMax); - inline TString GetPrimaryEventParticleName(int n) const { - if (fPrimaryParticleName.size() > n) { - return fPrimaryParticleName[n]; - } - return "Not defined"; - } + inline size_t GetNumberOfPrimaries() const { return fPrimaryParticleNames.size(); } + + inline TString GetPrimaryEventParticleName(size_t n = 0) const { return fPrimaryParticleNames[n]; } + inline TVector3 GetPrimaryEventDirection(size_t n = 0) const { return fPrimaryDirections[n]; } + inline TVector3 GetPrimaryEventOrigin() const { return fPrimaryPosition; } + inline Double_t GetPrimaryEventEnergy(size_t n = 0) const { return fPrimaryEnergies[n]; } - TVector3 GetPrimaryEventDirection(Int_t n = 0) const { return fPrimaryEventDirection[n]; } - TVector3 GetPrimaryEventOrigin() const { return fPrimaryEventOrigin; } - Double_t GetPrimaryEventEnergy(Int_t n = 0) const { return fPrimaryEventEnergy[n]; } + inline Bool_t IsSubEvent() const { return fSubEventID > 0; } + inline TString GetSubEventPrimaryEventParticleName() const { return fSubEventPrimaryParticleName; } + inline TVector3 GetSubEventPrimaryEventDirection() const { return fSubEventPrimaryDirection; } + inline TVector3 GetSubEventPrimaryEventOrigin() const { return fSubEventPrimaryPosition; } + inline Double_t GetSubEventPrimaryEventEnergy() const { return fSubEventPrimaryEnergy; } - Int_t GetNumberOfHits(Int_t volID = -1) const; - inline Int_t GetNumberOfTracks() const { return fNTracks; } - inline Int_t GetNumberOfPrimaries() const { return fPrimaryEventDirection.size(); } + size_t GetNumberOfHits(Int_t volID = -1) const; + size_t GetNumberOfPhysicalHits(Int_t volID = -1) const; + + inline const std::vector& GetTracks() const { return fTracks; } + inline size_t GetNumberOfTracks() const { return fTracks.size(); } inline Int_t GetNumberOfActiveVolumes() const { return fNVolumes; } inline Int_t isVolumeStored(int n) const { return fVolumeStored[n]; } - inline const TRestGeant4Track& GetTrack(int n) const { return fTrack[n]; } - inline TRestGeant4Track* GetTrackPointer(int n) { return &fTrack[n]; } - TRestGeant4Track* GetTrackByID(int id); - inline Int_t GetNumberOfSubEventIDTracks() const { return fMaxSubEventID + 1; } + inline const TRestGeant4Track& GetTrack(int n) const { return fTracks[n]; } + inline TRestGeant4Track* GetTrackPointer(int n) { return &fTracks[n]; } + TRestGeant4Track* GetTrackByID(Int_t trackID) const; inline Double_t GetTotalDepositedEnergy() const { return fTotalDepositedEnergy; } - Double_t GetTotalDepositedEnergyFromTracks() const; inline Double_t GetEnergyDepositedInVolume(Int_t volID) const { return fVolumeDepositedEnergy[volID]; } inline Double_t GetSensitiveVolumeEnergy() const { return fSensitiveVolumeEnergy; } TVector3 GetMeanPositionInVolume(Int_t volID) const; @@ -158,16 +165,25 @@ class TRestGeant4Event : public TRestEvent { TVector3 GetLastPositionInVolume(Int_t volID) const; TVector3 GetPositionDeviationInVolume(Int_t volID) const; + std::map>> + GetEnergyInVolumePerParticlePerProcessMap() const; + std::map> GetEnergyInVolumePerProcessMap() const; + std::map> GetEnergyInVolumePerParticleMap() const; + std::map GetEnergyPerProcessMap() const; + std::map GetEnergyPerParticleMap() const; + std::map GetEnergyInVolumeMap() const; + + inline Double_t GetEnergyInVolume(const std::string& volumeName) const { + auto energyMap = GetEnergyInVolumeMap(); + return energyMap[volumeName]; + } + TRestHits GetHits(Int_t volID = -1) const; inline TRestHits GetHitsInVolume(Int_t volID) const { return GetHits(volID); } Int_t GetNumberOfTracksForParticle(const TString& parName) const; Int_t GetEnergyDepositedByParticle(const TString& parName) const; - inline void SetPrimaryEventOrigin(const TVector3& pos) { fPrimaryEventOrigin = pos; } - inline void SetPrimaryEventDirection(const TVector3& dir) { fPrimaryEventDirection.push_back(dir); } - inline void SetPrimaryEventParticleName(const TString& pName) { fPrimaryParticleName.push_back(pName); } - inline void SetPrimaryEventEnergy(Double_t en) { fPrimaryEventEnergy.push_back(en); } inline void ActivateVolumeForStorage(Int_t n) { fVolumeStored[n] = 1; } inline void DisableVolumeForStorage(Int_t n) { fVolumeStored[n] = 0; } @@ -182,11 +198,11 @@ class TRestGeant4Event : public TRestEvent { inline Int_t GetLowestTrackID() const { Int_t lowestID = 0; - if (fNTracks > 0) { + if (GetNumberOfTracks() > 0) { lowestID = GetTrack(0).GetTrackID(); } - for (int i = 0; i < fNTracks; i++) { + for (int i = 0; i < GetNumberOfTracks(); i++) { auto tr = GetTrack(i); if (tr.GetTrackID() < lowestID) lowestID = tr.GetTrackID(); } @@ -194,8 +210,7 @@ class TRestGeant4Event : public TRestEvent { return lowestID; } - void SetTrackSubEventID(Int_t n, Int_t id); - void AddTrack(TRestGeant4Track trk); + const std::set GetUniqueParticles() const; Bool_t ContainsProcessInVolume(Int_t processID, Int_t volumeID = -1) const; inline Bool_t ContainsProcess(Int_t processID) const { return ContainsProcessInVolume(processID, -1); } @@ -226,6 +241,23 @@ class TRestGeant4Event : public TRestEvent { // Destructor virtual ~TRestGeant4Event(); - ClassDefOverride(TRestGeant4Event, 7); // REST event superclass + ClassDefOverride(TRestGeant4Event, 8); // REST event superclass + + // restG4 + public: + TRestGeant4Event(const G4Event*); //! // Implemented in restG4 + bool InsertTrack(const G4Track*); //! + void UpdateTrack(const G4Track*); //! + void InsertStep(const G4Step*); //! + + friend class OutputManager; + + private: + std::map fTrackIDToTrackIndex = {}; + TRestGeant4Hits fInitialStep; //! + + void AddEnergyInVolumeForParticleForProcess(Double_t energy, const std::string& volumeName, + const std::string& particleName, + const std::string& processName); }; #endif diff --git a/inc/TRestGeant4EventViewer.h b/inc/TRestGeant4EventViewer.h index 18f4fcb..e9d40a5 100644 --- a/inc/TRestGeant4EventViewer.h +++ b/inc/TRestGeant4EventViewer.h @@ -22,12 +22,13 @@ class TRestGeant4EventViewer : public TRestEveEventViewer { private: std::vector fHitConnectors; - TRestGeant4Event* fG4Event; + TRestGeant4Event* fG4Event = nullptr; + const TRestGeant4Metadata* fG4Metadata = nullptr; public: void Initialize(); void DeleteCurrentEvent(); - void AddEvent(TRestEvent* ev); + void AddEvent(TRestEvent* event); void NextTrackVertex(Int_t trkID, TVector3 to); void AddTrack(Int_t trkID, Int_t parentID, TVector3 from, TString name); @@ -41,6 +42,6 @@ class TRestGeant4EventViewer : public TRestEveEventViewer { // Destructor ~TRestGeant4EventViewer(); - ClassDef(TRestGeant4EventViewer, 1); // class inherited from TRestEventViewer + ClassDef(TRestGeant4EventViewer, 2); // class inherited from TRestEventViewer }; #endif diff --git a/inc/TRestGeant4GeometryInfo.h b/inc/TRestGeant4GeometryInfo.h index 26673f8..56b617a 100644 --- a/inc/TRestGeant4GeometryInfo.h +++ b/inc/TRestGeant4GeometryInfo.h @@ -20,6 +20,14 @@ class TRestGeant4GeometryInfo { std::map fVolumeNameMap = {}; std::map fVolumeNameReverseMap = {}; + void PopulateFromGeant4World(const G4VPhysicalVolume*); + + inline void InitializeOnDetectorConstruction(const TString& gdmlFilename, + const G4VPhysicalVolume* world) { + PopulateFromGdml(gdmlFilename); + PopulateFromGeant4World(world); + } + public: // Insertion order is important for GDML containers. These containers are filled from GDML only not Geant4 std::vector fGdmlNewPhysicalNames; @@ -44,10 +52,6 @@ class TRestGeant4GeometryInfo { TString GetAlternativeNameFromGeant4PhysicalName(const TString&) const; TString GetGeant4PhysicalNameFromAlternativeName(const TString&) const; - // Int_t GetIDFromVolumeName(const TString&) const; - - void PopulateFromGeant4World(const G4VPhysicalVolume*); - std::vector GetAllLogicalVolumes() const; std::vector GetAllPhysicalVolumes() const; std::vector GetAllAlternativePhysicalVolumes() const; @@ -77,6 +81,10 @@ class TRestGeant4GeometryInfo { return {}; } + inline TVector3 GetPosition(const TString& volume) const { + return fPhysicalToPositionInWorldMap.at(volume); + } + inline bool IsAssembly() const { return fIsAssembly; } void InsertVolumeName(Int_t id, const TString& volumeName); @@ -85,6 +93,8 @@ class TRestGeant4GeometryInfo { Int_t GetIDFromVolume(const TString& volumeName) const; void Print() const; + + friend class DetectorConstruction; }; #endif // REST_TRESTGEANT4GEOMETRYINFO_H diff --git a/inc/TRestGeant4Hits.h b/inc/TRestGeant4Hits.h index da6b351..b65f206 100644 --- a/inc/TRestGeant4Hits.h +++ b/inc/TRestGeant4Hits.h @@ -18,39 +18,45 @@ #ifndef RestCore_TRestGeant4Hits #define RestCore_TRestGeant4Hits -#include -#include #include #include -#include "TObject.h" +#include "TRestGeant4Metadata.h" + +class G4Step; +class TRestGeant4Track; +class TRestGeant4Event; class TRestGeant4Hits : public TRestHits { protected: - TArrayI fVolumeID; - TArrayI fProcessID; // [fNHits] - TArrayF fKineticEnergy; // [fNHits] + std::vector fProcessID = {}; + std::vector fVolumeID = {}; + std::vector fKineticEnergy = {}; + std::vector fMomentumDirection = {}; + + TRestGeant4Track* fTrack = nullptr; //! + TRestGeant4Event* fEvent = nullptr; //! public: - TArrayF fMomentumDirectionX; - TArrayF fMomentumDirectionY; - TArrayF fMomentumDirectionZ; + TRestGeant4Metadata* GetGeant4Metadata() const; + + inline const TRestGeant4Track* GetTrack() const { return fTrack; } + inline void SetTrack(TRestGeant4Track* track) { fTrack = track; } - TVector3 GetMomentumDirection(int n) const { - return {fMomentumDirectionX[n], fMomentumDirectionY[n], fMomentumDirectionZ[n]}; - } + inline const TRestGeant4Event* GetEvent() const { return fEvent; } + inline void SetEvent(TRestGeant4Event* event) { fEvent = event; } - Int_t GetProcess(int n) const { return fProcessID[n]; } + inline TVector3 GetMomentumDirection(size_t n) const { return fMomentumDirection[n]; } + + inline Int_t GetProcess(size_t n) const { return fProcessID[n]; } - void AddG4Hit(TVector3 pos, Double_t en, Double_t hit_global_time, Int_t process, Int_t volume, - Double_t eKin, TVector3 momentumDirection); void RemoveG4Hits(); - Int_t GetHitProcess(int n) const { return fProcessID[n]; } - Int_t GetHitVolume(int n) const { return fVolumeID[n]; } - Int_t GetVolumeId(int n) const { return fVolumeID[n]; } - Double_t GetKineticEnergy(int n) const { return fKineticEnergy[n]; } + inline Int_t GetHitProcess(size_t n) const { return fProcessID[n]; } + inline Int_t GetHitVolume(size_t n) const { return fVolumeID[n]; } + inline Int_t GetVolumeId(size_t n) const { return fVolumeID[n]; } + inline Double_t GetKineticEnergy(size_t n) const { return fKineticEnergy[n]; } Double_t GetEnergyInVolume(Int_t volumeID) const; @@ -65,6 +71,10 @@ class TRestGeant4Hits : public TRestHits { // Destructor virtual ~TRestGeant4Hits(); - ClassDef(TRestGeant4Hits, 6); // REST event superclass + ClassDef(TRestGeant4Hits, 7); // REST event superclass + + // restG4 + public: + void InsertStep(const G4Step*); }; #endif diff --git a/inc/TRestGeant4Metadata.h b/inc/TRestGeant4Metadata.h index fdc09ea..0d1037d 100644 --- a/inc/TRestGeant4Metadata.h +++ b/inc/TRestGeant4Metadata.h @@ -41,52 +41,7 @@ #include "TRestGeant4GeometryInfo.h" #include "TRestGeant4ParticleSource.h" #include "TRestGeant4PhysicsInfo.h" -//------------------------------------------------------------------------------------------------------------------------ -// -// * This section was added by Luis A. Obis (lobis@unizar.es) on 17/06/2019 -// -// Here we add all the possible options for different configurations such as all the types of generators, etc. -// We use a structure called 'enum' and a function to clean the strings so that we can easily implement case -// insensitivity or more options such as ignoring underscores. -// - -namespace g4_metadata_parameters { -std::string CleanString(std::string); - -enum class generator_types { - CUSTOM, - VOLUME, - SURFACE, - POINT, -}; -extern std::map generator_types_map; - -enum class generator_shapes { - GDML, - WALL, - CIRCLE, - BOX, - SPHERE, - CYLINDER, -}; -extern std::map generator_shapes_map; - -enum class energy_dist_types { - TH1D, - MONO, - FLAT, - LOG, -}; -extern std::map energy_dist_types_map; - -enum class angular_dist_types { - TH1D, - ISOTROPIC, - FLUX, - BACK_TO_BACK, -}; -extern std::map angular_dist_types_map; -} // namespace g4_metadata_parameters +#include "TRestGeant4PrimaryGeneratorInfo.h" /// The main class to store the *Geant4* simulation conditions that will be used by *restG4*. class TRestGeant4Metadata : public TRestMetadata { @@ -98,15 +53,20 @@ class TRestGeant4Metadata : public TRestMetadata { void ReadGenerator(); void ReadParticleSource(TRestGeant4ParticleSource* src, TiXmlElement* sourceDefinition); - void ReadStorage(); + void ReadDetector(); void ReadBiasing(); + bool fDetectorSectionInitialized = false; //! + /// Class used to store and retrieve geometry info TRestGeant4GeometryInfo fGeant4GeometryInfo; /// Class used to store and retrieve physics info such as process names or particle names TRestGeant4PhysicsInfo fGeant4PhysicsInfo; + /// Class used to store and retrieve Geant4 primary generator info + TRestGeant4PrimaryGeneratorInfo fGeant4PrimaryGeneratorInfo; + /// The version of Geant4 used to generate the data TString fGeant4Version; @@ -122,41 +82,9 @@ class TRestGeant4Metadata : public TRestMetadata { /// A GDML materials reference introduced in the header of the GDML of materials definition TString fMaterialsReference; - /// Type of spatial generator (surface, volume, custom) - TString fGenType; - - /// Shape of spatial generator (point, wall, gdml, sphere, etc) - TString fGenShape; - - /// The volume name where the events are generated, in case of volume or - /// surface generator types. - TString fGenFrom; - - /// The position of the generator for virtual, and point, generator types. - TVector3 fGenPosition; - - /// \brief A 3d-std::vector with the angles, measured in degrees, of a XYZ rotation - /// applied to the virtual generator. This rotation is used by virtualWall, - /// virtualCircleWall and virtualCylinder generators. - TVector3 fGenRotationAxis; - - /// \brief degrees of rotation for generator virtual shape along the axis - Double_t fGenRotationDegree; - - /// \brief The size of the generator. Stores up to three deminsions according to the shape - /// box: length, width, height - /// sphere: radius - /// wall: length, width - /// circle: radius - /// cylinder: radius, length - TVector3 fGenSize; - - /// \brief Defines density distribution of the generator shape. rho=F(x,y,z), in range 0~1 - TString fGenDensityFunction; - /// \brief A 2d-std::vector storing the energy range, in keV, to decide if a particular /// event should be written to disk or not. - TVector2 fEnergyRangeStored; + TVector2 fEnergyRangeStored = {0, 1E20}; /// \brief A std::vector to store the names of the active volumes. std::vector fActiveVolumes; @@ -174,7 +102,7 @@ class TRestGeant4Metadata : public TRestMetadata { /// \brief The number of biasing volumes used in the simulation. If zero, no biasing /// technique is used. - Int_t fNBiasingVolumes; + Int_t fNBiasingVolumes = 0; /// A std::vector containing the biasing volume properties. std::vector fBiasingVolumes; @@ -186,27 +114,46 @@ class TRestGeant4Metadata : public TRestMetadata { /// \brief A time gap, in us, determining if an energy hit should be considered (and /// stored) as an independent event. - Double_t fSubEventTimeDelay; + Double_t fSubEventTimeDelay = 100; /// \brief Defines if a radioactive isotope decay is simulated in full chain /// (fFullChain=true), or just a single decay (fFullChain=false). - Bool_t fFullChain; + Bool_t fFullChain = true; /// \brief The volume that serves as trigger for data storage. Only events that /// deposit a non-zero energy on this volume will be registered. - TString fSensitiveVolume; + std::vector fSensitiveVolumes; /// The number of events simulated, or to be simulated. - Int_t fNEvents; + Int_t fNEvents = 0; + + /// The number of events the user requested to be on the file + Int_t fNRequestedEntries = 0; + + /// Time before simulation is ended and saved + Int_t fSimulationMaxTimeSeconds = 0; - /// \brief The seed value used for Geant4 random event generator. If it is zero - /// its value will be assigned using the system timestamp. + /// \brief The seed value used for Geant4 random event generator. + /// If it is zero its value will be assigned using the system timestamp. Long_t fSeed = 0; /// \brief If this parameter is set to 'true' it will save all events even if they leave no energy in the - /// sensitive volume (used for debugging pourposes). It is set to 'false' by default. + /// sensitive volume (used for debugging purposes). It is set to 'false' by default. Bool_t fSaveAllEvents = false; + /// \brief Sets all volume as active without having to explicitly list them. + Bool_t fActivateAllVolumes = false; //! + + /// \brief If activated will remove tracks not present in volumes marked as "keep" or "sensitive". + Bool_t fRemoveUnwantedTracks = false; + + /// \brief Option for 'removeUnwantedTracks', if enabled tracks with hits in volumes will be kept event if + /// they deposit zero energy (such as neutron captures) + Bool_t fRemoveUnwantedTracksKeepZeroEnergyTracks = false; + + /// \brief A container related to fRemoveUnwantedTracks. + std::set fRemoveUnwantedTracksVolumesToKeep; + /// If this parameter is set to 'true' it will print out on screen every time 10k events are reached. Bool_t fPrintProgress = false; //! @@ -219,6 +166,8 @@ class TRestGeant4Metadata : public TRestMetadata { TVector3 fMagneticField = TVector3(0, 0, 0); public: + std::set fActiveVolumesSet = {}; //! // Used for faster lookup + /// \brief Returns the random seed that was used to generate the corresponding /// geant4 dataset. inline Long_t GetSeed() const { return fSeed; } @@ -229,6 +178,11 @@ class TRestGeant4Metadata : public TRestMetadata { /// \brief Returns an immutable reference to the physics info inline const TRestGeant4PhysicsInfo& GetGeant4PhysicsInfo() const { return fGeant4PhysicsInfo; } + /// \brief Returns an immutable reference to the primary generator info + inline const TRestGeant4PrimaryGeneratorInfo& GetGeant4PrimaryGeneratorInfo() const { + return fGeant4PrimaryGeneratorInfo; + } + /// \brief Returns a std::string with the version of Geant4 used on the event data /// simulation inline TString GetGeant4Version() const { return fGeant4Version; } @@ -247,43 +201,6 @@ class TRestGeant4Metadata : public TRestMetadata { /// Returns the reference provided at the materials file header inline TString GetMaterialsReference() const { return fMaterialsReference; } - /// \brief Returns a std::string specifying the generator type (volume, surface, point, - /// virtualWall, etc ) - inline TString GetGeneratorType() const { return fGenType; } - - /// \brief Returns a std::string specifying the generator shape (point, wall, box, etc ) - inline TString GetGeneratorShape() const { return fGenShape; } - - /// \brief Returns the name of the GDML volume where primary events are - /// produced. This value has meaning only when using volume or surface - /// generator types. - inline TString GetGeneratedFrom() const { return fGenFrom; } - - /// \brief Returns the name of the GDML volume where primary events are - /// produced. This value has meaning only when using volume or surface - /// generator types. - inline TString GetGDMLGeneratorVolume() const { return fGenFrom; } - - /// \brief Returns a 3d-std::vector with the position of the primary event - /// generator. This value has meaning only when using point and virtual - /// generator types. - inline TVector3 GetGeneratorPosition() const { return fGenPosition; } - - /// \brief Returns a 3d-std::vector, fGenRotation, with the XYZ rotation angle - /// values in degrees. This value is used by virtualWall, virtualCircleWall - /// and virtualCylinder generator types. - inline TVector3 GetGeneratorRotationAxis() const { return fGenRotationAxis; } - - /// \brief Returns the degree of rotation - inline Double_t GetGeneratorRotationDegree() const { return fGenRotationDegree; } - - /// \brief Returns the main spatial dimension of virtual generator. - /// It is the size of a virtualBox. - inline TVector3 GetGeneratorSize() const { return fGenSize; } - - /// \brief Returns the density function of the generator - inline TString GetGeneratorSpatialDensity() const { return fGenDensityFunction; } - /// \brief Returns true in case full decay chain simulation is enabled. inline Bool_t isFullChainActivated() const { return fFullChain; } @@ -316,26 +233,9 @@ class TRestGeant4Metadata : public TRestMetadata { /// Sets the value of the Geant4 version inline void SetGeant4Version(const TString& versionString) { fGeant4Version = versionString; } - /// \brief Sets the generator type. I.e. volume, surface, point, virtualWall, - /// virtualCylinder, etc. - inline void SetGeneratorType(TString type) { fGenType = type; } - - /// \brief Sets the generator main spatial dimension. In a virtual generator is the - /// radius of cylinder, size of wall, etc. - inline void SetGeneratorSize(const TVector3& size) { fGenSize = size; } - /// Enables/disables the full chain decay generation. inline void SetFullChain(Bool_t fullChain) { fFullChain = fullChain; } - /// Sets the position of the virtual generator using a TVector3. - inline void SetGeneratorPosition(const TVector3& pos) { fGenPosition = pos; } - - /// Sets the position of the virtual generator using x,y,z coordinates. - inline void SetGeneratorPosition(double x, double y, double z) { fGenPosition = TVector3(x, y, z); } - - /// Sets the number of events to be simulated. - inline void SetNEvents(Int_t n) { fNEvents = n; } - /// It sets the location of the geometry files inline void SetGeometryPath(std::string path) { fGeometryPath = path; } @@ -350,15 +250,19 @@ class TRestGeant4Metadata : public TRestMetadata { /// Returns the number of events to be simulated. inline Int_t GetNumberOfEvents() const { return fNEvents; } + + inline Int_t GetNumberOfRequestedEntries() const { return fNRequestedEntries; } + + inline Int_t GetSimulationMaxTimeSeconds() const { return fSimulationMaxTimeSeconds; } + /////////////////////////////////////////////////////////// // Direct access to sources definition in primary generator - /////////////////////////////////////////////////////////// /// Returns the number of primary event sources defined in the generator. inline Int_t GetNumberOfSources() const { return fParticleSource.size(); } /// Returns the name of the particle source with index n (Geant4 based names). - inline TRestGeant4ParticleSource* GetParticleSource(int n) { return fParticleSource[n]; } + inline TRestGeant4ParticleSource* GetParticleSource(size_t n = 0) const { return fParticleSource[n]; } /// Remove all the particle sources. void RemoveParticleSources(); @@ -381,20 +285,39 @@ class TRestGeant4Metadata : public TRestMetadata { inline Int_t isBiasingActive() const { return fBiasingVolumes.size(); } /// Returns a std::string with the name of the sensitive volume. - inline TString GetSensitiveVolume() const { return fSensitiveVolume; } + inline TString GetSensitiveVolume(int n = 0) const { return fSensitiveVolumes[n]; } + + inline size_t GetNumberOfSensitiveVolumes() const { return fSensitiveVolumes.size(); } + + inline const std::vector& GetSensitiveVolumes() const { return fSensitiveVolumes; } /// Sets the name of the sensitive volume - inline void SetSensitiveVolume(const TString& sensitiveVolume) { fSensitiveVolume = sensitiveVolume; } + inline void SetNumberOfEvents(Int_t n) { fNEvents = n; } + + inline void SetNumberOfRequestedEntries(Int_t n) { fNRequestedEntries = n; } + + inline void SetSimulationMaxTimeSeconds(Int_t seconds) { fSimulationMaxTimeSeconds = seconds; } + + /// Sets the name of the sensitive volume + inline void InsertSensitiveVolume(const TString& volume) { + for (const auto& sensitiveVolume : fSensitiveVolumes) { + // Do not add duplicate volumes + if (volume == sensitiveVolume) { + return; + } + } + fSensitiveVolumes.push_back(volume); + } /// \brief Returns the probability per event to register (write to disk) hits in the /// storage volume with index n. - inline Double_t GetStorageChance(Int_t n) { return fChance[n]; } + inline Double_t GetStorageChance(Int_t n) const { return fChance[n]; } /// Returns the probability per event to register (write to disk) hits in a /// GDML volume given its geometry name. - Double_t GetStorageChance(TString vol); + Double_t GetStorageChance(TString volume); - Double_t GetMaxStepSize(TString vol); + Double_t GetMaxStepSize(const TString& volume); /// Returns the minimum event energy required for an event to be stored. inline Double_t GetMinimumEnergyStored() const { return fEnergyRangeStored.X(); } @@ -406,15 +329,31 @@ class TRestGeant4Metadata : public TRestMetadata { /// selected for data storage. inline Int_t GetNumberOfActiveVolumes() const { return fActiveVolumes.size(); } + inline bool IsActiveVolume(const char* volumeName) const { + return fActiveVolumesSet.count(volumeName) > 0; + } //! + + inline bool IsKeepTracksVolume(const char* volumeName) const { + return fRemoveUnwantedTracksVolumesToKeep.count(volumeName) > 0; + } + /// Returns a std::string with the name of the active volume with index n - inline TString GetActiveVolumeName(Int_t n) { return fActiveVolumes[n]; } + inline TString GetActiveVolumeName(Int_t n) const { return fActiveVolumes[n]; } + + inline std::vector GetActiveVolumes() const { return fActiveVolumes; } + + inline bool GetRemoveUnwantedTracks() const { return fRemoveUnwantedTracks; } + + inline bool GetRemoveUnwantedTracksKeepZeroEnergyTracks() const { + return fRemoveUnwantedTracksKeepZeroEnergyTracks; + } /// Returns the world magnetic field in Tesla inline TVector3 GetMagneticField() const { return fMagneticField; } Int_t GetActiveVolumeID(TString name); - Bool_t isVolumeStored(TString volName); + Bool_t isVolumeStored(const TString& volume) const; void SetActiveVolume(const TString& name, Double_t chance, Double_t maxStep = 0); @@ -425,10 +364,11 @@ class TRestGeant4Metadata : public TRestMetadata { ~TRestGeant4Metadata(); - ClassDefOverride(TRestGeant4Metadata, 10); + ClassDefOverride(TRestGeant4Metadata, 11); // Allow modification of otherwise inaccessible / immutable members that shouldn't be modified by the user friend class SteppingAction; friend class DetectorConstruction; + friend class TRestGeant4Hits; }; #endif // RestCore_TRestGeant4Metadata diff --git a/inc/TRestGeant4Particle.h b/inc/TRestGeant4Particle.h index 2734d3d..3722d57 100644 --- a/inc/TRestGeant4Particle.h +++ b/inc/TRestGeant4Particle.h @@ -21,14 +21,12 @@ #include -#include "TObject.h" - -class TRestGeant4Particle : public TObject { +class TRestGeant4Particle { protected: TString fParticleName; Double_t fExcitationLevel = 0; - TVector3 fDirection; - Double_t fEnergy; + TVector3 fDirection = {1, 0, 0}; + Double_t fEnergy = 0; Int_t fCharge = 0; TVector3 fOrigin; @@ -40,17 +38,17 @@ class TRestGeant4Particle : public TObject { inline Int_t GetParticleCharge() const { return fCharge; } inline TVector3 GetOrigin() const { return fOrigin; } - void SetParticle(TRestGeant4Particle ptcle) { - fExcitationLevel = ptcle.GetExcitationLevel(); - fParticleName = ptcle.GetParticleName(); - fEnergy = ptcle.GetEnergy(); - fDirection = ptcle.GetMomentumDirection(); - fOrigin = ptcle.fOrigin; + void SetParticle(TRestGeant4Particle particle) { + fExcitationLevel = particle.GetExcitationLevel(); + fParticleName = particle.GetParticleName(); + fEnergy = particle.GetEnergy(); + fDirection = particle.GetMomentumDirection(); + fOrigin = particle.fOrigin; } - void SetParticleName(TString ptcle) { fParticleName = ptcle; } - void SetExcitationLevel(Double_t eenergy) { - fExcitationLevel = eenergy; + void SetParticleName(TString particle) { fParticleName = particle; } + void SetExcitationLevel(Double_t excitationEnergy) { + fExcitationLevel = excitationEnergy; if (fExcitationLevel < 0) fExcitationLevel = 0; } @@ -65,6 +63,6 @@ class TRestGeant4Particle : public TObject { // Destructor virtual ~TRestGeant4Particle(); - ClassDef(TRestGeant4Particle, 3); + ClassDef(TRestGeant4Particle, 4); }; #endif diff --git a/inc/TRestGeant4ParticleSource.h b/inc/TRestGeant4ParticleSource.h index 57dc2c5..ca18fe9 100644 --- a/inc/TRestGeant4ParticleSource.h +++ b/inc/TRestGeant4ParticleSource.h @@ -17,16 +17,16 @@ #ifndef RestCore_TRestGeant4ParticleSource #define RestCore_TRestGeant4ParticleSource +#include +#include +#include #include #include #include #include -#include "TObject.h" -// #include "TRestGeant4Particle.h" -#include "TRestMetadata.h" class TRestGeant4ParticleSource : public TRestGeant4Particle, public TRestMetadata { private: @@ -35,61 +35,112 @@ class TRestGeant4ParticleSource : public TRestGeant4Particle, public TRestMetada bool ReadOldDecay0File(TString fileName); protected: - TString fAngularDistType; - TString fEnergyDistType; - TVector2 fEnergyRange; + TString fAngularDistributionType = "Flux"; + TString fAngularDistributionFilename; + TString fAngularDistributionNameInFile; + TF1* fAngularDistributionFunction = nullptr; - TString fSpectrumFilename; - TString fSpectrumName; - - TString fAngularFilename; - TString fAngularName; + TString fEnergyDistributionType = "Mono"; + TString fEnergyDistributionFilename; + TString fEnergyDistributionNameInFile; + TF1* fEnergyDistributionFunction = nullptr; + TVector2 fEnergyDistributionRange; TString fGenFilename; // store a set of generated particles - std::vector fParticles; // + std::vector fParticles; // store a list of particle set templates that could be used to generate fParticles std::vector> fParticlesTemplate; //! // the random method to generate 0~1 equally distributed numbers - double (*fRndmMethod)(); //! + double (*fRandomMethod)(); //! public: virtual void Update(); virtual void InitFromConfigFile(); static TRestGeant4ParticleSource* instantiate(std::string model = ""); - inline TString GetParticleName() const { return fParticleName; } - inline TString GetAngularDistType() const { return fAngularDistType; } - inline TVector3 GetDirection() const { return fDirection; } - inline TString GetEnergyDistType() const { return fEnergyDistType; } - inline TVector2 GetEnergyRange() const { return fEnergyRange; } - inline Double_t GetMinEnergy() const { return fEnergyRange.X(); } - inline Double_t GetMaxEnergy() const { return fEnergyRange.Y(); } - inline TString GetSpectrumFilename() const { return fSpectrumFilename; } - inline TString GetSpectrumName() const { return fSpectrumName; } - inline TString GetAngularFilename() const { return fAngularFilename; } - inline TString GetAngularName() const { return fAngularName; } + inline TVector3 GetDirection() const { + if (fDirection.Mag() <= 0) { + std::cout << "TRestGeant4ParticleSource::GetDirection: " + << "Direction cannot be the zero vector" << std::endl; + exit(1); + } + return fDirection; + } + + inline TString GetEnergyDistributionType() const { return fEnergyDistributionType; } + inline TVector2 GetEnergyDistributionRange() const { return fEnergyDistributionRange; } + inline Double_t GetEnergyDistributionRangeMin() const { return fEnergyDistributionRange.X(); } + inline Double_t GetEnergyDistributionRangeMax() const { return fEnergyDistributionRange.Y(); } + inline TString GetEnergyDistributionFilename() const { return fEnergyDistributionFilename; } + inline TString GetEnergyDistributionNameInFile() const { return fEnergyDistributionNameInFile; } + inline const TF1* GetEnergyDistributionFunction() const { return fEnergyDistributionFunction; } + + inline TString GetAngularDistributionType() const { return fAngularDistributionType; } + inline TString GetAngularDistributionFilename() const { return fAngularDistributionFilename; } + inline TString GetAngularDistributionNameInFile() const { return fAngularDistributionNameInFile; } + inline const TF1* GetAngularDistributionFunction() const { return fAngularDistributionFunction; } + inline TString GetGenFilename() const { return fGenFilename; } inline std::vector GetParticles() const { return fParticles; } - inline void SetAngularDistType(const TString& type) { fAngularDistType = type; } - inline void SetEnergyDistType(const TString& type) { fEnergyDistType = type; } - inline void SetEnergyRange(const TVector2& range) { fEnergyRange = range; } - inline void SetSpectrumFilename(const TString& spectrumFilename) { fSpectrumFilename = spectrumFilename; } - inline void SetSpectrumName(const TString& spectrumName) { fSpectrumName = spectrumName; } - inline void SetAngularFilename(const TString& angFilename) { fAngularFilename = angFilename; } - inline void SetAngularName(const TString& angularName) { fAngularName = angularName; } + inline void SetAngularDistributionType(const TString& type) { + fAngularDistributionType = TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypesToString( + TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes(type.Data())); + } + inline void SetAngularDistributionFilename(const TString& filename) { + fAngularDistributionFilename = filename; + } + inline void SetAngularDistributionNameInFile(const TString& name) { + fAngularDistributionNameInFile = name; + } + inline void SetAngularDistributionFormula(const TString& formula) { + fAngularDistributionFunction = + (TF1*)TRestGeant4PrimaryGeneratorTypes::AngularDistributionFormulasToRootFormula( + TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionFormulas(formula.Data())) + .Clone(); + } + + inline void SetEnergyDistributionType(const TString& type) { + fEnergyDistributionType = TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypesToString( + TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionTypes(type.Data())); + } + inline void SetEnergyDistributionRange(const TVector2& range) { + fEnergyDistributionRange = range; + if (fEnergyDistributionRange.X() < 0) { + fEnergyDistributionRange.Set(0, fEnergyDistributionRange.Y()); + } + if (fEnergyDistributionRange.Y() < 0) { + fEnergyDistributionRange.Set(fEnergyDistributionRange.X(), 0); + } + if (fEnergyDistributionRange.X() > fEnergyDistributionRange.Y()) { + fEnergyDistributionRange.Set(fEnergyDistributionRange.Y(), fEnergyDistributionRange.X()); + } + } + inline void SetEnergyDistributionFilename(const TString& filename) { + fEnergyDistributionFilename = filename; + } + inline void SetEnergyDistributionNameInFile(const TString& name) { fEnergyDistributionNameInFile = name; } + inline void SetEnergyDistributionFormula(const TString& formula) { + fEnergyDistributionFunction = + (TF1*)TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToRootFormula( + TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionFormulas(formula.Data())) + .Clone(); + } + inline void SetGenFilename(const TString& name) { fGenFilename = name; } - inline void SetRndmMethod(double (*method)()) { fRndmMethod = method; } + + inline void SetRandomMethod(double (*method)()) { fRandomMethod = method; } inline void AddParticle(const TRestGeant4Particle& particle) { fParticles.push_back(particle); } + inline void RemoveParticles() { fParticles.clear(); } + inline void RemoveTemplates() { fParticlesTemplate.clear(); } inline void FlushParticlesTemplate() { fParticlesTemplate.push_back(fParticles); fParticles.clear(); } - inline void RemoveTemplates() { fParticlesTemplate.clear(); } virtual void PrintParticleSource(); @@ -98,6 +149,6 @@ class TRestGeant4ParticleSource : public TRestGeant4Particle, public TRestMetada // Destructor virtual ~TRestGeant4ParticleSource(); - ClassDef(TRestGeant4ParticleSource, 3); + ClassDef(TRestGeant4ParticleSource, 4); }; #endif diff --git a/inc/TRestGeant4PhysicsInfo.h b/inc/TRestGeant4PhysicsInfo.h index 50edc7d..d190497 100644 --- a/inc/TRestGeant4PhysicsInfo.h +++ b/inc/TRestGeant4PhysicsInfo.h @@ -4,12 +4,14 @@ #include #include +#include +#include #include class G4VProcess; class TRestGeant4PhysicsInfo { - ClassDef(TRestGeant4PhysicsInfo, 1); + ClassDef(TRestGeant4PhysicsInfo, 2); private: std::map fProcessNamesMap = {}; @@ -18,14 +20,22 @@ class TRestGeant4PhysicsInfo { std::map fParticleNamesMap = {}; std::map fParticleNamesReverseMap = {}; + std::map fProcessTypesMap = {}; // process name -> process type + + std::mutex fMutex; //! public: TString GetProcessName(Int_t id) const; Int_t GetProcessID(const TString& processName) const; - void InsertProcessName(Int_t id, const TString& processName); + void InsertProcessName(Int_t id, const TString& processName, const TString& processType); + std::set GetAllParticles() const; TString GetParticleName(Int_t id) const; Int_t GetParticleID(const TString& processName) const; void InsertParticleName(Int_t id, const TString& particleName); + std::set GetAllProcesses() const; + + TString GetProcessType(const TString& processName) const; + std::set GetAllProcessTypes() const; public: inline TRestGeant4PhysicsInfo() = default; diff --git a/inc/TRestGeant4PrimaryGeneratorInfo.h b/inc/TRestGeant4PrimaryGeneratorInfo.h new file mode 100644 index 0000000..06e24ff --- /dev/null +++ b/inc/TRestGeant4PrimaryGeneratorInfo.h @@ -0,0 +1,158 @@ + +#ifndef REST_TRESTGEANT4PRIMARYGENERATORINFO_H +#define REST_TRESTGEANT4PRIMARYGENERATORINFO_H + +#include +#include +#include + +#include + +namespace TRestGeant4PrimaryGeneratorTypes { + +enum class SpatialGeneratorTypes { + CUSTOM, + VOLUME, + SURFACE, + POINT, +}; + +std::string SpatialGeneratorTypesToString(const SpatialGeneratorTypes&); +SpatialGeneratorTypes StringToSpatialGeneratorTypes(const std::string&); + +enum class SpatialGeneratorShapes { + GDML, + WALL, + CIRCLE, + BOX, + SPHERE, + CYLINDER, +}; + +std::string SpatialGeneratorShapesToString(const SpatialGeneratorShapes&); +SpatialGeneratorShapes StringToSpatialGeneratorShapes(const std::string&); + +enum class EnergyDistributionTypes { + TH1D, + FORMULA, + MONO, + FLAT, + LOG, +}; + +std::string EnergyDistributionTypesToString(const EnergyDistributionTypes&); +EnergyDistributionTypes StringToEnergyDistributionTypes(const std::string&); + +enum class EnergyDistributionFormulas { + COSMIC_MUONS, + COSMIC_NEUTRONS, + COSMIC_GAMMAS, +}; + +std::string EnergyDistributionFormulasToString(const EnergyDistributionFormulas&); +EnergyDistributionFormulas StringToEnergyDistributionFormulas(const std::string&); +TF1 EnergyDistributionFormulasToRootFormula(const EnergyDistributionFormulas&); + +enum class AngularDistributionTypes { + TH1D, + FORMULA, + ISOTROPIC, + FLUX, + BACK_TO_BACK, +}; + +std::string AngularDistributionTypesToString(const AngularDistributionTypes&); +AngularDistributionTypes StringToAngularDistributionTypes(const std::string&); + +enum class AngularDistributionFormulas { + COS2, + COS3, +}; + +std::string AngularDistributionFormulasToString(const AngularDistributionFormulas&); +AngularDistributionFormulas StringToAngularDistributionFormulas(const std::string&); +TF1 AngularDistributionFormulasToRootFormula(const AngularDistributionFormulas&); + +} // namespace TRestGeant4PrimaryGeneratorTypes + +class TiXmlElement; + +class TRestGeant4PrimaryGeneratorInfo { + ClassDef(TRestGeant4PrimaryGeneratorInfo, 1); + + public: + void Print() const; + + private: + /* Members related to spatial generator */ + /// Type of spatial generator (point, surface, volume, custom) + TString fSpatialGeneratorType = "point"; + + /// Shape of spatial generator (wall, GDML, sphere, etc) + TString fSpatialGeneratorShape; + + /// The volume name where the events are generated, in case of volume or surface generator types + TString fSpatialGeneratorFrom; + + /// The position of the generator for virtual, and point, generator types + TVector3 fSpatialGeneratorPosition; + + /// \brief A 3d-std::vector with the angles, measured in degrees, of a XYZ rotation + /// applied to the virtual generator. This rotation is used by virtualWall, + /// virtualCircleWall and virtualCylinder generators. + TVector3 fSpatialGeneratorRotationAxis; + + /// \brief degrees of rotation for generator virtual shape along the axis + Double_t fSpatialGeneratorRotationValue = 0; + + /// \brief The size of the generator. Stores up to three dimensions according to the shape + /// box: length, width, height + /// sphere: radius + /// wall: length, width + /// circle: radius + /// cylinder: radius, length + TVector3 fSpatialGeneratorSize; + + /// \brief Defines density distribution of the generator shape. rho=F(x,y,z), in range 0~1 + TString fSpatialGeneratorSpatialDensityFunction; + + public: + /// \brief Returns a std::string specifying the generator type (volume, surface, point, + /// virtualWall, etc ) + inline TString GetSpatialGeneratorType() const { return fSpatialGeneratorType; } + + /// \brief Returns a std::string specifying the generator shape (point, wall, box, etc ) + inline TString GetSpatialGeneratorShape() const { return fSpatialGeneratorShape; } + + /// \brief Returns the name of the GDML volume where primary events are + /// produced. This value has meaning only when using volume or surface + /// generator types. + inline TString GetSpatialGeneratorFrom() const { return fSpatialGeneratorFrom; } + + /// \brief Returns a 3d-std::vector with the position of the primary event + /// generator. This value has meaning only when using point and virtual + /// generator types. + inline TVector3 GetSpatialGeneratorPosition() const { return fSpatialGeneratorPosition; } + + /// \brief Returns a 3d-std::vector, fGenRotation, with the XYZ rotation angle + /// values in degrees. This value is used by virtualWall, virtualCircleWall + /// and virtualCylinder generator types. + inline TVector3 GetSpatialGeneratorRotationAxis() const { return fSpatialGeneratorRotationAxis; } + + /// \brief Returns the degree of rotation + inline Double_t GetSpatialGeneratorRotationValue() const { return fSpatialGeneratorRotationValue; } + + /// \brief Returns the main spatial dimension of virtual generator. + /// It is the size of a virtualBox. + inline TVector3 GetSpatialGeneratorSize() const { return fSpatialGeneratorSize; } + + /// \brief Returns the density function of the generator + inline TString GetSpatialGeneratorSpatialDensityFunction() const { + return fSpatialGeneratorSpatialDensityFunction; + } + + friend class TRestGeant4Metadata; + friend class DetectorConstruction; +}; + +#endif // REST_TRESTGEANT4PRIMARYGENERATORINFO_H diff --git a/inc/TRestGeant4Track.h b/inc/TRestGeant4Track.h index a863a3f..0024bd5 100644 --- a/inc/TRestGeant4Track.h +++ b/inc/TRestGeant4Track.h @@ -16,71 +16,85 @@ #ifndef RestCore_TRestGeant4Track #define RestCore_TRestGeant4Track -#include #include -#include #include #include -#include -#include - -#include "TObject.h" +#include "TRestGeant4Hits.h" class TRestGeant4Event; class TRestGeant4Metadata; +class G4Track; +class G4Step; // Perhaps there might be need for a mother class TRestTrack (if there is future need) -class TRestGeant4Track : public TObject { +class TRestGeant4Track { protected: - Int_t fTrack_ID; - Int_t fParent_ID; - - Int_t fSubEventId; + Int_t fTrackID; + Int_t fParentID; - // We must change this to save space! (Might be not needed after all) TString fParticleName; + + TRestGeant4Hits fHits; + TString fCreatorProcess; - Double_t fGlobalTimestamp; // in seconds precision - Double_t fTrackTimestamp; // in ns precision (seconds have been removed) - Double_t fKineticEnergy; + std::vector fSecondaryTrackIDs; - TRestGeant4Hits fHits; + Double_t fGlobalTimestamp; + Double_t fTimeLength; + + Double_t fInitialKineticEnergy; + Double_t fLength; + + TVector3 fInitialPosition; - TVector3 fTrackOrigin; + Double_t fWeight = 1; // Used for biasing TRestGeant4Event* fEvent = nullptr; //! public: - inline void Initialize() { - RemoveHits(); - fSubEventId = 0.; - } - inline const TRestGeant4Hits& GetHits() const { return fHits; } - inline const TRestGeant4Event* GetEvent() const { return fEvent; } + const TRestGeant4Metadata* GetGeant4Metadata() const; + inline void SetEvent(TRestGeant4Event* event) { fEvent = event; } + inline void SetHits(const TRestGeant4Hits& hits) { + fHits = hits; + fHits.SetTrack(this); + } inline TString GetCreatorProcess() const { return fCreatorProcess; } - inline void SetCreatorProcess(const TString& processName) { fCreatorProcess = processName; } - inline Double_t GetEnergy() const { return fHits.GetEnergy(); } + inline void AddSecondaryTrackID(Int_t trackID) { fSecondaryTrackIDs.push_back(trackID); } size_t GetNumberOfHits(Int_t volID = -1) const; - inline Int_t GetTrackID() const { return fTrack_ID; } - inline Int_t GetParentID() const { return fParent_ID; } + size_t GetNumberOfPhysicalHits(Int_t volID = -1) const; + inline Int_t GetTrackID() const { return fTrackID; } + inline Int_t GetParentID() const { return fParentID; } inline TString GetParticleName() const { return fParticleName; } - EColor GetParticleColor() const; - inline Double_t GetGlobalTime() const { return fGlobalTimestamp; } - inline Double_t GetTrackTimeLength() const { return fTrackTimestamp; } - inline Double_t GetKineticEnergy() const { return fKineticEnergy; } + inline Double_t GetTimeLength() const { return fTimeLength; } + inline Double_t GetInitialKineticEnergy() const { return fInitialKineticEnergy; } inline Double_t GetTotalDepositedEnergy() const { return fHits.GetTotalDepositedEnergy(); } - inline TVector3 GetTrackOrigin() const { return fTrackOrigin; } - inline Int_t GetSubEventID() const { return fSubEventId; } + inline TVector3 GetInitialPosition() const { return fInitialPosition; } + inline Double_t GetWeight() const { return fWeight; } + inline Double_t GetEnergy() const { return fHits.GetEnergy(); } + inline Double_t GetLength() const { return fLength; } + + TString GetInitialVolume() const; + TString GetFinalVolume() const; + + inline std::vector GetSecondaryTrackIDs() const { return fSecondaryTrackIDs; } + std::vector GetSecondaryTracks() const; + inline std::vector GetChildrenTracks() const { return GetSecondaryTracks(); } + + TRestGeant4Track* GetParentTrack() const; + + inline TVector3 GetTrackOrigin() const { return GetInitialPosition(); } + + EColor GetParticleColor() const; inline Double_t GetEnergyInVolume(Int_t volID) const { return fHits.GetEnergyInVolume(volID); } inline TVector3 GetMeanPositionInVolume(Int_t volID) const { @@ -93,33 +107,6 @@ class TRestGeant4Track : public TObject { return fHits.GetLastPositionInVolume(volID); } - void SetSubEventID(Int_t id) { fSubEventId = id; } - void SetTrackID(Int_t id) { fTrack_ID = id; } - void SetParentID(Int_t id) { fParent_ID = id; } - // void SetParticleID( Int_t id ) { fParticle_ID = id; } - void SetParticleName(const TString& particleName) { fParticleName = particleName; } - void SetGlobalTrackTime(Double_t time) { fGlobalTimestamp = time; } - void SetTrackTimeLength(Double_t time) { fTrackTimestamp = time; } - void SetKineticEnergy(Double_t en) { fKineticEnergy = en; } - void SetTrackOrigin(const TVector3& pos) { fTrackOrigin = pos; } - void SetTrackOrigin(Double_t x, Double_t y, Double_t z) { fTrackOrigin.SetXYZ(x, y, z); } - - void AddG4Hit(const TVector3& pos, Double_t en, Double_t hit_global_time, Int_t pcs, Int_t vol, - Double_t eKin, const TVector3& momentumDirection) { - fHits.AddG4Hit(pos, en, hit_global_time, pcs, vol, eKin, momentumDirection); - } - - Double_t GetTrackLength() const; - - inline static Double_t GetDistance(const TVector3& v1, const TVector3& v2) { - return TMath::Sqrt((v1.X() - v2.X()) * (v1.X() - v2.X()) + (v1.Y() - v2.Y()) * (v1.Y() - v2.Y()) + - (v1.Z() - v2.Z()) * (v1.Z() - v2.Z())); - } - - inline void RemoveHits() { fHits.RemoveHits(); } - - const TRestGeant4Metadata* GetGeant4Metadata() const; - Int_t GetProcessID(const TString& processName) const; TString GetProcessName(Int_t id) const; @@ -131,15 +118,28 @@ class TRestGeant4Track : public TObject { return ContainsProcessInVolume(processName, -1); } + Double_t GetEnergyInVolume(const TString& volumeName, bool children = false) const; + + TString GetLastProcessName() const; + /// Prints the track information. N number of hits to print, 0 = all void PrintTrack(size_t maxHits = 0) const; + inline void RemoveHits() { fHits.RemoveHits(); } + // Constructor TRestGeant4Track(); + // Destructor virtual ~TRestGeant4Track(); - ClassDef(TRestGeant4Track, 4); // REST event superclass + ClassDef(TRestGeant4Track, 5); // REST event superclass + + // restG4 + public: + explicit TRestGeant4Track(const G4Track*); //! + void UpdateTrack(const G4Track*); //! + void InsertStep(const G4Step*); //! }; #endif diff --git a/inc/TRestGeant4VetoAnalysisProcess.h b/inc/TRestGeant4VetoAnalysisProcess.h index 0771be6..359f254 100644 --- a/inc/TRestGeant4VetoAnalysisProcess.h +++ b/inc/TRestGeant4VetoAnalysisProcess.h @@ -23,58 +23,56 @@ #ifndef RestCore_TRestGeant4VetoAnalysisProcess #define RestCore_TRestGeant4VetoAnalysisProcess -#include -#include +#include -#include "TRestEventProcess.h" - -using namespace std; +#include "TRestGeant4Event.h" +#include "TRestGeant4Metadata.h" class TRestGeant4VetoAnalysisProcess : public TRestEventProcess { private: - /// A pointer to the specific TRestGeant4Event input - TRestGeant4Event* fInputG4Event; //! - /// A pointer to the specific TRestGeant4Event output - TRestGeant4Event* fOutputG4Event; //! - /// A pointer to the simulation metadata information accessible to TRestRun - TRestGeant4Metadata* fG4Metadata; //! - - std::vector fVetoVolumeIds; //! - std::string fVetoKeyword = ""; //! - std::vector fVetoGroupKeywords; //! - std::map> fVetoGroupVolumeNames; //! - std::vector fQuenchingFactors; //! + TString fVetoVolumesExpression; + TString fVetoDetectorsExpression; + double fVetoDetectorOffsetSize = 0; + double fVetoLightAttenuation = 0; + double fVetoQuenchingFactor = 1.0; + + public: + inline TString GetVetoVolumesExpression() const { return fVetoVolumesExpression; } + inline TString GetVetoDetectorExpression() const { return fVetoDetectorsExpression; } + inline double GetVetoDetectorOffsetSize() const { return fVetoDetectorOffsetSize; } + inline double GetVetoLightAttenuation() const { return fVetoLightAttenuation; } + inline double GetVetoQuenchingFactor() const { return fVetoQuenchingFactor; } + + inline void SetVetoVolumesExpression(const TString& expression) { fVetoVolumesExpression = expression; } + inline void SetVetoDetectorsExpression(const TString& expression) { + fVetoDetectorsExpression = expression; + } + inline void SetVetoDetectorOffsetSize(double offset) { fVetoDetectorOffsetSize = offset; } + inline void SetVetoLightAttenuation(double attenuation) { fVetoLightAttenuation = attenuation; } + inline void SetVetoQuenchingFactor(double quenchingFactor) { fVetoQuenchingFactor = quenchingFactor; } + + private: + TRestGeant4Event* fInputEvent = nullptr; //! + TRestGeant4Event* fOutputEvent = nullptr; //! + const TRestGeant4Metadata* fGeant4Metadata = nullptr; //! + + std::vector fVetoVolumes; + std::vector fVetoDetectorVolumes; + std::map fVetoDetectorBoundaryPosition; + std::map fVetoDetectorBoundaryDirection; void InitFromConfigFile() override; void Initialize() override; void LoadDefaultConfig(); - // clean std::string (https://stackoverflow.com/questions/216823/whats-the-best-way-to-trim-stdstring) - inline std::string& rtrim(std::string& s, const char* t = " \t\n\r\f\v") { - s.erase(s.find_last_not_of(t) + 1); - return s; - } - // trim from beginning of std::string (left) - inline std::string& ltrim(std::string& s, const char* t = " \t\n\r\f\v") { - s.erase(0, s.find_first_not_of(t)); - return s; - } - // trim from both ends of std::string (right then left) - inline std::string& trim(std::string& s, const char* t = " \t\n\r\f\v") { return ltrim(rtrim(s, t), t); } - - // final clean std::string: trim and UPPER - inline std::string& clean_string(std::string& s) { - s = trim(s); - std::transform(s.begin(), s.end(), s.begin(), ::tolower); - return s; - } - - protected: - // add here the members of your event process - public: - any GetInputEvent() const override { return fInputG4Event; } - any GetOutputEvent() const override { return fOutputG4Event; } + any GetInputEvent() const override { return fInputEvent; } + any GetOutputEvent() const override { return fOutputEvent; } + + inline void SetGeant4Metadata(const TRestGeant4Metadata* metadata) { + fGeant4Metadata = metadata; + } // TODO: We should not need this! but `GetMetadata()` is not working early in the + // processing (look at the tests for more details) void InitProcess() override; TRestEvent* ProcessEvent(TRestEvent* inputEvent) override; @@ -82,52 +80,14 @@ class TRestGeant4VetoAnalysisProcess : public TRestEventProcess { void LoadConfig(const std::string& configFilename, const std::string& name = ""); - /// It prints out the process parameters stored in the metadata structure - void PrintMetadata() override { - BeginPrintProcess(); - - RESTDebug << "VETO KEYWORD: " << fVetoKeyword << RESTendl; - RESTDebug << RESTendl; - - RESTDebug << "VETO GROUP KEYWORDS:" << RESTendl; - for (unsigned int i = 0; i < fVetoGroupKeywords.size(); i++) { - RESTDebug << "\t" << fVetoGroupKeywords[i] << RESTendl; - } - RESTDebug << RESTendl; - - RESTDebug << "Found " << fVetoVolumeIds.size() << " veto volumes:" << RESTendl; - for (unsigned int i = 0; i < fVetoVolumeIds.size(); i++) { - RESTDebug << "\t" << fG4Metadata->GetActiveVolumeName(fVetoVolumeIds[i]) << RESTendl; - } - RESTDebug << RESTendl; - - RESTDebug << "GROUPS:" << RESTendl; - for (const auto& pair : fVetoGroupVolumeNames) { - RESTDebug << "GROUP " << pair.first << " (" << pair.second.size() << " volumes):\n"; - for (unsigned int i = 0; i < pair.second.size(); i++) { - RESTDebug << "\t" << pair.second[i] << RESTendl; - } - } - RESTDebug << RESTendl; - - RESTDebug << "QUENCHING FACTORS (" << fQuenchingFactors.size() << " Total)" << RESTendl; - for (unsigned int i = 0; i < fQuenchingFactors.size(); i++) { - RESTDebug << "\t" << fQuenchingFactors[i] << RESTendl; - } - RESTDebug << RESTendl; - - EndPrintProcess(); - } + void PrintMetadata() override; - /// Returns a new instance of this class - TRestEventProcess* Maker() { return new TRestGeant4VetoAnalysisProcess; } - /// Returns the name of this process const char* GetProcessName() const override { return "geant4VetoAnalysis"; } TRestGeant4VetoAnalysisProcess(); TRestGeant4VetoAnalysisProcess(const char* configFilename); ~TRestGeant4VetoAnalysisProcess(); - ClassDefOverride(TRestGeant4VetoAnalysisProcess, 1); + ClassDefOverride(TRestGeant4VetoAnalysisProcess, 2); }; #endif // RestCore_TRestGeant4VetoAnalysisProcess diff --git a/src/TRestGeant4Event.cxx b/src/TRestGeant4Event.cxx index 3daf0a7..80a1cb0 100644 --- a/src/TRestGeant4Event.cxx +++ b/src/TRestGeant4Event.cxx @@ -43,13 +43,9 @@ TRestGeant4Event::~TRestGeant4Event() { void TRestGeant4Event::Initialize() { TRestEvent::Initialize(); - fPrimaryParticleName.clear(); - fPrimaryEventDirection.clear(); - fPrimaryEventEnergy.clear(); - fPrimaryEventOrigin.SetXYZ(0, 0, 0); + fPrimaryParticleNames.clear(); - fTrack.clear(); - fNTracks = 0; + fTracks.clear(); // ClearVolumes(); fXZHitGraph = nullptr; @@ -76,7 +72,6 @@ void TRestGeant4Event::Initialize() { fTotalDepositedEnergy = 0; fSensitiveVolumeEnergy = 0; - fMaxSubEventID = 0; fMinX = 1e20; fMaxX = -1e20; @@ -109,29 +104,6 @@ void TRestGeant4Event::AddEnergyDepositToVolume(Int_t volID, Double_t eDep) { fVolumeDepositedEnergy[volID] += eDep; } -void TRestGeant4Event::SetTrackSubEventID(Int_t n, Int_t id) { - fTrack[n].SetSubEventID(id); - if (fMaxSubEventID < id) fMaxSubEventID = id; -} - -void TRestGeant4Event::AddTrack(TRestGeant4Track trk) { - fTrack.push_back(trk); - fNTracks = fTrack.size(); - fTotalDepositedEnergy += trk.GetTotalDepositedEnergy(); - for (int n = 0; n < GetNumberOfActiveVolumes(); n++) - fVolumeDepositedEnergy[n] += trk.GetEnergyInVolume(n); -} - -Double_t TRestGeant4Event::GetTotalDepositedEnergyFromTracks() const { - Double_t eDep = 0; - - for (int tk = 0; tk < GetNumberOfTracks(); tk++) { - eDep += GetTrack(tk).GetTotalDepositedEnergy(); - } - - return eDep; -} - TVector3 TRestGeant4Event::GetMeanPositionInVolume(Int_t volID) const { TVector3 pos; Double_t eDep = 0; @@ -216,10 +188,16 @@ TVector3 TRestGeant4Event::GetLastPositionInVolume(Int_t volID) const { return {nan, nan, nan}; } -TRestGeant4Track* TRestGeant4Event::GetTrackByID(int id) { - for (int i = 0; i < fNTracks; i++) - if (fTrack[i].GetTrackID() == id) return &fTrack[i]; - return nullptr; +TRestGeant4Track* TRestGeant4Event::GetTrackByID(Int_t trackID) const { + TRestGeant4Track* result = nullptr; + if (fTrackIDToTrackIndex.count(trackID) > 0) { + result = const_cast(&fTracks[fTrackIDToTrackIndex.at(trackID)]); + if (result == nullptr || result->GetTrackID() != trackID) { + cerr << "TRestGeant4Event::GetTrackByID - ERROR: trackIDToTrackIndex map is corrupted" << endl; + exit(1); + } + } + return result; } /////////////////////////////////////////////// @@ -227,12 +205,25 @@ TRestGeant4Track* TRestGeant4Event::GetTrackByID(int id) { /// a specific volume is given as argument only the hits of that specific volume /// will be counted. /// -Int_t TRestGeant4Event::GetNumberOfHits(Int_t volID) const { - Int_t hits = 0; - for (int i = 0; i < fNTracks; i++) { - hits += GetTrack(i).GetNumberOfHits(volID); +size_t TRestGeant4Event::GetNumberOfHits(Int_t volID) const { + size_t numberOfHits = 0; + for (const auto& track : fTracks) { + numberOfHits += track.GetNumberOfHits(volID); } - return hits; + return numberOfHits; +} + +/////////////////////////////////////////////// +/// \brief Function that returns the total number of hits with energy > 0 in the Geant4 event. If +/// a specific volume is given as argument only the hits of that specific volume +/// will be counted. +/// +size_t TRestGeant4Event::GetNumberOfPhysicalHits(Int_t volID) const { + size_t numberOfHits = 0; + for (const auto& track : fTracks) { + numberOfHits += track.GetNumberOfPhysicalHits(volID); + } + return numberOfHits; } /////////////////////////////////////////////// @@ -242,7 +233,7 @@ Int_t TRestGeant4Event::GetNumberOfHits(Int_t volID) const { /// TRestHits TRestGeant4Event::GetHits(Int_t volID) const { TRestHits hits; - for (int t = 0; t < fNTracks; t++) { + for (int t = 0; t < GetNumberOfTracks(); t++) { const auto& g4Hits = GetTrack(t).GetHits(); for (int n = 0; n < g4Hits.GetNumberOfHits(); n++) { if (volID != -1 && g4Hits.GetVolumeId(n) != volID) continue; @@ -336,7 +327,7 @@ void TRestGeant4Event::SetBoundaries() { } /* {{{ Get{XY,YZ,XZ}MultiGraph methods */ -TMultiGraph* TRestGeant4Event::GetXYMultiGraph(Int_t gridElement, std::vector pcsList, +TMultiGraph* TRestGeant4Event::GetXYMultiGraph(Int_t gridElement, vector pcsList, Double_t minPointSize, Double_t maxPointSize) { if (fXYHitGraph) { delete[] fXYHitGraph; @@ -428,7 +419,7 @@ TMultiGraph* TRestGeant4Event::GetXYMultiGraph(Int_t gridElement, std::vector pcsList, +TMultiGraph* TRestGeant4Event::GetYZMultiGraph(Int_t gridElement, vector pcsList, Double_t minPointSize, Double_t maxPointSize) { if (fYZHitGraph) { delete[] fYZHitGraph; @@ -520,7 +511,7 @@ TMultiGraph* TRestGeant4Event::GetYZMultiGraph(Int_t gridElement, std::vector pcsList, +TMultiGraph* TRestGeant4Event::GetXZMultiGraph(Int_t gridElement, vector pcsList, Double_t minPointSize, Double_t maxPointSize) { if (fXZHitGraph) { delete[] fXZHitGraph; @@ -614,7 +605,7 @@ TMultiGraph* TRestGeant4Event::GetXZMultiGraph(Int_t gridElement, std::vector optList) { +TH2F* TRestGeant4Event::GetXYHistogram(Int_t gridElement, vector optList) { if (fXYHisto) { delete fXYHisto; fXYHisto = nullptr; @@ -673,7 +664,7 @@ TH2F* TRestGeant4Event::GetXYHistogram(Int_t gridElement, std::vector o return fXYHisto; } -TH2F* TRestGeant4Event::GetXZHistogram(Int_t gridElement, std::vector optList) { +TH2F* TRestGeant4Event::GetXZHistogram(Int_t gridElement, vector optList) { if (fXZHisto) { delete fXZHisto; fXZHisto = nullptr; @@ -732,7 +723,7 @@ TH2F* TRestGeant4Event::GetXZHistogram(Int_t gridElement, std::vector o return fXZHisto; } -TH2F* TRestGeant4Event::GetYZHistogram(Int_t gridElement, std::vector optList) { +TH2F* TRestGeant4Event::GetYZHistogram(Int_t gridElement, vector optList) { if (fYZHisto) { delete fYZHisto; fYZHisto = nullptr; @@ -792,7 +783,7 @@ TH2F* TRestGeant4Event::GetYZHistogram(Int_t gridElement, std::vector o } /* }}} */ -TH1D* TRestGeant4Event::GetXHistogram(Int_t gridElement, std::vector optList) { +TH1D* TRestGeant4Event::GetXHistogram(Int_t gridElement, vector optList) { if (fXHisto) { delete fXHisto; fXHisto = nullptr; @@ -843,7 +834,7 @@ TH1D* TRestGeant4Event::GetXHistogram(Int_t gridElement, std::vector op return fXHisto; } -TH1D* TRestGeant4Event::GetZHistogram(Int_t gridElement, std::vector optList) { +TH1D* TRestGeant4Event::GetZHistogram(Int_t gridElement, vector optList) { if (fZHisto) { delete fZHisto; fZHisto = nullptr; @@ -894,7 +885,7 @@ TH1D* TRestGeant4Event::GetZHistogram(Int_t gridElement, std::vector op return fZHisto; } -TH1D* TRestGeant4Event::GetYHistogram(Int_t gridElement, std::vector optList) { +TH1D* TRestGeant4Event::GetYHistogram(Int_t gridElement, vector optList) { if (fYHisto) { delete fYHisto; fYHisto = nullptr; @@ -1025,7 +1016,7 @@ TPad* TRestGeant4Event::DrawEvent(const TString& option, Bool_t autoBoundaries) if (optList[n] == "print") this->PrintEvent(); } - optList.erase(std::remove(optList.begin(), optList.end(), "print"), optList.end()); + optList.erase(remove(optList.begin(), optList.end(), "print"), optList.end()); unsigned int nPlots = optList.size(); @@ -1126,49 +1117,37 @@ void TRestGeant4Event::PrintActiveVolumes() const { void TRestGeant4Event::PrintEvent(int maxTracks, int maxHits) const { TRestEvent::PrintEvent(); - cout.precision(4); + cout << "- Total deposited energy: " << ToEnergyString(fTotalDepositedEnergy) << endl; + cout << "- Sensitive detectors total energy: " << ToEnergyString(fSensitiveVolumeEnergy) << endl; - cout << "Total energy : " << fTotalDepositedEnergy << " keV" << endl; - cout << "Sensitive volume energy : " << fSensitiveVolumeEnergy << " keV" << endl; - cout << "Source origin : (" << fPrimaryEventOrigin.X() << "," << fPrimaryEventOrigin.Y() << "," - << fPrimaryEventOrigin.Z() << ") mm" << endl; + cout << "- Primary source position: " << VectorToString(fPrimaryPosition) << "{} mm" << endl; - for (int n = 0; n < GetNumberOfPrimaries(); n++) { - const auto dir = &fPrimaryEventDirection[n]; - cout << "Source " << n << " Particle name : " << GetPrimaryEventParticleName(n) << endl; - cout << "Source " << n << " direction : (" << dir->X() << "," << dir->Y() << "," << dir->Z() << ")" + for (int i = 0; i < GetNumberOfPrimaries(); i++) { + const char* sourceNumberString = + GetNumberOfPrimaries() <= 1 ? "" + : TString::Format(" %d of %zu", i + 1, GetNumberOfPrimaries()).Data(); + cout << " - Source" << sourceNumberString << " primary particle: " << fPrimaryParticleNames[i] << endl; - cout << "Source " << n << " energy : " << fPrimaryEventEnergy[n] << " keV" << endl; + cout << " Direction: " << VectorToString(fPrimaryDirections[i]) << endl; + cout << " Energy: " << ToEnergyString(fPrimaryEnergies[i]) << endl; } - cout << "Number of active volumes : " << GetNumberOfActiveVolumes() << endl; - for (int i = 0; i < GetNumberOfActiveVolumes(); i++) { - if (isVolumeStored(i)) { - cout << "Active volume " << i << ":" - << " has been stored." << endl; - cout << "Total energy deposit in volume " << i << ":" - << " : " << fVolumeDepositedEnergy[i] << " keV" << endl; - } else - cout << "Active volume " << i << ":" - << " has not been stored" << endl; - } - - cout << "---------------------------------------------------------------------------" << endl; - cout << "Total number of tracks : " << fNTracks << endl; + cout << endl; + cout << "Total number of tracks: " << GetNumberOfTracks() << endl; int nTracks = GetNumberOfTracks(); - if (maxTracks > 0) { - nTracks = min(maxTracks, GetNumberOfTracks()); - cout << " Printing only the first " << nTracks << " tracks" << endl; + if (maxTracks > 0 && maxTracks < GetNumberOfTracks()) { + nTracks = min(maxTracks, int(GetNumberOfTracks())); + cout << "Printing only the first " << nTracks << " tracks" << endl; } - for (int n = 0; n < nTracks; n++) { - GetTrack(n).PrintTrack(maxHits); + for (int i = 0; i < nTracks; i++) { + GetTrack(i).PrintTrack(maxHits); } } Bool_t TRestGeant4Event::ContainsProcessInVolume(Int_t processID, Int_t volumeID) const { - for (const auto& track : fTrack) { + for (const auto& track : fTracks) { if (track.ContainsProcessInVolume(processID, volumeID)) { return true; } @@ -1182,7 +1161,7 @@ Bool_t TRestGeant4Event::ContainsProcessInVolume(const TString& processName, Int return false; } const auto& processID = metadata->GetGeant4PhysicsInfo().GetProcessID(processName); - for (const auto& track : fTrack) { + for (const auto& track : fTracks) { if (track.ContainsProcessInVolume(processID, volumeID)) { return true; } @@ -1191,7 +1170,7 @@ Bool_t TRestGeant4Event::ContainsProcessInVolume(const TString& processName, Int } Bool_t TRestGeant4Event::ContainsParticle(const TString& particleName) const { - for (const auto& track : fTrack) { + for (const auto& track : fTracks) { if (track.GetParticleName() == particleName) { return true; } @@ -1200,7 +1179,7 @@ Bool_t TRestGeant4Event::ContainsParticle(const TString& particleName) const { } Bool_t TRestGeant4Event::ContainsParticleInVolume(const TString& particleName, Int_t volumeID) const { - for (const auto& track : fTrack) { + for (const auto& track : fTracks) { if (track.GetParticleName() != particleName) { continue; } @@ -1221,7 +1200,90 @@ void TRestGeant4Event::InitializeReferences(TRestRun* run) { This introduces overhead to event loading, but hopefully its small enough. If this is a problem, we could rework this approach */ - for (auto& track : fTrack) { + for (auto& track : fTracks) { track.SetEvent(this); } } + +const set TRestGeant4Event::GetUniqueParticles() const { + set result; + for (const auto& track : fTracks) { + result.insert(track.GetParticleName().Data()); + } + return result; +} + +map> TRestGeant4Event::GetEnergyInVolumePerProcessMap() const { + map> result; + for (const auto& [volume, particleProcessMap] : fEnergyInVolumePerParticlePerProcess) { + for (const auto& [particle, processMap] : particleProcessMap) { + for (const auto& [process, energy] : processMap) { + result[volume][process] += energy; + } + } + } + return result; +} + +map> TRestGeant4Event::GetEnergyInVolumePerParticleMap() const { + map> result; + for (const auto& [volume, particleProcessMap] : fEnergyInVolumePerParticlePerProcess) { + for (const auto& [particle, processMap] : particleProcessMap) { + for (const auto& [process, energy] : processMap) { + result[volume][particle] += energy; + } + } + } + return result; +} + +map TRestGeant4Event::GetEnergyPerProcessMap() const { + map result; + for (const auto& [volume, particleProcessMap] : fEnergyInVolumePerParticlePerProcess) { + for (const auto& [particle, processMap] : particleProcessMap) { + for (const auto& [process, energy] : processMap) { + result[process] += energy; + } + } + } + return result; +} + +map TRestGeant4Event::GetEnergyPerParticleMap() const { + map result; + for (const auto& [volume, particleProcessMap] : fEnergyInVolumePerParticlePerProcess) { + for (const auto& [particle, processMap] : particleProcessMap) { + for (const auto& [process, energy] : processMap) { + result[particle] += energy; + } + } + } + return result; +} + +map TRestGeant4Event::GetEnergyInVolumeMap() const { + map result; + for (const auto& [volume, particleProcessMap] : fEnergyInVolumePerParticlePerProcess) { + for (const auto& [particle, processMap] : particleProcessMap) { + for (const auto& [process, energy] : processMap) { + result[volume] += energy; + } + } + } + return result; +} + +map>> TRestGeant4Event::GetEnergyInVolumePerParticlePerProcessMap() + const { + return fEnergyInVolumePerParticlePerProcess; +} + +void TRestGeant4Event::AddEnergyInVolumeForParticleForProcess(Double_t energy, const string& volumeName, + const string& particleName, + const string& processName) { + if (energy <= 0) { + return; + } + fEnergyInVolumePerParticlePerProcess[volumeName][particleName][processName] += energy; + fTotalDepositedEnergy += energy; +} diff --git a/src/TRestGeant4EventViewer.cxx b/src/TRestGeant4EventViewer.cxx index 1d4e29f..095f0c0 100644 --- a/src/TRestGeant4EventViewer.cxx +++ b/src/TRestGeant4EventViewer.cxx @@ -14,7 +14,8 @@ #include "TRestGeant4EventViewer.h" -#include "TRestStringOutput.h" +#include +#include using namespace std; @@ -35,105 +36,199 @@ void TRestGeant4EventViewer::Initialize() { } void TRestGeant4EventViewer::DeleteCurrentEvent() { - // cout<<"Removing event"< 10 ? 10 : width); + width = (width < 1 ? 1 : width); + config.fLineWidth = TMath::Log10(track.GetInitialKineticEnergy() / 10); - Double_t eDepMin = 1.e6; - Double_t eDepMax = 0; - Double_t totalEDep = 0; + // line style - for (int i = 0; i < fG4Event->GetNumberOfTracks(); i++) { - const auto& g4Track = fG4Event->GetTrack(i); - const auto& g4Hits = g4Track.GetHits(); - Int_t nHits = g4Track.GetNumberOfHits(); - for (int j = 0; j < nHits; j++) { - Double_t eDep = g4Hits.GetEnergy(j); - if (eDep > eDepMax) eDepMax = eDep; - if (eDep < eDepMin) eDepMin = eDep; + if (track.GetCreatorProcess() == "nCapture") { + config.fLineStyle = 2; + } - // cout<<" track "<AddLine({static_cast(GEOM_SCALE * hits.GetPosition(i).x()), + static_cast(GEOM_SCALE * hits.GetPosition(i).y()), + static_cast(GEOM_SCALE * hits.GetPosition(i).z())}, // + {static_cast(GEOM_SCALE * hits.GetPosition(i + 1).x()), + static_cast(GEOM_SCALE * hits.GetPosition(i + 1).y()), + static_cast(GEOM_SCALE * hits.GetPosition(i + 1).z())}); + + const auto config = GetTrackVisualConfiguration(track); + lineSet->SetMainColor(config.fColor); + lineSet->SetLineColor(config.fColor); + lineSet->SetLineStyle(config.fLineStyle); + lineSet->SetLineWidth(config.fLineWidth); } - cout << "TRestGeant4EventViewer::AddEvent. Total EDep " << totalEDep << endl; - - Double_t slope = (fMaxRadius - fMinRadius) / (eDepMax - eDepMin); - Double_t bias = fMinRadius - slope * eDepMin; - - Int_t textAdded = 0; - for (int trkID = 1; trkID < fG4Event->GetNumberOfTracks() + 1; trkID++) { - const auto& g4Track = fG4Event->GetTrackByID(trkID); - - Int_t parentID = g4Track->GetParentID(); - TVector3 origin = g4Track->GetTrackOrigin(); - - // Building track name - Double_t eKin = g4Track->GetKineticEnergy(); - // cout << "eKin: " << eKin << endl; - TString ptlName = g4Track->GetParticleName(); - // cout << "ptlName: " << ptlName << endl; - char pcleStr[64]; - sprintf(pcleStr, "Track ID : %d %s (%6.2lf keV)", trkID, ptlName.Data(), eKin); - - char markerStr[256]; - sprintf(markerStr, " %s (%6.2lf keV). Position (%4.1lf,%4.1lf,%4.1lf)", ptlName.Data(), eKin, - origin.X(), origin.Y(), origin.Z()); - - if (parentID == 0) { - // cout << " Parent ID : 0" << endl; - char evInfoStr[256]; - sprintf(evInfoStr, "%s. EventID = %d at position (%4.2lf, %4.2lf, %4.2lf) mm", ptlName.Data(), - fG4Event->GetID(), origin.X(), origin.Y(), origin.Z()); - this->AddParentTrack(trkID, origin, pcleStr); - if (fG4Event->GetSubID() == 0) - this->AddText(evInfoStr, origin); - else if (!textAdded) { - textAdded = 1; - sprintf(evInfoStr, "%s. EventID = %d at position (%4.2lf, %4.2lf, %4.2lf) mm", - fG4Event->GetSubEventTag().Data(), fG4Event->GetID(), origin.X(), origin.Y(), - origin.Z()); - this->AddText(evInfoStr, origin); - } - } else { - // cout << "Adding track : " << trkID << " parent : " << parentID << " pt: " << pcleStr << endl; - this->AddTrack(trkID, parentID, origin, pcleStr); - } + return lineSet; +} + +struct HitsVisualConfiguration { + Color_t fColor = kBlack; + Style_t fMarkerStyle = 1; + Size_t fMarkerSize = 1.0; +}; + +HitsVisualConfiguration GetHitsVisualConfiguration(const TString& processNameOrType) { + auto config = HitsVisualConfiguration(); + // based on particle type + if (processNameOrType.EqualTo("Electromagnetic")) { + config.fColor = kRed; + } else if (processNameOrType.EqualTo("Init")) { + // custom process (not Geant4) + config.fColor = kWhite; + } + + // based on particle name + if (processNameOrType.EqualTo("Transportation")) { + config.fColor = kWhite; + } else if (processNameOrType.EqualTo("Init")) { + // custom process (not Geant4) + config.fColor = kWhite; + } else if (processNameOrType.EqualTo("hadElastic")) { + config.fColor = kOrange; + } else if (processNameOrType.EqualTo("neutronInelastic")) { + config.fColor = kGreen; + config.fMarkerStyle = 4; + } else if (processNameOrType.EqualTo("nCapture")) { + config.fColor = kBlue; + config.fMarkerStyle = 2; + config.fMarkerSize = 4.0; + } + return config; +} + +void TRestGeant4EventViewer::AddEvent(TRestEvent* event) { + DeleteCurrentEvent(); - // cout << "Adding marker" << endl; - this->AddMarker(trkID, origin, markerStr); + fG4Event = (TRestGeant4Event*)event; + fG4Metadata = fG4Event->GetGeant4Metadata(); + if (fG4Metadata == nullptr) { + cerr << "TRestGeant4EventViewer::Initialize. No TRestGeant4Metadata found in TRestGeant4Event" + << endl; + exit(1); + } + + size_t trackCounter = 0; + size_t hitsCounter = 0; - const auto& g4Hits = g4Track->GetHits(); - Int_t nHits = g4Track->GetNumberOfHits(); + map linesSet; + map hitsPoints; + map hitsType; - // cout << "Adding hits" << endl; - for (int i = 0; i < nHits; i++) { - Double_t x = g4Hits.GetX(i); - Double_t y = g4Hits.GetY(i); - Double_t z = g4Hits.GetZ(i); + auto trackList = new TEveElementList("Tracks"); + gEve->AddElement(trackList); + auto hitsList = new TEveElementList("Hits"); + gEve->AddElement(hitsList); - this->NextTrackVertex(trkID, TVector3(x, y, z)); + const auto& physicsInfo = fG4Metadata->GetGeant4PhysicsInfo(); - Double_t eDep = g4Hits.GetEnergy(i); + for (const auto& track : fG4Event->GetTracks()) { + if (track.GetInitialKineticEnergy() < 1.0 || track.GetLength() < 0.1) { + continue; + } - if (eDep > 0) { - Double_t radius = slope * eDep + bias; - AddSphericalHit(x, y, z, radius, eDep); + auto line = GetTrackEveDrawable(track); + linesSet[track.GetTrackID()] = line; + TEveElement* parentLine = trackList; + if (linesSet.count(track.GetParentID())) { + parentLine = linesSet.at(track.GetParentID()); + } + gEve->AddElement(line, parentLine); + trackCounter++; + + const auto& hits = track.GetHits(); + for (int i = 0; i < hits.GetNumberOfHits(); i++) { + const auto& processName = physicsInfo.GetProcessName(hits.GetProcess(i)); + const auto& processType = physicsInfo.GetProcessType(processName); + const auto& position = hits.GetPosition(i); + + if (hitsType.count(processType) == 0) { + hitsType[processType] = new TEveElementList(processType); + gEve->AddElement(hitsType[processType], hitsList); + } + if (hitsPoints.count(processName) == 0) { + hitsPoints[processName] = new TEvePointSet(processName); + auto hitPoints = hitsPoints.at(processName); + auto hitsVisualConfig = GetHitsVisualConfiguration(processName); + hitPoints->SetMarkerColor(hitsVisualConfig.fColor); + hitPoints->SetMarkerStyle(hitsVisualConfig.fMarkerStyle); + hitPoints->SetMarkerSize(hitsVisualConfig.fMarkerSize); + + gEve->AddElement(hitPoints, hitsType[processType]); } + hitsPoints.at(processName) + ->SetNextPoint(GEOM_SCALE * position.X(), GEOM_SCALE * position.Y(), + GEOM_SCALE * position.Z()); + hitsCounter++; } } - // cout << "Updating" << endl; + + // Add event text + const auto& firstTrack = fG4Event->GetTracks().front(); + TVector3 position = {firstTrack.GetInitialPosition().X(), firstTrack.GetInitialPosition().Y(), + firstTrack.GetInitialPosition().Z()}; + AddText( + TString::Format( + "Event ID: %d%s | Primary origin: (%4.2lf, %4.2lf, %4.2lf) mm", fG4Event->GetID(), + (fG4Event->GetSubID() > 0 ? TString::Format(" (SubID: %d)", fG4Event->GetSubID()) : "").Data(), + position.X(), position.Y(), position.Z()), + position); Update(); } @@ -182,7 +277,8 @@ void TRestGeant4EventViewer::AddTrack(Int_t trkID, Int_t parentID, TVector3 from fHitConnectors[parentID]->AddElement(fHitConnectors[trkID]); else { RESTWarning << "Parent ID: " << parentID << " of track " << trkID << " was not found!" << RESTendl; - RESTWarning << "This might be solved by enabling TRestGeant4Metadata::fRegisterEmptyTracks" << RESTendl; + RESTWarning << "This might be solved by enabling TRestGeant4Metadata::fRegisterEmptyTracks" + << RESTendl; } } diff --git a/src/TRestGeant4GeometryInfo.cxx b/src/TRestGeant4GeometryInfo.cxx index dda4260..48dfca9 100644 --- a/src/TRestGeant4GeometryInfo.cxx +++ b/src/TRestGeant4GeometryInfo.cxx @@ -37,12 +37,6 @@ TString GetNodeAttribute(TXMLEngine xml, XMLNodePointer_t node, const TString& a void AddVolumesRecursively(vector* physicalNames, vector* logicalNames, const vector& children, map& nameTable, map>& childrenTable, const TString& name = "") { - /* - cout << "called AddVolumesRecursively with name: " << name << endl; - for (const auto& child : children) { - cout << "\t" << child << endl; - } - */ if (children.empty()) { physicalNames->push_back(name); const auto logicalVolumeName = nameTable[name]; @@ -117,16 +111,21 @@ void TRestGeant4GeometryInfo::PopulateFromGdml(const TString& gdmlFilename) { xml.FreeDoc(xmldoc); + // Checks + if (fGdmlNewPhysicalNames.empty()) { + cout << "TRestGeant4GeometryInfo::PopulateFromGdml - ERROR - No physical volumes have been added!" + << endl; + exit(1); + } // Find duplicates - size_t n = fGdmlNewPhysicalNames.size(); set s; for (const auto& name : fGdmlNewPhysicalNames) { s.insert(name); } - if (s.size() != n) { - cout - << "TRestGeant4GeometryInfo::PopulateFromGdml - Generated a duplicate name, please check assembly" - << endl; + if (s.size() != fGdmlNewPhysicalNames.size()) { + cout << "TRestGeant4GeometryInfo::PopulateFromGdml - ERROR - Generated a duplicate name, please " + "check assembly" + << endl; exit(1); } } @@ -149,28 +148,6 @@ TString TRestGeant4GeometryInfo::GetGeant4PhysicalNameFromAlternativeName( return ""; } -/* - * DO NOT REMOVE! -Int_t TRestGeant4GeometryInfo::GetIDFromVolumeName(const TString& volumeName) const { - for (int i = 0; i < fGdmlNewPhysicalNames.size(); i++) { - if (volumeName.EqualTo(fGdmlNewPhysicalNames[i])) { - return i; - } - } - - int i = 0; - for (const auto& physical : GetAllPhysicalVolumes()) { - if (volumeName.EqualTo(physical)) { - return i; - } - i++; - } - - cout << "TRestGeant4GeometryInfo::GetIDFromPhysicalName - ID not found for " << volumeName << endl; - exit(1); -} -*/ - template U GetOrDefaultMapValueFromKey(const map* pMap, const T& key) { if (pMap->count(key) > 0) { @@ -184,6 +161,12 @@ TString TRestGeant4GeometryInfo::GetVolumeFromID(Int_t id) const { } Int_t TRestGeant4GeometryInfo::GetIDFromVolume(const TString& volumeName) const { + if (fVolumeNameReverseMap.count(volumeName) == 0) { + // if we do not find the volume we return -1 instead of default (which is 0 and may be confusing) + cout << "TRestGeant4GeometryInfo::GetIDFromVolume - volume '" << volumeName << "' not found in store!" + << endl; + return -1; + } return GetOrDefaultMapValueFromKey(&fVolumeNameReverseMap, volumeName); } @@ -198,11 +181,15 @@ void TRestGeant4GeometryInfo::Print() const { const auto physicalVolumes = GetAllPhysicalVolumes(); cout << "Physical volumes (" << physicalVolumes.size() << "):" << endl; for (const auto& physical : GetAllPhysicalVolumes()) { - auto newName = GetAlternativeNameFromGeant4PhysicalName(physical); - const auto logical = fPhysicalToLogicalVolumeMap.at(physical); - cout << "\t- " << (newName == physical ? physical : newName + " (" + physical + ")") + auto geant4Name = GetGeant4PhysicalNameFromAlternativeName(physical); + const auto& logical = fPhysicalToLogicalVolumeMap.at(physical); + const auto& position = GetPosition(physical); + cout << "\t- " << (geant4Name == physical ? physical : physical + " (" + geant4Name + ")") + << " - ID: " << GetIDFromVolume(physical) << " - Logical: " << fPhysicalToLogicalVolumeMap.at(physical) - << " - Material: " << fLogicalToMaterialMap.at(logical) << endl; + << " - Material: " << fLogicalToMaterialMap.at(logical) // + << " - Position: (" << position.X() << ", " << position.Y() << ", " << position.Z() << ") mm" // + << endl; } const auto logicalVolumes = GetAllLogicalVolumes(); diff --git a/src/TRestGeant4Hits.cxx b/src/TRestGeant4Hits.cxx index ffe4d2e..6178fa8 100644 --- a/src/TRestGeant4Hits.cxx +++ b/src/TRestGeant4Hits.cxx @@ -17,62 +17,31 @@ #include "TRestGeant4Hits.h" +#include "TRestGeant4Event.h" +#include "TRestGeant4Track.h" + using namespace std; ClassImp(TRestGeant4Hits); -TRestGeant4Hits::TRestGeant4Hits() : TRestHits() { - // TRestGeant4Hits default constructor -} - -TRestGeant4Hits::~TRestGeant4Hits() { - // TRestGeant4Hits destructor -} - -void TRestGeant4Hits::AddG4Hit(TVector3 pos, Double_t en, Double_t hit_global_time, Int_t process, - Int_t volume, Double_t eKin, TVector3 momentumDirection) { - AddHit(pos, en, hit_global_time); - - fProcessID.Set(fNHits); - fProcessID[fNHits - 1] = process; - - fVolumeID.Set(fNHits); - fVolumeID[fNHits - 1] = volume; - - fKineticEnergy.Set(fNHits); - fKineticEnergy[fNHits - 1] = eKin; - - momentumDirection = momentumDirection.Unit(); - - Float_t x = momentumDirection.X(); - Float_t y = momentumDirection.Y(); - Float_t z = momentumDirection.Z(); - - fMomentumDirectionX.Set(fNHits); - fMomentumDirectionX[fNHits - 1] = x; +TRestGeant4Hits::TRestGeant4Hits() : TRestHits() {} - fMomentumDirectionY.Set(fNHits); - fMomentumDirectionY[fNHits - 1] = y; - - fMomentumDirectionZ.Set(fNHits); - fMomentumDirectionZ[fNHits - 1] = z; -} +TRestGeant4Hits::~TRestGeant4Hits() {} void TRestGeant4Hits::RemoveG4Hits() { RemoveHits(); - fProcessID.Set(0); - - fVolumeID.Set(0); - - fKineticEnergy.Set(0); + fProcessID.clear(); + fVolumeID.clear(); + fKineticEnergy.clear(); } Double_t TRestGeant4Hits::GetEnergyInVolume(Int_t volumeID) const { Double_t en = 0; - for (int n = 0; n < fNHits; n++) + for (int n = 0; n < fNHits; n++) { if (fVolumeID[n] == volumeID) en += GetEnergy(n); + } return en; } @@ -120,3 +89,17 @@ size_t TRestGeant4Hits::GetNumberOfHitsInVolume(Int_t volumeID) const { } return result; } + +TRestGeant4Metadata* TRestGeant4Hits::GetGeant4Metadata() const { + const TRestGeant4Event* event = nullptr; + if (fTrack != nullptr) { + event = fTrack->GetEvent(); + } else { + event = fEvent; + } + if (event == nullptr) { + cout << "null event" << endl; + return nullptr; + } + return const_cast(event->GetGeant4Metadata()); +} diff --git a/src/TRestGeant4Metadata.cxx b/src/TRestGeant4Metadata.cxx index 8b431ec..1c4f185 100644 --- a/src/TRestGeant4Metadata.cxx +++ b/src/TRestGeant4Metadata.cxx @@ -107,8 +107,8 @@ /// *restG4*, as important as the detector geometry, or the number of events to /// be simulated. /// -/// * **Nevents**: The number of primary particles to be generated. The number -/// of registered events might differ from *Nevents* due to storage definitions. +/// * **nEvents**: The number of primary particles to be generated. The number +/// of registered events might differ from *nEvents* due to storage definitions. /// The number of events registered could be even larger than the number of /// primaries, as for example when launching full decay chain simulations, where /// different isotope decays are stored in different events. @@ -159,8 +159,8 @@ /// /// /// -/// -/// +/// +/// /// /// /// @@ -283,8 +283,8 @@ /// \code /// // 3. geant4 internal /// -/// -/// +/// +/// /// /// \endcode /// @@ -331,7 +331,7 @@ /// \code /// // We launch pre-generated 136Xe NLDBD events /// -/// // The energyDist and angularDist +/// // The energy and angular /// // definitions here will be just /// // ignored with this type of source /// @@ -341,7 +341,7 @@ /// /// A source, or particle, associated kinetic energy is defined by using. /// \code -/// +/// /// \endcode /// /// There are different energy distribution types we can use to define the @@ -354,7 +354,7 @@ /// `energy="E"`, where `E` is the kinetic energy of the particle in keV. /// \code /// // A mono energetic particle with 10keV. -/// +/// /// \endcode /// /// * **flat**: All the particles from this source will be launched in a @@ -362,7 +362,7 @@ /// define the parameter `range="(Ei,Ef)"`, where `Ei` is the minimum energy and /// `Ef` is the maximum energy (in keV). \code /// // A random uniform generation between 1 and 10 keV -/// +/// /// \endcode /// /// * **log**: All the particles from this source will be launched in a @@ -370,12 +370,12 @@ /// define the parameter `range="(Ei,Ef)"`, where `Ei` is the minimum energy and /// `Ef` is the maximum energy (in keV). `Ei` needs to be > 0. \code /// // A random logarithmic generation between 1 and 10 keV -/// +/// /// \endcode /// /// * **TH1D**: It will use a TH1D histogram from a ROOT file with a user /// defined spectrum. It requires to define the parameters -/// `file="mySpectrum.root"` `spctName="histName"` and `range="(Ei,Ef)"`. The +/// `file="mySpectrum.root"` `name="histName"` and `range="(Ei,Ef)"`. The /// ROOT file should contain a TH1D histogram with name `histName`. Only the /// region of the spectrum inside the range `Ei-Ef` will be considered. If `range` /// is not specified inside the RML, the full TH1D range definition will be used. @@ -385,7 +385,7 @@ /// \code /// // A TH1D input spectrum to produce underground muons in the range /// // between 150 and 400 GeV -/// /// \endcode /// @@ -393,7 +393,7 @@ /// /// The momentum direction of a particle is specified by using. /// \code -/// +/// /// \endcode /// /// We can use different types of angular distributions to define the @@ -406,7 +406,7 @@ /// normalized). /// \code /// // Particles will be launched in the positive y-axis direction. -/// +/// /// \endcode /// /// * **isotropic**: The momentum direction of each particle will be @@ -415,7 +415,7 @@ /// volume will be considered. /// \code /// // Particles will be launched without preferred direction. -/// +/// /// \endcode /// /// * **backtoback**: The source momentum direction will be opposite to @@ -423,18 +423,18 @@ /// angular distribution type will be re-defined to isotropic. /// \code /// // Particles will be launched without preferred direction. -/// +/// /// \endcode /// -/// * **TH1D** : It will use a TH1D histogram from a ROOT file with a +/// * **TH1D**: It will use a TH1D histogram from a ROOT file with a /// user defined angular distribution. It requires to define the /// additional parameters as `file="mySpectrum.root"` and -/// `spctName="histName"`. The file we give should be stored in +/// `name="histName"`. The file we give should be stored in /// `"data/distributions/"` and contain a TH1D histogram with /// name "histName". /// \code /// // A TH1D input angular distribution used for cosmic rays -/// +/// /// \endcode /// /// ## 3. The storage section definition @@ -648,52 +648,16 @@ #include "TRestGeant4Metadata.h" +#include #include #include +#include #include -#include "TRestGDMLParser.h" +#include "TRestGeant4PrimaryGeneratorInfo.h" using namespace std; - -namespace g4_metadata_parameters { -string CleanString(string s) { - // transform the string to lowercase - transform(s.begin(), s.end(), s.begin(), ::tolower); - // this is a temporary fix, TH1D name comparison is being done elsewhere and giving error - if (s == "th1d") { - s = "TH1D"; - } - return s; -} - -map generator_types_map = { - {CleanString("custom"), generator_types::CUSTOM}, - {CleanString("volume"), generator_types::VOLUME}, - {CleanString("surface"), generator_types::SURFACE}, - {CleanString("point"), generator_types::POINT}, -}; - -std::map generator_shapes_map = { - {CleanString("gdml"), generator_shapes::GDML}, {CleanString("wall"), generator_shapes::WALL}, - {CleanString("circle"), generator_shapes::CIRCLE}, {CleanString("box"), generator_shapes::BOX}, - {CleanString("sphere"), generator_shapes::SPHERE}, {CleanString("cylinder"), generator_shapes::CYLINDER}, -}; - -map energy_dist_types_map = { - {CleanString("TH1D"), energy_dist_types::TH1D}, - {CleanString("mono"), energy_dist_types::MONO}, - {CleanString("flat"), energy_dist_types::FLAT}, - {CleanString("log"), energy_dist_types::LOG}, -}; - -map angular_dist_types_map = { - {CleanString("TH1D"), angular_dist_types::TH1D}, - {CleanString("isotropic"), angular_dist_types::ISOTROPIC}, - {CleanString("flux"), angular_dist_types::FLUX}, - {CleanString("backtoback"), angular_dist_types::BACK_TO_BACK}, -}; -} // namespace g4_metadata_parameters +using namespace TRestGeant4PrimaryGeneratorTypes; ClassImp(TRestGeant4Metadata); @@ -721,8 +685,6 @@ TRestGeant4Metadata::TRestGeant4Metadata(const char* configFilename, const strin Initialize(); LoadConfigFromFile(fConfigFileName, name); - - PrintMetadata(); } /////////////////////////////////////////////// @@ -741,13 +703,7 @@ void TRestGeant4Metadata::Initialize() { fActiveVolumes.clear(); fBiasingVolumes.clear(); - fNBiasingVolumes = 0; - RemoveParticleSources(); - - fSensitiveVolume = "gas"; - - fEnergyRangeStored.Set(0, 1.E20); } /////////////////////////////////////////////// @@ -799,11 +755,7 @@ void TRestGeant4Metadata::InitFromConfigFile() { if (nEventsString == PARAMETER_NOT_FOUND_STR) { nEventsString = GetParameter("Nevents"); // old name } - if (nEventsString == PARAMETER_NOT_FOUND_STR) { - cout << "\"nEvents\" parameter is not defined!" << endl; - exit(1); - } - fNEvents = StringToInteger(nEventsString); + fNEvents = nEventsString == PARAMETER_NOT_FOUND_STR ? 0 : StringToInteger(nEventsString); fSaveAllEvents = ToUpper(GetParameter("saveAllEvents", "false")) == "TRUE" || ToUpper(GetParameter("saveAllEvents", "off")) == "ON"; @@ -815,9 +767,8 @@ void TRestGeant4Metadata::InitFromConfigFile() { ToUpper(GetParameter("registerEmptyTracks", "off")) == "ON"; ReadGenerator(); - - ReadStorage(); - + // Detector (old storage) section is processed after initializing geometry info in Detector Construction + // This allows to use regular expression to match logical or physical volumes etc. ReadBiasing(); fMaxTargetStepSize = GetDblParameterWithUnits("maxTargetStepSize", -1); @@ -874,7 +825,7 @@ void TRestGeant4Metadata::ReadBiasing() { TString biasEnabled = GetFieldValue("value", biasingDefinition); TString biasType = GetFieldValue("type", biasingDefinition); - RESTDebug << "Bias : " << biasEnabled << RESTendl; + RESTDebug << "Bias: " << biasEnabled << RESTendl; if (biasEnabled == "on" || biasEnabled == "ON" || biasEnabled == "On" || biasEnabled == "oN") { RESTInfo << "Biasing is enabled" << RESTendl; @@ -883,7 +834,7 @@ void TRestGeant4Metadata::ReadBiasing() { Int_t n = 0; while (biasVolumeDefinition) { TRestGeant4BiasingVolume biasVolume; - RESTDebug << "Def : " << biasVolumeDefinition << RESTendl; + RESTDebug << "Def: " << biasVolumeDefinition << RESTendl; biasVolume.SetBiasingVolumePosition( Get3DVectorParameterWithUnits("position", biasVolumeDefinition)); @@ -922,8 +873,8 @@ void TRestGeant4Metadata::ReadBiasing() { /// /// ///3. geant4 internal /// /// /// \endcode @@ -940,17 +891,27 @@ void TRestGeant4Metadata::ReadGenerator() { TiXmlElement* generatorDefinition = GetElement("generator"); - fGenType = GetParameter("type", generatorDefinition, "volume"); - fGenShape = GetParameter("shape", generatorDefinition, "box"); - fGenFrom = GetParameter("from", generatorDefinition); - if (fGenFrom != PARAMETER_NOT_FOUND_STR) { - fGenShape = "gdml"; + fGeant4PrimaryGeneratorInfo.fSpatialGeneratorType = + SpatialGeneratorTypesToString(StringToSpatialGeneratorTypes(GetParameter( + "type", generatorDefinition, SpatialGeneratorTypesToString(SpatialGeneratorTypes::VOLUME)))); + fGeant4PrimaryGeneratorInfo.fSpatialGeneratorShape = + SpatialGeneratorShapesToString(StringToSpatialGeneratorShapes(GetParameter( + "shape", generatorDefinition, SpatialGeneratorShapesToString(SpatialGeneratorShapes::BOX)))); + fGeant4PrimaryGeneratorInfo.fSpatialGeneratorFrom = GetParameter("from", generatorDefinition); + if (fGeant4PrimaryGeneratorInfo.fSpatialGeneratorFrom != PARAMETER_NOT_FOUND_STR) { + fGeant4PrimaryGeneratorInfo.fSpatialGeneratorShape = + SpatialGeneratorShapesToString(SpatialGeneratorShapes::GDML); } - fGenSize = Get3DVectorParameterWithUnits("size", generatorDefinition); - fGenPosition = Get3DVectorParameterWithUnits("position", generatorDefinition); - fGenRotationAxis = StringTo3DVector(GetParameter("rotationAxis", generatorDefinition, "(0,0,1)")); - fGenRotationDegree = GetDblParameterWithUnits("rotationAngle", generatorDefinition); - fGenDensityFunction = GetParameter("densityFunc", generatorDefinition, "1"); + fGeant4PrimaryGeneratorInfo.fSpatialGeneratorSize = + Get3DVectorParameterWithUnits("size", generatorDefinition, {0, 0, 0}); + fGeant4PrimaryGeneratorInfo.fSpatialGeneratorPosition = + Get3DVectorParameterWithUnits("position", generatorDefinition, {0, 0, 0}); + fGeant4PrimaryGeneratorInfo.fSpatialGeneratorRotationAxis = + Get3DVectorParameterWithUnits("rotationAxis", generatorDefinition, {0, 0, 1}); + fGeant4PrimaryGeneratorInfo.fSpatialGeneratorRotationValue = + GetDblParameterWithUnits("rotationAngle", generatorDefinition, 0); + fGeant4PrimaryGeneratorInfo.fSpatialGeneratorSpatialDensityFunction = + GetParameter("densityFunc", generatorDefinition, "1"); TiXmlElement* sourceDefinition = GetElement("source", generatorDefinition); while (sourceDefinition) { @@ -965,7 +926,7 @@ void TRestGeant4Metadata::ReadGenerator() { // check if the generator is valid. if (GetNumberOfSources() == 0) { - RESTError << "No sources are added inside geneartor!" << RESTendl; + RESTError << "No sources are added inside generator!" << RESTendl; exit(1); } } @@ -977,61 +938,90 @@ void TRestGeant4Metadata::ReadParticleSource(TRestGeant4ParticleSource* source, TiXmlElement* sourceDefinition = definition; source->SetParticleName(GetParameter("particle", sourceDefinition)); - source->SetExcitationLevel(StringToDouble(GetParameter("excitedLevel", sourceDefinition))); - - Int_t charge = 0; - if (GetParameter("charge", sourceDefinition) == PARAMETER_NOT_FOUND_STR) - charge = 0; - else - charge = StringToInteger(GetParameter("charge", sourceDefinition)); - - source->SetParticleCharge(charge); + source->SetParticleCharge(StringToInteger(GetParameter("charge", sourceDefinition, "0"))); TString fullChain = GetParameter("fullChain", sourceDefinition); - - if (fullChain == "on" || fullChain == "ON" || fullChain == "On" || fullChain == "oN") { + SetFullChain(false); + if (fullChain.EqualTo("on", TString::ECaseCompare::kIgnoreCase)) { SetFullChain(true); - } else { - SetFullChain(false); } // Angular distribution parameters - TiXmlElement* angularDefinition = GetElement("angularDist", sourceDefinition); - source->SetAngularDistType(GetParameter("type", angularDefinition, "flux")); - if (source->GetAngularDistType() == "TH1D") { - source->SetAngularFilename(SearchFile(GetParameter("file", angularDefinition))); - source->SetAngularName(GetParameter("spctName", angularDefinition)); + TiXmlElement* angularDefinition = GetElement("angular", sourceDefinition); + if (angularDefinition == nullptr) { + angularDefinition = GetElement("angularDist", sourceDefinition); // old name + } + source->SetAngularDistributionType(GetParameter( + "type", angularDefinition, AngularDistributionTypesToString(AngularDistributionTypes::FLUX))); + if (StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) == + AngularDistributionTypes::TH1D) { + source->SetAngularDistributionFilename(SearchFile(GetParameter("file", angularDefinition))); + auto name = GetParameter("name", angularDefinition); + if (name == "NO_SUCH_PARA") { + name = GetParameter("spctName", angularDefinition); // old name + } + source->SetAngularDistributionNameInFile(name); } - if (GetNumberOfSources() == 0 && source->GetAngularDistType() == "backtoback") { + if (StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) == + AngularDistributionTypes::FORMULA) { + source->SetAngularDistributionFormula(GetParameter("name", angularDefinition)); + } + if (GetNumberOfSources() == 0 && + StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) == + AngularDistributionTypes::BACK_TO_BACK) { cout << "WARNING: First source cannot be backtoback. Setting it to isotropic" << endl; - source->SetAngularDistType("isotropic"); + source->SetAngularDistributionType( + AngularDistributionTypesToString(AngularDistributionTypes::ISOTROPIC)); } source->SetDirection(StringTo3DVector(GetParameter("direction", angularDefinition, "(1,0,0)"))); // Energy distribution parameters - TiXmlElement* energyDefinition = GetElement("energyDist", sourceDefinition); - source->SetEnergyDistType(GetParameter("type", energyDefinition, "mono")); - if (source->GetEnergyDistType() == "TH1D") { - source->SetSpectrumFilename(SearchFile(GetParameter("file", energyDefinition))); - source->SetSpectrumName(GetParameter("spctName", energyDefinition)); + TiXmlElement* energyDefinition = GetElement("energy", sourceDefinition); + if (energyDefinition == nullptr) { + energyDefinition = GetElement("energyDist", sourceDefinition); // old name + } + source->SetEnergyDistributionType(GetParameter( + "type", energyDefinition, EnergyDistributionTypesToString(EnergyDistributionTypes::MONO))); + if (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) == + EnergyDistributionTypes::TH1D) { + source->SetEnergyDistributionFilename(SearchFile(GetParameter("file", energyDefinition))); + auto name = GetParameter("name", energyDefinition); + if (name == "NO_SUCH_PARA") { + name = GetParameter("spctName", energyDefinition); // old name + } + source->SetEnergyDistributionNameInFile(name); } - source->SetEnergyRange(Get2DVectorParameterWithUnits("range", energyDefinition)); - if (source->GetEnergyDistType() == "mono") { - Double_t en; - en = GetDblParameterWithUnits("energy", energyDefinition, 0); - source->SetEnergyRange(TVector2(en, en)); - source->SetEnergy(en); + source->SetEnergyDistributionRange(Get2DVectorParameterWithUnits("range", energyDefinition, {0, 1E20})); + if (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) == + EnergyDistributionTypes::MONO) { + Double_t energy; + energy = GetDblParameterWithUnits("energy", energyDefinition, 0); + source->SetEnergyDistributionRange({energy, energy}); + source->SetEnergy(energy); + } + if (StringToEnergyDistributionTypes(source->GetEnergyDistributionType().Data()) == + EnergyDistributionTypes::FORMULA) { + source->SetEnergyDistributionFormula(GetParameter("name", energyDefinition)); + // We cannot use an energy range bigger than the range of the formula + const auto function = source->GetEnergyDistributionFunction(); + if (source->GetEnergyDistributionRangeMin() < function->GetXaxis()->GetXmin()) { + source->SetEnergyDistributionRange( + {function->GetXaxis()->GetXmin(), source->GetEnergyDistributionRangeMax()}); + } + if (source->GetEnergyDistributionRangeMax() > function->GetXaxis()->GetXmax()) { + source->SetEnergyDistributionRange( + {source->GetEnergyDistributionRangeMin(), function->GetXaxis()->GetXmax()}); + } } - // allow custom configuration from the class source->LoadConfigFromElement(sourceDefinition, fElementGlobal, fVariables); // AddParticleSource(source); } void TRestGeant4Metadata::RemoveParticleSources() { - for (auto c : fParticleSource) { - delete c; + for (auto source : fParticleSource) { + delete source; } fParticleSource.clear(); } @@ -1045,76 +1035,162 @@ void TRestGeant4Metadata::AddParticleSource(TRestGeant4ParticleSource* src) { /// /// This section allows to define which hits will be stored to disk. /// Different volumes in the geometry can be tagged for hit storage. +/// At least one volume must be tagged as sensitive /// /// \code -/// +/// /// -/// -/// +/// +/// /// \endcode /// -void TRestGeant4Metadata::ReadStorage() { - TiXmlElement* storageDefinition = GetElement("storage"); - fSensitiveVolume = GetFieldValue("sensitiveVolume", storageDefinition); - if (fSensitiveVolume == "Not defined") { - RESTWarning << "Sensitive volume not defined. Setting it to 'gas'!" << RESTendl; - fSensitiveVolume = "gas"; +void TRestGeant4Metadata::ReadDetector() { + RESTInfo << "TRestGeant4Metadata: Processing detector section" << RESTendl; + TiXmlElement* detectorDefinition = GetElement("detector"); + if (detectorDefinition == nullptr) { + RESTError << "Detector section () is not present in metadata!" << RESTendl; + exit(1); } - Double_t defaultStep = GetDblParameterWithUnits("maxStepSize", storageDefinition); - if (defaultStep < 0) defaultStep = 0; - RESTInfo << "Sensitive volume: " << fSensitiveVolume << RESTendl; + const bool activateAllVolumes = + StringToBool(GetParameter("activateAllVolumes", detectorDefinition, "false")); + RESTInfo << "TRestGeant4Metadata: Setting 'fActivateAllVolumes' to " << fActivateAllVolumes << RESTendl; + fActivateAllVolumes = activateAllVolumes; + + TiXmlElement* removeUnwantedTracksDefinition = GetElement("removeUnwantedTracks", detectorDefinition); + if (removeUnwantedTracksDefinition != nullptr) { + fRemoveUnwantedTracks = StringToBool(GetParameter("enabled", removeUnwantedTracksDefinition, "true")); + fRemoveUnwantedTracksKeepZeroEnergyTracks = + StringToBool(GetParameter("keepZeroEnergyTracks", removeUnwantedTracksDefinition, "false")); + RESTInfo << "TRestGeant4Metadata: Setting 'fRemoveUnwantedTracks' to " << fRemoveUnwantedTracks + << " with 'keepZeroEnergyTracks' option set to " << fRemoveUnwantedTracksKeepZeroEnergyTracks + << RESTendl; + } - fEnergyRangeStored = Get2DVectorParameterWithUnits("energyRange", storageDefinition); + Double_t defaultStep = GetDblParameterWithUnits("maxStepSize", detectorDefinition); + if (defaultStep < 0) { + defaultStep = 0; + } - auto gdml = new TRestGDMLParser(); - gdml->Load(GetGdmlFilename().Data()); + TiXmlElement* volumeDefinition = GetElement("volume", detectorDefinition); + while (volumeDefinition != nullptr) { + string name = GetFieldValue("name", volumeDefinition); + if (name.empty() || name == "Not defined") { + RESTError << "volume inside detector section defined without name" << RESTendl; + exit(1); + } + vector physicalVolumes; + if (!fGeant4GeometryInfo.IsValidGdmlName(name)) { + if (fGeant4GeometryInfo.IsValidLogicalVolume(name)) { + for (const auto& physical : fGeant4GeometryInfo.GetAllPhysicalVolumesFromLogical(name)) { + physicalVolumes.emplace_back( + fGeant4GeometryInfo.GetAlternativeNameFromGeant4PhysicalName(physical)); + } + } + // does not match a explicit (gdml) physical name or a logical name, perhaps its a regular exp + if (physicalVolumes.empty()) { + for (const auto& physical : + fGeant4GeometryInfo.GetAllPhysicalVolumesMatchingExpression(name)) { + physicalVolumes.emplace_back( + fGeant4GeometryInfo.GetAlternativeNameFromGeant4PhysicalName(physical)); + } + } + if (physicalVolumes.empty()) { + for (const auto& logical : fGeant4GeometryInfo.GetAllLogicalVolumesMatchingExpression(name)) { + for (const auto& physical : + fGeant4GeometryInfo.GetAllPhysicalVolumesFromLogical(logical)) { + physicalVolumes.emplace_back( + fGeant4GeometryInfo.GetAlternativeNameFromGeant4PhysicalName(physical)); + } + } + } + } else { + physicalVolumes.push_back(name); + } - fGeant4GeometryInfo.PopulateFromGdml(gdml->GetOutputGDMLFile()); + if (physicalVolumes.empty()) { + RESTError + << "volume '" << name + << "' is not defined in the geometry. Could not match to a physical, logical volume or regex" + << RESTendl; + exit(1); + } - const auto& physicalVolumes = fGeant4GeometryInfo.fGdmlNewPhysicalNames; - TiXmlElement* volumeDefinition = GetElement("activeVolume", storageDefinition); - while (volumeDefinition) { - Double_t chance = StringToDouble(GetFieldValue("chance", volumeDefinition)); - if (chance == -1) chance = 1; + bool isSensitive = false; + const string isSensitiveValue = GetFieldValue("sensitive", volumeDefinition); + if (isSensitiveValue != "Not defined") { + isSensitive = StringToBool(isSensitiveValue); + } + bool isKeepTracks = isSensitive; + const string isKeepTracksValue = GetFieldValue("keepTracks", volumeDefinition); + if (isKeepTracksValue != "Not defined") { + isKeepTracks = StringToBool(isKeepTracksValue); + } + + Double_t chance = StringToDouble(GetFieldValue("chance", volumeDefinition)); + if (chance == -1 || chance < 0) { + chance = 1.0; + } Double_t maxStep = GetDblParameterWithUnits("maxStepSize", volumeDefinition); - if (maxStep < 0) maxStep = defaultStep; + if (maxStep < 0) { + maxStep = defaultStep; + } - const TString& name = GetFieldValue("name", volumeDefinition); - // first we verify it's in the list of valid volumes - if (!fGeant4GeometryInfo.IsValidGdmlName(name)) { - bool isValidLogical = false; - for (size_t i = 0; i < fGeant4GeometryInfo.fGdmlNewPhysicalNames.size(); i++) { - if (fGeant4GeometryInfo.fGdmlLogicalNames[i] == name) { - isValidLogical = true; - const auto& gdmlName = fGeant4GeometryInfo.fGdmlNewPhysicalNames[i]; - RESTInfo << "Adding active volume from RML: '" << gdmlName << "' from logical volume: '" - << name << "' with chance: " << chance << RESTendl; - SetActiveVolume(gdmlName, chance, maxStep); - } + for (const auto& physical : physicalVolumes) { + RESTInfo << "Adding " << (isSensitive ? "sensitive" : "active") << " volume from RML: '" + << physical << (chance != 1 ? " with change: " + to_string(chance) : " ") << RESTendl; + SetActiveVolume(physical, chance, maxStep); + if (isSensitive) { + InsertSensitiveVolume(physical); } - if (!isValidLogical) { - RESTError << "TRestGeant4Metadata: Problem reading storage section." << RESTendl; - RESTError << " - The volume '" << name << "' was not found in the GDML geometry." - << RESTendl; - exit(1); + if (fRemoveUnwantedTracks && isKeepTracks) { + fRemoveUnwantedTracksVolumesToKeep.insert(physical); } - } else { - RESTInfo << "Adding active volume from RML: '" << name << "' with chance: " << chance << RESTendl; - SetActiveVolume(name, chance, maxStep); } volumeDefinition = GetNextElement(volumeDefinition); } + if (fSensitiveVolumes.empty()) { + cout << "No sensitive volumes defined in detector section. Adding 'gas' as sensitive volume" << endl; + InsertSensitiveVolume("gas"); + } + + fEnergyRangeStored = Get2DVectorParameterWithUnits("energyRange", detectorDefinition, fEnergyRangeStored); + // TODO: Place this as a generic method + if (fEnergyRangeStored.X() < 0) { + fEnergyRangeStored.Set(0, fEnergyRangeStored.Y()); + } + if (fEnergyRangeStored.Y() < 0) { + fEnergyRangeStored.Set(fEnergyRangeStored.X(), 0); + } + if (fEnergyRangeStored.Y() > fEnergyRangeStored.Y()) { + fEnergyRangeStored.Set(fEnergyRangeStored.Y(), fEnergyRangeStored.X()); + } + if (fEnergyRangeStored.Y() <= 0) { + RESTError << "Energy range is not valid: (" << fEnergyRangeStored.X() << ", " + << fEnergyRangeStored.Y() << ") keV" << RESTendl; + exit(1); + } + + auto gdml = new TRestGDMLParser(); + gdml->Load(GetGdmlFilename().Data()); + + fGeant4GeometryInfo.PopulateFromGdml(gdml->GetOutputGDMLFile()); + // If the user did not add explicitly any volume to the storage section we understand // the user wants to register all the volumes - if (GetNumberOfActiveVolumes() == 0) { - for (auto& name : physicalVolumes) { - SetActiveVolume(name, 1, defaultStep); - RESTInfo << "Automatically adding active volume: '" << name << "' with chance: " << 1 << RESTendl; + if (fActivateAllVolumes || GetNumberOfActiveVolumes() == 0) { + for (auto& name : fGeant4GeometryInfo.fGdmlNewPhysicalNames) { + if (!IsActiveVolume(name)) { + SetActiveVolume(name, 1, defaultStep); + RESTInfo << "Automatically adding active volume: '" << name << "' with chance: " << 1 + << RESTendl; + } } } + + fDetectorSectionInitialized = true; } /////////////////////////////////////////////// @@ -1124,382 +1200,63 @@ void TRestGeant4Metadata::ReadStorage() { void TRestGeant4Metadata::PrintMetadata() { TRestMetadata::PrintMetadata(); - RESTMetadata << "Geant 4 version : " << GetGeant4Version() << RESTendl; - RESTMetadata << "Random seed : " << GetSeed() << RESTendl; - RESTMetadata << "GDML geometry : " << GetGdmlReference() << RESTendl; - RESTMetadata << "GDML materials reference : " << GetMaterialsReference() << RESTendl; - RESTMetadata << "Sub-event time delay : " << GetSubEventTimeDelay() << " us" << RESTendl; + RESTMetadata << "Geant4 version: " << GetGeant4Version() << RESTendl; + RESTMetadata << "Random seed: " << GetSeed() << RESTendl; + RESTMetadata << "GDML geometry: " << GetGdmlReference() << RESTendl; + RESTMetadata << "GDML materials reference: " << GetMaterialsReference() << RESTendl; + RESTMetadata << "Sub-event time delay: " << GetSubEventTimeDelay() << " us" << RESTendl; Double_t mx = GetMagneticField().X(); Double_t my = GetMagneticField().Y(); Double_t mz = GetMagneticField().Z(); - RESTMetadata << "Magnetic field : ( " << mx << ", " << my << ", " << mz << ") T" << RESTendl; + RESTMetadata << "Magnetic field: (" << mx << ", " << my << ", " << mz << ") T" << RESTendl; if (fSaveAllEvents) RESTMetadata << "Save all events was enabled!" << RESTendl; - if (fRegisterEmptyTracks) + if (fRegisterEmptyTracks) { RESTMetadata << "Register empty tracks was enabled" << RESTendl; - else - RESTMetadata << "Register empty tracks was NOT enabled" << RESTendl; - RESTMetadata << " ++++++++++ Generator +++++++++++ " << RESTendl; - RESTMetadata << "Number of generated events : " << GetNumberOfEvents() << RESTendl; - RESTMetadata << "Generator type : " << fGenType << RESTendl; - RESTMetadata << "Generator shape : " << fGenShape; - if (fGenShape == "gdml") { - RESTMetadata << "::" << GetGDMLGeneratorVolume() << RESTendl; } else { - if (fGenShape == "box") { - RESTMetadata << ", (length, width, height): "; - } else if (fGenShape == "sphere") { - RESTMetadata << ", (radius, , ): "; - } else if (fGenShape == "wall") { - RESTMetadata << ", (length, width, ): "; - } else if (fGenShape == "circle") { - RESTMetadata << ", (radius, , ): "; - } else if (fGenShape == "cylinder") { - RESTMetadata << ", (radius, length, ): "; - } - - if (fGenShape != "point") { - TVector3 s = GetGeneratorSize(); - RESTMetadata << s.X() << ", " << s.Y() << ", " << s.Z() << RESTendl; - } else { - RESTMetadata << RESTendl; - } + RESTMetadata << "Register empty tracks was NOT enabled" << RESTendl; } - TVector3 a = GetGeneratorPosition(); - RESTMetadata << "Generator center : (" << a.X() << "," << a.Y() << "," << a.Z() << ") mm" << RESTendl; - TVector3 b = GetGeneratorRotationAxis(); - RESTMetadata << "Generator rotation : (" << b.X() << "," << b.Y() << "," << b.Z() - << "), angle: " << GetGeneratorRotationDegree() << " rads" << RESTendl; + RESTMetadata << " ++++++++++ Generator +++++++++++ " << RESTendl; + RESTMetadata << "Number of generated events: " << GetNumberOfEvents() << RESTendl; + fGeant4PrimaryGeneratorInfo.Print(); for (int i = 0; i < GetNumberOfSources(); i++) GetParticleSource(i)->PrintParticleSource(); - RESTMetadata << " ++++++++++ Storage volumes +++++++++++ " << RESTendl; - RESTMetadata << "Energy range : Emin = " << GetMinimumEnergyStored() - << ", Emax : " << GetMaximumEnergyStored() << RESTendl; - RESTMetadata << "Sensitive volume : " << GetSensitiveVolume() << RESTendl; - RESTMetadata << "Active volumes : " << GetNumberOfActiveVolumes() << RESTendl; - RESTMetadata << "---------------------------------------" << RESTendl; - for (int n = 0; n < GetNumberOfActiveVolumes(); n++) { - RESTMetadata << "Name : " << GetActiveVolumeName(n) - << ", ID : " << GetActiveVolumeID(GetActiveVolumeName(n)) - << ", maxStep : " << GetMaxStepSize(GetActiveVolumeName(n)) << "mm " - << ", chance : " << GetStorageChance(GetActiveVolumeName(n)) << RESTendl; + if (fDetectorSectionInitialized) { + RESTMetadata << " ++++++++++ Detector +++++++++++ " << RESTendl; + RESTMetadata << "Energy range (keV): (" << GetMinimumEnergyStored() << ", " + << GetMaximumEnergyStored() << ")" << RESTendl; + RESTMetadata << "Number of sensitive volumes: " << GetNumberOfSensitiveVolumes() << RESTendl; + for (const auto& sensitiveVolume : fSensitiveVolumes) { + RESTMetadata << "Sensitive volume: " << sensitiveVolume << RESTendl; + } + RESTMetadata << "Number of active volumes: " << GetNumberOfActiveVolumes() << RESTendl; + for (int n = 0; n < GetNumberOfActiveVolumes(); n++) { + RESTMetadata << "Name: " << GetActiveVolumeName(n) + << ", ID: " << GetActiveVolumeID(GetActiveVolumeName(n)) + << ", maxStep: " << GetMaxStepSize(GetActiveVolumeName(n)) << "mm " + << ", chance: " << GetStorageChance(GetActiveVolumeName(n)) << RESTendl; + } } for (int n = 0; n < GetNumberOfBiasingVolumes(); n++) { GetBiasingVolume(n).PrintBiasingVolume(); } + + if (GetRemoveUnwantedTracks()) { + RESTMetadata << "removeUnwantedTracks is enabled " + << (fRemoveUnwantedTracksKeepZeroEnergyTracks ? "with" : "without") + << " keeping zero energy tracks" << RESTendl; + + /* + RESTMetadata << "keep volumes for removeUnwantedTracks:" << RESTendl; + for (const auto& volume : fRemoveUnwantedTracksVolumesToKeep) { + RESTMetadata << "\t" << volume << RESTendl; + } + */ + } + RESTMetadata << "+++++" << RESTendl; } -///////////////////////////////////////////////// -///// \brief Reads an input file produced by Decay0. -///// -///// The input file should contain the description of several -///// pre-generated events, providing the names (or ids) of -///// particles to be produced, their energy, and momentum. -///// The particles and their properties are stored in a -///// TRestG4ParticleCollection which will be randomly accessed -///// by the restG4 package. -///// -///// \param fName The Decay0 filename located at -///// REST_PATH/data/generator/ -///// -// void TRestGeant4Metadata::ReadEventDataFile(TString fName) { -// string fullPathName = SearchFile((string)fName); -// if (fullPathName == "") { -// ferr << "File not found : " << fName << endl; -// ferr << "Decay0 generator file could not be found!!" << endl; -// exit(1); -// } -// -// debug << "TRestGeant4Metadata::ReadGeneratorFile" << endl; -// debug << "Full path generator file : " << fullPathName << endl; -// -// if (!ReadOldDecay0File(fullPathName)) ReadNewDecay0File(fullPathName); -// } -// -///////////////////////////////////////////////// -///// \brief Reads particle information using the file format from newer Decay0 versions. -///// -///// This is an auxiliar method used in TRestGeant4Metadata::ReadEventDataFile that will read the Decay0 -/// files -///// produced with the newer Decay0 versions. -///// -// Int_t TRestGeant4Metadata::ReadNewDecay0File(TString fileName) { -// ifstream infile; -// infile.open(fileName); -// if (!infile.is_open()) { -// printf("Error when opening file %s\n", fileName.Data()); -// return 0; -// } -// -// int generatorEvents = 0; -// string s; -// // First lines are discarded. -// for (int i = 0; i < 24; i++) { -// getline(infile, s); -// int pos = s.find("@nevents="); -// if (pos != -1) generatorEvents = StringToInteger(s.substr(10, s.length() - 10)); -// } -// -// if (generatorEvents == 0) { -// ferr << "TRestGeant4Metadata::ReadNewDecay0File. Problem reading generator file" << endl; -// exit(1); -// } -// -// TRestGeant4ParticleSource* src = TRestGeant4ParticleSource::instantiate(); -// src->SetGenFilename(fileName); -// // src->SetAngularDistType("flux"); -// // src->SetEnergyDistType("mono"); -// -// TRestGeant4Particle particle; -// -// debug << "Reading generator file : " << fileName << endl; -// debug << "Total number of events : " << generatorEvents << endl; -// for (int n = 0; n < generatorEvents && !infile.eof(); n++) { -// int pos = -1; -// while (!infile.eof() && pos == -1) { -// getline(infile, s); -// pos = s.find("@event_start"); -// } -// -// // Time - nuclide is skipped -// getline(infile, s); -// -// Int_t nParticles; -// infile >> nParticles; -// debug << "Number of particles : " << nParticles << endl; -// -// // cout << evID <<" "<< time <<" "<< nParticles <<" "<< endl; -// for (int i = 0; i < nParticles && !infile.eof(); i++) { -// Int_t pID; -// Double_t momx, momy, momz, mass; -// Double_t energy = -1, momentum2; -// TString pName; -// -// Double_t time; -// infile >> pID >> time >> momx >> momy >> momz >> pName; -// -// debug << " ---- New particle found --- " << endl; -// debug << " Particle name : " << pName << endl; -// debug << " - pId : " << pID << endl; -// debug << " - Relative time : " << time << endl; -// debug << " - px: " << momx << " py: " << momy << " pz: " << momz << " " << endl; -// -// if (pID == 3) { -// momentum2 = (momx * momx) + (momy * momy) + (momz * momz); -// mass = 0.511; -// -// energy = TMath::Sqrt(momentum2 + mass * mass) - mass; -// particle.SetParticleName("e-"); -// particle.SetParticleCharge(-1); -// particle.SetExcitationLevel(0); -// -// } else if (pID == 1) { -// momentum2 = (momx * momx) + (momy * momy) + (momz * momz); -// -// energy = TMath::Sqrt(momentum2); -// particle.SetParticleName("gamma"); -// particle.SetParticleCharge(0); -// particle.SetExcitationLevel(0); -// } else { -// cout << "Particle id " << pID << " not recognized" << endl; -// } -// -// TVector3 momDirection(momx, momy, momz); -// momDirection = momDirection.Unit(); -// -// particle.SetEnergy(1000. * energy); -// particle.SetDirection(momDirection); -// -// src->AddParticle(particle); -// } -// src->FlushParticlesTemplate(); -// } -// -// AddParticleSource(src); -// -// return 1; -// } -// -///////////////////////////////////////////////// -///// \brief Reads particle information using the file format from older Decay0 versions. -///// -///// This is an auxiliar method used in TRestGeant4Metadata::ReadEventDataFile that will read the Decay0 -/// files -///// produced with the newer Decay0 versions. -///// -// Int_t TRestGeant4Metadata::ReadOldDecay0File(TString fileName) { -// ifstream infile; -// infile.open(fileName); -// if (!infile.is_open()) { -// printf("Error when opening file %s\n", fileName.Data()); -// return 0; -// } -// -// string s; -// // First lines are discarded. -// int headerFound = 0; -// for (int i = 0; i < 30; i++) { -// getline(infile, s); -// if (s.find("#!bxdecay0 1.0.0") != -1) return 0; -// if (s.find("First event and full number of events:") != -1) { -// headerFound = 1; -// break; -// } -// } -// if (!headerFound) { -// ferr -// << "TRestGeant4Metadata::ReadOldDecay0File. Problem reading generator file: no \"First event and -// " -// "full number of events:\" header.\n"; -// abort(); -// } -// int tmpInt; -// int fGeneratorEvents; -// infile >> tmpInt >> fGeneratorEvents; -// -// // cout << "i : " << tmpInt << " fN : " << fGeneratorEvents << endl; -// TRestGeant4ParticleSource* src = TRestGeant4ParticleSource::instantiate(); -// src->SetGenFilename(fileName); -// // src->SetAngularDistType("flux"); -// // src->SetEnergyDistType("mono"); -// -// TRestGeant4Particle particle; -// string type = (string)GetGeneratorType(); -// -// cout << "Reading generator file : " << fileName << endl; -// cout << "Total number of events : " << fGeneratorEvents << endl; -// for (int n = 0; n < fGeneratorEvents && !infile.eof(); n++) { -// Int_t nParticles; -// Int_t evID; -// Double_t time; -// -// infile >> evID >> time >> nParticles; -// -// // cout << evID <<" "<< time <<" "<< nParticles <<" "<< endl; -// for (int i = 0; i < nParticles && !infile.eof(); i++) { -// Int_t pID; -// Double_t momx, momy, momz, mass; -// Double_t energy = -1, momentum2; -// Double_t x, y, z; -// -// infile >> pID >> momx >> momy >> momz >> time; -// if (type == "file") infile >> x >> y >> z; -// -// // cout << momx << " " << momy << " " << momz << " " << endl; -// -// bool ise = 2 <= pID && pID <= 3, ismu = 5 <= pID && pID <= 6, isp = pID == 14, isg = pID == 1; -// if (ise || ismu || isp || isg) { -// momentum2 = (momx * momx) + (momy * momy) + (momz * momz); -// if (ise) { -// mass = 0.511; -// particle.SetParticleName(pID == 2 ? "e+" : "e-"); -// particle.SetParticleCharge(pID == 2 ? 1 : -1); -// } else if (ismu) { -// mass = 105.7; -// particle.SetParticleName(pID == 5 ? "mu+" : "mu-"); -// particle.SetParticleCharge(pID == 5 ? 1 : -1); -// } else if (isp) { -// mass = 938.3; -// particle.SetParticleName("proton"); -// particle.SetParticleCharge(1); -// } else { -// mass = 0; -// particle.SetParticleName("gamma"); -// particle.SetParticleCharge(0); -// } -// -// energy = TMath::Sqrt(momentum2 + mass * mass) - mass; -// particle.SetExcitationLevel(0); -// } else { -// cout << "Particle id " << pID << " not recognized" << endl; -// } -// -// TVector3 momDirection(momx, momy, momz); -// momDirection = momDirection.Unit(); -// -// particle.SetEnergy(1000. * energy); -// particle.SetDirection(momDirection); -// particle.SetOrigin(TVector3(x, y, z)); -// -// src->AddParticle(particle); -// } -// -// src->FlushParticlesTemplate(); -// } -// -// AddParticleSource(src); -// -// return 1; -// } -// -///////////////////////////////////////////////// -///// \brief It reads the SetParticleName(GetFieldValue("particle", sourceDefinition)); -// -// source->SetExcitationLevel(StringToDouble(GetFieldValue("excitedLevel", sourceDefinition))); -// -// Int_t charge = 0; -// if (GetFieldValue("charge", sourceDefinition) == "Not defined") -// charge = 0; -// else -// charge = StringToInteger(GetFieldValue("charge", sourceDefinition)); -// -// source->SetParticleCharge(charge); -// -// TString fullChain = GetFieldValue("fullChain", sourceDefinition); -// -// if (fullChain == "on" || fullChain == "ON" || fullChain == "On" || fullChain == "oN") { -// SetFullChain(true); -// } else { -// SetFullChain(false); -// } -// -// // Angular distribution parameters -// TiXmlElement* angularDefinition = GetElement("angularDist", sourceDefinition); -// -// source->SetAngularDistType(GetFieldValue("type", angularDefinition)); -// -// if (source->GetAngularDistType() == "TH1D") { -// source->SetAngularFilename(SearchFile(GetFieldValue("file", angularDefinition))); -// source->SetAngularName(GetFieldValue("spctName", angularDefinition)); -// } -// -// if (GetNumberOfSources() == 0 && source->GetAngularDistType() == "backtoback") { -// cout << "WARNING: First source cannot be backtoback. Setting it to isotropic" << endl; -// source->SetAngularDistType("isotropic"); -// } -// -// source->SetDirection(StringTo3DVector(GetFieldValue("direction", angularDefinition))); -// -// // Energy distribution parameters -// TiXmlElement* energyDefinition = GetElement("energyDist", sourceDefinition); -// -// source->SetEnergyDistType(GetFieldValue("type", energyDefinition)); -// -// if (source->GetEnergyDistType() == "TH1D") { -// source->SetSpectrumFilename(SearchFile(GetFieldValue("file", energyDefinition))); -// source->SetSpectrumName(GetFieldValue("spctName", energyDefinition)); -// } -// -// source->SetEnergyRange(Get2DVectorParameterWithUnits("range", energyDefinition)); -// -// if (source->GetEnergyDistType() == "mono") { -// Double_t en; -// en = GetDblParameterWithUnits("energy", energyDefinition); -// source->SetEnergyRange(TVector2(en, en)); -// source->SetEnergy(en); -// } -// -// AddParticleSource(source); -// } - /////////////////////////////////////////////// /// \brief Returns the id of an active volume giving as parameter its name. Int_t TRestGeant4Metadata::GetActiveVolumeID(TString name) { @@ -1536,15 +1293,16 @@ void TRestGeant4Metadata::SetActiveVolume(const TString& name, Double_t chance, fActiveVolumes.push_back(name); fChance.push_back(chance); fMaxStepSize.push_back(maxStep); + fActiveVolumesSet.insert(name.Data()); } /////////////////////////////////////////////// /// \brief Returns true if the volume named *volName* has been registered for /// data storage. /// -Bool_t TRestGeant4Metadata::isVolumeStored(TString volName) { +Bool_t TRestGeant4Metadata::isVolumeStored(const TString& volume) const { for (int n = 0; n < GetNumberOfActiveVolumes(); n++) - if (GetActiveVolumeName(n) == volName) return true; + if (GetActiveVolumeName(n) == volume) return true; return false; } @@ -1552,12 +1310,12 @@ Bool_t TRestGeant4Metadata::isVolumeStored(TString volName) { /////////////////////////////////////////////// /// \brief Returns the probability of an active volume being stored /// -Double_t TRestGeant4Metadata::GetStorageChance(TString vol) { +Double_t TRestGeant4Metadata::GetStorageChance(TString volume) { Int_t id; for (id = 0; id < (Int_t)fActiveVolumes.size(); id++) { - if (fActiveVolumes[id] == vol) return fChance[id]; + if (fActiveVolumes[id] == volume) return fChance[id]; } - RESTWarning << "TRestGeant4Metadata::GetStorageChance. Volume " << vol << " not found" << RESTendl; + RESTWarning << "TRestGeant4Metadata::GetStorageChance. Volume " << volume << " not found" << RESTendl; return 0; } @@ -1565,11 +1323,13 @@ Double_t TRestGeant4Metadata::GetStorageChance(TString vol) { /////////////////////////////////////////////// /// \brief Returns the maximum step at a particular active volume /// -Double_t TRestGeant4Metadata::GetMaxStepSize(TString vol) { - for (Int_t id = 0; id < (Int_t)fActiveVolumes.size(); id++) { - if (fActiveVolumes[id] == vol) return fMaxStepSize[id]; +Double_t TRestGeant4Metadata::GetMaxStepSize(const TString& volume) { + for (Int_t i = 0; i < (Int_t)fActiveVolumes.size(); i++) { + if (volume.EqualTo(fActiveVolumes[i], TString::kIgnoreCase)) { + return fMaxStepSize[i]; + } } - RESTWarning << "TRestGeant4Metadata::GetMaxStepSize. Volume " << vol << " not found" << RESTendl; + RESTWarning << "TRestGeant4Metadata::GetMaxStepSize. Volume " << volume << " not found" << RESTendl; return 0; } diff --git a/src/TRestGeant4NeutronTaggingProcess.cxx b/src/TRestGeant4NeutronTaggingProcess.cxx index 33e518a..d27071f 100644 --- a/src/TRestGeant4NeutronTaggingProcess.cxx +++ b/src/TRestGeant4NeutronTaggingProcess.cxx @@ -352,7 +352,7 @@ TRestEvent* TRestGeant4NeutronTaggingProcess::ProcessEvent(TRestEvent* inputEven } } fNeutronsCapturedIsCaptureVolume.push_back(isCaptureVolume); - fNeutronsCapturedProductionE.push_back(track.GetKineticEnergy()); + fNeutronsCapturedProductionE.push_back(track.GetInitialKineticEnergy()); // get energy deposited by neutron that undergoes capture and children double neutronsCapturedEDepByNeutron = 0; @@ -445,7 +445,7 @@ TRestEvent* TRestGeant4NeutronTaggingProcess::ProcessEvent(TRestEvent* inputEven } } fGammasNeutronCaptureIsCaptureVolume.push_back(isCaptureVolume); - fGammasNeutronCaptureProductionE.push_back(track.GetKineticEnergy()); + fGammasNeutronCaptureProductionE.push_back(track.GetInitialKineticEnergy()); // cout << "gamma capture" << endl; @@ -505,7 +505,7 @@ TRestEvent* TRestGeant4NeutronTaggingProcess::ProcessEvent(TRestEvent* inputEven fSecondaryNeutronsShieldingIsCapturedInCaptureVolume.push_back(0); } - fSecondaryNeutronsShieldingProductionE.push_back(track.GetKineticEnergy()); + fSecondaryNeutronsShieldingProductionE.push_back(track.GetInitialKineticEnergy()); fSecondaryNeutronsShieldingExitE.push_back(hits.GetKineticEnergy(j)); } } diff --git a/src/TRestGeant4ParticleSource.cxx b/src/TRestGeant4ParticleSource.cxx index 16b6c45..1dc80a4 100644 --- a/src/TRestGeant4ParticleSource.cxx +++ b/src/TRestGeant4ParticleSource.cxx @@ -3,7 +3,7 @@ ///_______________________________________________________________________________ /// /// -/// RESTSoft : Software for Rare Event Searches with TPCs +/// RESTSoft: Software for Rare Event Searches with TPCs /// /// TRestGeant4ParticleSource.cxx /// @@ -17,32 +17,20 @@ #include "TRestGeant4ParticleSource.h" +#include +#include +#include + #include "TRestGeant4Metadata.h" -#include "TRestReflector.h" -#include "TRestStringHelper.h" -#include "TRestStringOutput.h" using namespace std; - -// REST_Verbose_Level fLevel = REST_Essential; //! -//// TRestLeveledOutput(REST_Verbose_Level& vref, string _color = -//// COLOR_RESET, string BorderOrHeader = "", REST_Display_Format style = -//// kBorderedLeft) -// TRestLeveledOutput metadata = -// TRestLeveledOutput(fLevel, COLOR_BOLDGREEN, "||", -// (REST_Display_Format)kBorderedMiddle); //! +using namespace TRestGeant4PrimaryGeneratorTypes; ClassImp(TRestGeant4ParticleSource); -TRestGeant4ParticleSource::TRestGeant4ParticleSource() { - // TRestGeant4ParticleSource default constructor - fAngularDistType = "flux"; - fEnergyDistType = "mono"; -} +TRestGeant4ParticleSource::TRestGeant4ParticleSource() = default; -TRestGeant4ParticleSource::~TRestGeant4ParticleSource() { - // TRestGeant4ParticleSource destructor -} +TRestGeant4ParticleSource::~TRestGeant4ParticleSource() = default; void TRestGeant4ParticleSource::PrintParticleSource() { RESTMetadata << "---------------------------------------" << RESTendl; @@ -51,38 +39,46 @@ void TRestGeant4ParticleSource::PrintParticleSource() { RESTMetadata << "Generator file: " << GetGenFilename() << RESTendl; RESTMetadata << "Stored templates: " << fParticlesTemplate.size() << RESTendl; RESTMetadata << "Particles: "; - for (auto p : fParticles) RESTMetadata << p.GetParticleName() << ", "; + for (const auto& particle : fParticles) RESTMetadata << particle.GetParticleName() << ", "; RESTMetadata << RESTendl; } else { - RESTMetadata << "Charge : " << GetParticleCharge() << RESTendl; - RESTMetadata << "Angular distribution type : " << GetAngularDistType() << RESTendl; - if (GetAngularDistType() == "TH1D") { - RESTMetadata << "Angular distribution filename : " - << TRestTools::GetPureFileName((string)GetAngularFilename()) << RESTendl; - RESTMetadata << "Angular histogram name : " << GetAngularName() << RESTendl; + if (GetParticleCharge() != 0) { + RESTMetadata << "Particle charge: " << GetParticleCharge() << RESTendl; + } + RESTMetadata << "Angular distribution type: " << GetAngularDistributionType() << RESTendl; + if (StringToAngularDistributionTypes(GetAngularDistributionType().Data()) == + AngularDistributionTypes::TH1D) { + RESTMetadata << "Angular distribution filename: " + << TRestTools::GetPureFileName((string)GetAngularDistributionFilename()) << RESTendl; + RESTMetadata << "Angular distribution histogram name: " << GetAngularDistributionNameInFile() + << RESTendl; + } + RESTMetadata << "Angular distribution direction: (" << GetDirection().X() << "," << GetDirection().Y() + << "," << GetDirection().Z() << ")" << RESTendl; + RESTMetadata << "Energy distribution type: " << GetEnergyDistributionType() << RESTendl; + if (StringToEnergyDistributionTypes(GetEnergyDistributionType().Data()) == + EnergyDistributionTypes::TH1D) { + RESTMetadata << "Energy distribution filename: " + << TRestTools::GetPureFileName((string)GetEnergyDistributionFilename()) << RESTendl; + RESTMetadata << "Energy distribution histogram name: " << GetEnergyDistributionNameInFile() + << RESTendl; + } + if (GetEnergyDistributionRangeMin() == GetEnergyDistributionRangeMax()) { + RESTMetadata << "Energy distribution energy: " << GetEnergy() << " keV" << RESTendl; + } else { + RESTMetadata << "Energy distribution range (keV): (" << GetEnergyDistributionRangeMin() << ", " + << GetEnergyDistributionRangeMax() << ")" << RESTendl; } - RESTMetadata << "Direction : (" << GetDirection().X() << "," << GetDirection().Y() << "," - << GetDirection().Z() << ")" << RESTendl; - RESTMetadata << "Energy distribution : " << GetEnergyDistType() << RESTendl; - if (GetEnergyDistType() == "TH1D") { - RESTMetadata << "Energy distribution filename : " - << TRestTools::GetPureFileName((string)GetSpectrumFilename()) << RESTendl; - RESTMetadata << "Energy histogram name : " << GetSpectrumName() << RESTendl; - } else if (GetEnergyRange().X() == GetEnergyRange().Y()) - RESTMetadata << "Energy : " << GetEnergy() << " keV" << RESTendl; - else - RESTMetadata << "Energy range : (" << GetEnergyRange().X() << "," << GetEnergyRange().Y() << ") keV" - << RESTendl; } } TRestGeant4ParticleSource* TRestGeant4ParticleSource::instantiate(std::string model) { - if (model == "" || model == "geant4" || model.find(".dat") != -1) { + if (model.empty() || model == "geant4" || model.find(".dat") != -1) { // use default generator return new TRestGeant4ParticleSource(); } else { // use specific generator - // naming convension: TRestGeant4ParticleSourceXXX + // naming convention: TRestGeant4ParticleSourceXXX // currently supported generator: decay0 // in future we may add wrapper of generators: cry(for muon), pythia(for HEP), etc. model[0] = *REST_StringHelper::ToUpper(std::string(&model[0], 1)).c_str(); @@ -103,7 +99,7 @@ void TRestGeant4ParticleSource::InitFromConfigFile() { if (((string)modelUse).find(".dat") != -1) { string fullPathName = SearchFile((string)modelUse); if (fullPathName == "") { - RESTError << "File not found : " << modelUse << RESTendl; + RESTError << "File not found: " << modelUse << RESTendl; RESTError << "Decay0 generator file could not be found!!" << RESTendl; exit(1); } @@ -120,17 +116,14 @@ void TRestGeant4ParticleSource::InitFromConfigFile() { // base class's generator action: randomize the particle's energy/direction with distribution file void TRestGeant4ParticleSource::Update() { - if (fParticlesTemplate.size() > 0) { + if (!fParticlesTemplate.empty()) { // we use particle template to generate particles - Int_t rndCollection = (Int_t)(fRndmMethod() * fParticlesTemplate.size()); + Int_t rndCollection = (Int_t)(fRandomMethod() * fParticlesTemplate.size()); Int_t pCollectionID = rndCollection % fParticlesTemplate.size(); fParticles = fParticlesTemplate[pCollectionID]; } else { TRestGeant4Particle p(*this); - // Future: implement particle generation for toy simulation - // - // - + // TODO: implement particle generation for toy simulation fParticles = {p}; } } @@ -182,15 +175,10 @@ bool TRestGeant4ParticleSource::ReadNewDecay0File(TString fileName) { exit(1); } - // TRestGeant4ParticleSource* src = TRestGeant4ParticleSource::instantiate(); - // this->SetGenFilename(fileName); - // this->SetAngularDistType("flux"); - // this->SetEnergyDistType("mono"); - TRestGeant4Particle particle; - RESTDebug << "Reading generator file : " << fileName << RESTendl; - RESTDebug << "Total number of events : " << generatorEvents << RESTendl; + RESTDebug << "Reading generator file: " << fileName << RESTendl; + RESTDebug << "Total number of events: " << generatorEvents << RESTendl; for (int n = 0; n < generatorEvents && !infile.eof(); n++) { int pos = -1; while (!infile.eof() && pos == -1) { @@ -203,7 +191,7 @@ bool TRestGeant4ParticleSource::ReadNewDecay0File(TString fileName) { Int_t nParticles; infile >> nParticles; - RESTDebug << "Number of particles : " << nParticles << RESTendl; + RESTDebug << "Number of particles: " << nParticles << RESTendl; // cout << evID <<" "<< time <<" "<< nParticles <<" "<< std::endl; for (int i = 0; i < nParticles && !infile.eof(); i++) { @@ -216,9 +204,9 @@ bool TRestGeant4ParticleSource::ReadNewDecay0File(TString fileName) { infile >> pID >> time >> momx >> momy >> momz >> pName; RESTDebug << " ---- New particle found --- " << RESTendl; - RESTDebug << " Particle name : " << pName << RESTendl; - RESTDebug << " - pId : " << pID << RESTendl; - RESTDebug << " - Relative time : " << time << RESTendl; + RESTDebug << " Particle name: " << pName << RESTendl; + RESTDebug << " - pId: " << pID << RESTendl; + RESTDebug << " - Relative time: " << time << RESTendl; RESTDebug << " - px: " << momx << " py: " << momy << " pz: " << momz << " " << RESTendl; if (pID == 3) { @@ -281,24 +269,25 @@ bool TRestGeant4ParticleSource::ReadOldDecay0File(TString fileName) { } } if (!headerFound) { - RESTError << "TRestG4Metadata::ReadOldDecay0File. Problem reading generator file: no \"First event and " - "full number of events:\" header.\n"; + RESTError + << "TRestG4Metadata::ReadOldDecay0File. Problem reading generator file: no \"First event and " + "full number of events:\" header.\n"; abort(); } int tmpInt; int fGeneratorEvents; infile >> tmpInt >> fGeneratorEvents; - // cout << "i : " << tmpInt << " fN : " << fGeneratorEvents << std::endl; + // cout << "i: " << tmpInt << " fN: " << fGeneratorEvents << std::endl; // TRestGeant4ParticleSource* src = TRestGeant4ParticleSource::instantiate(); // this->SetGenFilename(fileName); - // this->SetAngularDistType("flux"); - // this->SetEnergyDistType("mono"); + // this->SetAngularDistributionType("flux"); + // this->SetEnergyDistributionType("mono"); TRestGeant4Particle particle; - cout << "Reading generator file : " << fileName << std::endl; - cout << "Total number of events : " << fGeneratorEvents << std::endl; + cout << "Reading generator file: " << fileName << std::endl; + cout << "Total number of events: " << fGeneratorEvents << std::endl; for (int n = 0; n < fGeneratorEvents && !infile.eof(); n++) { Int_t nParticles; Int_t evID; diff --git a/src/TRestGeant4PhysicsInfo.cxx b/src/TRestGeant4PhysicsInfo.cxx index 4f7e4da..0f833ae 100644 --- a/src/TRestGeant4PhysicsInfo.cxx +++ b/src/TRestGeant4PhysicsInfo.cxx @@ -7,29 +7,45 @@ using namespace std; ClassImp(TRestGeant4PhysicsInfo); -void TRestGeant4PhysicsInfo::PrintParticles() const { - vector ids = {}; - for (const auto& [id, _] : fParticleNamesMap) { - ids.push_back(id); +set TRestGeant4PhysicsInfo::GetAllParticles() const { + set particles = {}; + for (const auto& [_, name] : fParticleNamesMap) { + particles.insert(name); } - sort(ids.begin(), ids.end()); + return particles; +} - cout << "Particles:" << endl; - for (const auto& id : ids) { - cout << "\t" << id << " - " << GetParticleName(id) << endl; +std::set TRestGeant4PhysicsInfo::GetAllProcesses() const { + set processes = {}; + for (const auto& [_, name] : fProcessNamesMap) { + processes.insert(name); } + return processes; } -void TRestGeant4PhysicsInfo::PrintProcesses() const { - vector ids = {}; - for (const auto& [id, _] : fProcessNamesMap) { - ids.push_back(id); +std::set TRestGeant4PhysicsInfo::GetAllProcessTypes() const { + set types = {}; + for (const auto& [_, type] : fProcessTypesMap) { + types.insert(type); + } + return types; +} + +void TRestGeant4PhysicsInfo::PrintParticles() const { + const auto particleNames = GetAllParticles(); + cout << "Particles:" << endl; + for (const auto& name : particleNames) { + const auto id = GetParticleID(name); + cout << "\t" << name << " - " << id << endl; } - sort(ids.begin(), ids.end()); +} +void TRestGeant4PhysicsInfo::PrintProcesses() const { + const auto processNames = GetAllProcesses(); cout << "Processes:" << endl; - for (const auto& id : ids) { - cout << "\t" << id << " - " << GetProcessName(id) << endl; + for (const auto& name : processNames) { + const auto id = GetProcessID(name); + cout << "\t" << name << " - " << id << endl; } } @@ -38,12 +54,23 @@ void TRestGeant4PhysicsInfo::Print() const { PrintProcesses(); } -void TRestGeant4PhysicsInfo::InsertProcessName(Int_t id, const TString& processName) { +void TRestGeant4PhysicsInfo::InsertProcessName(Int_t id, const TString& processName, + const TString& processType) { + if (fProcessNamesMap.count(id) > 0) { + return; + } + std::lock_guard lock(fMutex); fProcessNamesMap[id] = processName; fProcessNamesReverseMap[processName] = id; + + fProcessTypesMap[processName] = processType; } void TRestGeant4PhysicsInfo::InsertParticleName(Int_t id, const TString& particleName) { + if (fParticleNamesMap.count(id) > 0) { + return; + } + std::lock_guard lock(fMutex); fParticleNamesMap[id] = particleName; fParticleNamesReverseMap[particleName] = id; } @@ -71,3 +98,7 @@ TString TRestGeant4PhysicsInfo::GetParticleName(Int_t id) const { Int_t TRestGeant4PhysicsInfo::GetParticleID(const TString& processName) const { return GetOrDefaultMapValueFromKey(&fParticleNamesReverseMap, processName); } + +TString TRestGeant4PhysicsInfo::GetProcessType(const TString& processName) const { + return GetOrDefaultMapValueFromKey(&fProcessTypesMap, processName); +} diff --git a/src/TRestGeant4PrimaryGeneratorInfo.cxx b/src/TRestGeant4PrimaryGeneratorInfo.cxx new file mode 100644 index 0000000..1433051 --- /dev/null +++ b/src/TRestGeant4PrimaryGeneratorInfo.cxx @@ -0,0 +1,356 @@ +// +// Created by lobis on 7/14/2022. +// + +#include "TRestGeant4PrimaryGeneratorInfo.h" + +#include +#include +#include +#include + +#include + +using namespace std; +using namespace TRestGeant4PrimaryGeneratorTypes; + +string TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypesToString(const SpatialGeneratorTypes& type) { + switch (type) { + case SpatialGeneratorTypes::CUSTOM: + return "Custom"; + case SpatialGeneratorTypes::VOLUME: + return "Volume"; + case SpatialGeneratorTypes::SURFACE: + return "Surface"; + case SpatialGeneratorTypes::POINT: + return "Point"; + } + cout << "TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypesToString - Error - Unknown " + "SpatialGeneratorTypes" + << endl; + exit(1); +} + +SpatialGeneratorTypes TRestGeant4PrimaryGeneratorTypes::StringToSpatialGeneratorTypes(const string& type) { + if (TString(type).EqualTo(SpatialGeneratorTypesToString(SpatialGeneratorTypes::CUSTOM), + TString::ECaseCompare::kIgnoreCase)) { + return SpatialGeneratorTypes::CUSTOM; + } else if (TString(type).EqualTo(SpatialGeneratorTypesToString(SpatialGeneratorTypes::VOLUME), + TString::ECaseCompare::kIgnoreCase)) { + return SpatialGeneratorTypes::VOLUME; + } else if (TString(type).EqualTo(SpatialGeneratorTypesToString(SpatialGeneratorTypes::SURFACE), + TString::ECaseCompare::kIgnoreCase)) { + return SpatialGeneratorTypes::SURFACE; + } else if (TString(type).EqualTo(SpatialGeneratorTypesToString(SpatialGeneratorTypes::POINT), + TString::ECaseCompare::kIgnoreCase)) { + return SpatialGeneratorTypes::POINT; + } else { + cout << "TRestGeant4PrimaryGeneratorTypes::StringToSpatialGeneratorTypes - Error - Unknown " + "SpatialGeneratorTypes: " + << type << endl; + exit(1); + } +} + +string TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorShapesToString(const SpatialGeneratorShapes& type) { + switch (type) { + case SpatialGeneratorShapes::GDML: + return "GDML"; + case SpatialGeneratorShapes::WALL: + return "Wall"; + case SpatialGeneratorShapes::CIRCLE: + return "Circle"; + case SpatialGeneratorShapes::BOX: + return "Box"; + case SpatialGeneratorShapes::SPHERE: + return "Sphere"; + case SpatialGeneratorShapes::CYLINDER: + return "Cylinder"; + } + cout << "TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorShapesToString - Error - Unknown " + "SpatialGeneratorShapes" + << endl; + exit(1); +} + +SpatialGeneratorShapes TRestGeant4PrimaryGeneratorTypes::StringToSpatialGeneratorShapes(const string& type) { + if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::GDML), + TString::ECaseCompare::kIgnoreCase)) { + return SpatialGeneratorShapes::GDML; + } else if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::WALL), + TString::ECaseCompare::kIgnoreCase)) { + return SpatialGeneratorShapes::WALL; + } else if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::CIRCLE), + TString::ECaseCompare::kIgnoreCase)) { + return SpatialGeneratorShapes::CIRCLE; + } else if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::BOX), + TString::ECaseCompare::kIgnoreCase)) { + return SpatialGeneratorShapes::BOX; + } else if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::SPHERE), + TString::ECaseCompare::kIgnoreCase)) { + return SpatialGeneratorShapes::SPHERE; + } else if (TString(type).EqualTo(SpatialGeneratorShapesToString(SpatialGeneratorShapes::CYLINDER), + TString::ECaseCompare::kIgnoreCase)) { + return SpatialGeneratorShapes::CYLINDER; + + } else { + cout << "TRestGeant4PrimaryGeneratorTypes::StringToSpatialGeneratorShapes - Error - Unknown " + "SpatialGeneratorShapes: " + << type << endl; + exit(1); + } +} + +string TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypesToString( + const EnergyDistributionTypes& type) { + switch (type) { + case EnergyDistributionTypes::TH1D: + return "TH1D"; + case EnergyDistributionTypes::FORMULA: + return "Formula"; + case EnergyDistributionTypes::MONO: + return "Mono"; + case EnergyDistributionTypes::FLAT: + return "Flat"; + case EnergyDistributionTypes::LOG: + return "Log"; + } + cout << "TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypesToString - Error - Unknown energy " + "distribution type" + << endl; + exit(1); +} + +EnergyDistributionTypes TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionTypes( + const string& type) { + if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::TH1D), + TString::ECaseCompare::kIgnoreCase)) { + return EnergyDistributionTypes::TH1D; + } else if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::FORMULA), + TString::ECaseCompare::kIgnoreCase)) { + return EnergyDistributionTypes::FORMULA; + } else if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::MONO), + TString::ECaseCompare::kIgnoreCase)) { + return EnergyDistributionTypes::MONO; + } else if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::FLAT), + TString::ECaseCompare::kIgnoreCase)) { + return EnergyDistributionTypes::FLAT; + } else if (TString(type).EqualTo(EnergyDistributionTypesToString(EnergyDistributionTypes::LOG), + TString::ECaseCompare::kIgnoreCase)) { + return EnergyDistributionTypes::LOG; + } else { + cout << "TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionTypes - Error - Unknown " + "EnergyDistributionTypes: " + << type << endl; + exit(1); + } +} + +string TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToString( + const EnergyDistributionFormulas& type) { + switch (type) { + case EnergyDistributionFormulas::COSMIC_MUONS: + return "CosmicMuons"; + case EnergyDistributionFormulas::COSMIC_NEUTRONS: + return "CosmicNeutrons"; + case EnergyDistributionFormulas::COSMIC_GAMMAS: + return "CosmicGammas"; + } + cout << "TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToString - Error - Unknown energy " + "distribution formula" + << endl; + exit(1); +} + +EnergyDistributionFormulas TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionFormulas( + const string& type) { + if (TString(type).EqualTo(EnergyDistributionFormulasToString(EnergyDistributionFormulas::COSMIC_MUONS), + TString::ECaseCompare::kIgnoreCase)) { + return EnergyDistributionFormulas::COSMIC_MUONS; + } else if (TString(type).EqualTo( + EnergyDistributionFormulasToString(EnergyDistributionFormulas::COSMIC_NEUTRONS), + TString::ECaseCompare::kIgnoreCase)) { + return EnergyDistributionFormulas::COSMIC_NEUTRONS; + } else if (TString(type).EqualTo( + EnergyDistributionFormulasToString(EnergyDistributionFormulas::COSMIC_GAMMAS), + TString::ECaseCompare::kIgnoreCase)) { + return EnergyDistributionFormulas::COSMIC_GAMMAS; + } else { + cout << "TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionFormulas - Error - Unknown " + "energyDistributionFormulas: " + << type << endl; + exit(1); + } +} + +TF1 TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToRootFormula( + const EnergyDistributionFormulas& type) { + switch (type) { + case EnergyDistributionFormulas::COSMIC_MUONS: + exit(1); + case EnergyDistributionFormulas::COSMIC_NEUTRONS: { + // Formula from https://ieeexplore.ieee.org/document/1369506 + const char* title = "Cosmic Neutrons at Sea Level"; + auto distribution = + TF1(title, + "1.006E-6 * TMath::Exp(-0.3500 * TMath::Power(TMath::Log(x * 1E-3), 2) + 2.1451 * " + "TMath::Log(x * 1E-3)) + " + "1.011E-3 * TMath::Exp(-0.4106 * TMath::Power(TMath::Log(x * 1E-3), 2) - 0.6670 * " + "TMath::Log(x * 1E-3))", + 1E2, 1E7); + distribution.SetNormalized(true); // Normalized, not really necessary + distribution.SetTitle(title); + distribution.GetXaxis()->SetTitle("Energy (keV)"); + return distribution; + } + case EnergyDistributionFormulas::COSMIC_GAMMAS: + exit(1); + } + cout << "TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToRootFormula - Error - Unknown " + "energy distribution formula" + << endl; + exit(1); +} + +string TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypesToString( + const AngularDistributionTypes& type) { + switch (type) { + case AngularDistributionTypes::TH1D: + return "TH1D"; + case AngularDistributionTypes::FORMULA: + return "Formula"; + case AngularDistributionTypes::ISOTROPIC: + return "Isotropic"; + case AngularDistributionTypes::FLUX: + return "Flux"; + case AngularDistributionTypes::BACK_TO_BACK: + return "Back to back"; + } + cout << "TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypesToString - Error - Unknown angular " + "distribution type" + << endl; + exit(1); +} + +AngularDistributionTypes TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes( + const string& type) { + if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::TH1D), + TString::ECaseCompare::kIgnoreCase)) { + return AngularDistributionTypes::TH1D; + } else if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::FORMULA), + TString::ECaseCompare::kIgnoreCase)) { + return AngularDistributionTypes::FORMULA; + } else if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::ISOTROPIC), + TString::ECaseCompare::kIgnoreCase)) { + return AngularDistributionTypes::ISOTROPIC; + } else if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::FLUX), + TString::ECaseCompare::kIgnoreCase)) { + return AngularDistributionTypes::FLUX; + } else if (TString(type).EqualTo(AngularDistributionTypesToString(AngularDistributionTypes::BACK_TO_BACK), + TString::ECaseCompare::kIgnoreCase)) { + return AngularDistributionTypes::BACK_TO_BACK; + } else { + cout << "TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes - Error - Unknown " + "AngularDistributionTypes: " + << type << endl; + exit(1); + } +} + +string TRestGeant4PrimaryGeneratorTypes::AngularDistributionFormulasToString( + const AngularDistributionFormulas& type) { + switch (type) { + case AngularDistributionFormulas::COS2: + return "Cos2"; + case AngularDistributionFormulas::COS3: + return "Cos3"; + } + cout << "TRestGeant4PrimaryGeneratorTypes::AngularDistributionFormulasToString - Error - Unknown angular " + "distribution formula" + << endl; + exit(1); +} + +AngularDistributionFormulas TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionFormulas( + const string& type) { + if (TString(type).EqualTo(AngularDistributionFormulasToString(AngularDistributionFormulas::COS2), + TString::ECaseCompare::kIgnoreCase)) { + return AngularDistributionFormulas::COS2; + } else if (TString(type).EqualTo(AngularDistributionFormulasToString(AngularDistributionFormulas::COS3), + TString::ECaseCompare::kIgnoreCase)) { + return AngularDistributionFormulas::COS3; + } else { + cout << "TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionFormulas - Error - Unknown " + "AngularDistributionFormulas: " + << type << endl; + exit(1); + } +} + +TF1 TRestGeant4PrimaryGeneratorTypes::AngularDistributionFormulasToRootFormula( + const AngularDistributionFormulas& type) { + switch (type) { + case AngularDistributionFormulas::COS2: { + auto cos2 = [](double* xs, double* ps) { + if (xs[0] >= 0 && xs[0] <= TMath::Pi() / 2) { + return TMath::Power(TMath::Cos(xs[0]), 2); + } + return 0.0; + }; + const char* title = "AngularDistribution: Cos2"; + auto f = TF1(title, cos2, 0.0, TMath::Pi()); + f.SetTitle(title); + return f; + } + case AngularDistributionFormulas::COS3: { + auto cos3 = [](double* xs, double* ps) { + if (xs[0] >= 0 && xs[0] <= TMath::Pi() / 2) { + return TMath::Power(TMath::Cos(xs[0]), 3); + } + return 0.0; + }; + const char* title = "AngularDistribution: Cos3"; + auto f = TF1(title, cos3, 0.0, TMath::Pi()); + f.SetTitle(title); + return f; + } + } + cout << "TRestGeant4PrimaryGeneratorTypes::AngularDistributionFormulasToRootFormula - Error - Unknown " + "angular distribution formula" + << endl; + exit(1); +} + +void TRestGeant4PrimaryGeneratorInfo::Print() const { + const auto typeEnum = StringToSpatialGeneratorTypes(fSpatialGeneratorType.Data()); + RESTMetadata << "Generator type: " << SpatialGeneratorTypesToString(typeEnum) << RESTendl; + + if (typeEnum != SpatialGeneratorTypes::POINT) { + const auto shapeEnum = StringToSpatialGeneratorShapes(fSpatialGeneratorShape.Data()); + RESTMetadata << "Generator shape: " << SpatialGeneratorShapesToString(shapeEnum); + if (shapeEnum == SpatialGeneratorShapes::GDML) { + RESTMetadata << "::" << fSpatialGeneratorFrom << RESTendl; + } else { + if (shapeEnum == SpatialGeneratorShapes::BOX) { + RESTMetadata << ", (length, width, height): "; + } else if (shapeEnum == SpatialGeneratorShapes::SPHERE) { + RESTMetadata << ", (radius, , ): "; + } else if (shapeEnum == SpatialGeneratorShapes::WALL) { + RESTMetadata << ", (length, width, ): "; + } else if (shapeEnum == SpatialGeneratorShapes::CIRCLE) { + RESTMetadata << ", (radius, , ): "; + } else if (shapeEnum == SpatialGeneratorShapes::CYLINDER) { + RESTMetadata << ", (radius, length, ): "; + } + + RESTMetadata << fSpatialGeneratorSize.X() << ", " << fSpatialGeneratorSize.Y() << ", " + << fSpatialGeneratorSize.Z() << RESTendl; + } + } + RESTMetadata << "Generator center : (" << fSpatialGeneratorPosition.X() << "," + << fSpatialGeneratorPosition.Y() << "," << fSpatialGeneratorPosition.Z() << ") mm" + << RESTendl; + RESTMetadata << "Generator rotation : (" << fSpatialGeneratorRotationAxis.X() << "," + << fSpatialGeneratorRotationAxis.Y() << "," << fSpatialGeneratorRotationAxis.Z() + << "), angle: " << fSpatialGeneratorRotationValue * TMath::RadToDeg() << "ยบ" << RESTendl; +} diff --git a/src/TRestGeant4Track.cxx b/src/TRestGeant4Track.cxx index 770d169..8f5a322 100644 --- a/src/TRestGeant4Track.cxx +++ b/src/TRestGeant4Track.cxx @@ -24,13 +24,9 @@ using namespace std; ClassImp(TRestGeant4Track); -TRestGeant4Track::TRestGeant4Track() { - // TRestGeant4Track default constructor -} +TRestGeant4Track::TRestGeant4Track() = default; -TRestGeant4Track::~TRestGeant4Track() { - // TRestGeant4Track destructor -} +TRestGeant4Track::~TRestGeant4Track() = default; Int_t TRestGeant4Track::GetProcessID(const TString& processName) const { const TRestGeant4Metadata* metadata = GetGeant4Metadata(); @@ -86,55 +82,63 @@ EColor TRestGeant4Track::GetParticleColor() const { /// the hits of that specific volume will be counted. /// size_t TRestGeant4Track::GetNumberOfHits(Int_t volID) const { - size_t hits = 0; + size_t numberOfHits = 0; for (int n = 0; n < fHits.GetNumberOfHits(); n++) { - if (volID != -1 && fHits.GetVolumeId(n) != volID) continue; - hits++; + if (volID != -1 && fHits.GetVolumeId(n) != volID) { + continue; + } + numberOfHits++; } - return hits; + return numberOfHits; } -Double_t TRestGeant4Track::GetTrackLength() const { - Double_t length = 0; - - length = GetDistance(fHits.GetPosition(0), GetTrackOrigin()); - - for (int i = 1; i < GetNumberOfHits(); i++) { - TVector3 prevHit = fHits.GetPosition(i - 1); - TVector3 hit = fHits.GetPosition(i); - length += GetDistance(hit, prevHit); +/////////////////////////////////////////////// +/// \brief Function that returns the number of hit depositions found inside +/// the TRestGeant4Track with energy > 0. If a specific volume id is given as argument only +/// the hits of that specific volume will be counted. +/// +size_t TRestGeant4Track::GetNumberOfPhysicalHits(Int_t volID) const { + size_t numberOfHits = 0; + for (int n = 0; n < fHits.GetNumberOfHits(); n++) { + if (volID != -1 && fHits.GetVolumeId(n) != volID) { + continue; + } + if (fHits.GetEnergy(n) <= 0) { + continue; + } + numberOfHits++; } - return length; + return numberOfHits; } void TRestGeant4Track::PrintTrack(size_t maxHits) const { - cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" << endl; - cout.precision(10); - cout << " SubEvent ID : " << fSubEventId << " Global timestamp : " << GetGlobalTime() << " seconds" - << endl; - cout.precision(5); - cout << " Track ID : " << GetTrackID() << " Parent ID : " << GetParentID() - << " Created by process: " << fCreatorProcess; - cout << " Particle : " << GetParticleName() << " Time track length : " << GetTrackTimeLength() << " us" - << endl; - cout << " Origin : X = " << GetTrackOrigin().X() << "mm Y = " << GetTrackOrigin().Y() - << "mm Z = " << GetTrackOrigin().Z() << "mm Ekin : " << GetKineticEnergy() << " keV" << endl; - cout << "++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++" - "++++++++++++" - << endl; - - int nHits = GetNumberOfHits(); - if (maxHits > 0) { - nHits = min(maxHits, GetNumberOfHits()); - cout << " Printing only the first " << nHits << " hits of the track" << endl; + cout + << " * TrackID: " << fTrackID << " - Particle: " << fParticleName << " - ParentID: " << fParentID + << "" + << (GetParentTrack() != nullptr + ? TString::Format(" - Parent particle: %s", GetParentTrack()->GetParticleName().Data()).Data() + : "") + << " - Created by '" << fCreatorProcess + << "' with initial " + "KE of " + << ToEnergyString(fInitialKineticEnergy) << "" << endl; + + cout << " Initial position " << VectorToString(fInitialPosition) << " mm at time " + << ToTimeString(fGlobalTimestamp) << " - Time length of " << ToTimeString(fTimeLength) + << " and spatial length of " << ToLengthString(fLength) << endl; + + size_t nHits = GetNumberOfHits(); + if (maxHits > 0 && maxHits < nHits) { + nHits = min(maxHits, nHits); + cout << "Printing only the first " << nHits << " hits of the track" << endl; } const TRestGeant4Metadata* metadata = GetGeant4Metadata(); for (int i = 0; i < nHits; i++) { TString processName = GetProcessName(fHits.GetHitProcess(i)); if (processName.IsNull()) { - // in case process name is not found, use ID - processName = TString(std::to_string(fHits.GetHitProcess(i))); + processName = + TString(std::to_string(fHits.GetHitProcess(i))); // in case process name is not found, use ID } TString volumeName = ""; @@ -145,15 +149,12 @@ void TRestGeant4Track::PrintTrack(size_t maxHits) const { // in case process name is not found, use ID volumeName = TString(std::to_string(fHits.GetHitVolume(i))); } - - cout << " # Hit " << i << " # process : " << processName << " volume : " << volumeName - << " X : " << fHits.GetX(i) << " Y : " << fHits.GetY(i) << " Z : " << fHits.GetZ(i) << " mm" - << endl; - cout << " Edep : " << fHits.GetEnergy(i) << " keV Ekin : " << fHits.GetKineticEnergy(i) << " keV" - << " Global time : " << fHits.GetTime(i) << " us" << endl; + cout << " - Hit " << i << " - Energy: " << ToEnergyString(fHits.GetEnergy(i)) + << " - Process: " << processName << " - Volume: " << volumeName + << " - Position: " << VectorToString(TVector3(fHits.GetX(i), fHits.GetY(i), fHits.GetZ(i))) + << " mm - Time: " << ToTimeString(fHits.GetTime(i)) + << " - KE: " << ToEnergyString(fHits.GetKineticEnergy(i)) << endl; } - cout << endl; - cout.precision(2); } Bool_t TRestGeant4Track::ContainsProcessInVolume(Int_t processID, Int_t volumeID) const { @@ -183,3 +184,79 @@ const TRestGeant4Metadata* TRestGeant4Track::GetGeant4Metadata() const { } return GetEvent()->GetGeant4Metadata(); } + +TRestGeant4Track* TRestGeant4Track::GetParentTrack() const { + if (fEvent == nullptr) { + return nullptr; + } + return fEvent->GetTrackByID(GetParentID()); +} + +vector TRestGeant4Track::GetSecondaryTracks() const { + vector secondaryTracks = {}; + for (const auto trackID : fSecondaryTrackIDs) { + const TRestGeant4Track* track = fEvent->GetTrackByID(trackID); + if (track != nullptr) { + secondaryTracks.push_back(track); + } + } + return secondaryTracks; +} + +TString TRestGeant4Track::GetInitialVolume() const { + const auto metadata = GetGeant4Metadata(); + if (metadata == nullptr) { + return ""; + } + const auto& hits = GetHits(); + return GetGeant4Metadata()->GetGeant4GeometryInfo().GetVolumeFromID(hits.GetVolumeId(0)); +} + +TString TRestGeant4Track::GetFinalVolume() const { + const auto metadata = GetGeant4Metadata(); + if (metadata == nullptr) { + return ""; + } + const auto& hits = GetHits(); + return GetGeant4Metadata()->GetGeant4GeometryInfo().GetVolumeFromID( + hits.GetVolumeId(hits.GetNumberOfHits() - 1)); +} + +Double_t TRestGeant4Track::GetEnergyInVolume(const TString& volumeName, bool children) const { + const auto metadata = GetGeant4Metadata(); + if (metadata == nullptr) { + return 0; + } + + const auto volumeId = metadata->GetGeant4GeometryInfo().GetIDFromVolume(volumeName); + + if (!children) { + return GetEnergyInVolume(volumeId); + } + + Double_t energy = 0; + vector tracks = {this}; + while (!tracks.empty()) { + const TRestGeant4Track* track = tracks.back(); + tracks.pop_back(); + if (track == nullptr) { + continue; + } + energy += track->GetEnergyInVolume(volumeId); + for (const TRestGeant4Track* secondaryTrack : track->GetSecondaryTracks()) { + tracks.push_back(secondaryTrack); + } + } + return energy; +} + +TString TRestGeant4Track::GetLastProcessName() const { + const auto metadata = GetGeant4Metadata(); + if (metadata == nullptr) { + return ""; + } + + const auto& hits = GetHits(); + return GetGeant4Metadata()->GetGeant4PhysicsInfo().GetProcessName( + hits.GetProcess(hits.GetNumberOfHits() - 1)); +} diff --git a/src/TRestGeant4VetoAnalysisProcess.cxx b/src/TRestGeant4VetoAnalysisProcess.cxx index a112743..f8e29db 100644 --- a/src/TRestGeant4VetoAnalysisProcess.cxx +++ b/src/TRestGeant4VetoAnalysisProcess.cxx @@ -72,14 +72,13 @@ ClassImp(TRestGeant4VetoAnalysisProcess); TRestGeant4VetoAnalysisProcess::TRestGeant4VetoAnalysisProcess() { Initialize(); } TRestGeant4VetoAnalysisProcess::TRestGeant4VetoAnalysisProcess(const char* configFilename) { - Initialize(); - if (LoadConfigFromFile(configFilename)) LoadDefaultConfig(); + TRestGeant4VetoAnalysisProcess(); + if (LoadConfigFromFile(configFilename)) { + LoadDefaultConfig(); + } } -/////////////////////////////////////////////// -/// \brief Default destructor -/// -TRestGeant4VetoAnalysisProcess::~TRestGeant4VetoAnalysisProcess() { delete fOutputG4Event; } +TRestGeant4VetoAnalysisProcess::~TRestGeant4VetoAnalysisProcess() { delete fOutputEvent; } /////////////////////////////////////////////// /// \brief Function to load the default config in absence of RML input @@ -91,13 +90,10 @@ void TRestGeant4VetoAnalysisProcess::LoadDefaultConfig() { SetTitle("Default con /// section name /// void TRestGeant4VetoAnalysisProcess::Initialize() { - fG4Metadata = nullptr; - SetSectionName(this->ClassName()); SetLibraryVersion(LIBRARY_VERSION); - fInputG4Event = nullptr; - fOutputG4Event = new TRestGeant4Event(); + fOutputEvent = new TRestGeant4Event(); } /////////////////////////////////////////////// @@ -120,34 +116,75 @@ void TRestGeant4VetoAnalysisProcess::LoadConfig(const string& configFilename, co /// \brief Process initialization. /// void TRestGeant4VetoAnalysisProcess::InitProcess() { - fG4Metadata = GetMetadata(); - // CAREFUL THIS METHOD IS CALLED TWICE! - // we need to reset these variables to zero - fVetoVolumeIds.clear(); - fVetoGroupVolumeNames.clear(); - - // get "veto" volumes - for (unsigned int i = 0; i < fG4Metadata->GetNumberOfActiveVolumes(); i++) { - string volume_name = (string)fG4Metadata->GetActiveVolumeName(i); - volume_name = clean_string(volume_name); - if (volume_name.find(clean_string(fVetoKeyword)) != string::npos) { - fVetoVolumeIds.push_back(i); + fVetoVolumes.clear(); + fVetoDetectorVolumes.clear(); + fVetoDetectorBoundaryDirection.clear(); + fVetoDetectorBoundaryPosition.clear(); + + if (fGeant4Metadata == nullptr) { + // maybe it was manually initialized + fGeant4Metadata = GetMetadata(); + } + if (fGeant4Metadata == nullptr) { + cerr << "TRestGeant4VetoAnalysisProcess::InitProcess: Geant4 metadata not found" << endl; + exit(1); + } + + const auto& geometryInfo = fGeant4Metadata->GetGeant4GeometryInfo(); + + fVetoVolumes = geometryInfo.GetAllPhysicalVolumesMatchingExpression(fVetoVolumesExpression); + if (fVetoVolumes.empty()) { + const auto logicalVolumes = + geometryInfo.GetAllLogicalVolumesMatchingExpression(fVetoVolumesExpression); + for (const auto& logicalVolume : logicalVolumes) { + for (const auto& physicalVolume : geometryInfo.GetAllPhysicalVolumesFromLogical(logicalVolume)) { + fVetoVolumes.push_back(geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolume)); + } } } + if (fVetoVolumes.empty()) { + cerr << "TRestGeant4VetoAnalysisProcess::InitProcess: No veto volumes found" << endl; + exit(1); + } - // veto groups (fill fVetoGroupVolumeNames) - for (unsigned int i = 0; i < fVetoGroupKeywords.size(); i++) { - string veto_group_keyword = clean_string(fVetoGroupKeywords[i]); - fVetoGroupVolumeNames[veto_group_keyword] = std::vector{}; - for (int& id : fVetoVolumeIds) { - string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); - volume_name = clean_string(volume_name); - if (volume_name.find(veto_group_keyword) != string::npos) { - fVetoGroupVolumeNames[veto_group_keyword].push_back( - (string)fG4Metadata->GetActiveVolumeName(id)); + // get detector volumes if requested + if (!fVetoDetectorsExpression.IsNull()) { + fVetoDetectorVolumes = geometryInfo.GetAllPhysicalVolumesMatchingExpression(fVetoDetectorsExpression); + if (fVetoDetectorVolumes.empty()) { + const auto logicalVolumes = + geometryInfo.GetAllLogicalVolumesMatchingExpression(fVetoDetectorsExpression); + for (const auto& logicalVolume : logicalVolumes) { + for (const auto& physicalVolume : + geometryInfo.GetAllPhysicalVolumesFromLogical(logicalVolume)) { + fVetoDetectorVolumes.push_back( + geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolume)); + } } } + if (fVetoDetectorVolumes.empty()) { + cerr << "TRestGeant4VetoAnalysisProcess::InitProcess: No detector volumes found" << endl; + exit(1); + } + if (fVetoDetectorVolumes.size() != fVetoVolumes.size()) { + cerr << "TRestGeant4VetoAnalysisProcess::InitProcess: Number of detector volumes " + << "does not match number of veto volumes" << endl; + exit(1); + } + } + + for (int i = 0; i < fVetoDetectorVolumes.size(); i++) { + const auto& vetoName = fVetoVolumes[i]; + const auto& vetoPosition = geometryInfo.GetPosition(vetoName); + + const auto& vetoDetectorName = fVetoDetectorVolumes[i]; + const auto& vetoDetectorPosition = geometryInfo.GetPosition(vetoDetectorName); + + const auto distance = vetoDetectorPosition - vetoPosition; + const auto direction = distance.Unit(); + + fVetoDetectorBoundaryDirection[vetoName] = direction; + fVetoDetectorBoundaryPosition[vetoName] = vetoDetectorPosition - direction * fVetoDetectorOffsetSize; } PrintMetadata(); @@ -157,14 +194,15 @@ void TRestGeant4VetoAnalysisProcess::InitProcess() { /// \brief The main processing event function /// TRestEvent* TRestGeant4VetoAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) { - fInputG4Event = (TRestGeant4Event*)inputEvent; - *fOutputG4Event = *((TRestGeant4Event*)inputEvent); + /* + fInputEvent = (TRestGeant4Event*)inputEvent; + *fOutputEvent = *((TRestGeant4Event*)inputEvent); std::map volume_energy_map; for (unsigned int i = 0; i < fVetoVolumeIds.size(); i++) { int id = fVetoVolumeIds[i]; - string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); + string volume_name = (string)fGeant4Metadata->GetActiveVolumeName(id); Double_t energy = fOutputG4Event->GetEnergyDepositedInVolume(id); volume_energy_map[volume_name] = energy; @@ -214,7 +252,7 @@ TRestEvent* TRestGeant4VetoAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) auto track = fOutputG4Event->GetTrack(i); string particle_name = (string)track.GetParticleName(); for (const auto& id : fVetoVolumeIds) { - string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); + string volume_name = (string)fGeant4Metadata->GetActiveVolumeName(id); if (particle_name == "e-" || particle_name == "e+" || particle_name == "gamma") { // no quenching factor @@ -261,6 +299,7 @@ TRestEvent* TRestGeant4VetoAnalysisProcess::ProcessEvent(TRestEvent* inputEvent) } return fOutputG4Event; + */ } /////////////////////////////////////////////// @@ -282,29 +321,59 @@ void TRestGeant4VetoAnalysisProcess::EndProcess() { /// void TRestGeant4VetoAnalysisProcess::InitFromConfigFile() { // word to identify active volume as veto (default = "veto" e.g. "vetoTop") - string veto_keyword = GetParameter("vetoKeyword", "veto"); - fVetoKeyword = clean_string(veto_keyword); - // comma separated tags: "top, bottom, ..." - string veto_group_keywords = GetParameter("vetoGroupKeywords", ""); - stringstream ss(veto_group_keywords); - while (ss.good()) { - string substr; - getline(ss, substr, ','); - fVetoGroupKeywords.push_back(clean_string(substr)); + fVetoVolumesExpression = GetParameter("vetoVolumesExpression", fVetoVolumesExpression); + fVetoDetectorsExpression = GetParameter("vetoDetectorsExpression", fVetoDetectorsExpression); + + fVetoDetectorOffsetSize = GetDblParameterWithUnits("vetoDetectorOffset", fVetoDetectorOffsetSize); + fVetoLightAttenuation = GetDblParameterWithUnits("vetoLightAttenuation", fVetoLightAttenuation); + fVetoQuenchingFactor = GetDblParameterWithUnits("quenchingFactor", fVetoQuenchingFactor); +} + +void TRestGeant4VetoAnalysisProcess::PrintMetadata() { + BeginPrintProcess(); + + cout << "Veto volume expression: " << fVetoVolumesExpression << endl; + if (!fVetoDetectorsExpression.IsNull()) { + cout << "Veto detector expression: " << fVetoDetectorsExpression << endl; + cout << "Veto detector offset: " << fVetoDetectorOffsetSize << endl; + cout << "Veto light attenuation: " << fVetoLightAttenuation << endl; + } else { + cout << "Veto detector expression: not set" << endl; + } + cout << "Veto quenching factor: " << fVetoQuenchingFactor << endl; + + RESTDebug << RESTendl; + + if (fVetoVolumes.empty()) { + cout << "Process not initialized yet" << endl; + return; } - // comma separated quenching factors: "0.15, 1.00, ..." - string quenching_factors = GetParameter("vetoQuenchingFactors", "-1"); - stringstream ss_qf(quenching_factors); - while (ss_qf.good()) { - string substr; - getline(ss_qf, substr, ','); - substr = clean_string(substr); - Float_t quenching_factor = (Float_t)std::atof(substr.c_str()); - if (quenching_factor > 1 || quenching_factor < 0) { - cout << "ERROR: quenching factor must be between 0 and 1" << endl; + cout << "Number of veto volumes: " << fVetoVolumes.size() << endl; + cout << "Number of veto detector volumes: " << fVetoDetectorVolumes.size() << endl; + + const auto& geometryInfo = fGeant4Metadata->GetGeant4GeometryInfo(); + for (int i = 0; i < fVetoVolumes.size(); i++) { + const auto& vetoName = fVetoVolumes[i]; + const auto& vetoPosition = geometryInfo.GetPosition(vetoName); + + cout << TString::Format(" - Veto volume: %d - name: '%s' - position: %s mm", i, vetoName.Data(), + VectorToString(vetoPosition).c_str()) + << endl; + + if (fVetoDetectorVolumes.empty()) { continue; } - fQuenchingFactors.push_back(quenching_factor); + + const auto& vetoDetectorName = fVetoDetectorVolumes[i]; + const auto& vetoDetectorPosition = geometryInfo.GetPosition(vetoDetectorName); + + cout << TString::Format(" Veto detector name: '%s' - position: %s mm", vetoDetectorName.Data(), + VectorToString(vetoDetectorPosition).c_str()) + << endl; + + cout << TString::Format(" Boundary position: %s mm - direction: %s", + VectorToString(fVetoDetectorBoundaryPosition.at(vetoName)).c_str(), + VectorToString(fVetoDetectorBoundaryDirection.at(vetoName)).c_str()); } } diff --git a/test/files/TRestGeant4Example.rml b/test/files/TRestGeant4Example.rml index 705c01d..175b79f 100644 --- a/test/files/TRestGeant4Example.rml +++ b/test/files/TRestGeant4Example.rml @@ -1,4 +1,4 @@ - + @@ -9,18 +9,25 @@ - + - - + + - - - - - + + + + + + + + + + + + diff --git a/test/files/TRestGeant4VetoAnalysisProcessExample.rml b/test/files/TRestGeant4VetoAnalysisProcessExample.rml new file mode 100644 index 0000000..c7754ce --- /dev/null +++ b/test/files/TRestGeant4VetoAnalysisProcessExample.rml @@ -0,0 +1,16 @@ + + + + + + + + + + + + + + + + diff --git a/test/src/Geant4AnalysisProcesses.cxx b/test/src/Geant4AnalysisProcesses.cxx new file mode 100644 index 0000000..ab68718 --- /dev/null +++ b/test/src/Geant4AnalysisProcesses.cxx @@ -0,0 +1,77 @@ + +#include +#include + +#include + +namespace fs = std::filesystem; + +using namespace std; + +const auto filesPath = fs::path(__FILE__).parent_path().parent_path() / "files"; +const auto vetoAnalysisRml = filesPath / "TRestGeant4VetoAnalysisProcessExample.rml"; +const auto vetoAnalysisRestG4Run = filesPath / "VetoAnalysisGeant4Run.root"; + +TEST(TRestGeant4VetoAnalysisProcess, TestFiles) { + GTEST_SKIP(); + + cout << "Test files path: " << filesPath << endl; + + // Check dir exists and is a directory + EXPECT_TRUE(fs::is_directory(filesPath)); + // Check it's not empty + EXPECT_TRUE(!fs::is_empty(filesPath)); + + // All used files in this tests + EXPECT_TRUE(fs::exists(vetoAnalysisRml)); + EXPECT_TRUE(fs::exists(vetoAnalysisRestG4Run)); +} + +TEST(TRestGeant4VetoAnalysisProcess, Default) { + TRestGeant4VetoAnalysisProcess process; + + process.PrintMetadata(); + + EXPECT_TRUE(process.GetVetoVolumesExpression().IsNull()); + EXPECT_TRUE(process.GetVetoDetectorExpression().IsNull()); + EXPECT_TRUE(process.GetVetoDetectorOffsetSize() == 0); + EXPECT_TRUE(process.GetVetoLightAttenuation() == 0); + EXPECT_TRUE(process.GetVetoQuenchingFactor() == 1.0); +} + +TEST(TRestGeant4VetoAnalysisProcess, FromRml) { + TRestGeant4VetoAnalysisProcess process(vetoAnalysisRml.c_str()); + + process.PrintMetadata(); + + EXPECT_TRUE(process.GetVetoVolumesExpression() == "^scintillatorVolume"); + EXPECT_TRUE(process.GetVetoDetectorExpression() == "^scintillatorLightGuideVolume"); + EXPECT_TRUE(process.GetVetoDetectorOffsetSize() == 0); + EXPECT_TRUE(process.GetVetoLightAttenuation() == 0); + EXPECT_TRUE(process.GetVetoQuenchingFactor() == 0); +} + +TEST(TRestGeant4VetoAnalysisProcess, Simulation) { + GTEST_SKIP(); + + TRestGeant4VetoAnalysisProcess process(vetoAnalysisRml.c_str()); + + EXPECT_TRUE(process.GetVetoVolumesExpression() == "^scintillatorVolume"); + EXPECT_TRUE(process.GetVetoDetectorExpression() == "^scintillatorLightGuideVolume"); + EXPECT_TRUE(process.GetVetoDetectorOffsetSize() == 0); + EXPECT_TRUE(process.GetVetoLightAttenuation() == 0); + EXPECT_TRUE(process.GetVetoQuenchingFactor() == 0); + + TRestRun run(vetoAnalysisRestG4Run.c_str()); + run.GetInputFile()->ls(); + + const auto metadata = + dynamic_cast(run.GetMetadataClass("TRestGeant4Metadata")); + EXPECT_TRUE(metadata != nullptr); + + process.SetGeant4Metadata(metadata); + + process.InitProcess(); + + process.PrintMetadata(); +} diff --git a/test/src/TRestGeant4Metadata.cxx b/test/src/TRestGeant4Metadata.cxx index 1bf09df..005b8f5 100644 --- a/test/src/TRestGeant4Metadata.cxx +++ b/test/src/TRestGeant4Metadata.cxx @@ -12,7 +12,7 @@ const auto filesPath = fs::path(__FILE__).parent_path().parent_path() / "files"; const auto geant4MetadataRml = filesPath / "TRestGeant4Example.rml"; TEST(TRestGeant4Metadata, TestFiles) { - cout << "FrameworkCore test files path: " << filesPath << endl; + cout << "Test files path: " << filesPath << endl; // Check dir exists and is a directory EXPECT_TRUE(fs::is_directory(filesPath)); @@ -29,14 +29,15 @@ TEST(TRestGeant4Metadata, Default) { restGeant4Metadata.PrintMetadata(); EXPECT_TRUE(restGeant4Metadata.GetSeed() == 0); - EXPECT_TRUE(restGeant4Metadata.GetSensitiveVolume() == "gas"); } TEST(TRestGeant4Metadata, FromRml) { - GTEST_SKIP_("Problem with paths..."); + cout << "Path: " << geant4MetadataRml << endl; TRestGeant4Metadata restGeant4Metadata(geant4MetadataRml.c_str()); + GTEST_SKIP_("TODO: fix this"); + restGeant4Metadata.PrintMetadata(); EXPECT_TRUE(restGeant4Metadata.GetSensitiveVolume() == "sensitiveVolume"); @@ -46,5 +47,5 @@ TEST(TRestGeant4Metadata, FromRml) { EXPECT_TRUE(restGeant4Metadata.GetNumberOfSources() == 1); const auto particleSource = restGeant4Metadata.GetParticleSource(0); EXPECT_TRUE(particleSource->GetParticleName() == "geantino"); - EXPECT_TRUE(particleSource->GetEnergyDistType() == "mono"); + EXPECT_TRUE(particleSource->GetEnergyDistributionType() == "mono"); }