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