]>
Commit | Line | Data |
---|---|---|
168c7b3c | 1 | // -*- mode: c++ -*- |
d655d7dd | 2 | /************************************************************************** |
3 | * This file is property of and copyright by * | |
4 | * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 * | |
5 | * * | |
e37e3c84 | 6 | * Primary Author: Per Thomas Hille <perthomas.hille@yale.edu> * |
d655d7dd | 7 | * * |
8 | * Contributors are mentioned in the code where appropriate. * | |
9 | * Please report bugs to p.t.hille@fys.uio.no * | |
10 | * * | |
11 | * Permission to use, copy, modify and distribute this software and its * | |
12 | * documentation strictly for non-commercial purposes is hereby granted * | |
13 | * without fee, provided that the above copyright notice appears in all * | |
14 | * copies and that both the copyright notice and this permission notice * | |
15 | * appear in the supporting documentation. The authors make no claims * | |
16 | * about the suitability of this software for any purpose. It is * | |
17 | * provided "as is" without express or implied warranty. * | |
18 | **************************************************************************/ | |
19 | ||
20 | ||
21 | // Base class for extraction | |
22 | // of signal amplitude and peak position | |
57839add | 23 | // From CALO Calorimeter RAW data (from the RCU) |
d655d7dd | 24 | // Contains some utilities for preparing / selecting |
25 | // Signals suitable for signal extraction | |
26 | // By derived classes | |
27 | ||
28 | #include "AliLog.h" | |
57839add | 29 | #include "AliCaloRawAnalyzer.h" |
30 | #include "AliCaloBunchInfo.h" | |
31 | #include "AliCaloFitResults.h" | |
32 | #include "TMath.h" | |
d655d7dd | 33 | #include <iostream> |
34 | using namespace std; | |
35 | ||
ce95bae9 | 36 | ClassImp(AliCaloRawAnalyzer) |
48a2e3eb | 37 | |
38 | AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char *name, const char *nameshort) : TObject(), | |
168c7b3c | 39 | fMinTimeIndex(-1), |
40 | fMaxTimeIndex(-1), | |
41 | fFitArrayCut(5), | |
42 | fAmpCut(4), | |
43 | fNsampleCut(5), | |
44 | fOverflowCut(950), | |
45 | fNsamplePed(3), | |
46 | fIsZerosupressed( false ), | |
47 | fVerbose( false ), | |
92d9f317 | 48 | fAlgo(Algo::kNONE), |
49 | // fFp(0), | |
50 | fL1Phase(0), | |
51 | fAmp(0), | |
396baaf6 | 52 | fTof(0), |
53 | fTau( EMCAL::TAU ) | |
54 | // fFixTau(true) | |
d655d7dd | 55 | { |
e37e3c84 | 56 | //Comment |
70e652a6 | 57 | |
58 | // snprintf(fName, 256,"%s", name); | |
59 | // snprintf(fNameShort,256, "%s", nameshort); | |
60 | ||
61 | sprintf(fName ,"%s", name); | |
62 | sprintf(fNameShort, "%s", nameshort); | |
63 | ||
92d9f317 | 64 | for(int i=0; i < ALTROMAXSAMPLES; i++ ) |
d655d7dd | 65 | { |
66 | fReversed[i] = 0; | |
67 | } | |
92d9f317 | 68 | |
69 | // fFp = fopen("amp2.txt", "w"); | |
70 | ||
d655d7dd | 71 | } |
72 | ||
57839add | 73 | AliCaloRawAnalyzer::~AliCaloRawAnalyzer() |
d655d7dd | 74 | { |
75 | ||
76 | } | |
77 | ||
78 | ||
79 | void | |
57839add | 80 | AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max ) |
d655d7dd | 81 | { |
82 | //Require that the bin if the maximum ADC value is between min and max (timebin) | |
92d9f317 | 83 | if( ( min > max ) || min > ALTROMAXSAMPLES || max > ALTROMAXSAMPLES ) |
d655d7dd | 84 | { |
85 | AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored", min, max ) ); | |
86 | } | |
87 | else | |
88 | { | |
89 | fMinTimeIndex = min; | |
90 | fMaxTimeIndex = max; | |
91 | } | |
92 | } | |
93 | ||
94 | ||
95 | UShort_t | |
92d9f317 | 96 | AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) |
d655d7dd | 97 | { |
98 | //------------ | |
99 | UShort_t tmpmax = data[0]; | |
100 | ||
101 | for(int i=0; i < length; i++) | |
102 | { | |
103 | if( tmpmax < data[i] ) | |
104 | { | |
105 | tmpmax = data[i]; | |
106 | } | |
107 | } | |
108 | return tmpmax; | |
109 | } | |
110 | ||
111 | ||
112 | void | |
92d9f317 | 113 | AliCaloRawAnalyzer::SelectSubarray( const Double_t *data, const int length, const short maxindex, int *const first, int *const last, const int cut) |
d655d7dd | 114 | { |
115 | //Selection of subset of data from one bunch that will be used for fitting or | |
116 | //Peak finding. Go to the left and right of index of the maximum time bin | |
e7acbc37 | 117 | //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump) |
d655d7dd | 118 | int tmpfirst = maxindex; |
119 | int tmplast = maxindex; | |
92d9f317 | 120 | Double_t prevFirst = data[maxindex]; |
121 | Double_t prevLast = data[maxindex]; | |
947c8e28 | 122 | bool firstJump = false; |
123 | bool lastJump = false; | |
124 | ||
92d9f317 | 125 | while( (tmpfirst >= 0) && (data[tmpfirst] >= cut ) && (!firstJump) ) |
d655d7dd | 126 | { |
947c8e28 | 127 | // jump check: |
128 | if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex | |
92d9f317 | 129 | if ( data[tmpfirst] >= prevFirst) { |
947c8e28 | 130 | firstJump = true; |
131 | } | |
132 | } | |
92d9f317 | 133 | prevFirst = data[tmpfirst]; |
d655d7dd | 134 | tmpfirst -- ; |
135 | } | |
136 | ||
92d9f317 | 137 | while( (tmplast < length) && (data[tmplast] >= cut ) && (!lastJump) ) |
d655d7dd | 138 | { |
947c8e28 | 139 | // jump check: |
140 | if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex | |
92d9f317 | 141 | if ( data[tmplast] >= prevLast) { |
947c8e28 | 142 | lastJump = true; |
143 | } | |
144 | } | |
92d9f317 | 145 | prevLast = data[tmplast]; |
d655d7dd | 146 | tmplast ++; |
147 | } | |
f57baa2d | 148 | |
947c8e28 | 149 | // we keep one pre- or post- sample if we can (as in online) |
150 | // check though if we ended up on a 'jump', or out of bounds: if so, back up | |
151 | if (firstJump || tmpfirst<0) tmpfirst ++; | |
152 | if (lastJump || tmplast>=length) tmplast --; | |
153 | ||
154 | *first = tmpfirst; | |
155 | *last = tmplast; | |
e7acbc37 | 156 | return; |
d655d7dd | 157 | } |
158 | ||
159 | ||
160 | ||
161 | Float_t | |
57839add | 162 | AliCaloRawAnalyzer::ReverseAndSubtractPed( const AliCaloBunchInfo *bunch, const UInt_t /*altrocfg1*/, const UInt_t /*altrocfg2*/, double *outarray ) const |
d655d7dd | 163 | { |
164 | //Time sample comes in reversed order, revers them back | |
165 | //Subtract the baseline based on content of altrocfg1 and altrocfg2. | |
166 | Int_t length = bunch->GetLength(); | |
167 | const UShort_t *sig = bunch->GetData(); | |
168 | ||
57839add | 169 | double ped = EvaluatePedestal( sig , length); |
170 | ||
171 | for( int i=0; i < length; i++ ) | |
172 | { | |
173 | outarray[i] = sig[length -i -1] - ped; | |
174 | } | |
175 | ||
176 | return ped; | |
177 | } | |
178 | ||
179 | ||
180 | ||
181 | Float_t | |
182 | AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * const data, const int /*length*/ ) const | |
183 | { | |
184 | // double ped = 0; | |
d655d7dd | 185 | double tmp = 0; |
186 | ||
187 | if( fIsZerosupressed == false ) | |
188 | { | |
f57baa2d | 189 | for(int i=0; i < fNsamplePed; i++ ) |
d655d7dd | 190 | { |
57839add | 191 | tmp += data[i]; |
d655d7dd | 192 | } |
f57baa2d | 193 | } |
e37e3c84 | 194 | |
f57baa2d | 195 | return tmp/fNsamplePed; |
196 | } | |
d655d7dd | 197 | |
198 | ||
199 | short | |
92d9f317 | 200 | AliCaloRawAnalyzer::Max( const AliCaloBunchInfo *const bunch , int *const maxindex ) |
d655d7dd | 201 | { |
202 | //comment | |
203 | short tmpmax = -1; | |
204 | int tmpindex = -1; | |
205 | const UShort_t *sig = bunch->GetData(); | |
206 | ||
207 | for(int i=0; i < bunch->GetLength(); i++ ) | |
208 | { | |
209 | if( sig[i] > tmpmax ) | |
210 | { | |
211 | tmpmax = sig[i]; | |
212 | tmpindex = i; | |
213 | } | |
214 | } | |
215 | ||
216 | if(maxindex != 0 ) | |
217 | { | |
48a2e3eb | 218 | // *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin(); |
219 | *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin(); | |
d655d7dd | 220 | } |
221 | ||
222 | return tmpmax; | |
223 | } | |
224 | ||
225 | ||
507751ce | 226 | bool |
227 | AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const | |
228 | { | |
229 | // a bunch is considered invalid if the maximum is in the first or last time-bin | |
230 | short tmpmax = -1; | |
231 | int tmpindex = -1; | |
232 | const UShort_t *sig = bunch->GetData(); | |
233 | ||
234 | for(int i=0; i < bunch->GetLength(); i++ ) | |
235 | { | |
236 | if( sig[i] > tmpmax ) | |
237 | { | |
238 | tmpmax = sig[i]; | |
239 | tmpindex = i; | |
240 | } | |
241 | } | |
242 | ||
243 | bool bunchOK = true; | |
244 | if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) ) | |
245 | { | |
246 | bunchOK = false; | |
247 | } | |
248 | ||
249 | return bunchOK; | |
250 | } | |
251 | ||
252 | ||
d655d7dd | 253 | int |
92d9f317 | 254 | AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,short *const maxampbin, short *const maxamplitude ) |
d655d7dd | 255 | { |
256 | //We select the bunch with the highest amplitude unless any time constraints is set | |
257 | short max = -1; | |
258 | short bunchindex = -1; | |
259 | short maxall = -1; | |
260 | int indx = -1; | |
261 | ||
262 | for(unsigned int i=0; i < bunchvector.size(); i++ ) | |
263 | { | |
48a2e3eb | 264 | max = Max( &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches |
92d9f317 | 265 | if( IsInTimeRange( indx, fMaxTimeIndex, fMinTimeIndex) ) |
d655d7dd | 266 | { |
267 | if( max > maxall ) | |
268 | { | |
269 | maxall = max; | |
270 | bunchindex = i; | |
48a2e3eb | 271 | *maxampbin = indx; |
272 | *maxamplitude = max; | |
d655d7dd | 273 | } |
274 | } | |
275 | } | |
276 | ||
507751ce | 277 | if (bunchindex >= 0) { |
278 | bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) ); | |
279 | if (! bunchOK) { | |
280 | bunchindex = -1; | |
281 | } | |
282 | } | |
283 | ||
d655d7dd | 284 | return bunchindex; |
285 | } | |
286 | ||
287 | ||
288 | bool | |
92d9f317 | 289 | AliCaloRawAnalyzer::IsInTimeRange( const int maxindex, const int maxtindx, const int mintindx ) |
d655d7dd | 290 | { |
291 | // Ckeck if the index of the max ADC vaue is consistent with trigger. | |
92d9f317 | 292 | if( ( mintindx < 0 && maxtindx < 0) ||maxtindx < 0 ) |
d655d7dd | 293 | { |
294 | return true; | |
295 | } | |
92d9f317 | 296 | return ( maxindex < maxtindx ) && ( maxindex > mintindx ) ? true : false; |
d655d7dd | 297 | } |
298 | ||
299 | ||
300 | void | |
57839add | 301 | AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) const |
d655d7dd | 302 | { |
303 | //comment | |
304 | cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl; | |
305 | cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl; | |
306 | ||
307 | for(unsigned int i=0; i < bvctr.size() ; i++ ) | |
308 | { | |
309 | PrintBunch( bvctr.at(i) ); | |
310 | cout << " bunch = " << i << endl; | |
311 | } | |
312 | cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl; | |
313 | } | |
314 | ||
315 | ||
316 | void | |
57839add | 317 | AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) const |
d655d7dd | 318 | { |
319 | //comment | |
320 | cout << __FILE__ << ":" << __LINE__ << endl; | |
321 | cout << __FILE__ << __LINE__ << ", startimebin = " << bunch.GetStartBin() << ", length =" << bunch.GetLength() << endl; | |
322 | const UShort_t *sig = bunch.GetData(); | |
323 | ||
324 | for ( Int_t j = 0; j < bunch.GetLength(); j++) | |
325 | { | |
326 | printf("%d\t", sig[j] ); | |
327 | } | |
328 | cout << endl; | |
329 | } | |
330 | ||
48a2e3eb | 331 | |
947c8e28 | 332 | Double_t |
333 | AliCaloRawAnalyzer::CalculateChi2(const Double_t amp, const Double_t time, | |
334 | const Int_t first, const Int_t last, | |
335 | const Double_t adcErr, | |
336 | const Double_t tau) | |
337 | { | |
338 | // Input: | |
339 | // amp - max amplitude; | |
340 | // time - time of max amplitude; | |
341 | // first, last - sample array indices to be used | |
342 | // adcErr - nominal error of amplitude measurement (one value for all channels) | |
343 | // if adcErr<0 that mean adcErr=1. | |
344 | // tau - filter time response (in timebin units) | |
345 | // Output: | |
346 | // chi2 - chi2 | |
347 | ||
348 | if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. | |
349 | // or, first is negative, the indices are not valid | |
168c7b3c | 350 | return Ret::kDummy; |
947c8e28 | 351 | } |
352 | ||
353 | int nsamples = last - first + 1; | |
2cd0ffda | 354 | // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i : amp %3.2f time %3.2f \n", first, last, nsamples, amp, time); |
947c8e28 | 355 | |
356 | Int_t x = 0; | |
357 | Double_t chi2 = 0; | |
358 | Double_t dy = 0.0, xx = 0.0, f=0.0; | |
359 | ||
360 | for (Int_t i=0; i<nsamples; i++) { | |
361 | x = first + i; // timebin | |
362 | xx = (x - time + tau) / tau; // help variable | |
363 | f = 0; | |
364 | if (xx > 0) { | |
365 | f = amp * xx*xx * TMath::Exp(2 * (1 - xx )) ; | |
366 | } | |
367 | dy = fReversed[x] - f; | |
368 | chi2 += dy*dy; | |
2cd0ffda | 369 | // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy); |
947c8e28 | 370 | } |
371 | ||
372 | if (adcErr>0.0) { // weight chi2 | |
373 | chi2 /= (adcErr*adcErr); | |
374 | } | |
375 | return chi2; | |
376 | } | |
377 | ||
378 | ||
379 | void | |
380 | AliCaloRawAnalyzer::CalculateMeanAndRMS(const Int_t first, const Int_t last, | |
381 | Double_t & mean, Double_t & rms) | |
382 | { | |
383 | // Input: | |
384 | // first, last - sample array indices to be used | |
385 | // Output: | |
386 | // mean and RMS of samples | |
387 | // | |
388 | // To possibly be used to differentiate good signals from bad before fitting | |
389 | // | |
168c7b3c | 390 | mean = Ret::kDummy; |
391 | rms = Ret::kDummy; | |
947c8e28 | 392 | |
393 | if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. | |
394 | // or, first is negative, the indices are not valid | |
395 | return; | |
396 | } | |
397 | ||
398 | int nsamples = last - first + 1; | |
399 | // printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples); | |
400 | ||
401 | int x = 0; | |
402 | Double_t sampleSum = 0; // sum of samples | |
403 | Double_t squaredSampleSum = 0; // sum of samples squared | |
404 | ||
405 | for (Int_t i=0; i<nsamples; i++) { | |
406 | x = first + i; | |
407 | sampleSum += fReversed[x]; | |
408 | squaredSampleSum += (fReversed[x] * fReversed[x]); | |
409 | } | |
410 | ||
411 | mean = sampleSum / nsamples; | |
412 | Double_t squaredMean = squaredSampleSum / nsamples; | |
413 | // The variance (rms squared) is equal to the mean of the squares minus the square of the mean.. | |
414 | rms = sqrt(squaredMean - mean*mean); | |
415 | ||
416 | return; | |
417 | } | |
418 | ||
168c7b3c | 419 | |
420 | ||
421 | // AliCaloFitResults | |
422 | // AliCaloRawAnalyzer::Evaluate( const vector<AliCaloBunchInfo> &/*bunchvector*/, const UInt_t /*altrocfg1*/, const UInt_t /*altrocfg2*/) | |
423 | // { // method to do the selection of what should possibly be fitted | |
424 | // // not implemented for base class | |
425 | // cout << __FILE__ << ":" << __LINE__ << " " << endl; | |
426 | ||
427 | // return AliCaloFitResults( 0, 0 ); | |
428 | // } | |
429 | ||
57839add | 430 | |
48a2e3eb | 431 | |
57839add | 432 | int |
92d9f317 | 433 | AliCaloRawAnalyzer::PreFitEvaluateSamples( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, |
434 | const UInt_t altrocfg2, Int_t & index, Float_t & maxf, short & maxamp, | |
435 | short & maxrev, Float_t & ped, int & first, int & last,const int acut ) | |
57839add | 436 | { // method to do the selection of what should possibly be fitted |
437 | int nsamples = 0; | |
f57baa2d | 438 | short maxampindex = 0; |
57839add | 439 | index = SelectBunch( bunchvector, &maxampindex, &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set |
440 | ||
441 | ||
92d9f317 | 442 | if( index >= 0 && maxamp >= acut ) // something valid was found, and non-zero amplitude |
57839add | 443 | { |
444 | // use more convenient numbering and possibly subtract pedestal | |
445 | ped = ReverseAndSubtractPed( &(bunchvector.at(index)), altrocfg1, altrocfg2, fReversed ); | |
446 | maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed ); | |
447 | ||
92d9f317 | 448 | if ( maxf >= acut ) // possibly significant signal |
57839add | 449 | { |
450 | // select array around max to possibly be used in fit | |
f57baa2d | 451 | maxrev = maxampindex - bunchvector.at(index).GetStartBin(); |
92d9f317 | 452 | SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, acut ); |
57839add | 453 | |
454 | // sanity check: maximum should not be in first or last bin | |
455 | // if we should do a fit | |
92d9f317 | 456 | if (first!=maxrev && last!=maxrev) |
457 | { | |
458 | // calculate how many samples we have | |
459 | nsamples = last - first + 1; | |
460 | } | |
57839add | 461 | } |
462 | } | |
463 | ||
464 | return nsamples; | |
465 | } | |
466 |