Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:51 

0001 // -*- C++ -*-
0002 //
0003 // Package:    RecoLocalCalo/HcalRecProducers
0004 // Class:      HBHEPhase1Reconstructor
0005 //
0006 /**\class HBHEPhase1Reconstructor HBHEPhase1Reconstructor.cc RecoLocalCalo/HcalRecProducers/plugins/HBHEPhase1Reconstructor.cc
0007 
0008  Description: Phase 1 reconstruction module for HB/HE
0009 
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Igor Volobouev
0015 //         Created:  Tue, 21 Jun 2016 00:56:40 GMT
0016 //
0017 //
0018 
0019 // system include files
0020 #include <cmath>
0021 #include <vector>
0022 #include <utility>
0023 #include <algorithm>
0024 
0025 // user include files
0026 #include "FWCore/Framework/interface/ConsumesCollector.h"
0027 #include "FWCore/Framework/interface/Frameworkfwd.h"
0028 #include "FWCore/Framework/interface/stream/EDProducer.h"
0029 #include "FWCore/Utilities/interface/ESGetToken.h"
0030 #include "FWCore/Framework/interface/Event.h"
0031 #include "FWCore/Framework/interface/MakerMacros.h"
0032 #include "FWCore/Utilities/interface/Exception.h"
0033 
0034 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0035 #include "FWCore/Utilities/interface/StreamID.h"
0036 
0037 #include "DataFormats/HcalDetId/interface/HcalGenericDetId.h"
0038 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0039 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0040 #include "DataFormats/HcalRecHit/interface/CaloRecHitAuxSetter.h"
0041 #include "DataFormats/METReco/interface/HcalPhase1FlagLabels.h"
0042 
0043 #include "CalibFormats/HcalObjects/interface/HcalDbService.h"
0044 #include "CalibFormats/HcalObjects/interface/HcalDbRecord.h"
0045 #include "CalibFormats/HcalObjects/interface/HcalCoderDb.h"
0046 #include "CondFormats/HcalObjects/interface/HFPhase1PMTParams.h"
0047 #include "CondFormats/DataRecord/interface/HFPhase1PMTParamsRcd.h"
0048 
0049 #include "CalibFormats/CaloObjects/interface/CaloSamples.h"
0050 
0051 #include "CalibCalorimetry/HcalAlgos/interface/HcalSiPMnonlinearity.h"
0052 
0053 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalChannelProperties.h"
0054 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalChannelPropertiesRecord.h"
0055 #include "RecoLocalCalo/HcalRecAlgos/interface/HBHEStatusBitSetter.h"
0056 #include "RecoLocalCalo/HcalRecAlgos/interface/HBHEPulseShapeFlag.h"
0057 
0058 #include "CondFormats/HcalObjects/interface/HBHENegativeEFilter.h"
0059 #include "CondFormats/DataRecord/interface/HBHENegativeEFilterRcd.h"
0060 
0061 // Parser for Phase 1 HB/HE reco algorithms
0062 #include "RecoLocalCalo/HcalRecAlgos/interface/parseHBHEPhase1AlgoDescription.h"
0063 
0064 // Some helper functions
0065 namespace {
0066   // Class for making SiPM/QIE11 look like HPD/QIE8. HPD/QIE8
0067   // needs only pedestal and gain to convert charge into energy.
0068   // Due to nonlinearities, response of SiPM/QIE11 is substantially
0069   // more complicated. It is possible to calculate all necessary
0070   // quantities from the charge and the info stored in the DB every
0071   // time the raw charge is needed. However, it does not make sense
0072   // to retrieve DB contents stored by channel for every time slice.
0073   // Because of this, we look things up only once, in the constructor.
0074   template <class DFrame>
0075   class RawChargeFromSample {
0076   public:
0077     inline RawChargeFromSample(const int sipmQTSShift,
0078                                const int sipmQNTStoSum,
0079                                const HcalDbService& cond,
0080                                const HcalChannelProperties& properties,
0081                                const CaloSamples& cs,
0082                                const int soi,
0083                                const DFrame& frame,
0084                                const int maxTS) {}
0085 
0086     inline double getRawCharge(const double decodedCharge, const double pedestal) const { return decodedCharge; }
0087   };
0088 
0089   template <>
0090   class RawChargeFromSample<QIE11DataFrame> {
0091   public:
0092     inline RawChargeFromSample(const int sipmQTSShift,
0093                                const int sipmQNTStoSum,
0094                                const HcalDbService& cond,
0095                                const HcalChannelProperties& properties,
0096                                const CaloSamples& cs,
0097                                const int soi,
0098                                const QIE11DataFrame& frame,
0099                                const int maxTS)
0100         : siPMParameter_(*properties.siPMParameter),
0101           fcByPE_(siPMParameter_.getFCByPE()),
0102           corr_(cond.getHcalSiPMCharacteristics()->getNonLinearities(siPMParameter_.getType())) {
0103       if (fcByPE_ <= 0.0)
0104         throw cms::Exception("HBHEPhase1BadDB") << "Invalid fC/PE conversion factor" << std::endl;
0105 
0106       const int firstTS = std::max(soi + sipmQTSShift, 0);
0107       const int lastTS = std::min(firstTS + sipmQNTStoSum, maxTS);
0108       double sipmQ = 0.0;
0109 
0110       for (int ts = firstTS; ts < lastTS; ++ts) {
0111         const float pedestal = properties.pedsAndGains[frame[ts].capid()].pedestal(false);
0112         sipmQ += (cs[ts] - pedestal);
0113       }
0114 
0115       const double effectivePixelsFired = sipmQ / fcByPE_;
0116       factor_ = corr_.getRecoCorrectionFactor(effectivePixelsFired);
0117     }
0118 
0119     inline double getRawCharge(const double decodedCharge, const double pedestal) const {
0120       return (decodedCharge - pedestal) * factor_ + pedestal;
0121 
0122       // Old version of TS-by-TS corrections looked as follows:
0123       // const double sipmQ = decodedCharge - pedestal;
0124       // const double nPixelsFired = sipmQ/fcByPE_;
0125       // return sipmQ*corr_.getRecoCorrectionFactor(nPixelsFired) + pedestal;
0126     }
0127 
0128   private:
0129     const HcalSiPMParameter& siPMParameter_;
0130     double fcByPE_;
0131     HcalSiPMnonlinearity corr_;
0132     double factor_;
0133   };
0134 
0135   float getTDCTimeFromSample(const QIE11DataFrame::Sample& s) { return HcalSpecialTimes::getTDCTime(s.tdc()); }
0136 
0137   float getTDCTimeFromSample(const HcalQIESample&) { return HcalSpecialTimes::UNKNOWN_T_NOTDC; }
0138 
0139   float getDifferentialChargeGain(const HcalQIECoder& coder,
0140                                   const HcalQIEShape& shape,
0141                                   const unsigned adc,
0142                                   const unsigned capid,
0143                                   const bool isQIE11) {
0144     // We have 5-bit ADC mantissa in QIE8 and 6-bit in QIE11
0145     static const unsigned mantissaMaskQIE8 = 0x1f;
0146     static const unsigned mantissaMaskQIE11 = 0x3f;
0147 
0148     const float q = coder.charge(shape, adc, capid);
0149     const unsigned mantissaMask = isQIE11 ? mantissaMaskQIE11 : mantissaMaskQIE8;
0150     const unsigned mantissa = adc & mantissaMask;
0151 
0152     // First, check if we are in the two lowest or two highest ADC
0153     // values for this range. Assume that they have the lowest and
0154     // the highest gain in the range, respectively.
0155     if (mantissa == 0U || mantissa == mantissaMask - 1U)
0156       return coder.charge(shape, adc + 1U, capid) - q;
0157     else if (mantissa == 1U || mantissa == mantissaMask)
0158       return q - coder.charge(shape, adc - 1U, capid);
0159     else {
0160       const float qup = coder.charge(shape, adc + 1U, capid);
0161       const float qdown = coder.charge(shape, adc - 1U, capid);
0162       const float upGain = qup - q;
0163       const float downGain = q - qdown;
0164       const float averageGain = (qup - qdown) / 2.f;
0165       if (std::abs(upGain - downGain) < 0.01f * averageGain)
0166         return averageGain;
0167       else {
0168         // We are in the gain transition region.
0169         // Need to determine if we are in the lower
0170         // gain ADC count or in the higher one.
0171         // This can be done by figuring out if the
0172         // "up" gain is more consistent then the
0173         // "down" gain.
0174         const float q2up = coder.charge(shape, adc + 2U, capid);
0175         const float q2down = coder.charge(shape, adc - 2U, capid);
0176         const float upGain2 = q2up - qup;
0177         const float downGain2 = qdown - q2down;
0178         if (std::abs(upGain2 - upGain) < std::abs(downGain2 - downGain))
0179           return upGain;
0180         else
0181           return downGain;
0182       }
0183     }
0184   }
0185 
0186   // The first element of the pair indicates presence of optical
0187   // link errors. The second indicated presence of capid errors.
0188   std::pair<bool, bool> findHWErrors(const HBHEDataFrame& df, const unsigned len) {
0189     bool linkErr = false;
0190     bool capidErr = false;
0191     if (len) {
0192       int expectedCapid = df[0].capid();
0193       for (unsigned i = 0; i < len; ++i) {
0194         if (df[i].er())
0195           linkErr = true;
0196         if (df[i].capid() != expectedCapid)
0197           capidErr = true;
0198         expectedCapid = (expectedCapid + 1) % 4;
0199       }
0200     }
0201     return std::pair<bool, bool>(linkErr, capidErr);
0202   }
0203 
0204   std::pair<bool, bool> findHWErrors(const QIE11DataFrame& df, const unsigned /* len */) {
0205     return std::pair<bool, bool>(df.linkError(), df.capidError());
0206   }
0207 
0208   std::unique_ptr<HBHEStatusBitSetter> parse_HBHEStatusBitSetter(const edm::ParameterSet& psdigi) {
0209     return std::make_unique<HBHEStatusBitSetter>(
0210         psdigi.getParameter<double>("nominalPedestal"),
0211         psdigi.getParameter<double>("hitEnergyMinimum"),
0212         psdigi.getParameter<int>("hitMultiplicityThreshold"),
0213         psdigi.getParameter<std::vector<edm::ParameterSet> >("pulseShapeParameterSets"));
0214   }
0215 
0216   std::unique_ptr<HBHEPulseShapeFlagSetter> parse_HBHEPulseShapeFlagSetter(const edm::ParameterSet& psPulseShape,
0217                                                                            const bool setLegacyFlags) {
0218     return std::make_unique<HBHEPulseShapeFlagSetter>(
0219         psPulseShape.getParameter<double>("MinimumChargeThreshold"),
0220         psPulseShape.getParameter<double>("TS4TS5ChargeThreshold"),
0221         psPulseShape.getParameter<double>("TS3TS4ChargeThreshold"),
0222         psPulseShape.getParameter<double>("TS3TS4UpperChargeThreshold"),
0223         psPulseShape.getParameter<double>("TS5TS6ChargeThreshold"),
0224         psPulseShape.getParameter<double>("TS5TS6UpperChargeThreshold"),
0225         psPulseShape.getParameter<double>("R45PlusOneRange"),
0226         psPulseShape.getParameter<double>("R45MinusOneRange"),
0227         psPulseShape.getParameter<unsigned int>("TrianglePeakTS"),
0228         psPulseShape.getParameter<std::vector<double> >("LinearThreshold"),
0229         psPulseShape.getParameter<std::vector<double> >("LinearCut"),
0230         psPulseShape.getParameter<std::vector<double> >("RMS8MaxThreshold"),
0231         psPulseShape.getParameter<std::vector<double> >("RMS8MaxCut"),
0232         psPulseShape.getParameter<std::vector<double> >("LeftSlopeThreshold"),
0233         psPulseShape.getParameter<std::vector<double> >("LeftSlopeCut"),
0234         psPulseShape.getParameter<std::vector<double> >("RightSlopeThreshold"),
0235         psPulseShape.getParameter<std::vector<double> >("RightSlopeCut"),
0236         psPulseShape.getParameter<std::vector<double> >("RightSlopeSmallThreshold"),
0237         psPulseShape.getParameter<std::vector<double> >("RightSlopeSmallCut"),
0238         psPulseShape.getParameter<std::vector<double> >("TS4TS5LowerThreshold"),
0239         psPulseShape.getParameter<std::vector<double> >("TS4TS5LowerCut"),
0240         psPulseShape.getParameter<std::vector<double> >("TS4TS5UpperThreshold"),
0241         psPulseShape.getParameter<std::vector<double> >("TS4TS5UpperCut"),
0242         psPulseShape.getParameter<bool>("UseDualFit"),
0243         psPulseShape.getParameter<bool>("TriangleIgnoreSlow"),
0244         setLegacyFlags);
0245   }
0246 
0247   // Determine array index shift so that a partial copy
0248   // of an array maps index "soi" into index "soiWanted"
0249   // under the constraint that there must be no out of
0250   // bounds access.
0251   const int determineIndexShift(const int soi, const int nRead, const int soiWanted, const int nReadWanted) {
0252     assert(nReadWanted <= nRead);
0253     if (soi < 0 || soi >= nRead)
0254       return 0;
0255     else
0256       return std::clamp(soi - soiWanted, 0, nRead - nReadWanted);
0257   }
0258 
0259   // Function for packing TDC data from the QIE11 data frame.
0260   // We want to pack five 6-bit TDC counts so that the soi
0261   // falls on index 3.
0262   uint32_t packTDCData(const QIE11DataFrame& frame, const int soi) {
0263     using namespace CaloRecHitAuxSetter;
0264 
0265     const unsigned six_bits_mask = 0x3f;
0266     const int soiWanted = 3;
0267     const int nRead = frame.size();
0268     const int nTSToCopy = std::min(5, nRead);
0269     const int tsShift = determineIndexShift(soi, nRead, soiWanted, nTSToCopy);
0270 
0271     uint32_t packed = 0;
0272     for (int ts = 0; ts < nTSToCopy; ++ts)
0273       setField(&packed, six_bits_mask, ts * 6, frame[ts + tsShift].tdc());
0274 
0275     // Set bit 30 indicating presence of TDC values
0276     setBit(&packed, 30, true);
0277     return packed;
0278   }
0279 }  // namespace
0280 
0281 //
0282 // class declaration
0283 //
0284 class HBHEPhase1Reconstructor : public edm::stream::EDProducer<> {
0285 public:
0286   explicit HBHEPhase1Reconstructor(const edm::ParameterSet&);
0287   ~HBHEPhase1Reconstructor() override;
0288 
0289   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0290 
0291 private:
0292   void beginRun(edm::Run const&, edm::EventSetup const&) override;
0293   void endRun(edm::Run const&, edm::EventSetup const&) override;
0294   void produce(edm::Event&, const edm::EventSetup&) override;
0295 
0296   // Configuration parameters
0297   std::string algoConfigClass_;
0298   bool processQIE8_;
0299   bool processQIE11_;
0300   bool saveInfos_;
0301   bool saveDroppedInfos_;
0302   bool makeRecHits_;
0303   bool dropZSmarkedPassed_;
0304   bool tsFromDB_;
0305   bool recoParamsFromDB_;
0306   bool saveEffectivePedestal_;
0307   bool use8ts_;
0308   int sipmQTSShift_;
0309   int sipmQNTStoSum_;
0310 
0311   // Parameters for turning status bit setters on/off
0312   bool setNegativeFlagsQIE8_;
0313   bool setNegativeFlagsQIE11_;
0314   bool setNoiseFlagsQIE8_;
0315   bool setNoiseFlagsQIE11_;
0316   bool setPulseShapeFlagsQIE8_;
0317   bool setPulseShapeFlagsQIE11_;
0318 
0319   // Other members
0320   edm::EDGetTokenT<HBHEDigiCollection> tok_qie8_;
0321   edm::EDGetTokenT<QIE11DigiCollection> tok_qie11_;
0322   std::unique_ptr<AbsHBHEPhase1Algo> reco_;
0323   std::unique_ptr<AbsHcalAlgoData> recoConfig_;
0324   edm::ESGetToken<HFPhase1PMTParams, HFPhase1PMTParamsRcd> recoConfigToken_;
0325 
0326   // Status bit setters
0327   const HBHENegativeEFilter* negEFilter_;  // We don't manage this pointer
0328   std::unique_ptr<HBHEStatusBitSetter> hbheFlagSetterQIE8_;
0329   std::unique_ptr<HBHEStatusBitSetter> hbheFlagSetterQIE11_;
0330   std::unique_ptr<HBHEPulseShapeFlagSetter> hbhePulseShapeFlagSetterQIE8_;
0331   std::unique_ptr<HBHEPulseShapeFlagSetter> hbhePulseShapeFlagSetterQIE11_;
0332 
0333   // ES tokens
0334   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> htopoToken_;
0335   edm::ESGetToken<HcalDbService, HcalDbRecord> conditionsToken_;
0336   edm::ESGetToken<HcalChannelPropertiesVec, HcalChannelPropertiesRecord> propertiesToken_;
0337   edm::ESGetToken<HBHENegativeEFilter, HBHENegativeEFilterRcd> negToken_;
0338   edm::ESGetToken<HcalFrontEndMap, HcalFrontEndMapRcd> feMapToken_;
0339 
0340   // For the function below, arguments "infoColl" and/or "rechits"
0341   // are allowed to be null.
0342   template <class DataFrame, class Collection>
0343   void processData(const Collection& coll,
0344                    const HcalTopology& htopo,
0345                    const HcalDbService& cond,
0346                    const HcalChannelPropertiesVec& prop,
0347                    const bool isRealData,
0348                    HBHEChannelInfo* info,
0349                    HBHEChannelInfoCollection* infoColl,
0350                    HBHERecHitCollection* rechits);
0351 
0352   // Methods for setting rechit status bits
0353   void setAsicSpecificBits(const HBHEDataFrame& frame,
0354                            const HcalCoder& coder,
0355                            const HBHEChannelInfo& info,
0356                            const HcalCalibrations& calib,
0357                            int soi,
0358                            HBHERecHit* rh);
0359   void setAsicSpecificBits(const QIE11DataFrame& frame,
0360                            const HcalCoder& coder,
0361                            const HBHEChannelInfo& info,
0362                            const HcalCalibrations& calib,
0363                            int soi,
0364                            HBHERecHit* rh);
0365   void setCommonStatusBits(const HBHEChannelInfo& info, const HcalCalibrations& calib, HBHERecHit* rh);
0366 
0367   void runHBHENegativeEFilter(const HBHEChannelInfo& info, HBHERecHit* rh);
0368 };
0369 
0370 //
0371 // constructors and destructor
0372 //
0373 HBHEPhase1Reconstructor::HBHEPhase1Reconstructor(const edm::ParameterSet& conf)
0374     : algoConfigClass_(conf.getParameter<std::string>("algoConfigClass")),
0375       processQIE8_(conf.getParameter<bool>("processQIE8")),
0376       processQIE11_(conf.getParameter<bool>("processQIE11")),
0377       saveInfos_(conf.getParameter<bool>("saveInfos")),
0378       saveDroppedInfos_(conf.getParameter<bool>("saveDroppedInfos")),
0379       makeRecHits_(conf.getParameter<bool>("makeRecHits")),
0380       dropZSmarkedPassed_(conf.getParameter<bool>("dropZSmarkedPassed")),
0381       tsFromDB_(conf.getParameter<bool>("tsFromDB")),
0382       recoParamsFromDB_(conf.getParameter<bool>("recoParamsFromDB")),
0383       saveEffectivePedestal_(conf.getParameter<bool>("saveEffectivePedestal")),
0384       use8ts_(conf.getParameter<bool>("use8ts")),
0385       sipmQTSShift_(conf.getParameter<int>("sipmQTSShift")),
0386       sipmQNTStoSum_(conf.getParameter<int>("sipmQNTStoSum")),
0387       setNegativeFlagsQIE8_(conf.getParameter<bool>("setNegativeFlagsQIE8")),
0388       setNegativeFlagsQIE11_(conf.getParameter<bool>("setNegativeFlagsQIE11")),
0389       setNoiseFlagsQIE8_(conf.getParameter<bool>("setNoiseFlagsQIE8")),
0390       setNoiseFlagsQIE11_(conf.getParameter<bool>("setNoiseFlagsQIE11")),
0391       setPulseShapeFlagsQIE8_(conf.getParameter<bool>("setPulseShapeFlagsQIE8")),
0392       setPulseShapeFlagsQIE11_(conf.getParameter<bool>("setPulseShapeFlagsQIE11")),
0393       reco_(parseHBHEPhase1AlgoDescription(conf.getParameter<edm::ParameterSet>("algorithm"), consumesCollector())),
0394       negEFilter_(nullptr) {
0395   // Check that the reco algorithm has been successfully configured
0396   if (!reco_.get())
0397     throw cms::Exception("HBHEPhase1BadConfig") << "Invalid HBHEPhase1Algo algorithm configuration" << std::endl;
0398   if (reco_->isConfigurable()) {
0399     if ("HFPhase1PMTParams" != algoConfigClass_) {
0400       throw cms::Exception("HBHEPhase1BadConfig")
0401           << "Invalid HBHEPhase1Reconstructor \"algoConfigClass\" parameter value \"" << algoConfigClass_ << '"'
0402           << std::endl;
0403     }
0404     recoConfigToken_ = esConsumes<edm::Transition::BeginRun>();
0405   }
0406 
0407   // Configure the status bit setters that have been turned on
0408   if (setNoiseFlagsQIE8_)
0409     hbheFlagSetterQIE8_ = parse_HBHEStatusBitSetter(conf.getParameter<edm::ParameterSet>("flagParametersQIE8"));
0410 
0411   if (setNoiseFlagsQIE11_)
0412     hbheFlagSetterQIE11_ = parse_HBHEStatusBitSetter(conf.getParameter<edm::ParameterSet>("flagParametersQIE11"));
0413 
0414   if (setPulseShapeFlagsQIE8_)
0415     hbhePulseShapeFlagSetterQIE8_ =
0416         parse_HBHEPulseShapeFlagSetter(conf.getParameter<edm::ParameterSet>("pulseShapeParametersQIE8"),
0417                                        conf.getParameter<bool>("setLegacyFlagsQIE8"));
0418 
0419   if (setPulseShapeFlagsQIE11_)
0420     hbhePulseShapeFlagSetterQIE11_ =
0421         parse_HBHEPulseShapeFlagSetter(conf.getParameter<edm::ParameterSet>("pulseShapeParametersQIE11"),
0422                                        conf.getParameter<bool>("setLegacyFlagsQIE11"));
0423 
0424   // Consumes and produces statements
0425   if (processQIE8_)
0426     tok_qie8_ = consumes<HBHEDigiCollection>(conf.getParameter<edm::InputTag>("digiLabelQIE8"));
0427 
0428   if (processQIE11_)
0429     tok_qie11_ = consumes<QIE11DigiCollection>(conf.getParameter<edm::InputTag>("digiLabelQIE11"));
0430 
0431   if (saveInfos_)
0432     produces<HBHEChannelInfoCollection>();
0433 
0434   if (makeRecHits_)
0435     produces<HBHERecHitCollection>();
0436 
0437   // ES tokens
0438   htopoToken_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
0439   conditionsToken_ = esConsumes<HcalDbService, HcalDbRecord>();
0440   propertiesToken_ = esConsumes<HcalChannelPropertiesVec, HcalChannelPropertiesRecord>();
0441   if (setNegativeFlagsQIE8_ || setNegativeFlagsQIE11_)
0442     negToken_ = esConsumes<HBHENegativeEFilter, HBHENegativeEFilterRcd>();
0443   if (setNoiseFlagsQIE8_ || setNoiseFlagsQIE11_)
0444     feMapToken_ = esConsumes<HcalFrontEndMap, HcalFrontEndMapRcd, edm::Transition::BeginRun>();
0445 }
0446 
0447 HBHEPhase1Reconstructor::~HBHEPhase1Reconstructor() {
0448   // do anything here that needs to be done at destruction time
0449   // (e.g. close files, deallocate resources etc.)
0450 }
0451 
0452 //
0453 // member functions
0454 //
0455 template <class DFrame, class Collection>
0456 void HBHEPhase1Reconstructor::processData(const Collection& coll,
0457                                           const HcalTopology& htopo,
0458                                           const HcalDbService& cond,
0459                                           const HcalChannelPropertiesVec& prop,
0460                                           const bool isRealData,
0461                                           HBHEChannelInfo* channelInfo,
0462                                           HBHEChannelInfoCollection* infos,
0463                                           HBHERecHitCollection* rechits) {
0464   // If "saveDroppedInfos_" flag is set, fill the info with something
0465   // meaningful even if the database tells us to drop this channel.
0466   // Note that this flag affects only "infos", the rechits are still
0467   // not going to be constructed from such channels.
0468   const bool skipDroppedChannels = !(infos && saveDroppedInfos_);
0469 
0470   // Iterate over the input collection
0471   for (typename Collection::const_iterator it = coll.begin(); it != coll.end(); ++it) {
0472     const DFrame& frame(*it);
0473     const HcalDetId cell(frame.id());
0474 
0475     // Protection against calibration channels which are not
0476     // in the database but can still come in the QIE11DataFrame
0477     // in the laser calibs, etc.
0478     const HcalSubdetector subdet = cell.subdet();
0479     if (!(subdet == HcalSubdetector::HcalBarrel || subdet == HcalSubdetector::HcalEndcap))
0480       continue;
0481 
0482     // Look up the channel properties. This lookup is O(1).
0483     const HcalChannelProperties& properties(prop.at(htopo.detId2denseId(cell)));
0484 
0485     // Check if the database tells us to drop this channel
0486     if (properties.taggedBadByDb && skipDroppedChannels)
0487       continue;
0488 
0489     // Check if the channel is zero suppressed
0490     bool dropByZS = false;
0491     if (dropZSmarkedPassed_)
0492       if (frame.zsMarkAndPass())
0493         dropByZS = true;
0494     if (dropByZS && skipDroppedChannels)
0495       continue;
0496 
0497     // ADC decoding tool
0498     const HcalCoderDb coder(*properties.channelCoder, *properties.shape);
0499 
0500     const bool saveEffectivePeds = channelInfo->hasEffectivePedestals();
0501     const double fcByPE = properties.siPMParameter->getFCByPE();
0502     double darkCurrent = 0.;
0503     double lambda = 0.;
0504     const double noisecorr = properties.siPMParameter->getauxi2();
0505     if (!saveEffectivePeds || saveInfos_) {
0506       // needed for the dark current in the M2 in alternative of the effectivePed
0507       darkCurrent = properties.siPMParameter->getDarkCurrent();
0508       lambda = cond.getHcalSiPMCharacteristics()->getCrossTalk(properties.siPMParameter->getType());
0509     }
0510 
0511     // ADC to fC conversion
0512     CaloSamples cs;
0513     coder.adc2fC(frame, cs);
0514 
0515     // Prepare to iterate over time slices
0516     const int nRead = cs.size();
0517     const int maxTS = std::min(nRead, static_cast<int>(HBHEChannelInfo::MAXSAMPLES));
0518     const int soi = tsFromDB_ ? properties.paramTs->firstSample() : frame.presamples();
0519     const bool badSOI = !(maxTS >= 3 && soi > 0 && soi < maxTS - 1);
0520     if (badSOI) {
0521       edm::LogWarning("HBHEDigi") << " bad SOI/maxTS in cell " << cell
0522                                   << "\n expect maxTS >= 3 && soi > 0 && soi < maxTS - 1"
0523                                   << "\n got maxTS = " << maxTS << ", SOI = " << soi;
0524     }
0525 
0526     const RawChargeFromSample<DFrame> rcfs(sipmQTSShift_, sipmQNTStoSum_, cond, properties, cs, soi, frame, maxTS);
0527     int soiCapid = 4;
0528 
0529     // Typical expected cases:
0530     //   maxTS = 10, SOI = 4 (typical 10-TS situation, in data or MC)
0531     //   maxTS = 10, SOI = 5 (new, for better modeling of SiPM noise in MC)
0532     //   maxTS = 8,  SOI = 3 (typical 8-TS situation in data)
0533     //
0534     // We want to fill the HBHEChannelInfo object with either
0535     // 8 or 10 time slices, depending on the maxTS value and
0536     // on the "use8ts_" parameter setting. If we want 8 time
0537     // slices in the HBHEChannelInfo, we want the SOI to fall
0538     // on the time slice number 3. For now, this number is not
0539     // expected to change and will be hardwired.
0540     //
0541     int nTSToCopy = maxTS;
0542     int tsShift = 0;
0543     if (maxTS > 8 && use8ts_) {
0544       const int soiWanted = 3;
0545 
0546       // We need to chop off the excess time slices
0547       // and configure the TS shift for filling
0548       // HBHEChannelInfo from the data frame.
0549       nTSToCopy = 8;
0550       tsShift = determineIndexShift(soi, nRead, soiWanted, nTSToCopy);
0551     }
0552 
0553     // Go over time slices and fill the samples
0554     for (int copyTS = 0; copyTS < nTSToCopy; ++copyTS) {
0555       const int inputTS = copyTS + tsShift;
0556       auto s(frame[inputTS]);
0557       const uint8_t adc = s.adc();
0558       const int capid = s.capid();
0559       const HcalPipelinePedestalAndGain& pAndGain(properties.pedsAndGains.at(capid));
0560 
0561       // Always use QIE-only pedestal for this computation
0562       const double rawCharge = rcfs.getRawCharge(cs[inputTS], pAndGain.pedestal(false));
0563       const float t = getTDCTimeFromSample(s);
0564       const float dfc = getDifferentialChargeGain(
0565           *properties.channelCoder, *properties.shape, adc, capid, channelInfo->hasTimeInfo());
0566       channelInfo->setSample(copyTS,
0567                              adc,
0568                              dfc,
0569                              rawCharge,
0570                              pAndGain.pedestal(saveEffectivePeds),
0571                              pAndGain.pedestalWidth(saveEffectivePeds),
0572                              pAndGain.gain(),
0573                              pAndGain.gainWidth(),
0574                              t);
0575       if (inputTS == soi)
0576         soiCapid = capid;
0577     }
0578 
0579     // Fill the overall channel info items
0580     const int fitSoi = soi - tsShift;
0581     const int pulseShapeID = properties.paramTs->pulseShapeID();
0582     const std::pair<bool, bool> hwerr = findHWErrors(frame, maxTS);
0583     channelInfo->setChannelInfo(cell,
0584                                 pulseShapeID,
0585                                 nTSToCopy,
0586                                 fitSoi,
0587                                 soiCapid,
0588                                 darkCurrent,
0589                                 fcByPE,
0590                                 lambda,
0591                                 noisecorr,
0592                                 hwerr.first,
0593                                 hwerr.second,
0594                                 properties.taggedBadByDb || dropByZS || badSOI);
0595 
0596     // If needed, add the channel info to the output collection
0597     const bool makeThisRechit = !channelInfo->isDropped();
0598     if (infos && (saveDroppedInfos_ || makeThisRechit))
0599       infos->push_back(*channelInfo);
0600 
0601     // Reconstruct the rechit
0602     if (rechits && makeThisRechit) {
0603       const HcalRecoParam* pptr = nullptr;
0604       if (recoParamsFromDB_)
0605         pptr = properties.paramTs;
0606       HBHERecHit rh = reco_->reconstruct(*channelInfo, pptr, *properties.calib, isRealData);
0607       if (rh.id().rawId()) {
0608         setAsicSpecificBits(frame, coder, *channelInfo, *properties.calib, soi, &rh);
0609         setCommonStatusBits(*channelInfo, *properties.calib, &rh);
0610         rechits->push_back(rh);
0611       }
0612     }
0613   }
0614 }
0615 
0616 void HBHEPhase1Reconstructor::setCommonStatusBits(const HBHEChannelInfo& /* info */,
0617                                                   const HcalCalibrations& /* calib */,
0618                                                   HBHERecHit* /* rh */) {}
0619 
0620 void HBHEPhase1Reconstructor::setAsicSpecificBits(const HBHEDataFrame& frame,
0621                                                   const HcalCoder& coder,
0622                                                   const HBHEChannelInfo& info,
0623                                                   const HcalCalibrations& calib,
0624                                                   int /* soi */,
0625                                                   HBHERecHit* rh) {
0626   if (setNoiseFlagsQIE8_)
0627     hbheFlagSetterQIE8_->rememberHit(*rh);
0628 
0629   if (setPulseShapeFlagsQIE8_)
0630     hbhePulseShapeFlagSetterQIE8_->SetPulseShapeFlags(*rh, frame, coder, calib);
0631 
0632   if (setNegativeFlagsQIE8_)
0633     runHBHENegativeEFilter(info, rh);
0634 
0635   rh->setAuxTDC(0U);
0636 }
0637 
0638 void HBHEPhase1Reconstructor::setAsicSpecificBits(const QIE11DataFrame& frame,
0639                                                   const HcalCoder& coder,
0640                                                   const HBHEChannelInfo& info,
0641                                                   const HcalCalibrations& calib,
0642                                                   const int soi,
0643                                                   HBHERecHit* rh) {
0644   if (setNoiseFlagsQIE11_)
0645     hbheFlagSetterQIE11_->rememberHit(*rh);
0646 
0647   if (setPulseShapeFlagsQIE11_)
0648     hbhePulseShapeFlagSetterQIE11_->SetPulseShapeFlags(*rh, frame, coder, calib);
0649 
0650   if (setNegativeFlagsQIE11_)
0651     runHBHENegativeEFilter(info, rh);
0652 
0653   rh->setAuxTDC(packTDCData(frame, soi));
0654 }
0655 
0656 void HBHEPhase1Reconstructor::runHBHENegativeEFilter(const HBHEChannelInfo& info, HBHERecHit* rh) {
0657   double ts[HBHEChannelInfo::MAXSAMPLES];
0658   const unsigned nRead = info.nSamples();
0659   for (unsigned i = 0; i < nRead; ++i)
0660     ts[i] = info.tsCharge(i);
0661   const bool passes = negEFilter_->checkPassFilter(info.id(), &ts[0], nRead);
0662   if (!passes)
0663     rh->setFlagField(1, HcalPhase1FlagLabels::HBHENegativeNoise);
0664 }
0665 
0666 // ------------ method called to produce the data  ------------
0667 void HBHEPhase1Reconstructor::produce(edm::Event& e, const edm::EventSetup& eventSetup) {
0668   using namespace edm;
0669 
0670   // Get the Hcal topology
0671   const HcalTopology* htopo = &eventSetup.getData(htopoToken_);
0672 
0673   // Fetch the calibrations
0674   const HcalDbService* conditions = &eventSetup.getData(conditionsToken_);
0675   const HcalChannelPropertiesVec* prop = &eventSetup.getData(propertiesToken_);
0676 
0677   // Configure the negative energy filter
0678   if (setNegativeFlagsQIE8_ || setNegativeFlagsQIE11_) {
0679     negEFilter_ = &eventSetup.getData(negToken_);
0680   }
0681 
0682   // Find the input data
0683   unsigned maxOutputSize = 0;
0684   Handle<HBHEDigiCollection> hbDigis;
0685   if (processQIE8_) {
0686     e.getByToken(tok_qie8_, hbDigis);
0687     maxOutputSize += hbDigis->size();
0688   }
0689 
0690   Handle<QIE11DigiCollection> heDigis;
0691   if (processQIE11_) {
0692     e.getByToken(tok_qie11_, heDigis);
0693     maxOutputSize += heDigis->size();
0694   }
0695 
0696   // Create new output collections
0697   std::unique_ptr<HBHEChannelInfoCollection> infos;
0698   if (saveInfos_) {
0699     infos = std::make_unique<HBHEChannelInfoCollection>();
0700     infos->reserve(maxOutputSize);
0701   }
0702 
0703   std::unique_ptr<HBHERecHitCollection> out;
0704   if (makeRecHits_) {
0705     out = std::make_unique<HBHERecHitCollection>();
0706     out->reserve(maxOutputSize);
0707   }
0708 
0709   // Process the input collections, filling the output ones
0710   const bool isData = e.isRealData();
0711   if (processQIE8_) {
0712     if (setNoiseFlagsQIE8_)
0713       hbheFlagSetterQIE8_->Clear();
0714 
0715     HBHEChannelInfo channelInfo(false, false);
0716     processData<HBHEDataFrame>(*hbDigis, *htopo, *conditions, *prop, isData, &channelInfo, infos.get(), out.get());
0717     if (setNoiseFlagsQIE8_)
0718       hbheFlagSetterQIE8_->SetFlagsFromRecHits(*out);
0719   }
0720 
0721   if (processQIE11_) {
0722     if (setNoiseFlagsQIE11_)
0723       hbheFlagSetterQIE11_->Clear();
0724 
0725     HBHEChannelInfo channelInfo(true, saveEffectivePedestal_);
0726     processData<QIE11DataFrame>(*heDigis, *htopo, *conditions, *prop, isData, &channelInfo, infos.get(), out.get());
0727     if (setNoiseFlagsQIE11_)
0728       hbheFlagSetterQIE11_->SetFlagsFromRecHits(*out);
0729   }
0730 
0731   // Add the output collections to the event record
0732   if (saveInfos_)
0733     e.put(std::move(infos));
0734   if (makeRecHits_)
0735     e.put(std::move(out));
0736 }
0737 
0738 // ------------ method called when starting to processes a run  ------------
0739 void HBHEPhase1Reconstructor::beginRun(edm::Run const& r, edm::EventSetup const& es) {
0740   if (reco_->isConfigurable()) {
0741     recoConfig_ = std::make_unique<HFPhase1PMTParams>(es.getData(recoConfigToken_));
0742     if (!reco_->configure(recoConfig_.get()))
0743       throw cms::Exception("HBHEPhase1BadConfig")
0744           << "Failed to configure HBHEPhase1Algo algorithm from EventSetup" << std::endl;
0745   }
0746 
0747   if (setNoiseFlagsQIE8_ || setNoiseFlagsQIE11_) {
0748     if (auto handle = es.getHandle(feMapToken_)) {
0749       const HcalFrontEndMap& hfemap = *handle;
0750       if (setNoiseFlagsQIE8_)
0751         hbheFlagSetterQIE8_->SetFrontEndMap(&hfemap);
0752       if (setNoiseFlagsQIE11_)
0753         hbheFlagSetterQIE11_->SetFrontEndMap(&hfemap);
0754     } else {
0755       edm::LogWarning("EventSetup") << "HBHEPhase1Reconstructor failed to get HcalFrontEndMap!" << std::endl;
0756     }
0757   }
0758 
0759   reco_->beginRun(r, es);
0760 }
0761 
0762 void HBHEPhase1Reconstructor::endRun(edm::Run const&, edm::EventSetup const&) { reco_->endRun(); }
0763 
0764 #define add_param_set(name) /**/     \
0765   edm::ParameterSetDescription name; \
0766   name.setAllowAnything();           \
0767   desc.add<edm::ParameterSetDescription>(#name, name)
0768 
0769 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0770 void HBHEPhase1Reconstructor::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0771   edm::ParameterSetDescription desc;
0772 
0773   desc.add<edm::InputTag>("digiLabelQIE8");
0774   desc.add<edm::InputTag>("digiLabelQIE11");
0775   desc.add<std::string>("algoConfigClass");
0776   desc.add<bool>("processQIE8");
0777   desc.add<bool>("processQIE11");
0778   desc.add<bool>("saveInfos");
0779   desc.add<bool>("saveDroppedInfos");
0780   desc.add<bool>("makeRecHits");
0781   desc.add<bool>("dropZSmarkedPassed");
0782   desc.add<bool>("tsFromDB");
0783   desc.add<bool>("recoParamsFromDB");
0784   desc.add<bool>("saveEffectivePedestal", false);
0785   desc.add<bool>("use8ts", false);
0786   desc.add<int>("sipmQTSShift", 0);
0787   desc.add<int>("sipmQNTStoSum", 3);
0788   desc.add<bool>("setNegativeFlagsQIE8");
0789   desc.add<bool>("setNegativeFlagsQIE11");
0790   desc.add<bool>("setNoiseFlagsQIE8");
0791   desc.add<bool>("setNoiseFlagsQIE11");
0792   desc.add<bool>("setPulseShapeFlagsQIE8");
0793   desc.add<bool>("setPulseShapeFlagsQIE11");
0794   desc.add<bool>("setLegacyFlagsQIE8");
0795   desc.add<bool>("setLegacyFlagsQIE11");
0796 
0797   desc.add<edm::ParameterSetDescription>("algorithm", fillDescriptionForParseHBHEPhase1Algo());
0798   add_param_set(flagParametersQIE8);
0799   add_param_set(flagParametersQIE11);
0800   add_param_set(pulseShapeParametersQIE8);
0801   add_param_set(pulseShapeParametersQIE11);
0802 
0803   descriptions.addDefault(desc);
0804 }
0805 
0806 //define this as a plug-in
0807 DEFINE_FWK_MODULE(HBHEPhase1Reconstructor);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!