#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(), funct::log(), random, s_, sigma_, mathSSE::sqrt(), ssquare_, v, and x.

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