Implement modular biasing options

inc/TRestGeant4BiasingInfo.h

@@ -0,0 +1,48 @@
+
+#ifndef REST_TRESTGEANT4BIASINGINFO_H
+#define REST_TRESTGEANT4BIASINGINFO_H
+
+#include <TVector3.h>
+
+#include <set>
+#include <string>
+
+class TRestGeant4BiasingInfo {
+    ClassDef(TRestGeant4BiasingInfo, 1);
+
+   private:
+    bool fEnabled = false;
+    int fSplittingFactor = 1;
+    bool fBiasOncePerTrack = true;
+    std::set<std::string> fBiasingVolumes;
+    TVector3 fBiasingCenter;
+
+   public:
+    inline TRestGeant4BiasingInfo() = default;
+
+    inline ~TRestGeant4BiasingInfo() = default;
+
+    inline void SetEnabled(bool enabled) { fEnabled = enabled; }
+
+    inline bool GetEnabled() const { return fEnabled; }
+
+    inline void SetSplittingFactor(unsigned int splittingFactor) { fSplittingFactor = splittingFactor; }
+
+    inline int GetSplittingFactor() const { return fSplittingFactor; }
+
+    inline void SetBiasOncePerTrack(bool biasOncePerTrack) { fBiasOncePerTrack = biasOncePerTrack; }
+
+    inline bool GetBiasOncePerTrack() const { return fBiasOncePerTrack; }
+
+    inline void AddBiasingVolume(const std::string& volumeName) { fBiasingVolumes.insert(volumeName); }
+
+    inline void ClearBiasingVolumes() { fBiasingVolumes.clear(); }
+
+    inline std::set<std::string> GetBiasingVolumes() const { return fBiasingVolumes; }
+
+    inline void SetBiasingCenter(const TVector3& biasingCenter) { fBiasingCenter = biasingCenter; }
+
+    inline TVector3 GetBiasingCenter() const { return fBiasingCenter; }
+};
+
+#endif  // REST_TRESTGEANT4BIASINGINFO_H

inc/TRestGeant4GeometryInfo.h

