00001 #ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecAnalFitAlgo_HH
00002 #define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecAnalFitAlgo_HH
00003
00014 #include "CLHEP/Matrix/Matrix.h"
00015 #include "CLHEP/Matrix/SymMatrix.h"
00016 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h"
00017 #include <vector>
00018 #include <string>
00019
00020 #include "TROOT.h"
00021 #include "TMinuit.h"
00022 #include "TGraph.h"
00023 #include "TF1.h"
00024
00025 template<class C> class EcalUncalibRecHitRecAnalFitAlgo : public EcalUncalibRecHitRecAbsAlgo<C>
00026 {
00027
00028
00029 private:
00030
00031 double pulseShapeFunction(double* var, double* par) {
00032 double x = var[0];
00033 double ampl = par[0];
00034 double tp = par[1];
00035 double alpha = par[2];
00036 double t0 = par[3];
00037
00038 double f = pow( (x-t0)/tp , alpha ) * exp( -alpha*(x-tp-t0)/tp );
00039 return ampl*f;
00040 };
00041
00042 double pedestalFunction(double* var, double* par) {
00043 double x = var[0];
00044 double ped = par[0];
00045 return ped;
00046 };
00047
00048 public:
00049
00050 virtual ~EcalUncalibRecHitRecAnalFitAlgo<C>() { };
00051
00052
00054 virtual EcalUncalibratedRecHit makeRecHit(const C& dataFrame, const double* pedestals,
00055 const double* gainRatios,
00056 const EcalWeightSet::EcalWeightMatrix** weights,
00057 const EcalWeightSet::EcalChi2WeightMatrix** chi2Matrix)
00058 {
00059 double amplitude_(-1.), pedestal_(-1.), jitter_(-1.), chi2_(-1.);
00060
00061
00062
00063 double frame[C::MAXSAMPLES];
00064
00065 int gainId0 = 1;
00066 int iGainSwitch = 0;
00067 double maxsample(-1);
00068 int imax(-1);
00069
00070 for(int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
00071 int gainId = dataFrame.sample(iSample).gainId();
00072 if ( gainId == 0 ) gainId = 3;
00073 if (gainId != gainId0) iGainSwitch++ ;
00074 if (!iGainSwitch)
00075 frame[iSample] = double(dataFrame.sample(iSample).adc());
00076 else
00077 frame[iSample] = double(((double)(dataFrame.sample(iSample).adc()) - pedestals[gainId-1]) * gainRatios[gainId-1]);
00078
00079 if( frame[iSample]>maxsample ) {
00080 maxsample= frame[iSample];
00081 imax=iSample;
00082 }
00083 }
00084
00085
00086
00087
00088
00089 double xarray[10]={0.,1.,2.,3.,4.,5.,6.,7.,8.,9.};
00090 TGraph graph(10,xarray,frame);
00091
00092
00093 TF1 pulseShape = TF1("pulseShape",
00094 "[0]*pow((x - [3])/[1],[2])*exp(-[2]*(x - [1] - [3])/[1])",
00095 imax-1.,imax+3.);
00096 TF1 pedestal = TF1("pedestal","[0]",0.,2.);
00097
00098
00099
00100 TF1 pluseAndPed = TF1("pulseAndPed","pedestal+pulseShape");
00101
00102
00103
00104 double FIT_A=(double)maxsample;
00105 pulseShape.SetParameter(0,FIT_A);
00106 pulseShape.SetParName(0,"Amplitude");
00107
00108 double FIT_Tp=(double)imax;
00109 pulseShape.SetParameter(1,FIT_Tp);
00110 pulseShape.SetParName(1,"t_{P}");
00111
00112 double FIT_ALFA=1.5;
00113 pulseShape.SetParameter(2,FIT_ALFA);
00114 pulseShape.SetParName(2,"\\alpha");
00115
00116 double FIT_To=3.;
00117 pulseShape.SetParameter(3,FIT_To);
00118 pulseShape.SetParName(3,"t_{0}");
00119
00120
00121 pedestal.SetParameter(0,frame[0]);
00122 pedestal.SetParName(0,"Pedestal");
00123
00124
00125
00126 graph.Fit("pulseShape","QRM");
00127
00128
00129 if ( std::string(gMinuit->fCstatu.Data()) == std::string("CONVERGED ") ) {
00130
00131 double amplitude_value=pulseShape.GetParameter(0);
00132
00133 graph.Fit("pedestal","QRL");
00134
00135 double pedestal_value=pedestal.GetParameter(0);
00136
00137 if (!iGainSwitch)
00138 amplitude_ = amplitude_value - pedestal_value;
00139 else
00140 amplitude_ = amplitude_value;
00141
00142 pedestal_ = pedestal_value;
00143 jitter_ = pulseShape.GetParameter(3);
00144 chi2_ = 1.;
00145
00146
00147
00148
00149
00150
00151
00152
00153 }
00154
00155 return EcalUncalibratedRecHit( dataFrame.id(), amplitude_, pedestal_, jitter_ - 6, chi2_);
00156 }
00157 };
00158 #endif