#include <EcalUncalibRecHitFixedAlphaBetaAlgo.h>

Inheritance diagram for EcalUncalibRecHitFixedAlphaBetaAlgo< 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.
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

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

template<class C>

EcalUncalibRecHitFixedAlphaBetaAlgo< C >::EcalUncalibRecHitFixedAlphaBetaAlgo ( ) [inline]

Definition at line 65 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

                                          :fAlpha_(0.),fBeta_(0.),fAmp_max_(-1.),fTim_max_(-1),fPed_max_(0),alfabeta_(0),fNb_iter_(4),fNum_samp_bef_max_(1),fNum_samp_after_max_(3),fSigma_ped(1.1),DM1_(3),temp_(3){
    un_sur_sigma = 1./double(fSigma_ped) ;
    for (int i=0;i<36;i++){
      for(int j=0;j<1701;j++){
        alpha_table_[i][j] = 1.2 ;
        beta_table_[i][j]  =  1.7 ;
      }
    }
    doFit_ = false;
    MinAmpl_ = 16;
    dyn_pedestal = true;
  }

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.

                                                                                                               :fAlpha_(0.),fBeta_(0.),fAmp_max_(-1.),fTim_max_(-1),fPed_max_(0),alfabeta_(0),DM1_(3),temp_(3){
    
   fNb_iter_ = n_iter ;
   fNum_samp_bef_max_   = n_bef_max ;
   fNum_samp_after_max_ = n_aft_max ;
   fSigma_ped = sigma_ped;
   un_sur_sigma = 1./double(fSigma_ped) ;
   for (int i=0;i<36;i++){
     for(int j=0;j<1701;j++){
       alpha_table_[i][j] = 1.2 ;
       beta_table_[i][j]  =  1.7 ;
     }
   }
   doFit_ = false;
   MinAmpl_ = 16;
   dyn_pedestal = true;
  };

template<class C>

virtual EcalUncalibRecHitFixedAlphaBetaAlgo< C >::~EcalUncalibRecHitFixedAlphaBetaAlgo ( ) [inline, virtual]

Definition at line 95 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

{ };

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.

                                                                                                                {

  // in a first attempt just use the value of the maximum sample 
  fAmp_max_ = samples[max_sample];
  fTim_max_ = max_sample;
  fPed_max_ = 0;

  // amplitude too low for fit to converge
  // timing set correctly is assumed
  if(fAmp_max_ <  MinAmpl_){
    fAmp_max_      = samples[5];
    double sumA    = samples[5]+samples[4]+samples[6];
    if(sumA != 0) { fTim_max_ = 5+(samples[6]-samples[4])/sumA; }
    else{ fTim_max_ = -993; }//-999+6
    doFit_ = false;
  }
  // if timing very badly off, that just use max sample
  else if(max_sample <1 || max_sample > 7)
    {    doFit_ = false;}
  else{
    //y=a*(x-xM)^2+b*(x-xM)+c
    doFit_ = true;
    float a = float(samples[max_sample-1]+samples[max_sample+1]-2*samples[max_sample])/2.;
    if(a==0){doFit_ =false; return;}
    float b = float(samples[max_sample+1]-samples[max_sample-1])/2.;
  
    fTim_max_ = max_sample - b/(2*a);
    fAmp_max_ =  samples[max_sample] - b*b/(4*a); 
  }  
  
}

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(), and EcalUncalibratedRecHit::kSaturated.

