NoiseGenerators/src/GaussianTailNoiseGenerator.cc

0001 #include "CLHEP/Random/RandFlat.h"
0002 #include "CLHEP/Random/RandGaussQ.h"
0003 #include "CLHEP/Random/RandPoissonQ.h"
0004 #include "SimGeneral/NoiseGenerators/interface/GaussianTailNoiseGenerator.h"
0005
0006 #include <cmath>
0007
0008 #include <gsl/gsl_sf_erf.h>
0009 #include <gsl/gsl_sf_result.h>
0010
0011 GaussianTailNoiseGenerator::GaussianTailNoiseGenerator() {
0012   // we have two cases: 512 and 768 channels
0013   // other cases are not allowed so far (performances issue)
0014   for (unsigned int i = 0; i < 512; ++i)
0015     channel512_[i] = i;
0016   for (unsigned int i = 0; i < 768; ++i)
0017     channel768_[i] = i;
0018 }
0019
0020 // this version is used by pixel
0021 void GaussianTailNoiseGenerator::generate(int NumberOfchannels,
0022                                           float threshold,
0023                                           float noiseRMS,
0024                                           std::map<int, float, std::less<int>> &theMap,
0025                                           CLHEP::HepRandomEngine *engine) {
0026   // Gaussian tail probability
0027   gsl_sf_result result;
0028   int status = gsl_sf_erf_Q_e(threshold, &result);
0029
0030   if (status != 0)
0031     std::cerr << "GaussianTailNoiseGenerator::could not compute gaussian tail "
0032                  "probability for the threshold chosen"
0033               << std::endl;
0034
0035   float probabilityLeft = result.val;
0036   float meanNumberOfNoisyChannels = probabilityLeft * NumberOfchannels;
0037
0038   CLHEP::RandPoissonQ randPoissonQ(*engine, meanNumberOfNoisyChannels);
0039   int numberOfNoisyChannels = randPoissonQ.fire();
0040
0041   float lowLimit = threshold * noiseRMS;
0042   for (int i = 0; i < numberOfNoisyChannels; i++) {
0043     // Find a random channel number
0044     int theChannelNumber = (int)CLHEP::RandFlat::shoot(engine, NumberOfchannels);
0045
0046     // Find random noise value
0047     double noise = generate_gaussian_tail(lowLimit, noiseRMS, engine);
0048
0049     // Fill in map
0050     theMap[theChannelNumber] = noise;
0051   }
0052 }
0053
0054 // this version is used by strips
0055 void GaussianTailNoiseGenerator::generate(int NumberOfchannels,
0056                                           float threshold,
0057                                           float noiseRMS,
0058                                           std::vector<std::pair<int, float>> &theVector,
0059                                           CLHEP::HepRandomEngine *engine) {
0060   // Compute number of channels with noise above threshold
0061   // Gaussian tail probability
0062   gsl_sf_result result;
0063   int status = gsl_sf_erf_Q_e(threshold, &result);
0064   if (status != 0)
0065     std::cerr << "GaussianTailNoiseGenerator::could not compute gaussian tail "
0066                  "probability for the threshold chosen"
0067               << std::endl;
0068   double probabilityLeft = result.val;
0069   double meanNumberOfNoisyChannels = probabilityLeft * NumberOfchannels;
0070
0071   CLHEP::RandPoissonQ randPoissonQ(*engine, meanNumberOfNoisyChannels);
0072   int numberOfNoisyChannels = randPoissonQ.fire();
0073
0074   if (numberOfNoisyChannels > NumberOfchannels)
0075     numberOfNoisyChannels = NumberOfchannels;
0076
0077   // Compute the list of noisy channels
0078   theVector.reserve(numberOfNoisyChannels);
0079   float lowLimit = threshold * noiseRMS;
0080   int *channels = getRandomChannels(numberOfNoisyChannels, NumberOfchannels, engine);
0081
0082   for (int i = 0; i < numberOfNoisyChannels; i++) {
0083     // Find random noise value
0084     double noise = generate_gaussian_tail(lowLimit, noiseRMS, engine);
0085     // Fill in the vector
0086     theVector.push_back(std::pair<int, float>(channels[i], noise));
0087   }
0088 }
0089
0090 /*
0091 // used by strips in VR mode
0092 void GaussianTailNoiseGenerator::generateRaw(int NumberOfchannels,
0093                                              float noiseRMS,
0094                                              std::vector<std::pair<int,float> >
0095 &theVector, CLHEP::HepRandomEngine* engine ) {
0096   theVector.reserve(NumberOfchannels);
0097   for (int i = 0; i < NumberOfchannels; i++) {
0098     // Find random noise value
0099     float noise = CLHEP::RandGaussQ::shoot(engine, 0., noiseRMS);
0100     // Fill in the vector
0101     theVector.push_back(std::pair<int, float>(i,noise));
0102   }
0103 }
0104 */
0105
0106 // used by strips in VR mode
0107 void GaussianTailNoiseGenerator::generateRaw(float noiseRMS,
0108                                              std::vector<double> &theVector,
0109                                              CLHEP::HepRandomEngine *engine) {
0110   // it was shown that a complex approach, inspired from the ZS case,
0111   // does not allow to gain much.
0112   // A cut at 2 sigmas only saves 25% of the processing time, while the cost
0113   // in terms of meaning is huge.
0114   // We therefore use here the trivial approach (as in the early 2XX cycle)
0115   unsigned int numberOfchannels = theVector.size();
0116   for (unsigned int i = 0; i < numberOfchannels; ++i) {
0117     if (theVector[i] == 0)
0118       theVector[i] = CLHEP::RandGaussQ::shoot(engine, 0., noiseRMS);
0119   }
0120 }
0121
0122 int *GaussianTailNoiseGenerator::getRandomChannels(int numberOfNoisyChannels,
0123                                                    int numberOfchannels,
0124                                                    CLHEP::HepRandomEngine *engine) {
0125   if (numberOfNoisyChannels > numberOfchannels)
0126     numberOfNoisyChannels = numberOfchannels;
0127   int *array = channel512_;
0128   if (numberOfchannels == 768)
0129     array = channel768_;
0130   int theChannelNumber;
0131   for (int j = 0; j < numberOfNoisyChannels; ++j) {
0132     theChannelNumber = (int)CLHEP::RandFlat::shoot(engine, numberOfchannels - j) + j;
0133     // swap the two array elements... this is optimized by the compiler
0134     int b = array[j];
0135     array[j] = array[theChannelNumber];
0136     array[theChannelNumber] = b;
0137   }
0138   return array;
0139 }
0140
0141 double GaussianTailNoiseGenerator::generate_gaussian_tail(const double a,
0142                                                           const double sigma,
0143                                                           CLHEP::HepRandomEngine *engine) {
0144   /* Returns a gaussian random variable larger than a
0145    * This implementation does one-sided upper-tailed deviates.
0146    */
0147
0148   double s = a / sigma;
0149
0150   if (s < 1) {
0151     /*
0152       For small s, use a direct rejection method. The limit s < 1
0153       can be adjusted to optimise the overall efficiency
0154     */
0155     double x;
0156
0157     do {
0158       x = CLHEP::RandGaussQ::shoot(engine, 0., 1.0);
0159     } while (x < s);
0160     return x * sigma;
0161
0162   } else {
0163     /* Use the "supertail" deviates from the last two steps
0164      * of Marsaglia's rectangle-wedge-tail method, as described
0165      * in Knuth, v2, 3rd ed, pp 123-128.  (See also exercise 11, p139,
0166      * and the solution, p586.)
0167      */
0168
0169     double u, v, x;
0170
0171     do {
0172       u = CLHEP::RandFlat::shoot(engine);
0173       do {
0174         v = CLHEP::RandFlat::shoot(engine);
0175       } while (v == 0.0);
0176       x = sqrt(s * s - 2 * log(v));
0177     } while (x * u > s);
0178     return x * sigma;
0179   }
0180 }