[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(sigLength); // array of sample amplitudes
  TArrayD times(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
  fTime = sigStart; 
  
  //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 ;
  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(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){
    if (aRMS < 2.) //sigle peak
      fQuality = 999. ;
    else
      fQuality =   0. ;
    //Evaluate time of signal arriving
    Int_t i=0;
    while(signal[sigLength-i-1]<pedestal+kBaseLine  && i<sigLength)
     i++ ;
    fTime+=i ;
    return kTRUE ;
  }
      
  // if sample has reasonable mean and RMS, try to fit it with gamma2
  //This method can not analyse overflow samples
  if(fOverflow){
    fQuality = 99. ;
    return kTRUE ;
  }
  // 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++){
      if(samples.At(i)<pedestal)
        continue ;
      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++){
      if(samples.At(i)<=pedestal)
        continue ;
      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++){
      if(samples.At(i)<=pedestal)
        continue ;
      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++){
      if(samples.At(i)<=pedestal)
        continue ;
      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 ; //If we have overflow the number of actually fitted points is smaller, but chi2 in this case is not important.
  }

  //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-sigStart,tau) ;
      fffit->SetLineColor(2) ;
      TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
      if(!can){
        can = new TCanvas("cSamples","cSamples",10,10,400,400) ;
        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);
      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() ;

      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",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();
 
      can->Update() ;
      printf("Press <enter> to continue\n") ;
      getchar();


      delete fffit ;
      delete spectrum_1 ;
      delete spectrum_2 ;
    }
  }
  
  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
ad308f96	112	TArrayD samples(sigLength); // array of sample amplitudes
ad308f96	113	TArrayD times(sigLength); // array of sample time stamps
c4cc99c3	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	}
ad308f96	120	samples.AddAt(signal[i],sigLength-i-1);
ad308f96	121	times.AddAt(i/tau ,i);
c4cc99c3	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
facd0d71	155	fTime = sigStart;
c4cc99c3	156
	157	//calculate time and energy
	158	Int_t maxBin=0 ;
	159	Int_t maxAmp=0 ;
	160	Int_t nMax = 0 ; //number of points in plato
	161	Double_t aMean =0. ;
	162	Double_t aRMS =0. ;
	163	Double_t wts =0 ;
c4cc99c3	164	for (Int_t i=0; i<sigLength; i++){
	165	if(signal[i] > pedestal){
	166	Double_t de = signal[i] - pedestal ;
	167	if(de > 1.) {
	168	aMean += de*i ;
	169	aRMS += deii ;
	170	wts += de;
	171	}
c4cc99c3	172	if(signal[i] > maxAmp){
	173	maxAmp = signal[i];
	174	nMax=0;
	175	maxBin = i ;
	176	}
	177	if(signal[i] == maxAmp)
	178	nMax++;
	179	}
	180	}
	181
	182	if (maxBin==sigLength-1){//bad "rising" sample
	183	fEnergy = 0. ;
	184	fTime = -999. ;
	185	fQuality= 999. ;
	186	return kTRUE ;
	187	}
	188
	189	fEnergy=Double_t(maxAmp)-pedestal ;
	190	fOverflow =0 ; //look for plato on the top of sample
	191	if (fEnergy>500 && //this is not fluctuation of soft sample
	192	nMax>2){ //and there is a plato
	193	fOverflow = kTRUE ;
	194	}
	195
	196	if (wts > 0) {
	197	aMean /= wts;
	198	aRMS = aRMS/wts - aMean*aMean;
	199	}
	200
	201	//do not take too small energies
	202	if (fEnergy < kBaseLine)
	203	fEnergy = 0;
	204
	205	//do not test quality of too soft samples
	206	if (fEnergy < maxEtoFit){
c4cc99c3	207	if (aRMS < 2.) //sigle peak
	208	fQuality = 999. ;
	209	else
	210	fQuality = 0. ;
facd0d71	211	//Evaluate time of signal arriving
	212	Int_t i=0;
	213	while(signal[sigLength-i-1]<pedestal+kBaseLine && i<sigLength)
	214	i++ ;
	215	fTime+=i ;
c4cc99c3	216	return kTRUE ;
	217	}
	218
	219	// if sample has reasonable mean and RMS, try to fit it with gamma2
