EcalSimAlgos/src/ESDigitizer.cc

0001 #include "SimCalorimetry/EcalSimAlgos/interface/ESDigitizer.h"
0002 #include "SimCalorimetry/EcalSimAlgos/interface/ESElectronicsSimFast.h"
0003 #include "SimCalorimetry/EcalSimAlgos/interface/ESHitResponse.h"
0004 #include "SimDataFormats/CaloHit/interface/PCaloHit.h"
0005
0006 #include "CLHEP/Random/RandGeneral.h"
0007 #include "CLHEP/Random/RandPoissonQ.h"
0008 #include "CLHEP/Random/RandFlat.h"
0009
0010 #include <gsl/gsl_sf_erf.h>
0011 #include <gsl/gsl_sf_result.h>
0012
0013 ESDigitizer::ESDigitizer(EcalHitResponse* hitResponse, ESElectronicsSimFast* electronicsSim, bool addNoise)
0014     : EcalTDigitizer<ESDigitizerTraits>(hitResponse, electronicsSim, addNoise),
0015       m_detIds(nullptr),
0016       m_ranGeneral(nullptr),
0017       m_ESGain(0),
0018       m_histoBin(0),
0019       m_histoInf(0),
0020       m_histoWid(0),
0021       m_meanNoisy(0),
0022       m_trip() {
0023   m_trip.reserve(2500);
0024 }
0025
0026 ESDigitizer::~ESDigitizer() { delete m_ranGeneral; }
0027
0028 /// tell the digitizer which cells exist; cannot change during a run
0029 void ESDigitizer::setDetIds(const std::vector<DetId>& detIds) {
0030   assert(nullptr == m_detIds || &detIds == m_detIds);  // sanity check; don't allow to change midstream
0031   m_detIds = &detIds;
0032 }
0033
0034 void ESDigitizer::setGain(const int gain) {
0035   if (0 != m_ESGain) {
0036     assert(gain == m_ESGain);  // only allow one value
0037   } else {
0038     assert(nullptr != m_detIds && !m_detIds->empty());  // detIds must already be set as we need size
0039
0040     assert(1 == gain || 2 == gain);  // legal values
0041
0042     m_ESGain = gain;
0043
0044     if (addNoise()) {
0045       double zsThresh(0.);
0046       std::string refFile;
0047
0048       if (1 == m_ESGain) {
0049         zsThresh = 3;
0050         refFile = "SimCalorimetry/EcalSimProducers/data/esRefHistosFile_LG.txt";
0051       } else {
0052         zsThresh = 4;
0053         refFile = "SimCalorimetry/EcalSimProducers/data/esRefHistosFile_HG.txt";
0054       }
0055
0056       gsl_sf_result result;
0057       int status = gsl_sf_erf_Q_e(zsThresh, &result);
0058       if (status != 0)
0059         std::cerr << "ESDigitizer::could not compute gaussian tail probability for the threshold chosen" << std::endl;
0060
0061       const double probabilityLeft(result.val);
0062       m_meanNoisy = probabilityLeft * m_detIds->size();
0063
0064       std::ifstream histofile(edm::FileInPath(refFile).fullPath().c_str());
0065       if (!histofile.good()) {
0066         throw edm::Exception(edm::errors::InvalidReference, "NullPointer") << "Reference histos file not opened";
0067       } else {
0068         // number of bins
0069         char buffer[200];
0070         int thisLine = 0;
0071         while (0 == thisLine) {
0072           histofile.getline(buffer, 200);
0073           if (!strstr(buffer, "#") && !(strspn(buffer, " ") == strlen(buffer))) {
0074             float histoBin;
0075             sscanf(buffer, "%f", &histoBin);
0076             m_histoBin = (double)histoBin;
0077             ++thisLine;
0078           }
0079         }
0080         const uint32_t histoBin1((int)m_histoBin);
0081         const uint32_t histoBin2(histoBin1 * histoBin1);
0082
0083         double t_histoSup(0);
0084
0085         std::vector<double> t_refHistos;
0086         t_refHistos.reserve(2500);
0087
0088         int thisBin = -2;
0089         while (!(histofile.eof())) {
0090           histofile.getline(buffer, 200);
0091           if (!strstr(buffer, "#") && !(strspn(buffer, " ") == strlen(buffer))) {
0092             if (-2 == thisBin) {
0093               float histoInf;
0094               sscanf(buffer, "%f", &histoInf);
0095               m_histoInf = (double)histoInf;
0096             }
0097             if (-1 == thisBin) {
0098               float histoSup;
0099               sscanf(buffer, "%f", &histoSup);
0100               t_histoSup = (double)histoSup;
0101             }
0102             if (0 <= thisBin) {
0103               float refBin;
0104               sscanf(buffer, "%f", &refBin);
0105               if (0.5 < refBin) {
0106                 t_refHistos.push_back((double)refBin);
0107                 const uint32_t i2(thisBin / histoBin2);
0108                 const uint32_t off(i2 * histoBin2);
0109                 const uint32_t i1((thisBin - off) / histoBin1);
0110                 const uint32_t i0(thisBin - off - i1 * histoBin1);
0111                 m_trip.emplace_back(i0, i1, i2);
0112               }
0113             }
0114             ++thisBin;
0115           }
0116         }
0117         m_histoWid = (t_histoSup - m_histoInf) / m_histoBin;
0118
0119         m_histoInf -= 1000.;
0120
0121         // creating the reference distribution to extract random numbers
0122         m_ranGeneral = new CLHEP::RandGeneral(nullptr, &t_refHistos.front(), t_refHistos.size(), 0);
0123         histofile.close();
0124       }
0125     }
0126   }
0127 }
0128
0129 /// turns hits into digis
0130 void ESDigitizer::run(ESDigiCollection& output, CLHEP::HepRandomEngine* engine) {
0131   assert(nullptr != m_detIds && !m_detIds->empty() && (!addNoise() || nullptr != m_ranGeneral));  // sanity check
0132
0133   // reserve space for how many digis we expect, with some cushion
0134   output.reserve(2 * ((int)m_meanNoisy) + hitResponse()->samplesSize());
0135
0136   EcalTDigitizer<ESDigitizerTraits>::run(output, engine);
0137
0138   // random generation of channel above threshold
0139   std::vector<DetId> abThreshCh;
0140   if (addNoise())
0141     createNoisyList(abThreshCh, engine);
0142
0143   // first make a raw digi for every cell where we have noise
0144   for (std::vector<DetId>::const_iterator idItr(abThreshCh.begin()); idItr != abThreshCh.end(); ++idItr) {
0145     if (hitResponse()->findDetId(*idItr)->zero())  // only if no true hit!
0146     {
0147       ESHitResponse::ESSamples analogSignal(*idItr, 3);  // space for the noise hit
0148       uint32_t myBin((uint32_t)m_trip.size() * m_ranGeneral->shoot(engine));
0149       if (myBin == m_trip.size())
0150         --myBin;  // guard against roundup
0151       assert(myBin < m_trip.size());
0152       const Triplet& trip(m_trip[myBin]);
0153       analogSignal[0] = m_histoInf + m_histoWid * trip.first;
0154       analogSignal[1] = m_histoInf + m_histoWid * trip.second;
0155       analogSignal[2] = m_histoInf + m_histoWid * trip.third;
0156       ESDataFrame digi(*idItr);
0157       const_cast<ESElectronicsSimFast*>(elecSim())->analogToDigital(engine, analogSignal, digi, true);
0158       output.push_back(digi);
0159     }
0160   }
0161 }
0162
0163 // preparing the list of channels where the noise has to be generated
0164 void ESDigitizer::createNoisyList(std::vector<DetId>& abThreshCh, CLHEP::HepRandomEngine* engine) {
0165   CLHEP::RandPoissonQ randPoissonQ(*engine, m_meanNoisy);
0166   const unsigned int nChan(randPoissonQ.fire());
0167   abThreshCh.reserve(nChan);
0168
0169   for (unsigned int i(0); i != nChan; ++i) {
0170     std::vector<DetId>::const_iterator idItr(abThreshCh.end());
0171     uint32_t iChan(0);
0172     DetId id;
0173     do {
0174       iChan = (uint32_t)CLHEP::RandFlat::shoot(engine, m_detIds->size());
0175       if (iChan == m_detIds->size())
0176         --iChan;                         //protect against roundup at end
0177       assert(m_detIds->size() > iChan);  // sanity check
0178       id = (*m_detIds)[iChan];
0179       idItr = find(abThreshCh.begin(), abThreshCh.end(), id);
0180     } while (idItr != abThreshCh.end());
0181
0182     abThreshCh.push_back(id);
0183   }
0184 }