DDG4/src/Geant4Random.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
0019 #include <CLHEP/Random/EngineFactory.h>
0020 #include <CLHEP/Random/RandGamma.h>
0021 #include <CLHEP/Random/Random.h>
0022
0023 // ROOT include files
0024 #include <TRandom1.h>
0025
0026 // C/C++ include files
0027 #include <cmath>
0028
0029 using namespace dd4hep::sim;
0030
0031 namespace CLHEP   {
0032   unsigned long crc32ul(const std::string& s);
0033 }
0034
0035 namespace    {
0036
0037   /// Helper class to redirect calls to gRandon to Geant4Random
0038   class RNDM : public TRandom  {
0039     /// Reference to the generator object
0040     Geant4Random* m_generator;
0041     /// Reference to HepRandomEngine
0042     CLHEP::HepRandomEngine* m_engine;
0043
0044   public:
0045     /// Initializing constructor
0046     RNDM(Geant4Random* r) : TRandom(), m_generator(r)    {
0047       m_engine = m_generator->engine();
0048     }
0049     /// Default destructor
0050     virtual ~RNDM() {      }
0051     /// Set new seed
0052     virtual  void SetSeed(UInt_t seed=0)  final  {
0053       fSeed = seed;
0054       m_generator->setSeed((long)seed);
0055     }
0056     /// Set new seed
0057     virtual  void SetSeed(ULong_t seed=0)  final  {
0058       fSeed = seed;
0059       m_generator->setSeed((long)seed);
0060     }
0061     /// Single shot random number creation
0062     virtual Double_t Rndm()  final {
0063       return m_engine->flat();
0064     }
0065     /// Single shot random number creation
0066     virtual Double_t Rndm(Int_t)  final  {
0067       return m_engine->flat();
0068     }
0069     /// Return an array of n random numbers uniformly distributed in ]0,1].
0070     virtual void RndmArray(Int_t size, Float_t *array)  final  {
0071       for (Int_t i=0;i<size;i++) array[i] = m_engine->flat();
0072     }
0073     /// Return an array of n random numbers uniformly distributed in ]0,1].
0074     virtual void RndmArray(Int_t size, Double_t *array)  final  {
0075       m_engine->flatArray(size,array);
0076     }
0077   };
0078   static Geant4Random* s_instance = 0;
0079 }
0080
0081 /// Default constructor
0082 Geant4Random::Geant4Random(Geant4Context* ctxt, const std::string& nam)
0083   : Geant4Action(ctxt,nam), m_engine(0), m_rootRandom(0), m_rootOLD(0),
0084     m_inited(false)
0085 {
0086   declareProperty("File",   m_file="");
0087   declareProperty("Type",   m_engineType="");
0088   declareProperty("Seed",   m_seed = 123456789);
0089   declareProperty("Luxury", m_luxury = 1);
0090   declareProperty("Replace_gRandom",  m_replace = true);
0091   // Default: static Geant4 random engine.
0092   m_engine = CLHEP::HepRandom::getTheEngine();
0093   InstanceCount::increment(this);
0094 }
0095
0096 /// Default destructor
0097 Geant4Random::~Geant4Random()  {
0098   // Only delete the engine if it is NOT the CLEP default one
0099   // BUT: Just cannot delete the engine. Causes havoc with static destructors!
0100   //CLHEP::HepRandomEngine* curr = CLHEP::HepRandom::getTheEngine();
0101   //if ( !m_engineType.empty() && m_engine != curr ) deletePtr(m_engine);
0102
0103   // Set gRandom to the old value
0104   if (  m_rootRandom == gRandom ) gRandom = m_rootOLD;
0105   // Reset instance pointer
0106   if (  s_instance == this ) s_instance = 0;
0107   // Finally delete the TRandom instance wrapper
0108   detail::deletePtr(m_rootRandom);
0109   InstanceCount::decrement(this);
0110 }
0111
0112 /// Access the main Geant4 random generator instance. Must be created before used!
0113 Geant4Random* Geant4Random::instance(bool throw_exception)   {
0114   if ( !s_instance && throw_exception )  {
0115     dd4hep::except("Geant4Random", "No global random number generator defined!");
0116   }
0117   return s_instance;
0118 }
0119
0120 /// Make this random generator instance the one used by Geant4
0121 Geant4Random* Geant4Random::setMainInstance(Geant4Random* ptr)   {
0122   if ( ptr && !ptr->m_inited )  {
0123     ptr->initialize();
0124   }
0125   if ( s_instance != ptr )   {
0126     if ( !ptr )  {
0127       dd4hep::except("Geant4Random","Attempt to declare invalid Geant4Random instance.");
0128     }
0129     if ( !ptr->m_inited )  {
0130       dd4hep::except("Geant4Random","Attempt to declare uninitialized Geant4Random instance.");
0131     }
0132     Geant4Random* old = s_instance;
0133     CLHEP::HepRandomEngine* curr = CLHEP::HepRandom::getTheEngine();
0134     if ( ptr->m_engine != curr )     {
0135       ptr->printP2("Moving CLHEP random instance from %p to %p",curr,ptr->m_engine);
0136       CLHEP::HepRandom::setTheEngine(ptr->m_engine);
0137     }
0138     if ( ptr->m_replace )   {
0139       ptr->m_rootOLD = gRandom;
0140       gRandom = ptr->m_rootRandom;
0141     }
0142     s_instance = ptr;
0143     return old;
0144   }
0145   return 0;
0146 }
0147
0148 #include <CLHEP/Random/DualRand.h>
0149 #include <CLHEP/Random/JamesRandom.h>
0150 #include <CLHEP/Random/MTwistEngine.h>
0151 #include <CLHEP/Random/RanecuEngine.h>
0152 #include <CLHEP/Random/Ranlux64Engine.h>
0153 #include <CLHEP/Random/RanluxEngine.h>
0154 #include <CLHEP/Random/RanshiEngine.h>
0155 #include <CLHEP/Random/NonRandomEngine.h>
0156
0157 /// Initialize the instance.
0158 void Geant4Random::initialize()   {
0159   if ( !m_file.empty() )  {
0160     std::ifstream in(m_file.c_str(), std::ios::in);
0161     m_engine = CLHEP::EngineFactory::newEngine(in);
0162     if ( !m_engine )    {
0163       except("Failed to create CLHEP random engine from file:%s.",m_file.c_str());
0164     }
0165     m_seed = m_engine->getSeed();
0166   }
0167   else if ( !m_engineType.empty() )    {
0168     /// Create new engine if a type is specified
0169     if ( m_engineType == CLHEP::HepJamesRandom::engineName() )
0170       m_engine = new CLHEP::HepJamesRandom();
0171     else if ( m_engineType == CLHEP::RanecuEngine::engineName() )
0172       m_engine = new CLHEP::RanecuEngine();
0173     else if ( m_engineType == CLHEP::Ranlux64Engine::engineName() )
0174       m_engine = new CLHEP::Ranlux64Engine();
0175     else if ( m_engineType == CLHEP::MTwistEngine::engineName() )
0176       m_engine = new CLHEP::MTwistEngine();
0177     else if ( m_engineType == CLHEP::DualRand::engineName() )
0178       m_engine = new CLHEP::DualRand();
0179     else if ( m_engineType == CLHEP::RanluxEngine::engineName() )
0180       m_engine = new CLHEP::RanluxEngine();
0181     else if ( m_engineType == CLHEP::RanshiEngine::engineName() )
0182       m_engine = new CLHEP::RanshiEngine();
0183     else if ( m_engineType == CLHEP::NonRandomEngine::engineName() )
0184       m_engine = new CLHEP::NonRandomEngine();
0185
0186     if ( !m_engine )    {
0187       except("Failed to create CLHEP random engine of type: %s.",m_engineType.c_str());
0188     }
0189   }
0190   m_engine->setSeed(m_seed,m_luxury);
0191   m_rootRandom = new RNDM(this);
0192   m_inited = true;
0193   if ( 0 == s_instance )   {
0194     setMainInstance(this);
0195   }
0196 }
0197
0198 /// Should initialise the status of the algorithm according to seed.
0199 void Geant4Random::setSeed(long seed)   {
0200   if ( !m_inited ) initialize();
0201   m_engine->setSeed(m_seed=seed,0);
0202 }
0203
0204 /// Should initialise the status of the algorithm
0205 /** Initialization according to the zero terminated
0206  *  array of seeds. It is allowed to ignore one or
0207  *  many seeds in this array.
0208  */
0209 void Geant4Random::setSeeds(const long* seeds, int size)   {
0210   if ( !m_inited ) initialize();
0211   m_seed = seeds[0];
0212   m_engine->setSeeds(seeds, size);
0213 }
0214
0215 /// Should save on a file specific to the instantiated engine in use the current status.
0216 void Geant4Random::saveStatus( const char filename[] ) const    {
0217   if ( !m_inited )  {
0218     except("Failed to save RandomGenerator status. [Not-inited]");
0219   }
0220   m_engine->saveStatus(filename);
0221 }
0222
0223 /// Should read from a file and restore the last saved engine configuration.
0224 void Geant4Random::restoreStatus( const char filename[] )    {
0225   if ( !m_inited ) initialize();
0226   m_engine->restoreStatus(filename);
0227 }
0228
0229 /// Should dump the current engine status on the screen.
0230 void Geant4Random::showStatus() const    {
0231   if ( !m_inited )  {
0232     error("Failed to show RandomGenerator status. [Not-inited]");
0233     return;
0234   }
0235   printP2("Random engine status of object of type Geant4Random @ 0x%p",this);
0236   if ( !m_file.empty() )
0237     printP2("   Created from file: %s",m_file.c_str());
0238   else if ( !m_engineType.empty() )
0239     printP2("   Special instance created of type:%s @ 0x%p",
0240             m_engineType.c_str(),m_engine);
0241   else
0242     printP2("   Reused HepRandom engine instance %s @ 0x%p",
0243             m_engine ? m_engine->name().c_str() : "???", m_engine);
0244
0245   if ( m_engine == CLHEP::HepRandom::getTheEngine() )
0246     printP2("   Instance is identical to Geant4's HepRandom instance.");
0247
0248   printP2("   Instance is %sidentical to ROOT's gRandom instance.",
0249           gRandom == m_rootRandom ? "" : "NOT ");
0250
0251   if ( gRandom != m_rootRandom )   {
0252     printP2("      Local TRandom: 0x%p  gRandom: 0x%p",m_rootRandom,gRandom);
0253   }
0254   if ( 0 == m_engine )   {
0255     error("   Geant4Random instance has not engine attached!");
0256     return;
0257   }
0258   m_engine->showStatus();
0259 }
0260
0261 /// Create flat distributed random numbers in the interval ]0,1]
0262 double Geant4Random::rndm_clhep()  {
0263   if ( !m_inited ) initialize();
0264   return m_engine->flat();
0265 }
0266
0267 /// Create flat distributed random numbers in the interval ]0,1]
0268 double Geant4Random::rndm(int i)  {
0269   if ( !m_inited ) initialize();
0270   return gRandom->Rndm(i);
0271 }
0272
0273 /// Create a float array of flat distributed random numbers in the interval ]0,1]
0274 void   Geant4Random::rndmArray(int n, float *array)  {
0275   if ( !m_inited ) initialize();
0276   gRandom->RndmArray(n,array);
0277 }
0278
0279 /// Create a double array of flat distributed random numbers in the interval ]0,1]
0280 void   Geant4Random::rndmArray(int n, double *array)  {
0281   if ( !m_inited ) initialize();
0282   gRandom->RndmArray(n,array);
0283 }
0284
0285 /// Create uniformly disributed random numbers in the interval ]0,x1]
0286 double Geant4Random::uniform(double x1)  {
0287   if ( !m_inited ) initialize();
0288   return gRandom->Uniform(x1);
0289 }
0290
0291 /// Create uniformly disributed random numbers in the interval ]x1,x2]
0292 double Geant4Random::uniform(double x1, double x2)  {
0293   if ( !m_inited ) initialize();
0294   return gRandom->Uniform(x1,x2);
0295 }
0296
0297 /// Create exponentially distributed random numbers
0298 double Geant4Random::exp(double tau)  {
0299   if ( !m_inited ) initialize();
0300   return gRandom->Exp(tau);
0301 }
0302
0303 /// Generates random vectors, uniformly distributed over a circle of given radius.
0304 double Geant4Random::gauss(double mean, double sigma)  {
0305   if ( !m_inited ) initialize();
0306   return gRandom->Gaus(mean,sigma);
0307 }
0308
0309 /// Create landau distributed random numbers
0310 double Geant4Random::landau(double mean, double sigma)  {
0311   if ( !m_inited ) initialize();
0312   return gRandom->Landau(mean,sigma);
0313 }
0314
0315 /// Create tuple of randum number around a circle with radius r
0316 void   Geant4Random::circle(double &x, double &y, double r)  {
0317   if ( !m_inited ) initialize();
0318   gRandom->Circle(x,y,r);
0319 }
0320
0321 /// Create tuple of randum number on a sphere with radius r
0322 void   Geant4Random::sphere(double &x, double &y, double &z, double r)  {
0323   if ( !m_inited ) initialize();
0324   gRandom->Sphere(x,y,z,r);
0325 }
0326
0327 /// Create poisson distributed random numbers
0328 double Geant4Random::poisson(double mean)   {
0329   if ( !m_inited ) initialize();
0330   return gRandom->PoissonD(mean);
0331 }
0332
0333 /// Create breit wigner distributed random numbers
0334 double Geant4Random::breit_wigner(double mean, double gamma)   {
0335   if ( !m_inited ) initialize();
0336   return gRandom->BreitWigner(mean, gamma);
0337 }
0338
0339 /// Create gamma distributed random numbers
0340 double Geant4Random::gamma(double k, double lambda)    {
0341   if ( !m_inited ) initialize();
0342   return CLHEP::RandGamma::shoot(this->m_engine, k, lambda);
0343 }