Copy TPCnclsF
[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawFitterv2.cxx
CommitLineData
379c5c09 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
6a265faf 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%
379c5c09 23//
24// Typical use case:
25// AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
379c5c09 26// fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
27// fitter->SetCalibData(fgCalibData) ;
1dfadc32 28// fitter->Eval(sig,sigStart,sigLength);
379c5c09 29// Double_t amplitude = fitter.GetEnergy();
30// Double_t time = fitter.GetTime();
31// Bool_t isLowGain = fitter.GetCaloFlag()==0;
32
6a265faf 33// Author: Dmitri Peressounko (after A.Pavlinov - see RAW/AliCaloFastAltroFitv0.cxx)
379c5c09 34
35// --- ROOT system ---
6a265faf 36#include "TArrayI.h"
379c5c09 37#include "TList.h"
38#include "TMath.h"
6a265faf 39#include "TH1I.h"
379c5c09 40#include "TF1.h"
6a265faf 41#include "TCanvas.h"
42#include "TFile.h"
379c5c09 43#include "TROOT.h"
44
45// --- AliRoot header files ---
46#include "AliLog.h"
47#include "AliPHOSCalibData.h"
48#include "AliPHOSRawFitterv2.h"
49#include "AliPHOSPulseGenerator.h"
50
51ClassImp(AliPHOSRawFitterv2)
52
53//-----------------------------------------------------------------------------
54AliPHOSRawFitterv2::AliPHOSRawFitterv2():
55 AliPHOSRawFitterv0(),
6a265faf 56 fAlpha(0.1),fBeta(0.035),fMax(0)
379c5c09 57{
58 //Default constructor.
379c5c09 59}
60
61//-----------------------------------------------------------------------------
62AliPHOSRawFitterv2::~AliPHOSRawFitterv2()
63{
64 //Destructor.
65}
66
67//-----------------------------------------------------------------------------
68AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ):
6a265faf 69 AliPHOSRawFitterv0(phosFitter),
70 fAlpha(0.1),fBeta(0.035),fMax(0)
379c5c09 71{
72 //Copy constructor.
379c5c09 73}
74
75//-----------------------------------------------------------------------------
6a265faf 76AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 & /*phosFitter*/)
379c5c09 77{
78 //Assignment operator.
379c5c09 79 return *this;
80}
81
82//-----------------------------------------------------------------------------
92236b27 83Bool_t AliPHOSRawFitterv2::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
379c5c09 84{
85 //Extract an energy deposited in the crystal,
86 //crystal' position (module,column,row),
87 //time and gain (high or low).
88 //First collects sample, then evaluates it and if it has
89 //reasonable shape, fits it with Gamma2 function and extracts
90 //energy and time.
91
6a265faf 92 const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
379c5c09 93 const Float_t kBaseLine = 1.0;
94 const Int_t kPreSamples = 10;
95
6a265faf 96 fOverflow = kFALSE ;
97 fEnergy=0 ;
98 if (fCaloFlag == 2 || fNBunches > 1) {
99 fQuality = 150;
100 return kTRUE;
101 }
102 if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
103 fQuality=200;
104 fEnergy=0 ;
105 return kTRUE;
106 }
379c5c09 107
6a265faf 108 //Evaluate pedestals
109 Float_t pedMean = 0;
110 Float_t pedRMS = 0;
111 Int_t nPed = 0;
112 for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) {
113 nPed++;
114 pedMean += signal[i];
115 pedRMS += signal[i]*signal[i] ;
116 }
379c5c09 117
6a265faf 118 fEnergy = -111;
119 fQuality= 999. ;
379c5c09 120 Double_t pedestal = 0;
379c5c09 121
122 if (fPedSubtract) {
123 if (nPed > 0) {
124 fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
125 if(fPedestalRMS > 0.)
126 fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
6a265faf 127 pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction
379c5c09 128 }
129 else
130 return kFALSE;
131 }
132 else {
133 //take pedestals from DB
134 pedestal = (Double_t) fAmpOffset ;
379c5c09 135 }
136
379c5c09 137
6a265faf 138 //calculate rough quality of the sample and check for overflow
139 Int_t maxBin=0 ;
140 Int_t maxAmp=0 ;
141 Int_t minAmp= signal[0] ;
142 Int_t nMax = 0 ; //number of points in plato
143 Double_t aMean =0. ;
144 Double_t aRMS =0. ;
145 Double_t wts =0 ;
146 Bool_t falling = kTRUE ; //Bad monotoneusly falling sample
147 Bool_t rising = kTRUE ; //Bad monotoneusly riging sample
148 for (Int_t i=0; i<sigLength; i++){
149 if(signal[i] > pedestal){
150 Double_t de = signal[i] - pedestal ;
151 if(de > 1.) {
152 aMean += de*i ;
153 aRMS += de*i*i ;
154 wts += de;
155 }
156 if(signal[i] > maxAmp){
157 maxAmp = signal[i];
158 nMax=0;
159 maxBin = i ;
160 }
161 if(signal[i] == maxAmp){
162 nMax++;
163 }
164 if(signal[i] < minAmp)
165 minAmp=signal[i] ;
166 if(falling && i>0 && signal[i]<signal[i-1])
167 falling=kFALSE ;
168 if(rising && i>0 && signal[i]>signal[i-1])
169 rising=kFALSE ;
170 }
171 }
172
173 if(rising || falling){//bad "rising" or falling sample
174 fEnergy = 0. ;
175 fTime = 0. ; //-999. ;
176 fQuality= 250. ;
177 return kTRUE ;
379c5c09 178 }
6a265faf 179 if(maxAmp-minAmp<3 && maxAmp>7 && sigLength>20){ //bad flat sample
180 fEnergy = 0. ;
181 fTime = 0; //-999. ;
182 fQuality= 260. ;
183 return kTRUE ;
379c5c09 184 }
185
6a265faf 186 fEnergy=Double_t(maxAmp)-pedestal ;
187 if (fEnergy < kBaseLine) fEnergy = 0;
188 fTime = sigStart-sigLength-3;
189
190 //do not test quality of too soft samples
191 if (wts > 0) {
192 aMean /= wts;
193 aRMS = aRMS/wts - aMean*aMean;
194 }
195 if (fEnergy <= maxEtoFit){
196 if (aRMS < 2.) //sigle peak
197 fQuality = 299. ;
198 else
199 fQuality = 0. ;
200 //Evaluate time of signal arriving
201 return kTRUE ;
202 }
379c5c09 203
6a265faf 204 //look for plato on the top of sample
205 if (fEnergy>500 && //this is not fluctuation of soft sample
206 nMax>2){ //and there is a plato
207 fOverflow = kTRUE ;
208 }
379c5c09 209
6a265faf 210
211 //do not fit High Gain samples with overflow
212 if(fCaloFlag==1 && fOverflow){
213 fQuality = 99. ;
214 return kTRUE;
215
216 }
217
218 //----Now fit sample with reasonable shape------
219 TArrayD samples(sigLength); // array of sample amplitudes
220 TArrayD times(sigLength); // array of sample time stamps
221 for (Int_t i=0; i<sigLength; i++) {
222 samples.AddAt(signal[i]-pedestal,sigLength-i-1);
223 times.AddAt(double(i),i);
224 }
225
226 if(fMax==0)
227 FindMax() ;
228 if(!FindAmpT(samples,times)){
229 if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv2")>3){
230 goto plot ;
231 }
232 else{
233 return kFALSE ;
379c5c09 234 }
235 }
6a265faf 236 fEnergy*=fMax ;
237 fTime += sigStart-sigLength-3;
238
239
240 //Impose cut on quality
241// fQuality/=4. ;
242 fQuality/=1.+0.005*fEnergy ;
243
244 //Draw corrupted samples
245 if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv2")>3){
246 if(fEnergy > 50. ){
247 plot:
248 printf("Sample par: amp=%f, t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
249 TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
250 if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
251 h->Reset() ;
252 for (Int_t i=0; i<sigLength; i++) {
253 h->SetBinContent(i+1,float(samples.At(i))) ;
254 }
255// TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
256 TF1 * fffit = new TF1("fffit","[0]*((x-[1])*(x-[1])*exp(-[2]*(x-[1]))+(x-[1])*exp(-[3]*(x-[1])))",0.,60.) ;
257 fffit->SetParameters(fEnergy/fMax,fTime-(sigStart-sigLength-3),fAlpha,fBeta) ;
258 fffit->SetLineColor(2) ;
259 TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
260 if(!can){
261 can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
262 can->SetFillColor(0) ;
263 can->SetFillStyle(0) ;
264 can->Range(0,0,1,1);
265 can->SetBorderSize(0);
266 }
267 can->cd() ;
379c5c09 268
6a265faf 269 TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
270 spectrum_1->Draw();
271 spectrum_1->cd();
272 spectrum_1->Range(0,0,1,1);
273 spectrum_1->SetFillColor(0);
274 spectrum_1->SetFillStyle(4000);
275 spectrum_1->SetBorderSize(1);
276 spectrum_1->SetBottomMargin(0.012);
277 spectrum_1->SetTopMargin(0.03);
278 spectrum_1->SetLeftMargin(0.10);
279 spectrum_1->SetRightMargin(0.05);
280
281 char title[155] ;
282 sprintf(title,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
283 h->SetTitle(title) ;
284// h->Fit(fffit,"","",0.,51.) ;
285 h->Draw() ;
286 fffit->Draw("same") ;
287/*
288 sprintf(title,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ;
289 TFile fout("samples_bad.root","update") ;
290 h->Write(title);
291 fout.Close() ;
292*/
293 can->cd() ;
294 TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
295 spectrum_2->SetFillColor(0) ;
296 spectrum_2->SetFillStyle(0) ;
297 spectrum_2->SetGridy() ;
298 spectrum_2->Draw();
299 spectrum_2->Range(0,0,1,1);
300 spectrum_2->SetFillColor(0);
301 spectrum_2->SetBorderSize(1);
302 spectrum_2->SetTopMargin(0.01);
303 spectrum_2->SetBottomMargin(0.25);
304 spectrum_2->SetLeftMargin(0.10);
305 spectrum_2->SetRightMargin(0.05);
306 spectrum_2->cd() ;
307
308 TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
309 if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
310 hd->Reset() ;
311 for (Int_t i=0; i<sigLength; i++) {
312 hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ;
313 }
314 hd->Draw();
315
316 can->Update() ;
317 printf("Press <enter> to continue\n") ;
318 getchar();
319
320
321 delete fffit ;
322 delete spectrum_1 ;
323 delete spectrum_2 ;
324 }
379c5c09 325 }
326
379c5c09 327 return kTRUE;
328}
6a265faf 329//------------------------------------------------------------------
330Bool_t AliPHOSRawFitterv2::FindAmpT(TArrayD samples, TArrayD times){
331// makes fit
332
333 const Int_t nMaxIter=50 ; //Maximal number of iterations
334 const Double_t epsdt = 1.e-3 ; //expected precision of t0 calculation
335
336 Double_t dTime=times.At(0)-0.5 ; //Most probable Initial approximation
337//printf(" start fit... \n") ;
338
339 Int_t nPoints = samples.GetSize() ;
340 Double_t dea=TMath::Exp(-fAlpha) ;
341 Double_t deb=TMath::Exp(-fBeta) ;
342 Double_t dt=1.,timeOld=dTime,dfOld=0. ;
343 for(Int_t iter=0; iter<nMaxIter; iter++){
344 Double_t yy=0.;
345 Double_t yf=0. ;
346 Double_t ydf=0. ;
347 Double_t yddf=0. ;
348 Double_t ff=0. ;
349 Double_t fdf=0. ;
350 Double_t dfdf=0. ;
351 Double_t fddf=0. ;
352 Int_t nfit=0 ;
353 Double_t aexp=TMath::Exp(-fAlpha*(times.At(0)-1.-dTime)) ;
354 Double_t bexp=TMath::Exp(-fBeta*(times.At(0)-1.-dTime)) ;
355 for(Int_t i=0; i<nPoints; i++){
356 Double_t t= times.At(i)-dTime ;
357 aexp*=dea ;
358 bexp*=deb ;
359 if(t<0.) continue ;
360 Double_t y=samples.At(i) ;
361 if(y<=fAmpThreshold)
362 continue ;
363 nfit++ ;
364 Double_t at=fAlpha*t ;
365 Double_t bt = fBeta*t ;
366 Double_t phi=t*(t*aexp+bexp) ;
367 Double_t dphi=t*aexp*(2.-at)+(1.-bt)*bexp ;
368 Double_t ddphi=aexp*(2.-4.*at+at*at)+bexp*fBeta*(bt-2.) ;
369 yy+=y*y ;
370 yf+=y*phi ;
371 ydf+=y*dphi ;
372 yddf+=y*ddphi ;
373 ff+=phi*phi ;
374 fdf+=phi*dphi ;
375 dfdf+=dphi*dphi ;
376 fddf+=phi*ddphi ;
377 }
378
379 if(ff==0.||nfit==0. ){
380 fQuality=199 ;
381 return kFALSE ;
382 }
383 Double_t f=ydf*ff-yf*fdf ; //d(chi2)/dt
384 Double_t df=yf*(dfdf+fddf)-yddf*ff-ydf*fdf;
385 if(df<=0.){ //we are too far from the root. In the wicinity of root df>0
386 if(iter!=0 && dfOld>0.){//If at previous step df was OK, just reduce step size
387 dt*=0.5 ;
388 dTime=timeOld+dt ;
389 continue ;
390 }
391 if(f<0){ //f<0 => dTime is too small and we still do not know root region
392 dTime+=2. ;
393 continue ;
394 }
395 else{ //dTime is too large, we are beyond the root region
396 dTime-=2. ;
397 continue ;
398 }
399 }
400 dt=-f/df ;
401 if(TMath::Abs(dt)<epsdt){
402 fQuality=(yy-yf*yf/ff)/nfit ;
403 fEnergy=yf/ff ; //ff!=0 already tested
404 fTime=dTime ;
405 return kTRUE ;
406 }
407 //In some cases time steps are huge (derivative ~0)
408 if(dt>10.) dt=10. ; //restrict step size
409 if(dt<-10.) dt=-5.3 ; //this restriction should be asimmetric to avoid jumping from one point to another
410 timeOld=dTime ; //remember current position for the case
411 dfOld=df ; //of reduction of dt step size
412 dTime+=dt ;
413
414 if(dTime>100. || dTime<-30.){ //this is corrupted sample, do not spend time improving accuracy.
415 fQuality=(yy-yf*yf/ff)/nfit ;
416 fEnergy=yf/ff ; //ff!=0 already tested
417 fTime=dTime ;
418 return kFALSE ;
419 }
420
421 }
422 //failed to find a root, too many iterations
423 fQuality=99.;
424 fEnergy=0 ;
425 return kFALSE ;
426}
427//_________________________________________
428void AliPHOSRawFitterv2::FindMax(){
429 //Finds maxumum of currecnt parameterization
430 Double_t t=2./fAlpha ;
431 fMax = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;
432 Double_t dt=15 ;
433 while(dt>0.01){
434 Double_t dfdt=(2.*t-fAlpha*t*t)*TMath::Exp(-fAlpha*t)+(1.-fBeta*t)*TMath::Exp(-fBeta*t) ;
435 if(dfdt>0.)
436 t+=dt ;
437 else
438 t-=dt ;
439 Double_t maxNew = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;
440 if(maxNew>fMax)
441 fMax=maxNew ;
442 else{
443 dt/=2 ;
444 if(dfdt<0.)
445 t+=dt ;
446 else
447 t-=dt ;
448 }
449 }
450}
451