File indexing completed on 2023-03-17 11:22:02
0001
0002
0003
0004
0005
0006
0007
0008
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(),
0048 testbeamEBShape(),
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();
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();
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();
0096 }
0097
0098
0099 const auto& grp = es.getData(weightXtalGroupsToken_);
0100
0101 const EcalXtalGroupsMap& grpMap = grp.getMap();
0102
0103
0104 const EcalGainRatioMap& gainMap = es.getData(gainRatiosToken_).getMap();
0105
0106
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
0117 #ifdef DEBUG
0118 LogDebug("EcalUncalibRecHitDebug") << "fetching pedestals....";
0119 #endif
0120 const EcalPedestalsMap& pedMap = es.getData(pedestalsToken_).getMap();
0121 #ifdef DEBUG
0122 LogDebug("EcalUncalibRecHitDebug") << "done.";
0123 #endif
0124
0125
0126 auto EBuncalibRechits = std::make_unique<EBUncalibratedRecHitCollection>();
0127 auto EEuncalibRechits = std::make_unique<EEUncalibratedRecHitCollection>();
0128
0129 EcalPedestalsMapIterator pedIter;
0130
0131 EcalGainRatioMap::const_iterator gainIter;
0132
0133 EcalXtalGroupsMap::const_iterator git;
0134
0135 EcalTBWeights::EcalTBWeightMap::const_iterator wit;
0136
0137
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
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
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
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
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 }
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
0226
0227
0228
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 }
0236
0237
0238 wit = wgts.getMap().find(std::make_pair(gid, tdcid));
0239 if (wit == wgts.getMap().end()) {
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;
0246
0247
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
0255
0256 const EcalWeightSet::EcalWeightMatrix* weights[2];
0257 weights[0] = &mat1;
0258 weights[1] = &mat2;
0259
0260
0261
0262
0263
0264
0265
0266
0267
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
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
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
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
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
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 }
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
0362
0363
0364
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 }
0372
0373
0374 wit = wgts.getMap().find(std::make_pair(gid, tdcid));
0375 if (wit == wgts.getMap().end()) {
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;
0382
0383
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
0391
0392 const EcalWeightSet::EcalWeightMatrix* weights[2];
0393 weights[0] = &mat1;
0394 weights[1] = &mat2;
0395
0396
0397
0398
0399
0400
0401
0402
0403
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
0416 evt.put(std::move(EEuncalibRechits), eeHitCollection_);
0417 }
0418
0419
0420
0421
0422
0423
0424
0425
0426
0427
0428
0429
0430
0431
0432
0433
0434
0435
0436
0437
0438
0439
0440
0441
0442
0443
0444