From 535efcab734639d3c8232ceb32f596fc1ae3fa29 Mon Sep 17 00:00:00 2001
From: Wouter Deconinck <wdconinc@gmail.com>
Date: Tue, 21 Sep 2021 16:52:10 -0500
Subject: [PATCH] Geant4IsotropeGenerator::getParticleMomentumUniform(double&
momentum)
---
DDG4/include/DDG4/Geant4IsotropeGenerator.h | 6 +++++
DDG4/src/Geant4IsotropeGenerator.cpp | 28 ++++++++++++++++++---
2 files changed, 31 insertions(+), 3 deletions(-)
diff --git a/DDG4/include/DDG4/Geant4IsotropeGenerator.h b/DDG4/include/DDG4/Geant4IsotropeGenerator.h
index 5bf0164b6..084ba36a8 100644
--- a/DDG4/include/DDG4/Geant4IsotropeGenerator.h
+++ b/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
diff --git a/DDG4/src/Geant4IsotropeGenerator.cpp b/DDG4/src/Geant4IsotropeGenerator.cpp
index 54588aafd..24f6176ce 100644
--- a/DDG4/src/Geant4IsotropeGenerator.cpp
+++ b/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;
}
--
GitLab