Fast ALTRO fit (A.Pavlinov)
[u/mrichter/AliRoot.git] / PHOS / AliPHOSFastAltroFit.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id:$ */
17
18 /* History of cvs commits:
19  * $Log$
20  */
21
22 //______________________________________________________
23 // Author : Aleksei Pavlinov; IHEP, Protvino, Russia
24 // Feb 17, 2009
25 // Implementation of fit procedure from ALICE-INT-2008-026:
26 // "Time and amplitude reconstruction from sampling 
27 //  measurements of the PHOS signal profile"
28 //  M.Yu.Bogolyubsky and ..
29 //
30 //  Fit by function en*x*x*exp(-2.*x): x = (t-t0)/tau.
31 //  Tme main goal is fast estimation of amplitude.
32 //
33
34 // --- AliRoot header files ---
35 #include "AliPHOSFastAltroFit.h"
36
37 #include <math.h>
38
39 ClassImp(AliPHOSFastAltroFit)
40
41 //____________________________________________________________________________
42   AliPHOSFastAltroFit:: AliPHOSFastAltroFit() 
43 : TNamed(), 
44   fSig(0),fTau(0),fN(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0) 
45 {
46 }
47
48 //____________________________________________________________________________
49 AliPHOSFastAltroFit::AliPHOSFastAltroFit(const char* name, const char* title, const Double_t tau)
50   : TNamed(name, title), 
51   fSig(0),fTau(0),fN(2), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0) 
52 {
53   if(strlen(name)==0) SetName("FastAltroFit");
54 }
55
56 //____________________________________________________________________________
57 AliPHOSFastAltroFit::~AliPHOSFastAltroFit() 
58 {  
59 }
60
61 void AliPHOSFastAltroFit::FastFit(Int_t* t, Int_t* y, Int_t n, Double_t sig, Double_t tau, Double_t ped)
62 {
63   fSig = sig;
64   fTau = tau;
65
66   Double_t* tn = new Double_t[n]; 
67   Double_t* yn = new Double_t[n]; 
68   Int_t nn=0;
69
70   DeductPedestal(t,y,n,  tau,ped,  tn,yn,nn);
71   //  printf(" n %i : nn %i : ped %f \n", n, nn, ped);
72   // for(int i=0; i<nn; i++) 
73   // printf(" i %i : yn %7.2f : tn %7.2f \n", i, yn[i], tn[i]); 
74
75   FastFit(tn,yn,nn,sig,tau, fAmp,fAmpErr, fT0,fT0Err,fChi2);
76
77   if(fChi2> 0.0) fNDF = nn - 2;
78   else           fNDF = 0; 
79
80   delete [] tn;
81   delete [] yn;
82 }
83
84 void AliPHOSFastAltroFit::DeductPedestal(Int_t* t, Int_t* y, Int_t n, Double_t tau, Double_t ped, 
85   Double_t* tn, Double_t* yn, Int_t &nn)
86 {
87   Int_t ymax=0, nmax=0;
88   for(Int_t i=0; i<n; i++){
89     if(y[i]>ymax) {
90       ymax = y[i];
91       nmax = i;
92     }
93   }
94   Int_t i1 = nmax - Int_t(tau);
95   Int_t i2 = n; // No rejection of points from right - should correct for EMCAL
96   i1 = i1<0?0:i1;
97
98   nn = 0;
99   Double_t yd=0.0;
100   for(Int_t i=i1; i<i2; i++) {
101     yd = Double_t(y[i]) - ped;
102     if(yd>1.) {
103       yn[nn] = yd;
104       tn[nn] = Double_t(t[i]);
105       nn++;
106     }
107   }
108 }
109
110 void AliPHOSFastAltroFit::FastFit(Double_t* t, Double_t* y, Int_t n, Double_t sig, Double_t tau,
111 Double_t &amp, Double_t &eamp, Double_t &t0, Double_t &et0, Double_t &chi2)
112 {
113   // It is case of n=k=2 : fnn = x*x*exp(2 - 2*x)
114   // Input: 
115   //     n  - number of points 
116   //   t[n] - array of time bins
117   //   y[n] - array of amplitudes after pedestal subtractions;
118   //   sig   - error of amplitude measurement (one value for all channels)
119   //   tau   - filter time response (in timebin units)
120   // Output:
121   //       amp - amplitude at t0;
122   //        t0 - timeof max amplitude; 
123   static Double_t xx; // t/tau
124   static Double_t a, b, c;
125   static Double_t f02, f12, f22;    // functions
126   static Double_t f02d, f12d, f22d; // functions derivations
127   chi2 = -1;
128   if(n<=0) {
129     printf("<I> FastFit : n<=0 \n"); 
130     return;
131   }
132
133   a = b = c = 0.0;
134   for(Int_t i=0; i<n; i++){
135     xx  = t[i]/tau;
136     f02 = exp(-2.*xx);
137     f12 = xx*f02;
138     f22 = xx*f12;
139     // Derivations
140     f02d = -2.*f02;
141     f12d = f02 - 2.*f12;
142     f22d = 2.*(f12 - f22); 
143     //
144     a += f02d * y[i];
145     b -= 2.*f12d * y[i];
146     c += f22d * y[i];
147   }
148   Double_t t0_1=0.0, t0_2=0.0;
149   Double_t amp_1=0.0, amp_2=0.0, chi2_1=0.0, chi2_2=0.0;
150   if(QuadraticRoots(a,b,c, t0_1,t0_2)) {
151     t0_1 *= tau;
152     t0_2 *= tau;
153     Amplitude(t,y,n, sig, tau, t0_1, amp_1, chi2_1);
154     Amplitude(t,y,n, sig, tau, t0_2, amp_2, chi2_2);
155     if(0) {
156       printf(" t0_1 %f : t0_2 %f \n", t0_1, t0_2);
157       printf(" amp_1 %f : amp_2 %f \n", amp_1, amp_2);
158       printf(" chi2_1 %f : chi2_2 %f \n", chi2_1, chi2_2);
159     }
160     // t0 less on one tau with comparing with value from  "canonical equation"
161     amp  = amp_1;
162     t0   = t0_1;
163     chi2 = chi2_1; 
164     if(chi2_1 > chi2_2) {
165       amp  = amp_2;
166       t0   = t0_2; 
167       chi2 = chi2_2; 
168     }
169     CalculateParsErrors(t, y, n, sig, tau, amp, t0, eamp, et0);
170
171     // Fill1();
172     
173     // DrawFastFunction(amp, t0, fUtils->GetPedestalValue(), "1");
174     //    DrawFastFunction(amp_1, t0_1, fUtils->GetPedestalValue(), "1");
175     // DrawFastFunction(amp_2, t0_2, fUtils->GetPedestalValue(), "2");
176   } else {
177     chi2 = -1.; // bad fit - negative chi2
178   }
179 }
180
181 Bool_t AliPHOSFastAltroFit::QuadraticRoots(Double_t a, Double_t b, Double_t c, Double_t &x1, Double_t &x2)
182 {
183   // Resolve quadratic equations a*x**2 + b*x + c
184   //printf(" a %12.5e b %12.5e c %12.5e \n", a, b, c);
185   static Double_t dtmp = 0.0;
186   dtmp = b*b - 4.*a*c;
187
188   if(dtmp>=-1.0e-7 && dtmp<0.0) {
189     printf("QuadraticRoots : small negative square : dtmp %f ", dtmp);
190     dtmp = 0.0;
191   }
192   if(dtmp>=0.0) {
193     dtmp = sqrt(dtmp);
194     x1   = (-b + dtmp) / (2.*a);
195     x2   = (-b - dtmp) / (2.*a);
196
197     //    printf(" x1 %f : x2 %f \n", x1, x2);
198     return kTRUE;
199   } else {
200     printf("QuadraticRoots : negative square : dtmp %f ", dtmp);
201     return kFALSE;
202   }
203 }
204
205 void AliPHOSFastAltroFit::Amplitude(Double_t* t,Double_t* y,Int_t n, Double_t sig, Double_t tau, 
206 Double_t t0, Double_t &amp, Double_t &chi2)
207 {  
208   // Calculate parameters error too - Mar 24,09
209   // sig is independent from points
210   amp = 0.;
211   Double_t x=0.0, f=0.0, den=0.0, f02;
212   for(Int_t i=0; i<n; i++){
213     x    = (t[i] - t0)/tau;
214     f02  = exp(-2.*x);
215     f    = x*x*f02;     
216     amp += f*y[i];
217     den += f*f;
218   }
219   amp /= den;
220   //
221   // chi2 calculation
222   //
223   Double_t dy=0.0;
224   chi2=0.;
225   for(Int_t i=0; i<n; i++){
226     x    = (t[i] - t0)/tau;
227     f02  = exp(-2.*x);
228     f    = amp*x*x*f02;
229     dy   = y[i]-f;
230     chi2 += dy*dy;
231     //printf(" %i : y %f -> f %f : dy %f \n", i, y[i], f, dy); 
232   }
233   chi2 /= (sig*sig);
234 }
235
236 void AliPHOSFastAltroFit::CalculateParsErrors(Double_t* t, Double_t* y, Int_t n, Double_t sig, Double_t tau,
237 Double_t &amp, Double_t &t0, Double_t &eamp, Double_t &et0)
238 {
239   // fmax_nk = (n/k)**n*exp(-n) => n=k=2 => exp(-n) = exp(-2.)
240   static Double_t cc = exp(-2.);
241
242   Double_t sumf2=0.0, sumfd2=0.0, x, f02, f12, f22, f22d;
243
244   for(Int_t i=0; i<n; i++){
245     x    = (t[i] - t0)/tau;
246     f02  = amp*exp(-2.*x);
247     f12  = x*f02;
248     f22  = x*f12;
249     sumf2 += f22 * f22;
250     //
251     f22d = 2.*(f12 - f22); 
252     sumfd2 += f22d * f22d;
253   }
254
255   et0  = sig/amp/sqrt(sumfd2);
256   amp *= cc;
257   eamp = sig*cc/sqrt(sumf2); 
258 }