1 /**************************************************************************
2 * This file is property of and copyright by *
3 * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 *
5 * Primary Author: Per Thomas Hille <perthomas.hille@yale.edu> *
7 * Contributors are mentioned in the code where appropriate. *
8 * Please report bugs to p.t.hille@fys.uio.no *
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 **************************************************************************/
20 // Base class for extraction
21 // of signal amplitude and peak position
22 // From CALO Calorimeter RAW data (from the RCU)
23 // Contains some utilities for preparing / selecting
24 // Signals suitable for signal extraction
28 #include "AliCaloRawAnalyzer.h"
29 #include "AliCaloBunchInfo.h"
30 #include "AliCaloFitResults.h"
35 ClassImp(AliCaloRawAnalyzer)
37 AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char *name, const char *nameshort) : TObject(),
44 fIsZerosupressed( false ),
48 sprintf(fName, "%s", name);
49 sprintf(fNameShort, "%s", nameshort);
51 for(int i=0; i < MAXSAMPLES; i++ )
57 AliCaloRawAnalyzer::~AliCaloRawAnalyzer()
64 AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max )
66 //Require that the bin if the maximum ADC value is between min and max (timebin)
67 if( ( min > max ) || min > MAXSAMPLES || max > MAXSAMPLES )
69 AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored", min, max ) );
80 AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) const
83 UShort_t tmpmax = data[0];
85 for(int i=0; i < length; i++)
87 if( tmpmax < data[i] )
97 AliCaloRawAnalyzer::SelectSubarray( const Double_t *fData, const int length, const short maxindex, int *const first, int *const last ) const
99 //Selection of subset of data from one bunch that will be used for fitting or
100 //Peak finding. Go to the left and right of index of the maximum time bin
101 //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump)
102 int tmpfirst = maxindex;
103 int tmplast = maxindex;
104 Double_t prevFirst = fData[maxindex];
105 Double_t prevLast = fData[maxindex];
106 bool firstJump = false;
107 bool lastJump = false;
109 while( (tmpfirst >= 0) && (fData[tmpfirst] >= fFitArrayCut) && (!firstJump) )
112 if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex
113 if (fData[tmpfirst] >= prevFirst) {
117 prevFirst = fData[tmpfirst];
121 while( (tmplast < length) && (fData[tmplast] >= fFitArrayCut) && (!lastJump) )
124 if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex
125 if (fData[tmplast] >= prevLast) {
129 prevLast = fData[tmplast];
133 // we keep one pre- or post- sample if we can (as in online)
134 // check though if we ended up on a 'jump', or out of bounds: if so, back up
135 if (firstJump || tmpfirst<0) tmpfirst ++;
136 if (lastJump || tmplast>=length) tmplast --;
146 AliCaloRawAnalyzer::ReverseAndSubtractPed( const AliCaloBunchInfo *bunch, const UInt_t /*altrocfg1*/, const UInt_t /*altrocfg2*/, double *outarray ) const
148 //Time sample comes in reversed order, revers them back
149 //Subtract the baseline based on content of altrocfg1 and altrocfg2.
150 Int_t length = bunch->GetLength();
151 const UShort_t *sig = bunch->GetData();
153 double ped = EvaluatePedestal( sig , length);
155 for( int i=0; i < length; i++ )
157 outarray[i] = sig[length -i -1] - ped;
166 AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * const data, const int /*length*/ ) const
171 if( fIsZerosupressed == false )
173 for(int i=0; i < fNsamplePed; i++ )
179 return tmp/fNsamplePed;
184 AliCaloRawAnalyzer::Max( const AliCaloBunchInfo *const bunch , int *const maxindex ) const
189 const UShort_t *sig = bunch->GetData();
191 for(int i=0; i < bunch->GetLength(); i++ )
193 if( sig[i] > tmpmax )
202 // *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
203 *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
211 AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const
213 // a bunch is considered invalid if the maximum is in the first or last time-bin
216 const UShort_t *sig = bunch->GetData();
218 for(int i=0; i < bunch->GetLength(); i++ )
220 if( sig[i] > tmpmax )
228 if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) )
238 AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,short *const maxampbin, short *const maxamplitude ) const
240 //We select the bunch with the highest amplitude unless any time constraints is set
242 short bunchindex = -1;
246 for(unsigned int i=0; i < bunchvector.size(); i++ )
248 max = Max( &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches
249 if( IsInTimeRange( indx) )
261 if (bunchindex >= 0) {
262 bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) );
273 AliCaloRawAnalyzer::IsInTimeRange( const int maxindex ) const
275 // Ckeck if the index of the max ADC vaue is consistent with trigger.
276 if( ( fMinTimeIndex < 0 && fMaxTimeIndex < 0) ||fMaxTimeIndex < 0 )
280 return ( maxindex < fMaxTimeIndex ) && ( maxindex > fMinTimeIndex ) ? true : false;
285 AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) const
288 cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl;
289 cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl;
291 for(unsigned int i=0; i < bvctr.size() ; i++ )
293 PrintBunch( bvctr.at(i) );
294 cout << " bunch = " << i << endl;
296 cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl;
301 AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) const
304 cout << __FILE__ << ":" << __LINE__ << endl;
305 cout << __FILE__ << __LINE__ << ", startimebin = " << bunch.GetStartBin() << ", length =" << bunch.GetLength() << endl;
306 const UShort_t *sig = bunch.GetData();
308 for ( Int_t j = 0; j < bunch.GetLength(); j++)
310 printf("%d\t", sig[j] );
317 AliCaloRawAnalyzer::CalculateChi2(const Double_t amp, const Double_t time,
318 const Int_t first, const Int_t last,
319 const Double_t adcErr,
323 // amp - max amplitude;
324 // time - time of max amplitude;
325 // first, last - sample array indices to be used
326 // adcErr - nominal error of amplitude measurement (one value for all channels)
327 // if adcErr<0 that mean adcErr=1.
328 // tau - filter time response (in timebin units)
332 if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined..
333 // or, first is negative, the indices are not valid
334 return AliCaloFitResults::kDummy;
337 int nsamples = last - first + 1;
338 // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i \n", first, last, nsamples);
342 Double_t dy = 0.0, xx = 0.0, f=0.0;
344 for (Int_t i=0; i<nsamples; i++) {
345 x = first + i; // timebin
346 xx = (x - time + tau) / tau; // help variable
349 f = amp * xx*xx * TMath::Exp(2 * (1 - xx )) ;
351 dy = fReversed[x] - f;
353 // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy);
356 if (adcErr>0.0) { // weight chi2
357 chi2 /= (adcErr*adcErr);
364 AliCaloRawAnalyzer::CalculateMeanAndRMS(const Int_t first, const Int_t last,
365 Double_t & mean, Double_t & rms)
368 // first, last - sample array indices to be used
370 // mean and RMS of samples
372 // To possibly be used to differentiate good signals from bad before fitting
374 mean = AliCaloFitResults::kDummy;
375 rms = AliCaloFitResults::kDummy;
377 if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined..
378 // or, first is negative, the indices are not valid
382 int nsamples = last - first + 1;
383 // printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples);
386 Double_t sampleSum = 0; // sum of samples
387 Double_t squaredSampleSum = 0; // sum of samples squared
389 for (Int_t i=0; i<nsamples; i++) {
391 sampleSum += fReversed[x];
392 squaredSampleSum += (fReversed[x] * fReversed[x]);
395 mean = sampleSum / nsamples;
396 Double_t squaredMean = squaredSampleSum / nsamples;
397 // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
398 rms = sqrt(squaredMean - mean*mean);
404 AliCaloRawAnalyzer::Evaluate( const vector<AliCaloBunchInfo> &/*bunchvector*/, const UInt_t /*altrocfg1*/, const UInt_t /*altrocfg2*/)
405 { // method to do the selection of what should possibly be fitted
406 // not implemented for base class
407 return AliCaloFitResults( 0, 0 );
412 AliCaloRawAnalyzer::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)
413 { // method to do the selection of what should possibly be fitted
415 short maxampindex = 0;
416 index = SelectBunch( bunchvector, &maxampindex, &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set
419 if( index >= 0 && maxamp >= fAmpCut) // something valid was found, and non-zero amplitude
421 // use more convenient numbering and possibly subtract pedestal
422 ped = ReverseAndSubtractPed( &(bunchvector.at(index)), altrocfg1, altrocfg2, fReversed );
423 maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
425 if ( maxf >= fAmpCut ) // possibly significant signal
427 // select array around max to possibly be used in fit
428 maxrev = maxampindex - bunchvector.at(index).GetStartBin();
429 SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last);
431 // sanity check: maximum should not be in first or last bin
432 // if we should do a fit
433 if (first!=maxrev && last!=maxrev) {
434 // calculate how many samples we have
435 nsamples = last - first + 1;