#include <GaussianTail.h>

Public Member Functions
	GaussianTail (double sigma=1., double threshold=2.)

void	setParameters (double sigma, double threshold)

double	shoot (RandomEngineAndDistribution const *) const

	~GaussianTail ()

Private Attributes
double	s_

double	sigma_

double	ssquare_

double	threshold_

Detailed Description

Definition at line 10 of file GaussianTail.h.

Constructor & Destructor Documentation

GaussianTail::GaussianTail	(	double	sigma = `1.`,
		double	threshold = `2.`
	)

Definition at line 5 of file GaussianTail.cc.

References s_, sigma_, ssquare_, and threshold_.

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

GaussianTail::~GaussianTail ( )

Definition at line 13 of file GaussianTail.cc.

 {
   ;
 }

Member Function Documentation

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

inline

Definition at line 15 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 ( RandomEngineAndDistribution const * random ) const

Definition at line 18 of file GaussianTail.cc.

References RandomEngineAndDistribution::flatShoot(), RandomEngineAndDistribution::gaussShoot(), log, s_, sigma_, mathSSE::sqrt(), ssquare_, findQualityFiles::v, and x.

Referenced by EcalPreshowerRecHitsMaker::noisifySubdet(), HcalRecHitsMaker::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_; 
     }
 }