Geant4IsotropeGenerator::getParticleMomentumUniform(double& momentum)

DDG4/include/DDG4/Geant4IsotropeGenerator.h

@@ -40,6 +40,10 @@ namespace dd4hep {
       double      m_thetaMin;
       /// Property: Maximal theta angular value
       double      m_thetaMax;
+      /// Property: Minimal momentum value
+      double      m_momentumMin;
+      /// Property: Maximal momentum value
+      double      m_momentumMax;
       
       /// Particle modification. Caller presets defaults to: ( direction = m_direction, momentum = m_energy)
       /** Use this function to implement isotrop guns, multiple guns etc. 
@@ -54,6 +58,8 @@ namespace dd4hep {
       void getParticleDirectionCosTheta(int num, ROOT::Math::XYZVector& direction, double& momentum) const;
       /// Uniform particle distribution
       void getParticleDirectionUniform(int num, ROOT::Math::XYZVector& direction, double& momentum) const;
+      /// Uniform particle momentum
+      void getParticleMomentumUniform(double& momentum) const;
 
     public:
       /// Inhibit default constructor

DDG4/src/Geant4IsotropeGenerator.cpp

@@ -29,6 +29,8 @@ Geant4IsotropeGenerator::Geant4IsotropeGenerator(Geant4Context* ctxt, const stri
   declareProperty("PhiMax", m_phiMax = 2.0*M_PI);
   declareProperty("ThetaMin", m_thetaMin = 0.0);
   declareProperty("ThetaMax", m_thetaMax = M_PI);
+  declareProperty("MomentumMin", m_momentumMin =  0.0);
+  declareProperty("MomentumMax", m_momentumMax = -1.0);
   declareProperty("Distribution", m_distribution = "uniform" );
 }
 
@@ -37,6 +39,16 @@ Geant4IsotropeGenerator::~Geant4IsotropeGenerator() {
   InstanceCount::decrement(this);
 }
 
+/// Uniform momentum distribution
+void Geant4IsotropeGenerator::getParticleMomentumUniform(double& momentum) const  {
+  Geant4Event&  evt = context()->event();
+  Geant4Random& rnd = evt.random();
+  if (m_momentumMax < m_momentumMin)
+    momentum = m_momentumMin+(momentum-m_momentumMin)*rnd.rndm();
+  else
+    momentum = m_momentumMin+(m_momentumMax-m_momentumMin)*rnd.rndm();
+}
+
 /// Uniform particle distribution
 void Geant4IsotropeGenerator::getParticleDirectionUniform(int, ROOT::Math::XYZVector& direction, double& momentum) const   {
   Geant4Event&  evt = context()->event();
@@ -48,7 +60,7 @@ void Geant4IsotropeGenerator::getParticleDirectionUniform(int, ROOT::Math::XYZVe
   double x3 = std::cos(theta);
 
   direction.SetXYZ(x1,x2,x3);
-  momentum = rnd.rndm()*momentum;
+  getParticleMomentumUniform(momentum);
 }
 
 /// Particle distribution ~ cos(theta)
@@ -63,7 +75,10 @@ void Geant4IsotropeGenerator::getParticleDirectionCosTheta(int, ROOT::Math::XYZV
   double x3 = cos_theta;
 
   direction.SetXYZ(x1,x2,x3);
-  momentum = rnd.rndm()*momentum;
+  if (m_momentumMax < m_momentumMin)
+    momentum = m_momentumMin+(momentum-m_momentumMin)*rnd.rndm();
+  else
+    momentum = m_momentumMin+(m_momentumMax-m_momentumMin)*rnd.rndm();
 }
 
 /// Particle distribution flat in eta (pseudo rapidity)
@@ -84,7 +99,10 @@ void Geant4IsotropeGenerator::getParticleDirectionEta(int, ROOT::Math::XYZVector
   double x3         = std::sinh(eta);
   double r          = std::sqrt(1.0+x3*x3);
   direction.SetXYZ(x1/r,x2/r,x3/r);
-  momentum = rnd.rndm()*momentum;
+  if (m_momentumMax < m_momentumMin)
+    momentum = m_momentumMin+(momentum-m_momentumMin)*rnd.rndm();
+  else
+    momentum = m_momentumMin+(m_momentumMax-m_momentumMin)*rnd.rndm();
 }
 
 /// e+e- --> ffbar particle distribution ~ 1 + cos^2(theta)
@@ -110,6 +128,10 @@ void Geant4IsotropeGenerator::getParticleDirectionFFbar(int, ROOT::Math::XYZVect
       double x2 = std::sin(theta)*std::sin(phi);
       double x3 = std::cos(theta);
       direction.SetXYZ(x1,x2,x3);
+      if (m_momentumMax < m_momentumMin)
+        momentum = m_momentumMin+(momentum-m_momentumMin)*rnd.rndm();
+      else
+        momentum = m_momentumMin+(m_momentumMax-m_momentumMin)*rnd.rndm();
       momentum = rnd.rndm()*momentum;
       return;
     }