EcalRecAlgos/interface/EcalUncalibRecHitFixedAlphaBetaAlgo.h

0001 #ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitFixedAlphaBetaAlgo_HH
0002 #define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitFixedAlphaBetaAlgo_HH
0003
0004 /** \class EcalUncalibRecHitFixedAlphaBetaAlgo
0005   *  Template used to compute amplitude, pedestal, time jitter, chi2 of a pulse
0006   *  using an analytical function fit, with the pulse parameters alpha and beta fixed.
0007   *  It follows a fast fit algorithms devolped on test beam data by P. Jarry
0008   *  If the pedestal is not given, it is calculated from the first 2 pre-samples;
0009   *  FIXME: conversion gainID (1,2,3) with gain factor (12,6,1) is hardcoded here !
0010   *
0011   *  \author A.Ghezzi
0012   */
0013
0014 #include "CLHEP/Matrix/Vector.h"
0015 #include "CLHEP/Matrix/SymMatrix.h"
0016 //#include "CLHEP/Matrix/Matrix.h"
0017 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h"
0018
0019 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0020 #include "FWCore/Utilities/interface/isFinite.h"
0021
0022 #include <vector>
0023 #include <string>
0024
0025 //#include "TROOT.h"
0026 //#include "TMinuit.h"
0027 //#include "TGraph.h"
0028 //#include "TF1.h"
0029 //#include "TMatrixD.h"
0030 //#include "TVectorD.h"
0031
0032 template <class C>
0033 class EcalUncalibRecHitFixedAlphaBetaAlgo : public EcalUncalibRecHitRecAbsAlgo<C> {
0034 private:
0035   double MinAmpl_;
0036   bool dyn_pedestal;
0037   double fAlpha_;    //parameter of the shape
0038   double fBeta_;     //parameter of the shape
0039   double fAmp_max_;  // peak amplitude
0040   double fTim_max_;  // time of the peak (in 25ns units)
0041   double fPed_max_;  // pedestal value
0042   double alfabeta_;
0043
0044   int fNb_iter_;
0045   int fNum_samp_bef_max_;
0046   int fNum_samp_after_max_;
0047
0048   float fSigma_ped;
0049   double un_sur_sigma;
0050   //temporary solution for different alpha and beta
0051   float alpha_table_[36][1701];
0052   float beta_table_[36][1701];
0053
0054   bool doFit_;
0055
0056   double pulseShapeFunction(double t);
0057   float PerformAnalyticFit(double* samples, int max_sample);
0058   void InitFitParameters(double* samples, int max_sample);
0059   CLHEP::HepSymMatrix DM1_;
0060   CLHEP::HepVector temp_;
0061
0062 public:
0063   EcalUncalibRecHitFixedAlphaBetaAlgo()
0064       : fAlpha_(0.),
0065         fBeta_(0.),
0066         fAmp_max_(-1.),
0067         fTim_max_(-1),
0068         fPed_max_(0),
0069         alfabeta_(0),
0070         fNb_iter_(4),
0071         fNum_samp_bef_max_(1),
0072         fNum_samp_after_max_(3),
0073         fSigma_ped(1.1),
0074         DM1_(3),
0075         temp_(3) {
0076     un_sur_sigma = 1. / double(fSigma_ped);
0077     for (int i = 0; i < 36; i++) {
0078       for (int j = 0; j < 1701; j++) {
0079         alpha_table_[i][j] = 1.2;
0080         beta_table_[i][j] = 1.7;
0081       }
0082     }
0083     doFit_ = false;
0084     MinAmpl_ = 16;
0085     dyn_pedestal = true;
0086   }
0087   EcalUncalibRecHitFixedAlphaBetaAlgo(int n_iter, int n_bef_max = 1, int n_aft_max = 3, float sigma_ped = 1.1)
0088       : fAlpha_(0.), fBeta_(0.), fAmp_max_(-1.), fTim_max_(-1), fPed_max_(0), alfabeta_(0), DM1_(3), temp_(3) {
0089     fNb_iter_ = n_iter;
0090     fNum_samp_bef_max_ = n_bef_max;
0091     fNum_samp_after_max_ = n_aft_max;
0092     fSigma_ped = sigma_ped;
0093     un_sur_sigma = 1. / double(fSigma_ped);
0094     for (int i = 0; i < 36; i++) {
0095       for (int j = 0; j < 1701; j++) {
0096         alpha_table_[i][j] = 1.2;
0097         beta_table_[i][j] = 1.7;
0098       }
0099     }
0100     doFit_ = false;
0101     MinAmpl_ = 16;
0102     dyn_pedestal = true;
0103   };
0104
0105   ~EcalUncalibRecHitFixedAlphaBetaAlgo() override {}
0106   EcalUncalibratedRecHit makeRecHit(const C& dataFrame,
0107                                     const double* pedestals,
0108                                     const double* gainRatios,
0109                                     const EcalWeightSet::EcalWeightMatrix** weights,
0110                                     const EcalWeightSet::EcalChi2WeightMatrix** chi2Matrix) override;
0111   void SetAlphaBeta(double alpha, double beta);
0112   void SetMinAmpl(double ampl);
0113   void SetDynamicPedestal(bool dyn_pede);
0114 };
0115
0116 /// Compute parameters
0117 template <class C>
0118 EcalUncalibratedRecHit EcalUncalibRecHitFixedAlphaBetaAlgo<C>::makeRecHit(
0119     const C& dataFrame,
0120     const double* pedestals,
0121     const double* gainRatios,
0122     const EcalWeightSet::EcalWeightMatrix** weights,
0123     const EcalWeightSet::EcalChi2WeightMatrix** chi2Matrix) {
0124   double chi2_(-1.);
0125
0126   //  double Gain12Equivalent[4]={0,1,2,12};
0127   double frame[C::MAXSAMPLES];  // will contain the ADC values
0128   double pedestal = 0;          // carries pedestal for gain12 i.e. gainId==1
0129
0130   int gainId0 = 1;       // expected gainId at the beginning of dataFrame
0131   int iGainSwitch = 0;   // flags whether there's any gainId other than gainId0
0132   int GainId = 0;        // stores gainId at every sample
0133   double maxsample(-1);  // ADC value of maximal ped-subtracted sample
0134   int imax(-1);          // sample number of maximal ped-subtracted sample
0135   bool external_pede = false;
0136   bool isSaturated = false;  // flag reporting whether gain0 has been found
0137
0138   // Get time samples checking for Gain Switch and pedestals
0139   if (pedestals) {
0140     external_pede = true;
0141     if (dyn_pedestal) {
0142       pedestal = (double(dataFrame.sample(0).adc()) + double(dataFrame.sample(1).adc())) / 2.;
0143     } else {
0144       pedestal = pedestals[0];
0145     }
0146     for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
0147       //create frame in adc gain 12 equivalent
0148       GainId = dataFrame.sample(iSample).gainId();
0149
0150       // FIX-ME: warning: the vector pedestal is supposed to have in the order G12, G6 and G1
0151       // if GainId is zero treat it as 3 temporarily to protect against undefined
0152       // frame will be set to ~max of gain1
0153       if (GainId == 0) {
0154         GainId = 3;
0155         isSaturated = true;
0156       }
0157
0158       if (GainId != gainId0)
0159         iGainSwitch = 1;
0160
0161       if (GainId == gainId0) {
0162         frame[iSample] = double(dataFrame.sample(iSample).adc()) - pedestal;
0163       } else {
0164         frame[iSample] = (double(dataFrame.sample(iSample).adc()) - pedestals[GainId - 1]) * gainRatios[GainId - 1];
0165       }
0166
0167       if (frame[iSample] > maxsample) {
0168         maxsample = frame[iSample];
0169         imax = iSample;
0170       }
0171     }
0172   } else {  // pedestal from pre-sample
0173     external_pede = false;
0174     pedestal = (double(dataFrame.sample(0).adc()) + double(dataFrame.sample(1).adc())) / 2.;
0175
0176     for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
0177       //create frame in adc gain 12 equivalent
0178       GainId = dataFrame.sample(iSample).gainId();
0179       //no gain switch forseen if there is no external pedestal
0180       if (GainId == 0) {
0181         GainId = 3;
0182         isSaturated = true;
0183       }
0184
0185       frame[iSample] = double(dataFrame.sample(iSample).adc()) - pedestal;
0186       // if gain has switched but no pedestals are available, no much good you can do...
0187       if (GainId > gainId0)
0188         iGainSwitch = 1;
0189       if (frame[iSample] > maxsample) {
0190         maxsample = frame[iSample];
0191         imax = iSample;
0192       }
0193     }
0194   }
0195
0196   if ((iGainSwitch == 1 && external_pede == false) ||  // ... thus you return dummy rechit
0197       imax == -1) {                                    // protect against all frames being <-1
0198     return EcalUncalibratedRecHit(dataFrame.id(), -1., -100., -1., -1.);
0199   }
0200
0201   InitFitParameters(frame, imax);
0202   chi2_ = PerformAnalyticFit(frame, imax);
0203   uint32_t flags = 0;
0204   if (isSaturated)
0205     flags = EcalUncalibratedRecHit::kSaturated;
0206
0207   /*    std::cout << "separate fits\nA: " << fAmp_max_  << ", ResidualPed: " <<  fPed_max_
0208               <<", pedestal: "<<pedestal << ", tPeak " << fTim_max_ << std::endl;
0209   */
0210   return EcalUncalibratedRecHit(dataFrame.id(), fAmp_max_, pedestal + fPed_max_, fTim_max_ - 5, chi2_, flags);
0211 }
0212
0213 template <class C>
0214 double EcalUncalibRecHitFixedAlphaBetaAlgo<C>::pulseShapeFunction(double t) {
0215   if (alfabeta_ <= 0)
0216     return ((double)0.);
0217   double dtsbeta, variable, puiss;
0218   double dt = t - fTim_max_;
0219   if (dt > -alfabeta_) {
0220     dtsbeta = dt / fBeta_;
0221     variable = 1. + dt / alfabeta_;
0222     puiss = pow(variable, fAlpha_);
0223     return fAmp_max_ * puiss * exp(-dtsbeta) + fPed_max_;
0224   }
0225   return fPed_max_;
0226 }
0227
0228 template <class C>
0229 void EcalUncalibRecHitFixedAlphaBetaAlgo<C>::InitFitParameters(double* samples, int max_sample) {
0230   // in a first attempt just use the value of the maximum sample
0231   fAmp_max_ = samples[max_sample];
0232   fTim_max_ = max_sample;
0233   fPed_max_ = 0;
0234
0235   // amplitude too low for fit to converge
0236   // timing set correctly is assumed
0237   if (fAmp_max_ < MinAmpl_) {
0238     fAmp_max_ = samples[5];
0239     double sumA = samples[5] + samples[4] + samples[6];
0240     if (sumA != 0) {
0241       fTim_max_ = 5 + (samples[6] - samples[4]) / sumA;
0242     } else {
0243       fTim_max_ = -993;
0244     }  //-999+6
0245     doFit_ = false;
0246   }
0247   // if timing very badly off, that just use max sample
0248   else if (max_sample < 1 || max_sample > 7) {
0249     doFit_ = false;
0250   } else {
0251     //y=a*(x-xM)^2+b*(x-xM)+c
0252     doFit_ = true;
0253     float a = float(samples[max_sample - 1] + samples[max_sample + 1] - 2 * samples[max_sample]) / 2.;
0254     if (a == 0) {
0255       doFit_ = false;
0256       return;
0257     }
0258     float b = float(samples[max_sample + 1] - samples[max_sample - 1]) / 2.;
0259
0260     fTim_max_ = max_sample - b / (2 * a);
0261     fAmp_max_ = samples[max_sample] - b * b / (4 * a);
0262   }
0263 }
0264
0265 template <class C>
0266 float EcalUncalibRecHitFixedAlphaBetaAlgo<C>::PerformAnalyticFit(double* samples, int max_sample) {
0267   //int fValue_tim_max = max_sample;
0268   //! fit electronic function from simulation
0269   //! parameters fAlpha_ and fBeta_ are fixed and fit is providing the 3 following parameters
0270   //! the maximum amplitude ( fAmp_max_ )
0271   //! the time of the maximum  ( fTim_max_)
0272   //| the pedestal (fPed_max_)
0273
0274   double chi2 = -1, db[3];
0275
0276   //HepSymMatrix DM1(3) ; CLHEP::HepVector temp(3) ;
0277
0278   int num_fit_min = (int)(max_sample - fNum_samp_bef_max_);
0279   int num_fit_max = (int)(max_sample + fNum_samp_after_max_);
0280
0281   if (num_fit_min < 0)
0282     num_fit_min = 0;
0283   //if (num_fit_max>=fNsamples-1) num_fit_max = fNsamples-2 ;
0284   if (num_fit_max >= C::MAXSAMPLES) {
0285     num_fit_max = C::MAXSAMPLES - 1;
0286   }
0287
0288   if (!doFit_) {
0289     LogDebug("EcalUncalibRecHitFixedAlphaBetaAlgo") << "No fit performed. The amplitude of sample 5 will be used";
0290     return -1;
0291   }
0292
0293   double func, delta;
0294   double variation_func_max = 0.;
0295   double variation_tim_max = 0.;
0296   double variation_ped_max = 0.;
0297   //!          Loop on iterations
0298   for (int iter = 0; iter < fNb_iter_; iter++) {
0299     //!          initialization inside iteration loop !
0300     chi2 = 0.;  //PROD.Zero() ;  DM1.Zero() ;
0301
0302     for (int i1 = 0; i1 < 3; i1++) {
0303       temp_[i1] = 0;
0304       for (int j1 = i1; j1 < 3; j1++) {
0305         DM1_.fast(j1 + 1, i1 + 1) = 0;
0306       }
0307     }
0308
0309     fAmp_max_ += variation_func_max;
0310     fTim_max_ += variation_tim_max;
0311     fPed_max_ += variation_ped_max;
0312
0313     //! Then we loop on samples to be fitted
0314     for (int i = num_fit_min; i <= num_fit_max; i++) {
0315       //if(i>fsamp_edge_fit && i<num_fit_min) continue ; // remove front edge samples
0316       //! calculate function to be fitted
0317       func = pulseShapeFunction((double)i);
0318       //! then calculate derivatives of function to be fitted
0319       double dt = (double)i - fTim_max_;
0320       if (dt > -alfabeta_) {
0321         double dt_sur_beta = dt / fBeta_;
0322         double variable = (double)1. + dt / alfabeta_;
0323         double expo = exp(-dt_sur_beta);
0324         double puissance = pow(variable, fAlpha_);
0325
0326         db[0] = un_sur_sigma * puissance * expo;
0327         db[1] = fAmp_max_ * db[0] * dt_sur_beta / (alfabeta_ * variable);
0328       } else {
0329         db[0] = 0.;
0330         db[1] = 0.;
0331       }
0332       db[2] = un_sur_sigma;
0333       //! compute matrix elements DM1
0334       for (int i1 = 0; i1 < 3; i1++) {
0335         for (int j1 = i1; j1 < 3; j1++) {
0336           //double & fast(int row, int col);
0337           DM1_.fast(j1 + 1, i1 + 1) += db[i1] * db[j1];
0338         }
0339       }
0340       //! compute delta
0341       delta = (samples[i] - func) * un_sur_sigma;
0342       //! compute vector elements PROD
0343       for (int ii = 0; ii < 3; ii++) {
0344         temp_[ii] += delta * db[ii];
0345       }
0346       chi2 += delta * delta;
0347     }  //! end of loop on samples
0348
0349     int fail = 0;
0350     DM1_.invert(fail);
0351     if (fail != 0.) {
0352       //just a guess from the value of the parameters in the previous interaction;
0353       //printf("wH4PulseFitWithFunction =====> determinant error --> No Fit Provided !\n") ;
0354       InitFitParameters(samples, max_sample);
0355       return -101;
0356     }
0357     /*     for(int i1=0 ; i1<3 ; i1++) { */
0358     /*       for(int j1=0 ; j1<3 ; j1++) {  */
0359     /*  //double & fast(int row, int col); */
0360     /*  std::cout<<"inverted: "<<DM1[j1][i1]<<std::endl;;} */
0361     /*     } */
0362     /*     std::cout<<"vector temp: "<< temp[0]<<" "<<temp[1]<<" "<<temp[2]<<std::endl; */
0363     //! compute variations of parameters fAmp_max and fTim_max
0364     CLHEP::HepVector PROD = DM1_ * temp_;
0365     //    std::cout<<"vector PROD: "<< PROD[0]<<" "<<PROD[1]<<" "<<PROD[2]<<std::endl;
0366
0367     // Probably the fastest way to protect against
0368     // +-inf value in the matrix DM1_ after inversion
0369     // (which is nevertheless flagged as successfull...)
0370     if (edm::isNotFinite(PROD[0])) {
0371       InitFitParameters(samples, max_sample);
0372       return -103;
0373     }
0374
0375     variation_func_max = PROD[0];
0376     variation_tim_max = PROD[1];
0377     variation_ped_max = PROD[2];
0378     //chi2 = chi2/((double)nsamp_used - 3.) ;
0379   }  //!end of loop on iterations
0380
0381   //!   protection again diverging/unstable fit
0382   if (variation_func_max > 2000. || variation_func_max < -1000.) {
0383     InitFitParameters(samples, max_sample);
0384     return -102;
0385   }
0386
0387   //!      results of the fit are calculated
0388   fAmp_max_ += variation_func_max;
0389   fTim_max_ += variation_tim_max;
0390   fPed_max_ += variation_ped_max;
0391
0392   // protection against unphysical results:
0393   // ampli mismatched to MaxSample, ampli largely negative, time off preselected range
0394   if (fAmp_max_ > 2 * samples[max_sample] || fAmp_max_ < -100 || fTim_max_ < 0 || 9 < fTim_max_) {
0395     InitFitParameters(samples, max_sample);
0396     return -104;
0397   }
0398
0399   //std::cout <<"chi2: "<<chi2<<" ampl: "<<fAmp_max_<<" time: "<<fTim_max_<<" pede: "<<fPed_max_<<std::endl;
0400   return chi2;
0401 }
0402
0403 template <class C>
0404 void EcalUncalibRecHitFixedAlphaBetaAlgo<C>::SetAlphaBeta(double alpha, double beta) {
0405   fAlpha_ = alpha;
0406   fBeta_ = beta;
0407   alfabeta_ = fAlpha_ * fBeta_;
0408 }
0409
0410 template <class C>
0411 void EcalUncalibRecHitFixedAlphaBetaAlgo<C>::SetMinAmpl(double ampl) {
0412   MinAmpl_ = ampl;
0413 }
0414 template <class C>
0415 void EcalUncalibRecHitFixedAlphaBetaAlgo<C>::SetDynamicPedestal(bool p) {
0416   dyn_pedestal = p;
0417 }
0418 #endif