[u/mrichter/AliRoot.git] / PHOS / AliPHOSFastAltroFit.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id:$ */

/* History of cvs commits:
 * $Log$
 */

//______________________________________________________
// Author : Aleksei Pavlinov; IHEP, Protvino, Russia
// Feb 17, 2009
// Implementation of fit procedure from ALICE-INT-2008-026:
// "Time and amplitude reconstruction from sampling 
//  measurements of the PHOS signal profile"
//  M.Yu.Bogolyubsky and ..
//
//  Fit by function en*x*x*exp(-2.*x): x = (t-t0)/tau.
//  Tme main goal is fast estimation of amplitude.
//

// --- AliRoot header files ---
#include "AliPHOSFastAltroFit.h"

#include <TH1.h>
#include <TF1.h>
#include <TCanvas.h>
#include <TGraphErrors.h>
#include <TMath.h>

#include <math.h>

ClassImp(AliPHOSFastAltroFit)

//____________________________________________________________________________
  AliPHOSFastAltroFit:: AliPHOSFastAltroFit() 
: TNamed(), 
  fSig(0),fTau(0),fN(0),fPed(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0)
,fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
{
}

//____________________________________________________________________________
AliPHOSFastAltroFit::AliPHOSFastAltroFit(const char* name, const char* title, const Double_t tau)
  : TNamed(name, title), 
  fSig(0),fTau(tau),fN(2),fPed(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0) 
 ,fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
{
  if(strlen(name)==0) SetName("FastAltroFit");
}

//____________________________________________________________________________
AliPHOSFastAltroFit::AliPHOSFastAltroFit(const AliPHOSFastAltroFit &obj)
  : TNamed(obj), 
  fSig(0),fTau(0),fN(2),fPed(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0) 
 ,fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
{
}

//____________________________________________________________________________
AliPHOSFastAltroFit::~AliPHOSFastAltroFit() 
{  
  if(fTfit) delete [] fTfit;
  if(fAmpfit) delete [] fAmpfit;
}

//____________________________________________________________________________
AliPHOSFastAltroFit& AliPHOSFastAltroFit::operator= (const AliPHOSFastAltroFit &/*obj*/)
{
  // Not implemented yet
  return *this;
}

//____________________________________________________________________________
void AliPHOSFastAltroFit::FastFit(TH1F* h, Double_t sig, Double_t tau, Double_t ped)
{
  // Service method for convinience only  
  Reset();

  if(h==0) return;
  Int_t n = h->GetNbinsX();
  if(n<=0) return;

  Int_t* t = new Int_t[n];
  Int_t* y = new Int_t[n];

  for(Int_t i=0; i<n; i++) {
    t[i] = Int_t(h->GetBinCenter(i+1));
    y[i] = Int_t(h->GetBinContent(i+1));
  }
  FastFit(t,y,n, sig,tau,ped);

  delete [] t;
  delete [] y;
}

void AliPHOSFastAltroFit::FastFit(Int_t* t, Int_t* y, Int_t n, Double_t sig, Double_t tau, Double_t ped)
{
  Reset();

  fSig = sig;
  fTau = tau;
  fPed = ped;

  if(fTfit) delete [] fTfit;
  if(fAmpfit) delete [] fAmpfit;

  fNfit   = 0;
  fTfit   = new Double_t[n]; 
  fAmpfit = new Double_t[n]; 

  DeductPedestal(t,y,n,  tau,ped,  fTfit,fAmpfit,fNfit);
  //  printf(" n %i : fNfit %i : ped %f \n", n, fNfit, ped);
  // for(int i=0; i<fNfit; i++) 
  // printf(" i %i : fAmpfit %7.2f : fTfit %7.2f \n", i, fAmpfit[i], fTfit[i]); 

  if(fNfit>=2) {
    FastFit(fTfit,fAmpfit,fNfit,sig,tau, fAmp,fAmpErr, fT0,fT0Err,fChi2);

    if(fChi2> 0.0) fNDF = fNfit - 2;
    else           fNDF = 0; 
  } else if(fNfit==1){
    Reset(); // What to do here => fT0 = fTfit[0]; fAmp = fAmpFit[0] ??
  } else {
    Reset();
  }
}

void AliPHOSFastAltroFit::Reset()
{
  fSig  = fTau = 0.0;
  fAmp  = fAmpErr = fT0 = fT0Err = 0.0;
  fChi2 = -1.;
  fNDF  = fNfit = 0;
  fTfit = fAmpfit = 0;
}


