CLHEP VERSION Reference Documentation
   
CLHEP Home Page     CLHEP Documentation     CLHEP Bug Reports

RandPoissonQ.cc
Go to the documentation of this file.
1 // $Id: RandPoissonQ.cc,v 1.7 2010/06/16 17:24:53 garren Exp $
2 // -*- C++ -*-
3 //
4 // -----------------------------------------------------------------------
5 // HEP Random
6 // --- RandPoissonQ ---
7 // class implementation file
8 // -----------------------------------------------------------------------
9 
10 // =======================================================================
11 // M. Fischler - Implemented new, much faster table-driven algorithm
12 // applicable for mu < 100
13 // - Implemented "quick()" methods, shich are the same as the
14 // new methods for mu < 100 and are a skew-corrected gaussian
15 // approximation for large mu.
16 // M. Fischler - Removed mean=100 from the table-driven set, since it
17 // uses a value just off the end of the table. (April 2004)
18 // M Fischler - put and get to/from streams 12/15/04
19 // M Fischler - Utilize RandGaussQ rather than RandGauss, as clearly
20 // intended by the inclusion of RandGaussQ.h. Using RandGauss
21 // introduces a subtle trap in that the state of RandPoissonQ
22 // can never be properly captured without also saveing the
23 // state of RandGauss! RandGaussQ is, on the other hand,
24 // stateless except for the engine used.
25 // M Fisculer - Modified use of wrong engine when shoot (anEngine, mean)
26 // is called. This flaw was preventing any hope of proper
27 // saving and restoring in the instance cases.
28 // M Fischler - fireArray using defaultMean 2/10/05
29 // M Fischler - put/get to/from streams uses pairs of ulongs when
30 // + storing doubles avoid problems with precision
31 // 4/14/05
32 // M Fisculer - Modified use of shoot (mean) instead of
33 // shoot(getLocalEngine(), mean) when fire(mean) is called.
34 // This flaw was causing bad "cross-talk" between modules
35 // in CMS, where one used its own engine, and the other
36 // used the static generator. 10/18/07
37 //
38 // =======================================================================
39 
40 #include "CLHEP/Random/defs.h"
41 #include "CLHEP/Random/RandPoissonQ.h"
42 #include "CLHEP/Random/RandGaussQ.h"
43 #include "CLHEP/Random/DoubConv.hh"
44 #include "CLHEP/Random/Stat.h"
45 #include <cmath> // for std::pow()
46 
47 namespace CLHEP {
48 
49 std::string RandPoissonQ::name() const {return "RandPoissonQ";}
51 
52 // Initialization of static data: Note that this is all const static data,
53 // so that saveEngineStatus properly saves all needed information.
54 
55  // The following MUST MATCH the corresponding values used (in
56  // poissonTables.cc) when poissonTables.cdat was created.
57 
58 const double RandPoissonQ::FIRST_MU = 10;// lowest mu value in table
59 const double RandPoissonQ::LAST_MU = 95;// highest mu value
60 const double RandPoissonQ::S = 5; // Spacing between mu values
61 const int RandPoissonQ::BELOW = 30; // Starting point for N is at mu - BELOW
62 const int RandPoissonQ::ENTRIES = 51; // Number of entries in each mu row
63 
64 const double RandPoissonQ::MAXIMUM_POISSON_DEVIATE = 2.0E9;
65  // Careful -- this is NOT the maximum number that can be held in
66  // a long. It actually should be some large number of sigma below
67  // that.
68 
69  // Here comes the big (9K bytes) table, kept in a file of
70  // ENTRIES * (FIRST_MU - LAST_MU + 1)/S doubles
71 
72 static const double poissonTables [ 51 * ( (95-10)/5 + 1 ) ] = {
73 #include "poissonTables.cdat"
74 };
75 
76 
77 //
78 // Constructors and destructors:
79 //
80 
82 }
83 
84 void RandPoissonQ::setupForDefaultMu() {
85 
86  // The following are useful for quick approximation, for large mu
87 
88  double sig2 = defaultMean * (.9998654 - .08346/defaultMean);
89  sigma = std::sqrt(sig2);
90  // sigma for the Guassian which approximates the Poisson -- naively
91  // std::sqrt (defaultMean).
92  //
93  // The multiplier corrects for fact that discretization of the form
94  // [gaussian+.5] increases the second moment by a small amount.
95 
96  double t = 1./(sig2);
97 
98  a2 = t/6 + t*t/324;
99  a1 = std::sqrt (1-2*a2*a2*sig2);
100  a0 = defaultMean + .5 - sig2 * a2;
101 
102  // The formula will be a0 + a1*x + a2*x*x where x has 2nd moment of sigma.
103  // The coeffeicients are chosen to match the first THREE moments of the
104  // true Poisson distribution.
105  //
106  // Actually, if the correction for discretization were not needed, then
107  // a2 could be taken one order higher by adding t*t*t/5832. However,
108  // the discretization correction is not perfect, leading to inaccuracy
109  // on the order to 1/mu**2, so adding a third term is overkill.
110 
111 } // setupForDefaultMu()
112 
113 
114 //
115 // fire, quick, operator(), and shoot methods:
116 //
117 
118 long RandPoissonQ::shoot(double xm) {
119  return shoot(getTheEngine(), xm);
120 }
121 
123  return (double) fire();
124 }
125 
126 double RandPoissonQ::operator()( double mean ) {
127  return (double) fire(mean);
128 }
129 
130 long RandPoissonQ::fire(double mean) {
131  return shoot(getLocalEngine(), mean);
132 }
133 
135  if ( defaultMean < LAST_MU + S ) {
136  return poissonDeviateSmall ( getLocalEngine(), defaultMean );
137  } else {
138  return poissonDeviateQuick ( getLocalEngine(), a0, a1, a2, sigma );
139  }
140 } // fire()
141 
142 long RandPoissonQ::shoot(HepRandomEngine* anEngine, double mean) {
143 
144  // The following variables, static to this method, apply to the
145  // last time a large mean was supplied; they obviate certain calculations
146  // if consecutive calls use the same mean.
147 
148  static double lastLargeMean = -1.; // Mean from previous shoot
149  // requiring poissonDeviateQuick()
150  static double lastA0;
151  static double lastA1;
152  static double lastA2;
153  static double lastSigma;
154 
155  if ( mean < LAST_MU + S ) {
156  return poissonDeviateSmall ( anEngine, mean );
157  } else {
158  if ( mean != lastLargeMean ) {
159  // Compute the coefficients defining the quadratic transformation from a
160  // Gaussian to a Poisson for this mean. Also save these for next time.
161  double sig2 = mean * (.9998654 - .08346/mean);
162  lastSigma = std::sqrt(sig2);
163  double t = 1./sig2;
164  lastA2 = t*(1./6.) + t*t*(1./324.);
165  lastA1 = std::sqrt (1-2*lastA2*lastA2*sig2);
166  lastA0 = mean + .5 - sig2 * lastA2;
167  }
168  return poissonDeviateQuick ( anEngine, lastA0, lastA1, lastA2, lastSigma );
169  }
170 
171 } // shoot (anEngine, mean)
172 
173 void RandPoissonQ::shootArray(const int size, long* vect, double m) {
174  for( long* v = vect; v != vect + size; ++v )
175  *v = shoot(m);
176  // Note: We could test for m > 100, and if it is, precompute a0, a1, a2,
177  // and sigma and call the appropriate form of poissonDeviateQuick.
178  // But since those are cached anyway, not much time would be saved.
179 }
180 
181 void RandPoissonQ::fireArray(const int size, long* vect, double m) {
182  for( long* v = vect; v != vect + size; ++v )
183  *v = fire( m );
184 }
185 
186 void RandPoissonQ::fireArray(const int size, long* vect) {
187  for( long* v = vect; v != vect + size; ++v )
188  *v = fire( defaultMean );
189 }
190 
191 
192 // Quick Poisson deviate algorithm used by quick for large mu:
193 
194 long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e, double mu ) {
195 
196  // Compute the coefficients defining the quadratic transformation from a
197  // Gaussian to a Poisson:
198 
199  double sig2 = mu * (.9998654 - .08346/mu);
200  double sig = std::sqrt(sig2);
201  // The multiplier corrects for fact that discretization of the form
202  // [gaussian+.5] increases the second moment by a small amount.
203 
204  double t = 1./sig2;
205 
206  double sa2 = t*(1./6.) + t*t*(1./324.);
207  double sa1 = std::sqrt (1-2*sa2*sa2*sig2);
208  double sa0 = mu + .5 - sig2 * sa2;
209 
210  // The formula will be sa0 + sa1*x + sa2*x*x where x has sigma of sq.
211  // The coeffeicients are chosen to match the first THREE moments of the
212  // true Poisson distribution.
213 
214  return poissonDeviateQuick ( e, sa0, sa1, sa2, sig );
215 }
216 
217 
218 long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e,
219  double A0, double A1, double A2, double sig) {
220 //
221 // Quick Poisson deviate algorithm used by quick for large mu:
222 //
223 // The principle: For very large mu, a poisson distribution can be approximated
224 // by a gaussian: return the integer part of mu + .5 + g where g is a unit
225 // normal. However, this yelds a miserable approximation at values as
226 // "large" as 100. The primary problem is that the poisson distribution is
227 // supposed to have a skew of 1/mu**2, and the zero skew of the Guassian
228 // leads to errors of order as big as 1/mu**2.
229 //
230 // We substitute for the gaussian a quadratic function of that gaussian random.
231 // The expression looks very nearly like mu + .5 - 1/6 + g + g**2/(6*mu).
232 // The small positive quadratic term causes the resulting variate to have
233 // a positive skew; the -1/6 constant term is there to correct for this bias
234 // in the mean. By adjusting these two and the linear term, we can match the
235 // first three moments to high accuracy in 1/mu.
236 //
237 // The sigma used is not precisely std::sqrt(mu) since a rounded-off Gaussian
238 // has a second moment which is slightly larger than that of the Gaussian.
239 // To compensate, sig is multiplied by a factor which is slightly less than 1.
240 
241  // double g = RandGauss::shootQuick( e ); // TEMPORARY MOD:
242  double g = RandGaussQ::shoot( e ); // Unit normal
243  g *= sig;
244  double p = A2*g*g + A1*g + A0;
245  if ( p < 0 ) return 0; // Shouldn't ever possibly happen since
246  // mean should not be less than 100, but
247  // we check due to paranoia.
249 
250  return long(p);
251 
252 } // poissonDeviateQuick ()
253 
254 
255 
256 long RandPoissonQ::poissonDeviateSmall (HepRandomEngine * e, double mean) {
257  long N1;
258  long N2;
259  // The following are for later use to form a secondary random s:
260  double rRange; // This will hold the interval between cdf for the
261  // computed N1 and cdf for N1+1.
262  double rRemainder = 0; // This will hold the length into that interval.
263 
264  // Coming in, mean should not be more than LAST_MU + S. However, we will
265  // be paranoid and test for this:
266 
267  if ( mean > LAST_MU + S ) {
268  return RandPoisson::shoot(e, mean);
269  }
270 
271  if (mean <= 0) {
272  return 0; // Perhaps we ought to balk harder here!
273  }
274 
275  // >>> 1 <<<
276  // Generate the first random, which we always will need.
277 
278  double r = e->flat();
279 
280  // >>> 2 <<<
281  // For small mean, below the start of the tables,
282  // do the series for cdf directly.
283 
284  // In this case, since we know the series will terminate relatively quickly,
285  // almost alwaye use a precomputed 1/N array without fear of overrunning it.
286 
287  static const double oneOverN[50] =
288  { 0, 1., 1/2., 1/3., 1/4., 1/5., 1/6., 1/7., 1/8., 1/9.,
289  1/10., 1/11., 1/12., 1/13., 1/14., 1/15., 1/16., 1/17., 1/18., 1/19.,
290  1/20., 1/21., 1/22., 1/23., 1/24., 1/25., 1/26., 1/27., 1/28., 1/29.,
291  1/30., 1/31., 1/32., 1/33., 1/34., 1/35., 1/36., 1/37., 1/38., 1/39.,
292  1/40., 1/41., 1/42., 1/43., 1/44., 1/45., 1/46., 1/47., 1/48., 1/49. };
293 
294 
295  if ( mean < FIRST_MU ) {
296 
297  long N = 0;
298  double term = std::exp(-mean);
299  double cdf = term;
300 
301  if ( r < (1 - 1.0E-9) ) {
302  //
303  // **** This is a normal path: ****
304  //
305  // Except when r is very close to 1, it is certain that we will exceed r
306  // before the 30-th term in the series, so a simple while loop is OK.
307  const double* oneOverNptr = oneOverN;
308  while( cdf <= r ) {
309  N++ ;
310  oneOverNptr++;
311  term *= ( mean * (*oneOverNptr) );
312  cdf += term;
313  }
314  return N;
315  //
316  // **** ****
317  //
318  } else { // r is almost 1...
319  // For r very near to 1 we would have to check that we don't fall
320  // off the end of the table of 1/N. Since this is very rare, we just
321  // ignore the table and do the identical while loop, using explicit
322  // division.
323  double cdf0;
324  while ( cdf <= r ) {
325  N++ ;
326  term *= ( mean / N );
327  cdf0 = cdf;
328  cdf += term;
329  if (cdf == cdf0) break; // Can't happen, but just in case...
330  }
331  return N;
332  } // end of if ( r compared to (1 - 1.0E-9) )
333 
334  } // End of the code for mean < FIRST_MU
335 
336  // >>> 3 <<<
337  // Find the row of the tables corresponding to the highest tabulated mu
338  // which is no greater than our actual mean.
339 
340  int rowNumber = int((mean - FIRST_MU)/S);
341  const double * cdfs = &poissonTables [rowNumber*ENTRIES];
342  double mu = FIRST_MU + rowNumber*S;
343  double deltaMu = mean - mu;
344  int Nmin = int(mu - BELOW);
345  if (Nmin < 1) Nmin = 1;
346  int Nmax = Nmin + (ENTRIES - 1);
347 
348 
349  // >>> 4 <<<
350  // If r is less that the smallest entry in the row, then
351  // generate the deviate directly from the series.
352 
353  if ( r < cdfs[0] ) {
354 
355  // In this case, we are tempted to use the actual mean, and not
356  // generate a second deviate to account for the leftover part mean - mu.
357  // That would be an error, generating a distribution with enough excess
358  // at Nmin + (mean-mu)/2 to be detectable in 4,000,000 trials.
359 
360  // Since this case is very rare (never more than .2% of the r values)
361  // and can happen where N will be large (up to 65 for the mu=95 row)
362  // we use explicit division so as to avoid having to worry about running
363  // out of oneOverN table.
364 
365  long N = 0;
366  double term = std::exp(-mu);
367  double cdf = term;
368  double cdf0;
369 
370  while(cdf <= r) {
371  N++ ;
372  term *= ( mu / N );
373  cdf0 = cdf;
374  cdf += term;
375  if (cdf == cdf0) break; // Can't happen, but just in case...
376  }
377  N1 = N;
378  // std::cout << r << " " << N << " ";
379  // DBG_small = true;
380  rRange = 0; // In this case there is always a second r needed
381 
382  } // end of small-r case
383 
384 
385  // >>> 5 <<<
386  // Assuming r lies within the scope of the row for this mu, find the
387  // largest entry not greater than r. N1 is the N corresponding to the
388  // index a.
389 
390  else if ( r < cdfs[ENTRIES-1] ) { // r is also >= cdfs[0]
391 
392  //
393  // **** This is the normal code path ****
394  //
395 
396  int a = 0; // Highest value of index such that cdfs[a]
397  // is known NOT to be greater than r.
398  int b = ENTRIES - 1; // Lowest value of index such that cdfs[b] is
399  // known to exeed r.
400 
401  while (b != (a+1) ) {
402  int c = (a+b+1)>>1;
403  if (r > cdfs[c]) {
404  a = c;
405  } else {
406  b = c;
407  }
408  }
409 
410  N1 = Nmin + a;
411  rRange = cdfs[a+1] - cdfs[a];
412  rRemainder = r - cdfs[a];
413 
414  //
415  // **** ****
416  //
417 
418  } // end of medium-r (normal) case
419 
420 
421  // >>> 6 <<<
422  // If r exceeds the greatest entry in the table for this mu, then start
423  // from that cdf, and use the series to compute from there until r is
424  // exceeded.
425 
426  else { // if ( r >= cdfs[ENTRIES-1] ) {
427 
428  // Here, division must be done explicitly, and we must also protect against
429  // roundoff preventing termination.
430 
431  //
432  //+++ cdfs[ENTRIES-1] is std::exp(-mu) sum (mu**m/m! , m=0 to Nmax)
433  //+++ (where Nmax = mu - BELOW + ENTRIES - 1)
434  //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is std::exp(-mu) mu**(Nmax)/(Nmax)!
435  //+++ If the sum up to k-1 <= r < sum up to k, then N = k-1
436  //+++ Consider k = Nmax in the above statement:
437  //+++ If cdfs[ENTRIES-2] <= r < cdfs[ENTRIES-1], N would be Nmax-1
438  //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax
439  //
440 
441  // Erroneous:
442  //+++ cdfs[ENTRIES-1] is std::exp(-mu) sum (mu**m/m! , m=0 to Nmax-1)
443  //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is std::exp(-mu) mu**(Nmax-1)/(Nmax-1)!
444  //+++ If a sum up to k-1 <= r < sum up to k, then N = k-1
445  //+++ So if cdfs[ENTRIES-1] were > r, N would be Nmax-1 (or less)
446  //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax
447  //
448 
449  // std::cout << "r = " << r << " mu = " << mu << "\n";
450  long N = Nmax -1;
451  double cdf = cdfs[ENTRIES-1];
452  double term = cdf - cdfs[ENTRIES-2];
453  double cdf0;
454  while(cdf <= r) {
455  N++ ;
456  // std::cout << " N " << N << " term " <<
457  // term << " cdf " << cdf << "\n";
458  term *= ( mu / N );
459  cdf0 = cdf;
460  cdf += term;
461  if (cdf == cdf0) break; // If term gets so small cdf stops increasing,
462  // terminate using that value of N since we
463  // would never reach r.
464  }
465  N1 = N;
466  rRange = 0; // We can't validly omit the second true random
467 
468  // N = Nmax -1;
469  // cdf = cdfs[ENTRIES-1];
470  // term = cdf - cdfs[ENTRIES-2];
471  // for (int isxz=0; isxz < 100; isxz++) {
472  // N++ ;
473  // term *= ( mu / N );
474  // cdf0 = cdf;
475  // cdf += term;
476  // }
477  // std::cout.precision(20);
478  // std::cout << "Final sum is " << cdf << "\n";
479 
480  } // end of large-r case
481 
482 
483 
484  // >>> 7 <<<
485  // Form a second random, s, based on the position of r within the range
486  // of this table entry to the next entry.
487 
488  // However, if this range is very small, then we lose too many bits of
489  // randomness. In that situation, we generate a second random for s.
490 
491  double s;
492 
493  static const double MINRANGE = .01; // Sacrifice up to two digits of
494  // randomness when using r to produce
495  // a second random s. Leads to up to
496  // .09 extra randoms each time.
497 
498  if ( rRange > MINRANGE ) {
499  //
500  // **** This path taken 90% of the time ****
501  //
502  s = rRemainder / rRange;
503  } else {
504  s = e->flat(); // extra true random needed about one time in 10.
505  }
506 
507  // >>> 8 <<<
508  // Use the direct summation method to form a second poisson deviate N2
509  // from deltaMu and s.
510 
511  N2 = 0;
512  double term = std::exp(-deltaMu);
513  double cdf = term;
514 
515  if ( s < (1 - 1.0E-10) ) {
516  //
517  // This is the normal path:
518  //
519  const double* oneOverNptr = oneOverN;
520  while( cdf <= s ) {
521  N2++ ;
522  oneOverNptr++;
523  term *= ( deltaMu * (*oneOverNptr) );
524  cdf += term;
525  }
526  } else { // s is almost 1...
527  while( cdf <= s ) {
528  N2++ ;
529  term *= ( deltaMu / N2 );
530  cdf += term;
531  }
532  } // end of if ( s compared to (1 - 1.0E-10) )
533 
534  // >>> 9 <<<
535  // The result is the sum of those two deviates
536 
537  // if (DBG_small) {
538  // std::cout << N2 << " " << N1+N2 << "\n";
539  // DBG_small = false;
540  // }
541 
542  return N1 + N2;
543 
544 } // poissonDeviate()
545 
546 std::ostream & RandPoissonQ::put ( std::ostream & os ) const {
547  int pr=os.precision(20);
548  std::vector<unsigned long> t(2);
549  os << " " << name() << "\n";
550  os << "Uvec" << "\n";
551  t = DoubConv::dto2longs(a0);
552  os << a0 << " " << t[0] << " " << t[1] << "\n";
553  t = DoubConv::dto2longs(a1);
554  os << a1 << " " << t[0] << " " << t[1] << "\n";
555  t = DoubConv::dto2longs(a2);
556  os << a2 << " " << t[0] << " " << t[1] << "\n";
557  t = DoubConv::dto2longs(sigma);
558  os << sigma << " " << t[0] << " " << t[1] << "\n";
559  RandPoisson::put(os);
560  os.precision(pr);
561  return os;
562 #ifdef REMOVED
563  int pr=os.precision(20);
564  os << " " << name() << "\n";
565  os << a0 << " " << a1 << " " << a2 << "\n";
566  os << sigma << "\n";
567  RandPoisson::put(os);
568  os.precision(pr);
569  return os;
570 #endif
571 }
572 
573 std::istream & RandPoissonQ::get ( std::istream & is ) {
574  std::string inName;
575  is >> inName;
576  if (inName != name()) {
577  is.clear(std::ios::badbit | is.rdstate());
578  std::cerr << "Mismatch when expecting to read state of a "
579  << name() << " distribution\n"
580  << "Name found was " << inName
581  << "\nistream is left in the badbit state\n";
582  return is;
583  }
584  if (possibleKeywordInput(is, "Uvec", a0)) {
585  std::vector<unsigned long> t(2);
586  is >> a0 >> t[0] >> t[1]; a0 = DoubConv::longs2double(t);
587  is >> a1 >> t[0] >> t[1]; a1 = DoubConv::longs2double(t);
588  is >> a2 >> t[0] >> t[1]; a2 = DoubConv::longs2double(t);
589  is >> sigma >> t[0] >> t[1]; sigma = DoubConv::longs2double(t);
591  return is;
592  }
593  // is >> a0 encompassed by possibleKeywordInput
594  is >> a1 >> a2 >> sigma;
596  return is;
597 }
598 
599 } // namespace CLHEP
600 
CLHEP::RandPoissonQ::~RandPoissonQ
virtual ~RandPoissonQ()
Definition: RandPoissonQ.cc:81
CLHEP::RandPoissonQ::shootArray
static void shootArray(const int size, long *vect, double m=1.0)
Definition: RandPoissonQ.cc:173
CLHEP::RandPoissonQ::fire
long fire()
Definition: RandPoissonQ.cc:134
a
@ a
Definition: testCategories.cc:125
CLHEP::RandPoisson::get
std::istream & get(std::istream &is)
Definition: RandPoisson.cc:311
CLHEP::HepRandomEngine
Definition: Matrix/CLHEP/Random/RandomEngine.h:55
g
int g(shared_ptr< X >)
Definition: testSharedPtrConvertible.cc:46
CLHEP::RandPoisson::put
std::ostream & put(std::ostream &os) const
Definition: RandPoisson.cc:282
CLHEP::RandPoisson::engine
HepRandomEngine & engine()
Definition: RandPoisson.cc:37
b
@ b
Definition: testCategories.cc:125
mu
the goal is to keep the overall false rejection probability down at the to level For each validated we discuss which of course is by necessity relative timing We take the time for a single random via one of the fastest good and at any rate the ratios will vary by around depending on the processor and memory configuration used A timing for a distribution of units would mean no time used beyond the uniform random Summary Distribution Validated Validation Rejected Past N RandGauss RandGaussT RandGaussQ bins stepwise bins RandPoisson RandPoissonT mu< 100 N=50, 000, 000 ------- mu > RandPoissonQ mu< 100 N=50, 000, 000 -------(same as RandPoissonT) mu > RandGauss shoot() etc 2.5 units Validation tests applied and a very accurate series is substituted for the table method shoot() etc 1.7 units Validation tests applied and below about **a series approximation is but we have applied this test with N up to and never get anywhere near the rejection point Analytical considerations indicate that the validation test would not reject until O(10 **24) samples were inspected. ---------------------------------------------------------- 2. RandGeneral Method since we wish to have good resolution power even if just one of the mu values is it would be unwise to dial this down any further Validation for several values o mu)
Definition: validation.doc:263
is
HepRotation and so forth isNear() norm2() rectify() static Rotation row1 row4(To avoid bloat in the code pulled in for programs which don 't use all these features, we split the implementation .cc files. Only isNear() goes into the original Rotation.cc) --------------------------------------- HepAxisAngle and HepEulerAngles classes --------------------------------------- These classes are very useful and simple structures for holding the result of a nice intuituve decomposition of a rotation there is no longer much content in the distinct ZOOM PhysicsVectors library The only content left in the library is the object files representing the various Exception objects When we build the CLHEP classes for the ZOOM we will set up so as to use ZOOM SpaceVector is(but we can disable namespace usage and most of our users do so at this point). What I do is leave Hep3Vector in the global namespace
CLHEP::RandPoissonQ::fireArray
void fireArray(const int size, long *vect)
Definition: RandPoissonQ.cc:186
CLHEP::RandGaussQ::shoot
static double shoot()
CLHEP::DoubConv::longs2double
static double longs2double(const std::vector< unsigned long > &v)
Definition: DoubConv.cc:106
CLHEP::RandPoissonQ::shoot
static long shoot(double m=1.0)
Definition: RandPoissonQ.cc:118
size
user code seldom needs to call this function directly ZMerrno whether or not they are still recorded ZMerrno size() Return the(integer) number of ZMthrow 'n exceptions currently recorded. 5) ZMerrno.clear() Set an internal counter to zero. This counter is available(see next function) to user code to track ZMthrow 'n exceptions that have occurred during any arbitrary time interval. 6) ZMerrno.countSinceCleared() Return the(integer) number of ZMthrow 'n exceptions that have been recorded via ZMerrno.write()
CLHEP::RandPoissonQ::engine
HepRandomEngine & engine()
Definition: RandPoissonQ.cc:50
N
the goal is to keep the overall false rejection probability down at the to level For each validated we discuss which of course is by necessity relative timing We take the time for a single random via one of the fastest good and at any rate the ratios will vary by around depending on the processor and memory configuration used A timing for a distribution of units would mean no time used beyond the uniform random Summary Distribution Validated Validation Rejected Past N RandGauss N
Definition: validation.doc:48
CLHEP
Definition: ClhepVersion.h:13
CLHEP::RandPoissonQ::MAXIMUM_POISSON_DEVIATE
static const double MAXIMUM_POISSON_DEVIATE
Definition: Matrix/CLHEP/Random/RandPoissonQ.h:110
v
they are gone ZOOM Features Discontinued The following features of the ZOOM package were felt to be extreme overkill These have been after checking that no existing user code was utilizing as in SpaceVector v
Definition: keyMergeIssues.doc:324
CLHEP::RandPoisson::defaultMean
double defaultMean
Definition: Matrix/CLHEP/Random/RandPoisson.h:100
CLHEP::RandPoissonQ::name
std::string name() const
Definition: RandPoissonQ.cc:49
CLHEP::possibleKeywordInput
bool possibleKeywordInput(IS &is, const std::string &key, T &t)
Definition: Matrix/CLHEP/Random/RandomEngine.h:168
CLHEP::RandPoissonQ::put
std::ostream & put(std::ostream &os) const
Definition: RandPoissonQ.cc:546
CLHEP::RandPoisson::getLocalEngine
HepRandomEngine * getLocalEngine()
s
Methods applicble to containers of as in std::list< LorentzVector > s
Definition: keyMergeIssues.doc:328
CLHEP::DoubConv::dto2longs
static std::vector< unsigned long > dto2longs(double d)
Definition: DoubConv.cc:90
CLHEP::HepRandom::getTheEngine
static HepRandomEngine * getTheEngine()
Definition: Random.cc:166
CLHEP::RandPoissonQ::get
std::istream & get(std::istream &is)
Definition: RandPoissonQ.cc:573
CLHEP::RandPoisson::shoot
static long shoot(double m=1.0)
Definition: RandPoisson.cc:92
CLHEP::RandPoissonQ::operator()
double operator()()
Definition: RandPoissonQ.cc:122