Referenced by EcalUncalibRecHitWorkerFixedAlphaBetaFit::run().

                                                                                                                                                  {
  double chi2_(-1.);
  
  //  double Gain12Equivalent[4]={0,1,2,12};
  double frame[C::MAXSAMPLES];// will contain the ADC values
  double pedestal =0;     // carries pedestal for gain12 i.e. gainId==1
  
  int gainId0 = 1;        // expected gainId at the beginning of dataFrame 
  int iGainSwitch = 0;    // flags whether there's any gainId other than gainId0
  int GainId= 0;          // stores gainId at every sample
  double maxsample(-1);   // ADC value of maximal ped-subtracted sample
  int imax(-1);           // sample number of maximal ped-subtracted sample
  bool external_pede = false;
  bool isSaturated = 0;   // flag reporting whether gain0 has been found

  // Get time samples checking for Gain Switch and pedestals
  if(pedestals){
    external_pede = true;
    if(dyn_pedestal) { pedestal = (double(dataFrame.sample(0).adc()) + double(dataFrame.sample(1).adc()))/2.;}
    else{ pedestal  = pedestals[0];}
    for(int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
        //create frame in adc gain 12 equivalent
        GainId = dataFrame.sample(iSample).gainId();

        // FIX-ME: warning: the vector pedestal is supposed to have in the order G12, G6 and G1
        // if GainId is zero treat it as 3 temporarily to protect against undefined
        // frame will be set to ~max of gain1
        if ( GainId == 0 )
          { 
            GainId = 3;
            isSaturated = 1;
          }

        if (GainId != gainId0) iGainSwitch = 1;

        if(GainId==gainId0){frame[iSample] = double(dataFrame.sample(iSample).adc())-pedestal ;}
        else {frame[iSample] = (double(dataFrame.sample(iSample).adc())-pedestals[GainId-1])*gainRatios[GainId-1];}

        if( frame[iSample]>maxsample ) {
          maxsample = frame[iSample];
          imax = iSample;
        }
      }
  }
  else {// pedestal from pre-sample
    external_pede = false;
    pedestal = (double(dataFrame.sample(0).adc()) + double(dataFrame.sample(1).adc()))/2.;

    for(int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
      //create frame in adc gain 12 equivalent
      GainId = dataFrame.sample(iSample).gainId();
      //no gain switch forseen if there is no external pedestal
      if ( GainId == 0 ) 
        {
          GainId = 3;
          isSaturated = 1;
        }

      frame[iSample] = double(dataFrame.sample(iSample).adc())-pedestal ;
      // if gain has switched but no pedestals are available, no much good you can do...
      if (GainId > gainId0) iGainSwitch = 1;
      if( frame[iSample]>maxsample ) {
        maxsample = frame[iSample];
        imax = iSample;
      }
    } 
  }

  if( (iGainSwitch==1 && external_pede==false) ||  // ... thus you return dummy rechit
      imax ==-1 ){                                 // protect against all frames being <-1
    return EcalUncalibratedRecHit( dataFrame.id(), -1., -100., -1. , -1.);
  }
  
  InitFitParameters(frame, imax);
  chi2_ = PerformAnalyticFit(frame,imax);
  uint32_t flags = 0;
  if (isSaturated) flags = EcalUncalibratedRecHit::kSaturated;

  /*    std::cout << "separate fits\nA: " << fAmp_max_  << ", ResidualPed: " <<  fPed_max_
              <<", pedestal: "<<pedestal << ", tPeak " << fTim_max_ << std::endl;
  */
  return EcalUncalibratedRecHit( dataFrame.id(),fAmp_max_, pedestal+fPed_max_, fTim_max_ - 5 , chi2_, flags );
}

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.

                                                                                                                 {
  
  //int fValue_tim_max = max_sample;  
  //| the pedestal (fPed_max_)  
  
  double chi2=-1 , db[3] ;
  

  //HepSymMatrix DM1(3) ; CLHEP::HepVector temp(3) ;

  int num_fit_min =(int)(max_sample - fNum_samp_bef_max_ ) ;
  int num_fit_max =(int)(max_sample + fNum_samp_after_max_) ;

  if (num_fit_min<0) num_fit_min = 0 ; 
  //if (num_fit_max>=fNsamples-1) num_fit_max = fNsamples-2 ;
  if (num_fit_max>= C::MAXSAMPLES ) {num_fit_max = C::MAXSAMPLES-1 ;}

  if(! doFit_ ) {
    LogDebug("EcalUncalibRecHitFixedAlphaBetaAlgo")<<"No fit performed. The amplitude of sample 5 will be used";
    return -1;
  }

  double func,delta ;
  double variation_func_max = 0. ;double variation_tim_max = 0. ; double variation_ped_max = 0. ;
  for (int iter=0 ; iter < fNb_iter_ ; iter ++) {
    chi2 = 0. ; //PROD.Zero() ;  DM1.Zero() ;

     for(int i1=0 ; i1<3 ; i1++) {
       temp_[i1]=0;
        for(int j1=i1 ; j1<3 ; j1++) { 
          DM1_.fast(j1+1,i1+1) = 0; }
      }

    fAmp_max_ += variation_func_max ;
    fTim_max_ += variation_tim_max ;
    fPed_max_ += variation_ped_max ;
    
    for( int i = num_fit_min ; i <= num_fit_max ; i++) {
      //if(i>fsamp_edge_fit && i<num_fit_min) continue ; // remove front edge samples
      func = pulseShapeFunction( (double)i  ) ;
      double dt =(double)i - fTim_max_ ;
      if(dt > -alfabeta_)  {      
        double dt_sur_beta = dt/fBeta_ ;
        double variable = (double)1. + dt/alfabeta_ ;
        double expo = exp(-dt_sur_beta) ;        
        double puissance = pow(variable,fAlpha_) ;
        
        db[0]=un_sur_sigma*puissance*expo ;
        db[1]=fAmp_max_*db[0]*dt_sur_beta/(alfabeta_*variable) ;
      }
      else {
        db[0]=0. ; db[1]=0. ; 
      }
      db[2]=un_sur_sigma ;
      for(int i1=0 ; i1<3 ; i1++) {
        for(int j1=i1 ; j1<3 ; j1++) { 
          //double & fast(int row, int col);
          DM1_.fast(j1+1,i1+1) += db[i1]*db[j1]; }
      }
      delta = (samples[i]-func)*un_sur_sigma ;
      for(int ii=0 ; ii<3 ;ii++) {temp_[ii] += delta*db[ii] ;}
      chi2 += delta * delta ;
    }
    
    int fail=0 ;
    DM1_.invert(fail) ;
      if(fail != 0.) {
      //just a guess from the value of the parameters in the previous interaction;
      //printf("wH4PulseFitWithFunction =====> determinant error --> No Fit Provided !\n") ;
      InitFitParameters(samples,max_sample);
      return -101 ;
    }
/*     for(int i1=0 ; i1<3 ; i1++) { */
/*       for(int j1=0 ; j1<3 ; j1++) {  */
/*      //double & fast(int row, int col); */
/*      std::cout<<"inverted: "<<DM1[j1][i1]<<std::endl;;} */
/*     } */
/*     std::cout<<"vector temp: "<< temp[0]<<" "<<temp[1]<<" "<<temp[2]<<std::endl; */
    CLHEP::HepVector PROD = DM1_*temp_ ;
    //    std::cout<<"vector PROD: "<< PROD[0]<<" "<<PROD[1]<<" "<<PROD[2]<<std::endl;

    // Probably the fastest way to protect against
    // +-inf value in the matrix DM1_ after inversion
    // (which is nevertheless flagged as successfull...)
    if ( edm::isNotFinite( PROD[0] ) ) {
            InitFitParameters(samples,max_sample);
            return -103 ;
    }

    variation_func_max = PROD[0] ;
    variation_tim_max = PROD[1] ;
    variation_ped_max = PROD[2] ;
    //chi2 = chi2/((double)nsamp_used - 3.) ;
  }
  

  if( variation_func_max > 2000. || variation_func_max < -1000. ) {
    InitFitParameters(samples,max_sample);
    return -102;
  }
  

  fAmp_max_ += variation_func_max ;
  fTim_max_ += variation_tim_max ;
  fPed_max_ += variation_ped_max ;


  // protection against unphysical results:
  // ampli mismatched to MaxSample, ampli largely negative, time off preselected range
  if( fAmp_max_>2*samples[max_sample]  ||  fAmp_max_<-100 ||  fTim_max_<0  ||  9<fTim_max_ ) {
    InitFitParameters(samples,max_sample);
    return -104;
  }
  

  //std::cout <<"chi2: "<<chi2<<" ampl: "<<fAmp_max_<<" time: "<<fTim_max_<<" pede: "<<fPed_max_<<std::endl;
  return chi2;
}

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().

                                                                                           {
  if( alfabeta_ <= 0 ) return((double)0.);
  double dtsbeta,variable,puiss;
  double dt = t-fTim_max_ ;
  if(dt > -alfabeta_)  {
    dtsbeta=dt/fBeta_ ;
    variable=1.+dt/alfabeta_ ;
    puiss=pow(variable,fAlpha_);
    return fAmp_max_*puiss*exp(-dtsbeta)+fPed_max_ ;
  }
  return  fPed_max_ ;
}