void AliPHOSFastAltroFit::GetFitResult(Double_t &amp,Double_t &eamp,Double_t &t0,Double_t &et0, 
Double_t &chi2, Int_t &ndf) 
{
  amp  = fAmp;
  eamp = fAmpErr;
  t0   = fT0;
  et0  = fT0Err;
  chi2 = fChi2;
  ndf  = fNDF;
}

void AliPHOSFastAltroFit::GetFittedPoints(Int_t &nfit, Double_t* ar[2])
{
  nfit  = fNfit;
  ar[0] = fTfit;
  ar[1] = fAmpfit;
}

void AliPHOSFastAltroFit::DeductPedestal(Int_t* t, Int_t* y, Int_t n, Double_t tau, Double_t ped, 
  Double_t* tn, Double_t* yn, Int_t &nn)
{
  static Double_t yMinUnderPed=2.; // should be tune
  Int_t ymax=0, nmax=0;
  for(Int_t i=0; i<n; i++){
    if(y[i]>ymax) {
      ymax = y[i];
      nmax = i;
    }
  }
  Int_t i1 = nmax - Int_t(tau);
  //i1 = 0;
  i1 = i1<0?0:i1;
  Int_t i2 = n;

  nn = 0;
  Double_t yd=0.0, tdiff=0.0;;
  for(Int_t i=i1; i<i2; i++) {
    if(ped>0.0) {
      yd = Double_t(y[i]) - ped;
    } else {
      yd = Double_t(y[i]);
    }
    if(yd < yMinUnderPed) continue;

    if(i>i1 && nn>0){
      tdiff = t[i] - tn[nn-1];
      //      printf(" i %i : nn %i : tdiff %6.2f : tn[nn] %6.2f \n", i,nn, tdiff, tn[nn-1]);
      if(tdiff>1.) {
     // discard previous points if its are before maximum point and with gap>1
        if(i<nmax ) {
          nn = 0; // nn--;
     // if point with gap after maximum - finish selection
        } else if(i>=nmax ) {
          break;
        }
      }
     // Far away from maximum
     //if(i-nmax > Int_t(5*tau))              break;
    }
    tn[nn] = Double_t(t[i]);
    yn[nn] = yd;
    //printf("i %i : nn %i : tn %6.2f : yn %6.2f \n", i, nn, tn[nn], yn[nn]); 
    nn++;
  }
  //printf(" nmax %i : n %i : nn %i i1 %i \n", nmax, n, nn, i1);
}

void AliPHOSFastAltroFit::FastFit(Double_t* t, Double_t* y, Int_t n, Double_t sig, Double_t tau,
Double_t &amp, Double_t &eamp, Double_t &t0, Double_t &et0, Double_t &chi2)
{
  // It is case of n=k=2 : fnn = x*x*exp(2 - 2*x)
  // Input: 
  //     n  - number of points 
  //   t[n] - array of time bins
  //   y[n] - array of amplitudes after pedestal subtractions;
  //   sig   - error of amplitude measurement (one value for all channels)
  //   tau   - filter time response (in timebin units)
  // Output:
  //       amp - amplitude at t0;
  //        t0 - time of max amplitude; 
  static Double_t xx; // t/tau
  static Double_t a, b, c;
  static Double_t f02, f12, f22;    // functions
  static Double_t f02d, f12d, f22d; // functions derivations

  chi2 = -1.;

  if(n<=0) {
    printf("<I> FastFit : n<=0 \n"); 
    return;
  }

  a = b = c = 0.0;
  for(Int_t i=0; i<n; i++){
    xx  = t[i]/tau;
    f02 = exp(-2.*xx);
    f12 = xx*f02;
    f22 = xx*f12;
    // Derivations
    f02d = -2.*f02;
    f12d = f02 - 2.*f12;
    f22d = 2.*(f12 - f22); 
    //
    a += f02d * y[i];
    b -= 2.*f12d * y[i];
    c += f22d * y[i];
  }
  Double_t t0_1=0.0, t0_2=0.0;
  Double_t amp_1=0.0, amp_2=0.0, chi2_1=0.0, chi2_2=0.0;
  if(QuadraticRoots(a,b,c, t0_1,t0_2)) {
    t0_1 *= tau;
    t0_2 *= tau;
    Amplitude(t,y,n, sig, tau, t0_1, amp_1, chi2_1);
    Amplitude(t,y,n, sig, tau, t0_2, amp_2, chi2_2);
    if(0) {
      printf(" t0_1 %f : t0_2 %f \n", t0_1, t0_2);
      printf(" amp_1 %f : amp_2 %f \n", amp_1, amp_2);
      printf(" chi2_1 %f : chi2_2 %f \n", chi2_1, chi2_2);
    }
    // t0 less on one tau with comparing with value from  "canonical equation"
    amp  = amp_1;
    t0   = t0_1;
    chi2 = chi2_1; 
    if(chi2_1 > chi2_2) {
      amp  = amp_2;
      t0   = t0_2; 
      chi2 = chi2_2; 
    }
    if(tau<3.) { // EMCAL case : small tau 
      t0 += -0.03; 
      Amplitude(t,y,n, sig, tau, t0, amp, chi2);
    }
    CalculateParsErrors(t, y, n, sig, tau, amp, t0, eamp, et0);

    // Fill1();
    
    // DrawFastFunction(amp, t0, fUtils->GetPedestalValue(), "1");
    //    DrawFastFunction(amp_1, t0_1, fUtils->GetPedestalValue(), "1");
    // DrawFastFunction(amp_2, t0_2, fUtils->GetPedestalValue(), "2");
  } else {
    chi2 = t0_1; // no roots, bad fit - negative chi2
  }
}

