EcalRecProducers/plugins/EcalRecHitProducer.cc

0001 /** \class EcalRecHitProducer
0002  *   produce ECAL rechits from uncalibrated rechits
0003  *
0004  *  \author Shahram Rahatlou, University of Rome & INFN, March 2006
0005  *
0006  **/
0007
0008 #include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
0009 #include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"
0010 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0011 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0012 #include "DataFormats/EcalDetId/interface/EcalElectronicsId.h"
0013 #include "DataFormats/EcalDetId/interface/EcalScDetId.h"
0014 #include "DataFormats/EcalDetId/interface/EcalTrigTowerDetId.h"
0015 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0016 #include "FWCore/Framework/interface/Event.h"
0017 #include "FWCore/Framework/interface/EventSetup.h"
0018 #include "FWCore/Framework/interface/stream/EDProducer.h"
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0021 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0023 #include "FWCore/Utilities/interface/ESGetToken.h"
0024 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalCleaningAlgo.h"
0025 #include "RecoLocalCalo/EcalRecProducers/interface/EcalRecHitWorkerBaseClass.h"
0026 #include "RecoLocalCalo/EcalRecProducers/interface/EcalRecHitWorkerFactory.h"
0027
0028 class EcalRecHitProducer : public edm::stream::EDProducer<> {
0029 public:
0030   explicit EcalRecHitProducer(const edm::ParameterSet& ps);
0031   ~EcalRecHitProducer() override;
0032   void produce(edm::Event& evt, const edm::EventSetup& es) override;
0033   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0034
0035 private:
0036   const bool doEB_;  // consume and use the EB uncalibrated RecHits. An EB collection is produced even if this is false
0037   const bool doEE_;  // consume and use the EE uncalibrated RecHits. An EE collection is produced even if this is false
0038   const bool recoverEBIsolatedChannels_;
0039   const bool recoverEEIsolatedChannels_;
0040   const bool recoverEBVFE_;
0041   const bool recoverEEVFE_;
0042   const bool recoverEBFE_;
0043   const bool recoverEEFE_;
0044   const bool killDeadChannels_;
0045
0046   std::unique_ptr<EcalRecHitWorkerBaseClass> worker_;
0047   std::unique_ptr<EcalRecHitWorkerBaseClass> workerRecover_;
0048
0049   std::unique_ptr<EcalCleaningAlgo> cleaningAlgo_;
0050
0051   edm::EDGetTokenT<EBUncalibratedRecHitCollection> ebUncalibRecHitToken_;
0052   edm::EDGetTokenT<EEUncalibratedRecHitCollection> eeUncalibRecHitToken_;
0053   edm::EDGetTokenT<std::set<EBDetId>> ebDetIdToBeRecoveredToken_;
0054   edm::EDGetTokenT<std::set<EEDetId>> eeDetIdToBeRecoveredToken_;
0055   edm::EDGetTokenT<std::set<EcalTrigTowerDetId>> ebFEToBeRecoveredToken_;
0056   edm::EDGetTokenT<std::set<EcalScDetId>> eeFEToBeRecoveredToken_;
0057   edm::ESGetToken<EcalChannelStatus, EcalChannelStatusRcd> ecalChannelStatusToken_;
0058   const edm::EDPutTokenT<EBRecHitCollection> ebRecHitToken_;
0059   const edm::EDPutTokenT<EERecHitCollection> eeRecHitToken_;
0060 };
0061
0062 EcalRecHitProducer::EcalRecHitProducer(const edm::ParameterSet& ps)
0063     : doEB_(!ps.getParameter<edm::InputTag>("EBuncalibRecHitCollection").label().empty()),
0064       doEE_(!ps.getParameter<edm::InputTag>("EEuncalibRecHitCollection").label().empty()),
0065       recoverEBIsolatedChannels_(ps.getParameter<bool>("recoverEBIsolatedChannels")),
0066       recoverEEIsolatedChannels_(ps.getParameter<bool>("recoverEEIsolatedChannels")),
0067       recoverEBVFE_(ps.getParameter<bool>("recoverEBVFE")),
0068       recoverEEVFE_(ps.getParameter<bool>("recoverEEVFE")),
0069       recoverEBFE_(ps.getParameter<bool>("recoverEBFE")),
0070       recoverEEFE_(ps.getParameter<bool>("recoverEEFE")),
0071       killDeadChannels_(ps.getParameter<bool>("killDeadChannels")),
0072       ebRecHitToken_(produces<EBRecHitCollection>(ps.getParameter<std::string>("EBrechitCollection"))),
0073       eeRecHitToken_(produces<EERecHitCollection>(ps.getParameter<std::string>("EErechitCollection"))) {
0074   if (doEB_) {
0075     ebUncalibRecHitToken_ =
0076         consumes<EBUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("EBuncalibRecHitCollection"));
0077
0078     if (recoverEBIsolatedChannels_ || recoverEBFE_ || killDeadChannels_) {
0079       ebDetIdToBeRecoveredToken_ = consumes<std::set<EBDetId>>(ps.getParameter<edm::InputTag>("ebDetIdToBeRecovered"));
0080     }
0081
0082     if (recoverEBFE_ || killDeadChannels_) {
0083       ebFEToBeRecoveredToken_ =
0084           consumes<std::set<EcalTrigTowerDetId>>(ps.getParameter<edm::InputTag>("ebFEToBeRecovered"));
0085     }
0086   }
0087
0088   if (doEE_) {
0089     eeUncalibRecHitToken_ =
0090         consumes<EEUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("EEuncalibRecHitCollection"));
0091
0092     if (recoverEEIsolatedChannels_ || recoverEEFE_ || killDeadChannels_) {
0093       eeDetIdToBeRecoveredToken_ = consumes<std::set<EEDetId>>(ps.getParameter<edm::InputTag>("eeDetIdToBeRecovered"));
0094     }
0095
0096     if (recoverEEFE_ || killDeadChannels_) {
0097       eeFEToBeRecoveredToken_ = consumes<std::set<EcalScDetId>>(ps.getParameter<edm::InputTag>("eeFEToBeRecovered"));
0098     }
0099   }
0100
0101   if (recoverEBIsolatedChannels_ || recoverEBFE_ || recoverEEIsolatedChannels_ || recoverEEFE_ || killDeadChannels_) {
0102     ecalChannelStatusToken_ = esConsumes<EcalChannelStatus, EcalChannelStatusRcd>();
0103   }
0104
0105   std::string componentType = ps.getParameter<std::string>("algo");
0106   edm::ConsumesCollector c{consumesCollector()};
0107   worker_ = EcalRecHitWorkerFactory::get()->create(componentType, ps, c);
0108
0109   // to recover problematic channels
0110   componentType = ps.getParameter<std::string>("algoRecover");
0111   workerRecover_ = EcalRecHitWorkerFactory::get()->create(componentType, ps, c);
0112
0113   edm::ParameterSet cleaningPs = ps.getParameter<edm::ParameterSet>("cleaningConfig");
0114   cleaningAlgo_ = std::make_unique<EcalCleaningAlgo>(cleaningPs);
0115 }
0116
0117 EcalRecHitProducer::~EcalRecHitProducer() = default;
0118
0119 void EcalRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
0120   using namespace edm;
0121
0122   // collection of rechits to put in the event
0123   auto ebRecHits = std::make_unique<EBRecHitCollection>();
0124   auto eeRecHits = std::make_unique<EERecHitCollection>();
0125
0126   worker_->set(es);
0127
0128   if (recoverEBIsolatedChannels_ || recoverEEIsolatedChannels_ || recoverEBFE_ || recoverEEFE_ || recoverEBVFE_ ||
0129       recoverEEVFE_ || killDeadChannels_) {
0130     workerRecover_->set(es);
0131   }
0132
0133   // Make EB rechits
0134   if (doEB_) {
0135     const auto& ebUncalibRecHits = evt.get(ebUncalibRecHitToken_);
0136     LogDebug("EcalRecHitDebug") << "total # EB uncalibrated rechits: " << ebUncalibRecHits.size();
0137
0138     // loop over uncalibrated rechits to make calibrated ones
0139     for (const auto& uncalibRecHit : ebUncalibRecHits) {
0140       worker_->run(evt, uncalibRecHit, *ebRecHits);
0141     }
0142   }
0143
0144   // Make EE rechits
0145   if (doEE_) {
0146     const auto& eeUncalibRecHits = evt.get(eeUncalibRecHitToken_);
0147     LogDebug("EcalRecHitDebug") << "total # EE uncalibrated rechits: " << eeUncalibRecHits.size();
0148
0149     // loop over uncalibrated rechits to make calibrated ones
0150     for (const auto& uncalibRecHit : eeUncalibRecHits) {
0151       worker_->run(evt, uncalibRecHit, *eeRecHits);
0152     }
0153   }
0154
0155   // sort collections before attempting recovery, to avoid insertion of double recHits
0156   ebRecHits->sort();
0157   eeRecHits->sort();
0158
0159   if (recoverEBIsolatedChannels_ || recoverEBFE_ || killDeadChannels_) {
0160     const auto& detIds = evt.get(ebDetIdToBeRecoveredToken_);
0161     const auto& chStatus = es.getData(ecalChannelStatusToken_);
0162
0163     for (const auto& detId : detIds) {
0164       // get channel status map to treat dead VFE separately
0165       EcalChannelStatusMap::const_iterator chit = chStatus.find(detId);
0166       EcalChannelStatusCode chStatusCode;
0167       if (chit != chStatus.end()) {
0168         chStatusCode = *chit;
0169       } else {
0170         edm::LogError("EcalRecHitProducerError") << "No channel status found for xtal " << detId.rawId()
0171                                                  << "! something wrong with EcalChannelStatus in your DB? ";
0172       }
0173       EcalUncalibratedRecHit urh;
0174       if (chStatusCode.getStatusCode() == EcalChannelStatusCode::kDeadVFE) {  // dead VFE (from DB info)
0175         // uses the EcalUncalibratedRecHit to pass the DetId info
0176         urh = EcalUncalibratedRecHit(detId, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_VFE);
0177         if (recoverEBVFE_ || killDeadChannels_)
0178           workerRecover_->run(evt, urh, *ebRecHits);
0179       } else {
0180         // uses the EcalUncalibratedRecHit to pass the DetId info
0181         urh = EcalUncalibratedRecHit(detId, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_single);
0182         if (recoverEBIsolatedChannels_ || killDeadChannels_)
0183           workerRecover_->run(evt, urh, *ebRecHits);
0184       }
0185     }
0186   }
0187
0188   if (recoverEEIsolatedChannels_ || recoverEEVFE_ || killDeadChannels_) {
0189     const auto& detIds = evt.get(eeDetIdToBeRecoveredToken_);
0190     const auto& chStatus = es.getData(ecalChannelStatusToken_);
0191
0192     for (const auto& detId : detIds) {
0193       // get channel status map to treat dead VFE separately
0194       EcalChannelStatusMap::const_iterator chit = chStatus.find(detId);
0195       EcalChannelStatusCode chStatusCode;
0196       if (chit != chStatus.end()) {
0197         chStatusCode = *chit;
0198       } else {
0199         edm::LogError("EcalRecHitProducerError") << "No channel status found for xtal " << detId.rawId()
0200                                                  << "! something wrong with EcalChannelStatus in your DB? ";
0201       }
0202       EcalUncalibratedRecHit urh;
0203       if (chStatusCode.getStatusCode() == EcalChannelStatusCode::kDeadVFE) {  // dead VFE (from DB info)
0204         // uses the EcalUncalibratedRecHit to pass the DetId info
0205         urh = EcalUncalibratedRecHit(detId, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EE_VFE);
0206         if (recoverEEVFE_ || killDeadChannels_)
0207           workerRecover_->run(evt, urh, *eeRecHits);
0208       } else {
0209         // uses the EcalUncalibratedRecHit to pass the DetId info
0210         urh = EcalUncalibratedRecHit(detId, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EE_single);
0211         if (recoverEEIsolatedChannels_ || killDeadChannels_)
0212           workerRecover_->run(evt, urh, *eeRecHits);
0213       }
0214     }
0215   }
0216
0217   if (recoverEBFE_ || killDeadChannels_) {
0218     const auto& ttIds = evt.get(ebFEToBeRecoveredToken_);
0219
0220     for (const auto& ttId : ttIds) {
0221       // uses the EcalUncalibratedRecHit to pass the DetId info
0222       int ieta = ((ttId.ietaAbs() - 1) * 5 + 1) * ttId.zside();  // from EcalTrigTowerConstituentsMap
0223       int iphi = ((ttId.iphi() - 1) * 5 + 11) % 360;             // from EcalTrigTowerConstituentsMap
0224       if (iphi <= 0)
0225         iphi += 360;  // from EcalTrigTowerConstituentsMap
0226       EcalUncalibratedRecHit urh(
0227           EBDetId(ieta, iphi, EBDetId::ETAPHIMODE), 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_FE);
0228       workerRecover_->run(evt, urh, *ebRecHits);
0229     }
0230   }
0231
0232   if (recoverEEFE_ || killDeadChannels_) {
0233     const auto& scIds = evt.get(eeFEToBeRecoveredToken_);
0234
0235     for (const auto& scId : scIds) {
0236       // uses the EcalUncalibratedRecHit to pass the DetId info
0237       if (EEDetId::validDetId((scId.ix() - 1) * 5 + 1, (scId.iy() - 1) * 5 + 1, scId.zside())) {
0238         EcalUncalibratedRecHit urh(EEDetId((scId.ix() - 1) * 5 + 1, (scId.iy() - 1) * 5 + 1, scId.zside()),
0239                                    0,
0240                                    0,
0241                                    0,
0242                                    0,
0243                                    EcalRecHitWorkerBaseClass::EE_FE);
0244         workerRecover_->run(evt, urh, *eeRecHits);
0245       }
0246     }
0247   }
0248
0249   // without re-sorting, find (used below in cleaning) will lead
0250   // to undefined results
0251   ebRecHits->sort();
0252   eeRecHits->sort();
0253
0254   // apply spike cleaning
0255   if (cleaningAlgo_) {
0256     cleaningAlgo_->setFlags(*ebRecHits);
0257     cleaningAlgo_->setFlags(*eeRecHits);
0258   }
0259
0260   // put the collection of reconstructed hits in the event
0261   LogInfo("EcalRecHitInfo") << "total # EB calibrated rechits: " << ebRecHits->size();
0262   LogInfo("EcalRecHitInfo") << "total # EE calibrated rechits: " << eeRecHits->size();
0263
0264   evt.put(ebRecHitToken_, std::move(ebRecHits));
0265   evt.put(eeRecHitToken_, std::move(eeRecHits));
0266 }
0267
0268 void EcalRecHitProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0269   edm::ParameterSetDescription desc;
0270   desc.add<bool>("recoverEEVFE", false);
0271   desc.add<std::string>("EErechitCollection", "EcalRecHitsEE");
0272   desc.add<bool>("recoverEBIsolatedChannels", false);
0273   desc.add<bool>("recoverEBVFE", false);
0274   desc.add<bool>("laserCorrection", true);
0275   desc.add<double>("EBLaserMIN", 0.5);
0276   desc.add<bool>("killDeadChannels", true);
0277   {
0278     std::vector<int> temp1;
0279     temp1.reserve(3);
0280     temp1.push_back(14);
0281     temp1.push_back(78);
0282     temp1.push_back(142);
0283     desc.add<std::vector<int>>("dbStatusToBeExcludedEB", temp1);
0284   }
0285   desc.add<edm::InputTag>("EEuncalibRecHitCollection",
0286                           edm::InputTag("ecalMultiFitUncalibRecHit", "EcalUncalibRecHitsEE"));
0287   {
0288     std::vector<int> temp1;
0289     temp1.reserve(3);
0290     temp1.push_back(14);
0291     temp1.push_back(78);
0292     temp1.push_back(142);
0293     desc.add<std::vector<int>>("dbStatusToBeExcludedEE", temp1);
0294   }
0295   desc.add<double>("EELaserMIN", 0.5);
0296   desc.add<edm::InputTag>("ebFEToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "ebFE"));
0297   {
0298     edm::ParameterSetDescription psd0;
0299     psd0.add<double>("e6e2thresh", 0.04);
0300     psd0.add<double>("tightenCrack_e6e2_double", 3);
0301     psd0.add<double>("e4e1Threshold_endcap", 0.3);
0302     psd0.add<double>("tightenCrack_e4e1_single", 3);
0303     psd0.add<double>("tightenCrack_e1_double", 2);
0304     psd0.add<double>("cThreshold_barrel", 4);
0305     psd0.add<double>("e4e1Threshold_barrel", 0.08);
0306     psd0.add<double>("tightenCrack_e1_single", 2);
0307     psd0.add<double>("e4e1_b_barrel", -0.024);
0308     psd0.add<double>("e4e1_a_barrel", 0.04);
0309     psd0.add<double>("ignoreOutOfTimeThresh", 1000000000.0);
0310     psd0.add<double>("cThreshold_endcap", 15);
0311     psd0.add<double>("e4e1_b_endcap", -0.0125);
0312     psd0.add<double>("e4e1_a_endcap", 0.02);
0313     psd0.add<double>("cThreshold_double", 10);
0314     desc.add<edm::ParameterSetDescription>("cleaningConfig", psd0);
0315   }
0316   desc.add<double>("logWarningEtThreshold_EE_FE", 50);
0317   desc.add<edm::InputTag>("eeDetIdToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "eeDetId"));
0318   desc.add<bool>("recoverEBFE", true);
0319   desc.add<edm::InputTag>("eeFEToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "eeFE"));
0320   desc.add<edm::InputTag>("ebDetIdToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "ebDetId"));
0321   desc.add<double>("singleChannelRecoveryThreshold", 8);
0322   desc.add<double>("sum8ChannelRecoveryThreshold", 0.);
0323   desc.add<edm::FileInPath>("bdtWeightFileNoCracks",
0324                             edm::FileInPath("RecoLocalCalo/EcalDeadChannelRecoveryAlgos/data/BDTWeights/"
0325                                             "bdtgAllRH_8GT700MeV_noCracks_ZskimData2017_v1.xml"));
0326   desc.add<edm::FileInPath>("bdtWeightFileCracks",
0327                             edm::FileInPath("RecoLocalCalo/EcalDeadChannelRecoveryAlgos/data/BDTWeights/"
0328                                             "bdtgAllRH_8GT700MeV_onlyCracks_ZskimData2017_v1.xml"));
0329   {
0330     std::vector<std::string> temp1;
0331     temp1.reserve(9);
0332     temp1.push_back("kNoisy");
0333     temp1.push_back("kNNoisy");
0334     temp1.push_back("kFixedG6");
0335     temp1.push_back("kFixedG1");
0336     temp1.push_back("kFixedG0");
0337     temp1.push_back("kNonRespondingIsolated");
0338     temp1.push_back("kDeadVFE");
0339     temp1.push_back("kDeadFE");
0340     temp1.push_back("kNoDataNoTP");
0341     desc.add<std::vector<std::string>>("ChannelStatusToBeExcluded", temp1);
0342   }
0343   desc.add<std::string>("EBrechitCollection", "EcalRecHitsEB");
0344   desc.add<edm::InputTag>("triggerPrimitiveDigiCollection", edm::InputTag("ecalDigis", "EcalTriggerPrimitives"));
0345   desc.add<bool>("recoverEEFE", true);
0346   desc.add<std::string>("singleChannelRecoveryMethod", "NeuralNetworks");
0347   desc.add<double>("EBLaserMAX", 3.0);
0348   {
0349     edm::ParameterSetDescription psd0;
0350     {
0351       std::vector<std::string> temp2;
0352       temp2.reserve(4);
0353       temp2.push_back("kOk");
0354       temp2.push_back("kDAC");
0355       temp2.push_back("kNoLaser");
0356       temp2.push_back("kNoisy");
0357       psd0.add<std::vector<std::string>>("kGood", temp2);
0358     }
0359     {
0360       std::vector<std::string> temp2;
0361       temp2.reserve(3);
0362       temp2.push_back("kFixedG0");
0363       temp2.push_back("kNonRespondingIsolated");
0364       temp2.push_back("kDeadVFE");
0365       psd0.add<std::vector<std::string>>("kNeighboursRecovered", temp2);
0366     }
0367     {
0368       std::vector<std::string> temp2;
0369       temp2.reserve(1);
0370       temp2.push_back("kNoDataNoTP");
0371       psd0.add<std::vector<std::string>>("kDead", temp2);
0372     }
0373     {
0374       std::vector<std::string> temp2;
0375       temp2.reserve(3);
0376       temp2.push_back("kNNoisy");
0377       temp2.push_back("kFixedG6");
0378       temp2.push_back("kFixedG1");
0379       psd0.add<std::vector<std::string>>("kNoisy", temp2);
0380     }
0381     {
0382       std::vector<std::string> temp2;
0383       temp2.reserve(1);
0384       temp2.push_back("kDeadFE");
0385       psd0.add<std::vector<std::string>>("kTowerRecovered", temp2);
0386     }
0387     desc.add<edm::ParameterSetDescription>("flagsMapDBReco", psd0);
0388   }
0389   desc.add<edm::InputTag>("EBuncalibRecHitCollection",
0390                           edm::InputTag("ecalMultiFitUncalibRecHit", "EcalUncalibRecHitsEB"));
0391   desc.add<std::string>("algoRecover", "EcalRecHitWorkerRecover");
0392   desc.add<std::string>("algo", "EcalRecHitWorkerSimple");
0393   desc.add<double>("EELaserMAX", 8.0);
0394   desc.add<double>("logWarningEtThreshold_EB_FE", 50);
0395   desc.add<bool>("recoverEEIsolatedChannels", false);
0396   desc.add<edm::ESInputTag>("timeCalibTag", edm::ESInputTag());
0397   desc.add<edm::ESInputTag>("timeOffsetTag", edm::ESInputTag());
0398   desc.add<bool>("skipTimeCalib", false);
0399   descriptions.add("ecalRecHit", desc);
0400 }
0401
0402 #include "FWCore/Framework/interface/MakerMacros.h"
0403 DEFINE_FWK_MODULE(EcalRecHitProducer);