CLHEP::RandPoissonQ Class Reference

#include <RandPoissonQ.h>

Inheritance diagram for CLHEP::RandPoissonQ:

Public Member Functions
HepRandomEngine &	engine ()
long	fire ()
long	fire (double mean)
void	fireArray (const int size, long *vect)
void	fireArray (const int size, long *vect, double mean)
double	flat ()
double	flat (HepRandomEngine *theNewEngine)
void	flatArray (const int size, double *vect)
void	flatArray (HepRandomEngine *theNewEngine, const int size, double *vect)
std::istream &	get (std::istream &is)
std::string	name () const
double	operator() ()
double	operator() (double mean)
std::ostream &	put (std::ostream &os) const
	RandPoissonQ (HepRandomEngine &anEngine, double b1=1.0)
	RandPoissonQ (HepRandomEngine *anEngine, double b1=1.0)
virtual	~RandPoissonQ ()

Static Public Member Functions
static int	createInstance ()
static std::string	distributionName ()
static HepRandomEngine *	getTheEngine ()
static HepRandom *	getTheGenerator ()
static long	getTheSeed ()
static const long *	getTheSeeds ()
static void	getTheTableSeeds (long *seeds, int index)
static std::istream &	restoreDistState (std::istream &is)
static void	restoreEngineStatus (const char filename[]="Config.conf")
static std::istream &	restoreFullState (std::istream &is)
static std::istream &	restoreStaticRandomStates (std::istream &is)
static std::ostream &	saveDistState (std::ostream &os)
static void	saveEngineStatus (const char filename[]="Config.conf")
static std::ostream &	saveFullState (std::ostream &os)
static std::ostream &	saveStaticRandomStates (std::ostream &os)
static void	setTheEngine (HepRandomEngine *theNewEngine)
static void	setTheSeed (long seed, int lxr=3)
static void	setTheSeeds (const long *seeds, int aux=-1)
static long	shoot (double mean=1.0)
static long	shoot (HepRandomEngine *anEngine, double mean=1.0)
static void	shootArray (const int size, long *vect, double mean=1.0)
static void	shootArray (HepRandomEngine *anEngine, const int size, long *vect, double mean=1.0)
static void	showEngineStatus ()
static int	tableBoundary ()

Static Public Attributes
static const double	MAXIMUM_POISSON_DEVIATE = 2.0E9

Protected Member Functions
HepRandomEngine *	getLocalEngine ()

Static Protected Member Functions
static double	getMaxMean ()
static double	getOldMean ()
static double *	getPStatus ()
static void	setOldMean (double val)
static void	setPStatus (double sq, double alxm, double g1)

Protected Attributes
double	defaultMean
double	meanMax

Static Protected Attributes
static const long	seedTable [215][2]

Private Member Functions
void	setupForDefaultMu ()

Static Private Member Functions
static long	poissonDeviateQuick (HepRandomEngine *e, double A0, double A1, double A2, double sig)
static long	poissonDeviateQuick (HepRandomEngine *e, double mean)
static long	poissonDeviateSmall (HepRandomEngine *e, double mean)

Private Attributes
double	a0
double	a1
double	a2
std::shared_ptr< HepRandomEngine >	localEngine
double	oldm
double	sigma
double	status [3]

Static Private Attributes
static const int	BELOW = 30
static const int	ENTRIES = 51
static const double	FIRST_MU = 10
static const double	LAST_MU = 95
static const double	meanMax_st = 2.0E9
static CLHEP_THREAD_LOCAL double	oldm_st = -1.0
static const double	S = 5
static CLHEP_THREAD_LOCAL double	status_st [3] = {0., 0., 0.}

Detailed Description

Author

Definition at line 31 of file RandPoissonQ.h.

Constructor & Destructor Documentation

RandPoissonQ() [1/2]

CLHEP::RandPoissonQ::RandPoissonQ ( HepRandomEngine &  anEngine, double  b1 = 1.0  )

inline

RandPoissonQ() [2/2]

CLHEP::RandPoissonQ::RandPoissonQ ( HepRandomEngine *  anEngine, double  b1 = 1.0  )

inline

~RandPoissonQ()

CLHEP::RandPoissonQ::~RandPoissonQ ( )

virtual

Definition at line 82 of file RandPoissonQ.cc.

                            {
}

Member Function Documentation

createInstance()

int CLHEP::HepRandom::createInstance ( )

staticinherited

Definition at line 311 of file Random.cc.

{
  return static_cast<int>( theDefaults().ensureInitialized() );
}

References CLHEP::anonymous_namespace{Random.cc}::defaults::ensureInitialized(), and CLHEP::anonymous_namespace{Random.cc}::theDefaults().

distributionName()

static std::string CLHEP::RandPoissonQ::distributionName ( )

inlinestatic

Definition at line 101 of file RandPoissonQ.h.

{return "RandPoissonQ";}  

engine()

HepRandomEngine & CLHEP::RandPoissonQ::engine ( )

virtual

Reimplemented from CLHEP::RandPoisson.

Definition at line 52 of file RandPoissonQ.cc.

{return RandPoisson::engine();}

References CLHEP::RandPoisson::engine().

