Updated DA from Brigitte. Steerable clock range for pedestal evaluation - default...
[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawFitterv3.cxx
CommitLineData
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
51ClassImp(AliPHOSRawFitterv3)
52
53//-----------------------------------------------------------------------------
54AliPHOSRawFitterv3::AliPHOSRawFitterv3():
55 AliPHOSRawFitterv0()
56{
57 //Default constructor.
58}
59
60//-----------------------------------------------------------------------------
61AliPHOSRawFitterv3::~AliPHOSRawFitterv3()
62{
63 //Destructor.
64}
65
66//-----------------------------------------------------------------------------
67AliPHOSRawFitterv3::AliPHOSRawFitterv3(const AliPHOSRawFitterv3 &phosFitter ):
68 AliPHOSRawFitterv0(phosFitter)
69{
70 //Copy constructor.
71}
72
73//-----------------------------------------------------------------------------
74AliPHOSRawFitterv3& AliPHOSRawFitterv3::operator = (const AliPHOSRawFitterv3 & /*phosFitter*/)
75{
76 //Assignment operator.
77 return *this;
78}
79
80//-----------------------------------------------------------------------------
81Bool_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 += de*i*i ;
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 = b*b - a*c;
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 = dt1*dt1*TMath::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- amp1*dt1*dt1*TMath::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 = dt*dt*TMath::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)- amp*t*t*TMath::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.004*fEnergy*fEnergy ;
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])^2*exp(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}