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

/**************************************************************************
 * Copyright(c) 2007, 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.                  *
 **************************************************************************/


// This class extracts the signal parameters (energy, time, quality)
// from ALTRO samples. Energy is in ADC counts, time is in time bin units.
// If sample is not in saturation, a coarse algorithm is used (a-la AliPHOSRawFitterv0)
// If sample is in saturation, the unsaturated part of the sample is fit a-la AliPHOSRawFitterv1
// 
//     AliPHOSRawFitterv4 *fitterv4=new AliPHOSRawFitterv4();
//     fitterv4->SetChannelGeo(module,cellX,cellZ,caloFlag);
//     fitterv4->SetCalibData(fgCalibData) ;
//     fitterv4->Eval(sig,sigStart,sigLength);
//     Double_t amplitude = fitterv4.GetEnergy();
//     Double_t time      = fitterv4.GetTime();
//     Bool_t   isLowGain = fitterv4.GetCaloFlag()==0;

// Author: Dmitri Peressounko

// --- ROOT system ---
#include "TArrayI.h"
#include "TMath.h"
#include "TObject.h"
#include "TArrayD.h"
#include "TList.h"
#include "TMinuit.h"


//Used for debug
//#include "TROOT.h"
//#include "TH1.h"
//#include "TCanvas.h"
//#include "TPad.h"
//#include "TF1.h"


// --- AliRoot header files ---
#include "AliPHOSRawFitterv4.h"
#include "AliPHOSCalibData.h"
#include "AliLog.h"

ClassImp(AliPHOSRawFitterv4)

//-----------------------------------------------------------------------------
AliPHOSRawFitterv4::AliPHOSRawFitterv4():
  AliPHOSRawFitterv1(),
  fFitHighGain(0) 
{
  //Default constructor
}

//-----------------------------------------------------------------------------
AliPHOSRawFitterv4::~AliPHOSRawFitterv4()
{
}

//-----------------------------------------------------------------------------

