[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawFitterv3.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.                  *
 **************************************************************************/

/* $Id: $ */

// This class extracts the signal parameters (energy, time, quality)
// from ALTRO samples. Energy is in ADC counts, time is in time bin units.
// Class uses FastFitting algorithm to fit sample and extract time and Amplitude 
// and evaluate sample quality = (chi^2/NDF)/some parameterization providing 
// efficiency close to 100%
// 
// Typical use case:
//     AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
//     fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
//     fitter->SetCalibData(fgCalibData) ;
//     fitter->Eval(sig,sigStart,sigLength);
//     Double_t amplitude = fitter.GetEnergy();
//     Double_t time      = fitter.GetTime();
//     Bool_t   isLowGain = fitter.GetCaloFlag()==0;

// Author: Dmitri Peressounko (after A.Pavlinov - see RAW/AliCaloFastAltroFitv0.cxx)

// --- ROOT system ---
#include "TArrayI.h"
#include "TList.h"
#include "TMath.h"
#include "TH1I.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TFile.h"
#include "TROOT.h"

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

ClassImp(AliPHOSRawFitterv3)

//-----------------------------------------------------------------------------
AliPHOSRawFitterv3::AliPHOSRawFitterv3():
  AliPHOSRawFitterv0()
{
  //Default constructor.
}

//-----------------------------------------------------------------------------
AliPHOSRawFitterv3::~AliPHOSRawFitterv3()
{
  //Destructor.
}

//-----------------------------------------------------------------------------
AliPHOSRawFitterv3::AliPHOSRawFitterv3(const AliPHOSRawFitterv3 &phosFitter ):
  AliPHOSRawFitterv0(phosFitter) 
{
  //Copy constructor.
}

//-----------------------------------------------------------------------------
AliPHOSRawFitterv3& AliPHOSRawFitterv3::operator = (const AliPHOSRawFitterv3 & /*phosFitter*/)
{
  //Assignment operator.
  return *this;
}

//-----------------------------------------------------------------------------
Bool_t AliPHOSRawFitterv3::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
{
  //Extract an energy deposited in the crystal,
  //crystal' position (module,column,row),
  //time and gain (high or low).
  //First collects sample, then evaluates it and if it has
  //reasonable shape, fits it with Gamma2 function and extracts 
  //energy and time.

  if (fCaloFlag == 2 || fNBunches > 1) {
    fQuality = 1000;
    return kTRUE;
  }
  if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
    fQuality=2000;
    fEnergy=0 ;
    return kTRUE;
  }

  Float_t pedMean = 0;
  Float_t pedRMS  = 0;
  Int_t   nPed    = 0;
  const Float_t kBaseLine   = 1.0;
  const Int_t   kPreSamples = 10;
  

  //We tryed individual taus for each channel, but
  //this approach seems to be unstable. Much better results are obtaned with
  //fixed decay time for all channels.
  const Double_t tau=21. ;

  TArrayD *samples = new TArrayD(sigLength); // array of sample amplitudes
  TArrayD *times   = new TArrayD(sigLength); // array of sample time stamps
  for (Int_t i=0; i<sigLength; i++) {
    if (i<kPreSamples) {
      nPed++;
      pedMean += signal[i];
      pedRMS  += signal[i]*signal[i] ;
    }
    samples->AddAt(signal[i],sigLength-i-1);
    times  ->AddAt(i/tau ,i);
  }

  fEnergy = -111;
  fQuality= 999. ;
  const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
  Double_t pedestal = 0;

  if (fPedSubtract) {
    if (nPed > 0) {
      fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
      if(fPedestalRMS > 0.) 
	fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
      fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction
    }
    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{
      AliWarning(Form("Can not read data from OCDB")) ;
    }
    fEnergy-=pedestal ;
  }

  if (fEnergy < kBaseLine) fEnergy = 0;
  //Evaluate time
  Int_t iStart = 0;
  while(iStart<sigLength && samples->At(iStart)-pedestal <kBaseLine) iStart++ ;
  fTime = sigStart+iStart; 
  
  //calculate time and energy
  Int_t    maxBin=0 ;
  Int_t    maxAmp=0 ;
  Int_t    nMax = 0 ; //number of points in plato
  Double_t aMean =0. ;
  Double_t aRMS  =0. ;
  Double_t wts   =0 ;
  Int_t    tStart=0 ;

  for (Int_t i=0; i<sigLength; i++){
    if(signal[i] > pedestal){
      Double_t de = signal[i] - pedestal ;
      if(de > 1.) {
	aMean += de*i ;
	aRMS  += de*i*i ;
	wts   += de; 
      }
      if(de > 2 && tStart==0) 
	tStart = i ;
      if(signal[i] >  maxAmp){
        maxAmp = signal[i]; 
        nMax=0;
	maxBin = i ;
      }
      if(signal[i] == maxAmp)
        nMax++;
    }
  }

  if (maxBin==sigLength-1){//bad "rising" sample
    fEnergy =    0. ;
    fTime   = -999. ;
    fQuality=  999. ;
    return kTRUE ;
  }

  fEnergy=Double_t(maxAmp)-pedestal ;
  fOverflow =0 ;  //look for plato on the top of sample
  if (fEnergy>500 &&  //this is not fluctuation of soft sample
     nMax>2){ //and there is a plato
    fOverflow = kTRUE ;
  }
  
  if (wts > 0) {
    aMean /= wts; 
    aRMS   = aRMS/wts - aMean*aMean;
  }

  //do not take too small energies
  if (fEnergy < kBaseLine) 
    fEnergy = 0;
  
  //do not test quality of too soft samples
  if (fEnergy < maxEtoFit){
    fTime = tStart;
    if (aRMS < 2.) //sigle peak
      fQuality = 999. ;
    else
      fQuality =   0. ;
    return kTRUE ;
  }
      
  // if sample has reasonable mean and RMS, try to fit it with gamma2
  // First estimate t0
  Double_t a=0,b=0,c=0 ;
  for(Int_t i=0; i<sigLength; i++){
    if(samples->At(i)<pedestal)
      continue ;
    Double_t t= times->At(i) ;
    Double_t f02 = TMath::Exp(-2.*t);
    Double_t f12 = t*f02;
    Double_t f22 = t*f12;
    // Derivatives
    Double_t f02d = -2.*f02;
    Double_t f12d = f02 - 2.*f12;
    Double_t f22d = 2.*(f12 - f22);
    a += f02d * (samples->At(i)-pedestal) ;
    b -= f12d * (samples->At(i)-pedestal) ;
    c += f22d * (samples->At(i)-pedestal) ;
  }
  
  //Find 2 roots
  Double_t det = b*b - a*c;
  if(det>=1.e-6 && det<0.0) {
    det = 0.0; //rounding error
  }
  if(det<0.){ //Problem
    fQuality = 1500. ;
    if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
      printf(" det=%e \n",det) ;
      goto plot ;
    }
    return kTRUE ;
  }

  if(det>0.0 && a!=0.) {
    det = TMath::Sqrt(det);
    Double_t t1 = (-b + det) / a;
//    Double_t t2 = (-b - det) / a; //second root is wrong one
    Double_t amp1=0., den1=0. ;
    for(Int_t i=0; i<sigLength; i++){
       Double_t dt1 = times->At(i) - t1;
       Double_t f01 = dt1*dt1*TMath::Exp(-2.*dt1);
       amp1 += f01*(samples->At(i)-pedestal);
       den1 += f01*f01;
     }
     if(den1>0.0) amp1 /= den1;
     Double_t chi1=0.; // chi2=0. ;
     for(Int_t i=0; i<sigLength; i++){
       Double_t dt1 = times->At(i)- t1;
       Double_t dy1 = samples->At(i)-pedestal- amp1*dt1*dt1*TMath::Exp(-2.*dt1) ;
       chi1 += dy1*dy1;
     }
     fEnergy=amp1*TMath::Exp(-2.) ; ; 
     fTime=t1*tau ;
     fQuality=chi1/sigLength ;
  } 
  else { 
    Double_t t1 ;
    if(a!=0.)
      t1=-b/a ;
    else
      t1=-c/b ;
    Double_t amp=0.,den=0.; ;
    for(Int_t i=0; i<sigLength; i++){
       Double_t dt = times->At(i) - t1;
       Double_t f = dt*dt*TMath::Exp(-2.*dt);
       amp += f*samples->At(i);
       den += f*f;
     }
     if(den>0.0) amp /= den;
     // chi2 calculation
     fQuality=0.;
     for(Int_t i=0; i<sigLength; i++){
       Double_t t = times->At(i)- t1;
       Double_t dy = samples->At(i)- amp*t*t*TMath::Exp(-2.*t) ;
       fQuality += dy*dy;
     }
     fTime=t1*tau ;
     fEnergy = amp*TMath::Exp(-2.);
     fQuality/= sigLength ;
  }

  //Impose cut on quality
  fQuality/=2.+0.004*fEnergy*fEnergy ;

  //Draw corrupted samples
  if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
    plot:
    if(fQuality >1. && !fOverflow ){ //Draw only bad samples
      printf("Sample par: amp=%f,  t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
      TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
      if(!h) h = new TH1I("h","Samples",50,0.,50.) ;
      h->Reset() ;
      for (Int_t i=0; i<sigLength; i++) {
        h->SetBinContent(i,samples->At(i)) ;
      }
      TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
      fffit->SetParameters(pedestal,fEnergy,fTime,tau) ;
      fffit->SetLineColor(2) ;
      TCanvas * c = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
      if(!c){
        c = new TCanvas("cSamples","cSamples",10,10,400,400) ;
        c->SetFillColor(0) ;
        c->SetFillStyle(0) ;
        c->Range(0,0,1,1);
        c->SetBorderSize(0);
      }
      c->cd() ;
  
      TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
      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] ;
      sprintf(title,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
      h->SetTitle(title) ;
      h->Draw() ;
      fffit->Draw("same") ;

      sprintf(title,"mod%d_x%d_z%d_HG",fModule,fCellX,fCellZ) ;
      TFile fout("samples_bad.root","update") ;
      h->Write(title);
      fout.Close() ;

      c->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",50,0.,50.) ;
      hd->Reset() ;
      for (Int_t i=0; i<sigLength; i++) {
        hd->SetBinContent(i,TMath::Max(-1023.,TMath::Min(1023.,samples->At(i)-fffit->Eval(i)))) ;
      }
      hd->Draw();
 
      c->Update() ;
      printf("Press <enter> to continue\n") ;
      getchar();


      delete fffit ;
      delete spectrum_1 ;
      delete spectrum_2 ;
    }
  }
  
  delete samples ;
  delete times ;
  return kTRUE;
}
Commit	Line	Data
c4cc99c3	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	/* $Id: $ */
	17
	18	// This class extracts the signal parameters (energy, time, quality)
	19	// from ALTRO samples. Energy is in ADC counts, time is in time bin units.
	20	// Class uses FastFitting algorithm to fit sample and extract time and Amplitude
	21	// and evaluate sample quality = (chi^2/NDF)/some parameterization providing
	22	// efficiency close to 100%
	23	//
	24	// Typical use case:
	25	// AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
	26	// fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
	27	// fitter->SetCalibData(fgCalibData) ;
	28	// fitter->Eval(sig,sigStart,sigLength);
	29	// Double_t amplitude = fitter.GetEnergy();
	30	// Double_t time = fitter.GetTime();
	31	// Bool_t isLowGain = fitter.GetCaloFlag()==0;
	32
	33	// Author: Dmitri Peressounko (after A.Pavlinov - see RAW/AliCaloFastAltroFitv0.cxx)
	34
	35	// --- ROOT system ---
	36	#include "TArrayI.h"
	37	#include "TList.h"
	38	#include "TMath.h"
	39	#include "TH1I.h"
	40	#include "TF1.h"
	41	#include "TCanvas.h"
	42	#include "TFile.h"
	43	#include "TROOT.h"
	44
	45	// --- AliRoot header files ---
	46	#include "AliLog.h"
	47	#include "AliPHOSCalibData.h"
	48	#include "AliPHOSRawFitterv3.h"
	49	#include "AliPHOSPulseGenerator.h"
	50
	51	ClassImp(AliPHOSRawFitterv3)
	52
	53	//-----------------------------------------------------------------------------
	54	AliPHOSRawFitterv3::AliPHOSRawFitterv3():
	55	AliPHOSRawFitterv0()
	56	{
	57	//Default constructor.
	58	}
	59
	60	//-----------------------------------------------------------------------------
	61	AliPHOSRawFitterv3::~AliPHOSRawFitterv3()
	62	{
	63	//Destructor.
	64	}
