RPDigiProducer/plugins/RPGaussianTailNoiseAdder.cc

0001 #include "SimPPS/RPDigiProducer/plugins/RPGaussianTailNoiseAdder.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include <iostream>
0004 #include "CLHEP/Random/RandGauss.h"
0005 #include "CLHEP/Random/RandFlat.h"
0006 #include "CLHEP/Random/RandBinomial.h"
0007
0008 #include <cmath>
0009
0010 using namespace std;
0011
0012 RPGaussianTailNoiseAdder::RPGaussianTailNoiseAdder(
0013     int numStrips, double theNoiseInElectrons, double theStripThresholdInE, CLHEP::HepRandomEngine &eng, int verbosity)
0014     : numStrips_(numStrips),
0015       theNoiseInElectrons(theNoiseInElectrons),
0016       theStripThresholdInE(theStripThresholdInE),
0017       rndEngine_(eng),
0018       verbosity_(verbosity) {
0019   strips_above_threshold_prob_ = std::erfc(theStripThresholdInE / sqrt(2.0) / theNoiseInElectrons) / 2;
0020 }
0021
0022 simromanpot::strip_charge_map RPGaussianTailNoiseAdder::addNoise(const simromanpot::strip_charge_map &theSignal) {
0023   simromanpot::strip_charge_map the_strip_charge_map;
0024
0025   // noise on strips with signal:
0026   for (simromanpot::strip_charge_map::const_iterator i = theSignal.begin(); i != theSignal.end(); ++i) {
0027     double noise = CLHEP::RandGauss::shoot(&(rndEngine_), 0.0, theNoiseInElectrons);
0028     the_strip_charge_map[i->first] = i->second + noise;
0029     if (verbosity_)
0030       edm::LogInfo("RPDigiProducer") << "noise added to signal strips: " << noise << "\n";
0031   }
0032
0033   // noise on the other strips
0034   int strips_no_above_threshold =
0035       CLHEP::RandBinomial::shoot(&(rndEngine_), (long)numStrips_, strips_above_threshold_prob_);
0036
0037   for (int j = 0; j < strips_no_above_threshold; j++) {
0038     int strip = CLHEP::RandFlat::shootInt(&(rndEngine_), numStrips_);
0039     if (the_strip_charge_map[strip] == 0) {
0040       the_strip_charge_map[strip] = 2 * theStripThresholdInE;
0041       //only binary decision later, no need to simulate the noise precisely,
0042       //enough to know that it exceeds the threshold
0043       if (verbosity_)
0044         edm::LogInfo("RPDigiProducer") << "nonsignal strips noise :" << strip << " " << the_strip_charge_map[strip]
0045                                        << "\n";
0046     }
0047   }
0048
0049   return the_strip_charge_map;
0050 }