Bool_t AliPHOSFastAltroFit::QuadraticRoots(Double_t a, Double_t b, Double_t c, Double_t &x1, Double_t &x2)
{
  // Resolve quadratic equations a*x**2 + b*x + c
  //printf(" a %12.5e b %12.5e c %12.5e \n", a, b, c);
  static Double_t dtmp = 0.0, dtmpCut = -1.e-6;
  static Int_t ierr=0;
  dtmp = b*b - 4.*a*c;

  if(dtmp>=dtmpCut && dtmp<0.0) {
    if(ierr<5 || ierr%1000==0)
      printf("QuadraticRoots : %i small negative square : dtmp %12.5e \n", ierr++, dtmp);
    dtmp = 0.0;
  }
  if(dtmp>=0.0) {
    dtmp = sqrt(dtmp);
    x1   = (-b + dtmp) / (2.*a);
    x2   = (-b - dtmp) / (2.*a);

    //    printf(" x1 %f : x2 %f \n", x1, x2);
    return kTRUE;
  } else {
    x1 = dtmp;
    printf("QuadraticRoots : negative square : dtmp %12.5e \n", dtmp);
    return kFALSE;
  }
}

void AliPHOSFastAltroFit::Amplitude(Double_t* t,Double_t* y,Int_t n, Double_t sig, Double_t tau, 
Double_t t0, Double_t &amp, Double_t &chi2)
{  
  // Calculate parameters error too - Mar 24,09
  // sig is independent from points
  amp = 0.;
  Double_t x=0.0, f=0.0, den=0.0, f02;
  for(Int_t i=0; i<n; i++){
    x    = (t[i] - t0)/tau;
    f02  = exp(-2.*x);
    f    = x*x*f02;     
    amp += f*y[i];
    den += f*f;
  }
  if(den>0.0) amp /= den;
  //
  // chi2 calculation
  //
  Double_t dy=0.0;
  chi2=0.;
  for(Int_t i=0; i<n; i++){
    x    = (t[i] - t0)/tau;
    f02  = exp(-2.*x);
    f    = amp*x*x*f02;
    dy   = y[i]-f;
    chi2 += dy*dy;
    //    printf(" %i : y %f -> f %f : dy %f \n", i, y[i], f, dy); 
  }
  chi2 /= (sig*sig);
}

void AliPHOSFastAltroFit::CalculateParsErrors(Double_t* t, Double_t* /*y*/, Int_t n, Double_t sig, 
					      Double_t tau, Double_t &amp, Double_t &t0, 
					      Double_t &eamp, Double_t &et0)
{
  // Remember that fmax = exp(-n);
  // fmax_nk = (n/k)**n*exp(-n) => n=k=2 => exp(-n) = exp(-2.)
  static Double_t cc = exp(-2.);
  //   static Double_t cc = exp(-fN); // mean(N)~1.5 ??

  Double_t sumf2=0.0, sumfd2=0.0, x, f02, f12, f22, f22d;

  for(Int_t i=0; i<n; i++){
    x    = (t[i] - t0)/tau;
    f02  = exp(-2.*x);
    f12  = x*f02;
    f22  = x*f12;
    sumf2 += f22 * f22;
    //
    f22d = 2.*(f12 - f22); 
    sumfd2 += f22d * f22d;
  }
  et0  = (sig/amp)/sqrt(sumfd2);
  eamp = sig/sqrt(sumf2);

  amp  *= cc;
  eamp *= cc;
}

