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