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 : // 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();
25 : // fitterv0->SetChannelGeo(module,cellX,cellZ,caloFlag);
26 : // fitterv0->SetCalibData(fgCalibData) ;
27 : // fitterv0->Eval(sig,sigStart,sigLength);
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 :
44 22 : ClassImp(AliPHOSRawFitterv0)
45 :
46 : //-----------------------------------------------------------------------------
47 : AliPHOSRawFitterv0::AliPHOSRawFitterv0():
48 4 : TObject(),
49 4 : fModule(0),
50 4 : fCellX(0),
51 4 : fCellZ(0),
52 4 : fCaloFlag(0),
53 4 : fNBunches(0),
54 4 : fPedSubtract(kFALSE),
55 4 : fEnergy(-111),
56 4 : fTime(-111),
57 4 : fQuality(0.),
58 4 : fPedestalRMS(0.),
59 4 : fAmpOffset(0),
60 4 : fAmpThreshold(0),
61 4 : fOverflow(kFALSE),
62 4 : fCalibData(0)
63 20 : {
64 : //Default constructor
65 8 : }
66 :
67 : //-----------------------------------------------------------------------------
68 : AliPHOSRawFitterv0::~AliPHOSRawFitterv0()
69 16 : {
70 : //Destructor
71 16 : }
72 :
73 : //-----------------------------------------------------------------------------
74 : AliPHOSRawFitterv0::AliPHOSRawFitterv0(const AliPHOSRawFitterv0 &phosFitter ):
75 0 : TObject(),
76 0 : fModule (phosFitter.fModule),
77 0 : fCellX (phosFitter.fCellX),
78 0 : fCellZ (phosFitter.fCellZ),
79 0 : fCaloFlag (phosFitter.fCaloFlag),
80 0 : fNBunches (phosFitter.fNBunches),
81 0 : fPedSubtract (phosFitter.fPedSubtract),
82 0 : fEnergy (phosFitter.fEnergy),
83 0 : fTime (phosFitter.fTime),
84 0 : fQuality (phosFitter.fQuality),
85 0 : fPedestalRMS (phosFitter.fPedestalRMS),
86 0 : fAmpOffset (phosFitter.fAmpOffset),
87 0 : fAmpThreshold(phosFitter.fAmpThreshold),
88 0 : fOverflow (phosFitter.fOverflow),
89 0 : fCalibData (phosFitter.fCalibData)
90 0 : {
91 : //Copy constructor
92 0 : }
93 :
94 : //-----------------------------------------------------------------------------
95 : AliPHOSRawFitterv0& AliPHOSRawFitterv0::operator = (const AliPHOSRawFitterv0 &phosFitter)
96 : {
97 : //Assignment operator.
98 :
99 0 : if(this != &phosFitter) {
100 0 : fModule = phosFitter.fModule;
101 0 : fCellX = phosFitter.fCellX;
102 0 : fCellZ = phosFitter.fCellZ;
103 0 : fCaloFlag = phosFitter.fCaloFlag;
104 0 : fNBunches = phosFitter.fNBunches;
105 0 : fPedSubtract = phosFitter.fPedSubtract;
106 0 : fEnergy = phosFitter.fEnergy;
107 0 : fTime = phosFitter.fTime;
108 0 : fQuality = phosFitter.fQuality;
109 0 : fPedestalRMS = phosFitter.fPedestalRMS;
110 0 : fAmpOffset = phosFitter.fAmpOffset;
111 0 : fAmpThreshold = phosFitter.fAmpThreshold;
112 0 : fOverflow = phosFitter.fOverflow;
113 0 : fCalibData = phosFitter.fCalibData;
114 0 : }
115 :
116 0 : return *this;
117 : }
118 :
119 : //-----------------------------------------------------------------------------
120 :
121 : void 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 308 : fModule = module;
128 154 : fCellX = cellX;
129 154 : fCellZ = cellZ;
130 154 : fCaloFlag = caloFlag;
131 154 : }
132 : //-----------------------------------------------------------------------------
133 :
134 : Bool_t AliPHOSRawFitterv0::Eval(const UShort_t *signal, Int_t sigStart, Int_t sigLength)
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 308 : fOverflow= kFALSE ;
142 154 : fEnergy = 0;
143 154 : if (fNBunches > 1) {
144 0 : fQuality = 1000;
145 0 : 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 :
157 17338 : for (Int_t i=0; i<sigLength; i++) {
158 8515 : if (i>sigLength-kPreSamples) { //inverse signal time order
159 1258 : nPed++;
160 1258 : pedMean += signal[i];
161 1258 : pedRMS += signal[i]*signal[i] ;
162 1258 : }
163 10667 : if(signal[i] > maxSample){ maxSample = signal[i]; nMax=0;}
164 15485 : if(signal[i] == maxSample) nMax++;
165 :
166 : }
167 :
168 :
169 154 : fEnergy = (Double_t)maxSample;
170 156 : if (maxSample > 900 && nMax > 2) fOverflow = kTRUE;
171 :
172 : Double_t pedestal = 0 ;
173 154 : if (fPedSubtract) {
174 0 : if (nPed > 0) {
175 0 : fPedestalRMS=(pedRMS - pedMean*pedMean/nPed)/nPed ;
176 0 : if(fPedestalRMS > 0.)
177 0 : fPedestalRMS = TMath::Sqrt(fPedestalRMS) ;
178 0 : pedestal = (Double_t)(pedMean/nPed);
179 : }
180 : else
181 0 : return kFALSE;
182 0 : }
183 : else {
184 : //take pedestals from DB
185 154 : pedestal = (Double_t) fAmpOffset ;
186 154 : if (fCalibData) {
187 154 : Float_t truePed = fCalibData->GetADCpedestalEmc(fModule, fCellZ, fCellX) ;
188 154 : Int_t altroSettings = fCalibData->GetAltroOffsetEmc(fModule, fCellZ, fCellX) ;
189 154 : pedestal += truePed - altroSettings ;
190 154 : }
191 : else{
192 0 : AliDebug(2,Form("Pedestal and offset are not read from OCDB. Use 0 for their values.")) ;
193 : }
194 : }
195 154 : fEnergy-=pedestal ;
196 170 : if (fEnergy < kBaseLine) fEnergy = 0;
197 :
198 : //Evaluate time
199 154 : fTime = sigStart-sigLength-3;
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 154 : if(fEnergy < eMinTOF){
205 101 : return kTRUE;
206 : }
207 :
208 : // Find index posK (kLevel is a level of "timestamp" point Tk):
209 53 : Int_t posK =sigLength-1 ; //last point before crossing k-level
210 53 : Double_t levelK = pedestal + kAmp*fEnergy;
211 612 : while(signal[posK] <= levelK && posK>=0){
212 253 : posK-- ;
213 : }
214 53 : posK++ ;
215 :
216 106 : if(posK == 0 || posK==sigLength-1){
217 0 : 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 265 : while(np<nLine){
228 : //point above crossing point
229 159 : if(iup>=0){
230 159 : x = sigLength-iup-1;
231 159 : y = signal[iup];
232 159 : sx += x;
233 159 : sy += y;
234 159 : sxx += (x*x);
235 159 : sxy += (x*y);
236 159 : np++ ;
237 159 : iup-- ;
238 159 : }
239 : //Point below crossing point
240 159 : if(idn<sigLength){
241 159 : if(signal[idn]<pedestal){
242 : idn=sigLength-1 ; //do not scan further
243 : idn++ ;
244 0 : continue ;
245 : }
246 159 : x = sigLength-idn-1;
247 159 : y = signal[idn];
248 159 : sx += x;
249 159 : sy += y;
250 159 : sxx += (x*x);
251 159 : sxy += (x*y);
252 159 : np++;
253 159 : idn++ ;
254 159 : }
255 318 : if(idn>=sigLength && iup<0){
256 : break ; //can not fit futher
257 : }
258 : }
259 :
260 53 : Double_t det = np*sxx - sx*sx;
261 53 : if(det == 0){
262 0 : return kTRUE;
263 : }
264 53 : if(np == 0){
265 0 : return kFALSE;
266 : }
267 53 : Double_t c1 = (np*sxy - sx*sy)/det; //slope
268 53 : Double_t c0 = (sy-c1*sx)/np; //offset
269 53 : if(c1 == 0){
270 0 : return kTRUE;
271 : }
272 :
273 : // Find where the line cross kLevel:
274 53 : fTime += (levelK - c0)/c1-5. ; //5: mean offset between k-Level and start times
275 53 : return kTRUE;
276 :
277 154 : }
|