Bool_t AliPHOSRawFitterv4::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
{
  // Calculate signal parameters (energy, time, quality) from array of samples
  // Energy is a maximum sample minus pedestal 9
  // Time is the first time bin
  // Signal overflows is there are at least 3 samples of the same amplitude above 900

  fOverflow= kFALSE ;
  fEnergy  = 0;
  if (fNBunches > 1) {
    fQuality = 1000;
    return kTRUE;
  }
  
  const Float_t kBaseLine   = 1.0;
  const Int_t   kPreSamples = 10;

  Float_t  pedMean   = 0;
  Float_t  pedRMS    = 0;
  Int_t    nPed      = 0;
  UShort_t maxSample = 0;
  Int_t    nMax      = 0;

  for (Int_t i=0; i<sigLength; i++) {
    if (i>sigLength-kPreSamples) { //inverse signal time order
      nPed++;
      pedMean += signal[i];
      pedRMS  += signal[i]*signal[i] ;
    }
    if(signal[i] > maxSample ){ maxSample = signal[i]; nMax=0;}
    if(signal[i] == maxSample) nMax++;

  }

  fEnergy = (Double_t)maxSample;
  if (maxSample > 850 && nMax > 2) fOverflow = kTRUE;
  
  Double_t pedestal = 0 ;
  if (fPedSubtract) {
    if (nPed > 0) {
      fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
      if(fPedestalRMS > 0.) 
	fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
      pedestal = (Double_t)(pedMean/nPed);
    }
    else
      return kFALSE;
  }
  else {
    //take pedestals from DB
    pedestal = (Double_t) fAmpOffset ;
    if (fCalibData) {
      Float_t truePed       = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
      Int_t   altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
      pedestal += truePed - altroSettings ;
    }
    else{
      AliDebug(2,Form("Pedestal and offset are not read from OCDB. Use 0 for their values.")) ;
    }
  }
  fEnergy-=pedestal ;
  if (fEnergy < kBaseLine) fEnergy = 0;
  
  
  if(fOverflow && ((fCaloFlag == 0) || fFitHighGain)){ //by default fit LowGain only
    TArrayI *samples = new TArrayI(sigLength); // array of sample amplitudes
    TArrayI *times   = new TArrayI(sigLength); // array of sample time stamps
    Bool_t result = kTRUE ;
    for (Int_t i=0; i<sigLength; i++) {
      samples->AddAt(signal[i]-pedestal,sigLength-i-1);
      times  ->AddAt(i ,i);
    }
    //Prepare fit parameters
    //Evaluate time
    Int_t iStart = 0;
    while(iStart<sigLength && samples->At(iStart) <kBaseLine) iStart++ ;
    if      (fCaloFlag == 0){ // Low gain
      fSampleParamsLow->AddAt(pedestal,4) ;
      fSampleParamsLow->AddAt(double(maxSample),5) ;
      fSampleParamsLow->AddAt(double(iStart),6) ;
    }
    else if (fCaloFlag == 1){ // High gain
      fSampleParamsHigh->AddAt(pedestal,4) ;
      fSampleParamsHigh->AddAt(double(maxSample),5) ;
      fSampleParamsHigh->AddAt(double(iStart),6) ;
    }
    result=EvalWithFitting(samples,times); 
    fToFit->Clear("nodelete") ;
    delete samples ;
    delete times ;
    
    return result; 
  
  }
  
  
  //Sample withour overflow - Evaluate time
  fTime = sigStart-sigLength-3; 
  const Int_t nLine= 6 ;        //Parameters of fitting
  const Float_t eMinTOF = 10. ; //Choosed from beam-test and cosmic analyis
  const Float_t kAmp=0.35 ;     //Result slightly depends on them, so no getters
  // Avoid too low peak:
  if(fEnergy < eMinTOF){
     return kTRUE;
  }
  // Find index posK (kLevel is a level of "timestamp" point Tk):
  Int_t posK =sigLength-1 ; //last point before crossing k-level
  Double_t levelK = pedestal + kAmp*fEnergy;
  while(signal[posK] <= levelK && posK>=0){
     posK-- ;
  }
  posK++ ;

  if(posK == 0 || posK==sigLength-1){
    return kTRUE; 
  }

  // Find crosing point by solving linear equation (least squares method)
  Int_t np = 0;
  Int_t iup=posK-1 ;
  Int_t idn=posK ;
  Double_t sx = 0., sy = 0., sxx = 0., sxy = 0.;
  Double_t x,y ;

  while(np<nLine){
    //point above crossing point
    if(iup>=0){
      x = sigLength-iup-1;
      y = signal[iup];
      sx += x;
      sy += y;
      sxx += (x*x);
      sxy += (x*y);
      np++ ;
      iup-- ;
    }
    //Point below crossing point
    if(idn<sigLength){
      if(signal[idn]<pedestal){
        idn=sigLength-1 ; //do not scan further
	idn++ ;
        continue ;
      }
      x = sigLength-idn-1;
      y = signal[idn];
      sx += x;
      sy += y;
      sxx += (x*x);
      sxy += (x*y);
      np++;
      idn++ ;
    }
    if(idn>=sigLength && iup<0){
      break ; //can not fit futher
    }
  }

  Double_t det = np*sxx - sx*sx;
  if(det == 0){
    return kTRUE;
  }
  Double_t c1 = (np*sxy - sx*sy)/det;  //slope
  Double_t c0 = 0;
  if (np>0)
    c0 = (sy-c1*sx)/np; //offset
  if(c1 == 0){
    return kTRUE;
  }

  // Find where the line cross kLevel:
  fTime += (levelK - c0)/c1-5. ; //5: mean offset between k-Level and start times
  return kTRUE;

}
//===================================================
Bool_t AliPHOSRawFitterv4::EvalWithFitting(TArrayI*samples, TArrayI * times){

  // if sample has reasonable mean and RMS, try to fit it with gamma2
  const Float_t sampleMaxHG=102.332 ; 
  const Float_t sampleMaxLG=277.196 ; 
  
  gMinuit->mncler();                     // Reset Minuit's list of paramters
  gMinuit->SetPrintLevel(-1) ;           // No Printout
  gMinuit->SetFCN(AliPHOSRawFitterv1::UnfoldingChiSquare) ;  
  // To set the address of the minimization function 
  
  fToFit->Clear("nodelete") ;
  Double_t b=0,bmin=0,bmax=0 ;
  if      (fCaloFlag == 0){ // Low gain
    b=fSampleParamsLow->At(2) ;
    bmin=0.5 ;
    bmax=10. ;
    fToFit->AddFirst((TObject*)fSampleParamsLow) ; 
  }
  else if (fCaloFlag == 1){ // High gain
    b=fSampleParamsHigh->At(2) ;
    bmin=0.05 ;
    bmax=0.4 ;
    fToFit->AddFirst((TObject*)fSampleParamsHigh) ; 
  }
  fToFit->AddLast((TObject*)samples) ;
  fToFit->AddLast((TObject*)times) ;
  

  gMinuit->SetObjectFit((TObject*)fToFit) ;         // To tranfer pointer to UnfoldingChiSquare
  Int_t ierflg ;
  gMinuit->mnparm(0, "t0",  1., 0.01, -50., 50., ierflg) ;
  if(ierflg != 0){
    //	  AliWarning(Form("Unable to set initial value for fit procedure : t0=%e\n",1.*tStart) ) ;
    fEnergy =   0. ;
    fTime   =-999. ;
    fQuality= 999. ;
    return kTRUE ; //will scan further
  }
  Double_t amp0=0; 
  if      (fCaloFlag == 0) // Low gain
    amp0 = fEnergy/sampleMaxLG;
  else if (fCaloFlag == 1) // High gain
    amp0 = fEnergy/sampleMaxHG;
  
  gMinuit->mnparm(1, "Energy", amp0 , 0.01*amp0, 0, 0, ierflg) ;
  if(ierflg != 0){
    //	  AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;
    fEnergy =   0. ;
    fTime   =-999. ;
    fQuality= 999. ;
    return kTRUE ; //will scan further
  }

  gMinuit->mnparm(2, "p2", b, 0.01*b, bmin, bmax, ierflg) ;
  if(ierflg != 0){                                         
    //        AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;  
    fEnergy =   0. ;
    fTime   =-999. ;
    fQuality= 999. ;
    return kTRUE ; //will scan further  
  }             
  
  Double_t p0 = 0.0001 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
  //  depends on it. 
  Double_t p1 = 1.0 ;
  Double_t p2 = 0.0 ;
  gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ;   // force TMinuit to reduce function calls  
  gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ;   // force TMinuit to use my gradient  
  //////	gMinuit->SetMaxIterations(100);
  gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ;  // No Warnings
  
  gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ;    // minimize 
  
  Double_t err=0.,t0err=0. ;
  Double_t t0=0.,efit=0. ;
  gMinuit->GetParameter(0,t0, t0err) ;
  gMinuit->GetParameter(1,efit, err) ;
  
  Double_t bfit=0., berr=0. ;
  gMinuit->GetParameter(2,bfit,berr) ;
  
   //Calculate total energy
  //this is parameterization of dependence of pulse height on parameter b
  if(fCaloFlag == 0) // Low gain
    fEnergy = efit*(99.54910 + 78.65038*bfit) ;
  else if(fCaloFlag == 1) // High gain
    fEnergy=efit*(80.33109 + 128.6433*bfit) ;
  
  if(fEnergy < 0. || fEnergy > 10000.){
    //set energy to previously found max
    fTime   =-999.;
    fQuality= 999 ;
    fEnergy=0. ;
    return kTRUE;
  }                                                                             
  
  //evaluate fit quality
  Double_t fmin=0.,fedm=0.,errdef=0. ;
  Int_t npari,nparx,istat;
  gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ;
  fQuality = fmin/samples->GetSize() ;
  //compare quality with some parameterization
  if      (fCaloFlag == 0) // Low gain
    fQuality /= 2.00 + 0.0020*fEnergy ;
  else if (fCaloFlag == 1) // High gain
    fQuality /= 0.75 + 0.0025*fEnergy ;
  
  fTime  += t0 - 4.024*bfit ;

  if(fQuality==0.){//no points to fit)
    fTime   =-999.;
    fQuality= 1999 ;
    fEnergy=0. ;  
  }
  
/*  
  if(1){
    TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
    if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
    h->Reset() ;
    for (Int_t i=0; i<samples->GetSize(); i++) {
        h->SetBinContent(i+1,float(samples->At(i))) ;
    }
    TF1 * fffit = new TF1("fffit","[0]*(abs(x-[1])^[3]*exp(-[2]*(x-[1]))+[4]*(x-[1])*(x-[1])*exp(-[5]*(x-[1])))",0.,60.) ;
    TArrayD * params=(TArrayD*)fToFit->At(0) ; 
    Double_t n=params->At(0) ;
    Double_t alpha=params->At(1) ;
    Double_t beta=params->At(3) ;
    fffit->SetParameters(efit,t0,alpha,n,bfit,beta) ;
    fffit->SetLineColor(2) ;
    TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
    if(!can){
        can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
        can->SetFillColor(0) ;
        can->SetFillStyle(0) ;
        can->Range(0,0,1,1);
        can->SetBorderSize(0);
    }
    can->cd() ;
  
//      TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
      TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.001,0.999,0.999);
      spectrum_1->Draw();
      spectrum_1->cd();
      spectrum_1->Range(0,0,1,1);
      spectrum_1->SetFillColor(0);
      spectrum_1->SetFillStyle(4000);
      spectrum_1->SetBorderSize(1);
      spectrum_1->SetBottomMargin(0.012);
      spectrum_1->SetTopMargin(0.03);
      spectrum_1->SetLeftMargin(0.10);
      spectrum_1->SetRightMargin(0.05);

      char title[155] ;
      snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
      h->SetTitle(title) ;
//      h->Fit(fffit,"","",0.,51.) ;
      h->Draw() ;
      fffit->Draw("same") ;


//       can->cd() ;
//       TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
//       spectrum_2->SetFillColor(0) ;
//       spectrum_2->SetFillStyle(0) ;
//       spectrum_2->SetGridy() ;
//       spectrum_2->Draw();
//       spectrum_2->Range(0,0,1,1);
//       spectrum_2->SetFillColor(0);
//       spectrum_2->SetBorderSize(1);
//       spectrum_2->SetTopMargin(0.01);
//       spectrum_2->SetBottomMargin(0.25);
//       spectrum_2->SetLeftMargin(0.10);
//       spectrum_2->SetRightMargin(0.05);
//       spectrum_2->cd() ;
// 
//       TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
//       if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
//       hd->Reset() ;
//       for (Int_t i=0; i<samples->GetSize(); i++) {
//         hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples->At(i)-fffit->Eval(i)))) ;
//       }
//       hd->Draw();

      can->Update() ;
      printf("Press <enter> to continue\n") ;
      getchar();


      delete fffit ;
      delete spectrum_1 ;
//      delete spectrum_2 ;
    }
*/
  
  return kTRUE;

}
Commit	Line	Data
de0cfd46	1	/**************************************************************************
	2	* Copyright(c) 2007, 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
	17	// This class extracts the signal parameters (energy, time, quality)
	18	// from ALTRO samples. Energy is in ADC counts, time is in time bin units.
	19	// If sample is not in saturation, a coarse algorithm is used (a-la AliPHOSRawFitterv0)
	20	// If sample is in saturation, the unsaturated part of the sample is fit a-la AliPHOSRawFitterv1
	21	//
	22	// AliPHOSRawFitterv4 *fitterv4=new AliPHOSRawFitterv4();
	23	// fitterv4->SetChannelGeo(module,cellX,cellZ,caloFlag);
	24	// fitterv4->SetCalibData(fgCalibData) ;
	25	// fitterv4->Eval(sig,sigStart,sigLength);
	26	// Double_t amplitude = fitterv4.GetEnergy();
	27	// Double_t time = fitterv4.GetTime();
	28	// Bool_t isLowGain = fitterv4.GetCaloFlag()==0;
	29
	30	// Author: Dmitri Peressounko
	31
	32	// --- ROOT system ---
	33	#include "TArrayI.h"
	34	#include "TMath.h"
	35	#include "TObject.h"
	36	#include "TArrayD.h"
	37	#include "TList.h"
	38	#include "TMinuit.h"
	39
	40
	41	//Used for debug
	42	//#include "TROOT.h"
	43	//#include "TH1.h"
	44	//#include "TCanvas.h"
	45	//#include "TPad.h"
	46	//#include "TF1.h"
	47
	48
	49
	50	// --- AliRoot header files ---
	51	#include "AliPHOSRawFitterv4.h"
	52	#include "AliPHOSCalibData.h"
	53	#include "AliLog.h"
	54
	55	ClassImp(AliPHOSRawFitterv4)
	56
	57	//-----------------------------------------------------------------------------
	58	AliPHOSRawFitterv4::AliPHOSRawFitterv4():
	59	AliPHOSRawFitterv1(),
	60	fFitHighGain(0)
	61	{
	62	//Default constructor
	63	}
	64
65	//-----------------------------------------------------------------------------
66	AliPHOSRawFitterv4::~AliPHOSRawFitterv4()
67	{
68	}
69
de0cfd46	70	//-----------------------------------------------------------------------------
	71
	72	Bool_t AliPHOSRawFitterv4::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
	73	{
	74	// Calculate signal parameters (energy, time, quality) from array of samples
	75	// Energy is a maximum sample minus pedestal 9
	76	// Time is the first time bin
	77	// Signal overflows is there are at least 3 samples of the same amplitude above 900
	78
	79	fOverflow= kFALSE ;
	80	fEnergy = 0;
	81	if (fNBunches > 1) {
	82	fQuality = 1000;
	83	return kTRUE;
	84	}
	85
	86	const Float_t kBaseLine = 1.0;
	87	const Int_t kPreSamples = 10;
	88
	89	Float_t pedMean = 0;
	90	Float_t pedRMS = 0;
	91	Int_t nPed = 0;
	92	UShort_t maxSample = 0;
	93	Int_t nMax = 0;
	94
	95	for (Int_t i=0; i<sigLength; i++) {
	96	if (i>sigLength-kPreSamples) { //inverse signal time order
	97	nPed++;
	98	pedMean += signal[i];
	99	pedRMS += signal[i]*signal[i] ;
	100	}
	101	if(signal[i] > maxSample ){ maxSample = signal[i]; nMax=0;}
	102	if(signal[i] == maxSample) nMax++;
	103
	104	}
	105
	106	fEnergy = (Double_t)maxSample;
	107	if (maxSample > 850 && nMax > 2) fOverflow = kTRUE;
	108
	109	Double_t pedestal = 0 ;
	110	if (fPedSubtract) {
	111	if (nPed > 0) {
	112	fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
	113	if(fPedestalRMS > 0.)
	114	fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
	115	pedestal = (Double_t)(pedMean/nPed);
	116	}
	117	else
	118	return kFALSE;
	119	}
	120	else {
	121	//take pedestals from DB
	122	pedestal = (Double_t) fAmpOffset ;
	123	if (fCalibData) {
	124	Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
	125	Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
	126	pedestal += truePed - altroSettings ;
	127	}
	128	else{
	129	AliDebug(2,Form("Pedestal and offset are not read from OCDB. Use 0 for their values.")) ;
	130	}
	131	}
	132	fEnergy-=pedestal ;
	133	if (fEnergy < kBaseLine) fEnergy = 0;
134
135
136	if(fOverflow && ((fCaloFlag == 0) \|\| fFitHighGain)){ //by default fit LowGain only
137	TArrayI *samples = new TArrayI(sigLength); // array of sample amplitudes
138	TArrayI *times = new TArrayI(sigLength); // array of sample time stamps
139	Bool_t result = kTRUE ;
140	for (Int_t i=0; i<sigLength; i++) {
141	samples->AddAt(signal[i]-pedestal,sigLength-i-1);
142	times ->AddAt(i ,i);
143	}
144	//Prepare fit parameters
145	//Evaluate time
146	Int_t iStart = 0;
147	while(iStart<sigLength && samples->At(iStart) <kBaseLine) iStart++ ;
148	if (fCaloFlag == 0){ // Low gain
149	fSampleParamsLow->AddAt(pedestal,4) ;
150	fSampleParamsLow->AddAt(double(maxSample),5) ;
151	fSampleParamsLow->AddAt(double(iStart),6) ;
de0cfd46	152	}
	153	else if (fCaloFlag == 1){ // High gain
	154	fSampleParamsHigh->AddAt(pedestal,4) ;
	155	fSampleParamsHigh->AddAt(double(maxSample),5) ;
	156	fSampleParamsHigh->AddAt(double(iStart),6) ;
de0cfd46	157	}
de0cfd46	158	result=EvalWithFitting(samples,times);
06532a37	159	fToFit->Clear("nodelete") ;
de0cfd46	160	delete samples ;
	161	delete times ;
	162
	163	return result;
	164
	165	}
	166
	167
	168	//Sample withour overflow - Evaluate time
	169	fTime = sigStart-sigLength-3;
	170	const Int_t nLine= 6 ; //Parameters of fitting
	171	const Float_t eMinTOF = 10. ; //Choosed from beam-test and cosmic analyis
	172	const Float_t kAmp=0.35 ; //Result slightly depends on them, so no getters
	173	// Avoid too low peak:
	174	if(fEnergy < eMinTOF){
	175	return kTRUE;
	176	}
	177	// Find index posK (kLevel is a level of "timestamp" point Tk):
	178	Int_t posK =sigLength-1 ; //last point before crossing k-level
	179	Double_t levelK = pedestal + kAmp*fEnergy;
	180	while(signal[posK] <= levelK && posK>=0){
	181	posK-- ;
	182	}
	183	posK++ ;
	184
	185	if(posK == 0 \|\| posK==sigLength-1){
	186	return kTRUE;
	187	}
	188
	189	// Find crosing point by solving linear equation (least squares method)
	190	Int_t np = 0;
	191	Int_t iup=posK-1 ;
	192	Int_t idn=posK ;
	193	Double_t sx = 0., sy = 0., sxx = 0., sxy = 0.;
	194	Double_t x,y ;
	195
	196	while(np<nLine){
	197	//point above crossing point
	198	if(iup>=0){
	199	x = sigLength-iup-1;
	200	y = signal[iup];
	201	sx += x;
	202	sy += y;
	203	sxx += (x*x);
	204	sxy += (x*y);
	205	np++ ;
	206	iup-- ;
	207	}
	208	//Point below crossing point
	209	if(idn<sigLength){
	210	if(signal[idn]<pedestal){
	211	idn=sigLength-1 ; //do not scan further
	212	idn++ ;
	213	continue ;
	214	}
	215	x = sigLength-idn-1;
	216	y = signal[idn];
	217	sx += x;
	218	sy += y;
	219	sxx += (x*x);
	220	sxy += (x*y);
	221	np++;
	222	idn++ ;
	223	}
224	if(idn>=sigLength && iup<0){
225	break ; //can not fit futher
226	}
227	}
228
229	Double_t det = npsxx - sxsx;
230	if(det == 0){
231	return kTRUE;
232	}
233	Double_t c1 = (npsxy - sxsy)/det; //slope
eaf6dbb0	234	Double_t c0 = 0;
	235	if (np>0)
	236	c0 = (sy-c1*sx)/np; //offset
de0cfd46	237	if(c1 == 0){
	238	return kTRUE;
	239	}
	240
	241	// Find where the line cross kLevel:
	242	fTime += (levelK - c0)/c1-5. ; //5: mean offset between k-Level and start times
	243	return kTRUE;
	244
	245	}
	246	//===================================================
	247	Bool_t AliPHOSRawFitterv4::EvalWithFitting(TArrayIsamples, TArrayI times){
	248
	249	// if sample has reasonable mean and RMS, try to fit it with gamma2
	250	const Float_t sampleMaxHG=102.332 ;
	251	const Float_t sampleMaxLG=277.196 ;
	252
	253	gMinuit->mncler(); // Reset Minuit's list of paramters
	254	gMinuit->SetPrintLevel(-1) ; // No Printout
	255	gMinuit->SetFCN(AliPHOSRawFitterv1::UnfoldingChiSquare) ;
	256	// To set the address of the minimization function
	257
	258	fToFit->Clear("nodelete") ;
	259	Double_t b=0,bmin=0,bmax=0 ;
	260	if (fCaloFlag == 0){ // Low gain
	261	b=fSampleParamsLow->At(2) ;
	262	bmin=0.5 ;
	263	bmax=10. ;
	264	fToFit->AddFirst((TObject*)fSampleParamsLow) ;
	265	}
	266	else if (fCaloFlag == 1){ // High gain
	267	b=fSampleParamsHigh->At(2) ;
	268	bmin=0.05 ;
	269	bmax=0.4 ;
	270	fToFit->AddFirst((TObject*)fSampleParamsHigh) ;
	271	}
	272	fToFit->AddLast((TObject*)samples) ;
	273	fToFit->AddLast((TObject*)times) ;
	274
	275
	276	gMinuit->SetObjectFit((TObject*)fToFit) ; // To tranfer pointer to UnfoldingChiSquare
	277	Int_t ierflg ;
	278	gMinuit->mnparm(0, "t0", 1., 0.01, -50., 50., ierflg) ;
	279	if(ierflg != 0){
	280	// AliWarning(Form("Unable to set initial value for fit procedure : t0=%e\n",1.*tStart) ) ;
	281	fEnergy = 0. ;
	282	fTime =-999. ;
	283	fQuality= 999. ;
	284	return kTRUE ; //will scan further
	285	}
	286	Double_t amp0=0;
	287	if (fCaloFlag == 0) // Low gain
	288	amp0 = fEnergy/sampleMaxLG;
	289	else if (fCaloFlag == 1) // High gain
	290	amp0 = fEnergy/sampleMaxHG;
	291
	292	gMinuit->mnparm(1, "Energy", amp0 , 0.01*amp0, 0, 0, ierflg) ;
	293	if(ierflg != 0){
	294	// AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;
	295	fEnergy = 0. ;
	296	fTime =-999. ;
	297	fQuality= 999. ;
	298	return kTRUE ; //will scan further
	299	}
	300
301	gMinuit->mnparm(2, "p2", b, 0.01*b, bmin, bmax, ierflg) ;
302	if(ierflg != 0){
303	// AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ;
304	fEnergy = 0. ;
305	fTime =-999. ;
306	fQuality= 999. ;
307	return kTRUE ; //will scan further
308	}
309
310	Double_t p0 = 0.0001 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly
311	// depends on it.
312	Double_t p1 = 1.0 ;
313	Double_t p2 = 0.0 ;
314	gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls
315	gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient
316	////// gMinuit->SetMaxIterations(100);
317	gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings
318
319	gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
320
321	Double_t err=0.,t0err=0. ;
322	Double_t t0=0.,efit=0. ;
323	gMinuit->GetParameter(0,t0, t0err) ;
324	gMinuit->GetParameter(1,efit, err) ;
325
326	Double_t bfit=0., berr=0. ;
327	gMinuit->GetParameter(2,bfit,berr) ;
328
329	//Calculate total energy
330	//this is parameterization of dependence of pulse height on parameter b
331	if(fCaloFlag == 0) // Low gain
332	fEnergy = efit(99.54910 + 78.65038bfit) ;
333	else if(fCaloFlag == 1) // High gain
334	fEnergy=efit(80.33109 + 128.6433bfit) ;
335
336	if(fEnergy < 0. \|\| fEnergy > 10000.){
337	//set energy to previously found max
338	fTime =-999.;
339	fQuality= 999 ;
340	fEnergy=0. ;
341	return kTRUE;
342	}
343
344	//evaluate fit quality
345	Double_t fmin=0.,fedm=0.,errdef=0. ;
346	Int_t npari,nparx,istat;
347	gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ;
348	fQuality = fmin/samples->GetSize() ;
349	//compare quality with some parameterization
350	if (fCaloFlag == 0) // Low gain
351	fQuality /= 2.00 + 0.0020*fEnergy ;
352	else if (fCaloFlag == 1) // High gain
353	fQuality /= 0.75 + 0.0025*fEnergy ;
354
355	fTime += t0 - 4.024*bfit ;
356
357	if(fQuality==0.){//no points to fit)
358	fTime =-999.;
359	fQuality= 1999 ;
360	fEnergy=0. ;
361	}
362
363	/*
364	if(1){
365	TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
366	if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
367	h->Reset() ;
368	for (Int_t i=0; i<samples->GetSize(); i++) {
369	h->SetBinContent(i+1,float(samples->At(i))) ;
370	}
371	TF1 * fffit = new TF1("fffit","[0](abs(x-[1])^[3]exp(-[2](x-[1]))+[4](x-[1])(x-[1])exp(-[5]*(x-[1])))",0.,60.) ;
372	TArrayD * params=(TArrayD*)fToFit->At(0) ;
373	Double_t n=params->At(0) ;
374	Double_t alpha=params->At(1) ;
375	Double_t beta=params->At(3) ;
376	fffit->SetParameters(efit,t0,alpha,n,bfit,beta) ;
377	fffit->SetLineColor(2) ;
378	TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
379	if(!can){
380	can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
381	can->SetFillColor(0) ;
382	can->SetFillStyle(0) ;
383	can->Range(0,0,1,1);
384	can->SetBorderSize(0);
385	}
386	can->cd() ;
387
388	// TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
389	TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.001,0.999,0.999);
390	spectrum_1->Draw();
391	spectrum_1->cd();
392	spectrum_1->Range(0,0,1,1);
393	spectrum_1->SetFillColor(0);
394	spectrum_1->SetFillStyle(4000);
395	spectrum_1->SetBorderSize(1);
396	spectrum_1->SetBottomMargin(0.012);
397	spectrum_1->SetTopMargin(0.03);
398	spectrum_1->SetLeftMargin(0.10);
399	spectrum_1->SetRightMargin(0.05);
400
401	char title[155] ;
402	snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
403	h->SetTitle(title) ;
404	// h->Fit(fffit,"","",0.,51.) ;
405	h->Draw() ;
406	fffit->Draw("same") ;
407
408
409	// can->cd() ;
410	// TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
411	// spectrum_2->SetFillColor(0) ;
412	// spectrum_2->SetFillStyle(0) ;
413	// spectrum_2->SetGridy() ;
414	// spectrum_2->Draw();
415	// spectrum_2->Range(0,0,1,1);
416	// spectrum_2->SetFillColor(0);
417	// spectrum_2->SetBorderSize(1);
418	// spectrum_2->SetTopMargin(0.01);
419	// spectrum_2->SetBottomMargin(0.25);
420	// spectrum_2->SetLeftMargin(0.10);
421	// spectrum_2->SetRightMargin(0.05);
422	// spectrum_2->cd() ;
423	//
424	// TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
425	// if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
426	// hd->Reset() ;
427	// for (Int_t i=0; i<samples->GetSize(); i++) {
428	// hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples->At(i)-fffit->Eval(i)))) ;
429	// }
430	// hd->Draw();
431
432	can->Update() ;
433	printf("Press <enter> to continue\n") ;
434	getchar();
435
436
437	delete fffit ;
438	delete spectrum_1 ;
439	// delete spectrum_2 ;
440	}
441	*/
442
443	return kTRUE;
444
445	}