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