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   declareProperty("Halton", m_halton = false);
0033   declareProperty("HaltonOffset", m_haltonOffset = 0ULL);
0034   m_haltonIndex = 0;
0035   m_haltonShift = {0.0, 0.0, 0.0};
0036 }
0037
0038 /// Default destructor
0039 Geant4IsotropeGenerator::~Geant4IsotropeGenerator() {
0040   InstanceCount::decrement(this);
0041 }
0042
0043 /// Initialize Cranley-Patterson shifts using rnd; set m_haltonIndex = m_haltonOffset
0044 void Geant4IsotropeGenerator::initHalton(Geant4Random& rnd) const  {
0045   for (auto& s : m_haltonShift) s = rnd.rndm();
0046   m_haltonIndex = m_haltonOffset;
0047 }
0048
0049 /// Compute the radical inverse of index in the given base (Halton value in [0,1))
0050 double Geant4IsotropeGenerator::haltonValue(uint64_t index, int base)  {
0051   uint64_t num   = 0;
0052   uint64_t denom = 1;
0053   while ( index > 0 )  {
0054     denom *= base;
0055     num    = num * base + (index % base);
0056     index /= base;
0057   }
0058   return static_cast<double>(num) / static_cast<double>(denom);
0059 }
0060
0061 /// Cranley-Patterson scrambled Halton sample for dimension dim (0=phi, 1=theta, 2=momentum)
0062 double Geant4IsotropeGenerator::haltonScrambled(uint64_t index, unsigned int dim) const  {
0063   static constexpr std::array<int,3> bases = {2, 3, 5};
0064   if ( dim >= bases.size() )  {
0065     except("haltonScrambled: dimension %u out of range [0,%zu].", dim, bases.size()-1);
0066   }
0067   return std::fmod(haltonValue(index, bases[dim]) + m_haltonShift[dim], 1.0);
0068 }
0069
0070 /// Build a per-particle sampler: PRNG or Halton depending on m_halton
0071 std::function<double(unsigned int)> Geant4IsotropeGenerator::makeSampler(Geant4Random& rnd) const  {
0072   if ( m_halton )  {
0073     std::call_once(m_haltonOnce, &Geant4IsotropeGenerator::initHalton, this, std::ref(rnd));
0074     uint64_t idx = m_haltonIndex++;
0075     return [this, idx](unsigned int dim) { return haltonScrambled(idx, dim); };
0076   }
0077   return [&rnd](unsigned int) { return rnd.rndm(); };
0078 }
0079
0080 /// Sample momentum from [m_momentumMin, m_momentumMax] using a pre-computed [0,1) value h
0081 void Geant4IsotropeGenerator::sampleMomentum(double h, double& momentum) const  {
0082   if ( m_energy != -1 )  {
0083     momentum = m_energy;
0084     return;
0085   }
0086   if ( m_momentumMax < m_momentumMin )
0087     momentum = m_momentumMin + (momentum - m_momentumMin) * h;
0088   else
0089     momentum = m_momentumMin + (m_momentumMax - m_momentumMin) * h;
0090 }
0091
0092 /// Uniform particle distribution
0093 void Geant4IsotropeGenerator::getParticleDirectionUniform(int, const std::function<double(unsigned int)>& sample,
0094                                                           ROOT::Math::XYZVector& direction, double& momentum) const  {
0095   double phi   = m_phiMin   + (m_phiMax   - m_phiMin)   * sample(0);
0096   double theta = m_thetaMin + (m_thetaMax - m_thetaMin) * sample(1);
0097   direction.SetXYZ(std::sin(theta)*std::cos(phi), std::sin(theta)*std::sin(phi), std::cos(theta));
0098   sampleMomentum(sample(2), momentum);
0099 }
0100
0101 /// Particle distribution ~ cos(theta)
0102 void Geant4IsotropeGenerator::getParticleDirectionCosTheta(int, const std::function<double(unsigned int)>& sample,
0103                                                            ROOT::Math::XYZVector& direction, double& momentum) const  {
0104   double phi       = m_phiMin + (m_phiMax - m_phiMin) * sample(0);
0105   double cos_theta = std::cos(m_thetaMin) + (std::cos(m_thetaMax) - std::cos(m_thetaMin)) * sample(1);
0106   double sin_theta = std::sqrt(1.0 - cos_theta*cos_theta);
0107   direction.SetXYZ(sin_theta*std::cos(phi), sin_theta*std::sin(phi), cos_theta);
0108   sampleMomentum(sample(2), momentum);
0109 }
0110
0111 /// Particle distribution flat in eta (pseudo rapidity)
0112 void Geant4IsotropeGenerator::getParticleDirectionEta(int, const std::function<double(unsigned int)>& sample,
0113                                                       ROOT::Math::XYZVector& direction, double& momentum) const  {
0114   struct Distribution {
0115     static double eta(double x)  { return -1.0*std::log(std::tan(x/2.0)); }
0116   };
0117   const double dmin = std::numeric_limits<double>::epsilon();
0118   double phi        = m_phiMin + (m_phiMax - m_phiMin) * sample(0);
0119   double eta_max    = Distribution::eta(m_thetaMin > dmin ? m_thetaMin : dmin);
0120   double eta_min    = Distribution::eta(m_thetaMax > (M_PI - dmin) ? M_PI - dmin : m_thetaMax);
0121   double eta        = eta_min + (eta_max - eta_min) * sample(1);
0122   double x1         = std::cos(phi);
0123   double x2         = std::sin(phi);
0124   double x3         = std::sinh(eta);
0125   double r          = std::sqrt(1.0+x3*x3);
0126   direction.SetXYZ(x1/r,x2/r,x3/r);
0127   sampleMomentum(sample(2), momentum);
0128 }
0129
0130 /// e+e- --> ffbar particle distribution ~ 1 + cos^2(theta) (PRNG only; incompatible with Halton)
0131 void Geant4IsotropeGenerator::getParticleDirectionFFbar(int, ROOT::Math::XYZVector& direction, double& momentum) const   {
0132   struct Distribution {
0133     static double ffbar(double x)  { double c = std::cos(x); return 1 + c*c; }
0134     //static double integral(double x)  { return 1.5*x + sin(2.*x)/4.0; }
0135   };
0136   Geant4Event&  evt = context()->event();
0137   Geant4Random& rnd = evt.random();
0138   double phi = m_phiMin+(m_phiMax-m_phiMin)*rnd.rndm();
0139
0140   // 1 + cos^2(theta) cannot be integrated and then inverted.
0141   // We have to throw the dice until it fits.....
0142   double v1 = Distribution::ffbar(m_thetaMin), v2 = Distribution::ffbar(m_thetaMax);
0143   double vmax = std::max(v1,v2); // ffbar symmetric and monotonic in |x|.
0144   while(1)  {
0145     double dice  = rnd.rndm()*vmax;
0146     double theta = m_thetaMin+(m_thetaMax-m_thetaMin)*rnd.rndm();
0147     if ( dice <= Distribution::ffbar(theta) )  {
0148       direction.SetXYZ(std::sin(theta)*std::cos(phi), std::sin(theta)*std::sin(phi), std::cos(theta));
0149       sampleMomentum(rnd.rndm(), momentum);
0150       return;
0151     }
0152   }
0153 }
0154
0155 /// Particle modification. Caller presets defaults to: ( direction = m_direction, momentum = [mMin, mMax])
0156 void Geant4IsotropeGenerator::getParticleDirection(int num, ROOT::Math::XYZVector& direction, double& momentum) const   {
0157   if ( m_halton && ::toupper(m_distribution[0]) == 'F' )  {
0158     except("Halton mode is incompatible with the 'ffbar' distribution: "
0159            "acceptance-rejection cannot be driven by a fixed per-particle Halton point. "
0160            "Use distribution='uniform', 'cos(theta)', or 'eta' with Halton mode.");
0161   }
0162   Geant4Event&  evt    = context()->event();
0163   Geant4Random& rnd    = evt.random();
0164   auto          sample = makeSampler(rnd);
0165   switch(::toupper(m_distribution[0]))  {
0166   case 'C':  // cos(theta)
0167     return getParticleDirectionCosTheta(num, sample, direction, momentum);
0168   case 'F':  // ffbar: ~ 1 + cos^2(theta)
0169     return getParticleDirectionFFbar(num, direction, momentum);
0170   case 'E':  // eta, rapidity, pseudorapidity
0171   case 'P':
0172   case 'R':
0173     return getParticleDirectionEta(num, sample, direction, momentum);
0174   case 'U':  // uniform
0175     return getParticleDirectionUniform(num, sample, direction, momentum);
0176   default:
0177     break;
0178   }
0179   except("Unknown distribution density: %s. Cannot generate primaries.",
0180          m_distribution.c_str());
0181 }