]>
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 |
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 |
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 += de*i*i ; | |
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 |
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 | } |
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 = 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. | |
c4cc99c3 | 277 | } |
a25f5856 | 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 | } | |
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++){ |
297 | if(samples.At(i)<=0.) | |
facd0d71 | 298 | continue ; |
ad308f96 | 299 | Double_t dt1 = times.At(i) - t1; |
facd0d71 | 300 | Double_t f01 = dt1*dt1*TMath::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++){ |
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)- amp1*dt1*dt1*TMath::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.004*fEnergy*fEnergy ; | |
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 |
325 | // if(!fOverflow ){ //Draw only bad samples | |
c4cc99c3 | 326 | printf("Sample par: amp=%f, t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ; |
327 | TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ; | |
a25f5856 | 328 | if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ; |
c4cc99c3 | 329 | h->Reset() ; |
330 | for (Int_t i=0; i<sigLength; i++) { | |
a25f5856 | 331 | h->SetBinContent(i+1,samples.At(i)+pedestal) ; |
c4cc99c3 | 332 | } |
333 | TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(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") ; |
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() ; |
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 | } |
390 | hd->Draw(); | |
a25f5856 | 391 | /* |
facd0d71 | 392 | can->Update() ; |
c4cc99c3 | 393 | printf("Press <enter> to continue\n") ; |
394 | getchar(); | |
a25f5856 | 395 | */ |
c4cc99c3 | 396 | |
397 | delete fffit ; | |
398 | delete spectrum_1 ; | |
399 | delete spectrum_2 ; | |
400 | } | |
401 | } | |
402 | ||
c4cc99c3 | 403 | return kTRUE; |
404 | } |