PHOS/AliPHOSFastAltroFit.cxx

   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
  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 ClassImp(AliPHOSFastAltroFit)
  46
  47 //____________________________________________________________________________
  48   AliPHOSFastAltroFit:: AliPHOSFastAltroFit()
  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 //____________________________________________________________________________
  56 AliPHOSFastAltroFit::AliPHOSFastAltroFit(const char* name, const char* title, const Double_t tau)
  57   : TNamed(name, title),
  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)
  60 {
  61   if(strlen(name)==0) SetName("FastAltroFit");
  62 }
  63
  64 //____________________________________________________________________________
  65 AliPHOSFastAltroFit::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 //____________________________________________________________________________
  73 AliPHOSFastAltroFit::~AliPHOSFastAltroFit()
  74 {
  75   if(fTfit) delete [] fTfit;
  76   if(fAmpfit) delete [] fAmpfit;
  77 }
  78
  79 //____________________________________________________________________________
  80 AliPHOSFastAltroFit& AliPHOSFastAltroFit::operator= (const AliPHOSFastAltroFit &/*obj*/)
  81 {
  82   // Not implemented yet
  83   return *this;
  84 }
  85
  86 //____________________________________________________________________________
  87 void 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;
 107 }
 108
 109 void AliPHOSFastAltroFit::FastFit(Int_t* t, Int_t* y, Int_t n, Double_t sig, Double_t tau, Double_t ped)
 110 {
 111   Reset();
 112
 113   fSig = sig;
 114   fTau = tau;
 115   fPed = ped;
 116
 117   if(fTfit) delete [] fTfit;
 118   if(fAmpfit) delete [] fAmpfit;
 119
 120   fNfit   = 0;
 121   fTfit   = new Double_t[n];
 122   fAmpfit = new Double_t[n];
 123
 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]);
 128
 129   if(fNfit>=2) {
 130     FastFit(fTfit,fAmpfit,fNfit,sig,tau, fAmp,fAmpErr, fT0,fT0Err,fChi2);
 131
 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 }
 140
 141 void 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
 151 void AliPHOSFastAltroFit::GetFitResult(Double_t &amp,Double_t &eamp,Double_t &t0,Double_t &et0,
 152 Double_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
 162 void AliPHOSFastAltroFit::GetFittedPoints(Int_t &nfit, Double_t* ar[2])
 163 {
 164   nfit  = fNfit;
 165   ar[0] = fTfit;
 166   ar[1] = fAmpfit;
 167 }
 168
 169 void 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 {
 172   static Double_t yMinUnderPed=2.; // should be tune
 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);
 181   //i1 = 0;
 182   i1 = i1<0?0:i1;
 183   Int_t i2 = n;
 184
 185   nn = 0;
 186   Double_t yd=0.0, tdiff=0.0;;
 187   for(Int_t i=i1; i<i2; i++) {
 188     if(ped>0.0) {
 189       yd = Double_t(y[i]) - ped;
 190     } else {
 191       yd = Double_t(y[i]);
 192     }
 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++;
 214   }
 215   //printf(" nmax %i : n %i : nn %i i1 %i \n", nmax, n, nn, i1);
 216 }
 217
 218 void AliPHOSFastAltroFit::FastFit(Double_t* t, Double_t* y, Int_t n, Double_t sig, Double_t tau,
 219 Double_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;
 230   //        t0 - time of max amplitude;
 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
 235
 236   chi2 = -1.;
 237
 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     }
 279     if(tau<3.) { // EMCAL case : small tau
 280       t0 += -0.03;
 281       Amplitude(t,y,n, sig, tau, t0, amp, chi2);
 282     }
 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 {
 291     chi2 = t0_1; // no roots, bad fit - negative chi2
 292   }
 293 }
 294
 295 Bool_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);
 299   static Double_t dtmp = 0.0, dtmpCut = -1.e-6;
 300   static Int_t ierr=0;
 301   dtmp = b*b - 4.*a*c;
 302
 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);
 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 {
 316     x1 = dtmp;
 317     printf("QuadraticRoots : negative square : dtmp %12.5e \n", dtmp);
 318     return kFALSE;
 319   }
 320 }
 321
 322 void AliPHOSFastAltroFit::Amplitude(Double_t* t,Double_t* y,Int_t n, Double_t sig, Double_t tau,
 323 Double_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   }
 336   if(den>0.0) amp /= den;
 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;
 348     //    printf(" %i : y %f -> f %f : dy %f \n", i, y[i], f, dy);
 349   }
 350   chi2 /= (sig*sig);
 351 }
 352
 353 void 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)
 356 {
 357   // Remember that fmax = exp(-n);
 358   // fmax_nk = (n/k)**n*exp(-n) => n=k=2 => exp(-n) = exp(-2.)
 359   static Double_t cc = exp(-2.);
 360   //   static Double_t cc = exp(-fN); // mean(N)~1.5 ??
 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;
 366     f02  = exp(-2.*x);
 367     f12  = x*f02;
 368     f22  = x*f12;
 369     sumf2 += f22 * f22;
 370     //
 371     f22d = 2.*(f12 - f22);
 372     sumfd2 += f22d * f22d;
 373   }
 374   et0  = (sig/amp)/sqrt(sumfd2);
 375   eamp = sig/sqrt(sumf2);
 376
 377   amp  *= cc;
 378   eamp *= cc;
 379 }
 380
 381 //
 382 // Drawing
 383 //
 384 TCanvas* 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
 433 Double_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 }