EcalSimAlgos/src/EcalLiteDTUCoder.cc

0001 #include "SimCalorimetry/EcalSimAlgos/interface/EcalLiteDTUCoder.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "SimGeneral/NoiseGenerators/interface/CorrelatedNoisifier.h"
0004 #include "DataFormats/EcalDigi/interface/EcalLiteDTUSample.h"
0005 #include "DataFormats/EcalDigi/interface/EcalDataFrame_Ph2.h"
0006
0007 #include <iostream>
0008
0009 EcalLiteDTUCoder::EcalLiteDTUCoder(bool addNoise,
0010                                    bool PreMix1,
0011                                    EcalLiteDTUCoder::Noisifier* ebCorrNoise0,
0012                                    EcalLiteDTUCoder::Noisifier* ebCorrNoise1)
0013     : m_peds(nullptr),
0014       m_gainRatios(0),
0015       m_intercals(nullptr),
0016       m_maxEneEB(ecalPh2::maxEneEB),  // Maximum for CATIA: LSB gain 10: 0.048 MeV
0017       m_addNoise(addNoise),
0018       m_PreMix1(PreMix1),
0019       m_ebCorrNoise{ebCorrNoise0, ebCorrNoise1}
0020
0021 {}
0022
0023 EcalLiteDTUCoder::~EcalLiteDTUCoder() {}
0024
0025 void EcalLiteDTUCoder::setFullScaleEnergy(double EBscale) { m_maxEneEB = ecalPh2::maxEneEB; }
0026
0027 void EcalLiteDTUCoder::setPedestals(const EcalLiteDTUPedestalsMap* pedestals) { m_peds = pedestals; }
0028
0029 void EcalLiteDTUCoder::setGainRatios(float gainRatios) { m_gainRatios = gainRatios; }
0030
0031 void EcalLiteDTUCoder::setIntercalibConstants(const EcalIntercalibConstantsMC* ical) { m_intercals = ical; }
0032
0033 double EcalLiteDTUCoder::fullScaleEnergy(const DetId& detId) const { return m_maxEneEB; }
0034
0035 void EcalLiteDTUCoder::analogToDigital(CLHEP::HepRandomEngine* engine,
0036                                        const EcalSamples& clf,
0037                                        EcalDataFrame_Ph2& df) const {
0038   df.setSize(clf.size());
0039   encode(clf, df, engine);
0040 }
0041
0042 void EcalLiteDTUCoder::encode(const EcalSamples& ecalSamples,
0043                               EcalDataFrame_Ph2& df,
0044                               CLHEP::HepRandomEngine* engine) const {
0045   const int nSamples(ecalSamples.size());
0046
0047   DetId detId = ecalSamples.id();
0048   double Emax = fullScaleEnergy(detId);
0049
0050   //N Gains set to 2 in the .h
0051   double pedestals[ecalPh2::NGAINS];
0052   double widths[ecalPh2::NGAINS];
0053   double LSB[ecalPh2::NGAINS];
0054   double trueRMS[ecalPh2::NGAINS];
0055   int nSaturatedSamples = 0;
0056   double icalconst = 1.;
0057   findIntercalibConstant(detId, icalconst);
0058
0059   for (unsigned int igain(0); igain < ecalPh2::NGAINS; ++igain) {
0060     // fill in the pedestal and width
0061     findPedestal(detId, igain, pedestals[igain], widths[igain]);
0062     // insert an absolute value in the trueRMS
0063     trueRMS[igain] = std::sqrt(std::abs(widths[igain] * widths[igain] - 1. / 12.));
0064
0065     LSB[igain] = Emax / (ecalPh2::MAXADC * ecalPh2::gains[igain]);
0066   }
0067
0068   CaloSamples noiseframe[ecalPh2::NGAINS] = {
0069       CaloSamples(detId, nSamples),
0070       CaloSamples(detId, nSamples),
0071   };
0072
0073   const Noisifier* noisy[ecalPh2::NGAINS] = {m_ebCorrNoise[0], m_ebCorrNoise[1]};
0074
0075   if (m_addNoise) {
0076     for (unsigned int ig = 0; ig < ecalPh2::NGAINS; ++ig) {
0077       noisy[ig]->noisify(noiseframe[ig], engine);
0078     }
0079   }
0080
0081   std::vector<std::vector<int>> adctrace(nSamples);
0082   int firstSaturatedSample[ecalPh2::NGAINS] = {0, 0};
0083
0084   for (int i(0); i != nSamples; ++i)
0085     adctrace[i].resize(ecalPh2::NGAINS);
0086
0087   for (unsigned int igain = 0; igain < ecalPh2::NGAINS; ++igain) {
0088     for (int i(0); i != nSamples; ++i) {
0089       adctrace[i][igain] = -1;
0090     }
0091   }
0092
0093   // fill ADC trace in gain 0 (x10) and gain 1 (x1)
0094   for (unsigned int igain = 0; igain < ecalPh2::NGAINS; ++igain) {
0095     for (int i(0); i != nSamples; ++i) {
0096       double asignal = 0;
0097       if (!m_PreMix1) {
0098         asignal = pedestals[igain] + ecalSamples[i] / (LSB[igain] * icalconst) + trueRMS[igain] * noiseframe[igain][i];
0099         //Analog signal value for each sample in ADC.
0100         //It is corrected by the intercalibration constants
0101
0102       } else {
0103         //  no noise nor pedestal when premixing
0104         asignal = ecalSamples[i] / (LSB[igain] * icalconst);
0105       }
0106       int isignal = asignal;
0107
0108       unsigned int adc = asignal - (double)isignal < 0.5 ? isignal : isignal + 1;
0109       // gain 0 (x10) channel is saturated, readout will use gain 1 (x10), but I count the number of saturated samples
0110       if (adc > ecalPh2::MAXADC) {
0111         adc = ecalPh2::MAXADC;
0112         if (nSaturatedSamples == 0)
0113           firstSaturatedSample[igain] = i;
0114         nSaturatedSamples++;
0115       }
0116       adctrace[i][igain] = adc;
0117     }
0118     if (nSaturatedSamples == 0) {
0119       break;  //  gain 0 (x10) is not saturated, so don't bother with gain 1
0120     }
0121   }  // for igain
0122
0123   int igain = 0;
0124
0125   //Limits of gain 1:
0126   //The Lite DTU sends 5 samples before the saturating one, and 10 after with gain 1.
0127   //we put the maximum in bin 5, but could happen that the system saturates before.
0128
0129   int previousSaturatedSamples = 5;
0130   int nextSaturatedSamples = 10;
0131   int startingLowerGainSample = 0;
0132   int endingLowerGainSample = (firstSaturatedSample[0] + nextSaturatedSamples + (nSaturatedSamples));
0133
0134   if (nSaturatedSamples != 0 and (firstSaturatedSample[0] - previousSaturatedSamples) < 0) {
0135     startingLowerGainSample = 0;
0136   } else {
0137     startingLowerGainSample = (firstSaturatedSample[0] - previousSaturatedSamples);
0138   }
0139
0140   //Setting values to the samples:
0141   for (int j = 0; j < nSamples; ++j) {
0142     if (nSaturatedSamples != 0 and j >= startingLowerGainSample and j < endingLowerGainSample) {
0143       igain = 1;
0144     } else {
0145       igain = 0;
0146     }
0147     df.setSample(j, EcalLiteDTUSample(adctrace[j][igain], igain));
0148   }
0149 }
0150
0151 void EcalLiteDTUCoder::findPedestal(const DetId& detId, int gainId, double& ped, double& width) const {
0152   EcalLiteDTUPedestalsMap::const_iterator itped = m_peds->getMap().find(detId);
0153   if (itped != m_peds->getMap().end()) {
0154     ped = (*itped).mean(gainId);
0155     width = (*itped).rms(gainId);
0156     LogDebug("EcalLiteDTUCoder") << "Pedestals for " << detId.rawId() << " gain range " << gainId << " : \n"
0157                                  << "Mean = " << ped << " rms = " << width;
0158   } else {
0159     LogDebug("EcalLiteDTUCoder") << "Pedestals not found, put default values (ped: 12; width: 2.5) \n";
0160     ped = 12.;
0161     width = 2.5;
0162   }
0163 }
0164
0165 void EcalLiteDTUCoder::findIntercalibConstant(const DetId& detId, double& icalconst) const {
0166   EcalIntercalibConstantMC thisconst = 1.;
0167   // find intercalib constant for this xtal
0168   const EcalIntercalibConstantMCMap& icalMap = m_intercals->getMap();
0169   EcalIntercalibConstantMCMap::const_iterator icalit = icalMap.find(detId);
0170   if (icalit != icalMap.end()) {
0171     thisconst = (*icalit);
0172   } else {
0173     LogDebug("EcalLiteDTUCoder") << "Intercalib Constant not found, put default value \n";
0174     thisconst = 1.;
0175   }
0176
0177   if (icalconst == 0.)
0178     thisconst = 1.;
0179
0180   icalconst = thisconst;
0181 }