GaussianTailNoiseGenerator Class Reference

#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

(

)

Definition at line 11 of file GaussianTailNoiseGenerator.cc.

References channel512_, channel768_, and 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

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

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

Definition at line 53 of file GaussianTailNoiseGenerator.cc.

References ecal_dqm_sourceclient-live_cfg::cerr, generate_gaussian_tail(), getRandomChannels(), i, HLT_25ns14e33_v1_cff::noise, and query::result.

                                                                         {
  // 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));
  }
}

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

protected

Definition at line 132 of file GaussianTailNoiseGenerator.cc.

References cmsBatch::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;
  }
}

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

Definition at line 101 of file GaussianTailNoiseGenerator.cc.

References 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);
  }
}

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

protected

Definition at line 116 of file GaussianTailNoiseGenerator.cc.

References b, channel512_, channel768_, and 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;
}

Member Data Documentation

int GaussianTailNoiseGenerator::channel512_[512]

private

Definition at line 62 of file GaussianTailNoiseGenerator.h.

Referenced by GaussianTailNoiseGenerator(), and getRandomChannels().

int GaussianTailNoiseGenerator::channel768_[768]

private

Definition at line 63 of file GaussianTailNoiseGenerator.h.

Referenced by GaussianTailNoiseGenerator(), and getRandomChannels().