facd0d71	220	//This method can not analyse overflow samples
	221	if(fOverflow){
	222	fQuality = 99. ;
	223	return kTRUE ;
	224	}
c4cc99c3	225	// First estimate t0
	226	Double_t a=0,b=0,c=0 ;
	227	for(Int_t i=0; i<sigLength; i++){
ad308f96	228	if(samples.At(i)<=pedestal)
c4cc99c3	229	continue ;
ad308f96	230	Double_t t= times.At(i) ;
c4cc99c3	231	Double_t f02 = TMath::Exp(-2.*t);
	232	Double_t f12 = t*f02;
	233	Double_t f22 = t*f12;
	234	// Derivatives
	235	Double_t f02d = -2.*f02;
	236	Double_t f12d = f02 - 2.*f12;
	237	Double_t f22d = 2.*(f12 - f22);
ad308f96	238	a += f02d * (samples.At(i)-pedestal) ;
	239	b -= f12d * (samples.At(i)-pedestal) ;
	240	c += f22d * (samples.At(i)-pedestal) ;
c4cc99c3	241	}
	242
	243	//Find 2 roots
	244	Double_t det = bb - ac;
	245	if(det>=1.e-6 && det<0.0) {
	246	det = 0.0; //rounding error
	247	}
	248	if(det<0.){ //Problem
	249	fQuality = 1500. ;
	250	if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
	251	printf(" det=%e \n",det) ;
	252	goto plot ;
	253	}
	254	return kTRUE ;
	255	}
	256
	257	if(det>0.0 && a!=0.) {
	258	det = TMath::Sqrt(det);
	259	Double_t t1 = (-b + det) / a;
	260	// Double_t t2 = (-b - det) / a; //second root is wrong one
	261	Double_t amp1=0., den1=0. ;
	262	for(Int_t i=0; i<sigLength; i++){
ad308f96	263	if(samples.At(i)<pedestal)
facd0d71	264	continue ;
ad308f96	265	Double_t dt1 = times.At(i) - t1;
facd0d71	266	Double_t f01 = dt1dt1TMath::Exp(-2.*dt1);
ad308f96	267	amp1 += f01*(samples.At(i)-pedestal);
facd0d71	268	den1 += f01*f01;
	269	}
	270	if(den1>0.0) amp1 /= den1;
	271	Double_t chi1=0.; // chi2=0. ;
	272	for(Int_t i=0; i<sigLength; i++){
ad308f96	273	if(samples.At(i)<=pedestal)
facd0d71	274	continue ;
ad308f96	275	Double_t dt1 = times.At(i)- t1;
ad308f96	276	Double_t dy1 = samples.At(i)-pedestal- amp1dt1dt1TMath::Exp(-2.dt1) ;
facd0d71	277	chi1 += dy1*dy1;
	278	}
	279	fEnergy=amp1*TMath::Exp(-2.) ; ;
	280	fTime+=t1*tau ;
	281	fQuality=chi1/sigLength ;
c4cc99c3	282	}
	283	else {
	284	Double_t t1 ;
	285	if(a!=0.)
	286	t1=-b/a ;
	287	else
	288	t1=-c/b ;
	289	Double_t amp=0.,den=0.; ;
	290	for(Int_t i=0; i<sigLength; i++){
ad308f96	291	if(samples.At(i)<=pedestal)
facd0d71	292	continue ;
ad308f96	293	Double_t dt = times.At(i) - t1;
facd0d71	294	Double_t f = dtdtTMath::Exp(-2.*dt);
ad308f96	295	amp += f*samples.At(i);
facd0d71	296	den += f*f;
	297	}
	298	if(den>0.0) amp /= den;
	299	// chi2 calculation
	300	fQuality=0.;
	301	for(Int_t i=0; i<sigLength; i++){
ad308f96	302	if(samples.At(i)<=pedestal)
facd0d71	303	continue ;
ad308f96	304	Double_t t = times.At(i)- t1;
ad308f96	305	Double_t dy = samples.At(i)- ampttTMath::Exp(-2.t) ;
facd0d71	306	fQuality += dy*dy;
	307	}
	308	fTime+=t1*tau ;
	309	fEnergy = amp*TMath::Exp(-2.);
	310	fQuality/= sigLength ; //If we have overflow the number of actually fitted points is smaller, but chi2 in this case is not important.
c4cc99c3	311	}
	312
	313	//Impose cut on quality
	314	fQuality/=2.+0.004fEnergyfEnergy ;
	315
	316	//Draw corrupted samples
	317	if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
	318	plot:
	319	if(fQuality >1. && !fOverflow ){ //Draw only bad samples
	320	printf("Sample par: amp=%f, t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
	321	TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
	322	if(!h) h = new TH1I("h","Samples",50,0.,50.) ;
	323	h->Reset() ;
	324	for (Int_t i=0; i<sigLength; i++) {
ad308f96	325	h->SetBinContent(i,samples.At(i)) ;
c4cc99c3	326	}
c4cc99c3	327	TF1 * fffit = new TF1("fffit","[0]+[1]((x-[2])/[3])^2exp(2.-2.*(x-[2])/[3])",0.,200.) ;
facd0d71	328	fffit->SetParameters(pedestal,fEnergy,fTime-sigStart,tau) ;
c4cc99c3	329	fffit->SetLineColor(2) ;
facd0d71	330	TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
	331	if(!can){
	332	can = new TCanvas("cSamples","cSamples",10,10,400,400) ;
	333	can->SetFillColor(0) ;
	334	can->SetFillStyle(0) ;
	335	can->Range(0,0,1,1);
	336	can->SetBorderSize(0);
c4cc99c3	337	}
facd0d71	338	can->cd() ;
c4cc99c3	339
	340	TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
	341	spectrum_1->Draw();
	342	spectrum_1->cd();
	343	spectrum_1->Range(0,0,1,1);
	344	spectrum_1->SetFillColor(0);
	345	spectrum_1->SetFillStyle(4000);
	346	spectrum_1->SetBorderSize(1);
	347	spectrum_1->SetBottomMargin(0.012);
	348	spectrum_1->SetTopMargin(0.03);
	349	spectrum_1->SetLeftMargin(0.10);
	350	spectrum_1->SetRightMargin(0.05);
	351
	352	char title[155] ;
	353	sprintf(title,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
	354	h->SetTitle(title) ;
	355	h->Draw() ;
	356	fffit->Draw("same") ;
	357
	358	sprintf(title,"mod%d_x%d_z%d_HG",fModule,fCellX,fCellZ) ;
	359	TFile fout("samples_bad.root","update") ;
	360	h->Write(title);
	361	fout.Close() ;
	362
facd0d71	363	can->cd() ;
c4cc99c3	364	TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
	365	spectrum_2->SetFillColor(0) ;
	366	spectrum_2->SetFillStyle(0) ;
	367	spectrum_2->SetGridy() ;
	368	spectrum_2->Draw();
	369	spectrum_2->Range(0,0,1,1);
	370	spectrum_2->SetFillColor(0);
	371	spectrum_2->SetBorderSize(1);
	372	spectrum_2->SetTopMargin(0.01);
	373	spectrum_2->SetBottomMargin(0.25);
	374	spectrum_2->SetLeftMargin(0.10);
	375	spectrum_2->SetRightMargin(0.05);
	376	spectrum_2->cd() ;
	377
	378	TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
	379	if(!hd) hd = new TH1I("hd","Samples",50,0.,50.) ;
	380	hd->Reset() ;
	381	for (Int_t i=0; i<sigLength; i++) {
ad308f96	382	hd->SetBinContent(i,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)-fffit->Eval(i)))) ;
c4cc99c3	383	}
	384	hd->Draw();
	385
facd0d71	386	can->Update() ;
c4cc99c3	387	printf("Press <enter> to continue\n") ;
	388	getchar();
	389
	390
	391	delete fffit ;
	392	delete spectrum_1 ;
	393	delete spectrum_2 ;
	394	}
	395	}
	396
c4cc99c3	397	return kTRUE;
c4cc99c3	398	}