fire() [1/2]

long CLHEP::RandPoissonQ::fire

(

)

Definition at line 135 of file RandPoissonQ.cc.

                        {
  if ( defaultMean < LAST_MU + S ) {
    return poissonDeviateSmall ( getLocalEngine(), defaultMean );
  } else {
    return poissonDeviateQuick ( getLocalEngine(), a0, a1, a2, sigma );
  }
} // fire()

References a0, a1, a2, CLHEP::RandPoisson::defaultMean, CLHEP::RandPoisson::getLocalEngine(), LAST_MU, poissonDeviateQuick(), poissonDeviateSmall(), S, and sigma.

Referenced by fireArray(), and operator()().

fire() [2/2]

long CLHEP::RandPoissonQ::fire

(

double

mean

)

Definition at line 131 of file RandPoissonQ.cc.

                                   {
  return shoot(getLocalEngine(), mean);
}

References CLHEP::RandPoisson::getLocalEngine(), and shoot().

fireArray() [1/2]

void CLHEP::RandPoissonQ::fireArray	(	const int	size,
		long *	vect
	)

Definition at line 187 of file RandPoissonQ.cc.

                                                       {
  for( long* v = vect; v != vect + size; ++v )
    *v = fire( defaultMean );
}

References CLHEP::RandPoisson::defaultMean, and fire().

fireArray() [2/2]

void CLHEP::RandPoissonQ::fireArray	(	const int	size,
		long *	vect,
		double	mean
	)

Definition at line 182 of file RandPoissonQ.cc.

                                                                 {
  for( long* v = vect; v != vect + size; ++v )
    *v = fire( m );
}

References fire(), and CLHEP::m.

flat() [1/2]

double CLHEP::HepRandom::flat ( )

inherited

Definition at line 199 of file Random.cc.

{
  return theDefaults().theEngine->flat();
}

References CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theEngine.

Referenced by CLHEP::RandExponential::shoot().

flat() [2/2]

double CLHEP::HepRandom::flat ( HepRandomEngine *  theNewEngine )

inlineinherited

flatArray() [1/2]

void CLHEP::HepRandom::flatArray ( const int  size, double *  vect  )

inherited

Definition at line 204 of file Random.cc.

{
  theDefaults().theEngine->flatArray(size,vect);
}

References CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theEngine.

flatArray() [2/2]

void CLHEP::HepRandom::flatArray ( HepRandomEngine *  theNewEngine, const int  size, double *  vect  )

inlineinherited

get()

std::istream & CLHEP::RandPoissonQ::get ( std::istream &  is )

virtual

Reimplemented from CLHEP::RandPoisson.

Definition at line 574 of file RandPoissonQ.cc.

                                                 {
  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;
  }
  if (possibleKeywordInput(is, "Uvec", a0)) {
    std::vector<unsigned long> t(2);
    is >> a0 >> t[0] >> t[1]; a0 = DoubConv::longs2double(t); 
    is >> a1 >> t[0] >> t[1]; a1 = DoubConv::longs2double(t); 
    is >> a2 >> t[0] >> t[1]; a2 = DoubConv::longs2double(t); 
    is >> sigma >> t[0] >> t[1]; sigma = DoubConv::longs2double(t); 
    RandPoisson::get(is);
    return is;
  }
  // is >> a0 encompassed by possibleKeywordInput
  is >> a1 >> a2 >> sigma;
  RandPoisson::get(is);
  return is;
}

References a0, a1, a2, CLHEP::RandPoisson::get(), CLHEP::DoubConv::longs2double(), name(), CLHEP::possibleKeywordInput(), and sigma.

getLocalEngine()

HepRandomEngine * CLHEP::RandPoisson::getLocalEngine ( )

inlineprotectedinherited

Referenced by fire().

getMaxMean()

static double CLHEP::RandPoisson::getMaxMean ( )

inlinestaticprotectedinherited

Definition at line 103 of file RandPoisson.h.

{return meanMax_st;}

References CLHEP::RandPoisson::meanMax_st.

Referenced by CLHEP::RandPoisson::fire(), and CLHEP::RandPoisson::shoot().

getOldMean()

static double CLHEP::RandPoisson::getOldMean ( )

inlinestaticprotectedinherited

Definition at line 101 of file RandPoisson.h.

{return oldm_st;}

References CLHEP::RandPoisson::oldm_st.

Referenced by CLHEP::RandPoisson::shoot().

getPStatus()

static double * CLHEP::RandPoisson::getPStatus ( )

inlinestaticprotectedinherited

Definition at line 107 of file RandPoisson.h.

{return status_st;}

References CLHEP::RandPoisson::status_st.

Referenced by CLHEP::RandPoisson::shoot().

getTheEngine()

HepRandomEngine * CLHEP::HepRandom::getTheEngine ( )

staticinherited

Definition at line 268 of file Random.cc.

{
  return theDefaults().theEngine.get();
}

References CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theEngine.

