EcalSimAlgos/src/EcalCoder.cc

0001 #include "SimCalorimetry/EcalSimAlgos/interface/EcalCoder.h"
0002 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0003 #include "SimGeneral/NoiseGenerators/interface/CorrelatedNoisifier.h"
0004 #include "DataFormats/EcalDigi/interface/EcalMGPASample.h"
0005 #include "DataFormats/EcalDigi/interface/EcalDataFrame.h"
0006 #include "CalibFormats/CaloObjects/interface/CaloTSamplesBase.icc"
0007
0008 #include <iostream>
0009
0010 //#include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0011
0012 EcalCoder::EcalCoder(bool addNoise,
0013                      bool PreMix1,
0014                      EcalCoder::Noisifier* ebCorrNoise0,
0015                      EcalCoder::Noisifier* eeCorrNoise0,
0016                      EcalCoder::Noisifier* ebCorrNoise1,
0017                      EcalCoder::Noisifier* eeCorrNoise1,
0018                      EcalCoder::Noisifier* ebCorrNoise2,
0019                      EcalCoder::Noisifier* eeCorrNoise2)
0020     : m_peds(nullptr),
0021       m_gainRatios(nullptr),
0022       m_intercals(nullptr),
0023       m_maxEneEB(1668.3),  // 4095(MAXADC)*12(gain 2)*0.035(GeVtoADC)*0.97
0024       m_maxEneEE(2859.9),  // 4095(MAXADC)*12(gain 2)*0.060(GeVtoADC)*0.97
0025       m_addNoise(addNoise),
0026       m_PreMix1(PreMix1) {
0027   m_ebCorrNoise[0] = ebCorrNoise0;
0028   assert(nullptr != m_ebCorrNoise[0]);
0029   m_eeCorrNoise[0] = eeCorrNoise0;
0030   m_ebCorrNoise[1] = ebCorrNoise1;
0031   m_eeCorrNoise[1] = eeCorrNoise1;
0032   m_ebCorrNoise[2] = ebCorrNoise2;
0033   m_eeCorrNoise[2] = eeCorrNoise2;
0034 }
0035
0036 EcalCoder::~EcalCoder() {}
0037
0038 void EcalCoder::setFullScaleEnergy(double EBscale, double EEscale) {
0039   m_maxEneEB = EBscale;
0040   m_maxEneEE = EEscale;
0041 }
0042
0043 void EcalCoder::setPedestals(const EcalPedestals* pedestals) { m_peds = pedestals; }
0044
0045 void EcalCoder::setGainRatios(const EcalGainRatios* gainRatios) { m_gainRatios = gainRatios; }
0046
0047 void EcalCoder::setIntercalibConstants(const EcalIntercalibConstantsMC* ical) { m_intercals = ical; }
0048
0049 double EcalCoder::fullScaleEnergy(const DetId& detId) const {
0050   return detId.subdetId() == EcalBarrel ? m_maxEneEB : m_maxEneEE;
0051 }
0052
0053 void EcalCoder::analogToDigital(CLHEP::HepRandomEngine* engine, const EcalSamples& clf, EcalDataFrame& df) const {
0054   df.setSize(clf.size());
0055   encode(clf, df, engine);
0056   /*   std::cout<<"   **Id=" ;
0057    if( CaloGenericDetId( clf.id() ).isEB() )
0058    {
0059       std::cout<<EBDetId( clf.id() ) ;
0060    }
0061    else
0062    {
0063       std::cout<<EEDetId( clf.id() ) ;
0064    }
0065    std::cout<<", size="<<df.size();
0066    for( int i ( 0 ) ; i != df.size() ; ++i )
0067    {
0068       std::cout<<", "<<df[i];
0069    }
0070    std::cout<<std::endl ;*/
0071 }
0072
0073 void EcalCoder::encode(const EcalSamples& ecalSamples, EcalDataFrame& df, CLHEP::HepRandomEngine* engine) const {
0074   assert(nullptr != m_peds);
0075
0076   const unsigned int csize(ecalSamples.size());
0077
0078   DetId detId = ecalSamples.id();
0079   double Emax = fullScaleEnergy(detId);
0080
0081   //....initialisation
0082   if (ecalSamples[5] > 0.)
0083     LogDebug("EcalCoder") << "Input caloSample"
0084                           << "\n"
0085                           << ecalSamples;
0086
0087   double LSB[NGAINS + 1];
0088   double pedestals[NGAINS + 1];
0089   double widths[NGAINS + 1];
0090   double gains[NGAINS + 1];
0091   int maxADC[NGAINS + 1];
0092   double trueRMS[NGAINS + 1];
0093
0094   double icalconst = 1.;
0095   findIntercalibConstant(detId, icalconst);
0096
0097   for (unsigned int igain(0); igain <= NGAINS; ++igain) {
0098     // fill in the pedestal and width
0099     findPedestal(detId, igain, pedestals[igain], widths[igain]);
0100
0101     if (0 < igain)
0102       trueRMS[igain] = std::sqrt(widths[igain] * widths[igain] - 1. / 12.);
0103
0104     // set nominal value first
0105     findGains(detId, gains);
0106
0107     LSB[igain] = 0.;
0108     if (igain > 0)
0109       LSB[igain] = Emax / (MAXADC * gains[igain]);
0110     maxADC[igain] = ADCGAINSWITCH;  // saturation at 4080 for middle and high gains x6 & x12
0111     if (igain == NGAINS)
0112       maxADC[igain] = MAXADC;  // saturation at 4095 for low gain x1
0113   }
0114
0115   CaloSamples noiseframe[] = {CaloSamples(detId, csize), CaloSamples(detId, csize), CaloSamples(detId, csize)};
0116
0117   const Noisifier* noisy[3] = {
0118       (nullptr == m_eeCorrNoise[0] || EcalBarrel == detId.subdetId() ? m_ebCorrNoise[0] : m_eeCorrNoise[0]),
0119       (EcalBarrel == detId.subdetId() ? m_ebCorrNoise[1] : m_eeCorrNoise[1]),
0120       (EcalBarrel == detId.subdetId() ? m_ebCorrNoise[2] : m_eeCorrNoise[2])};
0121
0122   if (m_addNoise) {
0123     noisy[0]->noisify(noiseframe[0], engine);  // high gain
0124     if (nullptr == noisy[1])
0125       noisy[0]->noisify(noiseframe[1], engine,
0126                         &noisy[0]->vecgau());  // med
0127     if (nullptr == noisy[2])
0128       noisy[0]->noisify(noiseframe[2], engine,
0129                         &noisy[0]->vecgau());  // low
0130   }
0131
0132   //   std::cout << " intercal, LSBs, gains " << icalconst << " " << LSB[0] << " " << LSB[1] << " " << gains[0] << " " << gains[1] << " " << Emax <<  std::endl;
0133
0134   int wait = 0;
0135   int gainId = 0;
0136   bool isSaturated = false;
0137
0138   for (unsigned int i(0); i != csize; ++i) {
0139     bool done(false);
0140     int adc = MAXADC;
0141     if (0 == wait)
0142       gainId = 1;
0143
0144     // see which gain bin it fits in
0145     int igain(gainId - 1);
0146     do {
0147       ++igain;
0148
0149       //     if( igain != gainId ) std::cout <<"$$$$ Gain switch from " << gainId
0150       //                     <<" to "<< igain << std::endl ;
0151
0152       if (1 != igain &&                     // not high gain
0153           m_addNoise &&                     // want to add noise
0154           nullptr != noisy[igain - 1] &&    // exists
0155           noiseframe[igain - 1].isBlank())  // not already done
0156       {
0157         noisy[igain - 1]->noisify(noiseframe[igain - 1], engine, &noisy[0]->vecgau());
0158         //std::cout<<"....noisifying gain level = "<<igain<<std::endl ;
0159       }
0160
0161       double signal;
0162
0163       if (!m_PreMix1) {
0164         // noiseframe filled with zeros if !m_addNoise
0165         const double asignal(pedestals[igain] + ecalSamples[i] / (LSB[igain] * icalconst) +
0166                              trueRMS[igain] * noiseframe[igain - 1][i]);
0167         signal = asignal;
0168       } else {  // Any changes made here must be reverse-engineered in EcalSignalGenerator!
0169
0170         if (igain == 1) {
0171           const double asignal(ecalSamples[i] * 1000.);  // save low level info
0172           signal = asignal;
0173         } else if (igain == 2) {
0174           const double asignal(ecalSamples[i] / (LSB[1] * icalconst));
0175           signal = asignal;
0176         } else if (igain == 3) {  // bet that no pileup hit has an energy over Emax/2
0177           const double asignal(ecalSamples[i] / (LSB[2] * icalconst));
0178           signal = asignal;
0179         } else {  //not sure we ever get here at gain=0, but hit wil be saturated anyway
0180           const double asignal(ecalSamples[i] / (LSB[3] * icalconst));  // just calculate something
0181           signal = asignal;
0182         }
0183         // old version
0184         //const double asignal ( // no pedestals for pre-mixing
0185         //           ecalSamples[i] /( LSB[igain]*icalconst ) );
0186         //signal = asignal;
0187       }
0188
0189       //     std::cout << " " << ecalSamples[i] << " " << noiseframe[igain-1][i] << std::endl;
0190
0191       const int isignal(signal);
0192       const int tmpadc(signal - (double)isignal < 0.5 ? isignal : isignal + 1);
0193       // LogDebug("EcalCoder") << "DetId " << detId.rawId() << " gain " << igain << " caloSample "
0194       //                       <<  ecalSamples[i] << " pededstal " << pedestals[igain]
0195       //                       << " noise " << widths[igain] << " conversion factor " << LSB[igain]
0196       //                       << " result (ped,tmpadc)= " << ped << " " << tmpadc;
0197
0198       if (tmpadc <= maxADC[igain]) {
0199         adc = tmpadc;
0200         done = true;
0201       }
0202     } while (!done && igain < 3);
0203
0204     if (igain == 1) {
0205       wait = 0;
0206       gainId = igain;
0207     } else {
0208       if (igain == gainId)
0209         --wait;
0210       else {
0211         wait = 5;
0212         gainId = igain;
0213       }
0214     }
0215
0216     // change the gain for saturation
0217     int storeGainId = gainId;
0218     if (gainId == 3 && adc == MAXADC) {
0219       storeGainId = 0;
0220       isSaturated = true;
0221     }
0222     // LogDebug("EcalCoder") << " Writing out frame " << i << " ADC = " << adc << " gainId = " << gainId << " storeGainId = " << storeGainId ;
0223
0224     df.setSample(i, EcalMGPASample(adc, storeGainId));
0225   }
0226   // handle saturation properly according to IN-2007/056
0227   // N.B. - FIXME
0228   // still missing the possibility for a signal to pass the saturation threshold
0229   // again within the 5 subsequent samples (higher order effect).
0230
0231   if (isSaturated) {
0232     for (unsigned int i = 0; i < ecalSamples.size(); ++i) {
0233       if (df.sample(i).gainId() == 0) {
0234         unsigned int hyst = i + 1 + 2;
0235         for (unsigned int j = i + 1; j < hyst && j < ecalSamples.size(); ++j) {
0236           df.setSample(j, EcalMGPASample(MAXADC, 0));
0237         }
0238       }
0239     }
0240   }
0241 }
0242
0243 void EcalCoder::findPedestal(const DetId& detId, int gainId, double& ped, double& width) const {
0244   /*
0245      EcalPedestalsMapIterator mapItr
0246      = m_peds->m_pedestals.find(detId.rawId());
0247      // should I care if it doesn't get found?
0248      if(mapItr == m_peds->m_pedestals.end()) {
0249      edm::LogError("EcalCoder") << "Could not find pedestal for " << detId.rawId() << " among the " << m_peds->m_pedestals.size();
0250      } else {
0251      EcalPedestals::Item item = mapItr->second;
0252    */
0253
0254   /*
0255     EcalPedestals::Item const & item = (*m_peds)(detId);
0256     ped = item.mean(gainId);
0257     width = item.rms(gainId);
0258    */
0259
0260   EcalPedestalsMap::const_iterator itped = m_peds->getMap().find(detId);
0261   ped = (*itped).mean(gainId);
0262   width = (*itped).rms(gainId);
0263
0264   if ((detId.subdetId() != EcalBarrel) && (detId.subdetId() != EcalEndcap)) {
0265     edm::LogError("EcalCoder") << "Could not find pedestal for " << detId.rawId() << " among the "
0266                                << m_peds->getMap().size();
0267   }
0268
0269   /*
0270     switch(gainId) {
0271     case 0:
0272       ped = 0.;
0273       width = 0.;
0274     case 1:
0275       ped = item.mean_x12;
0276       width = item.rms_x12;
0277       break;
0278     case 2:
0279       ped = item.mean_x6;
0280       width = item.rms_x6;
0281       break;
0282     case 3:
0283       ped = item.mean_x1;
0284       width = item.rms_x1;
0285       break;
0286     default:
0287       edm::LogError("EcalCoder") << "Bad Pedestal " << gainId;
0288       break;
0289     }
0290   */
0291
0292   LogDebug("EcalCoder") << "Pedestals for " << detId.rawId() << " gain range " << gainId << " : \n"
0293                         << "Mean = " << ped << " rms = " << width;
0294 }
0295
0296 void EcalCoder::findGains(const DetId& detId, double Gains[]) const {
0297   /*
0298     EcalGainRatios::EcalGainRatioMap::const_iterator grit=m_gainRatios->getMap().find(detId.rawId());
0299     EcalMGPAGainRatio mgpa;
0300     if( grit!=m_gainRatios->getMap().end() ){
0301     mgpa = grit->second;
0302     Gains[0] = 0.;
0303     Gains[3] = 1.;
0304     Gains[2] = mgpa.gain6Over1() ;
0305     Gains[1] = Gains[2]*(mgpa.gain12Over6()) ;
0306     LogDebug("EcalCoder") << "Gains for " << detId.rawId() << "\n" << " 1 = " << Gains[1] << "\n" << " 2 = " << Gains[2] << "\n" << " 3 = " << Gains[3];
0307     } else {
0308     edm::LogError("EcalCoder") << "Could not find gain ratios for " << detId.rawId() << " among the " << m_gainRatios->getMap().size();
0309     }
0310   */
0311
0312   EcalGainRatioMap::const_iterator grit = m_gainRatios->getMap().find(detId);
0313   Gains[0] = 0.;
0314   Gains[3] = 1.;
0315   Gains[2] = (*grit).gain6Over1();
0316   Gains[1] = Gains[2] * ((*grit).gain12Over6());
0317
0318   if ((detId.subdetId() != EcalBarrel) && (detId.subdetId() != EcalEndcap)) {
0319     edm::LogError("EcalCoder") << "Could not find gain ratios for " << detId.rawId() << " among the "
0320                                << m_gainRatios->getMap().size();
0321   }
0322 }
0323
0324 void EcalCoder::findIntercalibConstant(const DetId& detId, double& icalconst) const {
0325   EcalIntercalibConstantMC thisconst = 1.;
0326   // find intercalib constant for this xtal
0327   const EcalIntercalibConstantMCMap& icalMap = m_intercals->getMap();
0328   EcalIntercalibConstantMCMap::const_iterator icalit = icalMap.find(detId);
0329   if (icalit != icalMap.end()) {
0330     thisconst = (*icalit);
0331     if (icalconst == 0.) {
0332       thisconst = 1.;
0333     }
0334   } else {
0335     edm::LogError("EcalCoder") << "No intercalib const found for xtal " << detId.rawId()
0336                                << "! something wrong with EcalIntercalibConstants in your DB? ";
0337   }
0338   icalconst = thisconst;
0339 }