CLHEP 2.0.4.7 Reference Documentation
CLHEP Home Page CLHEP Documentation CLHEP Bug Reports |
00001 // $Id: RandPoissonQ.cc,v 1.4.4.4 2007/10/18 21:01:20 fischler Exp $ 00002 // -*- C++ -*- 00003 // 00004 // ----------------------------------------------------------------------- 00005 // HEP Random 00006 // --- RandPoissonQ --- 00007 // class implementation file 00008 // ----------------------------------------------------------------------- 00009 00010 // ======================================================================= 00011 // M. Fischler - Implemented new, much faster table-driven algorithm 00012 // applicable for mu < 100 00013 // - Implemented "quick()" methods, shich are the same as the 00014 // new methods for mu < 100 and are a skew-corrected gaussian 00015 // approximation for large mu. 00016 // M. Fischler - Removed mean=100 from the table-driven set, since it 00017 // uses a value just off the end of the table. (April 2004) 00018 // M Fischler - put and get to/from streams 12/15/04 00019 // M Fischler - Utilize RandGaussQ rather than RandGauss, as clearly 00020 // intended by the inclusion of RandGaussQ.h. Using RandGauss 00021 // introduces a subtle trap in that the state of RandPoissonQ 00022 // can never be properly captured without also saveing the 00023 // state of RandGauss! RandGaussQ is, on the other hand, 00024 // stateless except for the engine used. 00025 // M Fisculer - Modified use of wrong engine when shoot (anEngine, mean) 00026 // is called. This flaw was preventing any hope of proper 00027 // saving and restoring in the instance cases. 00028 // M Fischler - fireArray using defaultMean 2/10/05 00029 // M Fischler - put/get to/from streams uses pairs of ulongs when 00030 // + storing doubles avoid problems with precision 00031 // 4/14/05 00032 // M Fisculer - Modified use of shoot (mean) instead of 00033 // shoot(getLocalEngine(), mean) when fire(mean) is called. 00034 // This flaw was causing bad "cross-talk" between modules 00035 // in CMS, where one used its own engine, and the other 00036 // used the static generator. 10/18/07 00037 // 00038 // ======================================================================= 00039 00040 #include "CLHEP/Random/defs.h" 00041 #include "CLHEP/Random/RandPoissonQ.h" 00042 #include "CLHEP/Random/RandGaussQ.h" 00043 #include "CLHEP/Random/DoubConv.hh" 00044 #include "CLHEP/Random/Stat.h" 00045 #include <cmath> // for pow() 00046 00047 namespace CLHEP { 00048 00049 std::string RandPoissonQ::name() const {return "RandPoissonQ";} 00050 HepRandomEngine & RandPoissonQ::engine() {return RandPoisson::engine();} 00051 00052 // Initialization of static data: Note that this is all const static data, 00053 // so that saveEngineStatus properly saves all needed information. 00054 00055 // The following MUST MATCH the corresponding values used (in 00056 // poissonTables.cc) when poissonTables.cdat was created. 00057 00058 const double RandPoissonQ::FIRST_MU = 10;// lowest mu value in table 00059 const double RandPoissonQ::LAST_MU = 95;// highest mu value 00060 const double RandPoissonQ::S = 5; // Spacing between mu values 00061 const int RandPoissonQ::BELOW = 30; // Starting point for N is at mu - BELOW 00062 const int RandPoissonQ::ENTRIES = 51; // Number of entries in each mu row 00063 00064 const double RandPoissonQ::MAXIMUM_POISSON_DEVIATE = 2.0E9; 00065 // Careful -- this is NOT the maximum number that can be held in 00066 // a long. It actually should be some large number of sigma below 00067 // that. 00068 00069 // Here comes the big (9K bytes) table, kept in a file of 00070 // ENTRIES * (FIRST_MU - LAST_MU + 1)/S doubles 00071 00072 static const double poissonTables [ 51 * ( (95-10)/5 + 1 ) ] = { 00073 #include "poissonTables.cdat" 00074 }; 00075 00076 00077 // 00078 // Constructors and destructors: 00079 // 00080 00081 RandPoissonQ::RandPoissonQ(const RandPoissonQ& right) : RandPoisson(right) { 00082 a0 = right.a0; 00083 a1 = right.a1; 00084 a2 = right.a2; 00085 sigma = right.sigma; 00086 } 00087 00088 RandPoissonQ::~RandPoissonQ() { 00089 } 00090 00091 void RandPoissonQ::setupForDefaultMu() { 00092 00093 // The following are useful for quick approximation, for large mu 00094 00095 double sig2 = defaultMean * (.9998654 - .08346/defaultMean); 00096 sigma = sqrt(sig2); 00097 // sigma for the Guassian which approximates the Poisson -- naively 00098 // sqrt (defaultMean). 00099 // 00100 // The multiplier corrects for fact that discretization of the form 00101 // [gaussian+.5] increases the second moment by a small amount. 00102 00103 double t = 1./(sig2); 00104 00105 a2 = t/6 + t*t/324; 00106 a1 = sqrt (1-2*a2*a2*sig2); 00107 a0 = defaultMean + .5 - sig2 * a2; 00108 00109 // The formula will be a0 + a1*x + a2*x*x where x has 2nd moment of sigma. 00110 // The coeffeicients are chosen to match the first THREE moments of the 00111 // true Poisson distribution. 00112 // 00113 // Actually, if the correction for discretization were not needed, then 00114 // a2 could be taken one order higher by adding t*t*t/5832. However, 00115 // the discretization correction is not perfect, leading to inaccuracy 00116 // on the order to 1/mu**2, so adding a third term is overkill. 00117 00118 } // setupForDefaultMu() 00119 00120 00121 // 00122 // fire, quick, operator(), and shoot methods: 00123 // 00124 00125 long RandPoissonQ::shoot(double xm) { 00126 return shoot(getTheEngine(), xm); 00127 } 00128 00129 double RandPoissonQ::operator()() { 00130 return (double) fire(); 00131 } 00132 00133 double RandPoissonQ::operator()( double mean ) { 00134 return (double) fire(mean); 00135 } 00136 00137 long RandPoissonQ::fire(double mean) { 00138 return shoot(getLocalEngine(), mean); 00139 } 00140 00141 long RandPoissonQ::fire() { 00142 if ( defaultMean < LAST_MU + S ) { 00143 return poissonDeviateSmall ( getLocalEngine(), defaultMean ); 00144 } else { 00145 return poissonDeviateQuick ( getLocalEngine(), a0, a1, a2, sigma ); 00146 } 00147 } // fire() 00148 00149 long RandPoissonQ::shoot(HepRandomEngine* anEngine, double mean) { 00150 00151 // The following variables, static to this method, apply to the 00152 // last time a large mean was supplied; they obviate certain calculations 00153 // if consecutive calls use the same mean. 00154 00155 static double lastLargeMean = -1.; // Mean from previous shoot 00156 // requiring poissonDeviateQuick() 00157 static double lastA0; 00158 static double lastA1; 00159 static double lastA2; 00160 static double lastSigma; 00161 00162 if ( mean < LAST_MU + S ) { 00163 return poissonDeviateSmall ( anEngine, mean ); 00164 } else { 00165 if ( mean != lastLargeMean ) { 00166 // Compute the coefficients defining the quadratic transformation from a 00167 // Gaussian to a Poisson for this mean. Also save these for next time. 00168 double sig2 = mean * (.9998654 - .08346/mean); 00169 lastSigma = sqrt(sig2); 00170 double t = 1./sig2; 00171 lastA2 = t*(1./6.) + t*t*(1./324.); 00172 lastA1 = sqrt (1-2*lastA2*lastA2*sig2); 00173 lastA0 = mean + .5 - sig2 * lastA2; 00174 } 00175 return poissonDeviateQuick ( anEngine, lastA0, lastA1, lastA2, lastSigma ); 00176 } 00177 00178 } // shoot (anEngine, mean) 00179 00180 void RandPoissonQ::shootArray(const int size, long* vect, double m) { 00181 int i; 00182 for (i=0; i<size; ++i) { 00183 vect[i] = shoot(m); 00184 // Note: We could test for m > 100, and if it is, precompute a0, a1, a2, 00185 // and sigma and call the appropriate form of poissonDeviateQuick. 00186 // But since those are cached anyway, not much time would be saved. 00187 } 00188 } 00189 00190 void RandPoissonQ::fireArray(const int size, long* vect, double m) { 00191 int i; 00192 for (i=0; i<size; ++i) { 00193 vect[i] = fire( m ); 00194 } 00195 } 00196 00197 void RandPoissonQ::fireArray(const int size, long* vect) { 00198 int i; 00199 for (i=0; i<size; ++i) { 00200 vect[i] = fire( defaultMean ); 00201 } 00202 } 00203 00204 00205 // Quick Poisson deviate algorithm used by quick for large mu: 00206 00207 long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e, double mu ) { 00208 00209 // Compute the coefficients defining the quadratic transformation from a 00210 // Gaussian to a Poisson: 00211 00212 double sig2 = mu * (.9998654 - .08346/mu); 00213 double sig = sqrt(sig2); 00214 // The multiplier corrects for fact that discretization of the form 00215 // [gaussian+.5] increases the second moment by a small amount. 00216 00217 double t = 1./sig2; 00218 00219 double sa2 = t*(1./6.) + t*t*(1./324.); 00220 double sa1 = sqrt (1-2*sa2*sa2*sig2); 00221 double sa0 = mu + .5 - sig2 * sa2; 00222 00223 // The formula will be sa0 + sa1*x + sa2*x*x where x has sigma of sq. 00224 // The coeffeicients are chosen to match the first THREE moments of the 00225 // true Poisson distribution. 00226 00227 return poissonDeviateQuick ( e, sa0, sa1, sa2, sig ); 00228 } 00229 00230 00231 long RandPoissonQ::poissonDeviateQuick ( HepRandomEngine *e, 00232 double A0, double A1, double A2, double sig) { 00233 // 00234 // Quick Poisson deviate algorithm used by quick for large mu: 00235 // 00236 // The principle: For very large mu, a poisson distribution can be approximated 00237 // by a gaussian: return the integer part of mu + .5 + g where g is a unit 00238 // normal. However, this yelds a miserable approximation at values as 00239 // "large" as 100. The primary problem is that the poisson distribution is 00240 // supposed to have a skew of 1/mu**2, and the zero skew of the Guassian 00241 // leads to errors of order as big as 1/mu**2. 00242 // 00243 // We substitute for the gaussian a quadratic function of that gaussian random. 00244 // The expression looks very nearly like mu + .5 - 1/6 + g + g**2/(6*mu). 00245 // The small positive quadratic term causes the resulting variate to have 00246 // a positive skew; the -1/6 constant term is there to correct for this bias 00247 // in the mean. By adjusting these two and the linear term, we can match the 00248 // first three moments to high accuracy in 1/mu. 00249 // 00250 // The sigma used is not precisely sqrt(mu) since a rounded-off Gaussian 00251 // has a second moment which is slightly larger than that of the Gaussian. 00252 // To compensate, sig is multiplied by a factor which is slightly less than 1. 00253 00254 // double g = RandGauss::shootQuick( e ); // TEMPORARY MOD: 00255 double g = RandGaussQ::shoot( e ); // Unit normal 00256 g *= sig; 00257 double p = A2*g*g + A1*g + A0; 00258 if ( p < 0 ) return 0; // Shouldn't ever possibly happen since 00259 // mean should not be less than 100, but 00260 // we check due to paranoia. 00261 if ( p > MAXIMUM_POISSON_DEVIATE ) p = MAXIMUM_POISSON_DEVIATE; 00262 00263 return long(p); 00264 00265 } // poissonDeviateQuick () 00266 00267 00268 00269 long RandPoissonQ::poissonDeviateSmall (HepRandomEngine * e, double mean) { 00270 long N1; 00271 long N2; 00272 // The following are for later use to form a secondary random s: 00273 double rRange; // This will hold the interval between cdf for the 00274 // computed N1 and cdf for N1+1. 00275 double rRemainder = 0; // This will hold the length into that interval. 00276 00277 // Coming in, mean should not be more than LAST_MU + S. However, we will 00278 // be paranoid and test for this: 00279 00280 if ( mean > LAST_MU + S ) { 00281 return RandPoisson::shoot(e, mean); 00282 } 00283 00284 if (mean <= 0) { 00285 return 0; // Perhaps we ought to balk harder here! 00286 } 00287 00288 // >>> 1 <<< 00289 // Generate the first random, which we always will need. 00290 00291 double r = e->flat(); 00292 00293 // >>> 2 <<< 00294 // For small mean, below the start of the tables, 00295 // do the series for cdf directly. 00296 00297 // In this case, since we know the series will terminate relatively quickly, 00298 // almost alwaye use a precomputed 1/N array without fear of overrunning it. 00299 00300 static const double oneOverN[50] = 00301 { 0, 1., 1/2., 1/3., 1/4., 1/5., 1/6., 1/7., 1/8., 1/9., 00302 1/10., 1/11., 1/12., 1/13., 1/14., 1/15., 1/16., 1/17., 1/18., 1/19., 00303 1/20., 1/21., 1/22., 1/23., 1/24., 1/25., 1/26., 1/27., 1/28., 1/29., 00304 1/30., 1/31., 1/32., 1/33., 1/34., 1/35., 1/36., 1/37., 1/38., 1/39., 00305 1/40., 1/41., 1/42., 1/43., 1/44., 1/45., 1/46., 1/47., 1/48., 1/49. }; 00306 00307 00308 if ( mean < FIRST_MU ) { 00309 00310 long N = 0; 00311 double term = exp(-mean); 00312 double cdf = term; 00313 00314 if ( r < (1 - 1.0E-9) ) { 00315 // 00316 // **** This is a normal path: **** 00317 // 00318 // Except when r is very close to 1, it is certain that we will exceed r 00319 // before the 30-th term in the series, so a simple while loop is OK. 00320 const double* oneOverNptr = oneOverN; 00321 while( cdf <= r ) { 00322 N++ ; 00323 oneOverNptr++; 00324 term *= ( mean * (*oneOverNptr) ); 00325 cdf += term; 00326 } 00327 return N; 00328 // 00329 // **** **** 00330 // 00331 } else { // r is almost 1... 00332 // For r very near to 1 we would have to check that we don't fall 00333 // off the end of the table of 1/N. Since this is very rare, we just 00334 // ignore the table and do the identical while loop, using explicit 00335 // division. 00336 double cdf0; 00337 while ( cdf <= r ) { 00338 N++ ; 00339 term *= ( mean / N ); 00340 cdf0 = cdf; 00341 cdf += term; 00342 if (cdf == cdf0) break; // Can't happen, but just in case... 00343 } 00344 return N; 00345 } // end of if ( r compared to (1 - 1.0E-9) ) 00346 00347 } // End of the code for mean < FIRST_MU 00348 00349 // >>> 3 <<< 00350 // Find the row of the tables corresponding to the highest tabulated mu 00351 // which is no greater than our actual mean. 00352 00353 int rowNumber = int((mean - FIRST_MU)/S); 00354 const double * cdfs = &poissonTables [rowNumber*ENTRIES]; 00355 double mu = FIRST_MU + rowNumber*S; 00356 double deltaMu = mean - mu; 00357 int Nmin = int(mu - BELOW); 00358 if (Nmin < 1) Nmin = 1; 00359 int Nmax = Nmin + (ENTRIES - 1); 00360 00361 00362 // >>> 4 <<< 00363 // If r is less that the smallest entry in the row, then 00364 // generate the deviate directly from the series. 00365 00366 if ( r < cdfs[0] ) { 00367 00368 // In this case, we are tempted to use the actual mean, and not 00369 // generate a second deviate to account for the leftover part mean - mu. 00370 // That would be an error, generating a distribution with enough excess 00371 // at Nmin + (mean-mu)/2 to be detectable in 4,000,000 trials. 00372 00373 // Since this case is very rare (never more than .2% of the r values) 00374 // and can happen where N will be large (up to 65 for the mu=95 row) 00375 // we use explicit division so as to avoid having to worry about running 00376 // out of oneOverN table. 00377 00378 long N = 0; 00379 double term = exp(-mu); 00380 double cdf = term; 00381 double cdf0; 00382 00383 while(cdf <= r) { 00384 N++ ; 00385 term *= ( mu / N ); 00386 cdf0 = cdf; 00387 cdf += term; 00388 if (cdf == cdf0) break; // Can't happen, but just in case... 00389 } 00390 N1 = N; 00391 // std::cout << r << " " << N << " "; 00392 // DBG_small = true; 00393 rRange = 0; // In this case there is always a second r needed 00394 00395 } // end of small-r case 00396 00397 00398 // >>> 5 <<< 00399 // Assuming r lies within the scope of the row for this mu, find the 00400 // largest entry not greater than r. N1 is the N corresponding to the 00401 // index a. 00402 00403 else if ( r < cdfs[ENTRIES-1] ) { // r is also >= cdfs[0] 00404 00405 // 00406 // **** This is the normal code path **** 00407 // 00408 00409 int a = 0; // Highest value of index such that cdfs[a] 00410 // is known NOT to be greater than r. 00411 int b = ENTRIES - 1; // Lowest value of index such that cdfs[b] is 00412 // known to exeed r. 00413 00414 while (b != (a+1) ) { 00415 int c = (a+b+1)>>1; 00416 if (r > cdfs[c]) { 00417 a = c; 00418 } else { 00419 b = c; 00420 } 00421 } 00422 00423 N1 = Nmin + a; 00424 rRange = cdfs[a+1] - cdfs[a]; 00425 rRemainder = r - cdfs[a]; 00426 00427 // 00428 // **** **** 00429 // 00430 00431 } // end of medium-r (normal) case 00432 00433 00434 // >>> 6 <<< 00435 // If r exceeds the greatest entry in the table for this mu, then start 00436 // from that cdf, and use the series to compute from there until r is 00437 // exceeded. 00438 00439 else { // if ( r >= cdfs[ENTRIES-1] ) { 00440 00441 // Here, division must be done explicitly, and we must also protect against 00442 // roundoff preventing termination. 00443 00444 // 00445 //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m/m! , m=0 to Nmax) 00446 //+++ (where Nmax = mu - BELOW + ENTRIES - 1) 00447 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp(-mu) mu**(Nmax)/(Nmax)! 00448 //+++ If the sum up to k-1 <= r < sum up to k, then N = k-1 00449 //+++ Consider k = Nmax in the above statement: 00450 //+++ If cdfs[ENTRIES-2] <= r < cdfs[ENTRIES-1], N would be Nmax-1 00451 //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax 00452 // 00453 00454 // Erroneous: 00455 //+++ cdfs[ENTRIES-1] is exp(-mu) sum (mu**m/m! , m=0 to Nmax-1) 00456 //+++ cdfs[ENTRIES-1]-cdfs[ENTRIES-2] is exp(-mu) mu**(Nmax-1)/(Nmax-1)! 00457 //+++ If a sum up to k-1 <= r < sum up to k, then N = k-1 00458 //+++ So if cdfs[ENTRIES-1] were > r, N would be Nmax-1 (or less) 00459 //+++ But here r >= cdfs[ENTRIES-1] so N >= Nmax 00460 // 00461 00462 // std::cout << "r = " << r << " mu = " << mu << "\n"; 00463 long N = Nmax -1; 00464 double cdf = cdfs[ENTRIES-1]; 00465 double term = cdf - cdfs[ENTRIES-2]; 00466 double cdf0; 00467 while(cdf <= r) { 00468 N++ ; 00469 // std::cout << " N " << N << " term " << 00470 // term << " cdf " << cdf << "\n"; 00471 term *= ( mu / N ); 00472 cdf0 = cdf; 00473 cdf += term; 00474 if (cdf == cdf0) break; // If term gets so small cdf stops increasing, 00475 // terminate using that value of N since we 00476 // would never reach r. 00477 } 00478 N1 = N; 00479 rRange = 0; // We can't validly omit the second true random 00480 00481 // N = Nmax -1; 00482 // cdf = cdfs[ENTRIES-1]; 00483 // term = cdf - cdfs[ENTRIES-2]; 00484 // for (int isxz=0; isxz < 100; isxz++) { 00485 // N++ ; 00486 // term *= ( mu / N ); 00487 // cdf0 = cdf; 00488 // cdf += term; 00489 // } 00490 // std::cout.precision(20); 00491 // std::cout << "Final sum is " << cdf << "\n"; 00492 00493 } // end of large-r case 00494 00495 00496 00497 // >>> 7 <<< 00498 // Form a second random, s, based on the position of r within the range 00499 // of this table entry to the next entry. 00500 00501 // However, if this range is very small, then we lose too many bits of 00502 // randomness. In that situation, we generate a second random for s. 00503 00504 double s; 00505 00506 static const double MINRANGE = .01; // Sacrifice up to two digits of 00507 // randomness when using r to produce 00508 // a second random s. Leads to up to 00509 // .09 extra randoms each time. 00510 00511 if ( rRange > MINRANGE ) { 00512 // 00513 // **** This path taken 90% of the time **** 00514 // 00515 s = rRemainder / rRange; 00516 } else { 00517 s = e->flat(); // extra true random needed about one time in 10. 00518 } 00519 00520 // >>> 8 <<< 00521 // Use the direct summation method to form a second poisson deviate N2 00522 // from deltaMu and s. 00523 00524 N2 = 0; 00525 double term = exp(-deltaMu); 00526 double cdf = term; 00527 00528 if ( s < (1 - 1.0E-10) ) { 00529 // 00530 // This is the normal path: 00531 // 00532 const double* oneOverNptr = oneOverN; 00533 while( cdf <= s ) { 00534 N2++ ; 00535 oneOverNptr++; 00536 term *= ( deltaMu * (*oneOverNptr) ); 00537 cdf += term; 00538 } 00539 } else { // s is almost 1... 00540 while( cdf <= s ) { 00541 N2++ ; 00542 term *= ( deltaMu / N2 ); 00543 cdf += term; 00544 } 00545 } // end of if ( s compared to (1 - 1.0E-10) ) 00546 00547 // >>> 9 <<< 00548 // The result is the sum of those two deviates 00549 00550 // if (DBG_small) { 00551 // std::cout << N2 << " " << N1+N2 << "\n"; 00552 // DBG_small = false; 00553 // } 00554 00555 return N1 + N2; 00556 00557 } // poissonDeviate() 00558 00559 std::ostream & RandPoissonQ::put ( std::ostream & os ) const { 00560 int pr=os.precision(20); 00561 std::vector<unsigned long> t(2); 00562 os << " " << name() << "\n"; 00563 os << "Uvec" << "\n"; 00564 t = DoubConv::dto2longs(a0); 00565 os << a0 << " " << t[0] << " " << t[1] << "\n"; 00566 t = DoubConv::dto2longs(a1); 00567 os << a1 << " " << t[0] << " " << t[1] << "\n"; 00568 t = DoubConv::dto2longs(a2); 00569 os << a2 << " " << t[0] << " " << t[1] << "\n"; 00570 t = DoubConv::dto2longs(sigma); 00571 os << sigma << " " << t[0] << " " << t[1] << "\n"; 00572 RandPoisson::put(os); 00573 os.precision(pr); 00574 return os; 00575 #ifdef REMOVED 00576 int pr=os.precision(20); 00577 os << " " << name() << "\n"; 00578 os << a0 << " " << a1 << " " << a2 << "\n"; 00579 os << sigma << "\n"; 00580 RandPoisson::put(os); 00581 os.precision(pr); 00582 return os; 00583 #endif 00584 } 00585 00586 std::istream & RandPoissonQ::get ( std::istream & is ) { 00587 std::string inName; 00588 is >> inName; 00589 if (inName != name()) { 00590 is.clear(std::ios::badbit | is.rdstate()); 00591 std::cerr << "Mismatch when expecting to read state of a " 00592 << name() << " distribution\n" 00593 << "Name found was " << inName 00594 << "\nistream is left in the badbit state\n"; 00595 return is; 00596 } 00597 if (possibleKeywordInput(is, "Uvec", a0)) { 00598 std::vector<unsigned long> t(2); 00599 is >> a0 >> t[0] >> t[1]; a0 = DoubConv::longs2double(t); 00600 is >> a1 >> t[0] >> t[1]; a1 = DoubConv::longs2double(t); 00601 is >> a2 >> t[0] >> t[1]; a2 = DoubConv::longs2double(t); 00602 is >> sigma >> t[0] >> t[1]; sigma = DoubConv::longs2double(t); 00603 RandPoisson::get(is); 00604 return is; 00605 } 00606 // is >> a0 encompassed by possibleKeywordInput 00607 is >> a1 >> a2 >> sigma; 00608 RandPoisson::get(is); 00609 return is; 00610 } 00611 00612 } // namespace CLHEP 00613