CMS 3D CMS Logo

List of all members | Classes | Public Member Functions | Protected Attributes
EcalUncalibRecHitRatioMethodAlgo< C > Class Template Reference

#include <EcalUncalibRecHitRatioMethodAlgo.h>

Classes

struct  CalculatedRecHit
 
struct  Ratio
 
struct  Tmax
 

Public Member Functions

void computeAmplitude (std::vector< double > &amplitudeFitParameters)
 
double computeAmplitudeImpl (std::vector< double > &amplitudeFitParameters, double, double)
 
void computeTime (std::vector< double > &timeFitParameters, std::pair< double, double > &timeFitLimits, std::vector< double > &amplitudeFitParameters)
 
bool fixMGPAslew (const C &dataFrame)
 
CalculatedRecHit getCalculatedRecHit ()
 
void init (const C &dataFrame, const EcalSampleMask &sampleMask, const double *pedestals, const double *pedestalRMSes, const double *gainRatios)
 
EcalUncalibratedRecHit makeRecHit (const C &dataFrame, const EcalSampleMask &sampleMask, const double *pedestals, const double *pedestalRMSes, const double *gainRatios, std::vector< double > &timeFitParameters, std::vector< double > &amplitudeFitParameters, std::pair< double, double > &timeFitLimits)
 

Protected Attributes

std::array< double, C::MAXSAMPLES > amplitudeErrors_
 
std::array< double, C::MAXSAMPLES > amplitudeIE2_
 
std::array< double, C::MAXSAMPLES > amplitudes_
 
double ampMaxError_
 
CalculatedRecHit calculatedRechit_
 
int num_
 
double pedestal_
 
EcalSampleMask sampleMask_
 
DetId theDetId_
 
std::array< bool, C::MAXSAMPLES > useless_
 

Detailed Description

template<class C>
class EcalUncalibRecHitRatioMethodAlgo< C >

Template used to compute amplitude, pedestal, time jitter, chi2 of a pulse using a ratio method

Author
A. Ledovskoy (Design) - M. Balazs (Implementation)

Definition at line 31 of file EcalUncalibRecHitRatioMethodAlgo.h.

Member Function Documentation

template<class C >
void EcalUncalibRecHitRatioMethodAlgo< C >::computeAmplitude ( std::vector< double > &  amplitudeFitParameters)

Definition at line 560 of file EcalUncalibRecHitRatioMethodAlgo.h.

Referenced by EcalUncalibRecHitWorkerGlobal::run(), and EcalUncalibRecHitWorkerMultiFit::run().

561  {
562 
564  computeAmplitudeImpl(amplitudeFitParameters, 1., 1.);
565 
566 }
double computeAmplitudeImpl(std::vector< double > &amplitudeFitParameters, double, double)
template<class C >
double EcalUncalibRecHitRatioMethodAlgo< C >::computeAmplitudeImpl ( std::vector< double > &  amplitudeFitParameters,
double  corr4,
double  corr6 
)

Definition at line 569 of file EcalUncalibRecHitRatioMethodAlgo.h.

References funct::abs(), alpha, pfBoostedDoubleSVAK8TagInfos_cfi::beta, MillePedeFileConverter_cfg::e, f, myMath::fast_expf(), myMath::fast_logf(), and mps_fire::i.