@@ -50,13 +50,17 @@ class TRestGeant4GeometryInfo {
     void PopulateFromGdml(const TString&);
 
     TString GetAlternativeNameFromGeant4PhysicalName(const TString&) const;
+
     TString GetGeant4PhysicalNameFromAlternativeName(const TString&) const;
 
     std::vector<TString> GetAllLogicalVolumes() const;
+
     std::vector<TString> GetAllPhysicalVolumes() const;
+
     std::vector<TString> GetAllAlternativePhysicalVolumes() const;
 
     std::vector<TString> GetAllLogicalVolumesMatchingExpression(const TString&) const;
+
     std::vector<TString> GetAllPhysicalVolumesMatchingExpression(const TString&) const;
 
     inline bool IsValidGdmlName(const TString& volume) const {
@@ -71,16 +75,25 @@ class TRestGeant4GeometryInfo {
     inline bool IsValidPhysicalVolume(const TString& volume) const {
         return fPhysicalToLogicalVolumeMap.count(volume) > 0;
     }
+
     inline bool IsValidLogicalVolume(const TString& volume) const {
         return fLogicalToPhysicalMap.count(volume) > 0;
     }
+
     inline std::vector<TString> GetAllPhysicalVolumesFromLogical(const TString& logicalVolume) const {
         if (IsValidLogicalVolume(logicalVolume)) {
             return fLogicalToPhysicalMap.at(logicalVolume);
         }
         return {};
     }
 
+    inline TString GetLogicalVolumeFromPhysical(const TString& physicalVolume) const {
+        if (IsValidPhysicalVolume(physicalVolume)) {
+            return fPhysicalToLogicalVolumeMap.at(physicalVolume);
+        }
+        return {};
+    }
+
     inline TVector3 GetPosition(const TString& volume) const {
         return fPhysicalToPositionInWorldMap.at(volume);
     }
@@ -90,6 +103,7 @@ class TRestGeant4GeometryInfo {
     void InsertVolumeName(Int_t id, const TString& volumeName);
 
     TString GetVolumeFromID(Int_t id) const;
+
     Int_t GetIDFromVolume(const TString& volumeName) const;
 
     void Print() const;

inc/TRestGeant4Metadata.h

@@ -37,6 +37,7 @@
 #include <iostream>
 #include <vector>
 
+#include "TRestGeant4BiasingInfo.h"
 #include "TRestGeant4BiasingVolume.h"
 #include "TRestGeant4GeometryInfo.h"
 #include "TRestGeant4ParticleSource.h"
@@ -51,9 +52,11 @@ class TRestGeant4Metadata : public TRestMetadata {
     void InitFromConfigFile() override;
 
     void ReadGenerator();
+
     void ReadParticleSource(TRestGeant4ParticleSource* src, TiXmlElement* sourceDefinition);
 
     void ReadDetector();
+
     void ReadBiasing();
 
     bool fDetectorSectionInitialized = false;  //!
@@ -64,6 +67,9 @@ class TRestGeant4Metadata : public TRestMetadata {
     /// Class used to store and retrieve physics info such as process names or particle names
     TRestGeant4PhysicsInfo fGeant4PhysicsInfo;
 
+    /// Class used to store and retrieve biasing info
+    TRestGeant4BiasingInfo fGeant4BiasingInfo;
+
     /// Class used to store and retrieve Geant4 primary generator info
     TRestGeant4PrimaryGeneratorInfo fGeant4PrimaryGeneratorInfo;
 
@@ -152,6 +158,8 @@ class TRestGeant4Metadata : public TRestMetadata {
 
     std::set<std::string> fKillVolumes;
 
+    bool fKillVolumesComputeEnergy = false;
+
     /// If this parameter is set to 'true' it will print out on screen every time 10k events are reached.
     Bool_t fPrintProgress = false;  //!
 
@@ -175,6 +183,9 @@ 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 biasing info
+    inline const TRestGeant4BiasingInfo& GetGeant4BiasingInfo() const { return fGeant4BiasingInfo; }
+
     /// \brief Returns an immutable reference to the primary generator info
     inline const TRestGeant4PrimaryGeneratorInfo& GetGeant4PrimaryGeneratorInfo() const {
         return fGeant4PrimaryGeneratorInfo;
@@ -332,6 +343,10 @@ class TRestGeant4Metadata : public TRestMetadata {
 
     inline bool IsKillVolume(const char* volumeName) const { return fKillVolumes.count(volumeName) > 0; }
 
+    inline bool GetKillVolumesComputeEnergy() const { return fKillVolumesComputeEnergy; }
+
+    void SetKillVolumesComputeEnergy(bool computeEnergy) { fKillVolumesComputeEnergy = computeEnergy; }
+
     inline std::vector<std::string> GetKillVolumes() const {
         std::vector<std::string> result;
         for (const auto& volume : fKillVolumes) {
@@ -349,7 +364,9 @@ class TRestGeant4Metadata : public TRestMetadata {
     }
 
     Double_t GetCosmicFluxInCountsPerCm2PerSecond() const;
+
     Double_t GetCosmicIntensityInCountsPerSecond() const;
+
     Double_t GetEquivalentSimulatedTime() const;
 
     /// Returns a std::string with the name of the active volume with index n
@@ -375,15 +392,19 @@ class TRestGeant4Metadata : public TRestMetadata {
     void PrintMetadata() override;
 
     TRestGeant4Metadata();
+
     TRestGeant4Metadata(const char* configFilename, const std::string& name = "");
 
     ~TRestGeant4Metadata();
 
-    ClassDefOverride(TRestGeant4Metadata, 13);
+    ClassDefOverride(TRestGeant4Metadata, 14);
 
     // 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

inc/TRestGeant4PrimaryGeneratorInfo.h

@@ -48,7 +48,7 @@ EnergyDistributionTypes StringToEnergyDistributionTypes(const std::string&);
 
 enum class EnergyDistributionFormulas {
     COSMIC_NEUTRONS,
-    COSMIC_GAMMAS,
+    COSMIC,
 };
 
 std::string EnergyDistributionFormulasToString(const EnergyDistributionFormulas&);

src/TRestGeant4BiasingInfo.cxx

@@ -0,0 +1,2 @@
+
+#include "TRestGeant4BiasingInfo.h"

src/TRestGeant4Metadata.cxx

@@ -624,6 +624,8 @@
 /// to stop a track after entering the volume. The track will be immediately killed and its energy deposition
 /// won't be computed. This is useful to speed up certain kinds of simulations.
 ///
+/// If "killVolumesComputeEnergy" option is enabled,
+/// the energy deposition in the volumes marked as kill will be computed towards the sensitive volume.
 /// ## 4. The biasing volumes section (optional)
 ///
 /// The REST Geant4 toolkit (*restG4*) implements a particular biasing
@@ -846,13 +848,16 @@ void TRestGeant4Metadata::InitFromConfigFile() {
     fPrintProgress = ToUpper(GetParameter("printProgress", "true")) == "TRUE" ||
                      ToUpper(GetParameter("printProgress", "on")) == "ON";
 
+    fKillVolumesComputeEnergy = ToUpper(GetParameter("killVolumesComputeEnergy", "true")) == "TRUE" ||
+                                ToUpper(GetParameter("killVolumesComputeEnergy", "on")) == "ON";
+
     fRegisterEmptyTracks = ToUpper(GetParameter("registerEmptyTracks", "false")) == "TRUE" ||
                            ToUpper(GetParameter("registerEmptyTracks", "off")) == "ON";
 
+    ReadBiasing();
     ReadGenerator();
     // 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);
     if (fMaxTargetStepSize > 0) {
@@ -949,20 +954,62 @@ Double_t TRestGeant4Metadata::GetEquivalentSimulatedTime() const {
     return seconds;
 }
 
+/*
+ * OLD IMPLEMENTATION
+ *
+void TRestGeant4Metadata::ReadBiasing() {
+
+}
+*/
+
 void TRestGeant4Metadata::ReadBiasing() {
     TiXmlElement* biasingDefinition = GetElement("biasing");
     if (biasingDefinition == nullptr) {
         fNBiasingVolumes = 0;
+        fGeant4BiasingInfo.SetEnabled(false);
         return;
     }
-    TString biasEnabled = GetFieldValue("value", biasingDefinition);
-    TString biasType = GetFieldValue("type", biasingDefinition);
+    string biasEnabled = GetFieldValue("value", biasingDefinition);
+    // biasEnabled to lowercase
+    std::transform(biasEnabled.begin(), biasEnabled.end(), biasEnabled.begin(), ::tolower);
+    if (biasEnabled == "off") {
+        fNBiasingVolumes = 0;
+        fGeant4BiasingInfo.SetEnabled(false);
+        return;
+    }
+
+    RESTInfo << "Biasing is enabled" << RESTendl;
 
-    RESTDebug << "Bias: " << biasEnabled << RESTendl;
+    string biasType = GetFieldValue("type", biasingDefinition);
 
-    if (biasEnabled == "on" || biasEnabled == "ON" || biasEnabled == "On" || biasEnabled == "oN") {
-        RESTInfo << "Biasing is enabled" << RESTendl;
+    if (biasType == "split") {
+        fGeant4BiasingInfo.SetEnabled(true);
 
+        const int splittingFactor = StringToInteger(GetParameter("splittingFactor", biasingDefinition, "1"));
+
+        if (splittingFactor <= 1) {
+            RESTError << "Parameter 'splittingFactor' in biasing section must be set to an integer > 1"
+                      << RESTendl;
+            exit(1);
+        }
+
+        fGeant4BiasingInfo.SetSplittingFactor(splittingFactor);
+
+        const TVector3 center = StringTo3DVector(GetParameter("center", biasingDefinition, "(0,0,0)"));
+        fGeant4BiasingInfo.SetBiasingCenter(center);
+
+        TiXmlElement* volumeDefinition = GetElement("volume", biasingDefinition);
+        while (volumeDefinition != nullptr) {
+            string name = GetFieldValue("name", volumeDefinition);
+            if (name.empty() || name == "Not defined") {
+                RESTError << "volume inside biasing section defined without name" << RESTendl;
+                exit(1);
+            }
+            fGeant4BiasingInfo.AddBiasingVolume(name);
+            volumeDefinition = GetNextElement(volumeDefinition);
+        }
+    } else {
+        // default implementation
         TiXmlElement* biasVolumeDefinition = GetElement("biasingVolume", biasingDefinition);
         Int_t n = 0;
         while (biasVolumeDefinition) {
@@ -976,10 +1023,6 @@ void TRestGeant4Metadata::ReadBiasing() {
             biasVolume.SetEnergyRange(Get2DVectorParameterWithUnits("energyRange", biasVolumeDefinition));
             biasVolume.SetBiasingVolumeType(biasType);  // For the moment all the volumes should be same type
 
-            /* TODO check that values are right if not printBiasingVolume with
-            getchar() biasVolume.PrintBiasingVolume(); getchar();
-            */
-
             fBiasingVolumes.push_back(biasVolume);
             biasVolumeDefinition = GetNextElement(biasVolumeDefinition);
             n++;

src/TRestGeant4PrimaryGeneratorInfo.cxx

@@ -163,8 +163,8 @@ string TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToString(
     switch (type) {
         case EnergyDistributionFormulas::COSMIC_NEUTRONS:
             return "CosmicNeutrons";
-        case EnergyDistributionFormulas::COSMIC_GAMMAS:
-            return "CosmicGammas";
+        case EnergyDistributionFormulas::COSMIC:
+            return "Cosmic";
     }
     cout << "TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToString - Error - Unknown energy "
             "distribution formula"
@@ -177,10 +177,9 @@ EnergyDistributionFormulas TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistr
     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 if (TString(type).EqualTo(EnergyDistributionFormulasToString(EnergyDistributionFormulas::COSMIC),
+                                     TString::ECaseCompare::kIgnoreCase)) {
+        return EnergyDistributionFormulas::COSMIC;
     } else {
         cout << "TRestGeant4PrimaryGeneratorTypes::StringToEnergyDistributionFormulas - Error - Unknown "
                 "energyDistributionFormulas: "
@@ -207,8 +206,13 @@ TF1 TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToRootFormula(
             distribution.GetXaxis()->SetTitle("Energy (keV)");
             return distribution;
         }
-        case EnergyDistributionFormulas::COSMIC_GAMMAS:
-            exit(1);
+        case EnergyDistributionFormulas::COSMIC:
+            const char* title = "Cosmic distribution aproximation";
+            auto distribution = TF1(title, "TMath::Power(x, -2.7)", 1E2, 1E9);
+            distribution.SetNormalized(true);  // Normalized, not really necessary
+            distribution.SetTitle(title);
+            distribution.GetXaxis()->SetTitle("Energy (keV)");
+            return distribution;
     }
     cout << "TRestGeant4PrimaryGeneratorTypes::EnergyDistributionFormulasToRootFormula - Error - Unknown "
             "energy distribution formula"