File indexing completed on 2021-02-14 14:24:22
0001 #include "RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitWorkerRecover.h"
0002
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
0006 #include "CondFormats/DataRecord/interface/EcalTimeCalibConstantsRcd.h"
0007 #include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
0008
0009 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0010 #include "DataFormats/EcalDetId/interface/EcalScDetId.h"
0011
0012 #include "CondFormats/EcalObjects/interface/EcalTimeCalibConstants.h"
0013
0014 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0015 #include "FWCore/Framework/interface/EDProducer.h"
0016
0017 EcalRecHitWorkerRecover::EcalRecHitWorkerRecover(const edm::ParameterSet& ps, edm::ConsumesCollector& c)
0018 : EcalRecHitWorkerBaseClass(ps, c), ecalScaleTokens_(c), tpgscaleTokens_(c) {
0019 rechitMaker_ = std::make_unique<EcalRecHitSimpleAlgo>();
0020
0021 singleRecoveryMethod_ = ps.getParameter<std::string>("singleChannelRecoveryMethod");
0022 singleRecoveryThreshold_ = ps.getParameter<double>("singleChannelRecoveryThreshold");
0023 sum8RecoveryThreshold_ = ps.getParameter<double>("sum8ChannelRecoveryThreshold");
0024 killDeadChannels_ = ps.getParameter<bool>("killDeadChannels");
0025 recoverEBIsolatedChannels_ = ps.getParameter<bool>("recoverEBIsolatedChannels");
0026 recoverEEIsolatedChannels_ = ps.getParameter<bool>("recoverEEIsolatedChannels");
0027 recoverEBVFE_ = ps.getParameter<bool>("recoverEBVFE");
0028 recoverEEVFE_ = ps.getParameter<bool>("recoverEEVFE");
0029 recoverEBFE_ = ps.getParameter<bool>("recoverEBFE");
0030 recoverEEFE_ = ps.getParameter<bool>("recoverEEFE");
0031 laserToken_ = c.esConsumes<EcalLaserDbService, EcalLaserDbRecord>();
0032 caloTopologyToken_ = c.esConsumes<CaloTopology, CaloTopologyRecord>();
0033 pEcalMappingToken_ = c.esConsumes<EcalElectronicsMapping, EcalMappingRcd>();
0034 pEBGeomToken_ = c.esConsumes<CaloSubdetectorGeometry, EcalBarrelGeometryRecord>(edm::ESInputTag("", "EcalBarrel"));
0035 caloGeometryToken_ = c.esConsumes<CaloGeometry, CaloGeometryRecord>();
0036 chStatusToken_ = c.esConsumes<EcalChannelStatus, EcalChannelStatusRcd>();
0037 ttMapToken_ = c.esConsumes<EcalTrigTowerConstituentsMap, IdealGeometryRecord>();
0038
0039 dbStatusToBeExcludedEE_ = ps.getParameter<std::vector<int> >("dbStatusToBeExcludedEE");
0040 dbStatusToBeExcludedEB_ = ps.getParameter<std::vector<int> >("dbStatusToBeExcludedEB");
0041
0042 logWarningEtThreshold_EB_FE_ = ps.getParameter<double>("logWarningEtThreshold_EB_FE");
0043 logWarningEtThreshold_EE_FE_ = ps.getParameter<double>("logWarningEtThreshold_EE_FE");
0044
0045 tpDigiToken_ =
0046 c.consumes<EcalTrigPrimDigiCollection>(ps.getParameter<edm::InputTag>("triggerPrimitiveDigiCollection"));
0047
0048 if (recoverEBIsolatedChannels_ && singleRecoveryMethod_ == "BDTG")
0049 ebDeadChannelCorrector.setParameters(ps);
0050 }
0051
0052 void EcalRecHitWorkerRecover::set(const edm::EventSetup& es) {
0053 laser = es.getHandle(laserToken_);
0054 caloTopology_ = es.getHandle(caloTopologyToken_);
0055 pEcalMapping_ = es.getHandle(pEcalMappingToken_);
0056 ecalMapping_ = pEcalMapping_.product();
0057
0058 pEBGeom_ = es.getHandle(pEBGeomToken_);
0059 caloGeometry_ = es.getHandle(caloGeometryToken_);
0060 chStatus_ = es.getHandle(chStatusToken_);
0061 geo_ = caloGeometry_.product();
0062 ebGeom_ = pEBGeom_.product();
0063 ttMap_ = es.getHandle(ttMapToken_);
0064 recoveredDetIds_EB_.clear();
0065 recoveredDetIds_EE_.clear();
0066 eventSetup_ = &es;
0067 }
0068
0069 bool EcalRecHitWorkerRecover::run(const edm::Event& evt,
0070 const EcalUncalibratedRecHit& uncalibRH,
0071 EcalRecHitCollection& result) {
0072 DetId detId = uncalibRH.id();
0073 uint32_t flags = (0xF & uncalibRH.flags());
0074
0075
0076
0077
0078
0079
0080
0081
0082 if (killDeadChannels_) {
0083 if ((flags == EcalRecHitWorkerRecover::EB_single && !recoverEBIsolatedChannels_) ||
0084 (flags == EcalRecHitWorkerRecover::EE_single && !recoverEEIsolatedChannels_) ||
0085 (flags == EcalRecHitWorkerRecover::EB_VFE && !recoverEBVFE_) ||
0086 (flags == EcalRecHitWorkerRecover::EE_VFE && !recoverEEVFE_)) {
0087 EcalRecHit hit(detId, 0., 0., EcalRecHit::kDead);
0088 hit.setFlag(EcalRecHit::kDead);
0089 insertRecHit(hit, result);
0090 return true;
0091 }
0092 if (flags == EcalRecHitWorkerRecover::EB_FE && !recoverEBFE_) {
0093 EcalTrigTowerDetId ttDetId(((EBDetId)detId).tower());
0094 std::vector<DetId> vid = ttMap_->constituentsOf(ttDetId);
0095 for (std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); ++dit) {
0096 EcalRecHit hit((*dit), 0., 0., EcalRecHit::kDead);
0097 hit.setFlag(EcalRecHit::kDead);
0098 insertRecHit(hit, result);
0099 }
0100 if (logWarningEtThreshold_EB_FE_ < 0)
0101 return true;
0102 }
0103 if (flags == EcalRecHitWorkerRecover::EE_FE && !recoverEEFE_) {
0104 EEDetId id(detId);
0105 EcalScDetId sc(1 + (id.ix() - 1) / 5, 1 + (id.iy() - 1) / 5, id.zside());
0106 std::vector<DetId> eeC;
0107 for (int dx = 1; dx <= 5; ++dx) {
0108 for (int dy = 1; dy <= 5; ++dy) {
0109 int ix = (sc.ix() - 1) * 5 + dx;
0110 int iy = (sc.iy() - 1) * 5 + dy;
0111 int iz = sc.zside();
0112 if (EEDetId::validDetId(ix, iy, iz)) {
0113 eeC.push_back(EEDetId(ix, iy, iz));
0114 }
0115 }
0116 }
0117 for (size_t i = 0; i < eeC.size(); ++i) {
0118 EcalRecHit hit(eeC[i], 0., 0., EcalRecHit::kDead);
0119 hit.setFlag(EcalRecHit::kDead);
0120 insertRecHit(hit, result);
0121 }
0122 if (logWarningEtThreshold_EE_FE_ < 0)
0123 return true;
0124 }
0125 }
0126
0127 if (flags == EcalRecHitWorkerRecover::EB_single) {
0128
0129 ebDeadChannelCorrector.setCaloTopology(caloTopology_.product());
0130
0131
0132 bool AcceptRecHit = true;
0133 float ebEn = ebDeadChannelCorrector.correct(
0134 detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, sum8RecoveryThreshold_, &AcceptRecHit);
0135 EcalRecHit hit(detId, ebEn, 0., EcalRecHit::kDead);
0136
0137 if (hit.energy() != 0 and AcceptRecHit == true) {
0138 hit.setFlag(EcalRecHit::kNeighboursRecovered);
0139 } else {
0140
0141 hit.setFlag(EcalRecHit::kDead);
0142 }
0143 insertRecHit(hit, result);
0144
0145 } else if (flags == EcalRecHitWorkerRecover::EE_single) {
0146
0147 eeDeadChannelCorrector.setCaloTopology(caloTopology_.product());
0148
0149
0150 bool AcceptRecHit = true;
0151 float eeEn = eeDeadChannelCorrector.correct(
0152 detId, result, singleRecoveryMethod_, singleRecoveryThreshold_, sum8RecoveryThreshold_, &AcceptRecHit);
0153 EcalRecHit hit(detId, eeEn, 0., EcalRecHit::kDead);
0154 if (hit.energy() != 0 and AcceptRecHit == true) {
0155 hit.setFlag(EcalRecHit::kNeighboursRecovered);
0156 } else {
0157
0158 hit.setFlag(EcalRecHit::kDead);
0159 }
0160 insertRecHit(hit, result);
0161
0162 } else if (flags == EcalRecHitWorkerRecover::EB_VFE) {
0163
0164 EcalRecHit hit(detId, 0., 0.);
0165 hit.setFlag(EcalRecHit::kDead);
0166
0167 insertRecHit(hit, result);
0168 } else if (flags == EcalRecHitWorkerRecover::EB_FE) {
0169
0170
0171 EcalTrigTowerDetId ttDetId(((EBDetId)detId).tower());
0172 edm::Handle<EcalTrigPrimDigiCollection> pTPDigis;
0173 evt.getByToken(tpDigiToken_, pTPDigis);
0174 const EcalTrigPrimDigiCollection* tpDigis = nullptr;
0175 tpDigis = pTPDigis.product();
0176
0177 EcalTrigPrimDigiCollection::const_iterator tp = tpDigis->find(ttDetId);
0178
0179 if (tp != tpDigis->end()) {
0180 EcalTPGScale ecalScale{ecalScaleTokens_, *eventSetup_};
0181
0182 std::vector<DetId> vid = ttMap_->constituentsOf(ttDetId);
0183 float tpEt = ecalScale.getTPGInGeV(tp->compressedEt(), tp->id());
0184 float tpEtThreshEB = logWarningEtThreshold_EB_FE_;
0185 if (tpEt > tpEtThreshEB) {
0186 edm::LogWarning("EnergyInDeadEB_FE") << "TP energy in the dead TT = " << tpEt << " at " << ttDetId;
0187 }
0188 if (!killDeadChannels_ || recoverEBFE_) {
0189
0190
0191 for (std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); ++dit) {
0192 if (alreadyInserted(*dit))
0193 continue;
0194 float theta = ebGeom_->getGeometry(*dit)->getPosition().theta();
0195 float tpEt = ecalScale.getTPGInGeV(tp->compressedEt(), tp->id());
0196 if (checkChannelStatus(*dit, dbStatusToBeExcludedEB_)) {
0197 EcalRecHit hit(*dit, tpEt / ((float)vid.size()) / sin(theta), 0.);
0198 hit.setFlag(EcalRecHit::kTowerRecovered);
0199 if (tp->compressedEt() == 0xFF)
0200 hit.setFlag(EcalRecHit::kTPSaturated);
0201 if (tp->sFGVB())
0202 hit.setFlag(EcalRecHit::kL1SpikeFlag);
0203 insertRecHit(hit, result);
0204 }
0205 }
0206 } else {
0207
0208 std::vector<DetId> vid = ttMap_->constituentsOf(ttDetId);
0209 for (std::vector<DetId>::const_iterator dit = vid.begin(); dit != vid.end(); ++dit) {
0210 if (alreadyInserted(*dit))
0211 continue;
0212 EcalRecHit hit(*dit, 0., 0.);
0213 hit.setFlag(EcalRecHit::kDead);
0214 insertRecHit(hit, result);
0215 }
0216 }
0217 }
0218 } else if (flags == EcalRecHitWorkerRecover::EE_FE) {
0219
0220
0221
0222
0223
0224
0225
0226
0227 EEDetId eeId(detId);
0228 EcalScDetId sc((eeId.ix() - 1) / 5 + 1, (eeId.iy() - 1) / 5 + 1, eeId.zside());
0229 std::set<DetId> eeC;
0230 for (int dx = 1; dx <= 5; ++dx) {
0231 for (int dy = 1; dy <= 5; ++dy) {
0232 int ix = (sc.ix() - 1) * 5 + dx;
0233 int iy = (sc.iy() - 1) * 5 + dy;
0234 int iz = sc.zside();
0235 if (EEDetId::validDetId(ix, iy, iz)) {
0236 EEDetId id(ix, iy, iz);
0237 if (checkChannelStatus(id, dbStatusToBeExcludedEE_)) {
0238 eeC.insert(id);
0239 }
0240 }
0241 }
0242 }
0243
0244 edm::Handle<EcalTrigPrimDigiCollection> pTPDigis;
0245 evt.getByToken(tpDigiToken_, pTPDigis);
0246 const EcalTrigPrimDigiCollection* tpDigis = nullptr;
0247 tpDigis = pTPDigis.product();
0248
0249
0250 std::set<EcalTrigTowerDetId> aTT;
0251 for (std::set<DetId>::const_iterator it = eeC.begin(); it != eeC.end(); ++it) {
0252 aTT.insert(ttMap_->towerOf(*it));
0253 }
0254
0255 EcalTPGScale tpgscale(tpgscaleTokens_, *eventSetup_);
0256 EcalTPGScale ecalScale(ecalScaleTokens_, *eventSetup_);
0257
0258 float totE = 0;
0259
0260 std::set<DetId> aTTC;
0261 bool atLeastOneTPSaturated = false;
0262 for (std::set<EcalTrigTowerDetId>::const_iterator it = aTT.begin(); it != aTT.end(); ++it) {
0263
0264 EcalTrigPrimDigiCollection::const_iterator itTP = tpDigis->find(*it);
0265 if (itTP != tpDigis->end()) {
0266 std::vector<DetId> v = ttMap_->constituentsOf(*it);
0267
0268
0269 std::vector<DetId>::iterator ttcons = v.begin();
0270 while (ttcons != v.end()) {
0271 if (!checkChannelStatus(*ttcons, dbStatusToBeExcludedEE_)) {
0272 ttcons = v.erase(ttcons);
0273 } else {
0274 ++ttcons;
0275 }
0276 }
0277
0278 if (itTP->compressedEt() == 0xFF) {
0279 atLeastOneTPSaturated =
0280 true;
0281
0282
0283
0284
0285
0286 totE += estimateEnergy(itTP->id().ietaAbs(), &result, eeC, v, tpgscale);
0287
0288
0289
0290
0291
0292
0293
0294
0295
0296
0297
0298
0299
0300 } else {
0301 totE += ((it->ietaAbs() > 26) ? 2 : 1) * ecalScale.getTPGInGeV(itTP->compressedEt(), itTP->id());
0302 }
0303
0304
0305
0306 if (itTP->compressedEt() == 0) {
0307 for (size_t i = 0; i < v.size(); ++i)
0308 eeC.erase(v[i]);
0309 } else if (itTP->compressedEt() != 0xFF) {
0310
0311 for (size_t j = 0; j < v.size(); ++j) {
0312 aTTC.insert(v[j]);
0313 }
0314 }
0315 }
0316 }
0317
0318
0319
0320
0321 for (std::set<DetId>::const_iterator it = eeC.begin(); it != eeC.end(); ++it) {
0322 aTTC.erase(*it);
0323 }
0324
0325 const EcalRecHitCollection* hits = &result;
0326 for (std::set<DetId>::const_iterator it = aTTC.begin(); it != aTTC.end(); ++it) {
0327 EcalRecHitCollection::const_iterator jt = hits->find(*it);
0328 if (jt != hits->end()) {
0329 float energy = jt->energy();
0330 float eta = geo_->getPosition(jt->id()).eta();
0331 float pf = 1.0 / cosh(eta);
0332
0333 totE -= energy * pf;
0334 }
0335 }
0336
0337 float scEt = totE;
0338 float scEtThreshEE = logWarningEtThreshold_EE_FE_;
0339 if (scEt > scEtThreshEE) {
0340 edm::LogWarning("EnergyInDeadEE_FE") << "TP energy in the dead TT = " << scEt << " at " << sc;
0341 }
0342
0343
0344 if (!killDeadChannels_ || recoverEEFE_) {
0345
0346 for (std::set<DetId>::const_iterator it = eeC.begin(); it != eeC.end(); ++it) {
0347 float eta = geo_->getPosition(*it).eta();
0348 float pf = 1.0 / cosh(eta);
0349 EcalRecHit hit(*it, totE / ((float)eeC.size() * pf), 0);
0350
0351 if (atLeastOneTPSaturated)
0352 hit.setFlag(EcalRecHit::kTPSaturated);
0353 hit.setFlag(EcalRecHit::kTowerRecovered);
0354 insertRecHit(hit, result);
0355
0356 }
0357 }
0358 }
0359 return true;
0360 }
0361
0362 float EcalRecHitWorkerRecover::estimateEnergy(int ieta,
0363 EcalRecHitCollection* hits,
0364 const std::set<DetId>& sId,
0365 const std::vector<DetId>& vId,
0366 const EcalTPGScale& tpgscale) {
0367 float xtalE = 0;
0368 int count = 0;
0369 for (std::vector<DetId>::const_iterator vIdit = vId.begin(); vIdit != vId.end(); ++vIdit) {
0370 std::set<DetId>::const_iterator sIdit = sId.find(*vIdit);
0371 if (sIdit == sId.end()) {
0372 float energy = hits->find(*vIdit)->energy();
0373 float eta = geo_->getPosition(*vIdit).eta();
0374 float pf = 1.0 / cosh(eta);
0375 xtalE += energy * pf;
0376 count++;
0377 }
0378 }
0379
0380 if (count == 0) {
0381
0382 double etsat = tpgscale.getTPGInGeV(0xFF,
0383 ttMap_->towerOf(*vId.begin()));
0384
0385
0386 return etsat / cosh(ieta) * (ieta > 26 ? 2 : 1);
0387 } else
0388 return xtalE * ((vId.size() / (float)count) - 1) * (ieta > 26 ? 2 : 1);
0389 }
0390
0391 void EcalRecHitWorkerRecover::insertRecHit(const EcalRecHit& hit, EcalRecHitCollection& collection) {
0392
0393 if (alreadyInserted(hit.id())) {
0394 edm::LogWarning("EcalRecHitWorkerRecover") << "DetId already recovered! Skipping...";
0395 return;
0396 }
0397 EcalRecHitCollection::iterator it = collection.find(hit.id());
0398 if (it == collection.end()) {
0399
0400 collection.push_back(hit);
0401 } else {
0402
0403 *it = hit;
0404 }
0405 if (hit.id().subdetId() == EcalBarrel) {
0406 recoveredDetIds_EB_.insert(hit.id());
0407 } else if (hit.id().subdetId() == EcalEndcap) {
0408 recoveredDetIds_EE_.insert(hit.id());
0409 } else {
0410 edm::LogError("EcalRecHitWorkerRecover::InvalidDetId") << "Invalid DetId " << hit.id().rawId();
0411 }
0412 }
0413
0414 bool EcalRecHitWorkerRecover::alreadyInserted(const DetId& id) {
0415 bool res = false;
0416 if (id.subdetId() == EcalBarrel) {
0417 res = (recoveredDetIds_EB_.find(id) != recoveredDetIds_EB_.end());
0418 } else if (id.subdetId() == EcalEndcap) {
0419 res = (recoveredDetIds_EE_.find(id) != recoveredDetIds_EE_.end());
0420 } else {
0421 edm::LogError("EcalRecHitWorkerRecover::InvalidDetId") << "Invalid DetId " << id.rawId();
0422 }
0423 return res;
0424 }
0425
0426
0427
0428 float EcalRecHitWorkerRecover::recCheckCalib(float eTT, int ieta) { return eTT; }
0429
0430
0431
0432
0433 bool EcalRecHitWorkerRecover::checkChannelStatus(const DetId& id, const std::vector<int>& statusestoexclude) {
0434 if (!chStatus_.isValid())
0435 edm::LogError("ObjectNotFound") << "Channel Status not set";
0436
0437 EcalChannelStatus::const_iterator chIt = chStatus_->find(id);
0438 uint16_t dbStatus = 0;
0439 if (chIt != chStatus_->end()) {
0440 dbStatus = chIt->getEncodedStatusCode();
0441 } else {
0442 edm::LogError("ObjectNotFound") << "No channel status found for xtal " << id.rawId()
0443 << "! something wrong with EcalChannelStatus in your DB? ";
0444 }
0445
0446 for (std::vector<int>::const_iterator status = statusestoexclude.begin(); status != statusestoexclude.end();
0447 ++status) {
0448 if (*status == dbStatus)
0449 return false;
0450 }
0451
0452 return true;
0453 }
0454
0455 #include "FWCore/Framework/interface/MakerMacros.h"
0456 #include "RecoLocalCalo/EcalRecProducers/interface/EcalRecHitWorkerFactory.h"
0457 DEFINE_EDM_PLUGIN(EcalRecHitWorkerFactory, EcalRecHitWorkerRecover, "EcalRecHitWorkerRecover");