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