Referenced by CLHEP::StaticRandomStates::restore(), CLHEP::RandFlat::restoreEngineStatus(), CLHEP::RandGauss::restoreEngineStatus(), CLHEP::RandFlat::saveEngineStatus(), CLHEP::RandGauss::saveEngineStatus(), CLHEP::RandExponential::shoot(), CLHEP::RandExpZiggurat::shoot(), CLHEP::RandFlat::shoot(), CLHEP::RandGauss::shoot(), CLHEP::RandGaussZiggurat::shoot(), CLHEP::RandChiSquare::shoot(), CLHEP::RandStudentT::shoot(), CLHEP::RandBreitWigner::shoot(), CLHEP::RandGamma::shoot(), CLHEP::RandPoisson::shoot(), shoot(), CLHEP::RandBinomial::shoot(), and CLHEP::RandFlat::shootArray().

getTheGenerator()

HepRandom * CLHEP::HepRandom::getTheGenerator ( )

staticinherited

Definition at line 263 of file Random.cc.

{
  return theDefaults().theGenerator.get();
}

References CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theGenerator.

getTheSeed()

long CLHEP::HepRandom::getTheSeed ( )

staticinherited

Definition at line 239 of file Random.cc.

{
  return theDefaults().theEngine->getSeed();
}

References CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theEngine.

getTheSeeds()

const long * CLHEP::HepRandom::getTheSeeds ( )

staticinherited

Definition at line 249 of file Random.cc.

{
  return theDefaults().theEngine->getSeeds();
}

References CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theEngine.

getTheTableSeeds()

void CLHEP::HepRandom::getTheTableSeeds ( long *  seeds, int  index  )

staticinherited

Definition at line 254 of file Random.cc.

{
  if ((index >= 0) && (index < 215)) {
    seeds[0] = seedTable[index][0];
    seeds[1] = seedTable[index][1];
  }
  else seeds = NULL;
}

Referenced by CLHEP::HepJamesRandom::HepJamesRandom(), CLHEP::MTwistEngine::MTwistEngine(), CLHEP::RanecuEngine::RanecuEngine(), CLHEP::Ranlux64Engine::Ranlux64Engine(), CLHEP::RanluxEngine::RanluxEngine(), and CLHEP::RanecuEngine::setSeed().

name()

std::string CLHEP::RandPoissonQ::name ( ) const

virtual

Reimplemented from CLHEP::RandPoisson.

Definition at line 51 of file RandPoissonQ.cc.

{return "RandPoissonQ";}

Referenced by source.g4viscp.G4Scene::create_scene(), get(), mcscore.MCParticle::printout(), put(), and source.g4viscp.G4Scene::update_scene().

operator()() [1/2]

double CLHEP::RandPoissonQ::operator() ( )

virtual

Reimplemented from CLHEP::RandPoisson.

Definition at line 123 of file RandPoissonQ.cc.

                                {
  return (double) fire();
}

References fire().

operator()() [2/2]

double CLHEP::RandPoissonQ::operator()

(

double

mean

)

Definition at line 127 of file RandPoissonQ.cc.

                                             {
  return (double) fire(mean);
}

References fire().

poissonDeviateQuick() [1/2]

long CLHEP::RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *  e, double  A0, double  A1, double  A2, double  sig  )

staticprivate

Definition at line 219 of file RandPoissonQ.cc.

                                                     {
//
// Quick Poisson deviate algorithm used by quick for large mu:
//
// The principle:  For very large mu, a poisson distribution can be approximated
// by a gaussian: return the integer part of mu + .5 + g where g is a unit 
// normal.  However, this yelds a miserable approximation at values as
// "large" as 100.  The primary problem is that the poisson distribution is 
// supposed to have a skew of 1/mu**2, and the zero skew of the Guassian 
// leads to errors of order as big as 1/mu**2.
//
// We substitute for the gaussian a quadratic function of that gaussian random.
// The expression looks very nearly like mu + .5 - 1/6 + g + g**2/(6*mu).  
// The small positive quadratic term causes the resulting variate to have 
// a positive skew; the -1/6 constant term is there to correct for this bias 
// in the mean.  By adjusting these two and the linear term, we can match the
// first three moments to high accuracy in 1/mu.
//
// The sigma used is not precisely sqrt(mu) since a rounded-off Gaussian
// has a second moment which is slightly larger than that of the Gaussian.  
// To compensate, sig is multiplied by a factor which is slightly less than 1.
 
  //  double g = RandGauss::shootQuick( e );   // TEMPORARY MOD:
  double g = RandGaussQ::shoot( e );   // Unit normal
  g *= sig;
  double p = A2*g*g + A1*g + A0;
  if ( p < 0 ) return 0;    // Shouldn't ever possibly happen since 
                // mean should not be less than 100, but
                // we check due to paranoia.
  if ( p > MAXIMUM_POISSON_DEVIATE ) p = MAXIMUM_POISSON_DEVIATE;
 
  return long(p);
 
} // poissonDeviateQuick ()

References CLHEP::g, MAXIMUM_POISSON_DEVIATE, and CLHEP::RandGaussQ::shoot().

poissonDeviateQuick() [2/2]

long CLHEP::RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *  e, double  mean  )

staticprivate