65
66	//-----------------------------------------------------------------------------
67	AliPHOSRawFitterv3::AliPHOSRawFitterv3(const AliPHOSRawFitterv3 &phosFitter ):
68	AliPHOSRawFitterv0(phosFitter)
69	{
70	//Copy constructor.
71	}
72
73	//-----------------------------------------------------------------------------
74	AliPHOSRawFitterv3& AliPHOSRawFitterv3::operator = (const AliPHOSRawFitterv3 & /phosFitter/)
75	{
76	//Assignment operator.
77	return *this;
78	}
79
80	//-----------------------------------------------------------------------------
81	Bool_t AliPHOSRawFitterv3::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
82	{
83	//Extract an energy deposited in the crystal,
84	//crystal' position (module,column,row),
85	//time and gain (high or low).
86	//First collects sample, then evaluates it and if it has
87	//reasonable shape, fits it with Gamma2 function and extracts
88	//energy and time.
89
90	if (fCaloFlag == 2 \|\| fNBunches > 1) {
91	fQuality = 1000;
92	return kTRUE;
93	}
94	if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
95	fQuality=2000;
96	fEnergy=0 ;
97	return kTRUE;
98	}
99
100	Float_t pedMean = 0;
101	Float_t pedRMS = 0;
102	Int_t nPed = 0;
103	const Float_t kBaseLine = 1.0;
104	const Int_t kPreSamples = 10;
105
106
107	//We tryed individual taus for each channel, but
108	//this approach seems to be unstable. Much better results are obtaned with
109	//fixed decay time for all channels.
110	const Double_t tau=21. ;
111
112	TArrayD *samples = new TArrayD(sigLength); // array of sample amplitudes
113	TArrayD *times = new TArrayD(sigLength); // array of sample time stamps
114	for (Int_t i=0; i<sigLength; i++) {
115	if (i<kPreSamples) {
116	nPed++;
117	pedMean += signal[i];
118	pedRMS += signal[i]*signal[i] ;
119	}
120	samples->AddAt(signal[i],sigLength-i-1);
121	times ->AddAt(i/tau ,i);
122	}
123
124	fEnergy = -111;
125	fQuality= 999. ;
126	const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
127	Double_t pedestal = 0;
128
129	if (fPedSubtract) {
130	if (nPed > 0) {
131	fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
132	if(fPedestalRMS > 0.)
133	fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
134	fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction
135	}
136	else
137	return kFALSE;
138	}
139	else {
140	//take pedestals from DB
141	pedestal = (Double_t) fAmpOffset ;
142	if (fCalibData) {
143	Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
144	Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
145	pedestal += truePed - altroSettings ;
146	}
147	else{
148	AliWarning(Form("Can not read data from OCDB")) ;
149	}
150	fEnergy-=pedestal ;
151	}
152
153	if (fEnergy < kBaseLine) fEnergy = 0;
154	//Evaluate time
155	Int_t iStart = 0;
156	while(iStart<sigLength && samples->At(iStart)-pedestal <kBaseLine) iStart++ ;
157	fTime = sigStart+iStart;
158
159	//calculate time and energy
160	Int_t maxBin=0 ;
161	Int_t maxAmp=0 ;
162	Int_t nMax = 0 ; //number of points in plato
163	Double_t aMean =0. ;
164	Double_t aRMS =0. ;
165	Double_t wts =0 ;
166	Int_t tStart=0 ;
167
168	for (Int_t i=0; i<sigLength; i++){
169	if(signal[i] > pedestal){
170	Double_t de = signal[i] - pedestal ;
171	if(de > 1.) {
172	aMean += de*i ;
173	aRMS += deii ;
174	wts += de;
175	}
176	if(de > 2 && tStart==0)
177	tStart = i ;
178	if(signal[i] > maxAmp){
179	maxAmp = signal[i];
180	nMax=0;
181	maxBin = i ;
182	}
183	if(signal[i] == maxAmp)
184	nMax++;
185	}
186	}
187
188	if (maxBin==sigLength-1){//bad "rising" sample
189	fEnergy = 0. ;
190	fTime = -999. ;
191	fQuality= 999. ;
192	return kTRUE ;
193	}
194
195	fEnergy=Double_t(maxAmp)-pedestal ;
196	fOverflow =0 ; //look for plato on the top of sample
197	if (fEnergy>500 && //this is not fluctuation of soft sample
198	nMax>2){ //and there is a plato
199	fOverflow = kTRUE ;
200	}
201
202	if (wts > 0) {
203	aMean /= wts;
204	aRMS = aRMS/wts - aMean*aMean;
205	}
206
207	//do not take too small energies
208	if (fEnergy < kBaseLine)
209	fEnergy = 0;
210
211	//do not test quality of too soft samples
212	if (fEnergy < maxEtoFit){
213	fTime = tStart;
214	if (aRMS < 2.) //sigle peak
215	fQuality = 999. ;
216	else
217	fQuality = 0. ;
218	return kTRUE ;
219	}
220
221	// if sample has reasonable mean and RMS, try to fit it with gamma2
222	// First estimate t0
223	Double_t a=0,b=0,c=0 ;
224	for(Int_t i=0; i<sigLength; i++){
225	if(samples->At(i)<pedestal)
226	continue ;
227	Double_t t= times->At(i) ;
228	Double_t f02 = TMath::Exp(-2.*t);
229	Double_t f12 = t*f02;
230	Double_t f22 = t*f12;
231	// Derivatives
232	Double_t f02d = -2.*f02;
233	Double_t f12d = f02 - 2.*f12;
234	Double_t f22d = 2.*(f12 - f22);
235	a += f02d * (samples->At(i)-pedestal) ;
236	b -= f12d * (samples->At(i)-pedestal) ;
237	c += f22d * (samples->At(i)-pedestal) ;
238	}
239
240	//Find 2 roots
241	Double_t det = bb - ac;
242	if(det>=1.e-6 && det<0.0) {
243	det = 0.0; //rounding error
244	}
245	if(det<0.){ //Problem
246	fQuality = 1500. ;
247	if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
248	printf(" det=%e \n",det) ;
249	goto plot ;
250	}
251	return kTRUE ;
252	}
253
254	if(det>0.0 && a!=0.) {
255	det = TMath::Sqrt(det);
256	Double_t t1 = (-b + det) / a;
257	// Double_t t2 = (-b - det) / a; //second root is wrong one
258	Double_t amp1=0., den1=0. ;
259	for(Int_t i=0; i<sigLength; i++){
260	Double_t dt1 = times->At(i) - t1;
261	Double_t f01 = dt1dt1TMath::Exp(-2.*dt1);
262	amp1 += f01*(samples->At(i)-pedestal);
263	den1 += f01*f01;
264	}
265	if(den1>0.0) amp1 /= den1;
266	Double_t chi1=0.; // chi2=0. ;
267	for(Int_t i=0; i<sigLength; i++){
268	Double_t dt1 = times->At(i)- t1;
269	Double_t dy1 = samples->At(i)-pedestal- amp1dt1dt1TMath::Exp(-2.dt1) ;
270	chi1 += dy1*dy1;
271	}
272	fEnergy=amp1*TMath::Exp(-2.) ; ;
273	fTime=t1*tau ;
274	fQuality=chi1/sigLength ;
275	}
276	else {
277	Double_t t1 ;
278	if(a!=0.)
279	t1=-b/a ;
280	else
281	t1=-c/b ;
282	Double_t amp=0.,den=0.; ;
283	for(Int_t i=0; i<sigLength; i++){
284	Double_t dt = times->At(i) - t1;
285	Double_t f = dtdtTMath::Exp(-2.*dt);
286	amp += f*samples->At(i);
287	den += f*f;
288	}
289	if(den>0.0) amp /= den;
290	// chi2 calculation
291	fQuality=0.;
292	for(Int_t i=0; i<sigLength; i++){
293	Double_t t = times->At(i)- t1;
294	Double_t dy = samples->At(i)- ampttTMath::Exp(-2.t) ;
295	fQuality += dy*dy;
296	}
297	fTime=t1*tau ;
298	fEnergy = amp*TMath::Exp(-2.);
299	fQuality/= sigLength ;
300	}
301
302	//Impose cut on quality
303	fQuality/=2.+0.004fEnergyfEnergy ;
304
305	//Draw corrupted samples
306	if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
307	plot:
308	if(fQuality >1. && !fOverflow ){ //Draw only bad samples
309	printf("Sample par: amp=%f, t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
310	TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
311	if(!h) h = new TH1I("h","Samples",50,0.,50.) ;
312	h->Reset() ;
313	for (Int_t i=0; i<sigLength; i++) {
314	h->SetBinContent(i,samples->At(i)) ;
315	}
316	TF1 * fffit = new TF1("fffit","[0]+[1]((x-[2])/[3])^2exp(2.-2.*(x-[2])/[3])",0.,200.) ;
317	fffit->SetParameters(pedestal,fEnergy,fTime,tau) ;
318	fffit->SetLineColor(2) ;
319	TCanvas * c = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
320	if(!c){
321	c = new TCanvas("cSamples","cSamples",10,10,400,400) ;
322	c->SetFillColor(0) ;
323	c->SetFillStyle(0) ;
324	c->Range(0,0,1,1);
325	c->SetBorderSize(0);
326	}
327	c->cd() ;
328
329	TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
330	spectrum_1->Draw();
331	spectrum_1->cd();
332	spectrum_1->Range(0,0,1,1);
333	spectrum_1->SetFillColor(0);
334	spectrum_1->SetFillStyle(4000);
335	spectrum_1->SetBorderSize(1);
336	spectrum_1->SetBottomMargin(0.012);
337	spectrum_1->SetTopMargin(0.03);
338	spectrum_1->SetLeftMargin(0.10);
339	spectrum_1->SetRightMargin(0.05);
340
341	char title[155] ;
342	sprintf(title,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
343	h->SetTitle(title) ;
344	h->Draw() ;
345	fffit->Draw("same") ;
346
347	sprintf(title,"mod%d_x%d_z%d_HG",fModule,fCellX,fCellZ) ;
348	TFile fout("samples_bad.root","update") ;
349	h->Write(title);
350	fout.Close() ;
351
352	c->cd() ;
353	TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
354	spectrum_2->SetFillColor(0) ;
355	spectrum_2->SetFillStyle(0) ;
356	spectrum_2->SetGridy() ;
357	spectrum_2->Draw();
358	spectrum_2->Range(0,0,1,1);
359	spectrum_2->SetFillColor(0);
360	spectrum_2->SetBorderSize(1);
361	spectrum_2->SetTopMargin(0.01);
362	spectrum_2->SetBottomMargin(0.25);
363	spectrum_2->SetLeftMargin(0.10);
364	spectrum_2->SetRightMargin(0.05);
365	spectrum_2->cd() ;
366
367	TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
368	if(!hd) hd = new TH1I("hd","Samples",50,0.,50.) ;
369	hd->Reset() ;
370	for (Int_t i=0; i<sigLength; i++) {
371	hd->SetBinContent(i,TMath::Max(-1023.,TMath::Min(1023.,samples->At(i)-fffit->Eval(i)))) ;
372	}
373	hd->Draw();
374
375	c->Update() ;
376	printf("Press <enter> to continue\n") ;
377	getchar();
378
379
380	delete fffit ;
381	delete spectrum_1 ;
382	delete spectrum_2 ;
383	}
384	}
385
386	delete samples ;
387	delete times ;
388	return kTRUE;
389	}