CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
List of all members | Public Member Functions | Private Member Functions | Private Attributes
EcalUncalibRecHitFixedAlphaBetaAlgo< C > Class Template Reference

#include <EcalUncalibRecHitFixedAlphaBetaAlgo.h>

Inheritance diagram for EcalUncalibRecHitFixedAlphaBetaAlgo< C >:
EcalUncalibRecHitRecAbsAlgo< C >

Public Member Functions

 EcalUncalibRecHitFixedAlphaBetaAlgo ()
 
 EcalUncalibRecHitFixedAlphaBetaAlgo (int n_iter, int n_bef_max=1, int n_aft_max=3, float sigma_ped=1.1)
 
EcalUncalibratedRecHit makeRecHit (const C &dataFrame, const double *pedestals, const double *gainRatios, const EcalWeightSet::EcalWeightMatrix **weights, const EcalWeightSet::EcalChi2WeightMatrix **chi2Matrix) override
 Compute parameters. More...
 
void SetAlphaBeta (double alpha, double beta)
 
void SetDynamicPedestal (bool dyn_pede)
 
void SetMinAmpl (double ampl)
 
 ~EcalUncalibRecHitFixedAlphaBetaAlgo () override
 
- Public Member Functions inherited from EcalUncalibRecHitRecAbsAlgo< C >
virtual ~EcalUncalibRecHitRecAbsAlgo ()=default
 Constructor. More...
 

Private Member Functions

void InitFitParameters (double *samples, int max_sample)
 
float PerformAnalyticFit (double *samples, int max_sample)
 
double pulseShapeFunction (double t)
 

Private Attributes

double alfabeta_
 
float alpha_table_ [36][1701]
 
float beta_table_ [36][1701]
 
CLHEP::HepSymMatrix DM1_
 
bool doFit_
 
bool dyn_pedestal
 
double fAlpha_
 
double fAmp_max_
 
double fBeta_
 
int fNb_iter_
 
int fNum_samp_after_max_
 
int fNum_samp_bef_max_
 
double fPed_max_
 
float fSigma_ped
 
double fTim_max_
 
double MinAmpl_
 
CLHEP::HepVector temp_
 
double un_sur_sigma
 

Additional Inherited Members

- Public Types inherited from EcalUncalibRecHitRecAbsAlgo< C >
enum  { nWeightsRows = 3, iAmplitude = 0, iPedestal = 1, iTime = 2 }
 

Detailed Description

template<class C>
class EcalUncalibRecHitFixedAlphaBetaAlgo< C >

Template used to compute amplitude, pedestal, time jitter, chi2 of a pulse using an analytical function fit, with the pulse parameters alpha and beta fixed. It follows a fast fit algorithms devolped on test beam data by P. Jarry If the pedestal is not given, it is calculated from the first 2 pre-samples; FIXME: conversion gainID (1,2,3) with gain factor (12,6,1) is hardcoded here !

Author
A.Ghezzi

Definition at line 33 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Constructor & Destructor Documentation

Definition at line 63 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

64  : fAlpha_(0.),
65  fBeta_(0.),
66  fAmp_max_(-1.),
67  fTim_max_(-1),
68  fPed_max_(0),
69  alfabeta_(0),
70  fNb_iter_(4),
73  fSigma_ped(1.1),
74  DM1_(3),
75  temp_(3) {
76  un_sur_sigma = 1. / double(fSigma_ped);
77  for (int i = 0; i < 36; i++) {
78  for (int j = 0; j < 1701; j++) {
79  alpha_table_[i][j] = 1.2;
80  beta_table_[i][j] = 1.7;
81  }
82  }
83  doFit_ = false;
84  MinAmpl_ = 16;
85  dyn_pedestal = true;
86  }
template<class C>
EcalUncalibRecHitFixedAlphaBetaAlgo< C >::EcalUncalibRecHitFixedAlphaBetaAlgo ( int  n_iter,
int  n_bef_max = 1,
int  n_aft_max = 3,
float  sigma_ped = 1.1 
)
inline

Definition at line 87 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

