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