1 /**************************************************************************
2 * Copyright(c) 1998-2010 ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 /* History of cvs commits:
22 //______________________________________________________
23 // Author : Aleksei Pavlinov; IHEP, Protvino, Russia
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 ..
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.
34 // --- AliRoot header files ---
35 #include "AliCaloFastAltroFitv0.h"
40 #include <TGraphErrors.h>
45 ClassImp(AliCaloFastAltroFitv0)
47 //____________________________________________________________________________
48 AliCaloFastAltroFitv0::AliCaloFastAltroFitv0()
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)
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)
62 if(strlen(name)==0) SetName("CaloFastAltroFitv0");
65 //____________________________________________________________________________
66 AliCaloFastAltroFitv0::AliCaloFastAltroFitv0(const AliCaloFastAltroFitv0 &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)
73 //____________________________________________________________________________
74 AliCaloFastAltroFitv0::~AliCaloFastAltroFitv0()
76 if(fTfit) delete [] fTfit;
77 if(fAmpfit) delete [] fAmpfit;
80 //____________________________________________________________________________
81 AliCaloFastAltroFitv0& AliCaloFastAltroFitv0::operator= (const AliCaloFastAltroFitv0 &/*obj*/)
83 // Not implemented yet
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)
90 // n 2 here and unused
99 CutRightPart(t,y,nPoints,tMax, ii);
103 fTfit = new Double_t[nPoints];
104 fAmpfit = new Double_t[nPoints];
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]);
113 FastFit(fTfit,fAmpfit,fNfit,sig,tau, fAmp,fAmpErr, fT0,fT0Err,fChi2);
122 Reset(); // What to do here => fT0 = fTfit[0]; fAmp = fAmpFit[0] ??
128 //____________________________________________________________________________
129 void AliCaloFastAltroFitv0::FastFit(TH1F* h, Double_t sig, Double_t tau, Double_t n,
130 Double_t ped, Double_t tMax)
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
138 Int_t nPoints = h->GetNbinsX();
139 if(nPoints<2) return; // Sep 07, 09
141 Int_t* t = new Int_t[nPoints];
142 Int_t* y = new Int_t[nPoints];
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);
154 FastFit(t,y,nPoints, sig,tau,n,ped, tMax);
156 if(fChi2<=0.0) fNoFit++;
162 void AliCaloFastAltroFitv0::Reset()
166 fAmp = fAmpErr = fT0 = fT0Err = 0.0;
170 if(fTfit) delete [] fTfit;
171 if(fAmpfit) delete [] fAmpfit;
176 void AliCaloFastAltroFitv0::GetFitResult(Double_t &,Double_t &eamp,Double_t &t0,Double_t &et0,
177 Double_t &chi2, Int_t &ndf) const
179 // Return results of fitting procedure
188 void AliCaloFastAltroFitv0::GetFittedPoints(Int_t &nfit, Double_t* ar[2]) const
197 void AliCaloFastAltroFitv0::CutRightPart(Int_t *t,Int_t *y,Int_t nPoints,Double_t tMax, Int_t &ii)
199 // Cut right part of altro sample : static function
201 for(Int_t i=0; i<nPoints; i++) {
203 if(tMax && tt <= Int_t(tMax)) {
209 if(0) printf(" ii %i -> ii %i : tMax %7.2f \n", nPoints, ii, tMax);
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)
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++){
225 Int_t i1 = nmax - Int_t(tau);
231 Double_t yd=0.0, tdiff=0.0;;
232 for(Int_t i=i1; i<i2; i++) {
234 yd = Double_t(y[i]) - ped;
238 if(yd < yMinUnderPed) continue;
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]);
244 // discard previous points if its are before maximum point and with gap>1
246 nPointsOut = 0; // nPointsOut--;
247 // if point with gap after maximum - finish selection
248 } else if(i>=nmax ) {
252 // Far away from maximum
253 //if(i-nmax > Int_t(5*tau)) break;
255 tn[nPointsOut] = Double_t(t[i]);
257 //printf("i %i : nPointsOut %i : tn %6.2f : yn %6.2f \n", i, nPointsOut, tn[nPointsOut], yn[nPointsOut]);
260 //printf(" nmax %i : nPointsIn %i :nPointsOut %i i1 %i \n", nmax, nPointsIn, nPointsOut, i1);
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 &, Double_t &eamp, Double_t &t0, Double_t &et0, Double_t &chi2)
268 // It is case of n=k=2 : fnn = x*x*exp(2 - 2*x)
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)
276 // amp - amplitude at t0;
277 // eamp - error of amplitude;
278 // t0 - time of max amplitude;
279 // et0 - error of t0;
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
289 printf(" AliCaloFastAltroFitv0::FastFit : nPoints<=%i \n", nPoints);
294 for(Int_t i=0; i<nPoints; i++){
302 f22d = 2.*(f12 - f22);
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)) {
313 Amplitude(t,y,nPoints, sig, tau, t01, amp1, chi21);
314 Amplitude(t,y,nPoints, sig, tau, t02, amp2, chi22);
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);
320 // t0 less on one tau with comparing with value from "canonical equation"
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);
333 CalculateParsErrors(t, y, nPoints, sig, tau, amp, t0, eamp, et0);
337 // DrawFastFunction(amp, t0, fUtils->GetPedestalValue(), "1");
338 // DrawFastFunction(amp1, t01, fUtils->GetPedestalValue(), "1");
339 // DrawFastFunction(amp2, t02, fUtils->GetPedestalValue(), "2");
341 chi2 = t01; // no roots, bad fit - negative chi2
345 Bool_t AliCaloFastAltroFitv0::QuadraticRoots(const Double_t a, const Double_t b, const Double_t c,
346 Double_t &x1, Double_t &x2)
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;
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);
362 x1 = (-b + dtmp) / (2.*a);
363 x2 = (-b - dtmp) / (2.*a);
365 // printf(" x1 %f : x2 %f \n", x1, x2);
369 if(iNoSolution<5 || iNoSolution%1000==0)
370 printf("<No solution> %i neg. sq. : dtmp %12.5e \n", iNoSolution, dtmp);
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 &, Double_t &chi2)
380 // Calculate parameters error too - Mar 24,09
381 // sig is independent from points
383 Double_t x=0.0, f=0.0, den=0.0, f02;
384 for(Int_t i=0; i<nPoints; i++){
391 if(den>0.0) amp /= den;
397 for(Int_t i=0; i<nPoints; i++){
403 // printf(" %i : y %f -> f %f : dy %f \n", i, y[i], f, dy);
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 &, Double_t &t0, Double_t &eamp, Double_t &et0)
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 ??
417 Double_t sumf2=0.0, sumfd2=0.0, x, f02, f12, f22, f22d;
419 for(Int_t i=0; i<nPoints; i++){
426 f22d = 2.*(f12 - f22);
427 sumfd2 += f22d * f22d;
429 et0 = (sig/amp)/sqrt(sumfd2);
430 eamp = sig/sqrt(sumf2);
439 TCanvas* AliCaloFastAltroFitv0::DrawFastFunction()
442 if(fNfit<=0) return 0; // no points
444 static TCanvas *c = 0;
446 c = new TCanvas("fastFun","fastFun",800,600);
451 Double_t* eamp = new Double_t[fNfit];
452 Double_t* et = new Double_t[fNfit];
454 for(Int_t i=0; i<fNfit; i++) {
459 TGraphErrors *gr = new TGraphErrors(fNfit, fTfit,fAmpfit, et,eamp);
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 ");
466 fStdFun = new TF1("stdFun", StdResponseFunction, 0., fTfit[fNfit-1]+2., 5);
467 fStdFun->SetParNames("amp","t0","tau","N","ped");
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.); //
475 fStdFun->SetLineColor(kBlue);
476 fStdFun->SetLineWidth(1);
478 fStdFun->Draw("same");
488 Double_t AliCaloFastAltroFitv0::StdResponseFunction(const Double_t *x, const Double_t *par)
490 // Static Standard Response Function :
491 // look to Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
492 // Using for drawing only.
494 // Shape of the electronics raw reponse:
495 // It is a semi-gaussian, 2nd order Gamma function (n=2) of the general form
497 // t' = (t - t0 + tau) / tau
498 // F = A * t**N * exp( N * ( 1 - t) ) for t >= 0
502 // A: par[0] // Amplitude = peak value
508 static Double_t signal , tau, n, ped, xx;
513 xx = ( x[0] - par[1] + tau ) / tau ;
518 signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;