EcalTBRecProducers/src/EcalTBWeightUncalibRecHitProducer.cc

0001 /** \class EcalTBWeightUncalibRecHitProducer
0002  *   produce ECAL uncalibrated rechits from dataframes
0003  *
0004   *
0005   *  $Alex Zabi$
0006   *  Modification to detect first sample to switch gain.
0007   *  used for amplitude recontruction at high energy
0008   *  Add TDC convention option (P. Meridiani)
0009   *
0010   */
0011 #include "RecoTBCalo/EcalTBRecProducers/interface/EcalTBWeightUncalibRecHitProducer.h"
0012 #include "DataFormats/EcalDigi/interface/EcalMGPASample.h"
0013 #include "DataFormats/Common/interface/Handle.h"
0014
0015 #include <cmath>
0016
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0020 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0021
0022 #include "CondFormats/EcalObjects/interface/EcalXtalGroupId.h"
0023
0024 #include "CondFormats/EcalObjects/interface/EcalWeight.h"
0025 #include "CondFormats/EcalObjects/interface/EcalWeightSet.h"
0026
0027 #include "CondFormats/EcalObjects/interface/EcalMGPAGainRatio.h"
0028
0029 #include "CLHEP/Matrix/Matrix.h"
0030 #include "CLHEP/Matrix/SymMatrix.h"
0031 #include <vector>
0032
0033 #define DEBUG
0034 EcalTBWeightUncalibRecHitProducer::EcalTBWeightUncalibRecHitProducer(const edm::ParameterSet& ps)
0035     : ebDigiCollection_(ps.getParameter<edm::InputTag>("EBdigiCollection")),
0036       eeDigiCollection_(ps.getParameter<edm::InputTag>("EEdigiCollection")),
0037       tdcRecInfoCollection_(ps.getParameter<edm::InputTag>("tdcRecInfoCollection")),
0038       ebHitCollection_(ps.getParameter<std::string>("EBhitCollection")),
0039       eeHitCollection_(ps.getParameter<std::string>("EEhitCollection")),
0040       ebDigiToken_(consumes<EBDigiCollection>(ebDigiCollection_)),
0041       eeDigiToken_(consumes<EEDigiCollection>(eeDigiCollection_)),
0042       tbTDCRecInfoToken_(consumes<EcalTBTDCRecInfo>(tdcRecInfoCollection_)),
0043       weightXtalGroupsToken_(esConsumes()),
0044       gainRatiosToken_(esConsumes()),
0045       tbWeightsToken_(esConsumes()),
0046       pedestalsToken_(esConsumes()),
0047       testbeamEEShape(),  // Shapes have been updated in 2018 such as to be able to fetch shape from the DB if EBShape(consumesCollector())//EEShape(consumesCollector()) are used
0048       testbeamEBShape(),  // use default constructor if you would rather prefer to use Phase I hardcoded shapes (18.05.2018 K. Theofilatos)
0049       nbTimeBin_(ps.getParameter<int>("nbTimeBin")),
0050       use2004OffsetConvention_(ps.getUntrackedParameter<bool>("use2004OffsetConvention", false)) {
0051   produces<EBUncalibratedRecHitCollection>(ebHitCollection_);
0052   produces<EEUncalibratedRecHitCollection>(eeHitCollection_);
0053 }
0054
0055 EcalTBWeightUncalibRecHitProducer::~EcalTBWeightUncalibRecHitProducer() {}
0056
0057 void EcalTBWeightUncalibRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0058   using namespace edm;
0059
0060   Handle<EBDigiCollection> pEBDigis;
0061   const EBDigiCollection* EBdigis = nullptr;
0062   if (!ebDigiCollection_.label().empty() || !ebDigiCollection_.instance().empty()) {
0063     evt.getByToken(ebDigiToken_, pEBDigis);
0064     if (!pEBDigis.isValid()) {
0065       edm::LogError("EcalUncalibRecHitError") << "Error! can't get the product " << ebDigiCollection_;
0066     } else {
0067       EBdigis = pEBDigis.product();  // get a ptr to the produc
0068 #ifdef DEBUG
0069       LogDebug("EcalUncalibRecHitInfo") << "total # EBdigis: " << EBdigis->size();
0070 #endif
0071     }
0072   }
0073
0074   Handle<EEDigiCollection> pEEDigis;
0075   const EEDigiCollection* EEdigis = nullptr;
0076   if (!eeDigiCollection_.label().empty() || !eeDigiCollection_.instance().empty()) {
0077     evt.getByToken(eeDigiToken_, pEEDigis);
0078     if (!pEEDigis.isValid()) {
0079       edm::LogError("EcalUncalibRecHitError") << "Error! can't get the product " << eeDigiCollection_;
0080     } else {
0081       EEdigis = pEEDigis.product();  // get a ptr to the produc
0082 #ifdef DEBUG
0083       LogDebug("EcalUncalibRecHitInfo") << "total # EEdigis: " << EEdigis->size();
0084 #endif
0085     }
0086   }
0087
0088   if (!EBdigis && !EEdigis)
0089     return;
0090
0091   Handle<EcalTBTDCRecInfo> pRecTDC;
0092   const EcalTBTDCRecInfo* recTDC = nullptr;
0093   evt.getByToken(tbTDCRecInfoToken_, pRecTDC);
0094   if (pRecTDC.isValid()) {
0095     recTDC = pRecTDC.product();  // get a ptr to the product
0096   }
0097
0098   // fetch map of groups of xtals
0099   const auto& grp = es.getData(weightXtalGroupsToken_);
0100
0101   const EcalXtalGroupsMap& grpMap = grp.getMap();
0102
0103   // Gain Ratios
0104   const EcalGainRatioMap& gainMap = es.getData(gainRatiosToken_).getMap();  // map of gain ratios
0105
0106   // fetch TB weights
0107 #ifdef DEBUG
0108   LogDebug("EcalUncalibRecHitDebug") << "Fetching EcalTBWeights from DB ";
0109 #endif
0110   const auto& wgts = es.getData(tbWeightsToken_);
0111
0112 #ifdef DEBUG
0113   LogDebug("EcalUncalibRecHitDebug") << "EcalTBWeightMap.size(): " << std::setprecision(3) << wgts.getMap().size();
0114 #endif
0115
0116   // fetch the pedestals from the cond DB via EventSetup
0117 #ifdef DEBUG
0118   LogDebug("EcalUncalibRecHitDebug") << "fetching pedestals....";
0119 #endif
0120   const EcalPedestalsMap& pedMap = es.getData(pedestalsToken_).getMap();  // map of pedestals
0121 #ifdef DEBUG
0122   LogDebug("EcalUncalibRecHitDebug") << "done.";
0123 #endif
0124   // collection of reco'ed ampltudes to put in the event
0125
0126   auto EBuncalibRechits = std::make_unique<EBUncalibratedRecHitCollection>();
0127   auto EEuncalibRechits = std::make_unique<EEUncalibratedRecHitCollection>();
0128
0129   EcalPedestalsMapIterator pedIter;  // pedestal iterator
0130
0131   EcalGainRatioMap::const_iterator gainIter;  // gain iterator
0132
0133   EcalXtalGroupsMap::const_iterator git;  // group iterator
0134
0135   EcalTBWeights::EcalTBWeightMap::const_iterator wit;  //weights iterator
0136   // loop over EB digis
0137   //Getting the TDC bin
0138   EcalTBWeights::EcalTDCId tdcid(int(nbTimeBin_ / 2) + 1);
0139
0140   if (recTDC)
0141     if (recTDC->offset() == -999.) {
0142       edm::LogError("EcalUncalibRecHitError") << "TDC bin completely out of range. Returning";
0143       return;
0144     }
0145
0146   if (EBdigis) {
0147     for (unsigned int idig = 0; idig < EBdigis->size(); ++idig) {
0148       EBDataFrame itdg = (*EBdigis)[idig];
0149
0150       // find pedestals for this channel
0151 #ifdef DEBUG
0152       LogDebug("EcalUncalibRecHitDebug") << "looking up pedestal for crystal: " << EBDetId(itdg.id());
0153 #endif
0154       pedIter = pedMap.find(itdg.id().rawId());
0155       if (pedIter == pedMap.end()) {
0156         edm::LogError("EcalUncalibRecHitError")
0157             << "error!! could not find pedestals for channel: " << EBDetId(itdg.id())
0158             << "\n  no uncalib rechit will be made for this digi!";
0159         continue;
0160       }
0161       const EcalPedestals::Item& aped = (*pedIter);
0162       double pedVec[3];
0163       double pedRMSVec[3];
0164       pedVec[0] = aped.mean_x12;
0165       pedVec[1] = aped.mean_x6;
0166       pedVec[2] = aped.mean_x1;
0167       pedRMSVec[0] = aped.rms_x12;
0168       pedRMSVec[1] = aped.rms_x6;
0169       pedRMSVec[2] = aped.rms_x1;
0170
0171       // find gain ratios
0172 #ifdef DEBUG
0173       LogDebug("EcalUncalibRecHitDebug") << "looking up gainRatios for crystal: " << EBDetId(itdg.id());
0174 #endif
0175       gainIter = gainMap.find(itdg.id().rawId());
0176       if (gainIter == gainMap.end()) {
0177         edm::LogError("EcalUncalibRecHitError")
0178             << "error!! could not find gain ratios for channel: " << EBDetId(itdg.id())
0179             << "\n  no uncalib rechit will be made for this digi!";
0180         continue;
0181       }
0182       const EcalMGPAGainRatio& aGain = (*gainIter);
0183       double gainRatios[3];
0184       gainRatios[0] = 1.;
0185       gainRatios[1] = aGain.gain12Over6();
0186       gainRatios[2] = aGain.gain6Over1() * aGain.gain12Over6();
0187
0188       // lookup group ID for this channel
0189       git = grpMap.find(itdg.id().rawId());
0190       if (git == grpMap.end()) {
0191         edm::LogError("EcalUncalibRecHitError")
0192             << "No group id found for this crystal. something wrong with EcalWeightXtalGroups in your DB?"
0193             << "\n  no uncalib rechit will be made for digi with id: " << EBDetId(itdg.id());
0194         continue;
0195       }
0196       const EcalXtalGroupId& gid = (*git);
0197
0198       //GAIN SWITCHING DETECTION ///////////////////////////////////////////////////////////////////////////////////////////////////
0199       double sampleGainRef = 1;
0200       int sampleSwitch = 999;
0201       for (int sample = 0; sample < itdg.size(); ++sample) {
0202         double gainSample = itdg.sample(sample).gainId();
0203         if (gainSample != sampleGainRef) {
0204           sampleGainRef = gainSample;
0205           sampleSwitch = sample;
0206         }
0207       }  //loop sample
0208       ///////////////////////////////////////////////////////////////////////////////////////////////////
0209
0210       if (recTDC) {
0211         int tdcBin = 0;
0212         if (recTDC->offset() <= 0.)
0213           tdcBin = 1;
0214         else if (recTDC->offset() >= 1.)
0215           tdcBin = nbTimeBin_;
0216         else
0217           tdcBin = int(recTDC->offset() * float(nbTimeBin_)) + 1;
0218
0219         if (tdcBin < 1 || tdcBin > nbTimeBin_) {
0220           edm::LogError("EcalUncalibRecHitError")
0221               << "TDC bin out of range " << tdcBin << " offset " << recTDC->offset();
0222           continue;
0223         }
0224
0225         // In case gain switching happens at the sample 4 (5th sample)
0226         // (sample 5 (6th sample) in 2004 TDC convention) an extra
0227         // set of weights has to be used. This set of weights is assigned to
0228         // TDC values going from 25 and up.
0229         if (use2004OffsetConvention_ && sampleSwitch == 5)
0230           tdcid = EcalTBWeights::EcalTDCId(tdcBin + 25);
0231         else if (!use2004OffsetConvention_ && sampleSwitch == 4)
0232           tdcid = EcalTBWeights::EcalTDCId(tdcBin + 25);
0233         else
0234           tdcid = EcalTBWeights::EcalTDCId(tdcBin);
0235       }  //check TDC
0236
0237       // now lookup the correct weights in the map
0238       wit = wgts.getMap().find(std::make_pair(gid, tdcid));
0239       if (wit == wgts.getMap().end()) {  // no weights found for this group ID
0240         edm::LogError("EcalUncalibRecHitError")
0241             << "No weights found for EcalGroupId: " << gid.id() << " and  EcalTDCId: " << tdcid
0242             << "\n  skipping digi with id: " << EBDetId(itdg.id());
0243         continue;
0244       }
0245       const EcalWeightSet& wset = wit->second;  // this is the EcalWeightSet
0246
0247       // EcalWeightMatrix is vec<vec:double>>
0248
0249 #ifdef DEBUG
0250       LogDebug("EcalUncalibRecHitDebug") << "accessing matrices of weights...";
0251 #endif
0252       const EcalWeightSet::EcalWeightMatrix& mat1 = wset.getWeightsBeforeGainSwitch();
0253       const EcalWeightSet::EcalWeightMatrix& mat2 = wset.getWeightsAfterGainSwitch();
0254       //const EcalWeightSet::EcalChi2WeightMatrix& mat3 = wset.getChi2WeightsBeforeGainSwitch();
0255       //const EcalWeightSet::EcalChi2WeightMatrix& mat4 = wset.getChi2WeightsAfterGainSwitch();
0256       const EcalWeightSet::EcalWeightMatrix* weights[2];
0257       weights[0] = &mat1;
0258       weights[1] = &mat2;
0259       //     weights.push_back(clmat1);
0260       //     weights.push_back(clmat2);
0261       //     LogDebug("EcalUncalibRecHitDebug") << "weights before switch:\n" << clmat1 ;
0262       //     LogDebug("EcalUncalibRecHitDebug") << "weights after switch:\n" << clmat2 ;
0263
0264       // build CLHEP chi2  matrices
0265       //const EcalWeightSet::EcalChi2WeightMatrix* chi2mat[2];
0266       // chi2mat[0]=&mat3;
0267       // chi2mat[1]=&mat4;
0268
0269       EcalUncalibratedRecHit aHit = ebAlgo_.makeRecHit(itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEBShape);
0270       EBuncalibRechits->push_back(aHit);
0271 #ifdef DEBUG
0272       if (aHit.amplitude() > 0.) {
0273         LogDebug("EcalUncalibRecHitDebug") << "processed EBDataFrame with id: " << EBDetId(itdg.id()) << "\n"
0274                                            << "uncalib rechit amplitude: " << aHit.amplitude();
0275       }
0276 #endif
0277     }
0278   }
0279   // put the collection of reconstructed hits in the event
0280   evt.put(std::move(EBuncalibRechits), ebHitCollection_);
0281
0282   if (EEdigis) {
0283     for (unsigned int idig = 0; idig < EEdigis->size(); ++idig) {
0284       EEDataFrame itdg = (*EEdigis)[idig];
0285
0286       // find pedestals for this channel
0287 #ifdef DEBUG
0288       LogDebug("EcalUncalibRecHitDebug") << "looking up pedestal for crystal: " << EEDetId(itdg.id());
0289 #endif
0290       pedIter = pedMap.find(itdg.id().rawId());
0291       if (pedIter == pedMap.end()) {
0292         edm::LogError("EcalUncalibRecHitError")
0293             << "error!! could not find pedestals for channel: " << EEDetId(itdg.id())
0294             << "\n  no uncalib rechit will be made for this digi!";
0295         continue;
0296       }
0297       const EcalPedestals::Item& aped = (*pedIter);
0298       double pedVec[3];
0299       double pedRMSVec[3];
0300       pedVec[0] = aped.mean_x12;
0301       pedVec[1] = aped.mean_x6;
0302       pedVec[2] = aped.mean_x1;
0303       pedRMSVec[0] = aped.rms_x12;
0304       pedRMSVec[1] = aped.rms_x6;
0305       pedRMSVec[2] = aped.rms_x1;
0306
0307       // find gain ratios
0308 #ifdef DEBUG
0309       LogDebug("EcalUncalibRecHitDebug") << "looking up gainRatios for crystal: " << EEDetId(itdg.id());
0310 #endif
0311       gainIter = gainMap.find(itdg.id().rawId());
0312       if (gainIter == gainMap.end()) {
0313         edm::LogError("EcalUncalibRecHitError")
0314             << "error!! could not find gain ratios for channel: " << EEDetId(itdg.id())
0315             << "\n  no uncalib rechit will be made for this digi!";
0316         continue;
0317       }
0318       const EcalMGPAGainRatio& aGain = (*gainIter);
0319       double gainRatios[3];
0320       gainRatios[0] = 1.;
0321       gainRatios[1] = aGain.gain12Over6();
0322       gainRatios[2] = aGain.gain6Over1() * aGain.gain12Over6();
0323
0324       // lookup group ID for this channel
0325       git = grpMap.find(itdg.id().rawId());
0326       if (git == grpMap.end()) {
0327         edm::LogError("EcalUncalibRecHitError")
0328             << "No group id found for this crystal. something wrong with EcalWeightXtalGroups in your DB?"
0329             << "\n  no uncalib rechit will be made for digi with id: " << EEDetId(itdg.id());
0330         continue;
0331       }
0332       const EcalXtalGroupId& gid = (*git);
0333
0334       //GAIN SWITCHING DETECTION ///////////////////////////////////////////////////////////////////////////////////////////////////
0335       double sampleGainRef = 1;
0336       int sampleSwitch = 999;
0337       for (int sample = 0; sample < itdg.size(); ++sample) {
0338         double gainSample = itdg.sample(sample).gainId();
0339         if (gainSample != sampleGainRef) {
0340           sampleGainRef = gainSample;
0341           sampleSwitch = sample;
0342         }
0343       }  //loop sample
0344       ///////////////////////////////////////////////////////////////////////////////////////////////////
0345
0346       if (recTDC) {
0347         int tdcBin = 0;
0348         if (recTDC->offset() <= 0.)
0349           tdcBin = 1;
0350         else if (recTDC->offset() >= 1.)
0351           tdcBin = nbTimeBin_;
0352         else
0353           tdcBin = int(recTDC->offset() * float(nbTimeBin_)) + 1;
0354
0355         if (tdcBin < 1 || tdcBin > nbTimeBin_) {
0356           edm::LogError("EcalUncalibRecHitError")
0357               << "TDC bin out of range " << tdcBin << " offset " << recTDC->offset();
0358           continue;
0359         }
0360
0361         // In case gain switching happens at the sample 4 (5th sample)
0362         // (sample 5 (6th sample) in 2004 TDC convention) an extra
0363         // set of weights has to be used. This set of weights is assigned to
0364         // TDC values going from 25 and up.
0365         if (use2004OffsetConvention_ && sampleSwitch == 5)
0366           tdcid = EcalTBWeights::EcalTDCId(tdcBin + 25);
0367         else if (!use2004OffsetConvention_ && sampleSwitch == 4)
0368           tdcid = EcalTBWeights::EcalTDCId(tdcBin + 25);
0369         else
0370           tdcid = EcalTBWeights::EcalTDCId(tdcBin);
0371       }  //check TDC
0372
0373       // now lookup the correct weights in the map
0374       wit = wgts.getMap().find(std::make_pair(gid, tdcid));
0375       if (wit == wgts.getMap().end()) {  // no weights found for this group ID
0376         edm::LogError("EcalUncalibRecHitError")
0377             << "No weights found for EcalGroupId: " << gid.id() << " and  EcalTDCId: " << tdcid
0378             << "\n  skipping digi with id: " << EEDetId(itdg.id());
0379         continue;
0380       }
0381       const EcalWeightSet& wset = wit->second;  // this is the EcalWeightSet
0382
0383       // EcalWeightMatrix is vec<vec:double>>
0384
0385 #ifdef DEBUG
0386       LogDebug("EcalUncalibRecHitDebug") << "accessing matrices of weights...";
0387 #endif
0388       const EcalWeightSet::EcalWeightMatrix& mat1 = wset.getWeightsBeforeGainSwitch();
0389       const EcalWeightSet::EcalWeightMatrix& mat2 = wset.getWeightsAfterGainSwitch();
0390       //const EcalWeightSet::EcalChi2WeightMatrix& mat3 = wset.getChi2WeightsBeforeGainSwitch();
0391       //const EcalWeightSet::EcalChi2WeightMatrix& mat4 = wset.getChi2WeightsAfterGainSwitch();
0392       const EcalWeightSet::EcalWeightMatrix* weights[2];
0393       weights[0] = &mat1;
0394       weights[1] = &mat2;
0395       //     weights.push_back(clmat1);
0396       //     weights.push_back(clmat2);
0397       //     LogDebug("EcalUncalibRecHitDebug") << "weights before switch:\n" << clmat1 ;
0398       //     LogDebug("EcalUncalibRecHitDebug") << "weights after switch:\n" << clmat2 ;
0399
0400       // build CLHEP chi2  matrices
0401       //const EcalWeightSet::EcalChi2WeightMatrix* chi2mat[2];
0402       //chi2mat[0]=&mat3;
0403       //chi2mat[1]=&mat4;
0404
0405       EcalUncalibratedRecHit aHit = eeAlgo_.makeRecHit(itdg, pedVec, pedRMSVec, gainRatios, weights, testbeamEEShape);
0406       EEuncalibRechits->push_back(aHit);
0407 #ifdef DEBUG
0408       if (aHit.amplitude() > 0.) {
0409         LogDebug("EcalUncalibRecHitDebug") << "processed EEDataFrame with id: " << EEDetId(itdg.id()) << "\n"
0410                                            << "uncalib rechit amplitude: " << aHit.amplitude();
0411       }
0412 #endif
0413     }
0414   }
0415   // put the collection of reconstructed hits in the event
0416   evt.put(std::move(EEuncalibRechits), eeHitCollection_);
0417 }
0418
0419 // HepMatrix
0420 // EcalTBWeightUncalibRecHitProducer::makeMatrixFromVectors(const std::vector< std::vector<EcalWeight> >& vecvec) {
0421 //   int nrow = vecvec.size();
0422 //   int ncol = (vecvec[0]).size();
0423 //   HepMatrix clmat(nrow,ncol);
0424 //   //LogDebug("EcalUncalibRecHitDebug") << "created HepMatrix(" << nrow << "," << ncol << ")" ;
0425 //   for(int irow=0;irow<nrow;++irow) {
0426 //     for(int icol=0;icol<ncol;++icol) {
0427 //         clmat[irow][icol] = ((vecvec[irow])[icol]).value();
0428 //     }
0429 //   }
0430 //   return clmat;
0431 // }
0432
0433 // HepMatrix
0434 // EcalTBWeightUncalibRecHitProducer::makeDummySymMatrix(int size)
0435 // {
0436 //   HepMatrix clmat(10,10);
0437 //   //LogDebug("EcalUncalibRecHitDebug") << "created HepMatrix(" << nrow << "," << ncol << ")" ;
0438 //   for(int irow=0; irow<size; ++irow) {
0439 //     for(int icol=0 ; icol<size; ++icol) {
0440 //       clmat[irow][icol] = irow+icol;
0441 //     }
0442 //   }
0443 //   return clmat;
0444 // }