88  : fAlpha_(0.), fBeta_(0.), fAmp_max_(-1.), fTim_max_(-1), fPed_max_(0), alfabeta_(0), DM1_(3), temp_(3) {
89  fNb_iter_ = n_iter;
90  fNum_samp_bef_max_ = n_bef_max;
91  fNum_samp_after_max_ = n_aft_max;
92  fSigma_ped = sigma_ped;
93  un_sur_sigma = 1. / double(fSigma_ped);
94  for (int i = 0; i < 36; i++) {
95  for (int j = 0; j < 1701; j++) {
96  alpha_table_[i][j] = 1.2;
97  beta_table_[i][j] = 1.7;
98  }
99  }
100  doFit_ = false;
101  MinAmpl_ = 16;
102  dyn_pedestal = true;
103  };
template<class C>
EcalUncalibRecHitFixedAlphaBetaAlgo< C >::~EcalUncalibRecHitFixedAlphaBetaAlgo ( )
inlineoverride

Definition at line 105 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

105 {};

Member Function Documentation

template<class C >
void EcalUncalibRecHitFixedAlphaBetaAlgo< C >::InitFitParameters ( double *  samples,
int  max_sample 
)
private

Definition at line 229 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

References a, and b.

229  {
230  // in a first attempt just use the value of the maximum sample
231  fAmp_max_ = samples[max_sample];
232  fTim_max_ = max_sample;
233  fPed_max_ = 0;
234 
235  // amplitude too low for fit to converge
236  // timing set correctly is assumed
237  if (fAmp_max_ < MinAmpl_) {
238  fAmp_max_ = samples[5];
239  double sumA = samples[5] + samples[4] + samples[6];
240  if (sumA != 0) {
241  fTim_max_ = 5 + (samples[6] - samples[4]) / sumA;
242  } else {
243  fTim_max_ = -993;
244  } //-999+6
245  doFit_ = false;
246  }
247  // if timing very badly off, that just use max sample
248  else if (max_sample < 1 || max_sample > 7) {
249  doFit_ = false;
250  } else {
251  //y=a*(x-xM)^2+b*(x-xM)+c
252  doFit_ = true;
253  float a = float(samples[max_sample - 1] + samples[max_sample + 1] - 2 * samples[max_sample]) / 2.;
254  if (a == 0) {
255  doFit_ = false;
256  return;
257  }
258  float b = float(samples[max_sample + 1] - samples[max_sample - 1]) / 2.;
259 
260  fTim_max_ = max_sample - b / (2 * a);
261  fAmp_max_ = samples[max_sample] - b * b / (4 * a);
262  }
263 }
double b
Definition: hdecay.h:118
double a
Definition: hdecay.h:119
template<class C>
EcalUncalibratedRecHit EcalUncalibRecHitFixedAlphaBetaAlgo< C >::makeRecHit ( const C &  dataFrame,
const double *  pedestals,
const double *  gainRatios,
const EcalWeightSet::EcalWeightMatrix **  weights,
const EcalWeightSet::EcalChi2WeightMatrix **  chi2Matrix 
)
overridevirtual

Compute parameters.

Implements EcalUncalibRecHitRecAbsAlgo< C >.

Definition at line 118 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

References CastorSimpleRecAlgoImpl::isSaturated(), EcalUncalibratedRecHit::kSaturated, and EcalCondDBWriter_cfi::pedestal.

Referenced by EcalUncalibRecHitWorkerFixedAlphaBetaFit::run().

