CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros 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)
 
virtual EcalUncalibratedRecHit makeRecHit (const C &dataFrame, const double *pedestals, const double *gainRatios, const EcalWeightSet::EcalWeightMatrix **weights, const EcalWeightSet::EcalChi2WeightMatrix **chi2Matrix)
 Compute parameters. More...
 
void SetAlphaBeta (double alpha, double beta)
 
void SetDynamicPedestal (bool dyn_pede)
 
void SetMinAmpl (double ampl)
 
virtual ~EcalUncalibRecHitFixedAlphaBetaAlgo ()
 

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 32 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Constructor & Destructor Documentation

Definition at line 65 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

66  un_sur_sigma = 1./double(fSigma_ped) ;
67  for (int i=0;i<36;i++){
68  for(int j=0;j<1701;j++){
69  alpha_table_[i][j] = 1.2 ;
70  beta_table_[i][j] = 1.7 ;
71  }
72  }
73  doFit_ = false;
74  MinAmpl_ = 16;
75  dyn_pedestal = true;
76  }
int i
Definition: DBlmapReader.cc:9
int j
Definition: DBlmapReader.cc:9
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 77 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

77  :fAlpha_(0.),fBeta_(0.),fAmp_max_(-1.),fTim_max_(-1),fPed_max_(0),alfabeta_(0),DM1_(3),temp_(3){
78 
79  fNb_iter_ = n_iter ;
80  fNum_samp_bef_max_ = n_bef_max ;
81  fNum_samp_after_max_ = n_aft_max ;
82  fSigma_ped = sigma_ped;
83  un_sur_sigma = 1./double(fSigma_ped) ;
84  for (int i=0;i<36;i++){
85  for(int j=0;j<1701;j++){
86  alpha_table_[i][j] = 1.2 ;
87  beta_table_[i][j] = 1.7 ;
88  }
89  }
90  doFit_ = false;
91  MinAmpl_ = 16;
92  dyn_pedestal = true;
93  };
int i
Definition: DBlmapReader.cc:9
int j
Definition: DBlmapReader.cc:9
template<class C>
virtual EcalUncalibRecHitFixedAlphaBetaAlgo< C >::~EcalUncalibRecHitFixedAlphaBetaAlgo ( )
inlinevirtual

Definition at line 95 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

95 { };

Member Function Documentation

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

Definition at line 208 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

References a, and b.

208  {
209 
210  // in a first attempt just use the value of the maximum sample
211  fAmp_max_ = samples[max_sample];
212  fTim_max_ = max_sample;
213  fPed_max_ = 0;
214 
215  // amplitude too low for fit to converge
216  // timing set correctly is assumed
217  if(fAmp_max_ < MinAmpl_){
218  fAmp_max_ = samples[5];
219  double sumA = samples[5]+samples[4]+samples[6];
220  if(sumA != 0) { fTim_max_ = 5+(samples[6]-samples[4])/sumA; }
221  else{ fTim_max_ = -993; }//-999+6
222  doFit_ = false;
223  }
224  // if timing very badly off, that just use max sample
225  else if(max_sample <1 || max_sample > 7)
226  { doFit_ = false;}
227  else{
228  //y=a*(x-xM)^2+b*(x-xM)+c
229  doFit_ = true;
230  float a = float(samples[max_sample-1]+samples[max_sample+1]-2*samples[max_sample])/2.;
231  if(a==0){doFit_ =false; return;}
232  float b = float(samples[max_sample+1]-samples[max_sample-1])/2.;
233 
234  fTim_max_ = max_sample - b/(2*a);
235  fAmp_max_ = samples[max_sample] - b*b/(4*a);
236  }
237 
238 }
double b
Definition: hdecay.h:120
double a
Definition: hdecay.h:121
template<class C>
EcalUncalibratedRecHit EcalUncalibRecHitFixedAlphaBetaAlgo< C >::makeRecHit ( const C &  dataFrame,
const double *  pedestals,
const double *  gainRatios,
const EcalWeightSet::EcalWeightMatrix **  weights,
const EcalWeightSet::EcalChi2WeightMatrix **  chi2Matrix 
)
virtual

Compute parameters.

Implements EcalUncalibRecHitRecAbsAlgo< C >.

Definition at line 108 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

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

Referenced by EcalUncalibRecHitWorkerFixedAlphaBetaFit::run().

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

References EcalCondDB::db, delta, dt, create_public_lumi_plots::exp, cmsPerfPublish::fail(), i, cuy::ii, edm::isNotFinite(), LogDebug, funct::pow(), and PROD.

