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