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