Merge pull request #103 from rest-for-physics/lobis-isotropic

Add option for angular range for isotropic generator

Author: lobis

inc/TRestGeant4Particle.h

@@ -57,9 +57,9 @@ class TRestGeant4Particle {
 
     void SetParticleCharge(Int_t charge) { fCharge = charge; }
 
-    void SetDirection(TVector3 dir) { fDirection = dir; }
+    void SetDirection(const TVector3& dir) { fDirection = dir.Unit(); }
     void SetEnergy(Double_t en) { fEnergy = en; }
-    void SetOrigin(TVector3 pos) { fOrigin = pos; }
+    void SetOrigin(const TVector3& pos) { fOrigin = pos; }
 
     void Print() const;
 

inc/TRestGeant4ParticleSource.h

@@ -42,6 +42,7 @@ class TRestGeant4ParticleSource : public TRestGeant4Particle, public TRestMetada
     size_t fAngularDistributionFormulaNPoints = 500;
     TF1* fAngularDistributionFunction = nullptr;
     TVector2 fAngularDistributionRange;
+    double fAngularDistributionIsotropicConeHalfAngle = 0;
 
     TString fEnergyDistributionType = "Mono";
     TString fEnergyDistributionFilename;
@@ -65,14 +66,7 @@ class TRestGeant4ParticleSource : public TRestGeant4Particle, public TRestMetada
     virtual void InitFromConfigFile() override;
     static TRestGeant4ParticleSource* instantiate(std::string model = "");
 
-    inline TVector3 GetDirection() const {
-        if (fDirection.Mag() <= 0) {
-            std::cout << "TRestGeant4ParticleSource::GetDirection: "
-                      << "Direction cannot be the zero vector" << std::endl;
-            exit(1);
-        }
-        return fDirection;
-    }
+    TVector3 GetDirection() const;
 
     inline TString GetEnergyDistributionType() const { return fEnergyDistributionType; }
     inline TVector2 GetEnergyDistributionRange() const { return fEnergyDistributionRange; }
@@ -91,6 +85,9 @@ class TRestGeant4ParticleSource : public TRestGeant4Particle, public TRestMetada
     inline TString GetAngularDistributionFilename() const { return fAngularDistributionFilename; }
     inline TString GetAngularDistributionNameInFile() const { return fAngularDistributionNameInFile; }
     inline const TF1* GetAngularDistributionFunction() const { return fAngularDistributionFunction; }
+    inline double GetAngularDistributionIsotropicConeHalfAngle() const {
+        return fAngularDistributionIsotropicConeHalfAngle;
+    }
 
     inline const TF2* GetEnergyAndAngularDistributionFunction() const {
         return fEnergyAndAngularDistributionFunction;
@@ -100,6 +97,16 @@ class TRestGeant4ParticleSource : public TRestGeant4Particle, public TRestMetada
 
     inline std::vector<TRestGeant4Particle> GetParticles() const { return fParticles; }
 
+    inline void SetAngularDistributionIsotropicConeHalfAngle(double angle) {
+        if (angle < 0 || angle > TMath::Pi()) {
+            std::cerr << "TRestGeant4ParticleSource::SetAngularDistributionIsotropicConeHalfAngle: "
+                         "angle must be between 0 and Pi radians"
+                      << std::endl;
+            exit(1);
+        }
+        fAngularDistributionIsotropicConeHalfAngle = angle;
+    }
+
     inline void SetAngularDistributionType(const TString& type) {
         fAngularDistributionType = TRestGeant4PrimaryGeneratorTypes::AngularDistributionTypesToString(
             TRestGeant4PrimaryGeneratorTypes::StringToAngularDistributionTypes(type.Data()));
@@ -201,6 +208,6 @@ class TRestGeant4ParticleSource : public TRestGeant4Particle, public TRestMetada
     // Destructor
     virtual ~TRestGeant4ParticleSource();
 
-    ClassDefOverride(TRestGeant4ParticleSource, 5);
+    ClassDefOverride(TRestGeant4ParticleSource, 6);
 };
 #endif

src/TRestGeant4Metadata.cxx

@@ -461,6 +461,11 @@
 ///     <angular type="isotropic" />
 /// \endcode
 ///
+/// Additionally a direction and a cone angle can be specified to generate the particles in a cone around the
+/// direction. The angle specified is the half opening angle of the cone. \code
+///     <angular type="isotropic" direction="(0,1,0)" coneHalfAngle="30deg"/>
+/// \endcode
+///
 /// * **backtoback**: The source momentum direction will be opposite to
 /// the direction of the previous source. If it is the first source the
 /// angular distribution type will be re-defined to isotropic.
@@ -1146,6 +1151,12 @@ void TRestGeant4Metadata::ReadParticleSource(TRestGeant4ParticleSource* source,
     }
     source->SetDirection(StringTo3DVector(GetParameter("direction", angularDefinition, "(1,0,0)")));
 
+    if ((StringToAngularDistributionTypes(source->GetAngularDistributionType().Data()) ==
+         AngularDistributionTypes::ISOTROPIC)) {
+        source->SetAngularDistributionRange({0, TMath::Pi()});
+        const double coneHalfAngle = GetDblParameterWithUnits("coneHalfAngle", angularDefinition, 0);
+        source->SetAngularDistributionIsotropicConeHalfAngle(coneHalfAngle);
+    }
     // Energy distribution parameters
     TiXmlElement* energyDefinition = GetElement("energy", sourceDefinition);
     if (energyDefinition == nullptr) {

src/TRestGeant4ParticleSource.cxx

@@ -56,6 +56,17 @@ void TRestGeant4ParticleSource::PrintMetadata() {
         }
         RESTMetadata << "Angular distribution direction: (" << GetDirection().X() << "," << GetDirection().Y()
                      << "," << GetDirection().Z() << ")" << RESTendl;
+        if ((StringToAngularDistributionTypes(GetAngularDistributionType().Data()) ==
+             AngularDistributionTypes::ISOTROPIC)) {
+            if (GetAngularDistributionIsotropicConeHalfAngle() != 0) {
+                const double solidAngle =
+                    2 * TMath::Pi() * (1 - TMath::Cos(GetAngularDistributionIsotropicConeHalfAngle()));
+                RESTMetadata << "Angular distribution isotropic cone half angle (deg): "
+                             << GetAngularDistributionIsotropicConeHalfAngle() * TMath::RadToDeg()
+                             << " - solid angle: " << solidAngle << " ("
+                             << solidAngle / (4 * TMath::Pi()) * 100 << "%)" << RESTendl;
+            }
+        }
         if ((StringToAngularDistributionTypes(GetAngularDistributionType().Data()) ==
              AngularDistributionTypes::FORMULA) ||
             StringToAngularDistributionTypes(GetAngularDistributionType().Data()) ==
@@ -146,6 +157,18 @@ void TRestGeant4ParticleSource::Update() {
     }
 }
 
+TVector3 TRestGeant4ParticleSource::GetDirection() const {
+    // direction should be unit (normalized) vector with a tolerance of 0.001
+
+    const double magnitude = fDirection.Mag();
+    constexpr double tolerance = 0.001;
+    if (TMath::Abs(magnitude - 1) > tolerance) {
+        cerr << "TRestGeant4ParticleSource::GetDirection: Direction must be unit vector" << endl;
+        exit(1);
+    }
+
+    return fDirection;
+}
 ///////////////////////////////////////////////
 /// \brief Reads an input file produced by Decay0.
 ///