Geant4Random.cpp

Go to the documentation of this file.

//==========================================================================
//  AIDA Detector description implementation 
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author     : M.Frank
//
//==========================================================================
 
// Framework include files
#include <DD4hep/Printout.h>
#include <DD4hep/InstanceCount.h>
#include <DDG4/Geant4Random.h>
 
#include <CLHEP/Random/EngineFactory.h>
#include <CLHEP/Random/RandGamma.h>
#include <CLHEP/Random/Random.h>
 
// ROOT include files
#include <TRandom1.h>
 
// C/C++ include files
#include <cmath>
 
using namespace dd4hep::sim;
 
namespace CLHEP   {
  unsigned long crc32ul(const std::string& s);
}
 
namespace    {
 
  class RNDM : public TRandom  {
    Geant4Random* m_generator;
    CLHEP::HepRandomEngine* m_engine;
    
  public:
    RNDM(Geant4Random* r) : TRandom(), m_generator(r)    {
      m_engine = m_generator->engine();
    }
    virtual ~RNDM() {      }
    virtual  void SetSeed(UInt_t seed=0)  final  {
      fSeed = seed;
      m_generator->setSeed((long)seed);
    }
    virtual  void SetSeed(ULong_t seed=0)  final  {
      fSeed = seed;
      m_generator->setSeed((long)seed);
    }
    virtual Double_t Rndm()  final {
      return m_engine->flat();
    }
    virtual Double_t Rndm(Int_t)  final  {
      return m_engine->flat();
    }
    virtual void RndmArray(Int_t size, Float_t *array)  final  {
      for (Int_t i=0;i<size;i++) array[i] = m_engine->flat();
    }
    virtual void RndmArray(Int_t size, Double_t *array)  final  {
      m_engine->flatArray(size,array);
    }
  };
  // Thread-local so each worker thread has its own Geant4Random instance.
  static thread_local Geant4Random* s_instance = nullptr;
}
 
Geant4Random::Geant4Random(Geant4Context* ctxt, const std::string& nam)
  : Geant4Action(ctxt,nam), m_engine(0), m_rootRandom(0), m_rootOLD(0), 
    m_inited(false)
{
  declareProperty("File",   m_file="");
  declareProperty("Type",   m_engineType="");
  declareProperty("Seed",   m_seed = 123456789);
  declareProperty("Luxury", m_luxury = 1);
  declareProperty("Replace_gRandom",  m_replace = true);
  // Default: static Geant4 random engine.
  m_engine = CLHEP::HepRandom::getTheEngine();
  InstanceCount::increment(this);
}
 
Geant4Random::~Geant4Random()  {
  // Only delete the engine if it is NOT the CLEP default one
  // BUT: Just cannot delete the engine. Causes havoc with static destructors!
  //CLHEP::HepRandomEngine* curr = CLHEP::HepRandom::getTheEngine();
  //if ( !m_engineType.empty() && m_engine != curr ) deletePtr(m_engine);
 
  // Set gRandom to the old value
  if (  m_rootRandom == gRandom ) gRandom = m_rootOLD;
  // Reset instance pointer
  if (  s_instance == this ) s_instance = 0;
  // Finally delete the TRandom instance wrapper
  detail::deletePtr(m_rootRandom);
  InstanceCount::decrement(this);
}
 
Geant4Random* Geant4Random::instance(bool throw_exception)   {
  if ( !s_instance && throw_exception )  {
    dd4hep::except("Geant4Random", "No global random number generator defined!");
  }  
  return s_instance;
}
 
Geant4Random* Geant4Random::setMainInstance(Geant4Random* ptr)   {
  if ( ptr && !ptr->m_inited )  {
    ptr->initialize();
  }
  if ( s_instance != ptr )   {
    if ( !ptr )  {
      dd4hep::except("Geant4Random","Attempt to declare invalid Geant4Random instance.");
    }
    if ( !ptr->m_inited )  {  
      dd4hep::except("Geant4Random","Attempt to declare uninitialized Geant4Random instance.");
    }
    Geant4Random* old = s_instance;
    CLHEP::HepRandomEngine* curr = CLHEP::HepRandom::getTheEngine();
    if ( ptr->m_engine != curr )     {       
      ptr->printP2("Moving CLHEP random instance from %p to %p",curr,ptr->m_engine); 
      CLHEP::HepRandom::setTheEngine(ptr->m_engine);
    }
    if ( ptr->m_replace )   {
      ptr->m_rootOLD = gRandom;
      gRandom = ptr->m_rootRandom;
    }
    s_instance = ptr;
    return old;
  }
  return 0;
}
 
#include <CLHEP/Random/DualRand.h>
#include <CLHEP/Random/JamesRandom.h>
#include <CLHEP/Random/MTwistEngine.h>
#include <CLHEP/Random/RanecuEngine.h>
#include <CLHEP/Random/Ranlux64Engine.h>
#include <CLHEP/Random/RanluxEngine.h>
#include <CLHEP/Random/RanshiEngine.h>
#include <CLHEP/Random/NonRandomEngine.h>
 