template<class C >

void EcalUncalibRecHitFixedAlphaBetaAlgo< C >::SetAlphaBeta	(	double	alpha,
		double	beta
	)

Definition at line 374 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

References alpha, and beta.

Referenced by EcalUncalibRecHitWorkerFixedAlphaBetaFit::run().

                                                                                                     {
  fAlpha_ = alpha;
  fBeta_=  beta;
  alfabeta_ = fAlpha_*fBeta_;
}

template<class C >

void EcalUncalibRecHitFixedAlphaBetaAlgo< C >::SetDynamicPedestal ( bool dyn_pede )

Definition at line 383 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

References AlCaHLTBitMon_ParallelJobs::p.

Referenced by EcalUncalibRecHitWorkerFixedAlphaBetaFit::EcalUncalibRecHitWorkerFixedAlphaBetaFit().

                                                                                        {
  dyn_pedestal = p;
}

template<class C >

void EcalUncalibRecHitFixedAlphaBetaAlgo< C >::SetMinAmpl ( double ampl )

Definition at line 380 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Referenced by EcalUncalibRecHitWorkerFixedAlphaBetaFit::EcalUncalibRecHitWorkerFixedAlphaBetaFit().

                                                                                     {
  MinAmpl_ = 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]

Definition at line 54 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Referenced by EcalUncalibRecHitFixedAlphaBetaAlgo< EBDataFrame >::EcalUncalibRecHitFixedAlphaBetaAlgo().

