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