[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=22.18 ;

  TArrayD samples(sigLength); // array of sample amplitudes
  TArrayD times(sigLength); // array of sample time stamps
  for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) {
    nPed++;
    pedMean += signal[i];
    pedRMS  += signal[i]*signal[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) ;
      pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction
      fEnergy -= pedestal; // 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-sigLength-3;
  for (Int_t i=0; i<sigLength; i++) {
    samples.AddAt(signal[i]-pedestal,sigLength-i-1);
    times.AddAt(i/tau ,i);
  }
  
  //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
    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 ;
  Int_t minI=0 ;
  if (fPedSubtract) 
    minI=kPreSamples ;
  for(Int_t i=minI; i<sigLength; i++){
    if(samples.At(i)<=0.)
      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) ;
    b -= f12d * samples.At(i) ;
    c += f22d * samples.At(i) ;
  }
  
  //Find roots
  if(a==0.){
    if(b==0.){ //no roots
      fQuality = 2000. ;
      if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
        printf(" a=%f, b=%f, c=%f \n",a,b,c) ;
        goto plot ;
      }
      return kTRUE ;
    }
    Double_t t1=-c/b ;
    Double_t amp=0.,den=0.; ;
    for(Int_t i=minI; i<sigLength; i++){
      if(samples.At(i)<=0.)
        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=minI; i<sigLength; i++){
      if(samples.At(i)<=0.)
        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.
  }
  else{
    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 ;
    }

    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=minI; i<sigLength; i++){
      if(samples.At(i)<=0.)
        continue ;
      Double_t dt1 = times.At(i) - t1;
      Double_t f01 = dt1*dt1*TMath::Exp(-2.*dt1);
      amp1 += f01*samples.At(i);
      den1 += f01*f01;
    }
    if(den1>0.0) amp1 /= den1;
    Double_t chi1=0.; // chi2=0. ;
    for(Int_t i=minI; i<sigLength; i++){
      if(samples.At(i)<=0.)
        continue ;
      Double_t dt1 = times.At(i)- t1;
      Double_t dy1 = samples.At(i)- amp1*dt1*dt1*TMath::Exp(-2.*dt1) ;
      chi1 += dy1*dy1;
    }
    fEnergy=amp1*TMath::Exp(-2.) ; ; 
    fTime+=t1*tau ;
    fQuality=chi1/sigLength ;
  } 

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

  //Draw corrupted samples
  if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
    plot:
    if(fEnergy > 30. && fQuality >1. && !fOverflow ){ //Draw only bad samples
//    if(!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("hSamples","Samples",65,0.,65.) ;
      h->Reset() ;
      for (Int_t i=0; i<sigLength; i++) {
        h->SetBinContent(i+1,samples.At(i)+pedestal) ;
      }
      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 * 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);
      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->Draw() ;
      fffit->Draw("same") ;

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