HcalRecProducers/src/HFPhase1Reconstructor.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoLocalCalo/HcalRecProducers
0004 // Class:      HFPhase1Reconstructor
0005 //
0006 /**\class HFPhase1Reconstructor HFPhase1Reconstructor.cc RecoLocalCalo/HcalRecProducers/src/HFPhase1Reconstructor.cc
0007
0008  Description: Phase 1 HF reco with QIE 10 and split-anode readout
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Igor Volobouev
0015 //         Created:  Thu, 25 May 2016 00:17:51 GMT
0016 //
0017 //
0018
0019 // system include files
0020
0021 // user include files
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023 #include "FWCore/Framework/interface/stream/EDProducer.h"
0024 #include "FWCore/Utilities/interface/ESGetToken.h"
0025 #include "FWCore/Framework/interface/Event.h"
0026 #include "FWCore/Framework/interface/MakerMacros.h"
0027 #include "FWCore/Utilities/interface/Exception.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/Utilities/interface/StreamID.h"
0031
0032 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0033
0034 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0035 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
0036 #include "CondFormats/HcalObjects/interface/HFPhase1PMTParams.h"
0037 #include "CondFormats/DataRecord/interface/HFPhase1PMTParamsRcd.h"
0038
0039 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
0040 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
0041
0042 // Noise cleanup algos
0043 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalHF_PETalgorithm.h"
0044 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalHF_S9S1algorithm.h"
0045 #include "RecoLocalCalo/HcalRecAlgos/interface/HFStripFilter.h"
0046
0047 // Parser for Phase 1 HF reco algorithms
0048 #include "RecoLocalCalo/HcalRecAlgos/interface/parseHFPhase1AlgoDescription.h"
0049
0050 //
0051 // class declaration
0052 //
0053 class HFPhase1Reconstructor : public edm::stream::EDProducer<> {
0054 public:
0055   explicit HFPhase1Reconstructor(const edm::ParameterSet&);
0056   ~HFPhase1Reconstructor() override;
0057
0058   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0059
0060 private:
0061   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0062   void produce(edm::Event&, const edm::EventSetup&) override;
0063
0064   // Configuration parameters
0065   std::string algoConfigClass_;
0066   bool setNoiseFlags_;
0067   bool runHFStripFilter_;
0068   bool useChannelQualityFromDB_;
0069   bool checkChannelQualityForDepth3and4_;
0070
0071   // Other members
0072   edm::EDGetTokenT<HFPreRecHitCollection> tok_PreRecHit_;
0073   std::unique_ptr<AbsHFPhase1Algo> reco_;
0074   std::unique_ptr<AbsHcalAlgoData> recoConfig_;
0075   edm::ESGetToken<HFPhase1PMTParams, HFPhase1PMTParamsRcd> recoConfigToken_;
0076
0077   // Noise cleanup algos
0078   std::unique_ptr<HcalHF_S9S1algorithm> hfS9S1_;
0079   std::unique_ptr<HcalHF_S9S1algorithm> hfS8S1_;
0080   std::unique_ptr<HcalHF_PETalgorithm> hfPET_;
0081   std::unique_ptr<HFStripFilter> hfStripFilter_;
0082
0083   // ES tokens
0084   edm::ESGetToken<HcalDbService, HcalDbRecord> conditionsToken_;
0085   edm::ESGetToken<HcalChannelQuality, HcalChannelQualityRcd> qualToken_;
0086   edm::ESGetToken<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd> sevToken_;
0087 };
0088
0089 //
0090 // constructors and destructor
0091 //
0092 HFPhase1Reconstructor::HFPhase1Reconstructor(const edm::ParameterSet& conf)
0093     : algoConfigClass_(conf.getParameter<std::string>("algoConfigClass")),
0094       setNoiseFlags_(conf.getParameter<bool>("setNoiseFlags")),
0095       runHFStripFilter_(conf.getParameter<bool>("runHFStripFilter")),
0096       useChannelQualityFromDB_(conf.getParameter<bool>("useChannelQualityFromDB")),
0097       checkChannelQualityForDepth3and4_(conf.getParameter<bool>("checkChannelQualityForDepth3and4")),
0098       reco_(parseHFPhase1AlgoDescription(conf.getParameter<edm::ParameterSet>("algorithm"))) {
0099   // Check that the reco algorithm has been successfully configured
0100   if (!reco_.get())
0101     throw cms::Exception("HFPhase1BadConfig") << "Invalid HFPhase1Reconstructor algorithm configuration" << std::endl;
0102
0103   if (reco_->isConfigurable()) {
0104     if ("HFPhase1PMTParams" != algoConfigClass_) {
0105       throw cms::Exception("HBHEPhase1BadConfig")
0106           << "Invalid HBHEPhase1Reconstructor \"algoConfigClass\" parameter value \"" << algoConfigClass_ << '"'
0107           << std::endl;
0108     }
0109     recoConfigToken_ = esConsumes<edm::Transition::BeginRun>();
0110   }
0111
0112   // Configure the noise cleanup algorithms
0113   if (setNoiseFlags_) {
0114     const edm::ParameterSet& psS9S1 = conf.getParameter<edm::ParameterSet>("S9S1stat");
0115     hfS9S1_ = std::make_unique<HcalHF_S9S1algorithm>(psS9S1.getParameter<std::vector<double> >("short_optimumSlope"),
0116                                                      psS9S1.getParameter<std::vector<double> >("shortEnergyParams"),
0117                                                      psS9S1.getParameter<std::vector<double> >("shortETParams"),
0118                                                      psS9S1.getParameter<std::vector<double> >("long_optimumSlope"),
0119                                                      psS9S1.getParameter<std::vector<double> >("longEnergyParams"),
0120                                                      psS9S1.getParameter<std::vector<double> >("longETParams"),
0121                                                      psS9S1.getParameter<int>("HcalAcceptSeverityLevel"),
0122                                                      psS9S1.getParameter<bool>("isS8S1"));
0123
0124     const edm::ParameterSet& psS8S1 = conf.getParameter<edm::ParameterSet>("S8S1stat");
0125     hfS8S1_ = std::make_unique<HcalHF_S9S1algorithm>(psS8S1.getParameter<std::vector<double> >("short_optimumSlope"),
0126                                                      psS8S1.getParameter<std::vector<double> >("shortEnergyParams"),
0127                                                      psS8S1.getParameter<std::vector<double> >("shortETParams"),
0128                                                      psS8S1.getParameter<std::vector<double> >("long_optimumSlope"),
0129                                                      psS8S1.getParameter<std::vector<double> >("longEnergyParams"),
0130                                                      psS8S1.getParameter<std::vector<double> >("longETParams"),
0131                                                      psS8S1.getParameter<int>("HcalAcceptSeverityLevel"),
0132                                                      psS8S1.getParameter<bool>("isS8S1"));
0133
0134     const edm::ParameterSet& psPET = conf.getParameter<edm::ParameterSet>("PETstat");
0135     hfPET_ = std::make_unique<HcalHF_PETalgorithm>(psPET.getParameter<std::vector<double> >("short_R"),
0136                                                    psPET.getParameter<std::vector<double> >("shortEnergyParams"),
0137                                                    psPET.getParameter<std::vector<double> >("shortETParams"),
0138                                                    psPET.getParameter<std::vector<double> >("long_R"),
0139                                                    psPET.getParameter<std::vector<double> >("longEnergyParams"),
0140                                                    psPET.getParameter<std::vector<double> >("longETParams"),
0141                                                    psPET.getParameter<int>("HcalAcceptSeverityLevel"),
0142                                                    psPET.getParameter<std::vector<double> >("short_R_29"),
0143                                                    psPET.getParameter<std::vector<double> >("long_R_29"));
0144
0145     // Configure HFStripFilter
0146     if (runHFStripFilter_) {
0147       const edm::ParameterSet& psStripFilter = conf.getParameter<edm::ParameterSet>("HFStripFilter");
0148       hfStripFilter_ = HFStripFilter::parseParameterSet(psStripFilter);
0149     }
0150   }
0151
0152   // Describe consumed data
0153   tok_PreRecHit_ = consumes<HFPreRecHitCollection>(conf.getParameter<edm::InputTag>("inputLabel"));
0154
0155   // Register the product
0156   produces<HFRecHitCollection>();
0157
0158   // ES tokens
0159   conditionsToken_ = esConsumes<HcalDbService, HcalDbRecord>();
0160   qualToken_ = esConsumes<HcalChannelQuality, HcalChannelQualityRcd>(edm::ESInputTag("", "withTopo"));
0161   sevToken_ = esConsumes<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd>();
0162 }
0163
0164 HFPhase1Reconstructor::~HFPhase1Reconstructor() {
0165   // do anything here that needs to be done at destruction time
0166   // (e.g. close files, deallocate resources etc.)
0167 }
0168
0169 void HFPhase1Reconstructor::beginRun(const edm::Run& r, const edm::EventSetup& es) {
0170   if (reco_->isConfigurable()) {
0171     recoConfig_ = std::make_unique<HFPhase1PMTParams>(es.getData(recoConfigToken_));
0172     if (!reco_->configure(recoConfig_.get()))
0173       throw cms::Exception("HFPhase1BadConfig")
0174           << "Failed to configure HFPhase1Reconstructor algorithm from EventSetup" << std::endl;
0175   }
0176 }
0177
0178 // ------------ method called to produce the data  ------------
0179 void HFPhase1Reconstructor::produce(edm::Event& e, const edm::EventSetup& eventSetup) {
0180   using namespace edm;
0181
0182   // Fetch the calibrations
0183   const HcalDbService* conditions = &eventSetup.getData(conditionsToken_);
0184   const HcalChannelQuality* myqual = &eventSetup.getData(qualToken_);
0185   const HcalSeverityLevelComputer* mySeverity = &eventSetup.getData(sevToken_);
0186
0187   // Get the input data
0188   Handle<HFPreRecHitCollection> preRecHits;
0189   e.getByToken(tok_PreRecHit_, preRecHits);
0190
0191   // Create a new output collection
0192   std::unique_ptr<HFRecHitCollection> rec(std::make_unique<HFRecHitCollection>());
0193   rec->reserve(preRecHits->size());
0194
0195   // Iterate over the input and fill the output collection
0196   for (HFPreRecHitCollection::const_iterator it = preRecHits->begin(); it != preRecHits->end(); ++it) {
0197     // The check whether this PMT is single-anode one should go here.
0198     // Fix this piece of code if/when mixed-anode readout configurations
0199     // become available.
0200     const bool thisIsSingleAnodePMT = false;
0201
0202     // Check if the anodes were tagged bad in the database
0203     bool taggedBadByDb[2] = {false, false};
0204     if (useChannelQualityFromDB_) {
0205       if (checkChannelQualityForDepth3and4_ && !thisIsSingleAnodePMT) {
0206         HcalDetId anodeIds[2];
0207         anodeIds[0] = it->id();
0208         anodeIds[1] = anodeIds[0].secondAnodeId();
0209         for (unsigned i = 0; i < 2; ++i) {
0210           const HcalChannelStatus* mydigistatus = myqual->getValues(anodeIds[i].rawId());
0211           taggedBadByDb[i] = mySeverity->dropChannel(mydigistatus->getValue());
0212         }
0213       } else {
0214         const HcalChannelStatus* mydigistatus = myqual->getValues(it->id().rawId());
0215         const bool b = mySeverity->dropChannel(mydigistatus->getValue());
0216         taggedBadByDb[0] = b;
0217         taggedBadByDb[1] = b;
0218       }
0219     }
0220
0221     // Reconstruct the rechit
0222     const HFRecHit& rh =
0223         reco_->reconstruct(*it, conditions->getHcalCalibrations(it->id()), taggedBadByDb, thisIsSingleAnodePMT);
0224
0225     // The rechit will have the id of 0 if the algorithm
0226     // decides that it should be dropped
0227     if (rh.id().rawId())
0228       rec->push_back(rh);
0229   }
0230
0231   // At this point, the rechits contain energy, timing,
0232   // as well as proper auxiliary words. We do still need
0233   // to set certain flags using the noise cleanup algoritms.
0234
0235   // The following flags require the full set of rechits.
0236   // These flags need to be set consecutively.
0237   if (setNoiseFlags_) {
0238     // Step 1:  Set PET flag  (short fibers of |ieta|==29)
0239     for (HFRecHitCollection::iterator i = rec->begin(); i != rec->end(); ++i) {
0240       int depth = i->id().depth();
0241       int ieta = i->id().ieta();
0242       // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
0243       if (depth == 2 || abs(ieta) == 29)
0244         hfPET_->HFSetFlagFromPET(*i, *rec, myqual, mySeverity);
0245     }
0246
0247     // Step 2:  Set S8S1 flag (short fibers or |ieta|==29)
0248     for (HFRecHitCollection::iterator i = rec->begin(); i != rec->end(); ++i) {
0249       int depth = i->id().depth();
0250       int ieta = i->id().ieta();
0251       // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
0252       if (depth == 2 || abs(ieta) == 29)
0253         hfS8S1_->HFSetFlagFromS9S1(*i, *rec, myqual, mySeverity);
0254     }
0255
0256     // Step 3:  Set S9S1 flag (long fibers)
0257     for (HFRecHitCollection::iterator i = rec->begin(); i != rec->end(); ++i) {
0258       int depth = i->id().depth();
0259       int ieta = i->id().ieta();
0260       // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
0261       if (depth == 1 && abs(ieta) != 29)
0262         hfS9S1_->HFSetFlagFromS9S1(*i, *rec, myqual, mySeverity);
0263     }
0264
0265     // Step 4:  Run HFStripFilter if requested
0266     if (runHFStripFilter_)
0267       hfStripFilter_->runFilter(*rec, myqual);
0268   }
0269
0270   // Add the output collection to the event record
0271   e.put(std::move(rec));
0272 }
0273
0274 #define add_param_set(name) /**/     \
0275   edm::ParameterSetDescription name; \
0276   name.setAllowAnything();           \
0277   desc.add<edm::ParameterSetDescription>(#name, name)
0278
0279 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0280 void HFPhase1Reconstructor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0281   edm::ParameterSetDescription desc;
0282
0283   desc.add<edm::InputTag>("inputLabel");
0284   desc.add<std::string>("algoConfigClass");
0285   desc.add<bool>("setNoiseFlags");
0286   desc.add<bool>("runHFStripFilter", false);
0287   desc.add<bool>("useChannelQualityFromDB");
0288   desc.add<bool>("checkChannelQualityForDepth3and4");
0289   desc.add<edm::ParameterSetDescription>("algorithm", fillDescriptionForParseHFPhase1AlgoDescription());
0290   desc.add<edm::ParameterSetDescription>("HFStripFilter", HFStripFilter::fillDescription());
0291
0292   add_param_set(S9S1stat);
0293   add_param_set(S8S1stat);
0294   add_param_set(PETstat);
0295
0296   descriptions.addDefault(desc);
0297 }
0298
0299 //define this as a plug-in
0300 DEFINE_FWK_MODULE(HFPhase1Reconstructor);