Raw data decoder is migrated from AliCaloRawStream to AliCaloRawStreamV3
[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawFitterv2.cxx
CommitLineData
379c5c09 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 plots samples and qualify their quality according to their shape
19//
20// Typical use case:
21// AliPHOSRawFitter *fitter=new AliPHOSRawFitter();
22// fitter->SetSamples(sig,sigStart,sigLength);
23// fitter->SetNBunches(nBunches);
24// fitter->SetChannelGeo(module,cellX,cellZ,caloFlag);
25// fitter->SetCalibData(fgCalibData) ;
26// fitter->Eval();
27// Double_t amplitude = fitter.GetEnergy();
28// Double_t time = fitter.GetTime();
29// Bool_t isLowGain = fitter.GetCaloFlag()==0;
30
31// Author: Dmitri Peressounko (Oct.2008)
32// Modified: Yuri Kharlov (Jul.2009)
33
34// --- ROOT system ---
35#include "TList.h"
36#include "TMath.h"
37#include "TMinuit.h"
38#include "TCanvas.h"
39#include "TH1.h"
40#include "TH2.h"
41#include "TF1.h"
42#include "TROOT.h"
43
44// --- AliRoot header files ---
45#include "AliLog.h"
46#include "AliPHOSCalibData.h"
47#include "AliPHOSRawFitterv2.h"
48#include "AliPHOSPulseGenerator.h"
49
50ClassImp(AliPHOSRawFitterv2)
51
52//-----------------------------------------------------------------------------
53AliPHOSRawFitterv2::AliPHOSRawFitterv2():
54 AliPHOSRawFitterv0(),
55 fNtimeSamples(25),
56 fRMScut(11.)
57{
58 //Default constructor.
59 fLGpar[0]=0.971 ;
60 fLGpar[1]=0.0465;
61 fLGpar[2]=1.56 ;
62 fHGpar[0]=0.941 ;
63 fHGpar[1]=0.0436;
64 fHGpar[2]=1.65 ;
65}
66
67//-----------------------------------------------------------------------------
68AliPHOSRawFitterv2::~AliPHOSRawFitterv2()
69{
70 //Destructor.
71}
72
73//-----------------------------------------------------------------------------
74AliPHOSRawFitterv2::AliPHOSRawFitterv2(const AliPHOSRawFitterv2 &phosFitter ):
75 AliPHOSRawFitterv0(phosFitter),
76 fNtimeSamples(25),
77 fRMScut(11.)
78{
79 //Copy constructor.
80 fNtimeSamples=phosFitter.fNtimeSamples ;
81 for(Int_t i=0; i<3;i++){
82 fLGpar[i]=phosFitter.fLGpar[i] ;
83 fHGpar[i]=phosFitter.fHGpar[i] ;
84 }
85 fRMScut=phosFitter.fRMScut ;
86}
87
88//-----------------------------------------------------------------------------
89AliPHOSRawFitterv2& AliPHOSRawFitterv2::operator = (const AliPHOSRawFitterv2 &phosFitter)
90{
91 //Assignment operator.
92
93 fNtimeSamples=phosFitter.fNtimeSamples ;
94 for(Int_t i=0; i<3;i++){
95 fLGpar[i]=phosFitter.fLGpar[i] ;
96 fHGpar[i]=phosFitter.fHGpar[i] ;
97 }
98 fRMScut=phosFitter.fRMScut ;
99 return *this;
100}
101
102//-----------------------------------------------------------------------------
103Bool_t AliPHOSRawFitterv2::Eval()
104{
105 //Extract an energy deposited in the crystal,
106 //crystal' position (module,column,row),
107 //time and gain (high or low).
108 //First collects sample, then evaluates it and if it has
109 //reasonable shape, fits it with Gamma2 function and extracts
110 //energy and time.
111
112 if (fNBunches > 1) {
113 fQuality = 1000;
114 return kTRUE;
115 }
116
117 const Float_t kBaseLine = 1.0;
118 const Int_t kPreSamples = 10;
119
120 Float_t pedMean = 0;
121 Float_t pedRMS = 0;
122 Int_t nPed = 0;
123 UShort_t maxSample = 0;
124 Int_t nMax = 0;
125
126 TCanvas * cs = (TCanvas*)gROOT->FindObjectAny("CSample") ;
127 if(!cs)
128 cs = new TCanvas("CSample","CSample") ;
129
130 TH1D * h = (TH1D*)gROOT->FindObjectAny("hSample") ;
131 if(!h) h = new TH1D("hSample","",200,0.,200.) ;
132
133 Double_t pedestal = 0;
134 for (Int_t i=0; i<fLength; i++) {
135 if (i<kPreSamples) {
136 nPed++;
137 pedMean += fSignal[i];
138 pedRMS += fSignal[i]*fSignal[i] ;
139 }
140 if(fSignal[i] > maxSample) maxSample = fSignal[i];
141 if(fSignal[i] == maxSample) nMax++;
142 h->SetBinContent(i+1,fSignal[i]) ;
143 }
144 fEnergy = (Double_t)maxSample;
145 if (maxSample > 900 && nMax > 2) fOverflow = kTRUE;
146
147 if (fPedSubtract) {
148 if (nPed > 0) {
149 fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
150 if(fPedestalRMS > 0.)
151 fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
152 fEnergy -= (Double_t)(pedMean/nPed); // pedestal subtraction
153 }
154 else
155 return kFALSE;
156 }
157 else {
158 //take pedestals from DB
159 pedestal = (Double_t) fAmpOffset ;
160 if (fCalibData) {
161 Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
162 Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
163 pedestal += truePed - altroSettings ;
164 }
165 else{
166 AliWarning(Form("Can not read data from OCDB")) ;
167 }
168 fEnergy-=pedestal ;
169 }
170
171 if (fEnergy < kBaseLine) {
172 fEnergy = 0;
173 return kTRUE;
174 }
175
176 // calculate time
177 fTime=0. ;
178 Double_t tRMS = 0. ;
179 Double_t tW = 0. ;
180 Int_t cnts=0 ;
181 Double_t a=0,b=0,c=0 ;
182 if(fCaloFlag == 0){ // Low gain
183 a=fLGpar[0] ;
184 b=fLGpar[1] ;
185 c=fLGpar[2] ;
186 }
187 else if(fCaloFlag == 1){ // High gain
188 a=fHGpar[0] ;
189 b=fHGpar[1] ;
190 c=fHGpar[2] ;
191 }
192
193
194 fQuality = 0. ;
195
196 for(Int_t i=1; i<fLength && cnts<fNtimeSamples; i++){
197 if(fSignal[i] < pedestal)
198 continue ;
199 Double_t de = fSignal[i] - pedestal ;
200 Double_t av = fSignal[i-1] - pedestal + de;
201 if(av<=0.) //this is fluctuation around pedestal, skip it
202 continue ;
203 Double_t ds = fSignal[i] - fSignal[i-1] ;
204 Double_t ti = ds/av ; // calculate log. derivative
205 ti = a/(ti+b)-c*ti ; // and compare with parameterization
206 ti = i - ti ;
207 Double_t wi = TMath::Abs(ds) ;
208 fTime += ti*wi ;
209 tW += wi;
210 tRMS += ti*ti*wi ;
211 cnts++ ;
212 }
213
214 if(tW>0.){
215 fTime/=tW ;
216 fQuality = tRMS/tW-fTime*fTime ;
217 }
218 else{
219 fTime=-999. ;
220 fQuality=999. ;
221 }
222
223 Bool_t isBad = 0 ;
224 for(Int_t i=1; i<fLength-1&&!isBad; i++){
225 if(fSignal[i] > fSignal[i-1]+5 && fSignal[i] > fSignal[i+1]+5) { //single jump
226 isBad=1 ;
227 }
228 }
229 if(pedestal < 10.)
230 isBad=1 ;
231
232 if(fPedestalRMS > 0.1)
233 isBad=1 ;
234
235 for(Int_t i=1; i<fLength-1&&!isBad; i++){
236 if(fSignal[i] < pedestal-1)
237 isBad=1 ;
238 }
239
240 if(fEnergy>10. && !isBad ){
241 printf("fE=%f, ped=%f, fQuality=%f, pedRMS=%f \n",fEnergy,pedestal,fQuality,pedRMS) ;
242 if(fOverflow)printf(" Overflow \n") ;
243 if(isBad)printf("bad") ;
244 cs->cd() ;
245 h->Draw() ;
246 cs->Update() ;
247 getchar() ;
248 }
249
250 return kTRUE;
251}