File indexing completed on 2023-03-17 11:18:41
0001 #ifndef RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecAnalFitAlgo_HH
0002 #define RecoLocalCalo_EcalRecAlgos_EcalUncalibRecHitRecAnalFitAlgo_HH
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "CLHEP/Matrix/Matrix.h"
0012 #include "CLHEP/Matrix/SymMatrix.h"
0013 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalUncalibRecHitRecAbsAlgo.h"
0014 #include <vector>
0015 #include <string>
0016
0017 #include "TROOT.h"
0018 #include "TGraph.h"
0019 #include "TF1.h"
0020
0021 #include "TMinuitMinimizer.h"
0022
0023 template <class C>
0024 class EcalUncalibRecHitRecAnalFitAlgo : public EcalUncalibRecHitRecAbsAlgo<C> {
0025 private:
0026 double pulseShapeFunction(double* var, double* par) {
0027 double x = var[0];
0028 double ampl = par[0];
0029 double tp = par[1];
0030 double alpha = par[2];
0031 double t0 = par[3];
0032
0033 double f = pow((x - t0) / tp, alpha) * exp(-alpha * (x - tp - t0) / tp);
0034 return ampl * f;
0035 };
0036
0037 double pedestalFunction(double* var, double* par) {
0038 double ped = par[0];
0039 return ped;
0040 };
0041
0042 public:
0043 EcalUncalibRecHitRecAnalFitAlgo() {
0044
0045
0046 TMinuitMinimizer::UseStaticMinuit(false);
0047 }
0048
0049
0050 ~EcalUncalibRecHitRecAnalFitAlgo<C>() override{};
0051
0052
0053 EcalUncalibratedRecHit makeRecHit(const C& dataFrame,
0054 const double* pedestals,
0055 const double* gainRatios,
0056 const EcalWeightSet::EcalWeightMatrix** weights,
0057 const EcalWeightSet::EcalChi2WeightMatrix** chi2Matrix) override {
0058 double amplitude_(-1.), pedestal_(-1.), jitter_(-1.), chi2_(-1.);
0059
0060
0061
0062 double frame[C::MAXSAMPLES];
0063
0064 int gainId0 = 1;
0065 int iGainSwitch = 0;
0066 double maxsample(-1);
0067 int imax(-1);
0068 bool isSaturated = false;
0069 uint32_t flag = 0;
0070 for (int iSample = 0; iSample < C::MAXSAMPLES; iSample++) {
0071 int gainId = dataFrame.sample(iSample).gainId();
0072 if (dataFrame.isSaturated()) {
0073 gainId = 3;
0074 isSaturated = true;
0075 }
0076
0077 if (gainId != gainId0)
0078 ++iGainSwitch;
0079 if (!iGainSwitch)
0080 frame[iSample] = double(dataFrame.sample(iSample).adc());
0081 else
0082 frame[iSample] =
0083 double(((double)(dataFrame.sample(iSample).adc()) - pedestals[gainId - 1]) * gainRatios[gainId - 1]);
0084
0085 if (frame[iSample] > maxsample) {
0086 maxsample = frame[iSample];
0087 imax = iSample;
0088 }
0089 }
0090
0091
0092
0093
0094
0095 double xarray[10] = {0., 1., 2., 3., 4., 5., 6., 7., 8., 9.};
0096 TGraph graph(10, xarray, frame);
0097
0098
0099 TF1 pulseShape = TF1("pulseShape",
0100 "[0]*pow((x - [3])/[1],[2])*exp(-[2]*(x - [1] - [3])/[1])",
0101 imax - 1.,
0102 imax + 3.,
0103 TF1::EAddToList::kNo);
0104 TF1 pedestal = TF1("pedestal", "[0]", 0., 2., TF1::EAddToList::kNo);
0105
0106
0107
0108 double FIT_A = (double)maxsample;
0109 pulseShape.SetParameter(0, FIT_A);
0110 pulseShape.SetParName(0, "Amplitude");
0111
0112 double FIT_Tp = (double)imax;
0113 pulseShape.SetParameter(1, FIT_Tp);
0114 pulseShape.SetParName(1, "t_{P}");
0115
0116 double FIT_ALFA = 1.5;
0117 pulseShape.SetParameter(2, FIT_ALFA);
0118 pulseShape.SetParName(2, "\\alpha");
0119
0120 double FIT_To = 3.;
0121 pulseShape.SetParameter(3, FIT_To);
0122 pulseShape.SetParName(3, "t_{0}");
0123
0124
0125 pedestal.SetParameter(0, frame[0]);
0126 pedestal.SetParName(0, "Pedestal");
0127
0128 int result = graph.Fit(&pulseShape, "QRMN SERIAL");
0129
0130 if (0 == result) {
0131 double amplitude_value = pulseShape.GetParameter(0);
0132
0133 graph.Fit(&pedestal, "QRLN SERIAL");
0134 double pedestal_value = pedestal.GetParameter(0);
0135
0136 if (!iGainSwitch)
0137 amplitude_ = amplitude_value - pedestal_value;
0138 else
0139 amplitude_ = amplitude_value;
0140
0141 pedestal_ = pedestal_value;
0142 jitter_ = pulseShape.GetParameter(3);
0143 chi2_ = 1.;
0144 if (isSaturated)
0145 flag = EcalUncalibratedRecHit::kSaturated;
0146
0147
0148
0149
0150
0151
0152 }
0153
0154 return EcalUncalibratedRecHit(dataFrame.id(), amplitude_, pedestal_, jitter_ - 6, chi2_, flag);
0155 }
0156 };
0157 #endif