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/Event.h"
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/Framework/interface/stream/EDProducer.h"
0026 #include "FWCore/ParameterSet/interface/EmptyGroupDescription.h"
0027 #include "FWCore/Utilities/interface/ESGetToken.h"
0028 #include "FWCore/Utilities/interface/Exception.h"
0029
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/Utilities/interface/StreamID.h"
0032
0033 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0034
0035 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0036 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
0037 #include "CondFormats/HcalObjects/interface/HFPhase1PMTParams.h"
0038 #include "CondFormats/DataRecord/interface/HFPhase1PMTParamsRcd.h"
0039
0040 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputer.h"
0041 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalSeverityLevelComputerRcd.h"
0042
0043 // Noise cleanup algos
0044 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalHF_PETalgorithm.h"
0045 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalHF_S9S1algorithm.h"
0046 #include "RecoLocalCalo/HcalRecAlgos/interface/HFStripFilter.h"
0047
0048 // Parser for Phase 1 HF reco algorithms
0049 #include "RecoLocalCalo/HcalRecAlgos/interface/parseHFPhase1AlgoDescription.h"
0050
0051 //
0052 // class declaration
0053 //
0054 class HFPhase1Reconstructor : public edm::stream::EDProducer<> {
0055 public:
0056   explicit HFPhase1Reconstructor(const edm::ParameterSet&);
0057   ~HFPhase1Reconstructor() override = default;
0058
0059   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0060
0061 private:
0062   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0063   void produce(edm::Event&, const edm::EventSetup&) override;
0064
0065   // Configuration parameters
0066   std::string algoConfigClass_;
0067   bool setNoiseFlags_;
0068   bool runHFStripFilter_;
0069   bool useChannelQualityFromDB_;
0070   bool checkChannelQualityForDepth3and4_;
0071
0072   // Other members
0073   edm::EDGetTokenT<HFPreRecHitCollection> tok_PreRecHit_;
0074   std::unique_ptr<AbsHFPhase1Algo> reco_;
0075   std::unique_ptr<AbsHcalAlgoData> recoConfig_;
0076   edm::ESGetToken<HFPhase1PMTParams, HFPhase1PMTParamsRcd> recoConfigToken_;
0077
0078   // Noise cleanup algos
0079   std::unique_ptr<HcalHF_S9S1algorithm> hfS9S1_;
0080   std::unique_ptr<HcalHF_S9S1algorithm> hfS8S1_;
0081   std::unique_ptr<HcalHF_PETalgorithm> hfPET_;
0082   std::unique_ptr<HFStripFilter> hfStripFilter_;
0083
0084   // ES tokens
0085   edm::ESGetToken<HcalDbService, HcalDbRecord> conditionsToken_;
0086   edm::ESGetToken<HcalChannelQuality, HcalChannelQualityRcd> qualToken_;
0087   edm::ESGetToken<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd> sevToken_;
0088 };
0089
0090 //
0091 // constructors and destructor
0092 //
0093 HFPhase1Reconstructor::HFPhase1Reconstructor(const edm::ParameterSet& conf)
0094     : algoConfigClass_(conf.getParameter<std::string>("algoConfigClass")),
0095       setNoiseFlags_(conf.getParameter<bool>("setNoiseFlags")),
0096       runHFStripFilter_(conf.getParameter<bool>("runHFStripFilter")),
0097       useChannelQualityFromDB_(conf.getParameter<bool>("useChannelQualityFromDB")),
0098       checkChannelQualityForDepth3and4_(conf.getParameter<bool>("checkChannelQualityForDepth3and4")),
0099       reco_(parseHFPhase1AlgoDescription(conf.getParameter<edm::ParameterSet>("algorithm"))) {
0100   // Check that the reco algorithm has been successfully configured
0101   if (!reco_.get())
0102     throw cms::Exception("HFPhase1BadConfig") << "Invalid HFPhase1Reconstructor algorithm configuration" << std::endl;
0103
0104   if (reco_->isConfigurable()) {
0105     if ("HFPhase1PMTParams" != algoConfigClass_) {
0106       throw cms::Exception("HBHEPhase1BadConfig")
0107           << "Invalid HBHEPhase1Reconstructor \"algoConfigClass\" parameter value \"" << algoConfigClass_ << '"'
0108           << std::endl;
0109     }
0110     recoConfigToken_ = esConsumes<edm::Transition::BeginRun>();
0111   }
0112
0113   // Configure the noise cleanup algorithms
0114   if (setNoiseFlags_) {
0115     const edm::ParameterSet& psS9S1 = conf.getParameter<edm::ParameterSet>("S9S1stat");
0116     hfS9S1_ = std::make_unique<HcalHF_S9S1algorithm>(psS9S1.getParameter<std::vector<double>>("short_optimumSlope"),
0117                                                      psS9S1.getParameter<std::vector<double>>("shortEnergyParams"),
0118                                                      psS9S1.getParameter<std::vector<double>>("shortETParams"),
0119                                                      psS9S1.getParameter<std::vector<double>>("long_optimumSlope"),
0120                                                      psS9S1.getParameter<std::vector<double>>("longEnergyParams"),
0121                                                      psS9S1.getParameter<std::vector<double>>("longETParams"),
0122                                                      psS9S1.getParameter<int>("HcalAcceptSeverityLevel"),
0123                                                      psS9S1.getParameter<bool>("isS8S1"));
0124
0125     const edm::ParameterSet& psS8S1 = conf.getParameter<edm::ParameterSet>("S8S1stat");
0126     hfS8S1_ = std::make_unique<HcalHF_S9S1algorithm>(psS8S1.getParameter<std::vector<double>>("short_optimumSlope"),
0127                                                      psS8S1.getParameter<std::vector<double>>("shortEnergyParams"),
0128                                                      psS8S1.getParameter<std::vector<double>>("shortETParams"),
0129                                                      psS8S1.getParameter<std::vector<double>>("long_optimumSlope"),
0130                                                      psS8S1.getParameter<std::vector<double>>("longEnergyParams"),
0131                                                      psS8S1.getParameter<std::vector<double>>("longETParams"),
0132                                                      psS8S1.getParameter<int>("HcalAcceptSeverityLevel"),
0133                                                      psS8S1.getParameter<bool>("isS8S1"));
0134
0135     const edm::ParameterSet& psPET = conf.getParameter<edm::ParameterSet>("PETstat");
0136     hfPET_ = std::make_unique<HcalHF_PETalgorithm>(psPET.getParameter<std::vector<double>>("short_R"),
0137                                                    psPET.getParameter<std::vector<double>>("shortEnergyParams"),
0138                                                    psPET.getParameter<std::vector<double>>("shortETParams"),
0139                                                    psPET.getParameter<std::vector<double>>("long_R"),
0140                                                    psPET.getParameter<std::vector<double>>("longEnergyParams"),
0141                                                    psPET.getParameter<std::vector<double>>("longETParams"),
0142                                                    psPET.getParameter<int>("HcalAcceptSeverityLevel"),
0143                                                    psPET.getParameter<std::vector<double>>("short_R_29"),
0144                                                    psPET.getParameter<std::vector<double>>("long_R_29"));
0145
0146     // Configure HFStripFilter
0147     if (runHFStripFilter_) {
0148       const edm::ParameterSet& psStripFilter = conf.getParameter<edm::ParameterSet>("HFStripFilter");
0149       hfStripFilter_ = HFStripFilter::parseParameterSet(psStripFilter);
0150     }
0151   }
0152
0153   // Describe consumed data
0154   tok_PreRecHit_ = consumes<HFPreRecHitCollection>(conf.getParameter<edm::InputTag>("inputLabel"));
0155
0156   // Register the product
0157   produces<HFRecHitCollection>();
0158
0159   // ES tokens
0160   conditionsToken_ = esConsumes<HcalDbService, HcalDbRecord>();
0161   qualToken_ = esConsumes<HcalChannelQuality, HcalChannelQualityRcd>(edm::ESInputTag("", "withTopo"));
0162   sevToken_ = esConsumes<HcalSeverityLevelComputer, HcalSeverityLevelComputerRcd>();
0163 }
0164
0165 void HFPhase1Reconstructor::beginRun(const edm::Run& r, const edm::EventSetup& es) {
0166   if (reco_->isConfigurable()) {
0167     recoConfig_ = std::make_unique<HFPhase1PMTParams>(es.getData(recoConfigToken_));
0168     if (!reco_->configure(recoConfig_.get()))
0169       throw cms::Exception("HFPhase1BadConfig")
0170           << "Failed to configure HFPhase1Reconstructor algorithm from EventSetup" << std::endl;
0171   }
0172 }
0173
0174 // ------------ method called to produce the data  ------------
0175 void HFPhase1Reconstructor::produce(edm::Event& e, const edm::EventSetup& eventSetup) {
0176   using namespace edm;
0177
0178   // Fetch the calibrations
0179   const HcalDbService* conditions = &eventSetup.getData(conditionsToken_);
0180   const HcalChannelQuality* myqual = &eventSetup.getData(qualToken_);
0181   const HcalSeverityLevelComputer* mySeverity = &eventSetup.getData(sevToken_);
0182
0183   // Get the input data
0184   Handle<HFPreRecHitCollection> preRecHits;
0185   e.getByToken(tok_PreRecHit_, preRecHits);
0186
0187   // Create a new output collection
0188   std::unique_ptr<HFRecHitCollection> rec(std::make_unique<HFRecHitCollection>());
0189   rec->reserve(preRecHits->size());
0190
0191   // Iterate over the input and fill the output collection
0192   for (HFPreRecHitCollection::const_iterator it = preRecHits->begin(); it != preRecHits->end(); ++it) {
0193     // The check whether this PMT is single-anode one should go here.
0194     // Fix this piece of code if/when mixed-anode readout configurations
0195     // become available.
0196     const bool thisIsSingleAnodePMT = false;
0197
0198     // Check if the anodes were tagged bad in the database
0199     bool taggedBadByDb[2] = {false, false};
0200     if (useChannelQualityFromDB_) {
0201       if (checkChannelQualityForDepth3and4_ && !thisIsSingleAnodePMT) {
0202         HcalDetId anodeIds[2];
0203         anodeIds[0] = it->id();
0204         anodeIds[1] = anodeIds[0].secondAnodeId();
0205         for (unsigned i = 0; i < 2; ++i) {
0206           const HcalChannelStatus* mydigistatus = myqual->getValues(anodeIds[i].rawId());
0207           taggedBadByDb[i] = mySeverity->dropChannel(mydigistatus->getValue());
0208         }
0209       } else {
0210         const HcalChannelStatus* mydigistatus = myqual->getValues(it->id().rawId());
0211         const bool b = mySeverity->dropChannel(mydigistatus->getValue());
0212         taggedBadByDb[0] = b;
0213         taggedBadByDb[1] = b;
0214       }
0215     }
0216
0217     // Reconstruct the rechit
0218     const HFRecHit& rh =
0219         reco_->reconstruct(*it, conditions->getHcalCalibrations(it->id()), taggedBadByDb, thisIsSingleAnodePMT);
0220
0221     // The rechit will have the id of 0 if the algorithm
0222     // decides that it should be dropped
0223     if (rh.id().rawId())
0224       rec->push_back(rh);
0225   }
0226
0227   // At this point, the rechits contain energy, timing,
0228   // as well as proper auxiliary words. We do still need
0229   // to set certain flags using the noise cleanup algoritms.
0230
0231   // The following flags require the full set of rechits.
0232   // These flags need to be set consecutively.
0233   if (setNoiseFlags_) {
0234     // Step 1:  Set PET flag  (short fibers of |ieta|==29)
0235     for (HFRecHitCollection::iterator i = rec->begin(); i != rec->end(); ++i) {
0236       int depth = i->id().depth();
0237       int ieta = i->id().ieta();
0238       // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
0239       if (depth == 2 || abs(ieta) == 29)
0240         hfPET_->HFSetFlagFromPET(*i, *rec, myqual, mySeverity);
0241     }
0242
0243     // Step 2:  Set S8S1 flag (short fibers or |ieta|==29)
0244     for (HFRecHitCollection::iterator i = rec->begin(); i != rec->end(); ++i) {
0245       int depth = i->id().depth();
0246       int ieta = i->id().ieta();
0247       // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
0248       if (depth == 2 || abs(ieta) == 29)
0249         hfS8S1_->HFSetFlagFromS9S1(*i, *rec, myqual, mySeverity);
0250     }
0251
0252     // Step 3:  Set S9S1 flag (long fibers)
0253     for (HFRecHitCollection::iterator i = rec->begin(); i != rec->end(); ++i) {
0254       int depth = i->id().depth();
0255       int ieta = i->id().ieta();
0256       // Short fibers and all channels at |ieta|=29 use PET settings in Algo 3
0257       if (depth == 1 && abs(ieta) != 29)
0258         hfS9S1_->HFSetFlagFromS9S1(*i, *rec, myqual, mySeverity);
0259     }
0260
0261     // Step 4:  Run HFStripFilter if requested
0262     if (runHFStripFilter_)
0263       hfStripFilter_->runFilter(*rec, myqual);
0264   }
0265
0266   // Add the output collection to the event record
0267   e.put(std::move(rec));
0268 }
0269
0270 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0271 void HFPhase1Reconstructor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0272   edm::ParameterSetDescription desc;
0273
0274   desc.add<edm::InputTag>("inputLabel", edm::InputTag("hfprereco"))
0275       ->setComment("Label for the input HFPreRecHitCollection");
0276   desc.add<std::string>("algoConfigClass", "HFPhase1PMTParams")
0277       ->setComment("reconstruction algorithm data to fetch from DB, if any");
0278   desc.add<bool>("useChannelQualityFromDB", true)
0279       ->setComment("change the following to True in order to use the channel status from the DB");
0280   desc.add<bool>("checkChannelQualityForDepth3and4", true);
0281   desc.add<edm::ParameterSetDescription>("algorithm", fillDescriptionForParseHFPhase1AlgoDescription())
0282       ->setComment("configure the reconstruction algorithm");
0283
0284   desc.add<bool>("runHFStripFilter", true);
0285   desc.add<edm::ParameterSetDescription>("HFStripFilter", HFStripFilter::fillDescription());
0286   desc.add<bool>("setNoiseFlags", true);
0287
0288   {
0289     // Define common vectors
0290     std::vector<double> slopes_S9S1_run1 = {-99999,
0291                                             0.0164905,
0292                                             0.0238698,
0293                                             0.0321383,
0294                                             0.041296,
0295                                             0.0513428,
0296                                             0.0622789,
0297                                             0.0741041,
0298                                             0.0868186,
0299                                             0.100422,
0300                                             0.135313,
0301                                             0.136289,
0302                                             0.0589927};
0303     std::vector<double> coeffs = {1.0, 2.5, 2.2, 2.0, 1.8, 1.6, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0};
0304     std::vector<double> slopes_S9S1_run2(slopes_S9S1_run1.size());
0305     for (size_t i = 0; i < slopes_S9S1_run1.size(); ++i) {
0306       slopes_S9S1_run2[i] = slopes_S9S1_run1[i] * coeffs[i];
0307     }
0308
0309     // S9S1stat configuration
0310     edm::ParameterSetDescription S9S1statDesc;
0311     S9S1statDesc.add<std::vector<double>>("short_optimumSlope", slopes_S9S1_run2);
0312     S9S1statDesc.add<std::vector<double>>(
0313         "shortEnergyParams",
0314         {35.1773, 35.37, 35.7933, 36.4472, 37.3317, 38.4468, 39.7925, 41.3688, 43.1757, 45.2132, 47.4813, 49.98, 52.7093});
0315     S9S1statDesc.add<std::vector<double>>("shortETParams", std::vector<double>(13, 0));
0316     S9S1statDesc.add<std::vector<double>>("long_optimumSlope", slopes_S9S1_run2);
0317     S9S1statDesc.add<std::vector<double>>(
0318         "longEnergyParams", {43.5, 45.7, 48.32, 51.36, 54.82, 58.7, 63.0, 67.72, 72.86, 78.42, 84.4, 90.8, 97.62});
0319     S9S1statDesc.add<std::vector<double>>("longETParams", std::vector<double>(13, 0));
0320     S9S1statDesc.add<int>("HcalAcceptSeverityLevel", 9);
0321     S9S1statDesc.add<bool>("isS8S1", false);
0322     desc.add<edm::ParameterSetDescription>("S9S1stat", S9S1statDesc);
0323   }
0324
0325   {
0326     // S8S1stat configuration
0327     edm::ParameterSetDescription S8S1statDesc;
0328     S8S1statDesc.add<std::vector<double>>(
0329         "short_optimumSlope", {0.30, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10});
0330     S8S1statDesc.add<std::vector<double>>("shortEnergyParams",
0331                                           {40, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100});
0332     S8S1statDesc.add<std::vector<double>>("shortETParams", std::vector<double>(13, 0));
0333     S8S1statDesc.add<std::vector<double>>(
0334         "long_optimumSlope", {0.30, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10, 0.10});
0335     S8S1statDesc.add<std::vector<double>>("longEnergyParams",
0336                                           {40, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100, 100});
0337     S8S1statDesc.add<std::vector<double>>("longETParams", std::vector<double>(13, 0));
0338     S8S1statDesc.add<int>("HcalAcceptSeverityLevel", 9);
0339     S8S1statDesc.add<bool>("isS8S1", true);
0340     desc.add<edm::ParameterSetDescription>("S8S1stat", S8S1statDesc);
0341   }
0342
0343   {
0344     // PETstat configuration
0345     edm::ParameterSetDescription PETstatDesc;
0346     PETstatDesc.add<std::vector<double>>("short_R", {0.8});
0347     PETstatDesc.add<std::vector<double>>(
0348         "shortEnergyParams",
0349         {35.1773, 35.37, 35.7933, 36.4472, 37.3317, 38.4468, 39.7925, 41.3688, 43.1757, 45.2132, 47.4813, 49.98, 52.7093});
0350     PETstatDesc.add<std::vector<double>>("shortETParams", std::vector<double>(13, 0));
0351     PETstatDesc.add<std::vector<double>>("long_R", {0.98});
0352     PETstatDesc.add<std::vector<double>>(
0353         "longEnergyParams", {43.5, 45.7, 48.32, 51.36, 54.82, 58.7, 63.0, 67.72, 72.86, 78.42, 84.4, 90.8, 97.62});
0354     PETstatDesc.add<std::vector<double>>("longETParams", std::vector<double>(13, 0));
0355     PETstatDesc.add<std::vector<double>>("short_R_29", {0.8});
0356     PETstatDesc.add<std::vector<double>>("long_R_29", {0.8});
0357     PETstatDesc.add<int>("HcalAcceptSeverityLevel", 9);
0358     desc.add<edm::ParameterSetDescription>("PETstat", PETstatDesc);
0359   }
0360
0361   descriptions.addWithDefaultLabel(desc);
0362 }
0363
0364 //define this as a plug-in
0365 DEFINE_FWK_MODULE(HFPhase1Reconstructor);