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