570  {
572  // //
573  // CALCULATE AMPLITUDE //
574  // //
576 
577  double alpha = amplitudeFitParameters[0];
578  double beta = amplitudeFitParameters[1];
579 
580  // calculate pedestal, again
581 
582  double pedestalLimit = calculatedRechit_.timeMax - (alpha * beta) - 1.0;
583  double sumA = 0;
584  double sumF = 0;
585  double sumAF = 0;
586  double sumFF = 0;
587  double sum1 = 0;
588  for (unsigned int i = 0; i < amplitudes_.size(); i++) {
589  if (useless_[i]) continue;
590  double inverr2 = amplitudeIE2_[i];
591  double f = 0;
592  double termOne = 1 + (i - calculatedRechit_.timeMax) / (alpha * beta);
593  if (termOne > 1.e-5)
594  f = myMath::fast_expf(alpha * myMath::fast_logf(termOne) -
595  (i - calculatedRechit_.timeMax) / beta);
596 
597  // apply range of interesting samples
598 
599  if ((i < pedestalLimit) ||
600  (f > 0.6 * corr6 && i <= calculatedRechit_.timeMax) ||
601  (f > 0.4 * corr4 && i >= calculatedRechit_.timeMax)) {
602  sum1 += inverr2;
603  sumA += amplitudes_[i] * inverr2;
604  sumF += f * inverr2;
605  sumAF += (f * inverr2) * amplitudes_[i];
606  sumFF += f * (f * inverr2);
607  }
608  }
609 
610  double amplitudeMax = 0;
611  if (sum1 > 0) {
612  double denom = sumFF * sum1 - sumF * sumF;
613  if (std::abs(denom) > 1.0e-20) {
614  amplitudeMax = (sumAF * sum1 - sumA * sumF) / denom;
615  }
616  }
617  return amplitudeMax;
618 }
float alpha
Definition: AMPTWrapper.h:95
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
double f[11][100]
std::array< double, C::MAXSAMPLES > amplitudeIE2_
std::array< double, C::MAXSAMPLES > amplitudes_
std::array< bool, C::MAXSAMPLES > useless_
float fast_expf(float x)
float fast_logf(float x)
template<class C >
void EcalUncalibRecHitRatioMethodAlgo< C >::computeTime ( std::vector< double > &  timeFitParameters,
std::pair< double, double > &  timeFitLimits,
std::vector< double > &  amplitudeFitParameters 
)

Definition at line 224 of file EcalUncalibRecHitRatioMethodAlgo.h.

References alpha, pfBoostedDoubleSVAK8TagInfos_cfi::beta, vertices_cff::chi2, EcalUncalibRecHitRatioMethodAlgo< C >::Tmax::chi2, MillePedeFileConverter_cfg::e, f, myMath::fast_expf(), myMath::fast_logf(), mps_fire::i, EcalUncalibRecHitRatioMethodAlgo< C >::Tmax::index, gen::k, SiStripPI::max, PFRecoTauDiscriminationByIsolation_cfi::offset, Gflash::Rmax, Gflash::Rmin, mathSSE::sqrt(), trackingPlots::stat, and tmax.

Referenced by EcalUncalibRecHitWorkerGlobal::run(), and EcalUncalibRecHitWorkerMultiFit::run().