Definition at line 195 of file RandPoissonQ.cc.

                                                                       {
 
  // Compute the coefficients defining the quadratic transformation from a 
  // Gaussian to a Poisson:
 
  double sig2 = mu * (.9998654 - .08346/mu);
  double sig = std::sqrt(sig2);
    // The multiplier corrects for fact that discretization of the form
    // [gaussian+.5] increases the second moment by a small amount.
 
  double t = 1./sig2;
 
  double sa2 = t*(1./6.) + t*t*(1./324.);
  double sa1 = std::sqrt (1-2*sa2*sa2*sig2);
  double sa0 = mu + .5 - sig2 * sa2;
 
  // The formula will be sa0 + sa1*x + sa2*x*x where x has sigma of sq.
  // The coeffeicients are chosen to match the first THREE moments of the 
  // true Poisson distribution.
 
  return poissonDeviateQuick ( e, sa0, sa1, sa2, sig ); 
} 

References poissonDeviateQuick().

Referenced by fire(), poissonDeviateQuick(), and shoot().

poissonDeviateSmall()

long CLHEP::RandPoissonQ::poissonDeviateSmall ( HepRandomEngine *  e, double  mean  )

staticprivate

Definition at line 257 of file RandPoissonQ.cc.

                                                                        {
  long N1;
  long N2;
  // The following are for later use to form a secondary random s:
  double rRange;        // This will hold the interval between cdf for the
                // computed N1 and cdf for N1+1.
  double rRemainder = 0; // This will hold the length into that interval.
 
  // Coming in, mean should not be more than LAST_MU + S.  However, we will 
  // be paranoid and test for this:
 
  if ( mean > LAST_MU + S ) {
    return RandPoisson::shoot(e, mean);
  }
 
  if (mean <= 0) {
    return 0;           // Perhaps we ought to balk harder here!
  }
 
  // >>> 1 <<< 
  //    Generate the first random, which we always will need.
 
  double r = e->flat();
 
  // >>> 2 <<< 
  //    For small mean, below the start of the tables, 
  //    do the series for cdf directly.  
 
  // In this case, since we know the series will terminate relatively quickly, 
  // almost alwaye use a precomputed 1/N array without fear of overrunning it.
 
  static const double oneOverN[50] = 
  {    0,   1.,    1/2.,  1/3.,  1/4.,  1/5.,  1/6.,  1/7.,  1/8.,  1/9., 
   1/10.,  1/11.,  1/12., 1/13., 1/14., 1/15., 1/16., 1/17., 1/18., 1/19., 
   1/20.,  1/21.,  1/22., 1/23., 1/24., 1/25., 1/26., 1/27., 1/28., 1/29.,
   1/30.,  1/31.,  1/32., 1/33., 1/34., 1/35., 1/36., 1/37., 1/38., 1/39., 
   1/40.,  1/41.,  1/42., 1/43., 1/44., 1/45., 1/46., 1/47., 1/48., 1/49. };
 
 
  if ( mean < FIRST_MU ) {
 
    long N = 0;
    double term = std::exp(-mean);
    double cdf = term;
 
    if ( r < (1 - 1.0E-9) ) {
      //
      // **** This is a normal path: ****
      //
      // Except when r is very close to 1, it is certain that we will exceed r 
      // before the 30-th term in the series, so a simple while loop is OK.
      const double* oneOverNptr = oneOverN;
      while( cdf <= r ) {
        N++ ;  
        oneOverNptr++;
        term *= ( mean * (*oneOverNptr) );
        cdf  += term;
      }
      return N;
      //
      // **** ****
      //
    } else { // r is almost 1...
      // For r very near to 1 we would have to check that we don't fall
      // off the end of the table of 1/N.  Since this is very rare, we just
      // ignore the table and do the identical while loop, using explicit 
      // division.
      double cdf0;
      while ( cdf <= r ) {
        N++ ;  
        term *= ( mean / N );
        cdf0 = cdf;
        cdf  += term;
    if (cdf == cdf0) break; // Can't happen, but just in case...
      }
      return N;
    } // end of if ( r compared to (1 - 1.0E-9) )
 
  } // End of the code for mean < FIRST_MU
 
  // >>> 3 <<< 
  //    Find the row of the tables corresponding to the highest tabulated mu
  //    which is no greater than our actual mean.
 
  int rowNumber = int((mean - FIRST_MU)/S);
  const double * cdfs = &poissonTables [rowNumber*ENTRIES]; 
  double mu = FIRST_MU + rowNumber*S;
  double deltaMu = mean - mu;
  int Nmin = int(mu - BELOW);
  if (Nmin < 1) Nmin = 1;
  int Nmax = Nmin + (ENTRIES - 1);
 
 
  // >>> 4 <<< 
  //    If r is less that the smallest entry in the row, then 
  //    generate the deviate directly from the series.  
 
  if ( r < cdfs[0] ) {
  
    // In this case, we are tempted to use the actual mean, and not 
    // generate a second deviate to account for the leftover part mean - mu.
    // That would be an error, generating a distribution with enough excess
    // at Nmin + (mean-mu)/2 to be detectable in 4,000,000 trials.
 
    // Since this case is very rare (never more than .2% of the r values)
    // and can happen where N will be large (up to 65 for the mu=95 row)
    // we use explicit division so as to avoid having to worry about running
    // out of oneOverN table.
 
    long N = 0;
    double term = std::exp(-mu);
    double cdf = term;
    double cdf0;
 
    while(cdf <= r) {
        N++ ;  
        term *= ( mu / N );
        cdf0 = cdf;
        cdf  += term;
    if (cdf == cdf0) break; // Can't happen, but just in case...
    }
    N1 = N;
        // std::cout << r << "   " << N << "   ";
        // DBG_small = true;
    rRange = 0;     // In this case there is always a second r needed
 
  } // end of small-r case
 
 
  // >>> 5 <<< 
  //    Assuming r lies within the scope of the row for this mu, find the 
  //    largest entry not greater than r.  N1 is the N corresponding to the 
  //    index a.
 
  else if ( r < cdfs[ENTRIES-1] ) {         // r is also >= cdfs[0]
 
    //
    // **** This is the normal code path ****
    //
 
    int a = 0;                  // Highest value of index such that cdfs[a]
                // is known NOT to be greater than r.
    int b = ENTRIES - 1;    // Lowest value of index such that cdfs[b] is
                // known to exeed r.
 
    while (b != (a+1) ) {
      int c = (a+b+1)>>1;
      if (r > cdfs[c]) {
        a = c;
      } else {
        b = c;
      }
    }
 
    N1 = Nmin + a;
    rRange = cdfs[a+1] - cdfs[a];
    rRemainder = r - cdfs[a];
 
    //
    // **** ****
    //
 
  } // end of medium-r (normal) case
 
 
  // >>> 6 <<< 
  //    If r exceeds the greatest entry in the table for this mu, then start 
  //    from that cdf, and use the series to compute from there until r is 
  //    exceeded.  
 
  else { //       if ( r >= cdfs[ENTRIES-1] ) {
 
    // Here, division must be done explicitly, and we must also protect against
    // roundoff preventing termination.
 
    // 
    //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m/m! , m=0 to Nmax) 
    //+++ (where Nmax = mu - BELOW + ENTRIES - 1)
    //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp(-mu) mu**(Nmax)/(Nmax)!
    //+++ If the sum up to k-1 <= r < sum up to k, then N = k-1
    //+++ Consider k = Nmax in the above statement:
    //+++ If cdfs[ENTRIES-2] <= r < cdfs[ENTRIES-1], N would be Nmax-1 
    //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax
    //
 
    // Erroneous:
    //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m/m! , m=0 to Nmax-1) 
    //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp(-mu) mu**(Nmax-1)/(Nmax-1)!
    //+++ If a sum up to k-1 <= r < sum up to k, then N = k-1
    //+++ So if cdfs[ENTRIES-1] were > r, N would be Nmax-1 (or less)
    //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax
    //
 
        //  std::cout << "r = " << r << " mu = " << mu << "\n";
    long N = Nmax -1;
    double cdf = cdfs[ENTRIES-1];
    double term = cdf - cdfs[ENTRIES-2];
    double cdf0;
    while(cdf <= r) {
        N++ ;  
        //  std::cout << "  N " << N << " term " << 
        //  term << " cdf " << cdf << "\n";
        term *= ( mu / N );
    cdf0  = cdf;
        cdf  += term;
    if (cdf == cdf0) break; // If term gets so small cdf stops increasing,
                // terminate using that value of N since we
                // would never reach r.
    }
    N1 = N;
    rRange = 0;     // We can't validly omit the second true random
 
    //      N = Nmax -1;
    //      cdf = cdfs[ENTRIES-1];
    //      term = cdf - cdfs[ENTRIES-2];
    //      for (int isxz=0; isxz < 100; isxz++) {
    //          N++ ;  
    //          term *= ( mu / N );
    //      cdf0  = cdf;
    //          cdf  += term;
    //      }
    //      std::cout.precision(20);
    //      std::cout << "Final sum is " << cdf << "\n";
 
  } // end of large-r case
 
 
 
  // >>> 7 <<< 
  //    Form a second random, s, based on the position of r within the range
  //    of this table entry to the next entry.  
 
  // However, if this range is very small, then we lose too many bits of
  // randomness.  In that situation, we generate a second random for s.
 
  double s;
 
  static const double MINRANGE = .01;   // Sacrifice up to two digits of 
                        // randomness when using r to produce
                    // a second random s.  Leads to up to
                    // .09 extra randoms each time.
 
  if ( rRange > MINRANGE ) {
    //
    // **** This path taken 90% of the time ****
    //
    s = rRemainder / rRange;
  } else {
    s = e->flat();  // extra true random needed about one time in 10.
  }
 
  // >>> 8 <<< 
  //    Use the direct summation method to form a second poisson deviate N2 
  //    from deltaMu and s.
 
  N2 = 0;
  double term = std::exp(-deltaMu);
  double cdf = term;
 
  if ( s < (1 - 1.0E-10) ) {
      //
      // This is the normal path:
      //
      const double* oneOverNptr = oneOverN;
      while( cdf <= s ) {
        N2++ ;  
        oneOverNptr++;
        term *= ( deltaMu * (*oneOverNptr) );
        cdf  += term;
      }
  } else { // s is almost 1...
      while( cdf <= s ) {
        N2++ ;  
        term *= ( deltaMu / N2 );
        cdf  += term;
      }
  } // end of if ( s compared to (1 - 1.0E-10) )
 
  // >>> 9 <<< 
  //    The result is the sum of those two deviates
 
        // if (DBG_small) {
        //   std::cout << N2 << "   " << N1+N2 << "\n";
        //   DBG_small = false;
        // }
 
  return N1 + N2;
 
} // poissonDeviate()

