#include <GaussianTailNoiseGenerator.h>

Public Member Functions
	GaussianTailNoiseGenerator ()

void	generate (int NumberOfchannels, float threshold, float noiseRMS, std::map< int, float > &theMap, CLHEP::HepRandomEngine *)

void	generate (int NumberOfchannels, float threshold, float noiseRMS, std::vector< std::pair< int, float >> &, CLHEP::HepRandomEngine *)

void	generateRaw (float noiseRMS, std::vector< double > &, CLHEP::HepRandomEngine *)

Protected Member Functions
double	generate_gaussian_tail (const double, const double, CLHEP::HepRandomEngine *)

int *	getRandomChannels (int, int, CLHEP::HepRandomEngine *)

Private Attributes
int	channel512_ [512]

int	channel768_ [768]

Detailed Description

Generation of random noise on "numberOfChannels" channels with a given threshold. The generated noise :

corresponds to a Gaussian distribution of RMS = "noiseRMS"
is larger than threshold*noiseRMS.

Initial author : Veronique Lefebure 08.10.98
according to the FORTRAN code tgreset.F from Pascal Vanlaer
Modified by C. Delaere 01.10.09

Fills in a map < channel number, generated noise >

Definition at line 24 of file GaussianTailNoiseGenerator.h.

Constructor & Destructor Documentation

◆ GaussianTailNoiseGenerator()

GaussianTailNoiseGenerator::GaussianTailNoiseGenerator ( )

Definition at line 11 of file GaussianTailNoiseGenerator.cc.

References channel512_, channel768_, and mps_fire::i.

                                                        {
   // we have two cases: 512 and 768 channels
   // other cases are not allowed so far (performances issue)
   for (unsigned int i = 0; i < 512; ++i)
     channel512_[i] = i;
   for (unsigned int i = 0; i < 768; ++i)
     channel768_[i] = i;
 }

Member Function Documentation

◆ generate() [1/2]

void GaussianTailNoiseGenerator::generate	(	int	NumberOfchannels,
		float	threshold,
		float	noiseRMS,
		std::map< int, float > &	theMap,
		CLHEP::HepRandomEngine *
	)

◆ generate() [2/2]

void GaussianTailNoiseGenerator::generate	(	int	NumberOfchannels,
		float	threshold,
		float	noiseRMS,
		std::vector< std::pair< int, float >> &	theVector,
		CLHEP::HepRandomEngine *	engine
	)

Definition at line 55 of file GaussianTailNoiseGenerator.cc.

References DMR_cfg::cerr, ewkTauDQM_cfi::channels, generate_gaussian_tail(), getRandomChannels(), mps_fire::i, hgchebackDigitizer_cfi::noise, mps_fire::result, mps_update::status, and DiMuonV_cfg::threshold.

                                                                         {
   // Compute number of channels with noise above threshold
   // Gaussian tail probability
   gsl_sf_result result;
   int status = gsl_sf_erf_Q_e(threshold, &result);
   if (status != 0)
     std::cerr << "GaussianTailNoiseGenerator::could not compute gaussian tail "
                  "probability for the threshold chosen"
               << std::endl;
   double probabilityLeft = result.val;
   double meanNumberOfNoisyChannels = probabilityLeft * NumberOfchannels;
 
   CLHEP::RandPoissonQ randPoissonQ(*engine, meanNumberOfNoisyChannels);
   int numberOfNoisyChannels = randPoissonQ.fire();
 
   if (numberOfNoisyChannels > NumberOfchannels)
     numberOfNoisyChannels = NumberOfchannels;
 
   // Compute the list of noisy channels
   theVector.reserve(numberOfNoisyChannels);
   float lowLimit = threshold * noiseRMS;
   int *channels = getRandomChannels(numberOfNoisyChannels, NumberOfchannels, engine);
 
   for (int i = 0; i < numberOfNoisyChannels; i++) {
     // Find random noise value
     double noise = generate_gaussian_tail(lowLimit, noiseRMS, engine);
     // Fill in the vector
     theVector.push_back(std::pair<int, float>(channels[i], noise));
   }
 }

◆ generate_gaussian_tail()

double GaussianTailNoiseGenerator::generate_gaussian_tail	(	const double	a,
		const double	sigma,
		CLHEP::HepRandomEngine *	engine
	)

protected

Definition at line 141 of file GaussianTailNoiseGenerator.cc.

References a, dqm-mbProfile::log, alignCSCRings::s, mathSSE::sqrt(), findQualityFiles::v, and x.

Referenced by generate().

                                                                                         {
   /* Returns a gaussian random variable larger than a
    * This implementation does one-sided upper-tailed deviates.
    */
 
   double s = a / sigma;
 
   if (s < 1) {
     /*
       For small s, use a direct rejection method. The limit s < 1
       can be adjusted to optimise the overall efficiency
     */
     double x;
 
     do {
       x = CLHEP::RandGaussQ::shoot(engine, 0., 1.0);
     } while (x < s);
     return x * sigma;
 
   } else {
     /* 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 = CLHEP::RandFlat::shoot(engine);
       do {
         v = CLHEP::RandFlat::shoot(engine);
       } while (v == 0.0);
       x = sqrt(s * s - 2 * log(v));
     } while (x * u > s);
     return x * sigma;
   }
 }

◆ generateRaw()

void GaussianTailNoiseGenerator::generateRaw	(	float	noiseRMS,
		std::vector< double > &	theVector,
		CLHEP::HepRandomEngine *	engine
	)

Definition at line 107 of file GaussianTailNoiseGenerator.cc.

References mps_fire::i.

                                                                            {
   // it was shown that a complex approach, inspired from the ZS case,
   // does not allow to gain much.
   // A cut at 2 sigmas only saves 25% of the processing time, while the cost
   // in terms of meaning is huge.
   // We therefore use here the trivial approach (as in the early 2XX cycle)
   unsigned int numberOfchannels = theVector.size();
   for (unsigned int i = 0; i < numberOfchannels; ++i) {
     if (theVector[i] == 0)
       theVector[i] = CLHEP::RandGaussQ::shoot(engine, 0., noiseRMS);
   }
 }

◆ getRandomChannels()

int * GaussianTailNoiseGenerator::getRandomChannels	(	int	numberOfNoisyChannels,
		int	numberOfchannels,
		CLHEP::HepRandomEngine *	engine
	)

protected

Definition at line 122 of file GaussianTailNoiseGenerator.cc.

References mps_check::array, b, channel512_, channel768_, createfilelist::int, and dqmiolumiharvest::j.

Referenced by generate().

                                                                                  {
   if (numberOfNoisyChannels > numberOfchannels)
     numberOfNoisyChannels = numberOfchannels;
   int *array = channel512_;
   if (numberOfchannels == 768)
     array = channel768_;
   int theChannelNumber;
   for (int j = 0; j < numberOfNoisyChannels; ++j) {
     theChannelNumber = (int)CLHEP::RandFlat::shoot(engine, numberOfchannels - j) + j;
     // swap the two array elements... this is optimized by the compiler
     int b = array[j];
     array[j] = array[theChannelNumber];
     array[theChannelNumber] = b;
   }
   return array;
 }