227  {
229  // //
230  // RATIO METHOD FOR TIME STARTS HERE //
231  // //
233  double ampMaxAlphaBeta = 0;
234  double tMaxAlphaBeta = 5;
235  double tMaxErrorAlphaBeta = 999;
236  double tMaxRatio = 5;
237  double tMaxErrorRatio = 999;
238 
239  double sumAA = 0;
240  double sumA = 0;
241  double sum1 = 0;
242  double sum0 = 0;
243  double sumAf = 0;
244  double sumff = 0;
245  double NullChi2 = 0;
246 
247  // null hypothesis = no pulse, pedestal only
248  for (unsigned int i = 0; i < amplitudes_.size(); i++) {
249  if (useless_[i]) continue;
250  double inverr2 = amplitudeIE2_[i];
251  sum0 += 1;
252  sum1 += inverr2;
253  sumA += amplitudes_[i] * inverr2;
254  sumAA += amplitudes_[i] * (amplitudes_[i] * inverr2);
255  }
256  if (sum0 > 0) {
257  NullChi2 = (sumAA - sumA * sumA / sum1) / sum0;
258  } else {
259  // not enough samples to reconstruct the pulse
260  return;
261  }
262 
263  // Make all possible Ratio's based on any pair of samples i and j
264  // (j>i) with positive amplitudes_
265  //
266  // Ratio[k] = Amp[i]/Amp[j]
267  // where Amp[i] is pedestal subtracted ADC value in a time sample [i]
268  //
269  double alphabeta = amplitudeFitParameters[0] * amplitudeFitParameters[1];
270  double invalphabeta = 1. / alphabeta;
271  double alpha = amplitudeFitParameters[0];
272  double beta = amplitudeFitParameters[1];
273 
274  Ratio ratios_[C::MAXSAMPLES * (C::MAXSAMPLES - 1) / 2];
275  unsigned int ratios_size = 0;
276 
277  double Rlim[amplitudes_.size()];
278  for (unsigned int k = 1; k != amplitudes_.size(); ++k)
279  Rlim[k] = myMath::fast_expf(double(k) / beta) - 0.001;
280 
281  double relErr2[amplitudes_.size()];
282  double invampl[amplitudes_.size()];
283  for (unsigned int i = 0; i < amplitudes_.size(); i++) {
284  invampl[i] = (useless_[i]) ? 0 : 1. / amplitudes_[i];
285  relErr2[i] = (useless_[i]) ? 0 : (amplitudeErrors_[i] * invampl[i]) *
286  (amplitudeErrors_[i] * invampl[i]);
287  }
288 
289  for (unsigned int i = 0; i < amplitudes_.size() - 1; i++) {
290  if (useless_[i]) continue;
291  for (unsigned int j = i + 1; j < amplitudes_.size(); j++) {
292  if (useless_[j]) continue;
293 
294  if (amplitudes_[i] > 1 && amplitudes_[j] > 1) {
295 
296  // ratio
297  double Rtmp = amplitudes_[i] / amplitudes_[j];
298 
299  // error^2 due to stat fluctuations of time samples
300  // (uncorrelated for both samples)
301 
302  double err1 = Rtmp * Rtmp * (relErr2[i] + relErr2[j]);
303 
304  // error due to fluctuations of pedestal (common to both samples)
305  double stat;
306  if (num_ > 0)
307  stat = num_; // num presampeles used to compute pedestal
308  else
309  stat = 1; // pedestal from db
310  double err2 = amplitudeErrors_[j] * (amplitudes_[i] - amplitudes_[j]) *
311  (invampl[j] * invampl[j]); // /sqrt(stat);
312  err2 *= err2 / stat;
313 
314  //error due to integer round-down. It is relevant to low
315  //amplitudes_ in gainID=1 and negligible otherwise.
316  double err3 = (0.289 * 0.289) * (invampl[j] * invampl[j]);
317 
318  double totalError = std::sqrt(err1 + err2 + err3);
319 
320  // don't include useless ratios
321  if (totalError < 1.0 && Rtmp > 0.001 && Rtmp < Rlim[j - i]) {
322  Ratio currentRatio = { i, (j - i), Rtmp, totalError };
323  ratios_[ratios_size++] = currentRatio;
324  }
325  }
326  }
327  }
328 
329  // No useful ratios, return zero amplitude and no time measurement
330  if (0 == ratios_size) return;
331 
332  // std::array < Tmax, C::MAXSAMPLES*(C::MAXSAMPLES-1)/2 > times_;
333  Tmax timesAB_[C::MAXSAMPLES * (C::MAXSAMPLES - 1) / 2];
334  unsigned int timesAB_size = 0;
335 
336  // make a vector of Tmax measurements that correspond to each ratio
337  // and based on Alpha-Beta parameterization of the pulse shape
338 
339  for (unsigned int i = 0; i < ratios_size; i++) {
340 
341  double stepOverBeta = double(ratios_[i].step) / beta;
342  double offset = double(ratios_[i].index) + alphabeta;
343 
344  double Rmin = ratios_[i].value - ratios_[i].error;
345  if (Rmin < 0.001) Rmin = 0.001;
346 
347  double Rmax = ratios_[i].value + ratios_[i].error;
348  double RLimit = Rlim[ratios_[i].step];
349  if (Rmax > RLimit) Rmax = RLimit;
350 
351  double time1 =
352  offset -
353  ratios_[i].step /
354  (myMath::fast_expf((stepOverBeta - myMath::fast_logf(Rmin)) /
355  alpha) - 1.0);
356  double time2 =
357  offset -
358  ratios_[i].step /
359  (myMath::fast_expf((stepOverBeta - myMath::fast_logf(Rmax)) /
360  alpha) - 1.0);
361 
362  // this is the time measurement based on the ratios[i]
363  double tmax = 0.5 * (time1 + time2);
364  double tmaxerr = 0.5 * std::sqrt((time1 - time2) * (time1 - time2));
365 
366  // calculate chi2
367  sumAf = 0;
368  sumff = 0;
369  int itmin = std::max(-1, int(std::floor(tmax - alphabeta)));
370  double loffset = (double(itmin) - tmax) * invalphabeta;
371  for (unsigned int it = itmin + 1; it < amplitudes_.size(); it++) {
372  loffset += invalphabeta;
373  if (useless_[it]) continue;
374  double inverr2 = amplitudeIE2_[it];
375  // double offset = (double(it) - tmax)/alphabeta;
376  double term1 = 1.0 + loffset;
377  // assert(term1>1e-6);
378  double f =
379  (term1 > 1e-6)
380  ? myMath::fast_expf(alpha * (myMath::fast_logf(term1) - loffset))
381  : 0;
382  sumAf += amplitudes_[it] * (f * inverr2);
383  sumff += f * (f * inverr2);
384  }
385 
386  double chi2 = sumAA;
387  double amp = 0;
388  if (sumff > 0) {
389  chi2 = sumAA - sumAf * (sumAf / sumff);
390  amp = (sumAf / sumff);
391  }
392  chi2 /= sum0;
393 
394  // choose reasonable measurements. One might argue what is
395  // reasonable and what is not.
396  if (chi2 > 0 && tmaxerr > 0 && tmax > 0) {
397  Tmax currentTmax = { ratios_[i].index, ratios_[i].step, tmax, tmaxerr,
398  amp, chi2 };
399  timesAB_[timesAB_size++] = currentTmax;
400  }
401  }
402 
403  // no reasonable time measurements!
404  if (0 == timesAB_size) return;
405 
406  // find minimum chi2
407  double chi2min = 1.0e+9;
408  //double timeMinimum = 5;
409  //double errorMinimum = 999;
410  for (unsigned int i = 0; i < timesAB_size; i++) {
411  if (timesAB_[i].chi2 < chi2min) {
412  chi2min = timesAB_[i].chi2;
413  //timeMinimum = timesAB_[i].value;
414  //errorMinimum = timesAB_[i].error;
415  }
416  }
417 
418  // make a weighted average of tmax measurements with "small" chi2
419  // (within 1 sigma of statistical uncertainty :-)
420  double chi2Limit = chi2min + 1.0;
421  double time_max = 0;
422  double time_wgt = 0;
423  for (unsigned int i = 0; i < timesAB_size; i++) {
424  if (timesAB_[i].chi2 < chi2Limit) {
425  double inverseSigmaSquared =
426  1.0 / (timesAB_[i].error * timesAB_[i].error);
427  time_wgt += inverseSigmaSquared;
428  time_max += timesAB_[i].value * inverseSigmaSquared;
429  }
430  }
431 
432  tMaxAlphaBeta = time_max / time_wgt;
433  tMaxErrorAlphaBeta = 1.0 / sqrt(time_wgt);
434 
435  // find amplitude and chi2
436  sumAf = 0;
437  sumff = 0;
438  for (unsigned int i = 0; i < amplitudes_.size(); i++) {
439  if (useless_[i]) continue;
440  double inverr2 = amplitudeIE2_[i];
441  double offset = (double(i) - tMaxAlphaBeta) * invalphabeta;
442  double term1 = 1.0 + offset;
443  if (term1 > 1e-6) {
444  double f = myMath::fast_expf(alpha * (myMath::fast_logf(term1) - offset));
445  sumAf += amplitudes_[i] * (f * inverr2);
446  sumff += f * (f * inverr2);
447  }
448  }
449 
450  if (sumff > 0) {
451  ampMaxAlphaBeta = sumAf / sumff;
452  double chi2AlphaBeta = (sumAA - sumAf * sumAf / sumff) / sum0;
453  if (chi2AlphaBeta > NullChi2) {
454  // null hypothesis is better
455  return;
456  }
457 
458  } else {
459  // no visible pulse here
460  return;
461  }
462 
463  // if we got to this point, we have a reconstructied Tmax
464  // using RatioAlphaBeta Method. To summarize:
465  //
466  // tMaxAlphaBeta - Tmax value
467  // tMaxErrorAlphaBeta - error on Tmax, but I would not trust it
468  // ampMaxAlphaBeta - amplitude of the pulse
469  // ampMaxError_ - uncertainty of the time sample with max amplitude
470  //
471 
472  // Do Ratio's Method with "large" pulses
473  if (ampMaxAlphaBeta / ampMaxError_ > 5.0) {
474 
475  // make a vector of Tmax measurements that correspond to each
476  // ratio. Each measurement have it's value and the error
477 
478  double time_max = 0;
479  double time_wgt = 0;
480 
481  for (unsigned int i = 0; i < ratios_size; i++) {
482 
483  if (ratios_[i].step == 1 && ratios_[i].value >= timeFitLimits.first &&
484  ratios_[i].value <= timeFitLimits.second) {
485 
486  double time_max_i = ratios_[i].index;
487 
488  // calculate polynomial for Tmax
489 
490  double u = timeFitParameters[timeFitParameters.size() - 1];
491  for (int k = timeFitParameters.size() - 2; k >= 0; k--) {
492  u = u * ratios_[i].value + timeFitParameters[k];
493  }
494 
495  // calculate derivative for Tmax error
496  double du = (timeFitParameters.size() - 1) *
497  timeFitParameters[timeFitParameters.size() - 1];
498  for (int k = timeFitParameters.size() - 2; k >= 1; k--) {
499  du = du * ratios_[i].value + k * timeFitParameters[k];
500  }
501 
502  // running sums for weighted average
503  double errorsquared = ratios_[i].error * ratios_[i].error * du * du;
504  if (errorsquared > 0) {
505 
506  time_max += (time_max_i - u) / errorsquared;
507  time_wgt += 1.0 / errorsquared;
508  // Tmax currentTmax =
509  // { ratios_[i].index, 1, (time_max_i - u),
510  //sqrt(errorsquared),0,1 };
511  // times_.push_back(currentTmax);
512 
513  }
514  }
515  }
516 
517  // calculate weighted average of all Tmax measurements
518  if (time_wgt > 0) {
519  tMaxRatio = time_max / time_wgt;
520  tMaxErrorRatio = 1.0 / sqrt(time_wgt);
521 
522  // combine RatioAlphaBeta and Ratio Methods
523 
524  if (ampMaxAlphaBeta / ampMaxError_ > 10.0) {
525 
526  // use pure Ratio Method
527  calculatedRechit_.timeMax = tMaxRatio;
528  calculatedRechit_.timeError = tMaxErrorRatio;
529 
530  } else {
531 
532  // combine two methods
534  (tMaxAlphaBeta * (10.0 - ampMaxAlphaBeta / ampMaxError_) +
535  tMaxRatio * (ampMaxAlphaBeta / ampMaxError_ - 5.0)) / 5.0;
537  (tMaxErrorAlphaBeta * (10.0 - ampMaxAlphaBeta / ampMaxError_) +
538  tMaxErrorRatio * (ampMaxAlphaBeta / ampMaxError_ - 5.0)) / 5.0;
539 
540  }
541 
542  } else {
543 
544  // use RatioAlphaBeta Method
545  calculatedRechit_.timeMax = tMaxAlphaBeta;
546  calculatedRechit_.timeError = tMaxErrorAlphaBeta;
547 
548  }
549 
550  } else {
551 
552  // use RatioAlphaBeta Method
553  calculatedRechit_.timeMax = tMaxAlphaBeta;
554  calculatedRechit_.timeError = tMaxErrorAlphaBeta;
555 
556  }
557 }
float alpha
Definition: AMPTWrapper.h:95
std::array< double, C::MAXSAMPLES > amplitudeErrors_
T sqrt(T t)
Definition: SSEVec.h:18
double f[11][100]
Definition: value.py:1
std::array< double, C::MAXSAMPLES > amplitudeIE2_
std::array< double, C::MAXSAMPLES > amplitudes_
int k[5][pyjets_maxn]
static const double tmax[3]
std::array< bool, C::MAXSAMPLES > useless_
float fast_expf(float x)
const double Rmax[kNumberCalorimeter]
step
Definition: StallMonitor.cc:94
const double Rmin[kNumberCalorimeter]
float fast_logf(float x)
template<class C>
bool EcalUncalibRecHitRatioMethodAlgo< C >::fixMGPAslew ( const C &  dataFrame)