123  {
124  double chi2_(-1.);
125 
126  // double Gain12Equivalent[4]={0,1,2,12};
127  double frame[C::MAXSAMPLES]; // will contain the ADC values
128  double pedestal = 0; // carries pedestal for gain12 i.e. gainId==1
129 
130  int gainId0 = 1; // expected gainId at the beginning of dataFrame
131  int iGainSwitch = 0; // flags whether there's any gainId other than gainId0
132  int GainId = 0; // stores gainId at every sample
133  double maxsample(-1); // ADC value of maximal ped-subtracted sample
134  int imax(-1); // sample number of maximal ped-subtracted sample
135  bool external_pede = false;
136  bool isSaturated = false; // flag reporting whether gain0 has been found
137 
138  // Get time samples checking for Gain Switch and pedestals
139  if (pedestals) {
140  external_pede = true;
141  if (dyn_pedestal) {
142  pedestal = (double(dataFrame.sample(0).adc()) + double(dataFrame.sample(1).adc())) / 2.;
143  } else {
144  pedestal = pedestals[0];
145  }
146  for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
147  //create frame in adc gain 12 equivalent
148  GainId = dataFrame.sample(iSample).gainId();
149 
150  // FIX-ME: warning: the vector pedestal is supposed to have in the order G12, G6 and G1
151  // if GainId is zero treat it as 3 temporarily to protect against undefined
152  // frame will be set to ~max of gain1
153  if (GainId == 0) {
154  GainId = 3;
155  isSaturated = true;
156  }
157 
158  if (GainId != gainId0)
159  iGainSwitch = 1;
160 
161  if (GainId == gainId0) {
162  frame[iSample] = double(dataFrame.sample(iSample).adc()) - pedestal;
163  } else {
164  frame[iSample] = (double(dataFrame.sample(iSample).adc()) - pedestals[GainId - 1]) * gainRatios[GainId - 1];
165  }
166 
167  if (frame[iSample] > maxsample) {
168  maxsample = frame[iSample];
169  imax = iSample;
170  }
171  }
172  } else { // pedestal from pre-sample
173  external_pede = false;
174  pedestal = (double(dataFrame.sample(0).adc()) + double(dataFrame.sample(1).adc())) / 2.;
175 
176  for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
177  //create frame in adc gain 12 equivalent
178  GainId = dataFrame.sample(iSample).gainId();
179  //no gain switch forseen if there is no external pedestal
180  if (GainId == 0) {
181  GainId = 3;
182  isSaturated = true;
183  }
184 
185  frame[iSample] = double(dataFrame.sample(iSample).adc()) - pedestal;
186  // if gain has switched but no pedestals are available, no much good you can do...
187  if (GainId > gainId0)
188  iGainSwitch = 1;
189  if (frame[iSample] > maxsample) {
190  maxsample = frame[iSample];
191  imax = iSample;
192  }
193  }
194  }
195 
196  if ((iGainSwitch == 1 && external_pede == false) || // ... thus you return dummy rechit
197  imax == -1) { // protect against all frames being <-1
198  return EcalUncalibratedRecHit(dataFrame.id(), -1., -100., -1., -1.);
199  }
200 
201  InitFitParameters(frame, imax);
202  chi2_ = PerformAnalyticFit(frame, imax);
203  uint32_t flags = 0;
204  if (isSaturated)
206 
207  /* std::cout << "separate fits\nA: " << fAmp_max_ << ", ResidualPed: " << fPed_max_
208  <<", pedestal: "<<pedestal << ", tPeak " << fTim_max_ << std::endl;
209  */
210  return EcalUncalibratedRecHit(dataFrame.id(), fAmp_max_, pedestal + fPed_max_, fTim_max_ - 5, chi2_, flags);
211 }
void InitFitParameters(double *samples, int max_sample)
bool isSaturated(const Digi &digi, const int &maxADCvalue, int ifirst, int n)
float PerformAnalyticFit(double *samples, int max_sample)
template<class C >
float EcalUncalibRecHitFixedAlphaBetaAlgo< C >::PerformAnalyticFit ( double *  samples,
int  max_sample 
)
private

fit electronic function from simulation parameters fAlpha_ and fBeta_ are fixed and fit is providing the 3 following parameters the maximum amplitude ( fAmp_max_ ) the time of the maximum ( fTim_max_)

Loop on iterations

initialization inside iteration loop !

Then we loop on samples to be fitted

calculate function to be fitted

then calculate derivatives of function to be fitted

compute matrix elements DM1

compute delta

compute vector elements PROD

end of loop on samples

compute variations of parameters fAmp_max and fTim_max

end of loop on iterations

protection again diverging/unstable fit

results of the fit are calculated

Definition at line 266 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

References HLT_FULL_cff::chi2, EcalCondDB::db, CommonMethods::delta(), dt, funct::exp(), cms::cuda::func, mps_fire::i, cuy::ii, edm::isNotFinite(), LogDebug, funct::pow(), and PROD.