//
// Drawing
//
TCanvas* AliPHOSFastAltroFit::DrawFastFunction()
{
  // QA of fitting
  if(fNfit<=0) return 0; // no points

  static TCanvas *c = 0;
  if(c==0) {
    c =  new TCanvas("fastFun","fastFun",800,600);
  }

  c->cd();
  
  Double_t* eamp = new Double_t[fNfit];
  Double_t* et   = new Double_t[fNfit];

  for(Int_t i=0; i<fNfit; i++) {
    eamp[i] = fSig;
    et[i]   = 0.0;
  }

  TGraphErrors *gr = new TGraphErrors(fNfit, fTfit,fAmpfit, et,eamp);
  gr->Draw("Ap");
  gr->SetTitle(Form("Fast Fit : #chi^{2}/ndf = %8.2f / %i", GetChi2(), GetNDF()));
  gr->GetHistogram()->SetXTitle(" time bin ");
  gr->GetHistogram()->SetYTitle(" amplitude ");

  if(fStdFun==0) {
     fStdFun = new TF1("stdFun", StdResponseFunction, 0., fTfit[fNfit-1]+2., 5);    
     fStdFun->SetParNames("amp","t0","tau","N","ped");
  }
  fStdFun->SetParameter(0, GetEnergy());
  fStdFun->SetParameter(1, GetTime() + GetTau());
  fStdFun->SetParameter(2, GetTau());  // 
  fStdFun->SetParameter(3, GetN());    // 2
  fStdFun->SetParameter(4, 0.);  // 
  
  fStdFun->SetLineColor(kBlue);
  fStdFun->SetLineWidth(1);

  fStdFun->Draw("same");

  delete [] eamp;
  delete [] et;

  c->Update();

  return c;
}

