]>
Commit | Line | Data |
---|---|---|
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=22.18 ; | |
111 | ||
112 | TArrayD samples(sigLength); // array of sample amplitudes | |
113 | TArrayD times(sigLength); // array of sample time stamps | |
114 | for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) { | |
115 | nPed++; | |
116 | pedMean += signal[i]; | |
117 | pedRMS += signal[i]*signal[i] ; | |
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) ; | |
130 | pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction | |
131 | fEnergy -= pedestal; // pedestal subtraction | |
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{ | |
145 | // AliWarning(Form("Can not read data from OCDB")) ; | |
146 | } | |
147 | fEnergy-=pedestal ; | |
148 | } | |
149 | ||
150 | if (fEnergy < kBaseLine) fEnergy = 0; | |
151 | //Evaluate time | |
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 | } | |
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 ; | |
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 += de*i*i ; | |
171 | wts += de; | |
172 | } | |
173 | if(signal[i] > maxAmp){ | |
174 | maxAmp = signal[i]; | |
175 | nMax=0; | |
176 | maxBin = i ; | |
177 | } | |
178 | if(signal[i] == maxAmp){ | |
179 | nMax++; | |
180 | } | |
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){ | |
209 | if (aRMS < 2.) //sigle peak | |
210 | fQuality = 999. ; | |
211 | else | |
212 | fQuality = 0. ; | |
213 | //Evaluate time of signal arriving | |
214 | return kTRUE ; | |
215 | } | |
216 | ||
217 | // if sample has reasonable mean and RMS, try to fit it with gamma2 | |
218 | //This method can not analyse overflow samples | |
219 | if(fOverflow){ | |
220 | fQuality = 99. ; | |
221 | return kTRUE ; | |
222 | } | |
223 | // First estimate t0 | |
224 | Double_t a=0,b=0,c=0 ; | |
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.) | |
230 | continue ; | |
231 | Double_t t= times.At(i) ; | |
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); | |
239 | a += f02d * samples.At(i) ; | |
240 | b -= f12d * samples.At(i) ; | |
241 | c += f22d * samples.At(i) ; | |
242 | } | |
243 | ||
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 ; | |
253 | } | |
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 = dt*dt*TMath::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)- amp*t*t*TMath::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. | |
277 | } | |
278 | else{ | |
279 | Double_t det = b*b - a*c; | |
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 | } | |
291 | ||
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. ; | |
296 | for(Int_t i=minI; i<sigLength; i++){ | |
297 | if(samples.At(i)<=0.) | |
298 | continue ; | |
299 | Double_t dt1 = times.At(i) - t1; | |
300 | Double_t f01 = dt1*dt1*TMath::Exp(-2.*dt1); | |
301 | amp1 += f01*samples.At(i); | |
302 | den1 += f01*f01; | |
303 | } | |
304 | if(den1>0.0) amp1 /= den1; | |
305 | Double_t chi1=0.; // chi2=0. ; | |
306 | for(Int_t i=minI; i<sigLength; i++){ | |
307 | if(samples.At(i)<=0.) | |
308 | continue ; | |
309 | Double_t dt1 = times.At(i)- t1; | |
310 | Double_t dy1 = samples.At(i)- amp1*dt1*dt1*TMath::Exp(-2.*dt1) ; | |
311 | chi1 += dy1*dy1; | |
312 | } | |
313 | fEnergy=amp1*TMath::Exp(-2.) ; ; | |
314 | fTime+=t1*tau ; | |
315 | fQuality=chi1/sigLength ; | |
316 | } | |
317 | ||
318 | //Impose cut on quality | |
319 | fQuality/=2.+0.004*fEnergy*fEnergy ; | |
320 | ||
321 | //Draw corrupted samples | |
322 | if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){ | |
323 | plot: | |
324 | if(fEnergy > 30. && fQuality >1. && !fOverflow ){ //Draw only bad samples | |
325 | // if(!fOverflow ){ //Draw only bad samples | |
326 | printf("Sample par: amp=%f, t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ; | |
327 | TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ; | |
328 | if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ; | |
329 | h->Reset() ; | |
330 | for (Int_t i=0; i<sigLength; i++) { | |
331 | h->SetBinContent(i+1,samples.At(i)+pedestal) ; | |
332 | } | |
333 | TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ; | |
334 | fffit->SetParameters(pedestal,fEnergy,fTime,tau) ; | |
335 | fffit->SetLineColor(2) ; | |
336 | TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ; | |
337 | if(!can){ | |
338 | can = new TCanvas("cSamples","cSamples",10,10,600,600) ; | |
339 | can->SetFillColor(0) ; | |
340 | can->SetFillStyle(0) ; | |
341 | can->Range(0,0,1,1); | |
342 | can->SetBorderSize(0); | |
343 | } | |
344 | can->cd() ; | |
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] ; | |
359 | snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ; | |
360 | h->SetTitle(title) ; | |
361 | h->Draw() ; | |
362 | fffit->Draw("same") ; | |
363 | ||
364 | snprintf(title,155,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ; | |
365 | TFile fout("samples_bad.root","update") ; | |
366 | h->Write(title); | |
367 | fout.Close() ; | |
368 | ||
369 | can->cd() ; | |
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") ; | |
385 | if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ; | |
386 | hd->Reset() ; | |
387 | for (Int_t i=0; i<sigLength; i++) { | |
388 | hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ; | |
389 | } | |
390 | hd->Draw(); | |
391 | /* | |
392 | can->Update() ; | |
393 | printf("Press <enter> to continue\n") ; | |
394 | getchar(); | |
395 | */ | |
396 | ||
397 | delete fffit ; | |
398 | delete spectrum_1 ; | |
399 | delete spectrum_2 ; | |
400 | } | |
401 | } | |
402 | ||
403 | return kTRUE; | |
404 | } |