Undid setting LHC polarity convention in the constructor
[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
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// 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
141 fEnergy = 0;
142 if (fNBunches > 1) {
143 fQuality = 1000;
144 return kTRUE;
145 }
146
147 const Float_t kBaseLine = 1.0;
148 const Int_t kPreSamples = 10;
149
150 Float_t pedMean = 0;
151 Float_t pedRMS = 0;
152 Int_t nPed = 0;
153 UShort_t maxSample = 0;
154 Int_t nMax = 0;
155
92236b27 156 for (Int_t i=0; i<sigLength; i++) {
fc5cc4cd 157 if (i>sigLength-kPreSamples) { //inverse signal time order
379c5c09 158 nPed++;
92236b27 159 pedMean += signal[i];
160 pedRMS += signal[i]*signal[i] ;
379c5c09 161 }
facd0d71 162 if(signal[i] > maxSample){ maxSample = signal[i]; nMax=0;}
92236b27 163 if(signal[i] == maxSample) nMax++;
a85e6165 164
379c5c09 165 }
fc5cc4cd 166
a85e6165 167
379c5c09 168 fEnergy = (Double_t)maxSample;
169 if (maxSample > 900 && nMax > 2) fOverflow = kTRUE;
170
fc5cc4cd 171 Double_t pedestal = 0 ;
379c5c09 172 if (fPedSubtract) {
173 if (nPed > 0) {
174 fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
175 if(fPedestalRMS > 0.)
176 fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
fc5cc4cd 177 pedestal = (Double_t)(pedMean/nPed);
379c5c09 178 }
179 else
180 return kFALSE;
181 }
182 else {
183 //take pedestals from DB
fc5cc4cd 184 pedestal = (Double_t) fAmpOffset ;
379c5c09 185 if (fCalibData) {
186 Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
187 Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
188 pedestal += truePed - altroSettings ;
189 }
190 else{
e10f780e 191 AliDebug(2,Form("Pedestal and offset are not read from OCDB. Use 0 for their values.")) ;
379c5c09 192 }
379c5c09 193 }
fc5cc4cd 194 fEnergy-=pedestal ;
379c5c09 195 if (fEnergy < kBaseLine) fEnergy = 0;
fc5cc4cd 196
197 //Evaluate time
facd0d71 198 fTime = sigStart;
fc5cc4cd 199 const Int_t nLine= 6 ; //Parameters of fitting
200 const Float_t eMinTOF = 10. ; //Choosed from beam-test and cosmic analyis
201 const Float_t kAmp=0.35 ; //Result slightly depends on them, so no getters
202 // Avoid too low peak:
203 if(fEnergy < eMinTOF){
204 return kTRUE;
205 }
206
207 // Find index posK (kLevel is a level of "timestamp" point Tk):
208 Int_t posK =sigLength-1 ; //last point before crossing k-level
209 Double_t levelK = pedestal + kAmp*fEnergy;
210 while(signal[posK] <= levelK && posK>=0){
211 posK-- ;
212 }
213 posK++ ;
214
215 if(posK == 0 || posK==sigLength-1){
216 return kTRUE;
217 }
218
219 // Find crosing point by solving linear equation (least squares method)
220 Int_t np = 0;
221 Int_t iup=posK-1 ;
222 Int_t idn=posK ;
223 Double_t sx = 0., sy = 0., sxx = 0., sxy = 0.;
224 Double_t x,y ;
225
226 while(np<nLine){
227 //point above crossing point
228 if(iup>=0){
229 x = sigLength-iup-1;
230 y = signal[iup];
231 sx += x;
232 sy += y;
233 sxx += (x*x);
234 sxy += (x*y);
235 np++ ;
236 iup-- ;
237 }
238 //Point below crossing point
239 if(idn<sigLength){
240 if(signal[idn]<pedestal){
241 idn=sigLength-1 ; //do not scan further
258abb9c 242 idn++ ;
fc5cc4cd 243 continue ;
244 }
245 x = sigLength-idn-1;
246 y = signal[idn];
247 sx += x;
248 sy += y;
249 sxx += (x*x);
250 sxy += (x*y);
251 np++;
252 idn++ ;
253 }
254 if(idn>=sigLength && iup<0){
255 break ; //can not fit futher
256 }
257 }
258
259 Double_t det = np*sxx - sx*sx;
260 if(det == 0){
261 return kTRUE;
262 }
263 Double_t c1 = (np*sxy - sx*sy)/det; //slope
264 Double_t c0 = (sy-c1*sx)/np; //offset
265 if(c1 == 0){
266 return kTRUE;
267 }
268
269 // Find where the line cross kLevel:
270 fTime += (levelK - c0)/c1-5. ; //5: mean offset between k-Level and start times
379c5c09 271 return kTRUE;
272
273}