HcalRecProducers/src/HFPreReconstructor.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoLocalCalo/HcalRecProducers
0004 // Class:      HFPreReconstructor
0005 //
0006 /**\class HFPreReconstructor HFPreReconstructor.cc RecoLocalCalo/HcalRecProducers/src/HFPreReconstructor.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, 05 May 2016 00:17:51 GMT
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021 #include <utility>
0022 #include <cassert>
0023
0024 // user include files
0025 #include "FWCore/Framework/interface/Frameworkfwd.h"
0026 #include "FWCore/Framework/interface/stream/EDProducer.h"
0027 #include "FWCore/Utilities/interface/ESGetToken.h"
0028 #include "FWCore/Framework/interface/Event.h"
0029 #include "FWCore/Framework/interface/MakerMacros.h"
0030
0031 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0032 #include "FWCore/Utilities/interface/StreamID.h"
0033
0034 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0035 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0036
0037 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
0038
0039 #include "CondFormats/HcalObjects/interface/HcalRecoParam.h"
0040
0041 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0042 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0043
0044 #include "RecoLocalCalo/HcalRecAlgos/interface/HFPreRecAlgo.h"
0045 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalChannelProperties.h"
0046 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalChannelPropertiesRecord.h"
0047
0048 //
0049 // class declaration
0050 //
0051
0052 class HFPreReconstructor : public edm::stream::EDProducer<> {
0053 public:
0054   explicit HFPreReconstructor(const edm::ParameterSet&);
0055   ~HFPreReconstructor() override;
0056
0057   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0058
0059 private:
0060   typedef std::pair<HcalDetId, int> PmtAnodeId;
0061   typedef std::pair<PmtAnodeId, const HFQIE10Info*> QIE10InfoWithId;
0062
0063   void beginRun(const edm::Run&, const edm::EventSetup&) override;
0064   void produce(edm::Event&, const edm::EventSetup&) override;
0065
0066   // Module configuration parameters
0067   edm::InputTag inputLabel_;
0068   int forceSOI_;
0069   int soiShift_;
0070   bool dropZSmarkedPassed_;
0071   bool tsFromDB_;
0072
0073   // Other members
0074   HFPreRecAlgo reco_;
0075   edm::EDGetTokenT<QIE10DigiCollection> tok_hfQIE10_;
0076   std::vector<HFQIE10Info> qie10Infos_;
0077   std::vector<QIE10InfoWithId> sortedQIE10Infos_;
0078
0079   // Fill qie10Infos_ from the event data
0080   void fillInfos(const edm::Event& e, const edm::EventSetup& eventSetup);
0081
0082   // Fill out sortedQIE10Infos_ from qie10Infos_ and return the PMT count
0083   unsigned sortDataByPmt();
0084
0085   // ES tokens
0086   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> htopoToken_;
0087   edm::ESGetToken<HcalChannelPropertiesVec, HcalChannelPropertiesRecord> propertiesToken_;
0088 };
0089
0090 //
0091 // constructors and destructor
0092 //
0093 HFPreReconstructor::HFPreReconstructor(const edm::ParameterSet& conf)
0094     : inputLabel_(conf.getParameter<edm::InputTag>("digiLabel")),
0095       forceSOI_(conf.getParameter<int>("forceSOI")),
0096       soiShift_(conf.getParameter<int>("soiShift")),
0097       dropZSmarkedPassed_(conf.getParameter<bool>("dropZSmarkedPassed")),
0098       tsFromDB_(conf.getParameter<bool>("tsFromDB")),
0099       reco_(conf.getParameter<bool>("sumAllTimeSlices")) {
0100   // Describe consumed data
0101   tok_hfQIE10_ = consumes<QIE10DigiCollection>(inputLabel_);
0102
0103   // Register the product
0104   produces<HFPreRecHitCollection>();
0105
0106   // ES tokens
0107   htopoToken_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
0108   propertiesToken_ = esConsumes<HcalChannelPropertiesVec, HcalChannelPropertiesRecord>();
0109 }
0110
0111 HFPreReconstructor::~HFPreReconstructor() {
0112   // do anything here that needs to be done at destruction time
0113   // (e.g. close files, deallocate resources etc.)
0114 }
0115
0116 //
0117 // member functions
0118 //
0119 unsigned HFPreReconstructor::sortDataByPmt() {
0120   sortedQIE10Infos_.clear();
0121   unsigned pmtCount = 0;
0122   const unsigned sz = qie10Infos_.size();
0123   if (sz) {
0124     // Perform sorting
0125     sortedQIE10Infos_.reserve(sz);
0126     const HFQIE10Info* info = &qie10Infos_[0];
0127     for (unsigned i = 0; i < sz; ++i) {
0128       const HcalDetId id(info[i].id());
0129       sortedQIE10Infos_.push_back(QIE10InfoWithId(PmtAnodeId(id.baseDetId(), id.depth()), info + i));
0130     }
0131     std::sort(sortedQIE10Infos_.begin(), sortedQIE10Infos_.end());
0132
0133     // Count the PMTs
0134     HcalDetId previousBaseId(sortedQIE10Infos_[0].first.first);
0135     pmtCount = 1;
0136     for (unsigned i = 1; i < sz; ++i) {
0137       const HcalDetId baseId(sortedQIE10Infos_[i].first.first);
0138       if (baseId != previousBaseId) {
0139         previousBaseId = baseId;
0140         ++pmtCount;
0141       }
0142     }
0143   }
0144   return pmtCount;
0145 }
0146
0147 void HFPreReconstructor::fillInfos(const edm::Event& e, const edm::EventSetup& eventSetup) {
0148   using namespace edm;
0149
0150   // Clear the collection we want to fill in this method
0151   qie10Infos_.clear();
0152
0153   // Get the calibrations
0154   const HcalTopology& htopo = eventSetup.getData(htopoToken_);
0155   const HcalChannelPropertiesVec& prop = eventSetup.getData(propertiesToken_);
0156
0157   // Get the input collection
0158   Handle<QIE10DigiCollection> digi;
0159   e.getByToken(tok_hfQIE10_, digi);
0160
0161   const unsigned inputSize = digi->size();
0162   if (inputSize) {
0163     // Process the digis and fill out the HFQIE10Info vector
0164     qie10Infos_.reserve(inputSize);
0165
0166     for (QIE10DigiCollection::const_iterator it = digi->begin(); it != digi->end(); ++it) {
0167       const QIE10DataFrame& frame(*it);
0168       const HcalDetId cell(frame.id());
0169
0170       // Protection against calibration channels which are not
0171       // in the database but can still come in the QIE10DataFrame
0172       // in the laser calibs, etc.
0173       if (cell.subdet() != HcalSubdetector::HcalForward)
0174         continue;
0175
0176       // Check zero suppression
0177       if (dropZSmarkedPassed_)
0178         if (frame.zsMarkAndPass())
0179           continue;
0180
0181       // Look up the channel properties. This lookup is O(1).
0182       const HcalChannelProperties& properties(prop.at(htopo.detId2denseId(cell)));
0183
0184       // ADC decoding tool
0185       const HcalCoderDb coder(*properties.channelCoder, *properties.shape);
0186
0187       int tsToUse = forceSOI_;
0188       if (tsToUse < 0) {
0189         if (tsFromDB_) {
0190           tsToUse = properties.paramTs->firstSample();
0191         } else {
0192           // Get the "sample of interest" from the data frame itself
0193           tsToUse = frame.presamples();
0194         }
0195       }
0196
0197       // Reconstruct the charge, energy, etc
0198       const HFQIE10Info& info = reco_.reconstruct(frame, tsToUse + soiShift_, coder, properties);
0199       if (info.id().rawId())
0200         qie10Infos_.push_back(info);
0201     }
0202   }
0203 }
0204
0205 void HFPreReconstructor::beginRun(const edm::Run& r, const edm::EventSetup& es) {}
0206
0207 // ------------ method called to produce the data  ------------
0208 void HFPreReconstructor::produce(edm::Event& e, const edm::EventSetup& eventSetup) {
0209   // Process the input data
0210   fillInfos(e, eventSetup);
0211
0212   // Create a new output collection
0213   std::unique_ptr<HFPreRecHitCollection> out(std::make_unique<HFPreRecHitCollection>());
0214
0215   // Fill the output collection
0216   const unsigned pmtCount = sortDataByPmt();
0217   if (pmtCount) {
0218     out->reserve(pmtCount);
0219     const unsigned sz = sortedQIE10Infos_.size();
0220     HcalDetId previousBaseId(sortedQIE10Infos_[0].first.first);
0221     unsigned nFound = 1;
0222     for (unsigned i = 1; i <= sz; ++i) {
0223       bool appendData = i == sz;
0224       if (i < sz) {
0225         const HcalDetId baseId(sortedQIE10Infos_[i].first.first);
0226         if (baseId == previousBaseId)
0227           ++nFound;
0228         else {
0229           appendData = true;
0230           previousBaseId = baseId;
0231         }
0232       }
0233
0234       if (appendData) {
0235         // If we have found more than two QIE10 with the same base id,
0236         // there is a bug somewhere in the dataframe. We can't do
0237         // anything useful about it here. Once we make sure that
0238         // everything works as expected, this assertion can be removed.
0239         assert(nFound <= 2);
0240
0241         const HFQIE10Info* first = nullptr;
0242         const HFQIE10Info* second = sortedQIE10Infos_[i - 1].second;
0243
0244         if (nFound >= 2)
0245           first = sortedQIE10Infos_[i - 2].second;
0246         else if (sortedQIE10Infos_[i - 1].first.second < 3) {
0247           // Only one QIE10 readout found for this PMT.
0248           // Arrange for depth 1 and 2 to be "first".
0249           first = second;
0250           second = nullptr;
0251         }
0252
0253         out->push_back(HFPreRecHit(sortedQIE10Infos_[i - nFound].first.first, first, second));
0254
0255         // Reset the QIE find count for this base id
0256         nFound = 1;
0257       }
0258     }
0259
0260     assert(out->size() == pmtCount);
0261   }
0262
0263   // Add the output collection to the event record
0264   e.put(std::move(out));
0265 }
0266
0267 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0268 void HFPreReconstructor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0269   edm::ParameterSetDescription desc;
0270
0271   desc.add<edm::InputTag>("digiLabel");
0272   desc.add<int>("forceSOI", -1);
0273   desc.add<int>("soiShift", 0);
0274   desc.add<bool>("dropZSmarkedPassed");
0275   desc.add<bool>("tsFromDB");
0276   desc.add<bool>("sumAllTimeSlices");
0277
0278   descriptions.addDefault(desc);
0279 }
0280
0281 //define this as a plug-in
0282 DEFINE_FWK_MODULE(HFPreReconstructor);