Definition at line 189 of file EcalUncalibRecHitRatioMethodAlgo.h.

References MillePedeFileConverter_cfg::e, and mps_fire::result.

Referenced by EcalUncalibRecHitWorkerRatio::run(), EcalUncalibRecHitWorkerGlobal::run(), and EcalUncalibRecHitWorkerMultiFit::run().

189  {
190 
191  // This fuction finds sample(s) preceeding gain switching and
192  // inflates errors on this sample, therefore, making this sample
193  // invisible for Ratio Method. Only gain switching DOWN is
194  // considered Only gainID=1,2,3 are considered. In case of the
195  // saturation (gainID=0), we keep "distorted" sample because it is
196  // the only chance to make time measurement; the qualilty of it will
197  // be bad anyway.
198 
199  bool result = false;
200 
201  int GainIdPrev;
202  int GainIdNext;
203  for (int iSample = 1; iSample < C::MAXSAMPLES; iSample++) {
204 
205  // only use samples which are desired
206  if (!sampleMask_.useSample(iSample, theDetId_)) continue;
207 
208  GainIdPrev = dataFrame.sample(iSample - 1).gainId();
209  GainIdNext = dataFrame.sample(iSample).gainId();
210  if (GainIdPrev >= 1 && GainIdPrev <= 3 && GainIdNext >= 1 &&
211  GainIdNext <= 3 && GainIdPrev < GainIdNext) {
212  amplitudes_[iSample - 1] = 0;
213  amplitudeErrors_[iSample - 1] = 1e+9;
214  amplitudeIE2_[iSample - 1] = 0;
215  useless_[iSample - 1] = true;
216  result = true;
217  }
218  }
219  return result;
220 
221 }
std::array< double, C::MAXSAMPLES > amplitudeErrors_
bool useSample(const int sampleId, DetId &theCrystal) const
std::array< double, C::MAXSAMPLES > amplitudeIE2_
std::array< double, C::MAXSAMPLES > amplitudes_
std::array< bool, C::MAXSAMPLES > useless_
template<class C>
CalculatedRecHit EcalUncalibRecHitRatioMethodAlgo< C >::getCalculatedRecHit ( )
inline
template<class C>
void EcalUncalibRecHitRatioMethodAlgo< C >::init ( const C &  dataFrame,
const EcalSampleMask sampleMask,
const double *  pedestals,
const double *  pedestalRMSes,
const double *  gainRatios 
)