References BELOW, ENTRIES, FIRST_MU, CLHEP::HepRandomEngine::flat(), LAST_MU, CLHEP::poissonTables, S, CLHEP::s, and CLHEP::RandPoisson::shoot().

Referenced by fire(), and shoot().

put()

std::ostream & CLHEP::RandPoissonQ::put ( std::ostream &  os ) const

virtual

Reimplemented from CLHEP::RandPoisson.

Definition at line 547 of file RandPoissonQ.cc.

                                                       {
  int pr=os.precision(20);
  std::vector<unsigned long> t(2);
  os << " " << name() << "\n";
  os << "Uvec" << "\n";
  t = DoubConv::dto2longs(a0);
  os << a0 << " " << t[0] << " " << t[1] << "\n";
  t = DoubConv::dto2longs(a1);
  os << a1 << " " << t[0] << " " << t[1] << "\n";
  t = DoubConv::dto2longs(a2);
  os << a2 << " " << t[0] << " " << t[1] << "\n";
  t = DoubConv::dto2longs(sigma);
  os << sigma << " " << t[0] << " " << t[1] << "\n";
  RandPoisson::put(os);
  os.precision(pr);
  return os;
#ifdef REMOVED
  int pr=os.precision(20);
  os << " " << name() << "\n";
  os << a0 << " " << a1 << " " << a2 << "\n";
  os << sigma << "\n";
  RandPoisson::put(os);
  os.precision(pr);
  return os;
#endif
}

