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