Definition at line 95 of file EcalUncalibRecHitRatioMethodAlgo.h.

References funct::abs(), MillePedeFileConverter_cfg::e, and simplePhotonAnalyzer_cfi::sample.

Referenced by EcalUncalibRecHitWorkerGlobal::run(), and EcalUncalibRecHitWorkerMultiFit::run().

99  {
100  sampleMask_ = sampleMask;
101  theDetId_ = DetId(dataFrame.id().rawId());
102 
106 
107  // to obtain gain 12 pedestal:
108  // -> if it's in gain 12, use first sample
109  // --> average it with second sample if in gain 12 and 3-sigma-noise
110  // compatible (better LF noise cancellation)
111  // -> else use pedestal from database
112  pedestal_ = 0;
113  num_ = 0;
114  if (dataFrame.sample(0).gainId() == 1 &&
116  pedestal_ += double(dataFrame.sample(0).adc());
117  num_++;
118  }
119  if (num_ != 0 && dataFrame.sample(1).gainId() == 1 &&
121  std::abs(dataFrame.sample(1).adc() - dataFrame.sample(0).adc()) <
122  3 * pedestalRMSes[0]) {
123  pedestal_ += double(dataFrame.sample(1).adc());
124  num_++;
125  }
126  if (num_ != 0)
127  pedestal_ /= num_;
128  else
129  pedestal_ = pedestals[0];
130 
131  // fill vector of amplitudes, pedestal subtracted and vector
132  // of amplitude uncertainties Also, find the uncertainty of a
133  // sample with max amplitude. We will use it later.
134 
135  ampMaxError_ = 0;
136  double ampMaxValue = -1000;
137 
138  // ped-subtracted and gain-renormalized samples. It is VERY
139  // IMPORTANT to have samples one clock apart which means to
140  // have vector size equal to MAXSAMPLES
141  double sample;
142  double sampleError;
143  int GainId;
144  for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
145 
146  GainId = dataFrame.sample(iSample).gainId();
147 
148  bool bad = false;
149  // only use normally samples which are desired; if sample not to be used
150  // inflate error so won't generate ratio considered for the measurement
151  if (!sampleMask_.useSample(iSample, theDetId_)) {
152  sample = 0;
153  sampleError = 0;
154  bad = true;
155  } else if (GainId == 1) {
156  sample = double(dataFrame.sample(iSample).adc() - pedestal_);
157  sampleError = pedestalRMSes[0];
158  } else if (GainId == 2 || GainId == 3) {
159  sample =
160  (double(dataFrame.sample(iSample).adc() - pedestals[GainId - 1])) *
161  gainRatios[GainId - 1];
162  sampleError = pedestalRMSes[GainId - 1] * gainRatios[GainId - 1];
163  } else {
164  sample = 0; // GainId=0 case falls here, from saturation
165  sampleError = 0; // inflate error so won't generate ratio considered for
166  // the measurement
167  bad = true;
168  }
169 
170  useless_[iSample] = (sampleError <= 0) | bad;
171  amplitudes_[iSample] = sample;
172  // inflate error for useless samples
173  amplitudeErrors_[iSample] = useless_[iSample] ? 1e+9 : sampleError;
174  amplitudeIE2_[iSample] =
175  useless_[iSample]
176  ? 0
177  : 1 / (amplitudeErrors_[iSample] * amplitudeErrors_[iSample]);
178  if (sampleError > 0) {
179  if (ampMaxValue < sample) {
180  ampMaxValue = sample;
181  ampMaxError_ = sampleError;
182  }
183  }
184 
185  }
186 
187 }
std::array< double, C::MAXSAMPLES > amplitudeErrors_
bool useSample(const int sampleId, DetId &theCrystal) const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::array< double, C::MAXSAMPLES > amplitudeIE2_
std::array< double, C::MAXSAMPLES > amplitudes_
Definition: DetId.h:18
std::array< bool, C::MAXSAMPLES > useless_
template<class C>
EcalUncalibratedRecHit EcalUncalibRecHitRatioMethodAlgo< C >::makeRecHit ( const C &  dataFrame,
const EcalSampleMask sampleMask,
const double *  pedestals,
const double *  pedestalRMSes,
const double *  gainRatios,
std::vector< double > &  timeFitParameters,
std::vector< double > &  amplitudeFitParameters,
std::pair< double, double > &  timeFitLimits 
)

