File indexing completed on 2024-09-10 02:59:09
0001
0002
0003
0004
0005
0006
0007 #include "CondFormats/EcalObjects/interface/EcalMGPAGainRatio.h"
0008 #include "DataFormats/Common/interface/Handle.h"
0009 #include "FWCore/Framework/interface/ESHandle.h"
0010 #include "FWCore/Framework/interface/EventSetup.h"
0011 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0012 #include <cmath>
0013 #include <map>
0014 #include <memory>
0015
0016
0017
0018 #include "DataMixingEMDigiWorker.h"
0019
0020 using namespace std;
0021
0022 namespace edm {
0023
0024
0025
0026 DataMixingEMDigiWorker::DataMixingEMDigiWorker() {}
0027
0028
0029 DataMixingEMDigiWorker::DataMixingEMDigiWorker(const edm::ParameterSet &ps, edm::ConsumesCollector &&iC)
0030 : label_(ps.getParameter<std::string>("Label"))
0031
0032 {
0033
0034
0035
0036
0037
0038
0039 EBProducerSig_ = ps.getParameter<edm::InputTag>("EBdigiProducerSig");
0040 EEProducerSig_ = ps.getParameter<edm::InputTag>("EEdigiProducerSig");
0041 ESProducerSig_ = ps.getParameter<edm::InputTag>("ESdigiProducerSig");
0042
0043 EBDigiToken_ = iC.consumes<EBDigiCollection>(EBProducerSig_);
0044 EEDigiToken_ = iC.consumes<EEDigiCollection>(EEProducerSig_);
0045 ESDigiToken_ = iC.consumes<ESDigiCollection>(ESProducerSig_);
0046
0047 EBPileInputTag_ = ps.getParameter<edm::InputTag>("EBPileInputTag");
0048 EEPileInputTag_ = ps.getParameter<edm::InputTag>("EEPileInputTag");
0049 ESPileInputTag_ = ps.getParameter<edm::InputTag>("ESPileInputTag");
0050
0051 EBDigiPileToken_ = iC.consumes<EBDigiCollection>(EBPileInputTag_);
0052 EEDigiPileToken_ = iC.consumes<EEDigiCollection>(EEPileInputTag_);
0053 ESDigiPileToken_ = iC.consumes<ESDigiCollection>(ESPileInputTag_);
0054
0055 EBDigiCollectionDM_ = ps.getParameter<std::string>("EBDigiCollectionDM");
0056 EEDigiCollectionDM_ = ps.getParameter<std::string>("EEDigiCollectionDM");
0057 ESDigiCollectionDM_ = ps.getParameter<std::string>("ESDigiCollectionDM");
0058
0059 pedToken_ = iC.esConsumes<EcalPedestals, EcalPedestalsRcd>();
0060 grToken_ = iC.esConsumes<EcalGainRatios, EcalGainRatiosRcd>();
0061 }
0062
0063
0064 DataMixingEMDigiWorker::~DataMixingEMDigiWorker() {}
0065
0066 void DataMixingEMDigiWorker::addEMSignals(const edm::Event &e, const edm::EventSetup &ES) {
0067
0068
0069 LogInfo("DataMixingEMDigiWorker") << "===============> adding MC signals for " << e.id();
0070
0071
0072
0073 Handle<EBDigiCollection> pEBDigis;
0074
0075 const EBDigiCollection *EBDigis = nullptr;
0076
0077 if (e.getByToken(EBDigiToken_, pEBDigis)) {
0078 EBDigis = pEBDigis.product();
0079 LogDebug("DataMixingEMDigiWorker") << "total # EB digis: " << EBDigis->size();
0080 }
0081
0082 if (EBDigis) {
0083
0084
0085 for (EBDigiCollection::const_iterator it = EBDigis->begin(); it != EBDigis->end(); ++it) {
0086 EBDigiStorage_.insert(EBDigiMap::value_type((it->id()), *it));
0087 #ifdef DEBUG
0088
0089
0090
0091
0092 #endif
0093 }
0094 }
0095
0096
0097
0098 Handle<EEDigiCollection> pEEDigis;
0099
0100 const EEDigiCollection *EEDigis = nullptr;
0101
0102 if (e.getByToken(EEDigiToken_, pEEDigis)) {
0103 EEDigis = pEEDigis.product();
0104 LogDebug("DataMixingEMDigiWorker") << "total # EE digis: " << EEDigis->size();
0105 }
0106
0107 if (EEDigis) {
0108
0109 for (EEDigiCollection::const_iterator it = EEDigis->begin(); it != EEDigis->end(); ++it) {
0110 EEDigiStorage_.insert(EEDigiMap::value_type((it->id()), *it));
0111 #ifdef DEBUG
0112
0113
0114
0115
0116 #endif
0117 }
0118 }
0119
0120
0121 Handle<ESDigiCollection> pESDigis;
0122
0123 const ESDigiCollection *ESDigis = nullptr;
0124
0125 if (e.getByToken(ESDigiToken_, pESDigis)) {
0126 ESDigis = pESDigis.product();
0127 #ifdef DEBUG
0128 LogDebug("DataMixingEMDigiWorker") << "total # ES digis: " << ESDigis->size();
0129 #endif
0130 }
0131
0132 if (ESDigis) {
0133
0134 for (ESDigiCollection::const_iterator it = ESDigis->begin(); it != ESDigis->end(); ++it) {
0135 ESDigiStorage_.insert(ESDigiMap::value_type((it->id()), *it));
0136
0137 #ifdef DEBUG
0138
0139
0140
0141
0142 #endif
0143 }
0144 }
0145
0146 }
0147
0148 void DataMixingEMDigiWorker::addEMPileups(const int bcr,
0149 const EventPrincipal *ep,
0150 unsigned int eventNr,
0151 const edm::EventSetup &ES,
0152 ModuleCallingContext const *mcc) {
0153 LogInfo("DataMixingEMDigiWorker") << "\n===============> adding pileups from event " << ep->id()
0154 << " for bunchcrossing " << bcr;
0155
0156
0157
0158
0159
0160
0161 std::shared_ptr<Wrapper<EBDigiCollection> const> EBDigisPTR =
0162 getProductByTag<EBDigiCollection>(*ep, EBPileInputTag_, mcc);
0163
0164 if (EBDigisPTR) {
0165 const EBDigiCollection *EBDigis = const_cast<EBDigiCollection *>(EBDigisPTR->product());
0166
0167 LogDebug("DataMixingEMDigiWorker") << "total # EB digis: " << EBDigis->size();
0168
0169
0170 for (EBDigiCollection::const_iterator it = EBDigis->begin(); it != EBDigis->end(); ++it) {
0171 EBDigiStorage_.insert(EBDigiMap::value_type((it->id()), *it));
0172
0173 #ifdef DEBUG
0174
0175
0176
0177
0178 #endif
0179 }
0180 }
0181
0182
0183
0184 std::shared_ptr<Wrapper<EEDigiCollection> const> EEDigisPTR =
0185 getProductByTag<EEDigiCollection>(*ep, EEPileInputTag_, mcc);
0186
0187 if (EEDigisPTR) {
0188 const EEDigiCollection *EEDigis = const_cast<EEDigiCollection *>(EEDigisPTR->product());
0189
0190 LogDebug("DataMixingEMDigiWorker") << "total # EE digis: " << EEDigis->size();
0191
0192 for (EEDigiCollection::const_iterator it = EEDigis->begin(); it != EEDigis->end(); ++it) {
0193 EEDigiStorage_.insert(EEDigiMap::value_type((it->id()), *it));
0194
0195 #ifdef DEBUG
0196
0197
0198
0199
0200 #endif
0201 }
0202 }
0203
0204
0205 std::shared_ptr<Wrapper<ESDigiCollection> const> ESDigisPTR =
0206 getProductByTag<ESDigiCollection>(*ep, ESPileInputTag_, mcc);
0207
0208 if (ESDigisPTR) {
0209 const ESDigiCollection *ESDigis = const_cast<ESDigiCollection *>(ESDigisPTR->product());
0210
0211 LogDebug("DataMixingEMDigiWorker") << "total # ES digis: " << ESDigis->size();
0212
0213 for (ESDigiCollection::const_iterator it = ESDigis->begin(); it != ESDigis->end(); ++it) {
0214 ESDigiStorage_.insert(ESDigiMap::value_type((it->id()), *it));
0215
0216 #ifdef DEBUG
0217
0218
0219
0220
0221 #endif
0222 }
0223 }
0224 }
0225
0226 void DataMixingEMDigiWorker::putEM(edm::Event &e, const edm::EventSetup &ES) {
0227
0228 std::unique_ptr<EBDigiCollection> EBdigis(new EBDigiCollection);
0229 std::unique_ptr<EEDigiCollection> EEdigis(new EEDigiCollection);
0230 std::unique_ptr<ESDigiCollection> ESdigis(new ESDigiCollection);
0231
0232
0233
0234 DetId formerID = 0;
0235 DetId currentID;
0236
0237 EBDataFrame EB_old;
0238
0239 int gain_new = 0;
0240 int gain_old = 0;
0241 int gain_consensus = 0;
0242 int adc_new;
0243 int adc_old;
0244 int adc_sum;
0245 uint16_t data;
0246
0247
0248
0249 EBDigiMap::const_iterator iEBchk;
0250
0251 for (EBDigiMap::const_iterator iEB = EBDigiStorage_.begin(); iEB != EBDigiStorage_.end(); iEB++) {
0252 currentID = iEB->first;
0253
0254 if (currentID == formerID) {
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268 unsigned int sizenew = (iEB->second).size();
0269 unsigned int sizeold = EB_old.size();
0270
0271 unsigned int max_samp = std::max(sizenew, sizeold);
0272
0273
0274
0275
0276
0277 int sw_gain_consensus = 0;
0278
0279 for (unsigned int isamp = 0; isamp < max_samp; isamp++) {
0280 if (isamp < sizenew) {
0281 gain_new = (iEB->second)[isamp].gainId();
0282 adc_new = (iEB->second)[isamp].adc();
0283 } else {
0284 adc_new = 0;
0285 }
0286
0287 if (isamp < sizeold) {
0288 gain_old = EB_old[isamp].gainId();
0289 adc_old = EB_old[isamp].adc();
0290 } else {
0291 adc_old = 0;
0292 }
0293
0294 const std::vector<float> pedeStals = GetPedestals(ES, currentID);
0295 const std::vector<float> gainRatios = GetGainRatios(ES, currentID);
0296
0297 if (adc_new > 0 && adc_old > 0) {
0298 if (gain_old == gain_new) {
0299 gain_consensus = gain_old;
0300 } else {
0301
0302 if (gain_old < gain_new) {
0303
0304 float ratio = gainRatios[gain_new - 1] / gainRatios[gain_old - 1];
0305 adc_old = (int)round((adc_old - pedeStals[gain_old - 1]) / ratio + pedeStals[gain_new - 1]);
0306 gain_consensus = gain_new;
0307 } else {
0308 float ratio = gainRatios[gain_old - 1] / gainRatios[gain_new - 1];
0309 adc_new = (int)round((adc_new - pedeStals[gain_new - 1]) / ratio + pedeStals[gain_old - 1]);
0310 gain_consensus = gain_old;
0311 }
0312 }
0313 }
0314
0315
0316 adc_sum = adc_new + adc_old - (int)round(pedeStals[gain_consensus - 1]);
0317
0318
0319 if (adc_sum > 4096) {
0320 if (gain_consensus < 3) {
0321 double ratio = gainRatios[gain_consensus] / gainRatios[gain_consensus - 1];
0322 adc_sum = (int)round((adc_sum - pedeStals[gain_consensus - 1]) / ratio + pedeStals[gain_consensus]);
0323 sw_gain_consensus = ++gain_consensus;
0324 } else
0325 adc_sum = 4096;
0326 }
0327
0328
0329
0330 if (gain_consensus < sw_gain_consensus) {
0331 double ratio = gainRatios[sw_gain_consensus - 1] / gainRatios[gain_consensus - 1];
0332 adc_sum = (int)round((adc_sum - pedeStals[gain_consensus - 1]) / ratio + pedeStals[sw_gain_consensus - 1]);
0333 gain_consensus = sw_gain_consensus;
0334 }
0335
0336 EcalMGPASample sample(adc_sum, gain_consensus);
0337 EB_old.setSample(isamp,
0338 sample);
0339 }
0340
0341 }
0342 else {
0343 if (formerID > 0) {
0344 EBdigis->push_back(formerID, EB_old.frame().begin());
0345 }
0346
0347 formerID = currentID;
0348 EB_old = iEB->second;
0349 }
0350
0351 iEBchk = iEB;
0352 if ((++iEBchk) == EBDigiStorage_.end()) {
0353 EBdigis->push_back(currentID, (iEB->second).frame().begin());
0354 }
0355 }
0356
0357
0358
0359 formerID = 0;
0360 EEDataFrame EE_old;
0361
0362 EEDigiMap::const_iterator iEEchk;
0363
0364 for (EEDigiMap::const_iterator iEE = EEDigiStorage_.begin(); iEE != EEDigiStorage_.end(); iEE++) {
0365 currentID = iEE->first;
0366
0367 if (currentID == formerID) {
0368
0369
0370 unsigned int sizenew = (iEE->second).size();
0371 unsigned int sizeold = EE_old.size();
0372
0373 unsigned int max_samp = std::max(sizenew, sizeold);
0374
0375
0376
0377
0378
0379 for (unsigned int isamp = 0; isamp < max_samp; isamp++) {
0380 if (isamp < sizenew) {
0381 gain_new = (iEE->second)[isamp].gainId();
0382 adc_new = (iEE->second)[isamp].adc();
0383 } else {
0384 adc_new = 0;
0385 }
0386
0387 if (isamp < sizeold) {
0388 gain_old = EE_old[isamp].gainId();
0389 adc_old = EE_old[isamp].adc();
0390 } else {
0391 adc_old = 0;
0392 }
0393
0394 const std::vector<float> pedeStals = GetPedestals(ES, currentID);
0395 const std::vector<float> gainRatios = GetGainRatios(ES, currentID);
0396
0397 if (adc_new > 0 && adc_old > 0) {
0398 if (gain_old == gain_new) {
0399 gain_consensus = gain_old;
0400 } else {
0401
0402 if (gain_old < gain_new) {
0403
0404 float ratio = gainRatios[gain_new - 1] / gainRatios[gain_old - 1];
0405 adc_old = (int)round((adc_old - pedeStals[gain_old - 1]) / ratio + pedeStals[gain_new - 1]);
0406 gain_consensus = gain_new;
0407 } else {
0408 float ratio = gainRatios[gain_old - 1] / gainRatios[gain_new - 1];
0409 adc_new = (int)round((adc_new - pedeStals[gain_new - 1]) / ratio + pedeStals[gain_old - 1]);
0410 gain_consensus = gain_old;
0411 }
0412 }
0413 }
0414
0415
0416 adc_sum = adc_new + adc_old - (int)round(pedeStals[gain_consensus - 1]);
0417
0418
0419 if (adc_sum > 4096) {
0420 if (gain_consensus < 3) {
0421 double ratio = gainRatios[gain_consensus] / gainRatios[gain_consensus - 1];
0422 adc_sum = (int)round((adc_sum - pedeStals[gain_consensus - 1]) / ratio + pedeStals[gain_consensus]);
0423 ++gain_consensus;
0424 } else
0425 adc_sum = 4096;
0426 }
0427
0428 EcalMGPASample sample(adc_sum, gain_consensus);
0429 EE_old.setSample(isamp, sample);
0430 }
0431
0432 } else {
0433 if (formerID > 0) {
0434 EEdigis->push_back(formerID, EE_old.frame().begin());
0435 }
0436
0437 formerID = currentID;
0438 EE_old = iEE->second;
0439 }
0440
0441 iEEchk = iEE;
0442 if ((++iEEchk) == EEDigiStorage_.end()) {
0443 EEdigis->push_back(currentID, (iEE->second).frame().begin());
0444 }
0445 }
0446
0447
0448
0449 formerID = 0;
0450 ESDataFrame ES_old;
0451
0452 ESDigiMap::const_iterator iESchk;
0453
0454 for (ESDigiMap::const_iterator iES = ESDigiStorage_.begin(); iES != ESDigiStorage_.end(); iES++) {
0455 currentID = iES->first;
0456
0457 if (currentID == formerID) {
0458
0459
0460 unsigned int sizenew = (iES->second).size();
0461 unsigned int sizeold = ES_old.size();
0462 uint16_t rawdat = 0;
0463 unsigned int max_samp = std::max(sizenew, sizeold);
0464
0465
0466
0467
0468
0469 for (unsigned int isamp = 0; isamp < max_samp; isamp++) {
0470 if (isamp < sizenew) {
0471 adc_new = (iES->second)[isamp].adc();
0472 rawdat = (iES->second)[isamp].raw();
0473 } else {
0474 adc_new = 0;
0475 }
0476
0477 if (isamp < sizeold) {
0478 adc_old = ES_old[isamp].adc();
0479 rawdat = ES_old[isamp].raw();
0480 } else {
0481 adc_old = 0;
0482 }
0483
0484
0485 adc_sum = adc_new + adc_old;
0486
0487 adc_sum = std::min(adc_sum, 4095);
0488 data = adc_sum + (rawdat & 0xF000);
0489 ES_old.setSample(isamp, data);
0490 }
0491
0492 } else {
0493 if (formerID > 0) {
0494 ESdigis->push_back(ES_old);
0495 }
0496
0497 formerID = currentID;
0498 ES_old = iES->second;
0499 }
0500
0501 iESchk = iES;
0502 if ((++iESchk) == ESDigiStorage_.end()) {
0503 ESdigis->push_back(iES->second);
0504
0505
0506
0507
0508 }
0509 }
0510
0511
0512
0513
0514 LogInfo("DataMixingEMDigiWorker") << "total # EB Merged digis: " << EBdigis->size();
0515 LogInfo("DataMixingEMDigiWorker") << "total # EE Merged digis: " << EEdigis->size();
0516 LogInfo("DataMixingEMDigiWorker") << "total # ES Merged digis: " << ESdigis->size();
0517
0518 e.put(std::move(EBdigis), EBDigiCollectionDM_);
0519 e.put(std::move(EEdigis), EEDigiCollectionDM_);
0520 e.put(std::move(ESdigis), ESDigiCollectionDM_);
0521
0522
0523
0524 EBDigiStorage_.clear();
0525 EEDigiStorage_.clear();
0526 ESDigiStorage_.clear();
0527 }
0528 const std::vector<float> DataMixingEMDigiWorker::GetPedestals(const edm::EventSetup &ES, const DetId &detid) {
0529 std::vector<float> pedeStals(3);
0530
0531
0532 const auto &pedHandle = ES.getHandle(pedToken_);
0533
0534 const EcalPedestalsMap &pedMap = pedHandle.product()->getMap();
0535 EcalPedestalsMapIterator pedIter;
0536 EcalPedestals::Item aped;
0537
0538 pedIter = pedMap.find(detid);
0539 if (pedIter != pedMap.end()) {
0540 aped = (*pedIter);
0541 pedeStals[0] = aped.mean_x12;
0542 pedeStals[1] = aped.mean_x6;
0543 pedeStals[2] = aped.mean_x1;
0544 } else {
0545 edm::LogError("DataMixingMissingInput") << "Cannot find pedestals";
0546 pedeStals[0] = 0;
0547 pedeStals[1] = 0;
0548 pedeStals[2] = 0;
0549 }
0550
0551 return pedeStals;
0552 }
0553
0554 const std::vector<float> DataMixingEMDigiWorker::GetGainRatios(const edm::EventSetup &ES, const DetId &detid) {
0555 std::vector<float> gainRatios(3);
0556
0557 const auto &grHandle = ES.getHandle(grToken_);
0558 EcalMGPAGainRatio theRatio = (*grHandle)[detid];
0559
0560 gainRatios[0] = 1.;
0561 gainRatios[1] = theRatio.gain12Over6();
0562 gainRatios[2] = theRatio.gain6Over1() * theRatio.gain12Over6();
0563
0564 return gainRatios;
0565 }
0566
0567 }