HcalRecAlgos/src/SimpleHBHEPhase1Algo.cc

0001 #include <algorithm>
0002
0003 #include "CalibCalorimetry/HcalAlgos/interface/HcalTimeSlew.h"
0004
0005 #include "RecoLocalCalo/HcalRecAlgos/interface/SimpleHBHEPhase1Algo.h"
0006 #include "RecoLocalCalo/HcalRecAlgos/interface/HcalCorrectionFunctions.h"
0007
0008 #include "FWCore/Framework/interface/ConsumesCollector.h"
0009 #include "FWCore/Framework/interface/Run.h"
0010
0011 #include "DataFormats/HcalRecHit/interface/HBHERecHitAuxSetter.h"
0012 #include "DataFormats/METReco/interface/HcalPhase1FlagLabels.h"
0013 #include "CondFormats/DataRecord/interface/HcalTimeSlewRecord.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015
0016 namespace {
0017   // Maximum fractional error for calculating Method 0
0018   // pulse containment correction
0019   constexpr float PulseContainmentFractionalError = 0.002f;
0020 }  // namespace
0021
0022 SimpleHBHEPhase1Algo::SimpleHBHEPhase1Algo(const int firstSampleShift,
0023                                            const int samplesToAdd,
0024                                            const float phaseNS,
0025                                            const float timeShift,
0026                                            const bool correctForPhaseContainment,
0027                                            const bool applyLegacyHBMCorrection,
0028                                            const bool applyFixPCC,
0029                                            std::unique_ptr<PulseShapeFitOOTPileupCorrection> m2,
0030                                            std::unique_ptr<HcalDeterministicFit> detFit,
0031                                            std::unique_ptr<MahiFit> mahi,
0032                                            edm::ConsumesCollector iC)
0033     : delayToken_(iC.esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", "HBHE"))),
0034       pulseCorr_(PulseContainmentFractionalError, applyFixPCC, iC),
0035       firstSampleShift_(firstSampleShift),
0036       samplesToAdd_(samplesToAdd),
0037       phaseNS_(phaseNS),
0038       timeShift_(timeShift),
0039       runnum_(0),
0040       corrFPC_(correctForPhaseContainment),
0041       applyLegacyHBMCorrection_(applyLegacyHBMCorrection),
0042       psFitOOTpuCorr_(std::move(m2)),
0043       hltOOTpuCorr_(std::move(detFit)),
0044       mahiOOTpuCorr_(std::move(mahi)) {
0045   hcalTimeSlew_delay_ = nullptr;
0046 }
0047
0048 void SimpleHBHEPhase1Algo::beginRun(const edm::Run& r, const edm::EventSetup& es) {
0049   hcalTimeSlew_delay_ = &es.getData(delayToken_);
0050
0051   runnum_ = r.run();
0052   pulseCorr_.beginRun(es);
0053 }
0054
0055 void SimpleHBHEPhase1Algo::endRun() { runnum_ = 0; }
0056
0057 HBHERecHit SimpleHBHEPhase1Algo::reconstruct(const HBHEChannelInfo& info,
0058                                              const HcalRecoParam* params,
0059                                              const HcalCalibrations& calibs,
0060                                              const bool isData) {
0061   const HcalDetId channelId(info.id());
0062
0063   // Calculate "Method 0" quantities
0064   float m0t = 0.f, m0E = 0.f;
0065   {
0066     int ibeg = static_cast<int>(info.soi()) + firstSampleShift_;
0067     if (ibeg < 0)
0068       ibeg = 0;
0069     const int nSamplesToAdd = params ? params->samplesToAdd() : samplesToAdd_;
0070     const double fc_ampl = info.chargeInWindow(ibeg, ibeg + nSamplesToAdd);
0071     const bool applyContainment = params ? params->correctForPhaseContainment() : corrFPC_;
0072     const float phasens = params ? params->correctionPhaseNS() : phaseNS_;
0073     m0E = m0Energy(info, fc_ampl, applyContainment, phasens, nSamplesToAdd);
0074     m0E *= hbminusCorrectionFactor(channelId, m0E, isData);
0075     m0t = m0Time(info, fc_ampl, nSamplesToAdd);
0076   }
0077
0078   // Run "Method 2"
0079   float m2t = 0.f, m2E = 0.f, chi2 = -1.f;
0080   bool useTriple = false;
0081   const PulseShapeFitOOTPileupCorrection* method2 = psFitOOTpuCorr_.get();
0082   if (method2) {
0083     psFitOOTpuCorr_->setPulseShapeTemplate(
0084         theHcalPulseShapes_.getShape(info.recoShape()), !info.hasTimeInfo(), info.nSamples(), hcalTimeSlew_delay_);
0085     // "phase1Apply" call below sets m2E, m2t, useTriple, and chi2.
0086     // These parameters are pased by non-const reference.
0087     method2->phase1Apply(info, m2E, m2t, useTriple, chi2);
0088     m2E *= hbminusCorrectionFactor(channelId, m2E, isData);
0089   }
0090
0091   // Run "Method 3"
0092   float m3t = 0.f, m3E = 0.f;
0093   const HcalDeterministicFit* method3 = hltOOTpuCorr_.get();
0094   if (method3) {
0095     // "phase1Apply" sets m3E and m3t (pased by non-const reference)
0096     method3->phase1Apply(info, m3E, m3t, hcalTimeSlew_delay_);
0097     m3E *= hbminusCorrectionFactor(channelId, m3E, isData);
0098   }
0099
0100   // Run Mahi
0101   float m4E = 0.f, m4ESOIPlusOne = 0.f, m4chi2 = -1.f;
0102   float m4T = 0.f;
0103   bool m4UseTriple = false;
0104
0105   const MahiFit* mahi = mahiOOTpuCorr_.get();
0106
0107   if (mahi) {
0108     mahiOOTpuCorr_->setPulseShapeTemplate(
0109         info.recoShape(), theHcalPulseShapes_, info.hasTimeInfo(), hcalTimeSlew_delay_, info.nSamples(), info.tsGain(0));
0110     mahi->phase1Apply(info, m4E, m4ESOIPlusOne, m4T, m4UseTriple, m4chi2);
0111     m4E *= hbminusCorrectionFactor(channelId, m4E, isData);
0112     m4ESOIPlusOne *= hbminusCorrectionFactor(channelId, m4ESOIPlusOne, isData);
0113   }
0114
0115   // Finally, construct the rechit
0116   HBHERecHit rh;
0117
0118   float rhE = m0E;
0119   float rht = m0t;
0120   float rhX = -1.f;
0121   if (mahi) {
0122     rhE = m4E;
0123     rht = m4T;
0124     rhX = m4chi2;
0125   } else if (method2) {
0126     rhE = m2E;
0127     rht = m2t;
0128     rhX = chi2;
0129   } else if (method3) {
0130     rhE = m3E;
0131     rht = m3t;
0132   }
0133   float tdcTime = info.soiRiseTime();
0134   if (!HcalSpecialTimes::isSpecial(tdcTime))
0135     tdcTime += timeShift_;
0136   rh = HBHERecHit(channelId, rhE, rht, tdcTime);
0137   rh.setRawEnergy(m0E);
0138   if (method3) {
0139     rh.setAuxEnergy(m3E);
0140   } else {
0141     rh.setAuxEnergy(m4ESOIPlusOne);
0142   }
0143   rh.setChiSquared(rhX);
0144
0145   // Set rechit aux words
0146   HBHERecHitAuxSetter::setAux(info, &rh);
0147
0148   // Set some rechit flags (here, for Method 2/Mahi)
0149   if (useTriple || m4UseTriple)
0150     rh.setFlagField(1, HcalPhase1FlagLabels::HBHEPulseFitBit);
0151
0152   return rh;
0153 }
0154
0155 float SimpleHBHEPhase1Algo::hbminusCorrectionFactor(const HcalDetId& cell,
0156                                                     const float energy,
0157                                                     const bool isRealData) const {
0158   float corr = 1.f;
0159   if (applyLegacyHBMCorrection_ && isRealData && runnum_ > 0)
0160     if (cell.subdet() == HcalBarrel) {
0161       const int ieta = cell.ieta();
0162       const int iphi = cell.iphi();
0163       corr = hbminus_special_ecorr(ieta, iphi, energy, runnum_);
0164     }
0165   return corr;
0166 }
0167
0168 float SimpleHBHEPhase1Algo::m0Energy(const HBHEChannelInfo& info,
0169                                      const double fc_ampl,
0170                                      const bool applyContainmentCorrection,
0171                                      const double phaseNs,
0172                                      const int nSamplesToAdd) {
0173   int ibeg = static_cast<int>(info.soi()) + firstSampleShift_;
0174   if (ibeg < 0)
0175     ibeg = 0;
0176   double e = info.energyInWindow(ibeg, ibeg + nSamplesToAdd);
0177
0178   // Pulse containment correction
0179   {
0180     double corrFactor = 1.0;
0181     if (applyContainmentCorrection)
0182       corrFactor = pulseCorr_.get(info.id(), nSamplesToAdd, phaseNs)->getCorrection(fc_ampl);
0183     e *= corrFactor;
0184   }
0185
0186   return e;
0187 }
0188
0189 float SimpleHBHEPhase1Algo::m0Time(const HBHEChannelInfo& info,
0190                                    const double fc_ampl,
0191                                    const int nSamplesToExamine) const {
0192   float time = -9999.f;  // historic value
0193
0194   const unsigned nSamples = info.nSamples();
0195   if (nSamples > 2U) {
0196     const int soi = info.soi();
0197     int ibeg = soi + firstSampleShift_;
0198     if (ibeg < 0)
0199       ibeg = 0;
0200     const int iend = std::min(ibeg + nSamplesToExamine, (int)nSamples - 1);  // actual array
0201
0202     unsigned maxI = info.peakEnergyTS((unsigned)ibeg, (unsigned)iend);  // requires unsigned params
0203     if (maxI < HBHEChannelInfo::MAXSAMPLES) {
0204       if (maxI >= nSamples)
0205         maxI = nSamples - 1U;  // just in case
0206
0207       // Simplified evaluation for Phase1
0208       float emax0 = info.tsEnergy(maxI);
0209       float emax1 = 0.f;
0210       if (maxI < (nSamples - 1U))
0211         emax1 = info.tsEnergy(maxI + 1U);
0212
0213       // consider soi reference for collisions
0214       int position = (int)maxI;
0215       if (nSamplesToExamine < (int)nSamples)
0216         position -= soi;
0217
0218       time = 25.f * (float)position;
0219       if (emax0 > 0.f && emax1 > 0.f)
0220         time += 25.f * emax1 / (emax0 + emax1);  // 1st order corr.
0221
0222       // TimeSlew correction
0223       time -= hcalTimeSlew_delay_->delay(std::max(1.0, fc_ampl), HcalTimeSlew::Medium);
0224     }
0225   }
0226   return time;
0227 }