]>
Commit | Line | Data |
---|---|---|
de0cfd46 | 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 | ||
17 | // This class extracts the signal parameters (energy, time, quality) | |
18 | // from ALTRO samples. Energy is in ADC counts, time is in time bin units. | |
19 | // If sample is not in saturation, a coarse algorithm is used (a-la AliPHOSRawFitterv0) | |
20 | // If sample is in saturation, the unsaturated part of the sample is fit a-la AliPHOSRawFitterv1 | |
21 | // | |
22 | // AliPHOSRawFitterv4 *fitterv4=new AliPHOSRawFitterv4(); | |
23 | // fitterv4->SetChannelGeo(module,cellX,cellZ,caloFlag); | |
24 | // fitterv4->SetCalibData(fgCalibData) ; | |
25 | // fitterv4->Eval(sig,sigStart,sigLength); | |
26 | // Double_t amplitude = fitterv4.GetEnergy(); | |
27 | // Double_t time = fitterv4.GetTime(); | |
28 | // Bool_t isLowGain = fitterv4.GetCaloFlag()==0; | |
29 | ||
30 | // Author: Dmitri Peressounko | |
31 | ||
32 | // --- ROOT system --- | |
33 | #include "TArrayI.h" | |
34 | #include "TMath.h" | |
35 | #include "TObject.h" | |
36 | #include "TArrayD.h" | |
37 | #include "TList.h" | |
38 | #include "TMinuit.h" | |
39 | ||
40 | ||
41 | //Used for debug | |
42 | //#include "TROOT.h" | |
43 | //#include "TH1.h" | |
44 | //#include "TCanvas.h" | |
45 | //#include "TPad.h" | |
46 | //#include "TF1.h" | |
47 | ||
48 | ||
49 | ||
50 | // --- AliRoot header files --- | |
51 | #include "AliPHOSRawFitterv4.h" | |
52 | #include "AliPHOSCalibData.h" | |
53 | #include "AliLog.h" | |
54 | ||
55 | ClassImp(AliPHOSRawFitterv4) | |
56 | ||
57 | //----------------------------------------------------------------------------- | |
58 | AliPHOSRawFitterv4::AliPHOSRawFitterv4(): | |
59 | AliPHOSRawFitterv1(), | |
60 | fFitHighGain(0) | |
61 | { | |
62 | //Default constructor | |
63 | } | |
64 | ||
65 | //----------------------------------------------------------------------------- | |
66 | AliPHOSRawFitterv4::~AliPHOSRawFitterv4() | |
67 | { | |
68 | } | |
69 | ||
de0cfd46 | 70 | //----------------------------------------------------------------------------- |
71 | ||
72 | Bool_t AliPHOSRawFitterv4::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength) | |
73 | { | |
74 | // Calculate signal parameters (energy, time, quality) from array of samples | |
75 | // Energy is a maximum sample minus pedestal 9 | |
76 | // Time is the first time bin | |
77 | // Signal overflows is there are at least 3 samples of the same amplitude above 900 | |
78 | ||
79 | fOverflow= kFALSE ; | |
80 | fEnergy = 0; | |
81 | if (fNBunches > 1) { | |
82 | fQuality = 1000; | |
83 | return kTRUE; | |
84 | } | |
85 | ||
86 | const Float_t kBaseLine = 1.0; | |
87 | const Int_t kPreSamples = 10; | |
88 | ||
89 | Float_t pedMean = 0; | |
90 | Float_t pedRMS = 0; | |
91 | Int_t nPed = 0; | |
92 | UShort_t maxSample = 0; | |
93 | Int_t nMax = 0; | |
94 | ||
95 | for (Int_t i=0; i<sigLength; i++) { | |
96 | if (i>sigLength-kPreSamples) { //inverse signal time order | |
97 | nPed++; | |
98 | pedMean += signal[i]; | |
99 | pedRMS += signal[i]*signal[i] ; | |
100 | } | |
101 | if(signal[i] > maxSample ){ maxSample = signal[i]; nMax=0;} | |
102 | if(signal[i] == maxSample) nMax++; | |
103 | ||
104 | } | |
105 | ||
106 | fEnergy = (Double_t)maxSample; | |
107 | if (maxSample > 850 && nMax > 2) fOverflow = kTRUE; | |
108 | ||
109 | Double_t pedestal = 0 ; | |
110 | if (fPedSubtract) { | |
111 | if (nPed > 0) { | |
112 | fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ; | |
113 | if(fPedestalRMS > 0.) | |
114 | fPedestalRMS = TMath::Sqrt(fPedestalRMS) ; | |
115 | pedestal = (Double_t)(pedMean/nPed); | |
116 | } | |
117 | else | |
118 | return kFALSE; | |
119 | } | |
120 | else { | |
121 | //take pedestals from DB | |
122 | pedestal = (Double_t) fAmpOffset ; | |
123 | if (fCalibData) { | |
124 | Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ; | |
125 | Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ; | |
126 | pedestal += truePed - altroSettings ; | |
127 | } | |
128 | else{ | |
129 | AliDebug(2,Form("Pedestal and offset are not read from OCDB. Use 0 for their values.")) ; | |
130 | } | |
131 | } | |
132 | fEnergy-=pedestal ; | |
133 | if (fEnergy < kBaseLine) fEnergy = 0; | |
134 | ||
135 | ||
136 | if(fOverflow && ((fCaloFlag == 0) || fFitHighGain)){ //by default fit LowGain only | |
137 | TArrayI *samples = new TArrayI(sigLength); // array of sample amplitudes | |
138 | TArrayI *times = new TArrayI(sigLength); // array of sample time stamps | |
139 | Bool_t result = kTRUE ; | |
140 | for (Int_t i=0; i<sigLength; i++) { | |
141 | samples->AddAt(signal[i]-pedestal,sigLength-i-1); | |
142 | times ->AddAt(i ,i); | |
143 | } | |
144 | //Prepare fit parameters | |
145 | //Evaluate time | |
146 | Int_t iStart = 0; | |
147 | while(iStart<sigLength && samples->At(iStart) <kBaseLine) iStart++ ; | |
148 | if (fCaloFlag == 0){ // Low gain | |
149 | fSampleParamsLow->AddAt(pedestal,4) ; | |
150 | fSampleParamsLow->AddAt(double(maxSample),5) ; | |
151 | fSampleParamsLow->AddAt(double(iStart),6) ; | |
de0cfd46 | 152 | } |
153 | else if (fCaloFlag == 1){ // High gain | |
154 | fSampleParamsHigh->AddAt(pedestal,4) ; | |
155 | fSampleParamsHigh->AddAt(double(maxSample),5) ; | |
156 | fSampleParamsHigh->AddAt(double(iStart),6) ; | |
de0cfd46 | 157 | } |
158 | result=EvalWithFitting(samples,times); | |
06532a37 | 159 | fToFit->Clear("nodelete") ; |
de0cfd46 | 160 | delete samples ; |
161 | delete times ; | |
162 | ||
163 | return result; | |
164 | ||
165 | } | |
166 | ||
167 | ||
168 | //Sample withour overflow - Evaluate time | |
169 | fTime = sigStart-sigLength-3; | |
170 | const Int_t nLine= 6 ; //Parameters of fitting | |
171 | const Float_t eMinTOF = 10. ; //Choosed from beam-test and cosmic analyis | |
172 | const Float_t kAmp=0.35 ; //Result slightly depends on them, so no getters | |
173 | // Avoid too low peak: | |
174 | if(fEnergy < eMinTOF){ | |
175 | return kTRUE; | |
176 | } | |
177 | // Find index posK (kLevel is a level of "timestamp" point Tk): | |
178 | Int_t posK =sigLength-1 ; //last point before crossing k-level | |
179 | Double_t levelK = pedestal + kAmp*fEnergy; | |
180 | while(signal[posK] <= levelK && posK>=0){ | |
181 | posK-- ; | |
182 | } | |
183 | posK++ ; | |
184 | ||
185 | if(posK == 0 || posK==sigLength-1){ | |
186 | return kTRUE; | |
187 | } | |
188 | ||
189 | // Find crosing point by solving linear equation (least squares method) | |
190 | Int_t np = 0; | |
191 | Int_t iup=posK-1 ; | |
192 | Int_t idn=posK ; | |
193 | Double_t sx = 0., sy = 0., sxx = 0., sxy = 0.; | |
194 | Double_t x,y ; | |
195 | ||
196 | while(np<nLine){ | |
197 | //point above crossing point | |
198 | if(iup>=0){ | |
199 | x = sigLength-iup-1; | |
200 | y = signal[iup]; | |
201 | sx += x; | |
202 | sy += y; | |
203 | sxx += (x*x); | |
204 | sxy += (x*y); | |
205 | np++ ; | |
206 | iup-- ; | |
207 | } | |
208 | //Point below crossing point | |
209 | if(idn<sigLength){ | |
210 | if(signal[idn]<pedestal){ | |
211 | idn=sigLength-1 ; //do not scan further | |
212 | idn++ ; | |
213 | continue ; | |
214 | } | |
215 | x = sigLength-idn-1; | |
216 | y = signal[idn]; | |
217 | sx += x; | |
218 | sy += y; | |
219 | sxx += (x*x); | |
220 | sxy += (x*y); | |
221 | np++; | |
222 | idn++ ; | |
223 | } | |
224 | if(idn>=sigLength && iup<0){ | |
225 | break ; //can not fit futher | |
226 | } | |
227 | } | |
228 | ||
229 | Double_t det = np*sxx - sx*sx; | |
230 | if(det == 0){ | |
231 | return kTRUE; | |
232 | } | |
233 | Double_t c1 = (np*sxy - sx*sy)/det; //slope | |
eaf6dbb0 | 234 | Double_t c0 = 0; |
235 | if (np>0) | |
236 | c0 = (sy-c1*sx)/np; //offset | |
de0cfd46 | 237 | if(c1 == 0){ |
238 | return kTRUE; | |
239 | } | |
240 | ||
241 | // Find where the line cross kLevel: | |
242 | fTime += (levelK - c0)/c1-5. ; //5: mean offset between k-Level and start times | |
243 | return kTRUE; | |
244 | ||
245 | } | |
246 | //=================================================== | |
247 | Bool_t AliPHOSRawFitterv4::EvalWithFitting(TArrayI*samples, TArrayI * times){ | |
248 | ||
249 | // if sample has reasonable mean and RMS, try to fit it with gamma2 | |
250 | const Float_t sampleMaxHG=102.332 ; | |
251 | const Float_t sampleMaxLG=277.196 ; | |
252 | ||
253 | gMinuit->mncler(); // Reset Minuit's list of paramters | |
254 | gMinuit->SetPrintLevel(-1) ; // No Printout | |
255 | gMinuit->SetFCN(AliPHOSRawFitterv1::UnfoldingChiSquare) ; | |
256 | // To set the address of the minimization function | |
257 | ||
258 | fToFit->Clear("nodelete") ; | |
259 | Double_t b=0,bmin=0,bmax=0 ; | |
260 | if (fCaloFlag == 0){ // Low gain | |
261 | b=fSampleParamsLow->At(2) ; | |
262 | bmin=0.5 ; | |
263 | bmax=10. ; | |
264 | fToFit->AddFirst((TObject*)fSampleParamsLow) ; | |
265 | } | |
266 | else if (fCaloFlag == 1){ // High gain | |
267 | b=fSampleParamsHigh->At(2) ; | |
268 | bmin=0.05 ; | |
269 | bmax=0.4 ; | |
270 | fToFit->AddFirst((TObject*)fSampleParamsHigh) ; | |
271 | } | |
272 | fToFit->AddLast((TObject*)samples) ; | |
273 | fToFit->AddLast((TObject*)times) ; | |
274 | ||
275 | ||
276 | gMinuit->SetObjectFit((TObject*)fToFit) ; // To tranfer pointer to UnfoldingChiSquare | |
277 | Int_t ierflg ; | |
278 | gMinuit->mnparm(0, "t0", 1., 0.01, -50., 50., ierflg) ; | |
279 | if(ierflg != 0){ | |
280 | // AliWarning(Form("Unable to set initial value for fit procedure : t0=%e\n",1.*tStart) ) ; | |
281 | fEnergy = 0. ; | |
282 | fTime =-999. ; | |
283 | fQuality= 999. ; | |
284 | return kTRUE ; //will scan further | |
285 | } | |
286 | Double_t amp0=0; | |
287 | if (fCaloFlag == 0) // Low gain | |
288 | amp0 = fEnergy/sampleMaxLG; | |
289 | else if (fCaloFlag == 1) // High gain | |
290 | amp0 = fEnergy/sampleMaxHG; | |
291 | ||
292 | gMinuit->mnparm(1, "Energy", amp0 , 0.01*amp0, 0, 0, ierflg) ; | |
293 | if(ierflg != 0){ | |
294 | // AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ; | |
295 | fEnergy = 0. ; | |
296 | fTime =-999. ; | |
297 | fQuality= 999. ; | |
298 | return kTRUE ; //will scan further | |
299 | } | |
300 | ||
301 | gMinuit->mnparm(2, "p2", b, 0.01*b, bmin, bmax, ierflg) ; | |
302 | if(ierflg != 0){ | |
303 | // AliWarning(Form("Unable to set initial value for fit procedure : E=%e\n", amp0)) ; | |
304 | fEnergy = 0. ; | |
305 | fTime =-999. ; | |
306 | fQuality= 999. ; | |
307 | return kTRUE ; //will scan further | |
308 | } | |
309 | ||
310 | Double_t p0 = 0.0001 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly | |
311 | // depends on it. | |
312 | Double_t p1 = 1.0 ; | |
313 | Double_t p2 = 0.0 ; | |
314 | gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls | |
315 | gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient | |
316 | ////// gMinuit->SetMaxIterations(100); | |
317 | gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings | |
318 | ||
319 | gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize | |
320 | ||
321 | Double_t err=0.,t0err=0. ; | |
322 | Double_t t0=0.,efit=0. ; | |
323 | gMinuit->GetParameter(0,t0, t0err) ; | |
324 | gMinuit->GetParameter(1,efit, err) ; | |
325 | ||
326 | Double_t bfit=0., berr=0. ; | |
327 | gMinuit->GetParameter(2,bfit,berr) ; | |
328 | ||
329 | //Calculate total energy | |
330 | //this is parameterization of dependence of pulse height on parameter b | |
331 | if(fCaloFlag == 0) // Low gain | |
332 | fEnergy = efit*(99.54910 + 78.65038*bfit) ; | |
333 | else if(fCaloFlag == 1) // High gain | |
334 | fEnergy=efit*(80.33109 + 128.6433*bfit) ; | |
335 | ||
336 | if(fEnergy < 0. || fEnergy > 10000.){ | |
337 | //set energy to previously found max | |
338 | fTime =-999.; | |
339 | fQuality= 999 ; | |
340 | fEnergy=0. ; | |
341 | return kTRUE; | |
342 | } | |
343 | ||
344 | //evaluate fit quality | |
345 | Double_t fmin=0.,fedm=0.,errdef=0. ; | |
346 | Int_t npari,nparx,istat; | |
347 | gMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat) ; | |
348 | fQuality = fmin/samples->GetSize() ; | |
349 | //compare quality with some parameterization | |
350 | if (fCaloFlag == 0) // Low gain | |
351 | fQuality /= 2.00 + 0.0020*fEnergy ; | |
352 | else if (fCaloFlag == 1) // High gain | |
353 | fQuality /= 0.75 + 0.0025*fEnergy ; | |
354 | ||
355 | fTime += t0 - 4.024*bfit ; | |
356 | ||
357 | if(fQuality==0.){//no points to fit) | |
358 | fTime =-999.; | |
359 | fQuality= 1999 ; | |
360 | fEnergy=0. ; | |
361 | } | |
362 | ||
363 | /* | |
364 | if(1){ | |
365 | TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ; | |
366 | if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ; | |
367 | h->Reset() ; | |
368 | for (Int_t i=0; i<samples->GetSize(); i++) { | |
369 | h->SetBinContent(i+1,float(samples->At(i))) ; | |
370 | } | |
371 | TF1 * fffit = new TF1("fffit","[0]*(abs(x-[1])^[3]*exp(-[2]*(x-[1]))+[4]*(x-[1])*(x-[1])*exp(-[5]*(x-[1])))",0.,60.) ; | |
372 | TArrayD * params=(TArrayD*)fToFit->At(0) ; | |
373 | Double_t n=params->At(0) ; | |
374 | Double_t alpha=params->At(1) ; | |
375 | Double_t beta=params->At(3) ; | |
376 | fffit->SetParameters(efit,t0,alpha,n,bfit,beta) ; | |
377 | fffit->SetLineColor(2) ; | |
378 | TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ; | |
379 | if(!can){ | |
380 | can = new TCanvas("cSamples","cSamples",10,10,600,600) ; | |
381 | can->SetFillColor(0) ; | |
382 | can->SetFillStyle(0) ; | |
383 | can->Range(0,0,1,1); | |
384 | can->SetBorderSize(0); | |
385 | } | |
386 | can->cd() ; | |
387 | ||
388 | // TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99); | |
389 | TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.001,0.999,0.999); | |
390 | spectrum_1->Draw(); | |
391 | spectrum_1->cd(); | |
392 | spectrum_1->Range(0,0,1,1); | |
393 | spectrum_1->SetFillColor(0); | |
394 | spectrum_1->SetFillStyle(4000); | |
395 | spectrum_1->SetBorderSize(1); | |
396 | spectrum_1->SetBottomMargin(0.012); | |
397 | spectrum_1->SetTopMargin(0.03); | |
398 | spectrum_1->SetLeftMargin(0.10); | |
399 | spectrum_1->SetRightMargin(0.05); | |
400 | ||
401 | char title[155] ; | |
402 | snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ; | |
403 | h->SetTitle(title) ; | |
404 | // h->Fit(fffit,"","",0.,51.) ; | |
405 | h->Draw() ; | |
406 | fffit->Draw("same") ; | |
407 | ||
408 | ||
409 | // can->cd() ; | |
410 | // TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32); | |
411 | // spectrum_2->SetFillColor(0) ; | |
412 | // spectrum_2->SetFillStyle(0) ; | |
413 | // spectrum_2->SetGridy() ; | |
414 | // spectrum_2->Draw(); | |
415 | // spectrum_2->Range(0,0,1,1); | |
416 | // spectrum_2->SetFillColor(0); | |
417 | // spectrum_2->SetBorderSize(1); | |
418 | // spectrum_2->SetTopMargin(0.01); | |
419 | // spectrum_2->SetBottomMargin(0.25); | |
420 | // spectrum_2->SetLeftMargin(0.10); | |
421 | // spectrum_2->SetRightMargin(0.05); | |
422 | // spectrum_2->cd() ; | |
423 | // | |
424 | // TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ; | |
425 | // if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ; | |
426 | // hd->Reset() ; | |
427 | // for (Int_t i=0; i<samples->GetSize(); i++) { | |
428 | // hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples->At(i)-fffit->Eval(i)))) ; | |
429 | // } | |
430 | // hd->Draw(); | |
431 | ||
432 | can->Update() ; | |
433 | printf("Press <enter> to continue\n") ; | |
434 | getchar(); | |
435 | ||
436 | ||
437 | delete fffit ; | |
438 | delete spectrum_1 ; | |
439 | // delete spectrum_2 ; | |
440 | } | |
441 | */ | |
442 | ||
443 | return kTRUE; | |
444 | ||
445 | } |