]>
Commit | Line | Data |
---|---|---|
48a2e3eb | 1 | /************************************************************************** |
2 | * This file is property of and copyright by * | |
3 | * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 * | |
4 | * * | |
5 | * Primary Author: Per Thomas Hille <perthomas.hille@yale.edu> * | |
6 | * * | |
7 | * Contributors are mentioned in the code where appropriate. * | |
8 | * Please report bugs to p.t.hille@fys.uio.no * | |
9 | * * | |
10 | * Permission to use, copy, modify and distribute this software and its * | |
11 | * documentation strictly for non-commercial purposes is hereby granted * | |
12 | * without fee, provided that the above copyright notice appears in all * | |
13 | * copies and that both the copyright notice and this permission notice * | |
14 | * appear in the supporting documentation. The authors make no claims * | |
15 | * about the suitability of this software for any purpose. It is * | |
16 | * provided "as is" without express or implied warranty. * | |
17 | **************************************************************************/ | |
18 | ||
19 | // The Peak-Finder algorithm | |
20 | // The amplitude is extracted as a | |
21 | // weighted sum of the samples using the | |
22 | // best possible weights. | |
23 | // The wights is calculated only once and the | |
24 | // Actual extraction of amplitude and peak position | |
25 | // Is done with a simple vector multiplication, allowing for | |
26 | // Extreemely fast computations. | |
27 | ||
28 | #include "AliCaloRawAnalyzerPeakFinderV2.h" | |
29 | #include "AliCaloBunchInfo.h" | |
30 | #include "AliCaloFitResults.h" | |
31 | #include <iostream> | |
32 | #include "unistd.h" | |
33 | #include "TMath.h" | |
34 | #include "AliLog.h" | |
35 | ||
36 | using namespace std; | |
37 | ||
38 | ClassImp( AliCaloRawAnalyzerPeakFinderV2 ) | |
39 | ||
40 | AliCaloRawAnalyzerPeakFinderV2::AliCaloRawAnalyzerPeakFinderV2() :AliCaloRawAnalyzer("Peak-FinderV2", "PF") | |
41 | // fTof(0), | |
42 | // fAmp(0) | |
43 | { | |
44 | //comment | |
45 | ||
46 | fNsampleCut = 5; | |
47 | ||
48 | for(int i=0; i < MAXSTART; i++) | |
49 | { | |
50 | for(int j=0; j < SAMPLERANGE; j++ ) | |
51 | { | |
52 | fPFAmpVectors[i][j] = new double[100]; | |
53 | fPFTofVectors[i][j] = new double[100]; | |
54 | fPFAmpVectorsCoarse[i][j] = new double[100]; | |
55 | fPFTofVectorsCoarse[i][j] = new double[100]; | |
56 | ||
57 | ||
58 | for(int k=0; k < 100; k++ ) | |
59 | { | |
60 | fPFAmpVectors[i][j][k] = 0; | |
61 | fPFTofVectors[i][j][k] = 0; | |
62 | fPFAmpVectorsCoarse[i][j][k] = 0; | |
63 | fPFTofVectorsCoarse[i][j][k] = 0; | |
64 | } | |
65 | } | |
66 | } | |
67 | ||
68 | LoadVectors(); | |
69 | ||
70 | } | |
71 | ||
72 | ||
73 | AliCaloRawAnalyzerPeakFinderV2::~AliCaloRawAnalyzerPeakFinderV2() | |
74 | { | |
75 | //comment | |
76 | for(int i=0; i < MAXSTART; i++) | |
77 | { | |
78 | for(int j=0; j < SAMPLERANGE; j++ ) | |
79 | { | |
80 | delete[] fPFAmpVectors[i][j]; | |
81 | delete[] fPFTofVectors[i][j]; | |
82 | delete[] fPFAmpVectorsCoarse[i][j]; | |
83 | delete[] fPFTofVectorsCoarse[i][j]; | |
84 | } | |
85 | } | |
86 | } | |
87 | ||
88 | ||
89 | Double_t | |
90 | AliCaloRawAnalyzerPeakFinderV2::ScanCoarse(const Double_t *const array, const int length ) const | |
91 | { | |
92 | Double_t tmpTof = 0; | |
93 | Double_t tmpAmp= 0; | |
94 | ||
95 | for(int i=0; i < length; i++) | |
96 | { | |
97 | tmpTof += fPFTofVectorsCoarse[0][length][i]*array[i]; | |
98 | tmpAmp += fPFAmpVectorsCoarse[0][length][i]*array[i]; | |
99 | } | |
100 | ||
101 | tmpTof = tmpTof / tmpAmp ; | |
102 | return tmpTof; | |
103 | } | |
104 | ||
105 | ||
106 | AliCaloFitResults | |
107 | AliCaloRawAnalyzerPeakFinderV2::Evaluate( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2 ) | |
108 | { | |
109 | // Extracting the amplitude using the Peak-FinderV2 algorithm | |
110 | // The amplitude is a weighted sum of the samples using | |
111 | // optimum weights. | |
112 | ||
113 | short maxampindex; //index of maximum amplitude | |
114 | short maxamp; //Maximum amplitude | |
115 | // fAmp = 0; | |
116 | fAmpA[0] = 0; | |
117 | fAmpA[1] = 0; | |
118 | fAmpA[2] = 0; | |
119 | ||
120 | int index = SelectBunch( bunchvector, &maxampindex, &maxamp ); | |
121 | ||
122 | if( index >= 0) | |
123 | { | |
124 | Float_t ped = ReverseAndSubtractPed( &(bunchvector.at(index)) , altrocfg1, altrocfg2, fReversed ); | |
125 | Float_t maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed ); | |
126 | ||
127 | if( maxf < fAmpCut || ( maxamp - ped) > 900 ) // (maxamp - ped) > 900 = Close to saturation (use low gain then) | |
128 | { | |
129 | // cout << __FILE__ << __LINE__ <<":, maxamp = " << maxamp << ", ped = "<< ped << ",. maxf = "<< maxf << ", maxampindex = "<< maxampindex << endl; | |
130 | return AliCaloFitResults( maxamp, ped, -1, maxf, maxampindex, -1, -1 ); | |
131 | } | |
132 | ||
133 | int first; | |
134 | int last; | |
135 | ||
136 | if ( maxf > fAmpCut ) | |
137 | { | |
138 | SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxampindex - bunchvector.at(index).GetStartBin(), &first, &last); | |
139 | int nsamples = last - first; | |
140 | if( ( nsamples ) >= fNsampleCut ) | |
141 | { | |
142 | int startbin = bunchvector.at(index).GetStartBin(); | |
143 | int n = last - first; | |
144 | int pfindex = n - fNsampleCut; | |
145 | pfindex = pfindex > SAMPLERANGE ? SAMPLERANGE : pfindex; | |
146 | ||
147 | int dt = maxampindex - startbin -1; | |
148 | ||
149 | ||
150 | // cout << __FILE__ << __LINE__ <<"\t The coarse estimated t0 is " << ScanCoarse( &fReversed[dt] , n ) << endl; | |
151 | ||
152 | ||
153 | ||
154 | // Float_t tmptof = ScanCoarse( &fReversed[dt] , n ); | |
155 | // cout << __FILE__ << __LINE__ << "\ttmptof = " << tmptof << endl; | |
156 | ||
157 | ||
158 | // cout << __FILE__ << __LINE__ << ", dt = " << dt << ",\tmaxamindex = " << maxampindex << "\tstartbin = "<< startbin << endl; | |
159 | ||
160 | for( int i=0; i < SAMPLERANGE; i++ ) | |
161 | { | |
162 | for( int j = 0; j < 3; j++ ) | |
163 | { | |
164 | // fAmpA[j] += fPFAmpVectors[0][pfindex][i]*tmp[j]; | |
165 | fAmpA[j] += fPFAmpVectors[0][pfindex][i]*fReversed[ dt +i +j -1 ]; | |
166 | } | |
167 | } | |
168 | ||
169 | double diff = 9999; | |
170 | int tmpindex = 0; | |
171 | ||
172 | for(int k=0; k < 3; k ++) | |
173 | { | |
174 | // cout << __FILE__ << __LINE__ << "amp[="<< k <<"] = " << fAmpA[k] << endl; | |
175 | if( TMath::Abs(fAmpA[k] - ( maxamp - ped) ) < diff) | |
176 | { | |
177 | diff = TMath::Abs(fAmpA[k] - ( maxamp - ped)); | |
178 | tmpindex = k; | |
179 | } | |
180 | } | |
181 | ||
182 | Float_t tmptof = ScanCoarse( &fReversed[dt] , n ); | |
183 | ||
184 | // Double_t tofnew = PolTof(tmptof) + ( dt + startbin )*100 ; | |
185 | // int dt = maxampindex - startbin -1; | |
186 | ||
187 | // Double_t tofnew = PolTof(tmptof) + startbin*100 ; | |
188 | ||
189 | // Double_t tofnew = PolTof(tmptof) + maxampindex*100 ; | |
190 | ||
191 | ||
192 | // Double_t tofnew = PolTof(tmptof) + (dt + startbin + tmpindex )*100 ; | |
193 | Double_t tofnew = PolTof(tmptof) ; | |
194 | ||
195 | // tmptof= tofnew; | |
196 | ||
197 | // Double_t tofnew = PolTof(tmptof) + maxampindex ; | |
198 | ||
199 | ||
200 | if(tmptof < 0 ) | |
201 | { | |
202 | cout << __FILE__ << __LINE__ << "\ttmptof = " << tmptof << endl; | |
203 | } | |
204 | ||
205 | ||
206 | if( tmptof < -1 ) | |
207 | { | |
208 | tmpindex = 0; | |
209 | } | |
210 | else | |
211 | if( tmptof > -1 && tmptof < 100 ) | |
212 | { | |
213 | tmpindex =1; | |
214 | } | |
215 | else | |
216 | { | |
217 | tmpindex = 2; | |
218 | } | |
219 | ||
220 | ||
221 | ||
222 | double tof = 0; | |
223 | ||
224 | for(int k=0; k < SAMPLERANGE; k++ ) | |
225 | { | |
226 | tof += fPFTofVectors[0][pfindex][k]*fReversed[ dt +k + tmpindex -1 ]; | |
227 | } | |
228 | ||
229 | // cout << __FILE__ << __LINE__ << "tofRaw = "<< tof / fAmpA[tmpindex] << endl; | |
230 | ||
231 | // tof = tof / fAmpA[tmpindex] + (dt + startbin)*100; | |
232 | ||
233 | if( TMath::Abs( (maxf - fAmpA[tmpindex])/maxf ) > 0.1 ) | |
234 | { | |
235 | fAmpA[tmpindex] = maxf; | |
236 | } | |
237 | ||
238 | // tof = (dt + startbin + tmpindex )*100 - tof/fAmpA[tmpindex]; // ns | |
239 | tof = dt + startbin + tmpindex - 0.01*tof/fAmpA[tmpindex]; // clock ticks | |
240 | ||
241 | ||
242 | return AliCaloFitResults( maxamp, ped , -1, fAmpA[tmpindex], tof , -2, -3 ); | |
243 | } | |
244 | else | |
245 | { | |
246 | return AliCaloFitResults( maxamp, ped , -5, maxf, -6, -7, -8 ); | |
247 | } | |
248 | } | |
249 | } | |
250 | // cout << __FILE__ << __LINE__ << "WARNING, returning amp = -1 " << endl; | |
251 | return AliCaloFitResults(-1, -1); | |
252 | } | |
253 | ||
254 | ||
255 | Double_t | |
256 | AliCaloRawAnalyzerPeakFinderV2::PolTof( const double fx1 ) const | |
257 | { | |
258 | ||
259 | //Newtons method | |
260 | Double_t tolerance = 0.01; | |
261 | // cout << "************************" << endl; | |
262 | Double_t fx2 = PolValue( fx1 ); | |
263 | Double_t tmpfx1 = fx1; | |
264 | ||
265 | while( TMath::Abs( fx2 - fx1 ) > tolerance ) | |
266 | { | |
267 | Double_t der = PolDerivative( tmpfx1 ); | |
268 | // tmpx = der*( x - val) +x; | |
269 | tmpfx1 = ( fx1 - fx2)/der +tmpfx1; | |
270 | ||
271 | // tmpx = der*( val - tmpx ) +tmpx; | |
272 | fx2 = PolValue( tmpfx1 ); | |
273 | // cout << __FILE__ << __LINE__ << "Der =\t" << der << " x=\t"<<x<<"\tval="<<val << endl; | |
274 | // cout << __FILE__ << __LINE__ << "Der =\t" << der << " tmpx=\t"<< tmpfx1 <<"\tval="<< fx2 << endl; | |
275 | } | |
276 | ||
277 | // cout << __FILE__ << __LINE__ << "CONVERGED !! fx1 = "<< fx1 <<" tmpfx1 = "<< tmpfx1 <<" f(tmpfx1) = "<< fx2 << endl; | |
278 | ||
279 | // cout << "************************" << endl; | |
280 | ||
281 | return tmpfx1; | |
282 | ||
283 | } | |
284 | ||
285 | ||
286 | Double_t | |
287 | AliCaloRawAnalyzerPeakFinderV2::PolValue(const Double_t x) const | |
288 | { | |
289 | static Double_t p0 = -55.69; | |
290 | static Double_t p1 = 4.718; | |
291 | static Double_t p2 = -0.05587; | |
292 | static Double_t p3 = 0.0003185; | |
293 | static Double_t p4 = -7.91E-7; | |
294 | static Double_t p5 = 7.576E-10; | |
295 | ||
296 | // return p0 + p1*x + p2*TMath::Power(x, 2) + p3*TMath::Power(x, 3) + p4*TMath::Power(x, 4) + p5*TMath::Power(x, 5); | |
297 | ||
298 | return p0 + p1*x + p2*x*x + p3*x*x*x + p4*x*x*x*x + p5*x*x*x*x*x; | |
299 | ||
300 | } | |
301 | ||
302 | ||
303 | Double_t | |
304 | AliCaloRawAnalyzerPeakFinderV2::PolDerivative(const Double_t x) const | |
305 | { | |
306 | static Double_t dp0 = 0; | |
307 | static Double_t dp1 = 4.718; | |
308 | static Double_t dp2 = -0.11174; | |
309 | static Double_t dp3 = 0.0009555; | |
310 | static Double_t dp4 = -3.164E-6; | |
311 | static Double_t dp5 = 3.788E-9; | |
312 | ||
313 | // return dp0 + dp1 + dp2*x + dp3*TMath::Power(x, 2) + dp4*TMath::Power(x, 3) + dp5*TMath::Power(x, 4); | |
314 | ||
315 | ||
316 | ||
317 | return dp0 + dp1 + dp2*x + dp3*x*x + dp4*x*x*x + dp5*x*x*x*x; | |
318 | ||
319 | } | |
320 | ||
321 | ||
322 | void | |
323 | AliCaloRawAnalyzerPeakFinderV2::LoadVectors() | |
324 | { | |
325 | //Read in the Peak finder vecors from file | |
326 | for(int i = 0; i < MAXSTART ; i++) | |
327 | { | |
328 | for( int j=0; j < SAMPLERANGE; j++) | |
329 | { | |
330 | char filenameCoarse[256]; | |
331 | char filename[256]; | |
332 | ||
333 | int n = j+fNsampleCut; | |
334 | ||
335 | // double start = (double)i+0.5; | |
336 | double start = (double)i+0; | |
337 | ||
338 | sprintf(filename, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt1.0.txt", getenv("ALICE_ROOT"), start, n); | |
339 | sprintf(filenameCoarse, "%s/EMCAL/vectors-emcal/start%.1fN%dtau0.235fs10dt3.0.txt", getenv("ALICE_ROOT"), start, n); | |
340 | ||
341 | FILE *fp = fopen(filename, "r"); | |
342 | FILE *fpc = fopen(filenameCoarse, "r"); | |
343 | ||
344 | if( fp == 0 ) | |
345 | { | |
346 | AliFatal( Form( "could not open file: %s", filename ) ); | |
347 | } | |
348 | ||
349 | if(fpc == 0) | |
350 | { | |
351 | AliFatal( Form( "could not open file: %s", filenameCoarse ) ); | |
352 | } | |
353 | else | |
354 | { | |
355 | for(int m = 0; m < n ; m++ ) | |
356 | { | |
357 | cout << __FILE__ << __LINE__ << "i="<<i <<"\tj=" <<j << "\tm=" << m << endl; | |
358 | ||
359 | fscanf(fp, "%lf\t", &fPFAmpVectors[i][j][m] ); | |
360 | // fPFAmpVectorsCoarse[i][j][m] = 1; | |
361 | fscanf(fpc, "%lf\t", &fPFAmpVectorsCoarse[i][j][m] ); | |
362 | } | |
363 | ||
364 | fscanf(fp, "\n" ); | |
365 | fscanf(fpc, "\n" ); | |
366 | ||
367 | for(int m = 0; m < n ; m++ ) | |
368 | { | |
369 | // fPFTofVectors[i][j][m] = 1; | |
370 | ||
371 | fscanf(fp, "%lf\t", &fPFTofVectors[i][j][m] ); | |
372 | fscanf(fpc, "%lf\t", &fPFTofVectorsCoarse[i][j][m] ); | |
373 | // fPFTofVectorsCoarse[i][j][m] = 1; | |
374 | } | |
375 | fclose (fp); | |
376 | fclose (fpc); | |
377 | } | |
378 | } | |
379 | } | |
380 | } | |
381 | ||
382 |