References a0, a1, a2, CLHEP::DoubConv::dto2longs(), name(), CLHEP::RandPoisson::put(), and sigma.

restoreDistState()

static std::istream & CLHEP::HepRandom::restoreDistState ( std::istream &  is )

inlinestaticinherited

Definition at line 136 of file Random.h.

{return is;}

restoreEngineStatus()

void CLHEP::HepRandom::restoreEngineStatus ( const char  filename[] = "Config.conf" )

staticinherited

Definition at line 283 of file Random.cc.

{
  theDefaults().theEngine->restoreStatus( filename );
}

References CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theEngine.

restoreFullState()

std::istream & CLHEP::HepRandom::restoreFullState ( std::istream &  is )

staticinherited

Definition at line 293 of file Random.cc.

                                                          {
  is >> *getTheEngine();
  return is;
}

Referenced by CLHEP::RandFlat::restoreFullState(), and CLHEP::RandGauss::restoreFullState().

restoreStaticRandomStates()

std::istream & CLHEP::HepRandom::restoreStaticRandomStates ( std::istream &  is )

staticinherited

Definition at line 302 of file Random.cc.

                                                                   {
  return StaticRandomStates::restore(is);
}

References CLHEP::StaticRandomStates::restore().

saveDistState()

static std::ostream & CLHEP::HepRandom::saveDistState ( std::ostream &  os )

inlinestaticinherited

Definition at line 133 of file Random.h.

{return os;}

saveEngineStatus()

void CLHEP::HepRandom::saveEngineStatus ( const char  filename[] = "Config.conf" )

staticinherited

Definition at line 278 of file Random.cc.

{
  theDefaults().theEngine->saveStatus( filename );
}

References CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theEngine.

Referenced by G4CascadeInterface::ApplyYourself(), G4HadronicProcess::PostStepDoIt(), and G4CascadeInterface::Propagate().

saveFullState()

std::ostream & CLHEP::HepRandom::saveFullState ( std::ostream &  os )

staticinherited

Definition at line 288 of file Random.cc.

                                                       {
  os << *getTheEngine();
  return os;
}

Referenced by G4EventManager::DoProcessing(), G4EventManager::ProcessOneEvent(), CLHEP::RandFlat::saveFullState(), and CLHEP::RandGauss::saveFullState().

saveStaticRandomStates()

std::ostream & CLHEP::HepRandom::saveStaticRandomStates ( std::ostream &  os )

staticinherited

Definition at line 298 of file Random.cc.

                                                                {
  return StaticRandomStates::save(os);
}

References CLHEP::StaticRandomStates::save().

setOldMean()

static void CLHEP::RandPoisson::setOldMean ( double  val )

inlinestaticprotectedinherited

Definition at line 105 of file RandPoisson.h.

{oldm_st = val;}

References CLHEP::RandPoisson::oldm_st.

Referenced by CLHEP::RandPoisson::shoot().

setPStatus()

static void CLHEP::RandPoisson::setPStatus ( double  sq, double  alxm, double  g1  )

inlinestaticprotectedinherited

Definition at line 109 of file RandPoisson.h.

                                                            {
    status_st[0] = sq; status_st[1] = alxm; status_st[2] = g1;
  }

References CLHEP::RandPoisson::status_st.