266  {
267  //int fValue_tim_max = max_sample;
272  //| the pedestal (fPed_max_)
273 
274  double chi2 = -1, db[3];
275 
276  //HepSymMatrix DM1(3) ; CLHEP::HepVector temp(3) ;
277 
278  int num_fit_min = (int)(max_sample - fNum_samp_bef_max_);
279  int num_fit_max = (int)(max_sample + fNum_samp_after_max_);
280 
281  if (num_fit_min < 0)
282  num_fit_min = 0;
283  //if (num_fit_max>=fNsamples-1) num_fit_max = fNsamples-2 ;
284  if (num_fit_max >= C::MAXSAMPLES) {
285  num_fit_max = C::MAXSAMPLES - 1;
286  }
287 
288  if (!doFit_) {
289  LogDebug("EcalUncalibRecHitFixedAlphaBetaAlgo") << "No fit performed. The amplitude of sample 5 will be used";
290  return -1;
291  }
292 
293  double func, delta;
294  double variation_func_max = 0.;
295  double variation_tim_max = 0.;
296  double variation_ped_max = 0.;
298  for (int iter = 0; iter < fNb_iter_; iter++) {
300  chi2 = 0.; //PROD.Zero() ; DM1.Zero() ;
301 
302  for (int i1 = 0; i1 < 3; i1++) {
303  temp_[i1] = 0;
304  for (int j1 = i1; j1 < 3; j1++) {
305  DM1_.fast(j1 + 1, i1 + 1) = 0;
306  }
307  }
308 
309  fAmp_max_ += variation_func_max;
310  fTim_max_ += variation_tim_max;
311  fPed_max_ += variation_ped_max;
312 
314  for (int i = num_fit_min; i <= num_fit_max; i++) {
315  //if(i>fsamp_edge_fit && i<num_fit_min) continue ; // remove front edge samples
317  func = pulseShapeFunction((double)i);
319  double dt = (double)i - fTim_max_;
320  if (dt > -alfabeta_) {
321  double dt_sur_beta = dt / fBeta_;
322  double variable = (double)1. + dt / alfabeta_;
323  double expo = exp(-dt_sur_beta);
324  double puissance = pow(variable, fAlpha_);
325 
326  db[0] = un_sur_sigma * puissance * expo;
327  db[1] = fAmp_max_ * db[0] * dt_sur_beta / (alfabeta_ * variable);
328  } else {
329  db[0] = 0.;
330  db[1] = 0.;
331  }
332  db[2] = un_sur_sigma;
334  for (int i1 = 0; i1 < 3; i1++) {
335  for (int j1 = i1; j1 < 3; j1++) {
336  //double & fast(int row, int col);
337  DM1_.fast(j1 + 1, i1 + 1) += db[i1] * db[j1];
338  }
339  }
341  delta = (samples[i] - func) * un_sur_sigma;
343  for (int ii = 0; ii < 3; ii++) {
344  temp_[ii] += delta * db[ii];
345  }
346  chi2 += delta * delta;
347  }
348 
349  int fail = 0;
350  DM1_.invert(fail);
351  if (fail != 0.) {
352  //just a guess from the value of the parameters in the previous interaction;
353  //printf("wH4PulseFitWithFunction =====> determinant error --> No Fit Provided !\n") ;
354  InitFitParameters(samples, max_sample);
355  return -101;
356  }
357  /* for(int i1=0 ; i1<3 ; i1++) { */
358  /* for(int j1=0 ; j1<3 ; j1++) { */
359  /* //double & fast(int row, int col); */
360  /* std::cout<<"inverted: "<<DM1[j1][i1]<<std::endl;;} */
361  /* } */
362  /* std::cout<<"vector temp: "<< temp[0]<<" "<<temp[1]<<" "<<temp[2]<<std::endl; */
364  CLHEP::HepVector PROD = DM1_ * temp_;
365  // std::cout<<"vector PROD: "<< PROD[0]<<" "<<PROD[1]<<" "<<PROD[2]<<std::endl;
366 
367  // Probably the fastest way to protect against
368  // +-inf value in the matrix DM1_ after inversion
369  // (which is nevertheless flagged as successfull...)
370  if (edm::isNotFinite(PROD[0])) {
371  InitFitParameters(samples, max_sample);
372  return -103;
373  }
374 
375  variation_func_max = PROD[0];
376  variation_tim_max = PROD[1];
377  variation_ped_max = PROD[2];
378  //chi2 = chi2/((double)nsamp_used - 3.) ;
379  }
380 
382  if (variation_func_max > 2000. || variation_func_max < -1000.) {
383  InitFitParameters(samples, max_sample);
384  return -102;
385  }
386 
388  fAmp_max_ += variation_func_max;
389  fTim_max_ += variation_tim_max;
390  fPed_max_ += variation_ped_max;
391 
392  // protection against unphysical results:
393  // ampli mismatched to MaxSample, ampli largely negative, time off preselected range
394  if (fAmp_max_ > 2 * samples[max_sample] || fAmp_max_ < -100 || fTim_max_ < 0 || 9 < fTim_max_) {
395  InitFitParameters(samples, max_sample);
396  return -104;
397  }
398 
399  //std::cout <<"chi2: "<<chi2<<" ampl: "<<fAmp_max_<<" time: "<<fTim_max_<<" pede: "<<fPed_max_<<std::endl;
400  return chi2;
401 }
void InitFitParameters(double *samples, int max_sample)
float dt
Definition: AMPTWrapper.h:136
constexpr bool isNotFinite(T x)
Definition: isFinite.h:9
tuple db
Definition: EcalCondDB.py:153
uint32_t T const *__restrict__ uint32_t const *__restrict__ int32_t int Histo::index_type cudaStream_t Func __host__ __device__ V int Func func
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
int ii
Definition: cuy.py:589
#define PROD(A, B)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
#define LogDebug(id)
template<class C >
double EcalUncalibRecHitFixedAlphaBetaAlgo< C >::pulseShapeFunction ( double  t)
private

