Geant4Random.cpp

//==========================================================================
//  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 std;
using namespace dd4hep::sim;

namespace CLHEP   {
  unsigned long crc32ul(const std::string& s);
}

namespace    {

  /// Helper class to redirect calls to gRandon to Geant4Random
  class RNDM : public TRandom  {
    /// Reference to the generator object
    Geant4Random* m_generator;
    /// Reference to HepRandomEngine
    CLHEP::HepRandomEngine* m_engine;
    
  public:
    /// Initializing constructor
    RNDM(Geant4Random* r) : TRandom(), m_generator(r)    {
      m_engine = m_generator->engine();
    }
    /// Default destructor
    virtual ~RNDM() {      }
    /// Set new seed
    virtual  void SetSeed(UInt_t seed=0)  final  {
      fSeed = seed;
      m_generator->setSeed((long)seed);
    }
    /// Set new seed
    virtual  void SetSeed(ULong_t seed=0)  final  {
      fSeed = seed;
      m_generator->setSeed((long)seed);
    }
    /// Single shot random number creation
    virtual Double_t Rndm()  final {
      return m_engine->flat();
    }
    /// Single shot random number creation
    virtual Double_t Rndm(Int_t)  final  {
      return m_engine->flat();
    }
    /// Return an array of n random numbers uniformly distributed in ]0,1].
    virtual void RndmArray(Int_t size, Float_t *array)  final  {
      for (Int_t i=0;i<size;i++) array[i] = m_engine->flat();
    }
    /// Return an array of n random numbers uniformly distributed in ]0,1].
    virtual void RndmArray(Int_t size, Double_t *array)  final  {
      m_engine->flatArray(size,array);
    }
  };
  static Geant4Random* s_instance = 0;
}

/// Default constructor
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);
}

/// Default destructor
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);
}

/// Access the main Geant4 random generator instance. Must be created before used!
Geant4Random* Geant4Random::instance(bool throw_exception)   {
  if ( !s_instance && throw_exception )  {
    throw runtime_error("No global random number generator defined!");
  }  
  return s_instance;
}

/// Make this random generator instance the one used by Geant4
Geant4Random* Geant4Random::setMainInstance(Geant4Random* ptr)   {
  if ( ptr && !ptr->m_inited )  {
    ptr->initialize();
  }
  if ( s_instance != ptr )   {
    if ( !ptr )  {
      throw runtime_error("Attempt to declare invalid Geant4Random instance.");
    }
    if ( !ptr->m_inited )  {  
      throw runtime_error("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"

/// Initialize the instance. 
void Geant4Random::initialize()   {
  if ( !m_file.empty() )  {
    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() )    {
    /// Create new engine if a type is specified
    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);
  }
}

/// Should initialise the status of the algorithm according to seed.
void Geant4Random::setSeed(long seed)   {
  if ( !m_inited ) initialize();
  m_engine->setSeed(m_seed=seed,0);
}

/// Should initialise the status of the algorithm
/** Initialization according to the zero terminated
 *  array of seeds. It is allowed to ignore one or 
 *  many seeds in this array.
 */
void Geant4Random::setSeeds(const long* seeds, int size)   {
  if ( !m_inited ) initialize();
  m_seed = seeds[0];
  m_engine->setSeeds(seeds, size);
}

/// Should save on a file specific to the instantiated engine in use the current status.
void Geant4Random::saveStatus( const char filename[] ) const    {
  if ( !m_inited )  {
    except("Failed to save RandomGenerator status. [Not-inited]");
  }
  m_engine->saveStatus(filename);
}

/// Should read from a file and restore the last saved engine configuration.
void Geant4Random::restoreStatus( const char filename[] )    {
  if ( !m_inited ) initialize();
  m_engine->restoreStatus(filename);
}

/// Should dump the current engine status on the screen.
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();
}

/// Create flat distributed random numbers in the interval ]0,1]
double Geant4Random::rndm_clhep()  {
  if ( !m_inited ) initialize();
  return m_engine->flat();
}

/// Create flat distributed random numbers in the interval ]0,1]
double Geant4Random::rndm(int i)  {
  if ( !m_inited ) initialize();
  return gRandom->Rndm(i);
}

/// Create a float array of flat distributed random numbers in the interval ]0,1]
void   Geant4Random::rndmArray(int n, float *array)  {
  if ( !m_inited ) initialize();
  gRandom->RndmArray(n,array);
}

/// Create a double array of flat distributed random numbers in the interval ]0,1]
void   Geant4Random::rndmArray(int n, double *array)  {
  if ( !m_inited ) initialize();
  gRandom->RndmArray(n,array);
}

/// Create uniformly disributed random numbers in the interval ]0,x1]
double Geant4Random::uniform(double x1)  {
  if ( !m_inited ) initialize();
  return gRandom->Uniform(x1);
}

/// Create uniformly disributed random numbers in the interval ]x1,x2]
double Geant4Random::uniform(double x1, double x2)  {
  if ( !m_inited ) initialize();
  return gRandom->Uniform(x1,x2);
}

/// Create exponentially distributed random numbers
double Geant4Random::exp(double tau)  {
  if ( !m_inited ) initialize();
  return gRandom->Exp(tau);
}

/// Generates random vectors, uniformly distributed over a circle of given radius.
double Geant4Random::gauss(double mean, double sigma)  {
  if ( !m_inited ) initialize();
  return gRandom->Gaus(mean,sigma);
}

/// Create landau distributed random numbers
double Geant4Random::landau(double mean, double sigma)  {
  if ( !m_inited ) initialize();
  return gRandom->Landau(mean,sigma);
}

/// Create tuple of randum number around a circle with radius r
void   Geant4Random::circle(double &x, double &y, double r)  {
  if ( !m_inited ) initialize();  
  gRandom->Circle(x,y,r);
}

/// Create tuple of randum number on a sphere with radius r
void   Geant4Random::sphere(double &x, double &y, double &z, double r)  {
  if ( !m_inited ) initialize();
  gRandom->Sphere(x,y,z,r);
}

/// Create poisson distributed random numbers
double Geant4Random::poisson(double mean)   {
  if ( !m_inited ) initialize();
  return gRandom->PoissonD(mean);
}

/// Create breit wigner distributed random numbers
double Geant4Random::breit_wigner(double mean, double gamma)   {
  if ( !m_inited ) initialize();
  return gRandom->BreitWigner(mean, gamma);
}

/// Create gamma distributed random numbers
double Geant4Random::gamma(double k, double lambda)    {
  if ( !m_inited ) initialize();
  return CLHEP::RandGamma::shoot(this->m_engine, k, lambda);
}