coverity
[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 ),
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 66AliCaloRawAnalyzer::~AliCaloRawAnalyzer()
d655d7dd 67{
68
69}
70
71
72void
57839add 73AliCaloRawAnalyzer::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
88UShort_t
9dfd5c8d 89AliCaloRawAnalyzer::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
105void
9dfd5c8d 106AliCaloRawAnalyzer::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
154Float_t
57839add 155AliCaloRawAnalyzer::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
174Float_t
175AliCaloRawAnalyzer::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
192short
9dfd5c8d 193AliCaloRawAnalyzer::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 219bool
220AliCaloRawAnalyzer::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 246int
9dfd5c8d 247AliCaloRawAnalyzer::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
281bool
9dfd5c8d 282AliCaloRawAnalyzer::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
293void
f67ce0df 294AliCaloRawAnalyzer::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
309void
f67ce0df 310AliCaloRawAnalyzer::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 325Double_t
326AliCaloRawAnalyzer::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
372void
373AliCaloRawAnalyzer::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 414int
92d9f317 415AliCaloRawAnalyzer::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