make the low multiplicity parameters the default parameters
[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
ad308f96 112 TArrayD samples(sigLength); // array of sample amplitudes
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);
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 += de*i*i ;
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 = b*b - a*c;
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 = dt1*dt1*TMath::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;
276 Double_t dy1 = samples.At(i)-pedestal- amp1*dt1*dt1*TMath::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 = dt*dt*TMath::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;
305 Double_t dy = samples.At(i)- amp*t*t*TMath::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.004*fEnergy*fEnergy ;
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 }
327 TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(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;
398}