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 : /* $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 "AliPHOSRawFitterv2.h"
49 : #include "AliPHOSPulseGenerator.h"
50 :
51 22 : ClassImp(AliPHOSRawFitterv2)
52 :
53 : //-----------------------------------------------------------------------------
54 : AliPHOSRawFitterv2::AliPHOSRawFitterv2():
55 0 : AliPHOSRawFitterv0(),
56 0 : fAlpha(0.1),fBeta(0.035),fMax(0)
57 0 : {
58 : //Default constructor.
59 0 : }
60 :
61 : //-----------------------------------------------------------------------------
62 : AliPHOSRawFitterv2::~AliPHOSRawFitterv2()
63 0 : {
64 : //Destructor.
65 0 : }
66 :
67 : //-----------------------------------------------------------------------------
68 : AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ):
69 0 : AliPHOSRawFitterv0(phosFitter),
70 0 : fAlpha(0.1),fBeta(0.035),fMax(0)
71 0 : {
72 : //Copy constructor.
73 0 : }
74 :
75 : //-----------------------------------------------------------------------------
76 : AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 & /*phosFitter*/)
77 : {
78 : //Assignment operator.
79 0 : return *this;
80 : }
81 :
82 : //-----------------------------------------------------------------------------
83 : Bool_t AliPHOSRawFitterv2::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
84 : {
85 : //Extract an energy deposited in the crystal,
86 : //crystal' position (module,column,row),
87 : //time and gain (high or low).
88 : //First collects sample, then evaluates it and if it has
89 : //reasonable shape, fits it with Gamma2 function and extracts
90 : //energy and time.
91 :
92 : const Float_t maxEtoFit=5 ; //fit only samples above this energy, accept all samples (with good aRMS) below it
93 : const Float_t kBaseLine = 1.0;
94 : const Int_t kPreSamples = 10;
95 :
96 0 : fOverflow = kFALSE ;
97 0 : fEnergy=0 ;
98 0 : if (fCaloFlag == 2 || fNBunches > 1) {
99 0 : fQuality = 150;
100 0 : return kTRUE;
101 : }
102 0 : if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
103 0 : fQuality=200;
104 0 : fEnergy=0 ;
105 0 : return kTRUE;
106 : }
107 :
108 : //Evaluate pedestals
109 : Float_t pedMean = 0;
110 : Float_t pedRMS = 0;
111 : Int_t nPed = 0;
112 0 : for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) {
113 0 : nPed++;
114 0 : pedMean += signal[i];
115 0 : pedRMS += signal[i]*signal[i] ;
116 : }
117 :
118 0 : fEnergy = -111;
119 0 : fQuality= 999. ;
120 : Double_t pedestal = 0;
121 :
122 0 : if (fPedSubtract) {
123 0 : if (nPed > 0) {
124 0 : fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
125 0 : if(fPedestalRMS > 0.)
126 0 : fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
127 0 : pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction
128 : }
129 : else
130 0 : return kFALSE;
131 0 : }
132 : else {
133 : //take pedestals from DB
134 0 : pedestal = (Double_t) fAmpOffset ;
135 : }
136 :
137 :
138 : //calculate rough quality of the sample and check for overflow
139 : Int_t maxAmp=0 ;
140 0 : Int_t minAmp= signal[0] ;
141 : Int_t nMax = 0 ; //number of points in plato
142 : Double_t aMean =0. ;
143 : Double_t aRMS =0. ;
144 : Double_t wts =0 ;
145 : Bool_t falling = kTRUE ; //Bad monotoneusly falling sample
146 : Bool_t rising = kTRUE ; //Bad monotoneusly riging sample
147 0 : for (Int_t i=0; i<sigLength; i++){
148 0 : if(signal[i] > pedestal){
149 0 : Double_t de = signal[i] - pedestal ;
150 0 : if(de > 1.) {
151 0 : aMean += de*i ;
152 0 : aRMS += de*i*i ;
153 0 : wts += de;
154 0 : }
155 0 : if(signal[i] > maxAmp){
156 : maxAmp = signal[i];
157 : nMax=0;
158 0 : }
159 0 : if(signal[i] == maxAmp){
160 0 : nMax++;
161 0 : }
162 0 : if(signal[i] < minAmp)
163 0 : minAmp=signal[i] ;
164 0 : if(falling && i>0 && signal[i]<signal[i-1])
165 0 : falling=kFALSE ;
166 0 : if(rising && i>0 && signal[i]>signal[i-1])
167 0 : rising=kFALSE ;
168 0 : }
169 : }
170 :
171 0 : if(rising || falling){//bad "rising" or falling sample
172 0 : fEnergy = 0. ;
173 0 : fTime = 0. ; //-999. ;
174 0 : fQuality= 250. ;
175 0 : return kTRUE ;
176 : }
177 0 : if(maxAmp-minAmp<3 && maxAmp>7 && sigLength>20){ //bad flat sample
178 0 : fEnergy = 0. ;
179 0 : fTime = 0; //-999. ;
180 0 : fQuality= 260. ;
181 0 : return kTRUE ;
182 : }
183 :
184 0 : fEnergy=Double_t(maxAmp)-pedestal ;
185 0 : if (fEnergy < kBaseLine) fEnergy = 0;
186 0 : fTime = sigStart-sigLength-3;
187 :
188 : //do not test quality of too soft samples
189 0 : if (wts > 0) {
190 0 : aMean /= wts;
191 0 : aRMS = aRMS/wts - aMean*aMean;
192 0 : }
193 0 : if (fEnergy <= maxEtoFit){
194 0 : if (aRMS < 2.) //sigle peak
195 0 : fQuality = 299. ;
196 : else
197 0 : fQuality = 0. ;
198 : //Evaluate time of signal arriving
199 0 : return kTRUE ;
200 : }
201 :
202 : //look for plato on the top of sample
203 0 : if (fEnergy>500 && //this is not fluctuation of soft sample
204 0 : nMax>2){ //and there is a plato
205 0 : fOverflow = kTRUE ;
206 0 : }
207 :
208 :
209 : //do not fit High Gain samples with overflow
210 0 : if(fCaloFlag==1 && fOverflow){
211 0 : fQuality = 99. ;
212 0 : return kTRUE;
213 :
214 : }
215 :
216 : //----Now fit sample with reasonable shape------
217 0 : TArrayD samples(sigLength); // array of sample amplitudes
218 0 : TArrayD times(sigLength); // array of sample time stamps
219 0 : for (Int_t i=0; i<sigLength; i++) {
220 0 : samples.AddAt(signal[i]-pedestal,sigLength-i-1);
221 0 : times.AddAt(double(i),i);
222 : }
223 :
224 0 : if(fMax==0)
225 0 : FindMax() ;
226 0 : if(!FindAmpT(samples,times)){
227 0 : if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv2")>3){
228 : goto plot ;
229 : }
230 : else{
231 0 : return kFALSE ;
232 : }
233 : }
234 0 : fEnergy*=fMax ;
235 0 : fTime += sigStart-sigLength-3;
236 :
237 :
238 : //Impose cut on quality
239 : // fQuality/=4. ;
240 0 : fQuality/=1.+0.005*fEnergy ;
241 :
242 : //Draw corrupted samples
243 0 : if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv2")>3){
244 0 : if(fEnergy > 50. ){
245 : plot:
246 0 : printf("Sample par: amp=%f, t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
247 0 : TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
248 0 : if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
249 0 : h->Reset() ;
250 0 : for (Int_t i=0; i<sigLength; i++) {
251 0 : h->SetBinContent(i+1,float(samples.At(i))) ;
252 : }
253 : // TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
254 0 : TF1 * fffit = new TF1("fffit","[0]*((x-[1])*(x-[1])*exp(-[2]*(x-[1]))+(x-[1])*exp(-[3]*(x-[1])))",0.,60.) ;
255 0 : fffit->SetParameters(fEnergy/fMax,fTime-(sigStart-sigLength-3),fAlpha,fBeta) ;
256 0 : fffit->SetLineColor(2) ;
257 0 : TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
258 0 : if(!can){
259 0 : can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
260 0 : can->SetFillColor(0) ;
261 0 : can->SetFillStyle(0) ;
262 0 : can->Range(0,0,1,1);
263 0 : can->SetBorderSize(0);
264 : }
265 0 : can->cd() ;
266 :
267 0 : TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
268 0 : spectrum_1->Draw();
269 0 : spectrum_1->cd();
270 0 : spectrum_1->Range(0,0,1,1);
271 0 : spectrum_1->SetFillColor(0);
272 0 : spectrum_1->SetFillStyle(4000);
273 0 : spectrum_1->SetBorderSize(1);
274 0 : spectrum_1->SetBottomMargin(0.012);
275 0 : spectrum_1->SetTopMargin(0.03);
276 0 : spectrum_1->SetLeftMargin(0.10);
277 0 : spectrum_1->SetRightMargin(0.05);
278 :
279 0 : char title[155] ;
280 0 : snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
281 0 : h->SetTitle(title) ;
282 : // h->Fit(fffit,"","",0.,51.) ;
283 0 : h->Draw() ;
284 0 : fffit->Draw("same") ;
285 : /*
286 : sprintf(title,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ;
287 : TFile fout("samples_bad.root","update") ;
288 : h->Write(title);
289 : fout.Close() ;
290 : */
291 0 : can->cd() ;
292 0 : TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
293 0 : spectrum_2->SetFillColor(0) ;
294 0 : spectrum_2->SetFillStyle(0) ;
295 0 : spectrum_2->SetGridy() ;
296 0 : spectrum_2->Draw();
297 0 : spectrum_2->Range(0,0,1,1);
298 0 : spectrum_2->SetFillColor(0);
299 0 : spectrum_2->SetBorderSize(1);
300 0 : spectrum_2->SetTopMargin(0.01);
301 0 : spectrum_2->SetBottomMargin(0.25);
302 0 : spectrum_2->SetLeftMargin(0.10);
303 0 : spectrum_2->SetRightMargin(0.05);
304 0 : spectrum_2->cd() ;
305 :
306 0 : TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
307 0 : if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
308 0 : hd->Reset() ;
309 0 : for (Int_t i=0; i<sigLength; i++) {
310 0 : hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ;
311 : }
312 0 : hd->Draw();
313 :
314 0 : can->Update() ;
315 0 : printf("Press <enter> to continue\n") ;
316 0 : getchar();
317 :
318 :
319 0 : delete fffit ;
320 0 : delete spectrum_1 ;
321 0 : delete spectrum_2 ;
322 0 : }
323 : }
324 :
325 0 : return kTRUE;
326 0 : }
327 : //------------------------------------------------------------------
328 : Bool_t AliPHOSRawFitterv2::FindAmpT(TArrayD samples, TArrayD times){
329 : // makes fit
330 :
331 : const Int_t nMaxIter=50 ; //Maximal number of iterations
332 : const Double_t epsdt = 1.e-3 ; //expected precision of t0 calculation
333 :
334 0 : Double_t dTime=times.At(0)-0.5 ; //Most probable Initial approximation
335 : //printf(" start fit... \n") ;
336 :
337 0 : Int_t nPoints = samples.GetSize() ;
338 0 : Double_t dea=TMath::Exp(-fAlpha) ;
339 0 : Double_t deb=TMath::Exp(-fBeta) ;
340 : Double_t dt=1.,timeOld=dTime,dfOld=0. ;
341 0 : for(Int_t iter=0; iter<nMaxIter; iter++){
342 : Double_t yy=0.;
343 : Double_t yf=0. ;
344 : Double_t ydf=0. ;
345 : Double_t yddf=0. ;
346 : Double_t ff=0. ;
347 : Double_t fdf=0. ;
348 : Double_t dfdf=0. ;
349 : Double_t fddf=0. ;
350 : Int_t nfit=0 ;
351 0 : Double_t aexp=TMath::Exp(-fAlpha*(times.At(0)-1.-dTime)) ;
352 0 : Double_t bexp=TMath::Exp(-fBeta*(times.At(0)-1.-dTime)) ;
353 0 : for(Int_t i=0; i<nPoints; i++){
354 0 : Double_t t= times.At(i)-dTime ;
355 0 : aexp*=dea ;
356 0 : bexp*=deb ;
357 0 : if(t<0.) continue ;
358 0 : Double_t y=samples.At(i) ;
359 0 : if(y<=fAmpThreshold)
360 0 : continue ;
361 0 : nfit++ ;
362 0 : Double_t at=fAlpha*t ;
363 0 : Double_t bt = fBeta*t ;
364 0 : Double_t phi=t*(t*aexp+bexp) ;
365 0 : Double_t dphi=t*aexp*(2.-at)+(1.-bt)*bexp ;
366 0 : Double_t ddphi=aexp*(2.-4.*at+at*at)+bexp*fBeta*(bt-2.) ;
367 0 : yy+=y*y ;
368 0 : yf+=y*phi ;
369 0 : ydf+=y*dphi ;
370 0 : yddf+=y*ddphi ;
371 0 : ff+=phi*phi ;
372 0 : fdf+=phi*dphi ;
373 0 : dfdf+=dphi*dphi ;
374 0 : fddf+=phi*ddphi ;
375 0 : }
376 :
377 0 : if(ff<1.e-09||nfit==0 ){
378 0 : fQuality=199 ;
379 0 : return kFALSE ;
380 : }
381 0 : Double_t f=ydf*ff-yf*fdf ; //d(chi2)/dt
382 0 : Double_t df=yf*(dfdf+fddf)-yddf*ff-ydf*fdf;
383 0 : if(df<=0.){ //we are too far from the root. In the wicinity of root df>0
384 0 : if(iter!=0 && dfOld>0.){//If at previous step df was OK, just reduce step size
385 0 : dt*=0.5 ;
386 0 : dTime=timeOld+dt ;
387 0 : continue ;
388 : }
389 0 : if(f<0){ //f<0 => dTime is too small and we still do not know root region
390 0 : dTime+=2. ;
391 0 : continue ;
392 : }
393 : else{ //dTime is too large, we are beyond the root region
394 0 : dTime-=2. ;
395 0 : continue ;
396 : }
397 : }
398 0 : dt=-f/df ;
399 0 : if(TMath::Abs(dt)<epsdt){
400 0 : fQuality=(yy-yf*yf/ff)/nfit ;
401 0 : fEnergy=yf/ff ; //ff!=0 already tested
402 0 : fTime=dTime ;
403 0 : return kTRUE ;
404 : }
405 : //In some cases time steps are huge (derivative ~0)
406 0 : if(dt>10.) dt=10. ; //restrict step size
407 0 : if(dt<-10.) dt=-5.3 ; //this restriction should be asimmetric to avoid jumping from one point to another
408 : timeOld=dTime ; //remember current position for the case
409 : dfOld=df ; //of reduction of dt step size
410 0 : dTime+=dt ;
411 :
412 0 : if(dTime>100. || dTime<-30.){ //this is corrupted sample, do not spend time improving accuracy.
413 0 : fQuality=(yy-yf*yf/ff)/nfit ;
414 0 : fEnergy=yf/ff ; //ff!=0 already tested
415 0 : fTime=dTime ;
416 0 : return kFALSE ;
417 : }
418 :
419 0 : }
420 : //failed to find a root, too many iterations
421 0 : fQuality=99.;
422 0 : fEnergy=0 ;
423 0 : return kFALSE ;
424 0 : }
425 : //_________________________________________
426 : void AliPHOSRawFitterv2::FindMax(){
427 : //Finds maxumum of currecnt parameterization
428 0 : Double_t t=2./fAlpha ;
429 0 : fMax = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;
430 : Double_t dt=15 ;
431 0 : while(dt>0.01){
432 0 : Double_t dfdt=(2.*t-fAlpha*t*t)*TMath::Exp(-fAlpha*t)+(1.-fBeta*t)*TMath::Exp(-fBeta*t) ;
433 0 : if(dfdt>0.)
434 0 : t+=dt ;
435 : else
436 0 : t-=dt ;
437 0 : Double_t maxNew = t*t*TMath::Exp(-fAlpha*t)+t*TMath::Exp(-fBeta*t) ;
438 0 : if(maxNew>fMax)
439 0 : fMax=maxNew ;
440 : else{
441 0 : dt/=2 ;
442 0 : if(dfdt<0.)
443 0 : t+=dt ;
444 : else
445 0 : t-=dt ;
446 : }
447 : }
448 0 : }
449 :
|
