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