Definition at line 621 of file EcalUncalibRecHitRatioMethodAlgo.h.

References init.

Referenced by EcalUncalibRecHitWorkerRatio::run().

626  {
627 
628  init(dataFrame, sampleMask, pedestals, pedestalRMSes, gainRatios);
629  computeTime(timeFitParameters, timeFitLimits, amplitudeFitParameters);
630  computeAmplitude(amplitudeFitParameters);
631 
632  // 1st parameters is id
633  //
634  // 2nd parameters is amplitude. It is calculated by this method.
635  //
636  // 3rd parameter is pedestal. It is not calculated. This method
637  // relies on input parameters for pedestals and gain ratio. Return
638  // zero.
639  //
640  // 4th parameter is jitter which is a bad choice to call Tmax. It is
641  // calculated by this method (in 25 nsec clock units)
642  //
643  // GF subtract 5 so that jitter==0 corresponds to synchronous hit
644  //
645  //
646  // 5th parameter is chi2. It is possible to calculate chi2 for
647  // Tmax. It is possible to calculate chi2 for Amax. However, these
648  // values are not very useful without TmaxErr and AmaxErr. This
649  // method can return one value for chi2 but there are 4 different
650  // parameters that have useful information about the quality of Amax
651  // ans Tmax. For now we can return TmaxErr. The quality of Tmax and
652  // Amax can be judged from the magnitude of TmaxErr
653 
657 }
void computeAmplitude(std::vector< double > &amplitudeFitParameters)
void computeTime(std::vector< double > &timeFitParameters, std::pair< double, double > &timeFitLimits, std::vector< double > &amplitudeFitParameters)
void init(const C &dataFrame, const EcalSampleMask &sampleMask, const double *pedestals, const double *pedestalRMSes, const double *gainRatios)