Double_t AliPHOSFastAltroFit::StdResponseFunction(Double_t *x, Double_t *par)
{
  // Standard Response Function : 
  // look to Double_t AliEMCALRawUtils::RawResponseFunction(Double_t *x, Double_t *par)
  // Using for drawing only.
  // 
  // Shape of the electronics raw reponse:
  // It is a semi-gaussian, 2nd order Gamma function (n=2) of the general form
  //
  // t' = (t - t0 + tau) / tau
  // F = A * t**N * exp( N * ( 1 - t) )   for t >= 0
  // F = 0                                for t < 0 
  //
  // parameters:
  // A:   par[0]   // Amplitude = peak value
  // t0:  par[1]
  // tau: par[2]
  // N:   par[3]
  // ped: par[4]
  //
  static Double_t signal , tau, n, ped, xx;

  tau = par[2];
  n   = par[3];
  ped = par[4];
  xx = ( x[0] - par[1] + tau ) / tau ;

  if (xx <= 0) 
    signal = ped ;  
  else {  
    signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ; 
  }
  return signal ;  
}
Commit	Line	Data
0b272b9b	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 enxxexp(-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
8d0626c1	37	#include <TH1.h>
	38	#include <TF1.h>
	39	#include <TCanvas.h>
	40	#include <TGraphErrors.h>
	41	#include <TMath.h>
	42
0b272b9b	43	#include <math.h>
	44
	45	ClassImp(AliPHOSFastAltroFit)
	46
	47	//____________________________________________________________________________
	48	AliPHOSFastAltroFit:: AliPHOSFastAltroFit()
	49	: TNamed(),
8d0626c1	50	fSig(0),fTau(0),fN(0),fPed(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0)
8d0626c1	51	,fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
0b272b9b	52	{
	53	}
	54
	55	//____________________________________________________________________________
	56	AliPHOSFastAltroFit::AliPHOSFastAltroFit(const char* name, const char* title, const Double_t tau)
	57	: TNamed(name, title),
8d0626c1	58	fSig(0),fTau(tau),fN(2),fPed(0), fAmp(0),fAmpErr(0),fT0(0),fT0Err(0),fChi2(0.),fNDF(0)
8d0626c1	59	,fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
0b272b9b	60	{
	61	if(strlen(name)==0) SetName("FastAltroFit");
	62	}
	63
8d0626c1	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
0b272b9b	72	//____________________________________________________________________________
	73	AliPHOSFastAltroFit::~AliPHOSFastAltroFit()
	74	{
8d0626c1	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;
0b272b9b	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	{
8d0626c1	111	Reset();
8d0626c1	112
0b272b9b	113	fSig = sig;
0b272b9b	114	fTau = tau;
8d0626c1	115	fPed = ped;
	116
	117	if(fTfit) delete [] fTfit;
	118	if(fAmpfit) delete [] fAmpfit;
0b272b9b	119
8d0626c1	120	fNfit = 0;
	121	fTfit = new Double_t[n];
	122	fAmpfit = new Double_t[n];
0b272b9b	123
8d0626c1	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]);
0b272b9b	128
8d0626c1	129	if(fNfit>=2) {
8d0626c1	130	FastFit(fTfit,fAmpfit,fNfit,sig,tau, fAmp,fAmpErr, fT0,fT0Err,fChi2);
0b272b9b	131
8d0626c1	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	}
0b272b9b	140
8d0626c1	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;
0b272b9b	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	{
8d0626c1	172	static Double_t yMinUnderPed=2.; // should be tune
0b272b9b	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);
8d0626c1	181	//i1 = 0;
0b272b9b	182	i1 = i1<0?0:i1;
8d0626c1	183	Int_t i2 = n;
0b272b9b	184
0b272b9b	185	nn = 0;
8d0626c1	186	Double_t yd=0.0, tdiff=0.0;;
0b272b9b	187	for(Int_t i=i1; i<i2; i++) {
8d0626c1	188	if(ped>0.0) {
	189	yd = Double_t(y[i]) - ped;
	190	} else {
	191	yd = Double_t(y[i]);
0b272b9b	192	}
8d0626c1	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++;
0b272b9b	214	}
8d0626c1	215	//printf(" nmax %i : n %i : nn %i i1 %i \n", nmax, n, nn, i1);
0b272b9b	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 = xxexp(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;
8d0626c1	230	// t0 - time of max amplitude;
0b272b9b	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
8d0626c1	235
3312b1e5	236	chi2 = -1.;
3312b1e5	237
0b272b9b	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	}
8d0626c1	279	if(tau<3.) { // EMCAL case : small tau
	280	t0 += -0.03;
	281	Amplitude(t,y,n, sig, tau, t0, amp, chi2);
	282	}
0b272b9b	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 {
3312b1e5	291	chi2 = t0_1; // no roots, bad fit - negative chi2
0b272b9b	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 ax2 + bx + c
	298	//printf(" a %12.5e b %12.5e c %12.5e \n", a, b, c);
3312b1e5	299	static Double_t dtmp = 0.0, dtmpCut = -1.e-6;
3312b1e5	300	static Int_t ierr=0;
0b272b9b	301	dtmp = bb - 4.a*c;
0b272b9b	302
3312b1e5	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);
0b272b9b	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 {
3312b1e5	316	x1 = dtmp;
3312b1e5	317	printf("QuadraticRoots : negative square : dtmp %12.5e \n", dtmp);
0b272b9b	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 = xxf02;
	333	amp += f*y[i];
	334	den += f*f;
	335	}
8d0626c1	336	if(den>0.0) amp /= den;
0b272b9b	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 = ampxx*f02;
	346	dy = y[i]-f;
	347	chi2 += dy*dy;
8d0626c1	348	// printf(" %i : y %f -> f %f : dy %f \n", i, y[i], f, dy);
0b272b9b	349	}
	350	chi2 /= (sig*sig);
	351	}
	352
8d0626c1	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)
0b272b9b	356	{
3312b1e5	357	// Remember that fmax = exp(-n);
0b272b9b	358	// fmax_nk = (n/k)*nexp(-n) => n=k=2 => exp(-n) = exp(-2.)
0b272b9b	359	static Double_t cc = exp(-2.);
8d0626c1	360	// static Double_t cc = exp(-fN); // mean(N)~1.5 ??
0b272b9b	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;
3312b1e5	366	f02 = exp(-2.*x);
0b272b9b	367	f12 = x*f02;
	368	f22 = x*f12;
	369	sumf2 += f22 * f22;
	370	//
	371	f22d = 2.*(f12 - f22);
	372	sumfd2 += f22d * f22d;
	373	}
3312b1e5	374	et0 = (sig/amp)/sqrt(sumfd2);
3312b1e5	375	eamp = sig/sqrt(sumf2);
0b272b9b	376
3312b1e5	377	amp *= cc;
3312b1e5	378	eamp *= cc;
0b272b9b	379	}
8d0626c1	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	}