]>
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), |
92d9f317 | 49 | fL1Phase(0), |
50 | fAmp(0), | |
396baaf6 | 51 | fTof(0), |
52 | fTau( EMCAL::TAU ) | |
d655d7dd | 53 | { |
e37e3c84 | 54 | //Comment |
576a5057 | 55 | snprintf(fName, 256,"%s", name); |
56 | snprintf(fNameShort,256, "%s", nameshort); | |
57 | // sprintf(fName ,"%s", name); | |
58 | // sprintf(fNameShort, "%s", nameshort); | |
70e652a6 | 59 | |
92d9f317 | 60 | for(int i=0; i < ALTROMAXSAMPLES; i++ ) |
d655d7dd | 61 | { |
62 | fReversed[i] = 0; | |
63 | } | |
64 | } | |
65 | ||
57839add | 66 | AliCaloRawAnalyzer::~AliCaloRawAnalyzer() |
d655d7dd | 67 | { |
68 | ||
69 | } | |
70 | ||
71 | ||
72 | void | |
57839add | 73 | AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max ) |
d655d7dd | 74 | { |
75 | //Require that the bin if the maximum ADC value is between min and max (timebin) | |
92d9f317 | 76 | if( ( min > max ) || min > ALTROMAXSAMPLES || max > ALTROMAXSAMPLES ) |
d655d7dd | 77 | { |
78 | AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored", min, max ) ); | |
79 | } | |
80 | else | |
81 | { | |
82 | fMinTimeIndex = min; | |
83 | fMaxTimeIndex = max; | |
84 | } | |
85 | } | |
86 | ||
87 | ||
88 | UShort_t | |
9dfd5c8d | 89 | AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) const |
d655d7dd | 90 | { |
91 | //------------ | |
92 | UShort_t tmpmax = data[0]; | |
93 | ||
94 | for(int i=0; i < length; i++) | |
95 | { | |
96 | if( tmpmax < data[i] ) | |
97 | { | |
98 | tmpmax = data[i]; | |
99 | } | |
100 | } | |
101 | return tmpmax; | |
102 | } | |
103 | ||
104 | ||
105 | void | |
9dfd5c8d | 106 | AliCaloRawAnalyzer::SelectSubarray( const Double_t *data, const int length, const short maxindex,int *const first, int *const last, const int cut) const |
d655d7dd | 107 | { |
108 | //Selection of subset of data from one bunch that will be used for fitting or | |
109 | //Peak finding. Go to the left and right of index of the maximum time bin | |
e7acbc37 | 110 | //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump) |
d655d7dd | 111 | int tmpfirst = maxindex; |
112 | int tmplast = maxindex; | |
92d9f317 | 113 | Double_t prevFirst = data[maxindex]; |
114 | Double_t prevLast = data[maxindex]; | |
947c8e28 | 115 | bool firstJump = false; |
116 | bool lastJump = false; | |
117 | ||
92d9f317 | 118 | while( (tmpfirst >= 0) && (data[tmpfirst] >= cut ) && (!firstJump) ) |
d655d7dd | 119 | { |
947c8e28 | 120 | // jump check: |
121 | if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex | |
92d9f317 | 122 | if ( data[tmpfirst] >= prevFirst) { |
947c8e28 | 123 | firstJump = true; |
124 | } | |
125 | } | |
92d9f317 | 126 | prevFirst = data[tmpfirst]; |
d655d7dd | 127 | tmpfirst -- ; |
128 | } | |
129 | ||
92d9f317 | 130 | while( (tmplast < length) && (data[tmplast] >= cut ) && (!lastJump) ) |
d655d7dd | 131 | { |
947c8e28 | 132 | // jump check: |
133 | if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex | |
92d9f317 | 134 | if ( data[tmplast] >= prevLast) { |
947c8e28 | 135 | lastJump = true; |
136 | } | |
137 | } | |
92d9f317 | 138 | prevLast = data[tmplast]; |
d655d7dd | 139 | tmplast ++; |
140 | } | |
f57baa2d | 141 | |
947c8e28 | 142 | // we keep one pre- or post- sample if we can (as in online) |
143 | // check though if we ended up on a 'jump', or out of bounds: if so, back up | |
144 | if (firstJump || tmpfirst<0) tmpfirst ++; | |
145 | if (lastJump || tmplast>=length) tmplast --; | |
146 | ||
147 | *first = tmpfirst; | |
148 | *last = tmplast; | |
e7acbc37 | 149 | return; |
d655d7dd | 150 | } |
151 | ||
152 | ||
153 | ||
154 | Float_t | |
57839add | 155 | AliCaloRawAnalyzer::ReverseAndSubtractPed( const AliCaloBunchInfo *bunch, const UInt_t /*altrocfg1*/, const UInt_t /*altrocfg2*/, double *outarray ) const |
d655d7dd | 156 | { |
157 | //Time sample comes in reversed order, revers them back | |
158 | //Subtract the baseline based on content of altrocfg1 and altrocfg2. | |
159 | Int_t length = bunch->GetLength(); | |
160 | const UShort_t *sig = bunch->GetData(); | |
161 | ||
57839add | 162 | double ped = EvaluatePedestal( sig , length); |
163 | ||
164 | for( int i=0; i < length; i++ ) | |
165 | { | |
166 | outarray[i] = sig[length -i -1] - ped; | |
167 | } | |
168 | ||
169 | return ped; | |
170 | } | |
171 | ||
172 | ||
173 | ||
174 | Float_t | |
175 | AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * const data, const int /*length*/ ) const | |
176 | { | |
177 | // double ped = 0; | |
d655d7dd | 178 | double tmp = 0; |
179 | ||
180 | if( fIsZerosupressed == false ) | |
181 | { | |
f57baa2d | 182 | for(int i=0; i < fNsamplePed; i++ ) |
d655d7dd | 183 | { |
57839add | 184 | tmp += data[i]; |
d655d7dd | 185 | } |
f57baa2d | 186 | } |
e37e3c84 | 187 | |
f57baa2d | 188 | return tmp/fNsamplePed; |
189 | } | |
d655d7dd | 190 | |
191 | ||
192 | short | |
9dfd5c8d | 193 | AliCaloRawAnalyzer::Max( const AliCaloBunchInfo *const bunch , int *const maxindex ) const |
d655d7dd | 194 | { |
195 | //comment | |
196 | short tmpmax = -1; | |
197 | int tmpindex = -1; | |
198 | const UShort_t *sig = bunch->GetData(); | |
199 | ||
200 | for(int i=0; i < bunch->GetLength(); i++ ) | |
201 | { | |
202 | if( sig[i] > tmpmax ) | |
203 | { | |
204 | tmpmax = sig[i]; | |
205 | tmpindex = i; | |
206 | } | |
207 | } | |
208 | ||
209 | if(maxindex != 0 ) | |
210 | { | |
48a2e3eb | 211 | // *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin(); |
212 | *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin(); | |
d655d7dd | 213 | } |
214 | ||
215 | return tmpmax; | |
216 | } | |
217 | ||
218 | ||
507751ce | 219 | bool |
220 | AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const | |
221 | { | |
222 | // a bunch is considered invalid if the maximum is in the first or last time-bin | |
223 | short tmpmax = -1; | |
224 | int tmpindex = -1; | |
225 | const UShort_t *sig = bunch->GetData(); | |
226 | ||
227 | for(int i=0; i < bunch->GetLength(); i++ ) | |
228 | { | |
229 | if( sig[i] > tmpmax ) | |
230 | { | |
231 | tmpmax = sig[i]; | |
232 | tmpindex = i; | |
233 | } | |
234 | } | |
235 | ||
236 | bool bunchOK = true; | |
237 | if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) ) | |
238 | { | |
239 | bunchOK = false; | |
240 | } | |
241 | ||
242 | return bunchOK; | |
243 | } | |
244 | ||
245 | ||
d655d7dd | 246 | int |
9dfd5c8d | 247 | AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector, short *const maxampbin, short *const maxamplitude ) |
d655d7dd | 248 | { |
249 | //We select the bunch with the highest amplitude unless any time constraints is set | |
250 | short max = -1; | |
251 | short bunchindex = -1; | |
252 | short maxall = -1; | |
253 | int indx = -1; | |
254 | ||
255 | for(unsigned int i=0; i < bunchvector.size(); i++ ) | |
256 | { | |
48a2e3eb | 257 | max = Max( &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches |
92d9f317 | 258 | if( IsInTimeRange( indx, fMaxTimeIndex, fMinTimeIndex) ) |
d655d7dd | 259 | { |
260 | if( max > maxall ) | |
261 | { | |
262 | maxall = max; | |
263 | bunchindex = i; | |
48a2e3eb | 264 | *maxampbin = indx; |
265 | *maxamplitude = max; | |
d655d7dd | 266 | } |
267 | } | |
268 | } | |
269 | ||
507751ce | 270 | if (bunchindex >= 0) { |
271 | bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) ); | |
272 | if (! bunchOK) { | |
273 | bunchindex = -1; | |
274 | } | |
275 | } | |
276 | ||
d655d7dd | 277 | return bunchindex; |
278 | } | |
279 | ||
280 | ||
281 | bool | |
9dfd5c8d | 282 | AliCaloRawAnalyzer::IsInTimeRange( const int maxindex, const int maxtindx, const int mintindx ) const |
d655d7dd | 283 | { |
284 | // Ckeck if the index of the max ADC vaue is consistent with trigger. | |
92d9f317 | 285 | if( ( mintindx < 0 && maxtindx < 0) ||maxtindx < 0 ) |
d655d7dd | 286 | { |
287 | return true; | |
288 | } | |
92d9f317 | 289 | return ( maxindex < maxtindx ) && ( maxindex > mintindx ) ? true : false; |
d655d7dd | 290 | } |
291 | ||
292 | ||
293 | void | |
f67ce0df | 294 | AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) |
d655d7dd | 295 | { |
296 | //comment | |
297 | cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl; | |
298 | cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl; | |
299 | ||
300 | for(unsigned int i=0; i < bvctr.size() ; i++ ) | |
301 | { | |
302 | PrintBunch( bvctr.at(i) ); | |
303 | cout << " bunch = " << i << endl; | |
304 | } | |
305 | cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl; | |
306 | } | |
307 | ||
308 | ||
309 | void | |
f67ce0df | 310 | AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) |
d655d7dd | 311 | { |
312 | //comment | |
313 | cout << __FILE__ << ":" << __LINE__ << endl; | |
314 | cout << __FILE__ << __LINE__ << ", startimebin = " << bunch.GetStartBin() << ", length =" << bunch.GetLength() << endl; | |
315 | const UShort_t *sig = bunch.GetData(); | |
316 | ||
317 | for ( Int_t j = 0; j < bunch.GetLength(); j++) | |
318 | { | |
319 | printf("%d\t", sig[j] ); | |
320 | } | |
321 | cout << endl; | |
322 | } | |
323 | ||
48a2e3eb | 324 | |
947c8e28 | 325 | Double_t |
326 | AliCaloRawAnalyzer::CalculateChi2(const Double_t amp, const Double_t time, | |
327 | const Int_t first, const Int_t last, | |
328 | const Double_t adcErr, | |
9dfd5c8d | 329 | const Double_t tau) const |
947c8e28 | 330 | { |
331 | // Input: | |
332 | // amp - max amplitude; | |
333 | // time - time of max amplitude; | |
334 | // first, last - sample array indices to be used | |
335 | // adcErr - nominal error of amplitude measurement (one value for all channels) | |
336 | // if adcErr<0 that mean adcErr=1. | |
337 | // tau - filter time response (in timebin units) | |
338 | // Output: | |
339 | // chi2 - chi2 | |
340 | ||
341 | if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. | |
342 | // or, first is negative, the indices are not valid | |
168c7b3c | 343 | return Ret::kDummy; |
947c8e28 | 344 | } |
345 | ||
346 | int nsamples = last - first + 1; | |
2cd0ffda | 347 | // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i : amp %3.2f time %3.2f \n", first, last, nsamples, amp, time); |
947c8e28 | 348 | |
349 | Int_t x = 0; | |
350 | Double_t chi2 = 0; | |
351 | Double_t dy = 0.0, xx = 0.0, f=0.0; | |
352 | ||
353 | for (Int_t i=0; i<nsamples; i++) { | |
354 | x = first + i; // timebin | |
355 | xx = (x - time + tau) / tau; // help variable | |
356 | f = 0; | |
357 | if (xx > 0) { | |
358 | f = amp * xx*xx * TMath::Exp(2 * (1 - xx )) ; | |
359 | } | |
360 | dy = fReversed[x] - f; | |
361 | chi2 += dy*dy; | |
2cd0ffda | 362 | // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy); |
947c8e28 | 363 | } |
364 | ||
365 | if (adcErr>0.0) { // weight chi2 | |
366 | chi2 /= (adcErr*adcErr); | |
367 | } | |
368 | return chi2; | |
369 | } | |
370 | ||
371 | ||
372 | void | |
373 | AliCaloRawAnalyzer::CalculateMeanAndRMS(const Int_t first, const Int_t last, | |
374 | Double_t & mean, Double_t & rms) | |
375 | { | |
376 | // Input: | |
377 | // first, last - sample array indices to be used | |
378 | // Output: | |
379 | // mean and RMS of samples | |
380 | // | |
381 | // To possibly be used to differentiate good signals from bad before fitting | |
382 | // | |
168c7b3c | 383 | mean = Ret::kDummy; |
384 | rms = Ret::kDummy; | |
947c8e28 | 385 | |
386 | if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined.. | |
387 | // or, first is negative, the indices are not valid | |
388 | return; | |
389 | } | |
390 | ||
391 | int nsamples = last - first + 1; | |
392 | // printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples); | |
393 | ||
394 | int x = 0; | |
395 | Double_t sampleSum = 0; // sum of samples | |
396 | Double_t squaredSampleSum = 0; // sum of samples squared | |
397 | ||
398 | for (Int_t i=0; i<nsamples; i++) { | |
399 | x = first + i; | |
400 | sampleSum += fReversed[x]; | |
401 | squaredSampleSum += (fReversed[x] * fReversed[x]); | |
402 | } | |
403 | ||
404 | mean = sampleSum / nsamples; | |
405 | Double_t squaredMean = squaredSampleSum / nsamples; | |
406 | // The variance (rms squared) is equal to the mean of the squares minus the square of the mean.. | |
407 | rms = sqrt(squaredMean - mean*mean); | |
408 | ||
409 | return; | |
410 | } | |
411 | ||
168c7b3c | 412 | |
413 | ||
57839add | 414 | int |
92d9f317 | 415 | AliCaloRawAnalyzer::PreFitEvaluateSamples( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, |
416 | const UInt_t altrocfg2, Int_t & index, Float_t & maxf, short & maxamp, | |
417 | short & maxrev, Float_t & ped, int & first, int & last,const int acut ) | |
57839add | 418 | { // method to do the selection of what should possibly be fitted |
419 | int nsamples = 0; | |
f57baa2d | 420 | short maxampindex = 0; |
57839add | 421 | index = SelectBunch( bunchvector, &maxampindex, &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set |
422 | ||
423 | ||
92d9f317 | 424 | if( index >= 0 && maxamp >= acut ) // something valid was found, and non-zero amplitude |
57839add | 425 | { |
426 | // use more convenient numbering and possibly subtract pedestal | |
427 | ped = ReverseAndSubtractPed( &(bunchvector.at(index)), altrocfg1, altrocfg2, fReversed ); | |
428 | maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed ); | |
429 | ||
92d9f317 | 430 | if ( maxf >= acut ) // possibly significant signal |
57839add | 431 | { |
432 | // select array around max to possibly be used in fit | |
f57baa2d | 433 | maxrev = maxampindex - bunchvector.at(index).GetStartBin(); |
92d9f317 | 434 | SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, acut ); |
57839add | 435 | |
436 | // sanity check: maximum should not be in first or last bin | |
437 | // if we should do a fit | |
92d9f317 | 438 | if (first!=maxrev && last!=maxrev) |
439 | { | |
440 | // calculate how many samples we have | |
441 | nsamples = last - first + 1; | |
442 | } | |
57839add | 443 | } |
444 | } | |
445 | ||
446 | return nsamples; | |
447 | } | |
448 |