Referenced by CLHEP::RandPoisson::shoot().

setTheEngine()

void CLHEP::HepRandom::setTheEngine ( HepRandomEngine *  theNewEngine )

staticinherited

Definition at line 273 of file Random.cc.

{
  theDefaults().theEngine.reset( theNewEngine, do_nothing_deleter() );
}

References CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theEngine.

Referenced by CLHEP::StaticRandomStates::restore().

setTheSeed()

void CLHEP::HepRandom::setTheSeed ( long  seed, int  lxr = 3  )

staticinherited

Definition at line 234 of file Random.cc.

{
  theDefaults().theEngine->setSeed(seed,lux);
}

References CLHEP::lux, CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theEngine.

setTheSeeds()

void CLHEP::HepRandom::setTheSeeds ( const long *  seeds, int  aux = -1  )

staticinherited

Definition at line 244 of file Random.cc.

{
  theDefaults().theEngine->setSeeds(seeds,aux);
}

References CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theEngine.

setupForDefaultMu()

void CLHEP::RandPoissonQ::setupForDefaultMu ( )

private

Definition at line 85 of file RandPoissonQ.cc.

                                     {
 
  // The following are useful for quick approximation, for large mu
  
  double sig2 = defaultMean * (.9998654 - .08346/defaultMean);
  sigma = std::sqrt(sig2);
    // sigma for the Guassian which approximates the Poisson -- naively
    // sqrt (defaultMean).
    //
    // The multiplier corrects for fact that discretization of the form
    // [gaussian+.5] increases the second moment by a small amount.
 
  double t = 1./(sig2);
 
  a2 = t/6 + t*t/324;
  a1 = std::sqrt (1-2*a2*a2*sig2);
  a0 = defaultMean + .5 - sig2 * a2;
 
  // The formula will be a0 + a1*x + a2*x*x where x has 2nd moment of sigma.
  // The coeffeicients are chosen to match the first THREE moments of the 
  // true Poisson distribution.   
  // 
  // Actually, if the correction for discretization were not needed, then 
  // a2 could be taken one order higher by adding t*t*t/5832.  However,
  // the discretization correction is not perfect, leading to inaccuracy
  // on the order to 1/mu**2, so adding a third term is overkill.  
 
} // setupForDefaultMu() 

References a0, a1, a2, CLHEP::RandPoisson::defaultMean, and sigma.

shoot() [1/2]

long CLHEP::RandPoissonQ::shoot ( double  mean = 1.0 )

static

Definition at line 119 of file RandPoissonQ.cc.

                                  {
  return shoot(getTheEngine(), xm);
}

References CLHEP::HepRandom::getTheEngine(), and shoot().

Referenced by fire(), shoot(), and shootArray().

shoot() [2/2]

long CLHEP::RandPoissonQ::shoot ( HepRandomEngine *  anEngine, double  mean = 1.0  )

static

Definition at line 143 of file RandPoissonQ.cc.

                                                               {
 
  // The following variables, static to this method, apply to the 
  // last time a large mean was supplied; they obviate certain calculations
  // if consecutive calls use the same mean.
 
  static CLHEP_THREAD_LOCAL double lastLargeMean = -1.; // Mean from previous shoot 
                    // requiring poissonDeviateQuick()
  static CLHEP_THREAD_LOCAL double lastA0;      
  static CLHEP_THREAD_LOCAL double lastA1;      
  static CLHEP_THREAD_LOCAL double lastA2;      
  static CLHEP_THREAD_LOCAL double lastSigma;       
 
  if ( mean < LAST_MU + S ) {
    return poissonDeviateSmall ( anEngine, mean );
  } else {
    if ( mean != lastLargeMean ) {
      // Compute the coefficients defining the quadratic transformation from a 
      // Gaussian to a Poisson for this mean.  Also save these for next time.
      double sig2 = mean * (.9998654 - .08346/mean);
      lastSigma = std::sqrt(sig2);
      double t = 1./sig2;
      lastA2 = t*(1./6.) + t*t*(1./324.);
      lastA1 = std::sqrt (1-2*lastA2*lastA2*sig2);
      lastA0 = mean + .5 - sig2 * lastA2;
    }
    return poissonDeviateQuick ( anEngine, lastA0, lastA1, lastA2, lastSigma );
  }
 
} // shoot (anEngine, mean)

References CLHEP_THREAD_LOCAL, LAST_MU, poissonDeviateQuick(), poissonDeviateSmall(), and S.

shootArray() [1/2]

void CLHEP::RandPoissonQ::shootArray ( const int  size, long *  vect, double  mean = 1.0  )

static

Definition at line 174 of file RandPoissonQ.cc.

                                                                  {
  for( long* v = vect; v != vect + size; ++v )
    *v = shoot(m);
     // Note: We could test for m > 100, and if it is, precompute a0, a1, a2, 
     // and sigma and call the appropriate form of poissonDeviateQuick.  
     // But since those are cached anyway, not much time would be saved.
}

References CLHEP::m, and shoot().

shootArray() [2/2]

static void CLHEP::RandPoissonQ::shootArray ( HepRandomEngine *  anEngine, const int  size, long *  vect, double  mean = 1.0  )

