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 "AliPHOSRawFitterv3.h"
49 : #include "AliPHOSPulseGenerator.h"
50 :
51 22 : ClassImp(AliPHOSRawFitterv3)
52 :
53 : //-----------------------------------------------------------------------------
54 : AliPHOSRawFitterv3::AliPHOSRawFitterv3():
55 0 : AliPHOSRawFitterv0()
56 0 : {
57 : //Default constructor.
58 0 : }
59 :
60 : //-----------------------------------------------------------------------------
61 : AliPHOSRawFitterv3::~AliPHOSRawFitterv3()
62 0 : {
63 : //Destructor.
64 0 : }
65 :
66 : //-----------------------------------------------------------------------------
67 : AliPHOSRawFitterv3::AliPHOSRawFitterv3(const AliPHOSRawFitterv3 &phosFitter ):
68 0 : AliPHOSRawFitterv0(phosFitter)
69 0 : {
70 : //Copy constructor.
71 0 : }
72 :
73 : //-----------------------------------------------------------------------------
74 : AliPHOSRawFitterv3& AliPHOSRawFitterv3::operator = (const AliPHOSRawFitterv3 & /*phosFitter*/)
75 : {
76 : //Assignment operator.
77 0 : 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 0 : if (fCaloFlag == 2 || fNBunches > 1) {
91 0 : fQuality = 1000;
92 0 : return kTRUE;
93 : }
94 0 : if(fCaloFlag!=0 && fCaloFlag!=1){//Corrupted sample
95 0 : fQuality=2000;
96 0 : fEnergy=0 ;
97 0 : 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.
110 : const Double_t tau=22.18 ;
111 :
112 0 : TArrayD samples(sigLength); // array of sample amplitudes
113 0 : TArrayD times(sigLength); // array of sample time stamps
114 0 : for (Int_t i=sigLength-kPreSamples; i<sigLength; i++) {
115 0 : nPed++;
116 0 : pedMean += signal[i];
117 0 : pedRMS += signal[i]*signal[i] ;
118 : }
119 :
120 0 : fEnergy = -111;
121 0 : 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 0 : if (fPedSubtract) {
126 0 : if (nPed > 0) {
127 0 : fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
128 0 : if(fPedestalRMS > 0.)
129 0 : fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
130 0 : pedestal = (Double_t)(pedMean/nPed); // pedestal subtraction
131 0 : fEnergy -= pedestal; // pedestal subtraction
132 : }
133 : else
134 0 : return kFALSE;
135 0 : }
136 : else {
137 : //take pedestals from DB
138 0 : pedestal = (Double_t) fAmpOffset ;
139 0 : if (fCalibData) {
140 0 : Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
141 0 : Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
142 0 : pedestal += truePed - altroSettings ;
143 0 : }
144 : else{
145 : // AliWarning(Form("Can not read data from OCDB")) ;
146 : }
147 0 : fEnergy-=pedestal ;
148 : }
149 :
150 0 : if (fEnergy < kBaseLine) fEnergy = 0;
151 : //Evaluate time
152 0 : fTime = sigStart-sigLength-3;
153 0 : for (Int_t i=0; i<sigLength; i++) {
154 0 : samples.AddAt(signal[i]-pedestal,sigLength-i-1);
155 0 : times.AddAt(i/tau ,i);
156 : }
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 ;
165 0 : for (Int_t i=0; i<sigLength; i++){
166 0 : if(signal[i] > pedestal){
167 0 : Double_t de = signal[i] - pedestal ;
168 0 : if(de > 1.) {
169 0 : aMean += de*i ;
170 0 : aRMS += de*i*i ;
171 0 : wts += de;
172 0 : }
173 0 : if(signal[i] > maxAmp){
174 : maxAmp = signal[i];
175 : nMax=0;
176 : maxBin = i ;
177 0 : }
178 0 : if(signal[i] == maxAmp){
179 0 : nMax++;
180 0 : }
181 0 : }
182 : }
183 :
184 0 : if (maxBin==sigLength-1){//bad "rising" sample
185 0 : fEnergy = 0. ;
186 0 : fTime = -999. ;
187 0 : fQuality= 999. ;
188 0 : return kTRUE ;
189 : }
190 :
191 0 : fEnergy=Double_t(maxAmp)-pedestal ;
192 0 : fOverflow =0 ; //look for plato on the top of sample
193 0 : if (fEnergy>500 && //this is not fluctuation of soft sample
194 0 : nMax>2){ //and there is a plato
195 0 : fOverflow = kTRUE ;
196 0 : }
197 :
198 0 : if (wts > 0) {
199 0 : aMean /= wts;
200 0 : aRMS = aRMS/wts - aMean*aMean;
201 0 : }
202 :
203 : //do not take too small energies
204 0 : if (fEnergy < kBaseLine)
205 0 : fEnergy = 0;
206 :
207 : //do not test quality of too soft samples
208 0 : if (fEnergy < maxEtoFit){
209 0 : if (aRMS < 2.) //sigle peak
210 0 : fQuality = 999. ;
211 : else
212 0 : fQuality = 0. ;
213 : //Evaluate time of signal arriving
214 0 : return kTRUE ;
215 : }
216 :
217 : // if sample has reasonable mean and RMS, try to fit it with gamma2
218 : //This method can not analyse overflow samples
219 0 : if(fOverflow){
220 0 : fQuality = 99. ;
221 0 : return kTRUE ;
222 : }
223 : // First estimate t0
224 : Double_t a=0,b=0,c=0 ;
225 : Int_t minI=0 ;
226 0 : if (fPedSubtract)
227 : minI=kPreSamples ;
228 0 : for(Int_t i=minI; i<sigLength; i++){
229 0 : if(samples.At(i)<=0.)
230 : continue ;
231 0 : Double_t t= times.At(i) ;
232 0 : Double_t f02 = TMath::Exp(-2.*t);
233 0 : Double_t f12 = t*f02;
234 0 : Double_t f22 = t*f12;
235 : // Derivatives
236 0 : Double_t f02d = -2.*f02;
237 0 : Double_t f12d = f02 - 2.*f12;
238 0 : Double_t f22d = 2.*(f12 - f22);
239 0 : a += f02d * samples.At(i) ;
240 0 : b -= f12d * samples.At(i) ;
241 0 : c += f22d * samples.At(i) ;
242 0 : }
243 :
244 : //Find roots
245 0 : if(a==0.){
246 0 : if(b==0.){ //no roots
247 0 : fQuality = 2000. ;
248 0 : if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
249 0 : printf(" a=%f, b=%f, c=%f \n",a,b,c) ;
250 : goto plot ;
251 : }
252 0 : return kTRUE ;
253 : }
254 0 : Double_t t1=-c/b ;
255 : Double_t amp=0.,den=0.; ;
256 0 : for(Int_t i=minI; i<sigLength; i++){
257 0 : if(samples.At(i)<=0.)
258 : continue ;
259 0 : Double_t dt = times.At(i) - t1;
260 0 : Double_t f = dt*dt*TMath::Exp(-2.*dt);
261 0 : amp += f*samples.At(i);
262 0 : den += f*f;
263 0 : }
264 0 : if(den>0.0) amp /= den;
265 : // chi2 calculation
266 0 : fQuality=0.;
267 0 : for(Int_t i=minI; i<sigLength; i++){
268 0 : if(samples.At(i)<=0.)
269 : continue ;
270 0 : Double_t t = times.At(i)- t1;
271 0 : Double_t dy = samples.At(i)- amp*t*t*TMath::Exp(-2.*t) ;
272 0 : fQuality += dy*dy;
273 0 : }
274 0 : fTime+=t1*tau ;
275 0 : fEnergy = amp*TMath::Exp(-2.);
276 0 : fQuality/= sigLength ; //If we have overflow the number of actually fitted points is smaller, but chi2 in this case is not important.
277 0 : }
278 : else{
279 0 : Double_t det = b*b - a*c;
280 0 : if(det>=1.e-6 && det<0.0) {
281 : det = 0.0; //rounding error
282 : }
283 0 : if(det<0.){ //Problem
284 0 : fQuality = 1500. ;
285 0 : if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
286 0 : printf(" det=%e \n",det) ;
287 0 : goto plot ;
288 : }
289 0 : return kTRUE ;
290 : }
291 :
292 0 : det = TMath::Sqrt(det);
293 0 : 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. ;
296 0 : for(Int_t i=minI; i<sigLength; i++){
297 0 : if(samples.At(i)<=0.)
298 : continue ;
299 0 : Double_t dt1 = times.At(i) - t1;
300 0 : Double_t f01 = dt1*dt1*TMath::Exp(-2.*dt1);
301 0 : amp1 += f01*samples.At(i);
302 0 : den1 += f01*f01;
303 0 : }
304 0 : if(den1>0.0) amp1 /= den1;
305 : Double_t chi1=0.; // chi2=0. ;
306 0 : for(Int_t i=minI; i<sigLength; i++){
307 0 : if(samples.At(i)<=0.)
308 : continue ;
309 0 : Double_t dt1 = times.At(i)- t1;
310 0 : Double_t dy1 = samples.At(i)- amp1*dt1*dt1*TMath::Exp(-2.*dt1) ;
311 0 : chi1 += dy1*dy1;
312 0 : }
313 0 : fEnergy=amp1*TMath::Exp(-2.) ; ;
314 0 : fTime+=t1*tau ;
315 0 : fQuality=chi1/sigLength ;
316 0 : }
317 :
318 : //Impose cut on quality
319 0 : fQuality/=2.+0.004*fEnergy*fEnergy ;
320 :
321 : //Draw corrupted samples
322 0 : if(AliLog::GetDebugLevel("PHOS","AliPHOSRawFitterv3")>3){
323 : plot:
324 0 : if(fEnergy > 30. && fQuality >1. && !fOverflow ){ //Draw only bad samples
325 : // if(!fOverflow ){ //Draw only bad samples
326 0 : printf("Sample par: amp=%f, t0=%f, Quality=%f \n",fEnergy,fTime,fQuality) ;
327 0 : TH1I * h = (TH1I*)gROOT->FindObjectAny("hSamples") ;
328 0 : if(!h) h = new TH1I("hSamples","Samples",65,0.,65.) ;
329 0 : h->Reset() ;
330 0 : for (Int_t i=0; i<sigLength; i++) {
331 0 : h->SetBinContent(i+1,samples.At(i)+pedestal) ;
332 : }
333 0 : TF1 * fffit = new TF1("fffit","[0]+[1]*((x-[2])/[3])^2*exp(2.-2.*(x-[2])/[3])",0.,200.) ;
334 0 : fffit->SetParameters(pedestal,fEnergy,fTime,tau) ;
335 0 : fffit->SetLineColor(2) ;
336 0 : TCanvas * can = (TCanvas*)gROOT->FindObjectAny("cSamples") ;
337 0 : if(!can){
338 0 : can = new TCanvas("cSamples","cSamples",10,10,600,600) ;
339 0 : can->SetFillColor(0) ;
340 0 : can->SetFillStyle(0) ;
341 0 : can->Range(0,0,1,1);
342 0 : can->SetBorderSize(0);
343 : }
344 0 : can->cd() ;
345 :
346 0 : TPad * spectrum_1 = new TPad("spectrum_1", "spectrum_1",0.001,0.32,0.99,0.99);
347 0 : spectrum_1->Draw();
348 0 : spectrum_1->cd();
349 0 : spectrum_1->Range(0,0,1,1);
350 0 : spectrum_1->SetFillColor(0);
351 0 : spectrum_1->SetFillStyle(4000);
352 0 : spectrum_1->SetBorderSize(1);
353 0 : spectrum_1->SetBottomMargin(0.012);
354 0 : spectrum_1->SetTopMargin(0.03);
355 0 : spectrum_1->SetLeftMargin(0.10);
356 0 : spectrum_1->SetRightMargin(0.05);
357 :
358 0 : char title[155] ;
359 0 : snprintf(title,155,"Sample, mod=%d, x=%d, z=%d, Quality=%5.1f",fModule,fCellX,fCellZ,fQuality) ;
360 0 : h->SetTitle(title) ;
361 0 : h->Draw() ;
362 0 : fffit->Draw("same") ;
363 :
364 0 : snprintf(title,155,"mod%d_x%d_z%d_HG_qu%4.1f",fModule,fCellX,fCellZ,fQuality) ;
365 0 : TFile fout("samples_bad.root","update") ;
366 0 : h->Write(title);
367 0 : fout.Close() ;
368 :
369 0 : can->cd() ;
370 0 : TPad *spectrum_2 = new TPad("spectrum_2", "spectrum_2",0.001,0.01,0.99,0.32);
371 0 : spectrum_2->SetFillColor(0) ;
372 0 : spectrum_2->SetFillStyle(0) ;
373 0 : spectrum_2->SetGridy() ;
374 0 : spectrum_2->Draw();
375 0 : spectrum_2->Range(0,0,1,1);
376 0 : spectrum_2->SetFillColor(0);
377 0 : spectrum_2->SetBorderSize(1);
378 0 : spectrum_2->SetTopMargin(0.01);
379 0 : spectrum_2->SetBottomMargin(0.25);
380 0 : spectrum_2->SetLeftMargin(0.10);
381 0 : spectrum_2->SetRightMargin(0.05);
382 0 : spectrum_2->cd() ;
383 :
384 0 : TH1I * hd = (TH1I*)gROOT->FindObjectAny("hSamplesDif") ;
385 0 : if(!hd) hd = new TH1I("hd","Samples",65,0.,65.) ;
386 0 : hd->Reset() ;
387 0 : for (Int_t i=0; i<sigLength; i++) {
388 0 : hd->SetBinContent(i+1,TMath::Max(-1023.,TMath::Min(1023.,samples.At(i)+pedestal-fffit->Eval(i)))) ;
389 : }
390 0 : hd->Draw();
391 : /*
392 : can->Update() ;
393 : printf("Press <enter> to continue\n") ;
394 : getchar();
395 : */
396 :
397 0 : delete fffit ;
398 0 : delete spectrum_1 ;
399 0 : delete spectrum_2 ;
400 0 : }
401 : }
402 :
403 0 : return kTRUE;
404 0 : }
|
