CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
SimpleHBHEPhase1Algo.cc
Go to the documentation of this file.
1 #include <algorithm>
2 
4 
7 
9 
15 
16 // Maximum fractional error for calculating Method 0
17 // pulse containment correction
19 
21  const int firstSampleShift,
22  const int samplesToAdd,
23  const float phaseNS,
24  const float timeShift,
25  const bool correctForPhaseContainment,
26  std::unique_ptr<PulseShapeFitOOTPileupCorrection> m2,
27  std::unique_ptr<HcalDeterministicFit> detFit,
28  std::unique_ptr<MahiFit> mahi)
29  : pulseCorr_(PulseContainmentFractionalError),
30  firstSampleShift_(firstSampleShift),
31  samplesToAdd_(samplesToAdd),
32  phaseNS_(phaseNS),
33  timeShift_(timeShift),
34  runnum_(0),
35  corrFPC_(correctForPhaseContainment),
36  psFitOOTpuCorr_(std::move(m2)),
37  hltOOTpuCorr_(std::move(detFit)),
38  mahiOOTpuCorr_(std::move(mahi))
39 {
40  hcalTimeSlew_delay_ = nullptr;
41 }
42 
44  const edm::EventSetup& es)
45 {
47  es.get<HcalTimeSlewRecord>().get("HBHE", delay);
48  hcalTimeSlew_delay_ = &*delay;
49 
50  runnum_ = r.run();
51  pulseCorr_.beginRun(es);
52 }
53 
55 {
56  runnum_ = 0;
58 }
59 
61  const HcalRecoParam* params,
62  const HcalCalibrations& calibs,
63  const bool isData)
64 {
65  HBHERecHit rh;
66 
67  const HcalDetId channelId(info.id());
68 
69  // Calculate "Method 0" quantities
70  float m0t = 0.f, m0E = 0.f;
71  {
72  int ibeg = static_cast<int>(info.soi()) + firstSampleShift_;
73  if (ibeg < 0)
74  ibeg = 0;
75  const int nSamplesToAdd = params ? params->samplesToAdd() : samplesToAdd_;
76  const double fc_ampl = info.chargeInWindow(ibeg, ibeg + nSamplesToAdd);
77  const bool applyContainment = params ? params->correctForPhaseContainment() : corrFPC_;
78  const float phasens = params ? params->correctionPhaseNS() : phaseNS_;
79  m0E = m0Energy(info, fc_ampl, applyContainment, phasens, nSamplesToAdd);
80  m0E *= hbminusCorrectionFactor(channelId, m0E, isData);
81  m0t = m0Time(info, fc_ampl, calibs, nSamplesToAdd);
82  }
83 
84  // Run "Method 2"
85  float m2t = 0.f, m2E = 0.f, chi2 = -1.f;
86  bool useTriple = false;
88  if (method2)
89  {
90  psFitOOTpuCorr_->setPulseShapeTemplate(theHcalPulseShapes_.getShape(info.recoShape()),
91  !info.hasTimeInfo(),info.nSamples());
92  // "phase1Apply" call below sets m2E, m2t, useTriple, and chi2.
93  // These parameters are pased by non-const reference.
94  method2->phase1Apply(info, m2E, m2t, useTriple, chi2, hcalTimeSlew_delay_);
95  m2E *= hbminusCorrectionFactor(channelId, m2E, isData);
96  }
97 
98  // Run "Method 3"
99  float m3t = 0.f, m3E = 0.f;
100  const HcalDeterministicFit* method3 = hltOOTpuCorr_.get();
101  if (method3)
102  {
103  // "phase1Apply" sets m3E and m3t (pased by non-const reference)
104  method3->phase1Apply(info, m3E, m3t, hcalTimeSlew_delay_);
105  m3E *= hbminusCorrectionFactor(channelId, m3E, isData);
106  }
107 
108  // Run Mahi
109  float m4E = 0.f, m4chi2 = -1.f;
110  float m4T = 0.f;
111  bool m4UseTriple=false;
112 
113  const MahiFit* mahi = mahiOOTpuCorr_.get();
114 
115  if (mahi) {
116  mahiOOTpuCorr_->setPulseShapeTemplate(theHcalPulseShapes_.getShape(info.recoShape()));
117  mahi->phase1Apply(info,m4E,m4T,m4UseTriple,m4chi2,hcalTimeSlew_delay_);
118  m4E *= hbminusCorrectionFactor(channelId, m4E, isData);
119  }
120 
121  // Finally, construct the rechit
122  float rhE = m0E;
123  float rht = m0t;
124  float rhX = -1.f;
125  if (mahi)
126  {
127  rhE = m4E;
128  rht = m4T;
129  rhX = m4chi2;
130  }
131  else if (method2)
132  {
133  rhE = m2E;
134  rht = m2t;
135  rhX = chi2;
136  }
137  else if (method3)
138  {
139  rhE = m3E;
140  rht = m3t;
141  }
142  float tdcTime = info.soiRiseTime();
143  if (!HcalSpecialTimes::isSpecial(tdcTime))
144  tdcTime += timeShift_;
145  rh = HBHERecHit(channelId, rhE, rht, tdcTime);
146  rh.setRawEnergy(m0E);
147  rh.setAuxEnergy(m3E);
148  rh.setChiSquared(rhX);
149 
150  // Set rechit aux words
151  HBHERecHitAuxSetter::setAux(info, &rh);
152 
153  // Set some rechit flags (here, for Method 2/Mahi)
154  if (useTriple || m4UseTriple)
156 
157  return rh;
158 }
159 
161  const float energy,
162  const bool isRealData) const
163 {
164  float corr = 1.f;
165  if (isRealData && runnum_ > 0)
166  if (cell.subdet() == HcalBarrel)
167  {
168  const int ieta = cell.ieta();
169  const int iphi = cell.iphi();
170  corr = hbminus_special_ecorr(ieta, iphi, energy, runnum_);
171  }
172  return corr;
173 }
174 
176  const double fc_ampl,
177  const bool applyContainmentCorrection,
178  const double phaseNs,
179  const int nSamplesToAdd)
180 {
181  int ibeg = static_cast<int>(info.soi()) + firstSampleShift_;
182  if (ibeg < 0)
183  ibeg = 0;
184  double e = info.energyInWindow(ibeg, ibeg + nSamplesToAdd);
185 
186  // Pulse containment correction
187  {
188  double corrFactor = 1.0;
189  if (applyContainmentCorrection)
190  corrFactor = pulseCorr_.get(info.id(), nSamplesToAdd, phaseNs)->getCorrection(fc_ampl);
191  e *= corrFactor;
192  }
193 
194  return e;
195 }
196 
198  const double fc_ampl,
199  const HcalCalibrations& calibs,
200  const int nSamplesToExamine) const
201 {
202  float time = -9999.f; // historic value
203 
204  const unsigned nSamples = info.nSamples();
205  if (nSamples > 2U)
206  {
207  const int soi = info.soi();
208  int ibeg = soi + firstSampleShift_;
209  if (ibeg < 0)
210  ibeg = 0;
211  const int iend = ibeg + nSamplesToExamine;
212  unsigned maxI = info.peakEnergyTS(ibeg, iend);
213  if (maxI < HBHEChannelInfo::MAXSAMPLES)
214  {
215  if (!maxI)
216  maxI = 1U;
217  else if (maxI >= nSamples - 1U)
218  maxI = nSamples - 2U;
219 
220  // The remaining code in this scope emulates
221  // the historic algorithm
222  float t0 = info.tsEnergy(maxI - 1U);
223  float maxA = info.tsEnergy(maxI);
224  float t2 = info.tsEnergy(maxI + 1U);
225 
226  // Handle negative excursions by moving "zero"
227  float minA = t0;
228  if (maxA < minA) minA = maxA;
229  if (t2 < minA) minA=t2;
230  if (minA < 0.f) { maxA-=minA; t0-=minA; t2-=minA; }
231  float wpksamp = (t0 + maxA + t2);
232  if (wpksamp) wpksamp = (maxA + 2.f*t2) / wpksamp;
233  time = (maxI - soi)*25.f + timeshift_ns_hbheho(wpksamp);
234 
235  // Legacy QIE8 timing correction
236  time -= hcalTimeSlew_delay_->delay(std::max(1.0, fc_ampl), HcalTimeSlew::Medium);
237  // Time calibration
238  time -= calibs.timecorr();
239  }
240  }
241  return time;
242 }
const Shape & getShape(int shapeType) const
HBHERecHit reconstruct(const HBHEChannelInfo &info, const HcalRecoParam *params, const HcalCalibrations &calibs, bool isRealData) override
const HcalTimeSlew * hcalTimeSlew_delay_
static const TGPicture * info(bool iBackgroundIsBlack)
float timeshift_ns_hbheho(float wpksamp)
HcalSubdetector subdet() const
get the subdetector
Definition: HcalDetId.h:49
RunNumber_t run() const
Definition: RunBase.h:40
const HcalPulseContainmentCorrection * get(const HcalDetId &detId, int toAdd, float fixedphase_ns)
bool hasTimeInfo() const
std::unique_ptr< HcalDeterministicFit > hltOOTpuCorr_
void phase1Apply(const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &reconstructedTime, bool &useTriple, float &chi2, const HcalTimeSlew *hcalTimeSlew_delay) const
void beginRun(const edm::Run &, const edm::EventSetup &) override
float PulseContainmentFractionalError
bool correctForPhaseContainment() const
Definition: HcalRecoParam.h:29
void phase1Apply(const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &reconstructedTime, bool &useTriple, float &chi2, const HcalTimeSlew *hcalTimeSlew_delay) const
Definition: MahiFit.cc:45
void setFlagField(uint32_t value, int base, int width=1)
Definition: CaloRecHit.cc:20
double delay(double fC, BiasSetting bias=Medium) const
Returns the amount (ns) by which a pulse of the given number of fC will be delayed by the timeslew ef...
Definition: HcalTimeSlew.cc:14
unsigned soi() const
HcalPulseContainmentManager pulseCorr_
float hbminus_special_ecorr(int ieta, int iphi, double energy, int runnum)
Special energy correction for some HB- cells.
HcalDetId id() const
#define constexpr
static const unsigned MAXSAMPLES
bool isSpecial(const float t)
auto const T2 &decltype(t1.eta()) t2
Definition: deltaR.h:16
float m0Time(const HBHEChannelInfo &info, double reconstructedCharge, const HcalCalibrations &calibs, int nSamplesToExamine) const
float soiRiseTime() const
int ieta() const
get the cell ieta
Definition: HcalDetId.h:56
double energyInWindow(const unsigned begin, const unsigned end) const
double f[11][100]
double chargeInWindow(const unsigned begin, const unsigned end) const
SimpleHBHEPhase1Algo(int firstSampleShift, int samplesToAdd, float phaseNS, float timeShift, bool correctForPhaseContainment, std::unique_ptr< PulseShapeFitOOTPileupCorrection > m2, std::unique_ptr< HcalDeterministicFit > detFit, std::unique_ptr< MahiFit > mahi)
int recoShape() const
unsigned int samplesToAdd() const
Definition: HcalRecoParam.h:33
float hbminusCorrectionFactor(const HcalDetId &cell, float energy, bool isRealData) const
constexpr size_t nSamples
std::unique_ptr< MahiFit > mahiOOTpuCorr_
float correctionPhaseNS() const
Definition: HcalRecoParam.h:31
JetCorrectorParameters corr
Definition: classes.h:5
int iphi() const
get the cell iphi
Definition: HcalDetId.cc:103
HcalPulseShapes theHcalPulseShapes_
double tsEnergy(const unsigned ts) const
double timecorr() const
get time correction factor
const T & get() const
Definition: EventSetup.h:58
void beginRun(edm::EventSetup const &es)
float m0Energy(const HBHEChannelInfo &info, double reconstructedCharge, bool applyContainmentCorrection, double phaseNS, int nSamplesToAdd)
static void setAux(const HBHEChannelInfo &info, HBHERecHit *rechit)
void setAuxEnergy(const float en)
Definition: HBHERecHit.h:33
void setChiSquared(const float chi2)
Definition: HBHERecHit.h:27
unsigned nSamples() const
std::unique_ptr< PulseShapeFitOOTPileupCorrection > psFitOOTpuCorr_
void setRawEnergy(const float en)
Definition: HBHERecHit.h:30
unsigned peakEnergyTS(const unsigned begin, const unsigned end) const
def move(src, dest)
Definition: eostools.py:510
void phase1Apply(const HBHEChannelInfo &channelData, float &reconstructedEnergy, float &reconstructedTime, const HcalTimeSlew *hcalTimeSlew_delay) const
Definition: Run.h:43