Definition at line 214 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

References dt, funct::exp(), and funct::pow().

214  {
215  if (alfabeta_ <= 0)
216  return ((double)0.);
217  double dtsbeta, variable, puiss;
218  double dt = t - fTim_max_;
219  if (dt > -alfabeta_) {
220  dtsbeta = dt / fBeta_;
221  variable = 1. + dt / alfabeta_;
222  puiss = pow(variable, fAlpha_);
223  return fAmp_max_ * puiss * exp(-dtsbeta) + fPed_max_;
224  }
225  return fPed_max_;
226 }
float dt
Definition: AMPTWrapper.h:136
Exp< T >::type exp(const T &t)
Definition: Exp.h:22
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29
template<class C >
void EcalUncalibRecHitFixedAlphaBetaAlgo< C >::SetAlphaBeta ( double  alpha,
double  beta 
)
template<class C >
void EcalUncalibRecHitFixedAlphaBetaAlgo< C >::SetDynamicPedestal ( bool  dyn_pede)
template<class C >
void EcalUncalibRecHitFixedAlphaBetaAlgo< C >::SetMinAmpl ( double  ampl)

Member Data Documentation

template<class C>
double EcalUncalibRecHitFixedAlphaBetaAlgo< C >::alfabeta_
private

Definition at line 42 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

template<class C>
float EcalUncalibRecHitFixedAlphaBetaAlgo< C >::alpha_table_[36][1701]
private
template<class C>
float EcalUncalibRecHitFixedAlphaBetaAlgo< C >::beta_table_[36][1701]
private
template<class C>
CLHEP::HepSymMatrix EcalUncalibRecHitFixedAlphaBetaAlgo< C >::DM1_
private

Definition at line 59 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

template<class C>
bool EcalUncalibRecHitFixedAlphaBetaAlgo< C >::doFit_
private
template<class C>
bool EcalUncalibRecHitFixedAlphaBetaAlgo< C >::dyn_pedestal
private
template<class C>
double EcalUncalibRecHitFixedAlphaBetaAlgo< C >::fAlpha_
private

Definition at line 37 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

template<class C>
double EcalUncalibRecHitFixedAlphaBetaAlgo< C >::fAmp_max_
private

Definition at line 39 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

template<class C>
double EcalUncalibRecHitFixedAlphaBetaAlgo< C >::fBeta_
private

Definition at line 38 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

template<class C>
int EcalUncalibRecHitFixedAlphaBetaAlgo< C >::fNb_iter_
private
template<class C>
int EcalUncalibRecHitFixedAlphaBetaAlgo< C >::fNum_samp_after_max_
private
template<class C>
int EcalUncalibRecHitFixedAlphaBetaAlgo< C >::fNum_samp_bef_max_
private
template<class C>
double EcalUncalibRecHitFixedAlphaBetaAlgo< C >::fPed_max_
private

Definition at line 41 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

template<class C>
float EcalUncalibRecHitFixedAlphaBetaAlgo< C >::fSigma_ped
private
template<class C>
double EcalUncalibRecHitFixedAlphaBetaAlgo< C >::fTim_max_
private

Definition at line 40 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

template<class C>
double EcalUncalibRecHitFixedAlphaBetaAlgo< C >::MinAmpl_
private
template<class C>
CLHEP::HepVector EcalUncalibRecHitFixedAlphaBetaAlgo< C >::temp_
private

Definition at line 60 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

template<class C>
double EcalUncalibRecHitFixedAlphaBetaAlgo< C >::un_sur_sigma
private