AliRoot Core  b857a14 (b857a14)
AliCaloRawAnalyzer.cxx
Go to the documentation of this file.
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 // AliRoot/EMCal system
21 #include "AliCaloRawAnalyzer.h"
22 #include "AliCaloBunchInfo.h"
23 #include "AliCaloFitResults.h"
24 #include "AliLog.h"
25 
26 // ROOT sytem
27 #include "TMath.h"
28 
29 // Standard libraries
30 #include <iostream>
31 using namespace std;
32 
34 ClassImp(AliCaloRawAnalyzer) ;
36 
39 //_______________________________________________________________________
40 AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char *name, const char *nameshort) : TObject(),
41  fMinTimeIndex(-1),
42  fMaxTimeIndex(-1),
43  fFitArrayCut(5),
44  fAmpCut(4),
45  fNsampleCut(5),
46  fOverflowCut(950),
47  fNsamplePed(3),
48  fIsZerosupressed( false ),
49  fVerbose( false ),
50  fAlgo(Algo::kNONE),
51  fL1Phase(0),
52  fAmp(0),
53  fTof(0),
54  fTau( EMCAL::TAU ),
55  fFixTau( true )
56 {
57  snprintf(fName, 256, "%s", name);
58  snprintf(fNameShort,256, "%s", nameshort);
59 
60  for(int i=0; i < ALTROMAXSAMPLES; i++ )
61  {
62  fReversed[i] = 0;
63  }
64 }
65 
68 //_______________________________________________________________________
69 void
70 AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max )
71 {
72 
73  if( ( min > max ) || min > ALTROMAXSAMPLES || max > ALTROMAXSAMPLES )
74  {
75  AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored", min, max ) );
76  }
77  else
78  {
79  fMinTimeIndex = min;
80  fMaxTimeIndex = max;
81  }
82 }
83 
86 //_______________________________________________________________________
87 UShort_t
88 AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) const
89 {
90  UShort_t tmpmax = data[0];
91 
92  for(int i=0; i < length; i++)
93  {
94  if( tmpmax < data[i] )
95  {
96  tmpmax = data[i];
97  }
98  }
99  return tmpmax;
100 }
101 
106 //_______________________________________________________________________
107 void
108 AliCaloRawAnalyzer::SelectSubarray( const Double_t *data, int length, short maxindex,
109  int * first, int * last, int cut) const
110 {
111  int tmpfirst = maxindex;
112  int tmplast = maxindex;
113  Double_t prevFirst = data[maxindex];
114  Double_t prevLast = data[maxindex];
115  bool firstJump = false;
116  bool lastJump = false;
117 
118  while( (tmpfirst >= 0) && (data[tmpfirst] >= cut ) && (!firstJump) )
119  {
120  // jump check:
121  if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex
122  if ( data[tmpfirst] >= prevFirst) {
123  firstJump = true;
124  }
125  }
126  prevFirst = data[tmpfirst];
127  tmpfirst -- ;
128  }
129 
130  while( (tmplast < length) && (data[tmplast] >= cut ) && (!lastJump) )
131  {
132  // jump check:
133  if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex
134  if ( data[tmplast] >= prevLast) {
135  lastJump = true;
136  }
137  }
138  prevLast = data[tmplast];
139  tmplast ++;
140  }
141 
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;
149 
150  return;
151 }
152 
153 
157 //_______________________________________________________________________
158 Float_t
160  UInt_t /*altrocfg1*/, UInt_t /*altrocfg2*/,
161  double *outarray ) const
162 {
163  Int_t length = bunch->GetLength();
164  const UShort_t *sig = bunch->GetData();
165 
166  double ped = EvaluatePedestal( sig , length);
167 
168  for( int i=0; i < length; i++ )
169  {
170  outarray[i] = sig[length -i -1] - ped;
171  }
172 
173  return ped;
174 }
175 
178 //_______________________________________________________________________
179 Float_t
180 AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * data, int /*length*/ ) const
181 {
182  double tmp = 0;
183 
184  if( fIsZerosupressed == false )
185  {
186  for(int i=0; i < fNsamplePed; i++ )
187  {
188  tmp += data[i];
189  }
190  }
191 
192  return tmp/fNsamplePed;
193 }
194 
197 //_______________________________________________________________________
198 short
199 AliCaloRawAnalyzer::Max( const AliCaloBunchInfo * bunch , int * maxindex ) const
200 {
201  short tmpmax = -1;
202  int tmpindex = -1;
203  const UShort_t *sig = bunch->GetData();
204 
205  for(int i=0; i < bunch->GetLength(); i++ )
206  {
207  if( sig[i] > tmpmax )
208  {
209  tmpmax = sig[i];
210  tmpindex = i;
211  }
212  }
213 
214  if(maxindex != 0 )
215  {
216  // *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
217  *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
218  }
219 
220  return tmpmax;
221 }
222 
225 //_______________________________________________________________________
226 bool
228 {
229  short tmpmax = -1;
230  int tmpindex = -1;
231  const UShort_t *sig = bunch->GetData();
232 
233  for(int i=0; i < bunch->GetLength(); i++ )
234  {
235  if( sig[i] > tmpmax )
236  {
237  tmpmax = sig[i];
238  tmpindex = i;
239  }
240  }
241 
242  bool bunchOK = true;
243  if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) )
244  {
245  bunchOK = false;
246  }
247 
248  return bunchOK;
249 }
250 
253 //_______________________________________________________________________
254 int
255 AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,
256  short * maxampbin, short * maxamplitude )
257 {
258  short max = -1;
259  short bunchindex = -1;
260  short maxall = -1;
261  int indx = -1;
262 
263  for(unsigned int i=0; i < bunchvector.size(); i++ )
264  {
265  max = Max( &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches
267  {
268  if( max > maxall )
269  {
270  maxall = max;
271  bunchindex = i;
272  *maxampbin = indx;
273  *maxamplitude = max;
274  }
275  }
276  }
277 
278  if (bunchindex >= 0)
279  {
280  bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) );
281  if (! bunchOK)
282  {
283  bunchindex = -1;
284  }
285  }
286 
287  return bunchindex;
288 }
289 
292 //_______________________________________________________________________
293 bool
294 AliCaloRawAnalyzer::IsInTimeRange( const int maxindex, const int maxtindx, const int mintindx ) const
295 {
296  if( ( mintindx < 0 && maxtindx < 0) ||maxtindx < 0 )
297  {
298  return true;
299  }
300 
301  return ( maxindex < maxtindx ) && ( maxindex > mintindx ) ? true : false;
302 }
303 
306 //_______________________________________________________________________
307 void
308 AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr )
309 {
310  cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl;
311  cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl;
312 
313  for(unsigned int i=0; i < bvctr.size() ; i++ )
314  {
315  PrintBunch( bvctr.at(i) );
316  cout << " bunch = " << i << endl;
317  }
318  cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl;
319 }
320 
323 //_______________________________________________________________________
324 void
326 {
327  cout << __FILE__ << ":" << __LINE__ << endl;
328  cout << __FILE__ << __LINE__ << ", startimebin = " << bunch.GetStartBin() << ", length =" << bunch.GetLength() << endl;
329  const UShort_t *sig = bunch.GetData();
330 
331  for ( Int_t j = 0; j < bunch.GetLength(); j++)
332  {
333  printf("%d\t", sig[j] );
334  }
335  cout << endl;
336 }
337 
347 //_______________________________________________________________________
348 Double_t
349 AliCaloRawAnalyzer::CalculateChi2( Double_t amp, Double_t time,
350  Int_t first, Int_t last,
351  Double_t adcErr, Double_t tau) const
352 {
353  if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined..
354  // or, first is negative, the indices are not valid
355  return Ret::kDummy;
356  }
357 
358  int nsamples = last - first + 1;
359  // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i : amp %3.2f time %3.2f \n", first, last, nsamples, amp, time);
360 
361  Int_t x = 0;
362  Double_t chi2 = 0;
363  Double_t dy = 0.0, xx = 0.0, f=0.0;
364 
365  for (Int_t i=0; i<nsamples; i++)
366  {
367  x = first + i; // timebin
368  xx = (x - time + tau) / tau; // help variable
369  f = 0;
370 
371  if (xx > 0)
372  {
373  f = amp * xx*xx * TMath::Exp(2 * (1 - xx )) ;
374  }
375 
376  dy = fReversed[x] - f;
377  chi2 += dy*dy;
378  // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy);
379  }
380 
381  if (adcErr>0.0)
382  { // weight chi2
383  chi2 /= (adcErr*adcErr);
384  }
385 
386  return chi2;
387 }
388 
398 //_______________________________________________________________________
399 void
401  Double_t & mean, Double_t & rms)
402 {
403  mean = Ret::kDummy;
404  rms = Ret::kDummy;
405 
406  // signal consists of single sample, chi2 estimate (0) not too well defined..
407  // or, first is negative, the indices are not valid
408  if (first == last || first<0 )
409  return;
410 
411  int nsamples = last - first + 1;
412  // printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples);
413 
414  int x = 0;
415  Double_t sampleSum = 0; // sum of samples
416  Double_t squaredSampleSum = 0; // sum of samples squared
417 
418  for (Int_t i=0; i<nsamples; i++)
419  {
420  x = first + i;
421  sampleSum += fReversed[x];
422  squaredSampleSum += (fReversed[x] * fReversed[x]);
423  }
424 
425  mean = sampleSum / nsamples;
426  Double_t squaredMean = squaredSampleSum / nsamples;
427  // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
428  rms = sqrt(squaredMean - mean*mean);
429 
430  return;
431 }
432 
433 
436 //_______________________________________________________________________
437 int
438 AliCaloRawAnalyzer::PreFitEvaluateSamples( const vector<AliCaloBunchInfo> &bunchvector,
439  UInt_t altrocfg1, UInt_t altrocfg2, Int_t & index,
440  Float_t & maxf, short & maxamp,
441  short & maxrev, Float_t & ped,
442  int & first, int & last, int acut )
443 {
444  int nsamples = 0;
445  short maxampindex = 0;
446 
447  // select the bunch with the highest amplitude unless any time constraints is set
448  index = SelectBunch( bunchvector, &maxampindex, &maxamp );
449 
450  // something valid was found, and non-zero amplitude
451  if( index >= 0 && maxamp >= acut )
452  {
453  // use more convenient numbering and possibly subtract pedestal
454  ped = ReverseAndSubtractPed( &(bunchvector.at(index)), altrocfg1, altrocfg2, fReversed );
455  maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
456 
457  if ( maxf >= acut ) // possibly significant signal
458  {
459  // select array around max to possibly be used in fit
460  maxrev = maxampindex - bunchvector.at(index).GetStartBin();
461  SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, acut );
462 
463  // sanity check: maximum should not be in first or last bin
464  // if we should do a fit
465  if (first!=maxrev && last!=maxrev)
466  {
467  // calculate how many samples we have
468  nsamples = last - first + 1;
469  }
470  }
471  }
472 
473  return nsamples;
474 }
475 
int fMaxTimeIndex
The timebin of the max signal value must be between fMinTimeIndex and fMaxTimeIndex.
bool IsInTimeRange(const int maxindex, const int maxtime, const int mintime) const
Check if the index of the max ADC vaue is consistent with trigger.
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
Base class for extraction of signal amplitude and peak position.
Int_t GetLength() const
bool CheckBunchEdgesForMax(const AliCaloBunchInfo *const bunch) const
A bunch is considered invalid if the maximum is in the first or last time-bin.
char fName[256]
Name of the algorithm.
const UShort_t * GetData() const
void CalculateMeanAndRMS(const Int_t first, const Int_t last, Double_t &mean, Double_t &rms)
short Max(const AliCaloBunchInfo *const bunch, int *maxindex) const
Get maximum in bunch array.
AliCaloRawAnalyzer(const char *name="AliCaloRawAnalyzer", const char *nameshort="RawAna")
Constructor.
int SelectBunch(const std::vector< AliCaloBunchInfo > &bunchvector, short *maxampbin, short *maxamplitude)
We select the bunch with the highest amplitude unless any time constraints is set.
const double TAU
Approximate shaping time.
#define AliWarning(message)
Definition: AliLog.h:541
Container of ALTRO information.
Double_t chi2
Definition: AnalyzeLaser.C:7
bool fIsZerosupressed
Wether or not the data is zeros supressed, by default its assumed that the baseline is also subtracte...
static void PrintBunches(const std::vector< AliCaloBunchInfo > &bunchvector)
Print bunch vector infomation.
const Double_t fTau
char fNameShort[256]
Abbrevation for the name.
UInt_t GetStartBin() const
void SelectSubarray(const Double_t *date, int length, short maxindex, int *first, int *last, int cut) const
int fNsamplePed
Number of samples used for pedestal calculation (first in bunch)
Float_t ReverseAndSubtractPed(const AliCaloBunchInfo *bunch, UInt_t altrocfg1, UInt_t altrocfg2, double *outarray) const
int PreFitEvaluateSamples(const std::vector< AliCaloBunchInfo > &bunchvector, UInt_t altrocfg1, UInt_t altrocfg2, Int_t &index, Float_t &maxf, short &maxamp, short &maxampindex, Float_t &ped, int &first, int &last, int acut)
Method to do the selection of what should possibly be fitted.
const int ALTROMAXSAMPLES
The maximum number of samples of the ALTRO.
static void PrintBunch(const AliCaloBunchInfo &bunch)
Print bunch information.
TF1 * f
Definition: interpolTest.C:21
Double_t fReversed[ALTROMAXSAMPLES]
Reversed sequence of samples (pedestalsubtracted)
Double_t CalculateChi2(const Double_t amp, const Double_t time, const Int_t first, const Int_t last, const Double_t adcErr=1, const Double_t tau=2.35) const
Float_t EvaluatePedestal(const UShort_t *const data, const int length) const
Pedestal evaluation if not zero suppressed.
int fMinTimeIndex
The timebin of the max signal value must be between fMinTimeIndex and fMaxTimeIndex.
void SetTimeConstraint(int min, int max)
Require that the bin if the maximum ADC value is between min and max (timebin)
TCut cut
Definition: MakeGlobalFit.C:75
static double fAmp
amplitude