[u/mrichter/AliRoot.git] / RAW / AliCaloFastAltroFitv0.cxx

/**************************************************************************
 * Copyright(c) 1998-2010 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.
//  The main goal is fast estimation of amplitude and t0.
//

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

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

#include <math.h>

ClassImp(AliCaloFastAltroFitv0)

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

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

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

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

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

void AliCaloFastAltroFitv0::FastFit(Int_t* t, Int_t* y, Int_t nPoints, Double_t sig, Double_t tau, 
				  Double_t n, Double_t ped, Double_t tMax)
{
  // n 2 here and unused
  n = 2.;
  Reset();

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

  Int_t ii=0;
  CutRightPart(t,y,nPoints,tMax, ii);
  nPoints = ii;

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

  
  DeductPedestal(t,y,nPoints,  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;
      fNoFit++;
    }
  } else if(fNfit==1){
    Reset(); // What to do here => fT0 = fTfit[0]; fAmp = fAmpFit[0] ??
  } else {
    Reset();
  }
}

//____________________________________________________________________________
void AliCaloFastAltroFitv0::FastFit(TH1F* h, Double_t sig, Double_t tau, Double_t n,
Double_t ped, Double_t tMax)
{
  // Service method for convinience only
  // h - hist with altro response, could have empty bin
  // and center of bin could be different from bin number  
  Reset();

  if(h==0) return;
  Int_t nPoints = h->GetNbinsX();
  if(nPoints<2) return; // Sep 07, 09

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

  Int_t nPositive=0;
  for(Int_t i=0; i<nPoints; i++) {
    if(h->GetBinContent(i+1) > 0.0){ // Get only positive
      y[nPositive] = Int_t(h->GetBinContent(i+1));
      t[nPositive] = Int_t(h->GetBinCenter(i+1)+0.0001);
      nPositive++;
    }
  }

  if(nPositive >= 2) {
    FastFit(t,y,nPoints, sig,tau,n,ped, tMax);
  }
  if(fChi2<=0.0) fNoFit++;

  delete [] t;
  delete [] y;
}

void AliCaloFastAltroFitv0::Reset()
{
  // Reset variables
  fSig  = fTau = 0.0;
  fAmp  = fAmpErr = fT0 = fT0Err = 0.0;
  fChi2 = -1.;
  fNDF  = fNfit = 0;

  if(fTfit)   delete [] fTfit;
  if(fAmpfit) delete [] fAmpfit;
  fTfit = fAmpfit = 0;
}


void AliCaloFastAltroFitv0::GetFitResult(Double_t &amp,Double_t &eamp,Double_t &t0,Double_t &et0, 
Double_t &chi2, Int_t &ndf) const
{
  // Return results of fitting procedure
  amp  = fAmp;
  eamp = fAmpErr;
  t0   = fT0;
  et0  = fT0Err;
  chi2 = fChi2;
  ndf  = fNDF;
}

void AliCaloFastAltroFitv0::GetFittedPoints(Int_t &nfit, Double_t* ar[2]) const
{
  nfit  = fNfit;
  ar[0] = fTfit;
  ar[1] = fAmpfit;
}
//
// static functions
// 
void  AliCaloFastAltroFitv0::CutRightPart(Int_t *t,Int_t *y,Int_t nPoints,Double_t tMax, Int_t &ii)
{
  // Cut right part of altro sample : static function
  Int_t tt=0;
  for(Int_t i=0; i<nPoints; i++) {
    tt = t[i];
    if(tMax && tt <= Int_t(tMax)) {
      t[ii] = tt;
      y[ii] = y[i];
      ii++;
    }
  }
  if(0) printf(" ii %i -> ii %i : tMax %7.2f \n", nPoints, ii, tMax);
}

void AliCaloFastAltroFitv0::DeductPedestal(Int_t* t, Int_t* y, Int_t nPoints, Double_t tau, Double_t ped, 
  Double_t* tn, Double_t* yn, Int_t &nPointsOut)
{
  // Pedestal deduction if ped is positive : static function
  // Discard part od samle if it is not compact.
  static Double_t yMinUnderPed=2.; // should be tune
  Int_t ymax=0, nmax=0;
  for(Int_t i=0; i<nPoints; 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 = nPoints;

  nPointsOut = 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 && nPointsOut>0){
      tdiff = t[i] - tn[nPointsOut-1];
      //      printf(" i %i : nPointsOut %i : tdiff %6.2f : tn[nPointsOut] %6.2f \n", i,nPointsOut, tdiff, tn[nPointsOut-1]);
      if(tdiff>1.) {
     // discard previous points if its are before maximum point and with gap>1
        if(i<nmax ) {
          nPointsOut = 0; // nPointsOut--;
     // 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[nPointsOut] = Double_t(t[i]);
    yn[nPointsOut] = yd;
    //printf("i %i : nPointsOut %i : tn %6.2f : yn %6.2f \n", i, nPointsOut, tn[nPointsOut], yn[nPointsOut]); 
    nPointsOut++;
  }
  //printf(" nmax %i : nPointsIn %i :nPointsOut %i i1 %i \n", nmax, nPointsIn, nPointsOut, i1);
}

void AliCaloFastAltroFitv0::FastFit(const Double_t* t, const Double_t* y, const Int_t nPoints, 
                                    const Double_t sig, const Double_t tau,
                                    Double_t &amp, Double_t &eamp, Double_t &t0, Double_t &et0, Double_t &chi2)
{
  // Static function
  // It is case of n=k=2 : fnn = x*x*exp(2 - 2*x)
  // Input: 
  //nPoints  - number of points 
  //   t[]   - array of time bins
  //   y[]   - 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;
  //      eamp - error of amplitude; 
  //        t0 - time of max amplitude; 
  //       et0 - error of t0;
  //      chi2 - chi2
  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(nPoints<2) {
    printf(" AliCaloFastAltroFitv0::FastFit : nPoints<=%i \n", nPoints); 
    return;
  }

  a = b = c = 0.0;
  for(Int_t i=0; i<nPoints; 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 t01=0.0, t02=0.0;
  Double_t amp1=0.0, amp2=0.0, chi21=0.0, chi22=0.0;
  if(QuadraticRoots(a,b,c, t01,t02)) {
    t01 *= tau;
    t02 *= tau;
    Amplitude(t,y,nPoints, sig, tau, t01, amp1, chi21);
    Amplitude(t,y,nPoints, sig, tau, t02, amp2, chi22);
    if(0) {
      printf(" t01 %f : t02 %f \n", t01, t02);
      printf(" amp1 %f : amp2 %f \n", amp1, amp2);
      printf(" chi21 %f : chi22 %f \n", chi21, chi22);
    }
    // t0 less on one tau with comparing with value from  "canonical equation"
    amp  = amp1;
    t0   = t01;
    chi2 = chi21; 
    if(chi21 > chi22) {
      amp  = amp2;
      t0   = t02; 
      chi2 = chi22; 
    }
    if(tau<3.) { // EMCAL case : small tau 
      t0 += -0.03; // Discard bias in t0
      Amplitude(t,y,nPoints, sig, tau, t0, amp, chi2);
    }
    CalculateParsErrors(t, y, nPoints, sig, tau, amp, t0, eamp, et0);

    // Fill1();
    
    // DrawFastFunction(amp, t0, fUtils->GetPedestalValue(), "1");
    //    DrawFastFunction(amp1, t01, fUtils->GetPedestalValue(), "1");
    // DrawFastFunction(amp2, t02, fUtils->GetPedestalValue(), "2");
  } else {
    chi2 = t01; // no roots, bad fit - negative chi2
  }
}

Bool_t AliCaloFastAltroFitv0::QuadraticRoots(const Double_t a, const Double_t b, const 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 iWarning=0, iNoSolution=0;
  dtmp = b*b - 4.*a*c;

  if(dtmp>=dtmpCut && dtmp<0.0) {
    if(iWarning<5 || iWarning%1000==0)
      printf("<W> %i small neg. sq. : dtmp %12.5e \n", iWarning, dtmp);
    iWarning++;
    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;
    if(iNoSolution<5 || iNoSolution%1000==0)
      printf("<No solution> %i neg. sq. : dtmp %12.5e \n", iNoSolution, dtmp);
    iNoSolution++;
    return kFALSE;
  }
}

void AliCaloFastAltroFitv0::Amplitude(const Double_t* t,const Double_t* y,const Int_t nPoints, 
                                      const Double_t sig, const Double_t tau, const 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<nPoints; 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<nPoints; 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 AliCaloFastAltroFitv0::CalculateParsErrors(const Double_t* t, const Double_t* /*y*/, const Int_t nPoints, 
                                                const Double_t sig, const 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<nPoints; 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* AliCaloFastAltroFitv0::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 AliCaloFastAltroFitv0::StdResponseFunction(const Double_t *x, const Double_t *par)
{
  // Static 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
40df175f	1	/**************************************************************************
276d61fd	2	* Copyright(c) 1998-2010 ALICE Experiment at CERN, All rights reserved. *
40df175f	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	// 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)
276d61fd	51	,fNoFit(0),fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
40df175f	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)
276d61fd	60	,fNoFit(0),fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
40df175f	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)
276d61fd	69	,fNoFit(0),fNfit(0),fTfit(0),fAmpfit(0), fStdFun(0)
40df175f	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
40df175f	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
276d61fd	115	if(fChi2> 0.0) {
	116	fNDF = fNfit - 2;
	117	} else {
	118	fNDF = 0;
	119	fNoFit++;
	120	}
40df175f	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	{
276d61fd	132	// Service method for convinience only
bb09d1e6	133	// h - hist with altro response, could have empty bin
bb09d1e6	134	// and center of bin could be different from bin number
40df175f	135	Reset();
	136
	137	if(h==0) return;
	138	Int_t nPoints = h->GetNbinsX();
61f3c0f2	139	if(nPoints<2) return; // Sep 07, 09
40df175f	140
	141	Int_t* t = new Int_t[nPoints];
	142	Int_t* y = new Int_t[nPoints];
	143
bb09d1e6	144	Int_t nPositive=0;
40df175f	145	for(Int_t i=0; i<nPoints; i++) {
bb09d1e6	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);
	149	nPositive++;
	150	}
40df175f	151	}
40df175f	152
bb09d1e6	153	if(nPositive >= 2) {
40df175f	154	FastFit(t,y,nPoints, sig,tau,n,ped, tMax);
40df175f	155	}
276d61fd	156	if(fChi2<=0.0) fNoFit++;
40df175f	157
	158	delete [] t;
	159	delete [] y;
	160	}
	161
	162	void AliCaloFastAltroFitv0::Reset()
	163	{
	164	// Reset variables
	165	fSig = fTau = 0.0;
	166	fAmp = fAmpErr = fT0 = fT0Err = 0.0;
	167	fChi2 = -1.;
	168	fNDF = fNfit = 0;
61f3c0f2	169
	170	if(fTfit) delete [] fTfit;
	171	if(fAmpfit) delete [] fAmpfit;
40df175f	172	fTfit = fAmpfit = 0;
	173	}
	174
	175
	176	void AliCaloFastAltroFitv0::GetFitResult(Double_t &amp,Double_t &eamp,Double_t &t0,Double_t &et0,
	177	Double_t &chi2, Int_t &ndf) const
	178	{
	179	// Return results of fitting procedure
	180	amp = fAmp;
	181	eamp = fAmpErr;
	182	t0 = fT0;
	183	et0 = fT0Err;
	184	chi2 = fChi2;
	185	ndf = fNDF;
	186	}
	187
	188	void AliCaloFastAltroFitv0::GetFittedPoints(Int_t &nfit, Double_t* ar[2]) const
	189	{
	190	nfit = fNfit;
	191	ar[0] = fTfit;
	192	ar[1] = fAmpfit;
	193	}
	194	//
	195	// static functions
	196	//
	197	void AliCaloFastAltroFitv0::CutRightPart(Int_t t,Int_t y,Int_t nPoints,Double_t tMax, Int_t &ii)
	198	{
61f3c0f2	199	// Cut right part of altro sample : static function
40df175f	200	Int_t tt=0;
	201	for(Int_t i=0; i<nPoints; i++) {
	202	tt = t[i];
	203	if(tMax && tt <= Int_t(tMax)) {
	204	t[ii] = tt;
	205	y[ii] = y[i];
	206	ii++;
	207	}
	208	}
	209	if(0) printf(" ii %i -> ii %i : tMax %7.2f \n", nPoints, ii, tMax);
	210	}
	211
276d61fd	212	void AliCaloFastAltroFitv0::DeductPedestal(Int_t* t, Int_t* y, Int_t nPoints, Double_t tau, Double_t ped,
276d61fd	213	Double_t* tn, Double_t* yn, Int_t &nPointsOut)
40df175f	214	{
61f3c0f2	215	// Pedestal deduction if ped is positive : static function
40df175f	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;
276d61fd	219	for(Int_t i=0; i<nPoints; i++){
40df175f	220	if(y[i]>ymax) {
	221	ymax = y[i];
	222	nmax = i;
	223	}
	224	}
	225	Int_t i1 = nmax - Int_t(tau);
	226	//i1 = 0;
	227	i1 = i1<0?0:i1;
276d61fd	228	Int_t i2 = nPoints;
40df175f	229
276d61fd	230	nPointsOut = 0;
40df175f	231	Double_t yd=0.0, tdiff=0.0;;
	232	for(Int_t i=i1; i<i2; i++) {
	233	if(ped>0.0) {
	234	yd = Double_t(y[i]) - ped;
	235	} else {
	236	yd = Double_t(y[i]);
	237	}
	238	if(yd < yMinUnderPed) continue;
	239
276d61fd	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]);
40df175f	243	if(tdiff>1.) {
	244	// discard previous points if its are before maximum point and with gap>1
	245	if(i<nmax ) {
276d61fd	246	nPointsOut = 0; // nPointsOut--;
40df175f	247	// if point with gap after maximum - finish selection
	248	} else if(i>=nmax ) {
	249	break;
	250	}
	251	}
	252	// Far away from maximum
	253	//if(i-nmax > Int_t(5*tau)) break;
	254	}
276d61fd	255	tn[nPointsOut] = Double_t(t[i]);
	256	yn[nPointsOut] = yd;
	257	//printf("i %i : nPointsOut %i : tn %6.2f : yn %6.2f \n", i, nPointsOut, tn[nPointsOut], yn[nPointsOut]);
	258	nPointsOut++;
40df175f	259	}
276d61fd	260	//printf(" nmax %i : nPointsIn %i :nPointsOut %i i1 %i \n", nmax, nPointsIn, nPointsOut, i1);
40df175f	261	}
40df175f	262
276d61fd	263	void AliCaloFastAltroFitv0::FastFit(const Double_t* t, const Double_t* y, const Int_t nPoints,
40df175f	264	const Double_t sig, const Double_t tau,
	265	Double_t &amp, Double_t &eamp, Double_t &t0, Double_t &et0, Double_t &chi2)
	266	{
276d61fd	267	// Static function
40df175f	268	// It is case of n=k=2 : fnn = xxexp(2 - 2*x)
40df175f	269	// Input:
276d61fd	270	//nPoints - number of points
	271	// t[] - array of time bins
	272	// y[] - array of amplitudes after pedestal subtractions;
40df175f	273	// sig - error of amplitude measurement (one value for all channels)
	274	// tau - filter time response (in timebin units)
	275	// Output:
	276	// amp - amplitude at t0;
276d61fd	277	// eamp - error of amplitude;
40df175f	278	// t0 - time of max amplitude;
276d61fd	279	// et0 - error of t0;
276d61fd	280	// chi2 - chi2
40df175f	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
	285
	286	chi2 = -1.;
	287
276d61fd	288	if(nPoints<2) {
276d61fd	289	printf(" AliCaloFastAltroFitv0::FastFit : nPoints<=%i \n", nPoints);
40df175f	290	return;
	291	}
	292
	293	a = b = c = 0.0;
276d61fd	294	for(Int_t i=0; i<nPoints; i++){
40df175f	295	xx = t[i]/tau;
	296	f02 = exp(-2.*xx);
	297	f12 = xx*f02;
	298	f22 = xx*f12;
	299	// Derivations
	300	f02d = -2.*f02;
	301	f12d = f02 - 2.*f12;
	302	f22d = 2.*(f12 - f22);
	303	//
	304	a += f02d * y[i];
	305	b -= 2.f12d y[i];
	306	c += f22d * y[i];
	307	}
	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)) {
	311	t01 *= tau;
	312	t02 *= tau;
276d61fd	313	Amplitude(t,y,nPoints, sig, tau, t01, amp1, chi21);
276d61fd	314	Amplitude(t,y,nPoints, sig, tau, t02, amp2, chi22);
40df175f	315	if(0) {
	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);
	319	}
	320	// t0 less on one tau with comparing with value from "canonical equation"
	321	amp = amp1;
	322	t0 = t01;
	323	chi2 = chi21;
	324	if(chi21 > chi22) {
	325	amp = amp2;
	326	t0 = t02;
	327	chi2 = chi22;
	328	}
	329	if(tau<3.) { // EMCAL case : small tau
276d61fd	330	t0 += -0.03; // Discard bias in t0
276d61fd	331	Amplitude(t,y,nPoints, sig, tau, t0, amp, chi2);
40df175f	332	}
276d61fd	333	CalculateParsErrors(t, y, nPoints, sig, tau, amp, t0, eamp, et0);
40df175f	334
	335	// Fill1();
	336
	337	// DrawFastFunction(amp, t0, fUtils->GetPedestalValue(), "1");
	338	// DrawFastFunction(amp1, t01, fUtils->GetPedestalValue(), "1");
	339	// DrawFastFunction(amp2, t02, fUtils->GetPedestalValue(), "2");
	340	} else {
	341	chi2 = t01; // no roots, bad fit - negative chi2
	342	}
	343	}
	344
	345	Bool_t AliCaloFastAltroFitv0::QuadraticRoots(const Double_t a, const Double_t b, const Double_t c,
	346	Double_t &x1, Double_t &x2)
	347	{
	348	// Resolve quadratic equations ax2 + bx + 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;
276d61fd	351	static Int_t iWarning=0, iNoSolution=0;
40df175f	352	dtmp = bb - 4.a*c;
	353
	354	if(dtmp>=dtmpCut && dtmp<0.0) {
276d61fd	355	if(iWarning<5 \|\| iWarning%1000==0)
	356	printf("<W> %i small neg. sq. : dtmp %12.5e \n", iWarning, dtmp);
	357	iWarning++;
40df175f	358	dtmp = 0.0;
	359	}
	360	if(dtmp>=0.0) {
	361	dtmp = sqrt(dtmp);
	362	x1 = (-b + dtmp) / (2.*a);
	363	x2 = (-b - dtmp) / (2.*a);
	364
	365	// printf(" x1 %f : x2 %f \n", x1, x2);
	366	return kTRUE;
	367	} else {
	368	x1 = dtmp;
276d61fd	369	if(iNoSolution<5 \|\| iNoSolution%1000==0)
	370	printf("<No solution> %i neg. sq. : dtmp %12.5e \n", iNoSolution, dtmp);
	371	iNoSolution++;
40df175f	372	return kFALSE;
	373	}
	374	}
	375
276d61fd	376	void AliCaloFastAltroFitv0::Amplitude(const Double_t* t,const Double_t* y,const Int_t nPoints,
40df175f	377	const Double_t sig, const Double_t tau, const Double_t t0,
	378	Double_t &amp, Double_t &chi2)
	379	{
	380	// Calculate parameters error too - Mar 24,09
	381	// sig is independent from points
	382	amp = 0.;
	383	Double_t x=0.0, f=0.0, den=0.0, f02;
276d61fd	384	for(Int_t i=0; i<nPoints; i++){
40df175f	385	x = (t[i] - t0)/tau;
	386	f02 = exp(-2.*x);
	387	f = xxf02;
	388	amp += f*y[i];
	389	den += f*f;
	390	}
	391	if(den>0.0) amp /= den;
	392	//
	393	// chi2 calculation
	394	//
	395	Double_t dy=0.0;
	396	chi2=0.;
276d61fd	397	for(Int_t i=0; i<nPoints; i++){
40df175f	398	x = (t[i] - t0)/tau;
	399	f02 = exp(-2.*x);
	400	f = ampxx*f02;
	401	dy = y[i]-f;
	402	chi2 += dy*dy;
	403	// printf(" %i : y %f -> f %f : dy %f \n", i, y[i], f, dy);
	404	}
	405	chi2 /= (sig*sig);
	406	}
	407
276d61fd	408	void AliCaloFastAltroFitv0::CalculateParsErrors(const Double_t* t, const Double_t* /y/, const Int_t nPoints,
40df175f	409	const Double_t sig, const Double_t tau,
	410	Double_t &amp, Double_t &t0, Double_t &eamp, Double_t &et0)
	411	{
	412	// Remember that fmax = exp(-n);
	413	// fmax_nk = (n/k)*nexp(-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 ??
	416
	417	Double_t sumf2=0.0, sumfd2=0.0, x, f02, f12, f22, f22d;
	418
276d61fd	419	for(Int_t i=0; i<nPoints; i++){
40df175f	420	x = (t[i] - t0)/tau;
	421	f02 = exp(-2.*x);
	422	f12 = x*f02;
	423	f22 = x*f12;
	424	sumf2 += f22 * f22;
	425	//
	426	f22d = 2.*(f12 - f22);
	427	sumfd2 += f22d * f22d;
	428	}
	429	et0 = (sig/amp)/sqrt(sumfd2);
	430	eamp = sig/sqrt(sumf2);
	431
	432	amp *= cc;
	433	eamp *= cc;
	434	}
	435
	436	//
	437	// Drawing
	438	//
	439	TCanvas* AliCaloFastAltroFitv0::DrawFastFunction()
	440	{
	441	// QA of fitting
	442	if(fNfit<=0) return 0; // no points
	443
	444	static TCanvas *c = 0;
	445	if(c==0) {
	446	c = new TCanvas("fastFun","fastFun",800,600);
	447	}
	448
	449	c->cd();
	450
	451	Double_t* eamp = new Double_t[fNfit];
	452	Double_t* et = new Double_t[fNfit];
	453
	454	for(Int_t i=0; i<fNfit; i++) {
	455	eamp[i] = fSig;
	456	et[i] = 0.0;
	457	}
	458
	459	TGraphErrors *gr = new TGraphErrors(fNfit, fTfit,fAmpfit, et,eamp);
	460	gr->Draw("Ap");
	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 ");
	464
	465	if(fStdFun==0) {
	466	fStdFun = new TF1("stdFun", StdResponseFunction, 0., fTfit[fNfit-1]+2., 5);
	467	fStdFun->SetParNames("amp","t0","tau","N","ped");
	468	}
	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.); //
	474
	475	fStdFun->SetLineColor(kBlue);
	476	fStdFun->SetLineWidth(1);
	477
	478	fStdFun->Draw("same");
	479
	480	delete [] eamp;
	481	delete [] et;
	482
	483	c->Update();
484
485	return c;
486	}
487
488	Double_t AliCaloFastAltroFitv0::StdResponseFunction(const Double_t x, const Double_t par)
489	{
61f3c0f2	490	// Static Standard Response Function :
40df175f	491	// look to Double_t AliEMCALRawUtils::RawResponseFunction(Double_t x, Double_t par)
	492	// Using for drawing only.
	493	//
	494	// Shape of the electronics raw reponse:
	495	// It is a semi-gaussian, 2nd order Gamma function (n=2) of the general form
	496	//
	497	// t' = (t - t0 + tau) / tau
	498	// F = A * t*N exp( N * ( 1 - t) ) for t >= 0
	499	// F = 0 for t < 0
	500	//
	501	// parameters:
	502	// A: par[0] // Amplitude = peak value
	503	// t0: par[1]
	504	// tau: par[2]
	505	// N: par[3]
	506	// ped: par[4]
	507	//
	508	static Double_t signal , tau, n, ped, xx;
	509
	510	tau = par[2];
	511	n = par[3];
	512	ped = par[4];
	513	xx = ( x[0] - par[1] + tau ) / tau ;
	514
	515	if (xx <= 0)
	516	signal = ped ;
	517	else {
	518	signal = ped + par[0] * TMath::Power(xx , n) * TMath::Exp(n * (1 - xx )) ;
	519	}
	520	return signal ;
	521	}