Line data Source code
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 42 : ClassImp(AliCaloFastAltroFitv0)
46 :
47 : //____________________________________________________________________________
48 : AliCaloFastAltroFitv0::AliCaloFastAltroFitv0()
49 0 : : TNamed(),
50 0 : fSig(0),fTau(0),fN(0),fPed(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0)
51 0 : ,fNoFit(0),fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
52 0 : {
53 0 : }
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 0 : : TNamed(name, title),
59 0 : fSig(sig),fTau(tau),fN(n),fPed(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0)
60 0 : ,fNoFit(0),fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
61 0 : {
62 0 : if(strlen(name)==0) SetName("CaloFastAltroFitv0");
63 0 : }
64 :
65 : //____________________________________________________________________________
66 : AliCaloFastAltroFitv0::AliCaloFastAltroFitv0(const AliCaloFastAltroFitv0 &obj)
67 0 : : TNamed(obj),
68 0 : fSig(0),fTau(0),fN(2.),fPed(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0)
69 0 : ,fNoFit(0),fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
70 0 : {
71 0 : }
72 :
73 : //____________________________________________________________________________
74 : AliCaloFastAltroFitv0::~AliCaloFastAltroFitv0()
75 0 : {
76 0 : if(fTfit) delete [] fTfit;
77 0 : if(fAmpfit) delete [] fAmpfit;
78 0 : }
79 :
80 : //____________________________________________________________________________
81 : AliCaloFastAltroFitv0& AliCaloFastAltroFitv0::operator= (const AliCaloFastAltroFitv0 &/*obj*/)
82 : {
83 : // Not implemented yet
84 0 : 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 0 : Reset();
91 :
92 0 : fSig = sig;
93 0 : fTau = tau;
94 0 : fPed = ped;
95 :
96 0 : Int_t ii=0;
97 0 : CutRightPart(t,y,nPoints,tMax, ii);
98 0 : nPoints = ii;
99 :
100 0 : fNfit = 0;
101 0 : fTfit = new Double_t[nPoints];
102 0 : fAmpfit = new Double_t[nPoints];
103 :
104 :
105 0 : 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 0 : if(fNfit>=2) {
111 0 : FastFit(fTfit,fAmpfit,fNfit,sig,tau, fAmp,fAmpErr, fT0,fT0Err,fChi2);
112 :
113 0 : if(fChi2> 0.0) {
114 0 : fNDF = fNfit - 2;
115 0 : } else {
116 0 : fNDF = 0;
117 0 : fNoFit++;
118 : }
119 : } else if(fNfit==1){
120 0 : Reset(); // What to do here => fT0 = fTfit[0]; fAmp = fAmpFit[0] ??
121 : } else {
122 : Reset();
123 : }
124 0 : }
125 :
126 : //____________________________________________________________________________
127 : void AliCaloFastAltroFitv0::FastFit(TH1F* h, Double_t sig, Double_t tau, Double_t n,
128 : Double_t ped, Double_t tMax)
129 : {
130 : // Service method for convinience only
131 : // h - hist with altro response, could have empty bin
132 : // and center of bin could be different from bin number
133 0 : Reset();
134 :
135 0 : if(h==0) return;
136 0 : Int_t nPoints = h->GetNbinsX();
137 0 : if(nPoints<2) return; // Sep 07, 09
138 :
139 0 : Int_t* t = new Int_t[nPoints];
140 0 : Int_t* y = new Int_t[nPoints];
141 :
142 : Int_t nPositive=0;
143 0 : for(Int_t i=0; i<nPoints; i++) {
144 0 : if(h->GetBinContent(i+1) > 0.0){ // Get only positive
145 0 : y[nPositive] = Int_t(h->GetBinContent(i+1));
146 0 : t[nPositive] = Int_t(h->GetBinCenter(i+1)+0.0001);
147 0 : nPositive++;
148 0 : }
149 : }
150 :
151 0 : if(nPositive >= 2) {
152 0 : FastFit(t,y,nPoints, sig,tau,n,ped, tMax);
153 0 : }
154 0 : if(fChi2<=0.0) fNoFit++;
155 :
156 0 : delete [] t;
157 0 : delete [] y;
158 0 : }
159 :
160 : void AliCaloFastAltroFitv0::Reset()
161 : {
162 : // Reset variables
163 0 : fSig = fTau = 0.0;
164 0 : fAmp = fAmpErr = fT0 = fT0Err = 0.0;
165 0 : fChi2 = -1.;
166 0 : fNDF = fNfit = 0;
167 :
168 0 : if(fTfit) delete [] fTfit;
169 0 : if(fAmpfit) delete [] fAmpfit;
170 0 : fTfit = fAmpfit = 0;
171 0 : }
172 :
173 :
174 : void AliCaloFastAltroFitv0::GetFitResult(Double_t &,Double_t &eamp,Double_t &t0,Double_t &et0,
175 : Double_t &chi2, Int_t &ndf) const
176 : {
177 : // Return results of fitting procedure
178 0 : amp = fAmp;
179 0 : eamp = fAmpErr;
180 0 : t0 = fT0;
181 0 : et0 = fT0Err;
182 0 : chi2 = fChi2;
183 0 : ndf = fNDF;
184 0 : }
185 :
186 : void AliCaloFastAltroFitv0::GetFittedPoints(Int_t &nfit, Double_t* ar[2]) const
187 : {
188 0 : nfit = fNfit;
189 0 : ar[0] = fTfit;
190 0 : ar[1] = fAmpfit;
191 0 : }
192 : //
193 : // static functions
194 : //
195 : void AliCaloFastAltroFitv0::CutRightPart(Int_t *t,Int_t *y,Int_t nPoints,Double_t tMax, Int_t &ii)
196 : {
197 : // Cut right part of altro sample : static function
198 : Int_t tt=0;
199 0 : for(Int_t i=0; i<nPoints; i++) {
200 0 : tt = t[i];
201 0 : if(tMax && tt <= Int_t(tMax)) {
202 0 : t[ii] = tt;
203 0 : y[ii] = y[i];
204 0 : ii++;
205 0 : }
206 : }
207 : if(0) printf(" ii %i -> ii %i : tMax %7.2f \n", nPoints, ii, tMax);
208 0 : }
209 :
210 : void 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)
212 : {
213 : // Pedestal deduction if ped is positive : static function
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;
217 0 : for(Int_t i=0; i<nPoints; i++){
218 0 : if(y[i]>ymax) {
219 : ymax = y[i];
220 : nmax = i;
221 0 : }
222 : }
223 0 : Int_t i1 = nmax - Int_t(tau);
224 : //i1 = 0;
225 0 : i1 = i1<0?0:i1;
226 : Int_t i2 = nPoints;
227 :
228 0 : nPointsOut = 0;
229 : Double_t yd=0.0, tdiff=0.0;;
230 0 : for(Int_t i=i1; i<i2; i++) {
231 0 : if(ped>0.0) {
232 0 : yd = Double_t(y[i]) - ped;
233 0 : } else {
234 : yd = Double_t(y[i]);
235 : }
236 0 : if(yd < yMinUnderPed) continue;
237 :
238 0 : if(i>i1 && nPointsOut>0){
239 0 : 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]);
241 0 : if(tdiff>1.) {
242 : // discard previous points if its are before maximum point and with gap>1
243 0 : if(i<nmax ) {
244 0 : nPointsOut = 0; // nPointsOut--;
245 : // if point with gap after maximum - finish selection
246 0 : } else if(i>=nmax ) {
247 0 : break;
248 : }
249 : }
250 : // Far away from maximum
251 : //if(i-nmax > Int_t(5*tau)) break;
252 : }
253 0 : tn[nPointsOut] = Double_t(t[i]);
254 0 : yn[nPointsOut] = yd;
255 : //printf("i %i : nPointsOut %i : tn %6.2f : yn %6.2f \n", i, nPointsOut, tn[nPointsOut], yn[nPointsOut]);
256 0 : nPointsOut++;
257 0 : }
258 : //printf(" nmax %i : nPointsIn %i :nPointsOut %i i1 %i \n", nmax, nPointsIn, nPointsOut, i1);
259 0 : }
260 :
261 : void AliCaloFastAltroFitv0::FastFit(const Double_t* t, const Double_t* y, const Int_t nPoints,
262 : const Double_t sig, const Double_t tau,
263 : Double_t &, Double_t &eamp, Double_t &t0, Double_t &et0, Double_t &chi2)
264 : {
265 : // Static function
266 : // It is case of n=k=2 : fnn = x*x*exp(2 - 2*x)
267 : // Input:
268 : //nPoints - number of points
269 : // t[] - array of time bins
270 : // y[] - array of amplitudes after pedestal subtractions;
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;
275 : // eamp - error of amplitude;
276 : // t0 - time of max amplitude;
277 : // et0 - error of t0;
278 : // chi2 - chi2
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 0 : chi2 = -1.;
285 :
286 0 : if(nPoints<2) {
287 0 : printf(" AliCaloFastAltroFitv0::FastFit : nPoints<=%i \n", nPoints);
288 0 : return;
289 : }
290 :
291 0 : a = b = c = 0.0;
292 0 : for(Int_t i=0; i<nPoints; i++){
293 0 : xx = t[i]/tau;
294 0 : f02 = exp(-2.*xx);
295 0 : f12 = xx*f02;
296 0 : f22 = xx*f12;
297 : // Derivations
298 0 : f02d = -2.*f02;
299 0 : f12d = f02 - 2.*f12;
300 0 : f22d = 2.*(f12 - f22);
301 : //
302 0 : a += f02d * y[i];
303 0 : b -= 2.*f12d * y[i];
304 0 : c += f22d * y[i];
305 : }
306 0 : Double_t t01=0.0, t02=0.0;
307 0 : Double_t amp1=0.0, amp2=0.0, chi21=0.0, chi22=0.0;
308 0 : if(QuadraticRoots(a,b,c, t01,t02)) {
309 0 : t01 *= tau;
310 0 : t02 *= tau;
311 0 : Amplitude(t,y,nPoints, sig, tau, t01, amp1, chi21);
312 0 : Amplitude(t,y,nPoints, sig, tau, t02, amp2, chi22);
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 0 : amp = amp1;
320 0 : t0 = t01;
321 0 : chi2 = chi21;
322 0 : if(chi21 > chi22) {
323 0 : amp = amp2;
324 0 : t0 = t02;
325 0 : chi2 = chi22;
326 0 : }
327 0 : if(tau<3.) { // EMCAL case : small tau
328 0 : t0 += -0.03; // Discard bias in t0
329 0 : Amplitude(t,y,nPoints, sig, tau, t0, amp, chi2);
330 0 : }
331 0 : CalculateParsErrors(t, y, nPoints, sig, tau, amp, t0, eamp, et0);
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 0 : } else {
339 0 : chi2 = t01; // no roots, bad fit - negative chi2
340 : }
341 0 : }
342 :
343 : Bool_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;
349 : static Int_t iWarning=0, iNoSolution=0;
350 0 : dtmp = b*b - 4.*a*c;
351 :
352 0 : if(dtmp>=dtmpCut && dtmp<0.0) {
353 0 : if(iWarning<5 || iWarning%1000==0)
354 0 : printf("<W> %i small neg. sq. : dtmp %12.5e \n", iWarning, dtmp);
355 0 : iWarning++;
356 0 : dtmp = 0.0;
357 0 : }
358 0 : if(dtmp>=0.0) {
359 0 : dtmp = sqrt(dtmp);
360 0 : x1 = (-b + dtmp) / (2.*a);
361 0 : x2 = (-b - dtmp) / (2.*a);
362 :
363 : // printf(" x1 %f : x2 %f \n", x1, x2);
364 0 : return kTRUE;
365 : } else {
366 0 : x1 = dtmp;
367 0 : if(iNoSolution<5 || iNoSolution%1000==0)
368 0 : printf("<No solution> %i neg. sq. : dtmp %12.5e \n", iNoSolution, dtmp);
369 0 : iNoSolution++;
370 0 : return kFALSE;
371 : }
372 0 : }
373 :
374 : void AliCaloFastAltroFitv0::Amplitude(const Double_t* t,const Double_t* y,const Int_t nPoints,
375 : const Double_t sig, const Double_t tau, const Double_t t0,
376 : Double_t &, Double_t &chi2)
377 : {
378 : // Calculate parameters error too - Mar 24,09
379 : // sig is independent from points
380 0 : amp = 0.;
381 : Double_t x=0.0, f=0.0, den=0.0, f02;
382 0 : for(Int_t i=0; i<nPoints; i++){
383 0 : x = (t[i] - t0)/tau;
384 0 : f02 = exp(-2.*x);
385 0 : f = x*x*f02;
386 0 : amp += f*y[i];
387 0 : den += f*f;
388 : }
389 0 : if(den>0.0) amp /= den;
390 : //
391 : // chi2 calculation
392 : //
393 : Double_t dy=0.0;
394 0 : chi2=0.;
395 0 : for(Int_t i=0; i<nPoints; i++){
396 0 : x = (t[i] - t0)/tau;
397 0 : f02 = exp(-2.*x);
398 0 : f = amp*x*x*f02;
399 0 : dy = y[i]-f;
400 0 : chi2 += dy*dy;
401 : // printf(" %i : y %f -> f %f : dy %f \n", i, y[i], f, dy);
402 : }
403 0 : chi2 /= (sig*sig);
404 0 : }
405 :
406 : void AliCaloFastAltroFitv0::CalculateParsErrors(const Double_t* t, const Double_t* /*y*/, const Int_t nPoints,
407 : const Double_t sig, const Double_t tau,
408 : Double_t &, 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 0 : 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 :
417 0 : for(Int_t i=0; i<nPoints; i++){
418 0 : x = (t[i] - t0)/tau;
419 0 : f02 = exp(-2.*x);
420 0 : f12 = x*f02;
421 0 : f22 = x*f12;
422 0 : sumf2 += f22 * f22;
423 : //
424 0 : f22d = 2.*(f12 - f22);
425 0 : sumfd2 += f22d * f22d;
426 : }
427 0 : et0 = (sig/amp)/sqrt(sumfd2);
428 0 : eamp = sig/sqrt(sumf2);
429 :
430 0 : amp *= cc;
431 0 : eamp *= cc;
432 0 : }
433 :
434 : //
435 : // Drawing
436 : //
437 : TCanvas* AliCaloFastAltroFitv0::DrawFastFunction()
438 : {
439 : // QA of fitting
440 0 : if(fNfit<=0) return 0; // no points
441 :
442 : static TCanvas *c = 0;
443 0 : if(c==0) {
444 0 : c = new TCanvas("fastFun","fastFun",800,600);
445 0 : }
446 :
447 0 : c->cd();
448 :
449 0 : Double_t* eamp = new Double_t[fNfit];
450 0 : Double_t* et = new Double_t[fNfit];
451 :
452 0 : for(Int_t i=0; i<fNfit; i++) {
453 0 : eamp[i] = fSig;
454 0 : et[i] = 0.0;
455 : }
456 :
457 0 : TGraphErrors *gr = new TGraphErrors(fNfit, fTfit,fAmpfit, et,eamp);
458 0 : gr->Draw("Ap");
459 0 : gr->SetTitle(Form("Fast Fit : #chi^{2}/ndf = %8.2f / %i", GetChi2(), GetNDF()));
460 0 : gr->GetHistogram()->SetXTitle(" time bin ");
461 0 : gr->GetHistogram()->SetYTitle(" amplitude ");
462 :
463 0 : if(fStdFun==0) {
464 0 : fStdFun = new TF1("stdFun", StdResponseFunction, 0., fTfit[fNfit-1]+2., 5);
465 0 : fStdFun->SetParNames("amp","t0","tau","N","ped");
466 0 : }
467 0 : fStdFun->SetParameter(0, GetEnergy());
468 0 : fStdFun->SetParameter(1, GetTime() + GetTau());
469 0 : fStdFun->SetParameter(2, GetTau()); //
470 0 : fStdFun->SetParameter(3, GetN()); // 2
471 0 : fStdFun->SetParameter(4, 0.); //
472 :
473 0 : fStdFun->SetLineColor(kBlue);
474 0 : fStdFun->SetLineWidth(1);
475 :
476 0 : fStdFun->Draw("same");
477 :
478 0 : delete [] eamp;
479 0 : delete [] et;
480 :
481 0 : c->Update();
482 :
483 0 : return c;
484 0 : }
485 :
486 : Double_t AliCaloFastAltroFitv0::StdResponseFunction(const Double_t *x, const Double_t *par)
487 : {
488 : // Static Standard Response Function :
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 0 : tau = par[2];
509 0 : n = par[3];
510 0 : ped = par[4];
511 0 : xx = ( x[0] - par[1] + tau ) / tau ;
512 :
513 0 : if (xx <= 0)
514 0 : signal = ped ;
515 : else {
516 0 : signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;
517 : }
518 0 : return signal ;
519 : }
|
