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