Utilities/src/GaussianTail.cc

0001 #include "FastSimulation/Utilities/interface/GaussianTail.h"
0002 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0003
0004 #include <cmath>
0005 GaussianTail::GaussianTail(double sigma, double threshold) : sigma_(sigma), threshold_(threshold) {
0006   s_ = threshold_ / sigma_;
0007   ssquare_ = s_ * s_;
0008 }
0009
0010 GaussianTail::~GaussianTail() { ; }
0011
0012 double GaussianTail::shoot(RandomEngineAndDistribution const* random) const {
0013   // in the zero suppresion case, s is usually >2
0014   if (s_ > 1.) {
0015     /* Use the "supertail" deviates from the last two steps
0016        * of Marsaglia's rectangle-wedge-tail method, as described
0017        * in Knuth, v2, 3rd ed, pp 123-128.  (See also exercise 11, p139,
0018        * and the solution, p586.)
0019        */
0020
0021     double u, v, x;
0022
0023     do {
0024       u = random->flatShoot();
0025       do {
0026         v = random->flatShoot();
0027       } while (v == 0.0);
0028       x = std::sqrt(ssquare_ - 2 * std::log(v));
0029     } while (x * u > s_);
0030     return x * sigma_;
0031   } else {
0032     /* For small s, use a direct rejection method. The limit s < 1
0033          can be adjusted to optimise the overall efficiency
0034        */
0035
0036     double x;
0037
0038     do {
0039       x = random->gaussShoot();
0040     } while (x < s_);
0041     return x * sigma_;
0042   }
0043 }