void Geant4Random::initialize()   {
  if ( !m_file.empty() )  {
    std::ifstream in(m_file.c_str(), std::ios::in);
    m_engine = CLHEP::EngineFactory::newEngine(in);
    if ( !m_engine )    {
      except("Failed to create CLHEP random engine from file:%s.",m_file.c_str());
    }
    m_seed = m_engine->getSeed();
  }
  else if ( !m_engineType.empty() )    {
    if ( m_engineType == CLHEP::HepJamesRandom::engineName() )
      m_engine = new CLHEP::HepJamesRandom();
    else if ( m_engineType == CLHEP::RanecuEngine::engineName() )
      m_engine = new CLHEP::RanecuEngine();
    else if ( m_engineType == CLHEP::Ranlux64Engine::engineName() )
      m_engine = new CLHEP::Ranlux64Engine();
    else if ( m_engineType == CLHEP::MTwistEngine::engineName() )
      m_engine = new CLHEP::MTwistEngine();
    else if ( m_engineType == CLHEP::DualRand::engineName() )
      m_engine = new CLHEP::DualRand();
    else if ( m_engineType == CLHEP::RanluxEngine::engineName() )
      m_engine = new CLHEP::RanluxEngine();
    else if ( m_engineType == CLHEP::RanshiEngine::engineName() )
      m_engine = new CLHEP::RanshiEngine();
    else if ( m_engineType == CLHEP::NonRandomEngine::engineName() )
      m_engine = new CLHEP::NonRandomEngine();
 
    if ( !m_engine )    {
      except("Failed to create CLHEP random engine of type: %s.",m_engineType.c_str());
    }
  }
  m_engine->setSeed(m_seed,m_luxury);
  m_rootRandom = new RNDM(this);
  m_inited = true;
  if ( 0 == s_instance )   {
    setMainInstance(this);
  }
}
 
void Geant4Random::setSeed(long seed)   {
  if ( !m_inited ) initialize();
  m_engine->setSeed(m_seed=seed,0);
}
 
 
void Geant4Random::setSeeds(const long* seeds, int size)   {
  if ( !m_inited ) initialize();
  m_seed = seeds[0];
  m_engine->setSeeds(seeds, size);
}
 
void Geant4Random::saveStatus( const char filename[] ) const    {
  if ( !m_inited )  {
    except("Failed to save RandomGenerator status. [Not-inited]");
  }
  m_engine->saveStatus(filename);
}
 
void Geant4Random::restoreStatus( const char filename[] )    {
  if ( !m_inited ) initialize();
  m_engine->restoreStatus(filename);
}
 
void Geant4Random::showStatus() const    {
  if ( !m_inited )  {
    error("Failed to show RandomGenerator status. [Not-inited]");
    return;
  }
  printP2("Random engine status of object of type Geant4Random @ 0x%p",this);
  if ( !m_file.empty() )
    printP2("   Created from file: %s",m_file.c_str());
  else if ( !m_engineType.empty() )
    printP2("   Special instance created of type:%s @ 0x%p",
            m_engineType.c_str(),m_engine);
  else
    printP2("   Reused HepRandom engine instance %s @ 0x%p",
            m_engine ? m_engine->name().c_str() : "???", m_engine);
  
  if ( m_engine == CLHEP::HepRandom::getTheEngine() )
    printP2("   Instance is identical to Geant4's HepRandom instance.");
 
  printP2("   Instance is %sidentical to ROOT's gRandom instance.",
          gRandom == m_rootRandom ? "" : "NOT ");
 
  if ( gRandom != m_rootRandom )   {
    printP2("      Local TRandom: 0x%p  gRandom: 0x%p",m_rootRandom,gRandom);
  }
  if ( 0 == m_engine )   {
    error("   Geant4Random instance has not engine attached!");
    return;
  }
  m_engine->showStatus();
}
 
double Geant4Random::rndm_clhep()  {
  if ( !m_inited ) initialize();
  return m_engine->flat();
}
 
double Geant4Random::rndm(int i)  {
  if ( !m_inited ) initialize();
  return m_rootRandom->Rndm(i);
}
 
void   Geant4Random::rndmArray(int n, float *array)  {
  if ( !m_inited ) initialize();
  m_rootRandom->RndmArray(n,array);
}
 
void   Geant4Random::rndmArray(int n, double *array)  {
  if ( !m_inited ) initialize();
  m_rootRandom->RndmArray(n,array);
}
 
double Geant4Random::uniform(double x1)  {
  if ( !m_inited ) initialize();
  return m_rootRandom->Uniform(x1);
}
 
double Geant4Random::uniform(double x1, double x2)  {
  if ( !m_inited ) initialize();
  return m_rootRandom->Uniform(x1,x2);
}
 
double Geant4Random::exp(double tau)  {
  if ( !m_inited ) initialize();
  return m_rootRandom->Exp(tau);
}
 
double Geant4Random::gauss(double mean, double sigma)  {
  if ( !m_inited ) initialize();
  return m_rootRandom->Gaus(mean,sigma);
}
 
double Geant4Random::landau(double mean, double sigma)  {
  if ( !m_inited ) initialize();
  return m_rootRandom->Landau(mean,sigma);
}
 
void   Geant4Random::circle(double &x, double &y, double r)  {
  if ( !m_inited ) initialize();
  m_rootRandom->Circle(x,y,r);
}
 
void   Geant4Random::sphere(double &x, double &y, double &z, double r)  {
  if ( !m_inited ) initialize();
  m_rootRandom->Sphere(x,y,z,r);
}
 
double Geant4Random::poisson(double mean)   {
  if ( !m_inited ) initialize();
  return m_rootRandom->PoissonD(mean);
}
 
double Geant4Random::breit_wigner(double mean, double gamma)   {
  if ( !m_inited ) initialize();
  return m_rootRandom->BreitWigner(mean, gamma);
}
 
double Geant4Random::gamma(double k, double lambda)    {
  if ( !m_inited ) initialize();
  return CLHEP::RandGamma::shoot(this->m_engine, k, lambda);
}