#include <GaussianTail.h>

Public Member Functions
	GaussianTail (const RandomEngine *engine, double sigma=1., double threshold=2.)

void	setParameters (double sigma, double threshold)

double	shoot () const

	~GaussianTail ()

Private Attributes
const RandomEngine *	random

double	s_

double	sigma_

double	ssquare_

double	threshold_

Detailed Description

Definition at line 10 of file GaussianTail.h.

Constructor & Destructor Documentation

GaussianTail::GaussianTail	(	const RandomEngine *	engine,
		double	sigma = `1.`,
		double	threshold = `2.`
	)

Definition at line 5 of file GaussianTail.cc.

References s_, sigma_, ssquare_, and threshold_.

                                               : 
   random(engine),
   sigma_(sigma),
   threshold_(threshold)
 {
   s_=threshold_/sigma_;
   ssquare_ = s_ * s_;
 }

GaussianTail::~GaussianTail ( )

Definition at line 15 of file GaussianTail.cc.

 {
   ;
 }

Member Function Documentation

void GaussianTail::setParameters	(	double	sigma,
		double	threshold
	)

inline

Definition at line 16 of file GaussianTail.h.

References s_, sigma_, ssquare_, dtDQMClient_cfg::threshold, and threshold_.

     {
       sigma_= sigma; 
       threshold_ = threshold;
       s_=threshold_/sigma_;
       ssquare_ = s_ * s_;
     };

double GaussianTail::shoot ( ) const

Definition at line 20 of file GaussianTail.cc.

References RandomEngine::flatShoot(), RandomEngine::gaussShoot(), create_public_lumi_plots::log, random, s_, sigma_, mathSSE::sqrt(), ssquare_, v, and vdt::x.

Referenced by HcalRecHitsMaker::noisifySubdet(), EcalPreshowerRecHitsMaker::noisifySubdet(), EcalBarrelRecHitsMaker::randomNoisifier(), and EcalEndcapRecHitsMaker::randomNoisifier().

 {
   // in the zero suppresion case, s is usually >2 
   if(s_>1.)
     {
       /* Use the "supertail" deviates from the last two steps
        * of Marsaglia's rectangle-wedge-tail method, as described
        * in Knuth, v2, 3rd ed, pp 123-128.  (See also exercise 11, p139,
        * and the solution, p586.)
        */
       
       double u, v, x;
       
       do
         {
           u = random->flatShoot();
           do
             {
               v = random->flatShoot();
             }
           while (v == 0.0);
           x = std::sqrt (ssquare_ - 2 * std::log (v));
         }
       while (x * u > s_);
       return x * sigma_;
     }
   else
     {
        /* For small s, use a direct rejection method. The limit s < 1
          can be adjusted to optimise the overall efficiency 
        */
       
       double x;
 
       do
         {
           x = random->gaussShoot();
         }
       while (x < s_);
       return x * sigma_; 
     }
 }