MUON + CheckCompiler
[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),
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
67void
57839add 68AliCaloRawAnalyzer::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
84UShort_t
9dfd5c8d 85AliCaloRawAnalyzer::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
102void
8a3ad886 103AliCaloRawAnalyzer::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
153Float_t
8a3ad886 154AliCaloRawAnalyzer::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
176Float_t
8a3ad886 177AliCaloRawAnalyzer::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
195short
8a3ad886 196AliCaloRawAnalyzer::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 223bool
224AliCaloRawAnalyzer::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 250int
8a3ad886 251AliCaloRawAnalyzer::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
286bool
9dfd5c8d 287AliCaloRawAnalyzer::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
298void
f67ce0df 299AliCaloRawAnalyzer::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
314void
f67ce0df 315AliCaloRawAnalyzer::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 330Double_t
8a3ad886 331AliCaloRawAnalyzer::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
376void
8a3ad886 377AliCaloRawAnalyzer::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 418int
8a3ad886 419AliCaloRawAnalyzer::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