+// -*- mode: c++ -*-
/**************************************************************************
* This file is property of and copyright by *
* the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009 *
* *
- * Primary Author: Per Thomas Hille <p.t.hille@fys.uio.no> *
+ * Primary Author: Per Thomas Hille <perthomas.hille@yale.edu> *
* *
* Contributors are mentioned in the code where appropriate. *
* Please report bugs to p.t.hille@fys.uio.no *
#include <iostream>
using namespace std;
-
-AliCaloRawAnalyzer::AliCaloRawAnalyzer() : TObject(),
- fMinTimeIndex(-1),
- fMaxTimeIndex(-1),
- fFitArrayCut(5),
- fAmpCut(4),
- fNsampleCut(5),
- fIsZerosupressed( false ),
- fVerbose( false )
+ClassImp(AliCaloRawAnalyzer)
+
+AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char *name, const char *nameshort) : TObject(),
+ fMinTimeIndex(-1),
+ fMaxTimeIndex(-1),
+ fFitArrayCut(5),
+ fAmpCut(4),
+ fNsampleCut(5),
+ fOverflowCut(950),
+ fNsamplePed(3),
+ fIsZerosupressed( false ),
+ fVerbose( false ),
+ fAlgo(Algo::kNONE),
+ fL1Phase(0),
+ fAmp(0),
+ fTof(0),
+ fTau( EMCAL::TAU )
{
- //Comment
- for(int i=0; i < MAXSAMPLES; i++ )
+ //Comment
+ snprintf(fName, 256,"%s", name);
+ snprintf(fNameShort,256, "%s", nameshort);
+ // sprintf(fName ,"%s", name);
+ // sprintf(fNameShort, "%s", nameshort);
+
+ for(int i=0; i < ALTROMAXSAMPLES; i++ )
{
fReversed[i] = 0;
}
AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max )
{
//Require that the bin if the maximum ADC value is between min and max (timebin)
- if( ( min > max ) || min > MAXSAMPLES || max > MAXSAMPLES )
+ if( ( min > max ) || min > ALTROMAXSAMPLES || max > ALTROMAXSAMPLES )
{
AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored", min, max ) );
}
void
-AliCaloRawAnalyzer::SelectSubarray( const Double_t *fData, const int length, const short maxindex, int *const first, int *const last ) const
+AliCaloRawAnalyzer::SelectSubarray( const Double_t *data, const int length, const short maxindex,int *const first, int *const last, const int cut) const
{
//Selection of subset of data from one bunch that will be used for fitting or
//Peak finding. Go to the left and right of index of the maximum time bin
- //Untile the ADC value is less that fFitArrayCut
+ //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump)
int tmpfirst = maxindex;
int tmplast = maxindex;
-
- while(( tmpfirst ) > 0 && ( fData[tmpfirst] > fFitArrayCut ))
+ Double_t prevFirst = data[maxindex];
+ Double_t prevLast = data[maxindex];
+ bool firstJump = false;
+ bool lastJump = false;
+
+ while( (tmpfirst >= 0) && (data[tmpfirst] >= cut ) && (!firstJump) )
{
+ // jump check:
+ if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex
+ if ( data[tmpfirst] >= prevFirst) {
+ firstJump = true;
+ }
+ }
+ prevFirst = data[tmpfirst];
tmpfirst -- ;
}
- while(( tmplast ) < length && ( fData [tmplast] > fFitArrayCut ))
+ while( (tmplast < length) && (data[tmplast] >= cut ) && (!lastJump) )
{
+ // jump check:
+ if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex
+ if ( data[tmplast] >= prevLast) {
+ lastJump = true;
+ }
+ }
+ prevLast = data[tmplast];
tmplast ++;
}
-
- *first = tmpfirst +1;
- *last = tmplast -1;
+
+ // we keep one pre- or post- sample if we can (as in online)
+ // check though if we ended up on a 'jump', or out of bounds: if so, back up
+ if (firstJump || tmpfirst<0) tmpfirst ++;
+ if (lastJump || tmplast>=length) tmplast --;
+
+ *first = tmpfirst;
+ *last = tmplast;
+ return;
}
if( fIsZerosupressed == false )
{
- for(int i=0; i < 5; i++ )
+ for(int i=0; i < fNsamplePed; i++ )
{
tmp += data[i];
}
- }
- return tmp/5;
-}
+ }
+
+ return tmp/fNsamplePed;
+ }
short
if(maxindex != 0 )
{
- *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
+ // *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
+ *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
}
return tmpmax;
}
+bool
+AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const
+{
+ // a bunch is considered invalid if the maximum is in the first or last time-bin
+ short tmpmax = -1;
+ int tmpindex = -1;
+ const UShort_t *sig = bunch->GetData();
+
+ for(int i=0; i < bunch->GetLength(); i++ )
+ {
+ if( sig[i] > tmpmax )
+ {
+ tmpmax = sig[i];
+ tmpindex = i;
+ }
+ }
+
+ bool bunchOK = true;
+ if (tmpindex == 0 || tmpindex == (bunch->GetLength()-1) )
+ {
+ bunchOK = false;
+ }
+
+ return bunchOK;
+}
+
+
int
-AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,short *const maxampbin, short *const maxamplitude ) const
+AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector, short *const maxampbin, short *const maxamplitude )
{
//We select the bunch with the highest amplitude unless any time constraints is set
short max = -1;
for(unsigned int i=0; i < bunchvector.size(); i++ )
{
- max = Max( &bunchvector.at(i), &indx );
- if( IsInTimeRange( indx) )
+ max = Max( &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches
+ if( IsInTimeRange( indx, fMaxTimeIndex, fMinTimeIndex) )
{
if( max > maxall )
{
maxall = max;
bunchindex = i;
+ *maxampbin = indx;
+ *maxamplitude = max;
}
}
}
- *maxampbin = indx;
- *maxamplitude = max;
+ if (bunchindex >= 0) {
+ bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) );
+ if (! bunchOK) {
+ bunchindex = -1;
+ }
+ }
+
return bunchindex;
}
bool
-AliCaloRawAnalyzer::IsInTimeRange( const int maxindex ) const
+AliCaloRawAnalyzer::IsInTimeRange( const int maxindex, const int maxtindx, const int mintindx ) const
{
// Ckeck if the index of the max ADC vaue is consistent with trigger.
- if( ( fMinTimeIndex < 0 && fMaxTimeIndex < 0) ||fMaxTimeIndex < 0 )
+ if( ( mintindx < 0 && maxtindx < 0) ||maxtindx < 0 )
{
return true;
}
- return ( maxindex < fMaxTimeIndex ) && ( maxindex > fMinTimeIndex ) ? true : false;
+ return ( maxindex < maxtindx ) && ( maxindex > mintindx ) ? true : false;
}
void
-AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) const
+AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr )
{
//comment
cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl;
void
-AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) const
+AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch )
{
//comment
cout << __FILE__ << ":" << __LINE__ << endl;
}
-AliCaloFitResults
-AliCaloRawAnalyzer::Evaluate( const vector<AliCaloBunchInfo> &/*bunchvector*/, const UInt_t /*altrocfg1*/, const UInt_t /*altrocfg2*/)
-{ // method to do the selection of what should possibly be fitted
- // not implemented for base class
- return AliCaloFitResults( 0, 0, 0, 0, 0, 0, 0 );
+Double_t
+AliCaloRawAnalyzer::CalculateChi2(const Double_t amp, const Double_t time,
+ const Int_t first, const Int_t last,
+ const Double_t adcErr,
+ const Double_t tau) const
+{
+ // Input:
+ // amp - max amplitude;
+ // time - time of max amplitude;
+ // first, last - sample array indices to be used
+ // adcErr - nominal error of amplitude measurement (one value for all channels)
+ // if adcErr<0 that mean adcErr=1.
+ // tau - filter time response (in timebin units)
+ // Output:
+ // chi2 - chi2
+
+ if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined..
+ // or, first is negative, the indices are not valid
+ return Ret::kDummy;
+ }
+
+ int nsamples = last - first + 1;
+ // printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i : amp %3.2f time %3.2f \n", first, last, nsamples, amp, time);
+
+ Int_t x = 0;
+ Double_t chi2 = 0;
+ Double_t dy = 0.0, xx = 0.0, f=0.0;
+
+ for (Int_t i=0; i<nsamples; i++) {
+ x = first + i; // timebin
+ xx = (x - time + tau) / tau; // help variable
+ f = 0;
+ if (xx > 0) {
+ f = amp * xx*xx * TMath::Exp(2 * (1 - xx )) ;
+ }
+ dy = fReversed[x] - f;
+ chi2 += dy*dy;
+ // printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy);
+ }
+
+ if (adcErr>0.0) { // weight chi2
+ chi2 /= (adcErr*adcErr);
+ }
+ return chi2;
}
+void
+AliCaloRawAnalyzer::CalculateMeanAndRMS(const Int_t first, const Int_t last,
+ Double_t & mean, Double_t & rms)
+{
+ // Input:
+ // first, last - sample array indices to be used
+ // Output:
+ // mean and RMS of samples
+ //
+ // To possibly be used to differentiate good signals from bad before fitting
+ //
+ mean = Ret::kDummy;
+ rms = Ret::kDummy;
+
+ if (first == last || first<0 ) { // signal consists of single sample, chi2 estimate (0) not too well defined..
+ // or, first is negative, the indices are not valid
+ return;
+ }
+
+ int nsamples = last - first + 1;
+ // printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples);
+
+ int x = 0;
+ Double_t sampleSum = 0; // sum of samples
+ Double_t squaredSampleSum = 0; // sum of samples squared
+
+ for (Int_t i=0; i<nsamples; i++) {
+ x = first + i;
+ sampleSum += fReversed[x];
+ squaredSampleSum += (fReversed[x] * fReversed[x]);
+ }
+
+ mean = sampleSum / nsamples;
+ Double_t squaredMean = squaredSampleSum / nsamples;
+ // The variance (rms squared) is equal to the mean of the squares minus the square of the mean..
+ rms = sqrt(squaredMean - mean*mean);
+
+ return;
+}
+
+
+
int
-AliCaloRawAnalyzer::PreFitEvaluateSamples( const vector<AliCaloBunchInfo> &bunchvector, const UInt_t altrocfg1, const UInt_t altrocfg2, Int_t & index, Float_t & maxf, short & maxamp, short & maxampindex, Float_t & ped, int & first, int & last)
+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,const int acut )
{ // method to do the selection of what should possibly be fitted
int nsamples = 0;
+ short maxampindex = 0;
index = SelectBunch( bunchvector, &maxampindex, &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set
- if( index >= 0 && maxamp > fAmpCut) // something valid was found, and non-zero amplitude
+ if( index >= 0 && maxamp >= acut ) // something valid was found, and non-zero amplitude
{
// use more convenient numbering and possibly subtract pedestal
ped = ReverseAndSubtractPed( &(bunchvector.at(index)), altrocfg1, altrocfg2, fReversed );
maxf = TMath::MaxElement( bunchvector.at(index).GetLength(), fReversed );
- if ( maxf > fAmpCut ) // possibly significant signal
+ if ( maxf >= acut ) // possibly significant signal
{
// select array around max to possibly be used in fit
- maxampindex -= bunchvector.at(index).GetStartBin(); // PTH - why isn't this index also reversed for call below?
- SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxampindex , &first, &last);
+ maxrev = maxampindex - bunchvector.at(index).GetStartBin();
+ SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last, acut );
// sanity check: maximum should not be in first or last bin
// if we should do a fit
- if (first!=maxampindex && last!=maxampindex) {
- // calculate how many samples we have
- nsamples = last - first + 1;
- }
+ if (first!=maxrev && last!=maxrev)
+ {
+ // calculate how many samples we have
+ nsamples = last - first + 1;
+ }
}
}