Merge pull request #129 from rest-for-physics/cosmics-benchmark

Comparison of cosmic generator with disk

CMakeLists.txt

@@ -2,16 +2,8 @@ set_library_version(LibraryVersion)
 add_definitions(-DLIBRARY_VERSION="${LibraryVersion}")
 
 if (${REST_DECAY0} MATCHES "ON")
-    add_definitions(-DUSE_Decay0)
-    # TODO Issue #6 at rest-framework
-    # ----------------------------------------------------------------------------
-    # Find package Decay0 and ROOT find_package(BxDecay0 1.0.9 CONFIG COMPONENTS
-    # manager REQUIRED)
-    #
-    # find_package was not working properly for me in 1.0.9. So I was using
-    # bxdecay0-config. But in case Decay0 flag is enabled and we do not have a
-    # bxdecay0 installation we might run into problems
-    #
+    add_definitions(-DUSE_Decay0
+    )# https://github.com/rest-for-physics/framework/issues/6
 
     execute_process(
         COMMAND bxdecay0-config --version

inc/TRestGeant4Event.h

@@ -119,6 +119,9 @@ class TRestGeant4Event : public TRestEvent {
     Double_t fTotalDepositedEnergy = 0;
     Double_t fSensitiveVolumeEnergy = 0;
 
+    double fEventTimeWall = 0;
+    double fEventTimeWallPrimaryGeneration = 0;
+
     Int_t fNVolumes;
     std::vector<Int_t> fVolumeStored;
     std::vector<std::string> fVolumeStoredNames;
@@ -150,6 +153,9 @@ class TRestGeant4Event : public TRestEvent {
     size_t GetNumberOfHits(Int_t volID = -1) const;
     size_t GetNumberOfPhysicalHits(Int_t volID = -1) const;
 
+    double GetEventTimeWall() const { return fEventTimeWall; }
+    double GetEventTimeWallPrimaryGeneration() const { return fEventTimeWallPrimaryGeneration; }
+
     inline const std::vector<TRestGeant4Track>& GetTracks() const { return fTracks; }
     inline size_t GetNumberOfTracks() const { return fTracks.size(); }
     inline Int_t GetNumberOfActiveVolumes() const { return fNVolumes; }
@@ -180,6 +186,8 @@ class TRestGeant4Event : public TRestEvent {
         return energyMap[volumeName];
     }
 
+    inline void ClearTracks() { fTracks.clear(); }
+
     TRestHits GetHits(Int_t volID = -1) const;
     inline TRestHits GetHitsInVolume(Int_t volID) const { return GetHits(volID); }
 
@@ -244,7 +252,7 @@ class TRestGeant4Event : public TRestEvent {
     // Destructor
     virtual ~TRestGeant4Event();
 
-    ClassDefOverride(TRestGeant4Event, 8);  // REST event superclass
+    ClassDefOverride(TRestGeant4Event, 9);  // REST event superclass
 
     // restG4
    public:

inc/TRestGeant4Metadata.h

@@ -133,7 +133,7 @@ class TRestGeant4Metadata : public TRestMetadata {
     /// \brief Time before simulation is ended and saved
     Double_t fSimulationMaxTimeSeconds = 0;
 
-    /// \brief The time, in seconds, that the simulation took to complete.
+    /// \brief The time, in seconds, that the simulation took to complete (wall time)
     Double_t fSimulationTime = 0;
 
     /// \brief The seed value used for Geant4 random event generator.
@@ -153,6 +153,9 @@ class TRestGeant4Metadata : public TRestMetadata {
     /// \brief If activated will remove tracks not present in volumes marked as "keep" or "sensitive".
     Bool_t fRemoveUnwantedTracks = false;
 
+    /// \brief Store event tracks (default is true)
+    Bool_t fStoreTracks = 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;
@@ -269,8 +272,6 @@ class TRestGeant4Metadata : public TRestMetadata {
 
     inline Double_t GetSimulationMaxTimeSeconds() const { return fSimulationMaxTimeSeconds; }
 
-    inline Double_t GetSimulationTime() const { return fSimulationTime; }
-
     // Direct access to sources definition in primary generator
     ///////////////////////////////////////////////////////////
     /// Returns the number of primary event sources defined in the generator.
@@ -369,17 +370,26 @@ class TRestGeant4Metadata : public TRestMetadata {
         return result;
     }
 
+    double GetGeneratorSurfaceCm2() const;
+
     Double_t GetCosmicFluxInCountsPerCm2PerSecond() const;
     Double_t GetCosmicIntensityInCountsPerSecond() const;
+
+    /// Returns the equivalent simulated time in seconds (physical time simulated)
     Double_t GetEquivalentSimulatedTime() const;
 
+    /// Returns the total time of the simulation in seconds (wall time)
+    inline Double_t GetSimulationWallTime() const { return fSimulationTime; }
+
     /// Returns a std::string with the name of the active volume with index n
     inline TString GetActiveVolumeName(Int_t n) const { return fActiveVolumes[n]; }
 
     inline std::vector<TString> GetActiveVolumes() const { return fActiveVolumes; }
 
     inline bool GetRemoveUnwantedTracks() const { return fRemoveUnwantedTracks; }
 
+    bool GetStoreTracks() const { return fStoreTracks; }
+
     inline bool GetRemoveUnwantedTracksKeepZeroEnergyTracks() const {
         return fRemoveUnwantedTracksKeepZeroEnergyTracks;
     }
@@ -393,7 +403,7 @@ class TRestGeant4Metadata : public TRestMetadata {
 
     void SetActiveVolume(const TString& name, Double_t chance, Double_t maxStep = 0);
 
-    void SetSimulationTime(Double_t time) { fSimulationTime = time; }
+    void SetSimulationWallTime(Double_t time) { fSimulationTime = time; }
 
     void PrintMetadata() override;
 
@@ -407,7 +417,7 @@ class TRestGeant4Metadata : public TRestMetadata {
     TRestGeant4Metadata(const TRestGeant4Metadata& metadata);
     TRestGeant4Metadata& operator=(const TRestGeant4Metadata& metadata);
 
-    ClassDefOverride(TRestGeant4Metadata, 18);
+    ClassDefOverride(TRestGeant4Metadata, 19);
 
     // Allow modification of otherwise inaccessible / immutable members that shouldn't be modified by the user
     friend class SteppingAction;

inc/TRestGeant4ParticleSourceCosmics.h

@@ -30,6 +30,7 @@ class TRestGeant4ParticleSourceCosmics : public TRestGeant4ParticleSource {
     const char* GetName() const override { return "TRestGeant4ParticleSourceCosmics"; }
 
     std::map<std::string, TH2D*> GetHistogramsTransformed() const { return fHistogramsTransformed; }
+    std::set<std::string> GetParticleNames() const { return fParticleNames; }
 
     ClassDefOverride(TRestGeant4ParticleSourceCosmics, 2);
 };

inc/TRestGeant4PrimaryGeneratorInfo.h

@@ -164,8 +164,8 @@ class TRestGeant4PrimaryGeneratorInfo {
 
     /// \brief Returns cosmic surface term (cm2) for simulation time computation
     inline Double_t GetSpatialGeneratorCosmicSurfaceTermCm2() const {
-        const auto radius = GetSpatialGeneratorCosmicRadius() / 10.;
-        return M_PI * radius * radius;
+        const auto radius = GetSpatialGeneratorCosmicRadius();
+        return M_PI * radius * radius * 0.01;  // cm2
     }
 
     /// \brief Returns the density function of the generator

macros/REST_Geant4_MergeRestG4Files.C

@@ -50,15 +50,15 @@ void REST_Geant4_MergeRestG4Files(const char* outputFilename, const char* inputF
     // open the first file
     TRestGeant4Metadata mergeMetadata;
 
-    auto mergeRun = new TRestRun();
-    mergeRun->SetName("run");
-    mergeRun->SetOutputFileName(outputFilename);
-    mergeRun->FormOutputFile();
-    mergeRun->GetOutputFile()->cd();
-    mergeRun->SetRunType("restG4");
+    TRestRun mergeRun;
+    mergeRun.SetName("run");
+    mergeRun.SetOutputFileName(outputFilename);
+    mergeRun.FormOutputFile();
+    mergeRun.GetOutputFile()->cd();
+    mergeRun.SetRunType("restG4");
 
     TRestGeant4Event* mergeEvent = nullptr;
-    auto mergeEventTree = mergeRun->GetEventTree();
+    auto mergeEventTree = mergeRun.GetEventTree();
     mergeEventTree->Branch("TRestGeant4EventBranch", "TRestGeant4Event", &mergeEvent);
 
     set<Int_t> eventIds;  // std::set is sorted from lower to higher automatically
@@ -71,8 +71,8 @@ void REST_Geant4_MergeRestG4Files(const char* outputFilename, const char* inputF
         map<Int_t, Int_t>
             eventIdUpdates;  // repeatedId -> newId. Make sure if there are repeated event ids in a file
                              // (because of sub-events) they keep the same event id after modification
-        auto run = TRestRun(inputFiles[i].c_str());
-        auto metadata = (TRestGeant4Metadata*)run.GetMetadataClass("TRestGeant4Metadata");
+        TRestRun run(inputFiles[i].c_str());
+        auto metadata = dynamic_cast<TRestGeant4Metadata*>(run.GetMetadataClass("TRestGeant4Metadata"));
         if (i == 0) {
             mergeMetadata = *metadata;
         } else {
@@ -106,23 +106,23 @@ void REST_Geant4_MergeRestG4Files(const char* outputFilename, const char* inputF
             eventIds.insert(mergeEvent->GetID());
 
             mergeEventTree->Fill();
-            mergeRun->GetAnalysisTree()->Fill();
+            mergeRun.GetAnalysisTree()->Fill();
 
             eventCounter++;
         }
     }
 
-    cout << "Output filename: " << mergeRun->GetOutputFileName() << endl;
-    cout << "Output file: " << mergeRun->GetOutputFile() << endl;
+    cout << "Output filename: " << mergeRun.GetOutputFileName() << endl;
+    cout << "Output file: " << mergeRun.GetOutputFile() << endl;
 
-    mergeRun->GetOutputFile()->cd();
+    mergeRun.GetOutputFile()->cd();
 
     gGeoManager->Write("Geometry", TObject::kOverwrite);
 
     mergeMetadata.SetName("geant4Metadata");
     mergeMetadata.Write();
-    mergeRun->UpdateOutputFile();
-    mergeRun->CloseFile();
+    mergeRun.UpdateOutputFile();
+    mergeRun.CloseFile();
 
     // Open the file again to check the number of events
     TRestRun runCheck(outputFilename);

src/TRestGeant4Event.cxx

@@ -181,9 +181,11 @@ TVector3 TRestGeant4Event::GetFirstPositionInVolume(Int_t volID) const {
 /// \param volID Int_t specifying volume ID
 ///
 TVector3 TRestGeant4Event::GetLastPositionInVolume(Int_t volID) const {
-    for (int t = GetNumberOfTracks() - 1; t >= 0; t--)
-        if (GetTrack(t).GetEnergyInVolume(volID) > 0) return GetTrack(t).GetLastPositionInVolume(volID);
-
+    for (int t = GetNumberOfTracks() - 1; t >= 0; t--) {
+        if (GetTrack(t).GetEnergyInVolume(volID) > 0) {
+            return GetTrack(t).GetLastPositionInVolume(volID);
+        }
+    }
     Double_t nan = TMath::QuietNaN();
     return {nan, nan, nan};
 }

src/TRestGeant4Metadata.cxx

@@ -856,6 +856,9 @@ void TRestGeant4Metadata::InitFromConfigFile() {
     }
     fNEvents = nEventsString == PARAMETER_NOT_FOUND_STR ? 0 : StringToInteger(nEventsString);
 
+    fStoreTracks = ToUpper(GetParameter("storeTracks", "true")) == "TRUE" ||
+                   ToUpper(GetParameter("storeTracks", "on")) == "ON";
+
     const auto fNRequestedEntriesString = GetParameter("nRequestedEntries");
     fNRequestedEntries =
         fNRequestedEntriesString == PARAMETER_NOT_FOUND_STR ? 0 : StringToInteger(fNRequestedEntriesString);
@@ -921,44 +924,59 @@ void TRestGeant4Metadata::InitFromConfigFile() {
 ///
 
 Double_t TRestGeant4Metadata::GetCosmicFluxInCountsPerCm2PerSecond() const {
-    if (TRestGeant4PrimaryGeneratorTypes::StringToSpatialGeneratorTypes(
-            fGeant4PrimaryGeneratorInfo.GetSpatialGeneratorType().Data()) !=
-        TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::COSMIC) {
-        RESTError
-            << "TRestGeant4Metadata::GetEquivalentSimulatedTime can only be called for 'cosmic' generator"
-            << RESTendl;
-        exit(1);
-    }
     const auto source = GetParticleSource();
+
+    double countsPerSecondPerCm2 = 0;
     if (TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionTypes(
-            source->GetEnergyDistributionType().Data()) !=
-        TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypes::FORMULA2) {
-        RESTError << "TRestGeant4Metadata::GetEquivalentSimulatedTime can only be called for 'formula2' "
-                     "energy distribution"
-                  << RESTendl;
-        exit(1);
+            source->GetEnergyDistributionType().Data()) ==
+            TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypes::FORMULA2 &&
+        TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes(
+            source->GetAngularDistributionType().Data()) ==
+            TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypes::FORMULA2) {
+        const auto energyRange = source->GetEnergyDistributionRange();
+        const auto angularRange = source->GetAngularDistributionRange();
+        auto function = (TF2*)source->GetEnergyAndAngularDistributionFunction()->Clone();
+        // counts per second per cm2 (distribution is already integrated over uniform phi)
+        countsPerSecondPerCm2 =
+            function->Integral(energyRange.X(), energyRange.Y(), angularRange.X(), angularRange.Y(), 1E-9);
     }
-    if (TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes(
-            source->GetAngularDistributionType().Data()) !=
-        TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypes::FORMULA2) {
-        RESTError << "TRestGeant4Metadata::GetEquivalentSimulatedTime can only be called for 'formula' "
-                     "angular distribution"
-                  << RESTendl;
-        exit(1);
+
+    else if (std::string(source->GetName()) == "TRestGeant4ParticleSourceCosmics") {
+        auto cosmicSource = dynamic_cast<TRestGeant4ParticleSourceCosmics*>(source);
+        const auto names = cosmicSource->GetParticleNames();
+        const auto histograms = cosmicSource->GetHistogramsTransformed();
+        for (const auto& name : names) {
+            countsPerSecondPerCm2 += histograms.at(name)->Integral();
+        }
+    } else if (TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionTypes(
+                   source->GetEnergyDistributionType().Data()) ==
+                   TRestGeant4PrimaryGeneratorTypes::EnergyDistributionTypes::TH2D &&
+               TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes(
+                   source->GetAngularDistributionType().Data()) ==
+                   TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypes::TH2D) {
+        const auto filename = source->GetEnergyDistributionFilename();
+        const auto name = source->GetEnergyDistributionNameInFile();
+        TFile* file = TFile::Open(filename);
+        auto energyAndAngularDistributionHistogram =
+            file->Get<TH2D>(name);  // it's the same for both angular and energy
+        if (energyAndAngularDistributionHistogram == nullptr) {
+            throw std::runtime_error("Histogram not found in file");
+        }
+        countsPerSecondPerCm2 = energyAndAngularDistributionHistogram->Integral();
+        file->Close();
+        delete file;
+    }
+
+    else {
+        throw std::runtime_error("Cosmic flux calculation is only supported for TFormula2 or TH2D sources");
     }
 
-    const auto energyRange = source->GetEnergyDistributionRange();
-    const auto angularRange = source->GetAngularDistributionRange();
-    auto function = (TF2*)source->GetEnergyAndAngularDistributionFunction()->Clone();
-    // counts per second per cm2 (distribution is already integrated over uniform phi)
-    const auto countsPerSecondPerCm2 =
-        function->Integral(energyRange.X(), energyRange.Y(), angularRange.X(), angularRange.Y(), 1E-9);
     return countsPerSecondPerCm2;
 }
 
 Double_t TRestGeant4Metadata::GetCosmicIntensityInCountsPerSecond() const {
     const auto countsPerSecondPerCm2 = GetCosmicFluxInCountsPerCm2PerSecond();
-    const auto surface = fGeant4PrimaryGeneratorInfo.GetSpatialGeneratorCosmicSurfaceTermCm2();
+    const auto surface = GetGeneratorSurfaceCm2();
     const auto countsPerSecond = countsPerSecondPerCm2 * surface;
     return countsPerSecond;
 }
@@ -1565,13 +1583,32 @@ void TRestGeant4Metadata::SetActiveVolume(const TString& name, Double_t chance,
     fActiveVolumesSet.insert(name.Data());
 }
 
+double TRestGeant4Metadata::GetGeneratorSurfaceCm2() const {
+    const auto type = ToLower(fGeant4PrimaryGeneratorInfo.GetSpatialGeneratorType());
+    const auto shape = ToLower(fGeant4PrimaryGeneratorInfo.GetSpatialGeneratorShape());
+
+    if (type == "surface" && shape == "circle") {
+        const auto radius = fGeant4PrimaryGeneratorInfo.GetSpatialGeneratorSize().X();
+        return TMath::Pi() * radius * radius * 0.01;  // cm2
+    } else if (type == "cosmic") {
+        return fGeant4PrimaryGeneratorInfo.GetSpatialGeneratorCosmicSurfaceTermCm2();
+    }
+
+    throw std::runtime_error(
+        "TRestGeant4Metadata::GetGeneratorSurfaceCm2: Can only be called for 'cosmic' and 'surface/circle' "
+        "generators");
+}
+
 ///////////////////////////////////////////////
 /// \brief Returns true if the volume named *volName* has been registered for
 /// data storage.
 ///
 Bool_t TRestGeant4Metadata::isVolumeStored(const TString& volume) const {
-    for (unsigned int n = 0; n < GetNumberOfActiveVolumes(); n++)
-        if (GetActiveVolumeName(n) == volume) return true;
+    for (unsigned int n = 0; n < GetNumberOfActiveVolumes(); n++) {
+        if (GetActiveVolumeName(n) == volume) {
+            return true;
+        }
+    }
 
     return false;
 }
@@ -1580,9 +1617,12 @@ Bool_t TRestGeant4Metadata::isVolumeStored(const TString& volume) const {
 /// \brief Returns the probability of an active volume being stored
 ///
 Double_t TRestGeant4Metadata::GetStorageChance(const TString& volume) {
-    Int_t id;
-    for (id = 0; id < (Int_t)fActiveVolumes.size(); id++) {
-        if (fActiveVolumes[id] == volume) return fChance[id];
+    for (auto id = 0; id < (Int_t)fActiveVolumes.size(); id++) {
+        {
+            if (fActiveVolumes[id] == volume) {
+                return fChance[id];
+            }
+        }
     }
     RESTWarning << "TRestGeant4Metadata::GetStorageChance. Volume " << volume << " not found" << RESTendl;
 
@@ -1624,7 +1664,9 @@ void TRestGeant4Metadata::Merge(const TRestGeant4Metadata& metadata) {
     fSimulationTime += metadata.fSimulationTime;
 }
 
-TRestGeant4Metadata::TRestGeant4Metadata(const TRestGeant4Metadata& metadata) { *this = metadata; }
+TRestGeant4Metadata::TRestGeant4Metadata(const TRestGeant4Metadata& metadata) : TRestMetadata(metadata) {
+    *this = metadata;
+}
 
 TRestGeant4Metadata& TRestGeant4Metadata::operator=(const TRestGeant4Metadata& metadata) {
     fIsMerge = metadata.fIsMerge;
@@ -1647,6 +1689,7 @@ TRestGeant4Metadata& TRestGeant4Metadata::operator=(const TRestGeant4Metadata& m
     fSensitiveVolumes = metadata.fSensitiveVolumes;
     fNEvents = metadata.fNEvents;
     fNRequestedEntries = metadata.fNRequestedEntries;
+    fStoreTracks = metadata.fStoreTracks;
     fSimulationMaxTimeSeconds = metadata.fSimulationMaxTimeSeconds;
     fSeed = metadata.fSeed;
     fSaveAllEvents = metadata.fSaveAllEvents;