240  {
241 
242  //int fValue_tim_max = max_sample;
247  //| the pedestal (fPed_max_)
248 
249  double chi2=-1 , db[3] ;
250 
251 
252  //HepSymMatrix DM1(3) ; CLHEP::HepVector temp(3) ;
253 
254  int num_fit_min =(int)(max_sample - fNum_samp_bef_max_ ) ;
255  int num_fit_max =(int)(max_sample + fNum_samp_after_max_) ;
256 
257  if (num_fit_min<0) num_fit_min = 0 ;
258  //if (num_fit_max>=fNsamples-1) num_fit_max = fNsamples-2 ;
259  if (num_fit_max>= C::MAXSAMPLES ) {num_fit_max = C::MAXSAMPLES-1 ;}
260 
261  if(! doFit_ ) {
262  LogDebug("EcalUncalibRecHitFixedAlphaBetaAlgo")<<"No fit performed. The amplitude of sample 5 will be used";
263  return -1;
264  }
265 
266  double func,delta ;
267  double variation_func_max = 0. ;double variation_tim_max = 0. ; double variation_ped_max = 0. ;
269  for (int iter=0 ; iter < fNb_iter_ ; iter ++) {
271  chi2 = 0. ; //PROD.Zero() ; DM1.Zero() ;
272 
273  for(int i1=0 ; i1<3 ; i1++) {
274  temp_[i1]=0;
275  for(int j1=i1 ; j1<3 ; j1++) {
276  DM1_.fast(j1+1,i1+1) = 0; }
277  }
278 
279  fAmp_max_ += variation_func_max ;
280  fTim_max_ += variation_tim_max ;
281  fPed_max_ += variation_ped_max ;
282 
284  for( int i = num_fit_min ; i <= num_fit_max ; i++) {
285  //if(i>fsamp_edge_fit && i<num_fit_min) continue ; // remove front edge samples
287  func = pulseShapeFunction( (double)i ) ;
289  double dt =(double)i - fTim_max_ ;
290  if(dt > -alfabeta_) {
291  double dt_sur_beta = dt/fBeta_ ;
292  double variable = (double)1. + dt/alfabeta_ ;
293  double expo = exp(-dt_sur_beta) ;
294  double puissance = pow(variable,fAlpha_) ;
295 
296  db[0]=un_sur_sigma*puissance*expo ;
297  db[1]=fAmp_max_*db[0]*dt_sur_beta/(alfabeta_*variable) ;
298  }
299  else {
300  db[0]=0. ; db[1]=0. ;
301  }
302  db[2]=un_sur_sigma ;
304  for(int i1=0 ; i1<3 ; i1++) {
305  for(int j1=i1 ; j1<3 ; j1++) {
306  //double & fast(int row, int col);
307  DM1_.fast(j1+1,i1+1) += db[i1]*db[j1]; }
308  }
310  delta = (samples[i]-func)*un_sur_sigma ;
312  for(int ii=0 ; ii<3 ;ii++) {temp_[ii] += delta*db[ii] ;}
313  chi2 += delta * delta ;
314  }
315 
316  int fail=0 ;
317  DM1_.invert(fail) ;
318  if(fail != 0.) {
319  //just a guess from the value of the parameters in the previous interaction;
320  //printf("wH4PulseFitWithFunction =====> determinant error --> No Fit Provided !\n") ;
321  InitFitParameters(samples,max_sample);
322  return -101 ;
323  }
324 /* for(int i1=0 ; i1<3 ; i1++) { */
325 /* for(int j1=0 ; j1<3 ; j1++) { */
326 /* //double & fast(int row, int col); */
327 /* std::cout<<"inverted: "<<DM1[j1][i1]<<std::endl;;} */
328 /* } */
329 /* std::cout<<"vector temp: "<< temp[0]<<" "<<temp[1]<<" "<<temp[2]<<std::endl; */
331  CLHEP::HepVector PROD = DM1_*temp_ ;
332  // std::cout<<"vector PROD: "<< PROD[0]<<" "<<PROD[1]<<" "<<PROD[2]<<std::endl;
333 
334  // Probably the fastest way to protect against
335  // +-inf value in the matrix DM1_ after inversion
336  // (which is nevertheless flagged as successfull...)
337  if ( edm::isNotFinite( PROD[0] ) ) {
338  InitFitParameters(samples,max_sample);
339  return -103 ;
340  }
341 
342  variation_func_max = PROD[0] ;
343  variation_tim_max = PROD[1] ;
344  variation_ped_max = PROD[2] ;
345  //chi2 = chi2/((double)nsamp_used - 3.) ;
346  }
347 
348 
350  if( variation_func_max > 2000. || variation_func_max < -1000. ) {
351  InitFitParameters(samples,max_sample);
352  return -102;
353  }
354 
355 
357  fAmp_max_ += variation_func_max ;
358  fTim_max_ += variation_tim_max ;
359  fPed_max_ += variation_ped_max ;
360 
361 
362  // protection against unphysical results:
363  // ampli mismatched to MaxSample, ampli largely negative, time off preselected range
364  if( fAmp_max_>2*samples[max_sample] || fAmp_max_<-100 || fTim_max_<0 || 9<fTim_max_ ) {
365  InitFitParameters(samples,max_sample);
366  return -104;
367  }
368 
369 
370  //std::cout <<"chi2: "<<chi2<<" ampl: "<<fAmp_max_<<" time: "<<fTim_max_<<" pede: "<<fPed_max_<<std::endl;
371  return chi2;
372 }
#define LogDebug(id)
dbl * delta
Definition: mlp_gen.cc:36
void InitFitParameters(double *samples, int max_sample)
float dt
Definition: AMPTWrapper.h:126
int i
Definition: DBlmapReader.cc:9
tuple db
Definition: EcalCondDB.py:151
int ii
Definition: cuy.py:588
bool isNotFinite(T x)
Definition: isFinite.h:10
#define PROD(A, B)
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
template<class C >
double EcalUncalibRecHitFixedAlphaBetaAlgo< C >::pulseShapeFunction ( double  t)
private

Definition at line 195 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

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

195  {
196  if( alfabeta_ <= 0 ) return((double)0.);
197  double dtsbeta,variable,puiss;
198  double dt = t-fTim_max_ ;
199  if(dt > -alfabeta_) {
200  dtsbeta=dt/fBeta_ ;
201  variable=1.+dt/alfabeta_ ;
202  puiss=pow(variable,fAlpha_);
203  return fAmp_max_*puiss*exp(-dtsbeta)+fPed_max_ ;
204  }
205  return fPed_max_ ;
206 }
tuple t
Definition: tree.py:139
float dt
Definition: AMPTWrapper.h:126
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
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 44 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 62 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 39 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

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

Definition at line 41 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

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

Definition at line 40 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 43 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 42 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

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

Definition at line 62 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

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