// $Id: RandLandau.cc,v 1.3.4.1 2005/03/18 22:26:48 garren Exp $ // -*- C++ -*- // // ----------------------------------------------------------------------- // HEP Random // --- RandLandau --- // class implementation file // ----------------------------------------------------------------------- // ======================================================================= // M Fischler - Created 1/6/2000. // // The key transform() method uses the algorithm in CERNLIB. // This is because I trust that RANLAN routine more than // I trust the Bukin-Grozina inverseLandau, which is not // claimed to be better than 1% accurate. // // M Fischler - put and get to/from streams 12/13/04 // ======================================================================= #include "CLHEP/Random/defs.h" #include "CLHEP/Random/RandLandau.h" #include #include // for log() namespace CLHEP { std::string RandLandau::name() const {return "RandLandau";} HepRandomEngine & RandLandau::engine() {return *localEngine;} RandLandau::~RandLandau() { if ( deleteEngine ) delete localEngine; } RandLandau::RandLandau(const RandLandau& right) {} void RandLandau::shootArray( const int size, double* vect ) { int i; for (i=0; i= 70 && index <= 800 ) { // (A) double f0 = inverseLandau [index]; double f1 = inverseLandau [index+1]; return f0 + du * (f1 - f0); } else if ( index >= 7 && index <= 980 ) { // (B) double f_1 = inverseLandau [index-1]; double f0 = inverseLandau [index]; double f1 = inverseLandau [index+1]; double f2 = inverseLandau [index+2]; return f0 + du * (f1 - f0 - .25*(1-du)* (f2 -f1 - f0 + f_1) ); } else if ( index < 7 ) { // (C) const double n0 = 0.99858950; const double n1 = 34.5213058; const double d1 = 34.1760202; const double n2 = 17.0854528; const double d2 = 4.01244582; double logr = log(r); double x = 1/logr; double x2 = x*x; double pade = (n0 + n1*x + n2*x2) / (1.0 + d1*x + d2*x2); return ( - log ( -.91893853 - logr ) -1 ) * pade; } else if ( index <= 999 ) { // (D) const double n0 = 1.00060006; const double n1 = 263.991156; const double d1 = 257.368075; const double n2 = 4373.20068; const double d2 = 3414.48018; double x = 1-r; double x2 = x*x; return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2)); } else { // (E) const double n0 = 1.00001538; const double n1 = 6075.14119; const double d1 = 6065.11919; const double n2 = 734266.409; const double d2 = 694021.044; double x = 1-r; double x2 = x*x; return (n0 + n1*x + n2*x2) / (x * (1.0 + d1*x + d2*x2)); } } // transform() std::ostream & RandLandau::put ( std::ostream & os ) const { int pr=os.precision(20); os << " " << name() << "\n"; os.precision(pr); return os; } std::istream & RandLandau::get ( std::istream & is ) { std::string inName; is >> inName; if (inName != name()) { is.clear(std::ios::badbit | is.rdstate()); std::cerr << "Mismatch when expecting to read state of a " << name() << " distribution\n" << "Name found was " << inName << "\nistream is left in the badbit state\n"; return is; } return is; } } // namespace CLHEP