File indexing completed on 2024-04-06 12:29:26
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),
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
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
0061 findPedestal(detId, igain, pedestals[igain], widths[igain]);
0062
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
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
0100
0101
0102 } else {
0103
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
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;
0120 }
0121 }
0122
0123 int igain = 0;
0124
0125
0126
0127
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
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
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 }