static

showEngineStatus()

void CLHEP::HepRandom::showEngineStatus ( )

staticinherited

Definition at line 306 of file Random.cc.

{
  theDefaults().theEngine->showStatus();
}

References CLHEP::anonymous_namespace{Random.cc}::theDefaults(), and CLHEP::anonymous_namespace{Random.cc}::defaults::theEngine.

tableBoundary()

static int CLHEP::RandPoissonQ::tableBoundary ( )

inlinestatic

Field Documentation

a0

double CLHEP::RandPoissonQ::a0

private

Definition at line 127 of file RandPoissonQ.h.

Referenced by fire(), get(), put(), and setupForDefaultMu().

a1

double CLHEP::RandPoissonQ::a1

private

Definition at line 128 of file RandPoissonQ.h.

Referenced by fire(), get(), put(), and setupForDefaultMu().

a2

double CLHEP::RandPoissonQ::a2

private

Definition at line 129 of file RandPoissonQ.h.

Referenced by fire(), get(), put(), and setupForDefaultMu().

BELOW

const int CLHEP::RandPoissonQ::BELOW = 30

staticprivate

Definition at line 145 of file RandPoissonQ.h.

Referenced by poissonDeviateSmall().

defaultMean

double CLHEP::RandPoisson::defaultMean

protectedinherited

Definition at line 99 of file RandPoisson.h.

Referenced by CLHEP::RandPoisson::fire(), fire(), CLHEP::RandPoisson::fireArray(), fireArray(), CLHEP::RandPoisson::get(), CLHEP::RandPoisson::operator()(), CLHEP::RandPoisson::put(), and setupForDefaultMu().

ENTRIES

const int CLHEP::RandPoissonQ::ENTRIES = 51

staticprivate

Definition at line 146 of file RandPoissonQ.h.

Referenced by poissonDeviateSmall().

FIRST_MU

const double CLHEP::RandPoissonQ::FIRST_MU = 10

staticprivate

Definition at line 142 of file RandPoissonQ.h.

Referenced by poissonDeviateSmall().

LAST_MU

const double CLHEP::RandPoissonQ::LAST_MU = 95

staticprivate

Definition at line 143 of file RandPoissonQ.h.

Referenced by fire(), poissonDeviateSmall(), and shoot().

localEngine

std::shared_ptr<HepRandomEngine> CLHEP::RandPoisson::localEngine

privateinherited

Definition at line 117 of file RandPoisson.h.

Referenced by CLHEP::RandPoisson::engine(), and CLHEP::RandPoisson::fire().

MAXIMUM_POISSON_DEVIATE

const double CLHEP::RandPoissonQ::MAXIMUM_POISSON_DEVIATE = 2.0E9

static

Definition at line 108 of file RandPoissonQ.h.

Referenced by poissonDeviateQuick().

meanMax

double CLHEP::RandPoisson::meanMax

protectedinherited

Definition at line 98 of file RandPoisson.h.

Referenced by CLHEP::RandPoisson::fire(), CLHEP::RandPoisson::get(), and CLHEP::RandPoisson::put().

meanMax_st

const double CLHEP::RandPoisson::meanMax_st = 2.0E9

staticprivateinherited

Definition at line 123 of file RandPoisson.h.

Referenced by CLHEP::RandPoisson::getMaxMean().

oldm

double CLHEP::RandPoisson::oldm

privateinherited

Definition at line 118 of file RandPoisson.h.

Referenced by CLHEP::RandPoisson::fire(), CLHEP::RandPoisson::get(), and CLHEP::RandPoisson::put().

oldm_st

CLHEP_THREAD_LOCAL double CLHEP::RandPoisson::oldm_st = -1.0

staticprivateinherited

Definition at line 122 of file RandPoisson.h.

Referenced by CLHEP::RandPoisson::getOldMean(), and CLHEP::RandPoisson::setOldMean().

S

const double CLHEP::RandPoissonQ::S = 5

staticprivate

Definition at line 144 of file RandPoissonQ.h.

Referenced by fire(), poissonDeviateSmall(), and shoot().

seedTable

const long CLHEP::HepRandom::seedTable

staticprotectedinherited

Definition at line 156 of file Random.h.

sigma

double CLHEP::RandPoissonQ::sigma

private

Definition at line 130 of file RandPoissonQ.h.

Referenced by fire(), get(), put(), and setupForDefaultMu().

status

double CLHEP::RandPoisson::status[3]

privateinherited

Definition at line 118 of file RandPoisson.h.

Referenced by CLHEP::RandPoisson::fire(), CLHEP::RandPoisson::get(), and CLHEP::RandPoisson::put().

status_st

CLHEP_THREAD_LOCAL double CLHEP::RandPoisson::status_st = {0., 0., 0.}

staticprivateinherited

Definition at line 121 of file RandPoisson.h.

Referenced by CLHEP::RandPoisson::getPStatus(), and CLHEP::RandPoisson::setPStatus().

---

The documentation for this class was generated from the following files:

• source/externals/clhep/include/CLHEP/Random/RandPoissonQ.h
• source/externals/clhep/src/RandPoissonQ.cc