DDG4/src/Geant4IsotropeGenerator.cpp

0001 //==========================================================================
0002 //  AIDA Detector description implementation
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013
0014 // Framework include files
0015 #include <DD4hep/Printout.h>
0016 #include <DD4hep/InstanceCount.h>
0017 #include <DDG4/Geant4Random.h>
0018 #include <DDG4/Geant4IsotropeGenerator.h>
0019
0020 using namespace dd4hep::sim;
0021
0022 /// Standard constructor
0023 Geant4IsotropeGenerator::Geant4IsotropeGenerator(Geant4Context* ctxt, const std::string& nam)
0024   : Geant4ParticleGenerator(ctxt, nam)
0025 {
0026   InstanceCount::increment(this);
0027   declareProperty("PhiMin", m_phiMin = 0.0);
0028   declareProperty("PhiMax", m_phiMax = 2.0*M_PI);
0029   declareProperty("ThetaMin", m_thetaMin = 0.0);
0030   declareProperty("ThetaMax", m_thetaMax = M_PI);
0031   declareProperty("Distribution", m_distribution = "uniform" );
0032 }
0033
0034 /// Default destructor
0035 Geant4IsotropeGenerator::~Geant4IsotropeGenerator() {
0036   InstanceCount::decrement(this);
0037 }
0038
0039 /// Uniform particle distribution
0040 void Geant4IsotropeGenerator::getParticleDirectionUniform(int, ROOT::Math::XYZVector& direction, double& momentum) const   {
0041   Geant4Event&  evt = context()->event();
0042   Geant4Random& rnd = evt.random();
0043   double phi   = m_phiMin+(m_phiMax-m_phiMin)*rnd.rndm();
0044   double theta = m_thetaMin+(m_thetaMax-m_thetaMin)*rnd.rndm();
0045   double x1 = std::sin(theta)*std::cos(phi);
0046   double x2 = std::sin(theta)*std::sin(phi);
0047   double x3 = std::cos(theta);
0048
0049   direction.SetXYZ(x1,x2,x3);
0050   getParticleMomentumUniform(momentum);
0051 }
0052
0053 /// Particle distribution ~ cos(theta)
0054 void Geant4IsotropeGenerator::getParticleDirectionCosTheta(int, ROOT::Math::XYZVector& direction, double& momentum) const   {
0055   Geant4Event&  evt = context()->event();
0056   Geant4Random& rnd = evt.random();
0057   double phi       = m_phiMin+(m_phiMax-m_phiMin)*rnd.rndm();
0058   double cos_theta = std::cos(m_thetaMin)+(std::cos(m_thetaMax)-std::cos(m_thetaMin))*rnd.rndm();
0059   double sin_theta = std::sqrt(1.0-cos_theta*cos_theta);
0060   double x1 = sin_theta*std::cos(phi);
0061   double x2 = sin_theta*std::sin(phi);
0062   double x3 = cos_theta;
0063
0064   direction.SetXYZ(x1,x2,x3);
0065   getParticleMomentumUniform(momentum);
0066 }
0067
0068 /// Particle distribution flat in eta (pseudo rapidity)
0069 void Geant4IsotropeGenerator::getParticleDirectionEta(int, ROOT::Math::XYZVector& direction, double& momentum) const   {
0070   struct Distribution {
0071     static double eta(double x)  { return -1.0*std::log(std::tan(x/2.0)); }
0072   };
0073   Geant4Event&  evt = context()->event();
0074   Geant4Random& rnd = evt.random();
0075   // See https://en.wikipedia.org/wiki/Pseudorapidity
0076   const double dmin = std::numeric_limits<double>::epsilon();
0077   double phi        = m_phiMin+(m_phiMax-m_phiMin)*rnd.rndm();
0078   double eta_max    = Distribution::eta(m_thetaMin>dmin ? m_thetaMin : dmin);
0079   double eta_min    = Distribution::eta(m_thetaMax>(M_PI-dmin) ? M_PI-dmin : m_thetaMax);
0080   double eta        = eta_min + (eta_max-eta_min)*rnd.rndm();
0081   double x1         = std::cos(phi);
0082   double x2         = std::sin(phi);
0083   double x3         = std::sinh(eta);
0084   double r          = std::sqrt(1.0+x3*x3);
0085   direction.SetXYZ(x1/r,x2/r,x3/r);
0086   getParticleMomentumUniform(momentum);
0087 }
0088
0089 /// e+e- --> ffbar particle distribution ~ 1 + cos^2(theta)
0090 void Geant4IsotropeGenerator::getParticleDirectionFFbar(int, ROOT::Math::XYZVector& direction, double& momentum) const   {
0091   struct Distribution {
0092     static double ffbar(double x)  { double c = std::cos(x); return 1 + c*c; }
0093     //static double integral(double x)  { return 1.5*x + sin(2.*x)/4.0; }
0094   };
0095   Geant4Event&  evt = context()->event();
0096   Geant4Random& rnd = evt.random();
0097   double phi = m_phiMin+(m_phiMax-m_phiMin)*rnd.rndm();
0098
0099   // 1 + cos^2(theta) cannot be integrated and then inverted.
0100   // We have to throw the dice until it fits.....
0101   double v1 = Distribution::ffbar(m_thetaMin), v2 = Distribution::ffbar(m_thetaMax);
0102   double vmax = std::max(v1,v2); // ffbar symmetric and monotonic in |x|.
0103   while(1)  {
0104     double dice  = rnd.rndm()*vmax;
0105     double theta = m_thetaMin+(m_thetaMax-m_thetaMin)*rnd.rndm();
0106     double dist  = Distribution::ffbar(theta);
0107     if ( dice <= dist )   {
0108       double x1 = std::sin(theta)*std::cos(phi);
0109       double x2 = std::sin(theta)*std::sin(phi);
0110       double x3 = std::cos(theta);
0111       direction.SetXYZ(x1,x2,x3);
0112       getParticleMomentumUniform(momentum);
0113       return;
0114     }
0115   }
0116 }
0117
0118 /// Particle modification. Caller presets defaults to: ( direction = m_direction, momentum = [mMin, mMax])
0119 void Geant4IsotropeGenerator::getParticleDirection(int num, ROOT::Math::XYZVector& direction, double& momentum) const   {
0120   switch(::toupper(m_distribution[0]))  {
0121   case 'C':  // cos(theta)
0122     return getParticleDirectionCosTheta(num, direction, momentum);
0123   case 'F':  // ffbar: ~ 1 + cos^2(theta)
0124     return getParticleDirectionFFbar(num, direction, momentum);
0125   case 'E':  // eta, rapidity, pseudorapidity
0126   case 'P':
0127   case 'R':
0128     return getParticleDirectionEta(num, direction, momentum);
0129   case 'U':  // uniform
0130     return getParticleDirectionUniform(num, direction, momentum);
0131   default:
0132     break;
0133   }
0134   except("Unknown distribution densitiy: %s. Cannot generate primaries.",
0135      m_distribution.c_str());
0136 }