File indexing completed on 2024-04-06 12:30:39
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
0013
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
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
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
0044 int theChannelNumber = (int)CLHEP::RandFlat::shoot(engine, NumberOfchannels);
0045
0046
0047 double noise = generate_gaussian_tail(lowLimit, noiseRMS, engine);
0048
0049
0050 theMap[theChannelNumber] = noise;
0051 }
0052 }
0053
0054
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
0061
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
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
0084 double noise = generate_gaussian_tail(lowLimit, noiseRMS, engine);
0085
0086 theVector.push_back(std::pair<int, float>(channels[i], noise));
0087 }
0088 }
0089
0090
0091
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101
0102
0103
0104
0105
0106
0107 void GaussianTailNoiseGenerator::generateRaw(float noiseRMS,
0108 std::vector<double> &theVector,
0109 CLHEP::HepRandomEngine *engine) {
0110
0111
0112
0113
0114
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
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
0145
0146
0147
0148 double s = a / sigma;
0149
0150 if (s < 1) {
0151
0152
0153
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
0164
0165
0166
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 }