Updates for the ROOT trunk
[u/mrichter/AliRoot.git] / PHOS / AliPHOSRawFitterv0.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
f18d93e0 16/* $Id$ */
379c5c09 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// A coarse algorithm takes the energy as the maximum
21// sample, time is the first sample index, and the quality is the number
22// of bunches in the signal.
23//
24// AliPHOSRawFitterv0 *fitterv0=new AliPHOSRawFitterv0();
379c5c09 25// fitterv0->SetChannelGeo(module,cellX,cellZ,caloFlag);
26// fitterv0->SetCalibData(fgCalibData) ;
1dfadc32 27// fitterv0->Eval(sig,sigStart,sigLength);
379c5c09 28// Double_t amplitude = fitterv0.GetEnergy();
29// Double_t time = fitterv0.GetTime();
30// Bool_t isLowGain = fitterv0.GetCaloFlag()==0;
31
32// Author: Yuri Kharlov
33
34// --- ROOT system ---
35#include "TArrayI.h"
36#include "TMath.h"
37#include "TObject.h"
38
39// --- AliRoot header files ---
40#include "AliPHOSRawFitterv0.h"
41#include "AliPHOSCalibData.h"
42#include "AliLog.h"
43
44ClassImp(AliPHOSRawFitterv0)
45
46//-----------------------------------------------------------------------------
47AliPHOSRawFitterv0::AliPHOSRawFitterv0():
48 TObject(),
379c5c09 49 fModule(0),
50 fCellX(0),
51 fCellZ(0),
52 fCaloFlag(0),
379c5c09 53 fNBunches(0),
54 fPedSubtract(kFALSE),
55 fEnergy(-111),
56 fTime(-111),
57 fQuality(0.),
58 fPedestalRMS(0.),
59 fAmpOffset(0),
60 fAmpThreshold(0),
61 fOverflow(kFALSE),
62 fCalibData(0)
63{
64 //Default constructor
65}
66
67//-----------------------------------------------------------------------------
68AliPHOSRawFitterv0::~AliPHOSRawFitterv0()
69{
70 //Destructor
379c5c09 71}
72
73//-----------------------------------------------------------------------------
74AliPHOSRawFitterv0::AliPHOSRawFitterv0(const AliPHOSRawFitterv0 &phosFitter ):
75 TObject(),
379c5c09 76 fModule (phosFitter.fModule),
77 fCellX (phosFitter.fCellX),
78 fCellZ (phosFitter.fCellZ),
79 fCaloFlag (phosFitter.fCaloFlag),
379c5c09 80 fNBunches (phosFitter.fNBunches),
81 fPedSubtract (phosFitter.fPedSubtract),
82 fEnergy (phosFitter.fEnergy),
83 fTime (phosFitter.fTime),
84 fQuality (phosFitter.fQuality),
85 fPedestalRMS (phosFitter.fPedestalRMS),
86 fAmpOffset (phosFitter.fAmpOffset),
87 fAmpThreshold(phosFitter.fAmpThreshold),
88 fOverflow (phosFitter.fOverflow),
89 fCalibData (phosFitter.fCalibData)
90{
91 //Copy constructor
92}
93
94//-----------------------------------------------------------------------------
95AliPHOSRawFitterv0& AliPHOSRawFitterv0::operator = (const AliPHOSRawFitterv0 &phosFitter)
96{
97 //Assignment operator.
98
99 if(this != &phosFitter) {
379c5c09 100 fModule = phosFitter.fModule;
101 fCellX = phosFitter.fCellX;
102 fCellZ = phosFitter.fCellZ;
103 fCaloFlag = phosFitter.fCaloFlag;
379c5c09 104 fNBunches = phosFitter.fNBunches;
105 fPedSubtract = phosFitter.fPedSubtract;
106 fEnergy = phosFitter.fEnergy;
107 fTime = phosFitter.fTime;
108 fQuality = phosFitter.fQuality;
109 fPedestalRMS = phosFitter.fPedestalRMS;
110 fAmpOffset = phosFitter.fAmpOffset;
111 fAmpThreshold = phosFitter.fAmpThreshold;
112 fOverflow = phosFitter.fOverflow;
113 fCalibData = phosFitter.fCalibData;
114 }
115
116 return *this;
117}
118
119//-----------------------------------------------------------------------------
120
379c5c09 121void AliPHOSRawFitterv0::SetChannelGeo(const Int_t module, const Int_t cellX,
122 const Int_t cellZ, const Int_t caloFlag)
123{
124 // Set geometry address of the channel
125 // (for a case if fitting parameters are different for different channels)
126
127 fModule = module;
128 fCellX = cellX;
129 fCellZ = cellZ;
130 fCaloFlag = caloFlag;
131}
132//-----------------------------------------------------------------------------
133
fc5cc4cd 134Bool_t AliPHOSRawFitterv0::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
379c5c09 135{
136 // Calculate signal parameters (energy, time, quality) from array of samples
137 // Energy is a maximum sample minus pedestal 9
138 // Time is the first time bin
139 // Signal overflows is there are at least 3 samples of the same amplitude above 900
140
61aa7aac 141 fOverflow= kFALSE ;
379c5c09 142 fEnergy = 0;
143 if (fNBunches > 1) {
144 fQuality = 1000;
145 return kTRUE;
146 }
147
148 const Float_t kBaseLine = 1.0;
149 const Int_t kPreSamples = 10;
150
151 Float_t pedMean = 0;
152 Float_t pedRMS = 0;
153 Int_t nPed = 0;
154 UShort_t maxSample = 0;
155 Int_t nMax = 0;
156
92236b27 157 for (Int_t i=0; i<sigLength; i++) {
fc5cc4cd 158 if (i>sigLength-kPreSamples) { //inverse signal time order
379c5c09 159 nPed++;
92236b27 160 pedMean += signal[i];
161 pedRMS += signal[i]*signal[i] ;
379c5c09 162 }
facd0d71 163 if(signal[i] > maxSample){ maxSample = signal[i]; nMax=0;}
92236b27 164 if(signal[i] == maxSample) nMax++;
a85e6165 165
379c5c09 166 }
fc5cc4cd 167
a85e6165 168
379c5c09 169 fEnergy = (Double_t)maxSample;
170 if (maxSample > 900 && nMax > 2) fOverflow = kTRUE;
171
fc5cc4cd 172 Double_t pedestal = 0 ;
379c5c09 173 if (fPedSubtract) {
174 if (nPed > 0) {
175 fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
176 if(fPedestalRMS > 0.)
177 fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
fc5cc4cd 178 pedestal = (Double_t)(pedMean/nPed);
379c5c09 179 }
180 else
181 return kFALSE;
182 }
183 else {
184 //take pedestals from DB
fc5cc4cd 185 pedestal = (Double_t) fAmpOffset ;
379c5c09 186 if (fCalibData) {
187 Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
188 Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
189 pedestal += truePed - altroSettings ;
190 }
191 else{
e10f780e 192 AliDebug(2,Form("Pedestal and offset are not read from OCDB. Use 0 for their values.")) ;
379c5c09 193 }
379c5c09 194 }
fc5cc4cd 195 fEnergy-=pedestal ;
379c5c09 196 if (fEnergy < kBaseLine) fEnergy = 0;
fc5cc4cd 197
198 //Evaluate time
df28635d 199 fTime = sigStart-sigLength-3;
fc5cc4cd 200 const Int_t nLine= 6 ; //Parameters of fitting
201 const Float_t eMinTOF = 10. ; //Choosed from beam-test and cosmic analyis
202 const Float_t kAmp=0.35 ; //Result slightly depends on them, so no getters
203 // Avoid too low peak:
204 if(fEnergy < eMinTOF){
205 return kTRUE;
206 }
207
208 // Find index posK (kLevel is a level of "timestamp" point Tk):
209 Int_t posK =sigLength-1 ; //last point before crossing k-level
210 Double_t levelK = pedestal + kAmp*fEnergy;
211 while(signal[posK] <= levelK && posK>=0){
212 posK-- ;
213 }
214 posK++ ;
215
216 if(posK == 0 || posK==sigLength-1){
217 return kTRUE;
218 }
219
220 // Find crosing point by solving linear equation (least squares method)
221 Int_t np = 0;
222 Int_t iup=posK-1 ;
223 Int_t idn=posK ;
224 Double_t sx = 0., sy = 0., sxx = 0., sxy = 0.;
225 Double_t x,y ;
226
227 while(np<nLine){
228 //point above crossing point
229 if(iup>=0){
230 x = sigLength-iup-1;
231 y = signal[iup];
232 sx += x;
233 sy += y;
234 sxx += (x*x);
235 sxy += (x*y);
236 np++ ;
237 iup-- ;
238 }
239 //Point below crossing point
240 if(idn<sigLength){
241 if(signal[idn]<pedestal){
242 idn=sigLength-1 ; //do not scan further
258abb9c 243 idn++ ;
fc5cc4cd 244 continue ;
245 }
246 x = sigLength-idn-1;
247 y = signal[idn];
248 sx += x;
249 sy += y;
250 sxx += (x*x);
251 sxy += (x*y);
252 np++;
253 idn++ ;
254 }
255 if(idn>=sigLength && iup<0){
256 break ; //can not fit futher
257 }
258 }
259
260 Double_t det = np*sxx - sx*sx;
261 if(det == 0){
262 return kTRUE;
263 }
929af5d7 264 if(np == 0){
265 return kFALSE;
266 }
fc5cc4cd 267 Double_t c1 = (np*sxy - sx*sy)/det; //slope
268 Double_t c0 = (sy-c1*sx)/np; //offset
269 if(c1 == 0){
270 return kTRUE;
271 }
272
273 // Find where the line cross kLevel:
274 fTime += (levelK - c0)/c1-5. ; //5: mean offset between k-Level and start times
379c5c09 275 return kTRUE;
276
277}