template<class C>

float EcalUncalibRecHitFixedAlphaBetaAlgo< C >::beta_table_[36][1701] [private]

Definition at line 55 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Referenced by EcalUncalibRecHitFixedAlphaBetaAlgo< EBDataFrame >::EcalUncalibRecHitFixedAlphaBetaAlgo().

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]

Definition at line 57 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Referenced by EcalUncalibRecHitFixedAlphaBetaAlgo< EBDataFrame >::EcalUncalibRecHitFixedAlphaBetaAlgo().

template<class C>

bool EcalUncalibRecHitFixedAlphaBetaAlgo< C >::dyn_pedestal [private]

Definition at line 38 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Referenced by EcalUncalibRecHitFixedAlphaBetaAlgo< EBDataFrame >::EcalUncalibRecHitFixedAlphaBetaAlgo().

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]

Definition at line 47 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Referenced by EcalUncalibRecHitFixedAlphaBetaAlgo< EBDataFrame >::EcalUncalibRecHitFixedAlphaBetaAlgo().

template<class C>

int EcalUncalibRecHitFixedAlphaBetaAlgo< C >::fNum_samp_after_max_ [private]

Definition at line 49 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Referenced by EcalUncalibRecHitFixedAlphaBetaAlgo< EBDataFrame >::EcalUncalibRecHitFixedAlphaBetaAlgo().

template<class C>

int EcalUncalibRecHitFixedAlphaBetaAlgo< C >::fNum_samp_bef_max_ [private]

Definition at line 48 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Referenced by EcalUncalibRecHitFixedAlphaBetaAlgo< EBDataFrame >::EcalUncalibRecHitFixedAlphaBetaAlgo().

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]

Definition at line 51 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Referenced by EcalUncalibRecHitFixedAlphaBetaAlgo< EBDataFrame >::EcalUncalibRecHitFixedAlphaBetaAlgo().

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]

Definition at line 37 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Referenced by EcalUncalibRecHitFixedAlphaBetaAlgo< EBDataFrame >::EcalUncalibRecHitFixedAlphaBetaAlgo().

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]

Definition at line 52 of file EcalUncalibRecHitFixedAlphaBetaAlgo.h.

Referenced by EcalUncalibRecHitFixedAlphaBetaAlgo< EBDataFrame >::EcalUncalibRecHitFixedAlphaBetaAlgo().

EcalUncalibRecHitFixedAlphaBetaAlgo< C > Class Template Reference

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

template<class C> class EcalUncalibRecHitFixedAlphaBetaAlgo< C >

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation

template<class C>
class EcalUncalibRecHitFixedAlphaBetaAlgo< C >