Member Data Documentation

template<class C>
std::array<double, C::MAXSAMPLES> EcalUncalibRecHitRatioMethodAlgo< C >::amplitudeErrors_
protected

Definition at line 83 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>
std::array<double, C::MAXSAMPLES> EcalUncalibRecHitRatioMethodAlgo< C >::amplitudeIE2_
protected

Definition at line 84 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>
std::array<double, C::MAXSAMPLES> EcalUncalibRecHitRatioMethodAlgo< C >::amplitudes_
protected

Definition at line 82 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>
double EcalUncalibRecHitRatioMethodAlgo< C >::ampMaxError_
protected

Definition at line 89 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>
CalculatedRecHit EcalUncalibRecHitRatioMethodAlgo< C >::calculatedRechit_
protected

Definition at line 91 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>
int EcalUncalibRecHitRatioMethodAlgo< C >::num_
protected

Definition at line 88 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>
double EcalUncalibRecHitRatioMethodAlgo< C >::pedestal_
protected

Definition at line 87 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>
EcalSampleMask EcalUncalibRecHitRatioMethodAlgo< C >::sampleMask_
protected

Definition at line 80 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>
DetId EcalUncalibRecHitRatioMethodAlgo< C >::theDetId_
protected

Definition at line 81 of file EcalUncalibRecHitRatioMethodAlgo.h.

template<class C>
std::array<bool, C::MAXSAMPLES> EcalUncalibRecHitRatioMethodAlgo< C >::useless_
protected

Definition at line 85 of file EcalUncalibRecHitRatioMethodAlgo.h.