RAW/AliCaloFastAltroFitv0.cxx

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