bug fixed
[u/mrichter/AliRoot.git] / RAW / AliCaloFastAltroFitv0.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2010 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 //  The main goal is fast estimation of amplitude and t0.
32 //
33
34 // --- AliRoot header files ---
35 #include "AliCaloFastAltroFitv0.h"
36
37 #include <TH1.h>
38 #include <TF1.h>
39 #include <TCanvas.h>
40 #include <TGraphErrors.h>
41 #include <TMath.h>
42
43 #include <math.h>
44
45 ClassImp(AliCaloFastAltroFitv0)
46
47 //____________________________________________________________________________
48   AliCaloFastAltroFitv0::AliCaloFastAltroFitv0() 
49 : TNamed(), 
50   fSig(0),fTau(0),fN(0),fPed(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0)
51 ,fNoFit(0),fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
52 {
53 }
54
55 //____________________________________________________________________________
56 AliCaloFastAltroFitv0::AliCaloFastAltroFitv0(const char* name, const char* title, 
57   const Double_t sig, const Double_t tau, const Double_t n)
58   : TNamed(name, title), 
59   fSig(sig),fTau(tau),fN(n),fPed(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0) 
60  ,fNoFit(0),fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
61 {
62   if(strlen(name)==0) SetName("CaloFastAltroFitv0");
63 }
64
65 //____________________________________________________________________________
66 AliCaloFastAltroFitv0::AliCaloFastAltroFitv0(const AliCaloFastAltroFitv0 &obj)
67   : TNamed(obj), 
68   fSig(0),fTau(0),fN(2.),fPed(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0) 
69  ,fNoFit(0),fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
70 {
71 }
72
73 //____________________________________________________________________________
74 AliCaloFastAltroFitv0::~AliCaloFastAltroFitv0() 
75 {  
76   if(fTfit) delete [] fTfit;
77   if(fAmpfit) delete [] fAmpfit;
78 }
79
80 //____________________________________________________________________________
81 AliCaloFastAltroFitv0& AliCaloFastAltroFitv0::operator= (const AliCaloFastAltroFitv0 &/*obj*/)
82 {
83   // Not implemented yet
84   return (*this);
85 }
86
87 void AliCaloFastAltroFitv0::FastFit(Int_t* t, Int_t* y, Int_t nPoints, Double_t sig, Double_t tau, 
88                                   Double_t n, Double_t ped, Double_t tMax)
89 {
90   // n 2 here and unused
91   n = 2.;
92   Reset();
93
94   fSig = sig;
95   fTau = tau;
96   fPed = ped;
97
98   Int_t ii=0;
99   CutRightPart(t,y,nPoints,tMax, ii);
100   nPoints = ii;
101
102   fNfit   = 0;
103   fTfit   = new Double_t[nPoints]; 
104   fAmpfit = new Double_t[nPoints]; 
105
106   
107   DeductPedestal(t,y,nPoints,  tau,ped,  fTfit,fAmpfit,fNfit);
108   //  printf(" n %i : fNfit %i : ped %f \n", n, fNfit, ped);
109   // for(int i=0; i<fNfit; i++) 
110   // printf(" i %i : fAmpfit %7.2f : fTfit %7.2f \n", i, fAmpfit[i], fTfit[i]); 
111
112   if(fNfit>=2) {
113     FastFit(fTfit,fAmpfit,fNfit,sig,tau, fAmp,fAmpErr, fT0,fT0Err,fChi2);
114
115     if(fChi2> 0.0) {
116       fNDF = fNfit - 2;
117     } else {
118       fNDF = 0;
119       fNoFit++;
120     }
121   } else if(fNfit==1){
122     Reset(); // What to do here => fT0 = fTfit[0]; fAmp = fAmpFit[0] ??
123   } else {
124     Reset();
125   }
126 }
127
128 //____________________________________________________________________________
129 void AliCaloFastAltroFitv0::FastFit(TH1F* h, Double_t sig, Double_t tau, Double_t n,
130 Double_t ped, Double_t tMax)
131 {
132   // Service method for convinience only
133   // h - hist with altro response, could have empty bin
134   // and center of bin could be different from bin number  
135   Reset();
136
137   if(h==0) return;
138   Int_t nPoints = h->GetNbinsX();
139   if(nPoints<2) return; // Sep 07, 09
140
141   Int_t* t = new Int_t[nPoints];
142   Int_t* y = new Int_t[nPoints];
143
144   Int_t nPositive=0;
145   for(Int_t i=0; i<nPoints; i++) {
146     if(h->GetBinContent(i+1) > 0.0){ // Get only positive
147       y[nPositive] = Int_t(h->GetBinContent(i+1));
148       t[nPositive] = Int_t(h->GetBinCenter(i+1)+0.0001);
149       nPositive++;
150     }
151   }
152
153   if(nPositive >= 2) {
154     FastFit(t,y,nPoints, sig,tau,n,ped, tMax);
155   }
156   if(fChi2<=0.0) fNoFit++;
157
158   delete [] t;
159   delete [] y;
160 }
161
162 void AliCaloFastAltroFitv0::Reset()
163 {
164   // Reset variables
165   fSig  = fTau = 0.0;
166   fAmp  = fAmpErr = fT0 = fT0Err = 0.0;
167   fChi2 = -1.;
168   fNDF  = fNfit = 0;
169
170   if(fTfit)   delete [] fTfit;
171   if(fAmpfit) delete [] fAmpfit;
172   fTfit = fAmpfit = 0;
173 }
174
175
176 void AliCaloFastAltroFitv0::GetFitResult(Double_t &amp,Double_t &eamp,Double_t &t0,Double_t &et0, 
177 Double_t &chi2, Int_t &ndf) const
178 {
179   // Return results of fitting procedure
180   amp  = fAmp;
181   eamp = fAmpErr;
182   t0   = fT0;
183   et0  = fT0Err;
184   chi2 = fChi2;
185   ndf  = fNDF;
186 }
187
188 void AliCaloFastAltroFitv0::GetFittedPoints(Int_t &nfit, Double_t* ar[2]) const
189 {
190   nfit  = fNfit;
191   ar[0] = fTfit;
192   ar[1] = fAmpfit;
193 }
194 //
195 // static functions
196 // 
197 void  AliCaloFastAltroFitv0::CutRightPart(Int_t *t,Int_t *y,Int_t nPoints,Double_t tMax, Int_t &ii)
198 {
199   // Cut right part of altro sample : static function
200   Int_t tt=0;
201   for(Int_t i=0; i<nPoints; i++) {
202     tt = t[i];
203     if(tMax && tt <= Int_t(tMax)) {
204       t[ii] = tt;
205       y[ii] = y[i];
206       ii++;
207     }
208   }
209   if(0) printf(" ii %i -> ii %i : tMax %7.2f \n", nPoints, ii, tMax);
210 }
211
212 void AliCaloFastAltroFitv0::DeductPedestal(Int_t* t, Int_t* y, Int_t nPoints, Double_t tau, Double_t ped, 
213   Double_t* tn, Double_t* yn, Int_t &nPointsOut)
214 {
215   // Pedestal deduction if ped is positive : static function
216   // Discard part od samle if it is not compact.
217   static Double_t yMinUnderPed=2.; // should be tune
218   Int_t ymax=0, nmax=0;
219   for(Int_t i=0; i<nPoints; i++){
220     if(y[i]>ymax) {
221       ymax = y[i];
222       nmax = i;
223     }
224   }
225   Int_t i1 = nmax - Int_t(tau);
226   //i1 = 0;
227   i1 = i1<0?0:i1;
228   Int_t i2 = nPoints;
229
230   nPointsOut = 0;
231   Double_t yd=0.0, tdiff=0.0;;
232   for(Int_t i=i1; i<i2; i++) {
233     if(ped>0.0) {
234       yd = Double_t(y[i]) - ped;
235     } else {
236       yd = Double_t(y[i]);
237     }
238     if(yd < yMinUnderPed) continue;
239
240     if(i>i1 && nPointsOut>0){
241       tdiff = t[i] - tn[nPointsOut-1];
242       //      printf(" i %i : nPointsOut %i : tdiff %6.2f : tn[nPointsOut] %6.2f \n", i,nPointsOut, tdiff, tn[nPointsOut-1]);
243       if(tdiff>1.) {
244      // discard previous points if its are before maximum point and with gap>1
245         if(i<nmax ) {
246           nPointsOut = 0; // nPointsOut--;
247      // if point with gap after maximum - finish selection
248         } else if(i>=nmax ) {
249           break;
250         }
251       }
252      // Far away from maximum
253      //if(i-nmax > Int_t(5*tau))              break;
254     }
255     tn[nPointsOut] = Double_t(t[i]);
256     yn[nPointsOut] = yd;
257     //printf("i %i : nPointsOut %i : tn %6.2f : yn %6.2f \n", i, nPointsOut, tn[nPointsOut], yn[nPointsOut]); 
258     nPointsOut++;
259   }
260   //printf(" nmax %i : nPointsIn %i :nPointsOut %i i1 %i \n", nmax, nPointsIn, nPointsOut, i1);
261 }
262
263 void AliCaloFastAltroFitv0::FastFit(const Double_t* t, const Double_t* y, const Int_t nPoints, 
264                                     const Double_t sig, const Double_t tau,
265                                     Double_t &amp, Double_t &eamp, Double_t &t0, Double_t &et0, Double_t &chi2)
266 {
267   // Static function
268   // It is case of n=k=2 : fnn = x*x*exp(2 - 2*x)
269   // Input: 
270   //nPoints  - number of points 
271   //   t[]   - array of time bins
272   //   y[]   - array of amplitudes after pedestal subtractions;
273   //   sig   - error of amplitude measurement (one value for all channels)
274   //   tau   - filter time response (in timebin units)
275   // Output:
276   //       amp - amplitude at t0;
277   //      eamp - error of amplitude; 
278   //        t0 - time of max amplitude; 
279   //       et0 - error of t0;
280   //      chi2 - chi2
281   static Double_t xx; // t/tau
282   static Double_t a, b, c;
283   static Double_t f02, f12, f22;    // functions
284   static Double_t f02d, f12d, f22d; // functions derivations
285
286   chi2 = -1.;
287
288   if(nPoints<2) {
289     printf(" AliCaloFastAltroFitv0::FastFit : nPoints<=%i \n", nPoints); 
290     return;
291   }
292
293   a = b = c = 0.0;
294   for(Int_t i=0; i<nPoints; i++){
295     xx  = t[i]/tau;
296     f02 = exp(-2.*xx);
297     f12 = xx*f02;
298     f22 = xx*f12;
299     // Derivations
300     f02d = -2.*f02;
301     f12d = f02 - 2.*f12;
302     f22d = 2.*(f12 - f22); 
303     //
304     a += f02d * y[i];
305     b -= 2.*f12d * y[i];
306     c += f22d * y[i];
307   }
308   Double_t t01=0.0, t02=0.0;
309   Double_t amp1=0.0, amp2=0.0, chi21=0.0, chi22=0.0;
310   if(QuadraticRoots(a,b,c, t01,t02)) {
311     t01 *= tau;
312     t02 *= tau;
313     Amplitude(t,y,nPoints, sig, tau, t01, amp1, chi21);
314     Amplitude(t,y,nPoints, sig, tau, t02, amp2, chi22);
315     if(0) {
316       printf(" t01 %f : t02 %f \n", t01, t02);
317       printf(" amp1 %f : amp2 %f \n", amp1, amp2);
318       printf(" chi21 %f : chi22 %f \n", chi21, chi22);
319     }
320     // t0 less on one tau with comparing with value from  "canonical equation"
321     amp  = amp1;
322     t0   = t01;
323     chi2 = chi21; 
324     if(chi21 > chi22) {
325       amp  = amp2;
326       t0   = t02; 
327       chi2 = chi22; 
328     }
329     if(tau<3.) { // EMCAL case : small tau 
330       t0 += -0.03; // Discard bias in t0
331       Amplitude(t,y,nPoints, sig, tau, t0, amp, chi2);
332     }
333     CalculateParsErrors(t, y, nPoints, sig, tau, amp, t0, eamp, et0);
334
335     // Fill1();
336     
337     // DrawFastFunction(amp, t0, fUtils->GetPedestalValue(), "1");
338     //    DrawFastFunction(amp1, t01, fUtils->GetPedestalValue(), "1");
339     // DrawFastFunction(amp2, t02, fUtils->GetPedestalValue(), "2");
340   } else {
341     chi2 = t01; // no roots, bad fit - negative chi2
342   }
343 }
344
345 Bool_t AliCaloFastAltroFitv0::QuadraticRoots(const Double_t a, const Double_t b, const Double_t c, 
346                                              Double_t &x1, Double_t &x2)
347 {
348   // Resolve quadratic equations a*x**2 + b*x + c
349   //printf(" a %12.5e b %12.5e c %12.5e \n", a, b, c);
350   static Double_t dtmp = 0.0, dtmpCut = -1.e-6;
351   static Int_t iWarning=0, iNoSolution=0;
352   dtmp = b*b - 4.*a*c;
353
354   if(dtmp>=dtmpCut && dtmp<0.0) {
355     if(iWarning<5 || iWarning%1000==0)
356       printf("<W> %i small neg. sq. : dtmp %12.5e \n", iWarning, dtmp);
357     iWarning++;
358     dtmp = 0.0;
359   }
360   if(dtmp>=0.0) {
361     dtmp = sqrt(dtmp);
362     x1   = (-b + dtmp) / (2.*a);
363     x2   = (-b - dtmp) / (2.*a);
364
365     //    printf(" x1 %f : x2 %f \n", x1, x2);
366     return kTRUE;
367   } else {
368     x1 = dtmp;
369     if(iNoSolution<5 || iNoSolution%1000==0)
370       printf("<No solution> %i neg. sq. : dtmp %12.5e \n", iNoSolution, dtmp);
371     iNoSolution++;
372     return kFALSE;
373   }
374 }
375
376 void AliCaloFastAltroFitv0::Amplitude(const Double_t* t,const Double_t* y,const Int_t nPoints, 
377                                       const Double_t sig, const Double_t tau, const Double_t t0, 
378                                       Double_t &amp, Double_t &chi2)
379 {  
380   // Calculate parameters error too - Mar 24,09
381   // sig is independent from points
382   amp = 0.;
383   Double_t x=0.0, f=0.0, den=0.0, f02;
384   for(Int_t i=0; i<nPoints; i++){
385     x    = (t[i] - t0)/tau;
386     f02  = exp(-2.*x);
387     f    = x*x*f02;     
388     amp += f*y[i];
389     den += f*f;
390   }
391   if(den>0.0) amp /= den;
392   //
393   // chi2 calculation
394   //
395   Double_t dy=0.0;
396   chi2=0.;
397   for(Int_t i=0; i<nPoints; i++){
398     x    = (t[i] - t0)/tau;
399     f02  = exp(-2.*x);
400     f    = amp*x*x*f02;
401     dy   = y[i]-f;
402     chi2 += dy*dy;
403     //    printf(" %i : y %f -> f %f : dy %f \n", i, y[i], f, dy); 
404   }
405   chi2 /= (sig*sig);
406 }
407
408 void AliCaloFastAltroFitv0::CalculateParsErrors(const Double_t* t, const Double_t* /*y*/, const Int_t nPoints, 
409                                                 const Double_t sig, const Double_t tau, 
410                                                 Double_t &amp, Double_t &t0, Double_t &eamp, Double_t &et0)
411 {
412   // Remember that fmax = exp(-n);
413   // fmax_nk = (n/k)**n*exp(-n) => n=k=2 => exp(-n) = exp(-2.)
414   static Double_t cc = exp(-2.);
415   //   static Double_t cc = exp(-fN); // mean(N)~1.5 ??
416
417   Double_t sumf2=0.0, sumfd2=0.0, x, f02, f12, f22, f22d;
418
419   for(Int_t i=0; i<nPoints; i++){
420     x    = (t[i] - t0)/tau;
421     f02  = exp(-2.*x);
422     f12  = x*f02;
423     f22  = x*f12;
424     sumf2 += f22 * f22;
425     //
426     f22d = 2.*(f12 - f22); 
427     sumfd2 += f22d * f22d;
428   }
429   et0  = (sig/amp)/sqrt(sumfd2);
430   eamp = sig/sqrt(sumf2);
431
432   amp  *= cc;
433   eamp *= cc;
434 }
435
436 //
437 // Drawing
438 //
439 TCanvas* AliCaloFastAltroFitv0::DrawFastFunction()
440 {
441   // QA of fitting
442   if(fNfit<=0) return 0; // no points
443
444   static TCanvas *c = 0;
445   if(c==0) {
446     c =  new TCanvas("fastFun","fastFun",800,600);
447   }
448
449   c->cd();
450   
451   Double_t* eamp = new Double_t[fNfit];
452   Double_t* et   = new Double_t[fNfit];
453
454   for(Int_t i=0; i<fNfit; i++) {
455     eamp[i] = fSig;
456     et[i]   = 0.0;
457   }
458
459   TGraphErrors *gr = new TGraphErrors(fNfit, fTfit,fAmpfit, et,eamp);
460   gr->Draw("Ap");
461   gr->SetTitle(Form("Fast Fit : #chi^{2}/ndf = %8.2f / %i", GetChi2(), GetNDF()));
462   gr->GetHistogram()->SetXTitle(" time bin ");
463   gr->GetHistogram()->SetYTitle(" amplitude ");
464
465   if(fStdFun==0) {
466      fStdFun = new TF1("stdFun", StdResponseFunction, 0., fTfit[fNfit-1]+2., 5);    
467      fStdFun->SetParNames("amp","t0","tau","N","ped");
468   }
469   fStdFun->SetParameter(0, GetEnergy());
470   fStdFun->SetParameter(1, GetTime() + GetTau());
471   fStdFun->SetParameter(2, GetTau());  // 
472   fStdFun->SetParameter(3, GetN());    // 2
473   fStdFun->SetParameter(4, 0.);  // 
474   
475   fStdFun->SetLineColor(kBlue);
476   fStdFun->SetLineWidth(1);
477
478   fStdFun->Draw("same");
479
480   delete [] eamp;
481   delete [] et;
482
483   c->Update();
484
485   return c;
486 }
487
488 Double_t AliCaloFastAltroFitv0::StdResponseFunction(const Double_t *x, const Double_t *par)
489 {
490   // Static Standard Response Function : 
491   // look to Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
492   // Using for drawing only.
493   // 
494   // Shape of the electronics raw reponse:
495   // It is a semi-gaussian, 2nd order Gamma function (n=2) of the general form
496   //
497   // t' = (t - t0 + tau) / tau
498   // F = A * t**N * exp( N * ( 1 - t) )   for t >= 0
499   // F = 0                                for t < 0 
500   //
501   // parameters:
502   // A:   par[0]   // Amplitude = peak value
503   // t0:  par[1]
504   // tau: par[2]
505   // N:   par[3]
506   // ped: par[4]
507   //
508   static Double_t signal , tau, n, ped, xx;
509
510   tau = par[2];
511   n   = par[3];
512   ped = par[4];
513   xx = ( x[0] - par[1] + tau ) / tau ;
514
515   if (xx <= 0) 
516     signal = ped ;  
517   else {  
518     signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ; 
519   }
520   return signal ;  
521 }