[u/mrichter/AliRoot.git] / EMCAL / AliCaloRawAnalyzer.cxx

/**************************************************************************
 * This file is property of and copyright by                              *
 * the Relativistic Heavy Ion Group (RHIG), Yale University, US, 2009     *
 *                                                                        *
 * 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                             *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/


// Base class for extraction 
// of signal amplitude and peak position
// From CALO Calorimeter RAW data (from the RCU)
// Contains some utilities for preparing / selecting
// Signals suitable for signal extraction
// By derived classes

#include "AliLog.h"
#include "AliCaloRawAnalyzer.h"
#include "AliCaloBunchInfo.h"
#include "AliCaloFitResults.h"
#include "TMath.h"
#include <iostream>
using namespace std;

ClassImp(AliCaloRawAnalyzer)  

AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char *name, const char *nameshort) :  TObject(),
										   fMinTimeIndex(-1),
										   fMaxTimeIndex(-1),
										   fFitArrayCut(5),
										   fAmpCut(4),
										   fNsampleCut(5),
										   fNsamplePed(3),
										   fIsZerosupressed( false ),
										   fVerbose( false )
{
  //Comment 
  sprintf(fName, "%s", name);
  sprintf(fNameShort, "%s", nameshort);
    
  for(int i=0; i < MAXSAMPLES; i++ )
    {
      fReversed[i] = 0;
    }
}

AliCaloRawAnalyzer::~AliCaloRawAnalyzer()
{

}


void 
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  )
    {
      AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored",  min, max ) ); 
    }
  else
    {
      fMinTimeIndex  = min;
      fMaxTimeIndex  = max;
    }
}


UShort_t 
AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) const
{
  //------------
  UShort_t tmpmax  = data[0];

  for(int i=0; i < length; i++)
    {
      if( tmpmax  <  data[i] )
	{
	  tmpmax = data[i];
	}
    }
  return tmpmax;
}


void 
AliCaloRawAnalyzer::SelectSubarray( const Double_t *fData, const int length, const short maxindex, int *const first,  int *const last ) 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
  //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump)
  int tmpfirst =  maxindex;
  int tmplast  =  maxindex;
  Double_t prevFirst =  fData[maxindex];
  Double_t prevLast  =  fData[maxindex];  
  bool firstJump = false;
  bool lastJump = false;

  while( (tmpfirst >= 0) && (fData[tmpfirst] >= fFitArrayCut) && (!firstJump) ) 
    {
      // jump check:
      if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex
	if (fData[tmpfirst] >= prevFirst) {
	  firstJump = true;
	}
      }
      prevFirst = fData[tmpfirst];
      tmpfirst -- ;
    }
  
  while( (tmplast < length) && (fData[tmplast] >= fFitArrayCut) && (!lastJump) ) 
    {
      // jump check:
      if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex
	if (fData[tmplast] >= prevLast) {
	  lastJump = true;
	}
      }
      prevLast = fData[tmplast];
      tmplast ++;
    }

  // 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;
}


Float_t 
AliCaloRawAnalyzer::ReverseAndSubtractPed( const AliCaloBunchInfo *bunch, const UInt_t /*altrocfg1*/,  const UInt_t /*altrocfg2*/, double *outarray ) const
{
  //Time sample comes in reversed order, revers them back
  //Subtract the baseline based on content of altrocfg1 and altrocfg2.
  Int_t length =  bunch->GetLength(); 
  const UShort_t *sig = bunch->GetData();

  double ped = EvaluatePedestal( sig , length);

  for( int i=0; i < length; i++ )
    {
      outarray[i] = sig[length -i -1] - ped;
    }

  return ped;
}


Float_t 
AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * const data, const int /*length*/ ) const
{
  //  double ped = 0;
  double tmp = 0;

  if( fIsZerosupressed == false ) 
    {
      for(int i=0; i < fNsamplePed; i++ )
	{
	  tmp += data[i];
	}
   }

  return  tmp/fNsamplePed;
 }


short  
AliCaloRawAnalyzer::Max( const AliCaloBunchInfo *const bunch , int *const maxindex ) const
{
  //comment
  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; 
	}
    }
  
  if(maxindex != 0 )
    {
      //   *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
{
  //We select the bunch with the highest amplitude unless any time constraints is set
  short max = -1;
  short bunchindex = -1;
  short maxall = -1;
  int indx = -1;

  for(unsigned int i=0; i < bunchvector.size(); i++ )
    {
      max = Max(  &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches  
      if( IsInTimeRange( indx) )
	{
	  if( max > maxall )
	    {
	      maxall = max;
	      bunchindex = i;
	      *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
{
  // Ckeck if the index of the max ADC vaue is consistent with trigger.
  if( ( fMinTimeIndex  < 0 && fMaxTimeIndex  < 0) ||fMaxTimeIndex  < 0 )
    {
      return true; 
    }
  return ( maxindex < fMaxTimeIndex ) && ( maxindex > fMinTimeIndex  ) ? true : false;
}


void 
AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) const
{
  //comment
  cout << __FILE__ << __LINE__<< "*************** Printing Bunches *******************" << endl;
  cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl;

  for(unsigned int i=0; i < bvctr.size() ; i++ )
    {
      PrintBunch( bvctr.at(i) );
      cout << " bunch = "  <<  i  << endl;
    }
  cout << __FILE__ << __LINE__<< "*************** Done ... *******************" << endl;
}


void 
AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) const
{
  //comment
  cout << __FILE__ << ":" << __LINE__ << endl;
  cout << __FILE__ << __LINE__   << ", startimebin = " << bunch.GetStartBin() << ", length =" <<  bunch.GetLength() << endl;
  const UShort_t *sig =  bunch.GetData();  
  
  for ( Int_t j = 0; j <  bunch.GetLength();  j++) 
    {
      printf("%d\t", sig[j] );
    }
  cout << endl; 
}


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)
{
  //   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 AliCaloFitResults::kDummy;
  }

  int nsamples =  last - first + 1;
  //  printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i \n", first, last, nsamples); 

  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 = AliCaloFitResults::kDummy;
  rms = AliCaloFitResults::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;
}

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 );
}


int
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)
{ // 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
    {
      // 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
	{
	  // select array around max to possibly be used in fit
	  maxrev = maxampindex - bunchvector.at(index).GetStartBin(); 
	  SelectSubarray( fReversed,  bunchvector.at(index).GetLength(),  maxrev, &first, &last);

	  // sanity check: maximum should not be in first or last bin
	  // if we should do a fit
	  if (first!=maxrev && last!=maxrev) {
	    // calculate how many samples we have 
	    nsamples =  last - first + 1;	  
	  }
	}
    }

  return nsamples;
}
Commit	Line	Data
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>
	33	using namespace std;
	34
ce95bae9	35	ClassImp(AliCaloRawAnalyzer)
48a2e3eb	36
	37	AliCaloRawAnalyzer::AliCaloRawAnalyzer(const char name, const char nameshort) : TObject(),
	38	fMinTimeIndex(-1),
	39	fMaxTimeIndex(-1),
	40	fFitArrayCut(5),
	41	fAmpCut(4),
	42	fNsampleCut(5),
f57baa2d	43	fNsamplePed(3),
48a2e3eb	44	fIsZerosupressed( false ),
48a2e3eb	45	fVerbose( false )
d655d7dd	46	{
e37e3c84	47	//Comment
e37e3c84	48	sprintf(fName, "%s", name);
48a2e3eb	49	sprintf(fNameShort, "%s", nameshort);
48a2e3eb	50
d655d7dd	51	for(int i=0; i < MAXSAMPLES; i++ )
	52	{
	53	fReversed[i] = 0;
	54	}
	55	}
	56
57839add	57	AliCaloRawAnalyzer::~AliCaloRawAnalyzer()
d655d7dd	58	{
	59
	60	}
	61
	62
	63	void
57839add	64	AliCaloRawAnalyzer::SetTimeConstraint(const int min, const int max )
d655d7dd	65	{
	66	//Require that the bin if the maximum ADC value is between min and max (timebin)
	67	if( ( min > max ) \|\| min > MAXSAMPLES \|\| max > MAXSAMPLES )
	68	{
	69	AliWarning( Form( "Attempt to set Invalid time bin range (Min , Max) = (%d, %d), Ingored", min, max ) );
	70	}
	71	else
	72	{
	73	fMinTimeIndex = min;
	74	fMaxTimeIndex = max;
	75	}
	76	}
	77
	78
	79	UShort_t
57839add	80	AliCaloRawAnalyzer::Max(const UShort_t *data, const int length ) const
d655d7dd	81	{
	82	//------------
	83	UShort_t tmpmax = data[0];
	84
	85	for(int i=0; i < length; i++)
	86	{
	87	if( tmpmax < data[i] )
	88	{
	89	tmpmax = data[i];
	90	}
	91	}
	92	return tmpmax;
	93	}
	94
	95
	96	void
57839add	97	AliCaloRawAnalyzer::SelectSubarray( const Double_t fData, const int length, const short maxindex, int const first, int *const last ) const
d655d7dd	98	{
	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
e7acbc37	101	//Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump)
d655d7dd	102	int tmpfirst = maxindex;
d655d7dd	103	int tmplast = maxindex;
947c8e28	104	Double_t prevFirst = fData[maxindex];
	105	Double_t prevLast = fData[maxindex];
	106	bool firstJump = false;
	107	bool lastJump = false;
	108
	109	while( (tmpfirst >= 0) && (fData[tmpfirst] >= fFitArrayCut) && (!firstJump) )
d655d7dd	110	{
947c8e28	111	// jump check:
	112	if (tmpfirst != maxindex) { // neighbor to maxindex can share peak with maxindex
	113	if (fData[tmpfirst] >= prevFirst) {
	114	firstJump = true;
	115	}
	116	}
	117	prevFirst = fData[tmpfirst];
d655d7dd	118	tmpfirst -- ;
	119	}
	120
947c8e28	121	while( (tmplast < length) && (fData[tmplast] >= fFitArrayCut) && (!lastJump) )
d655d7dd	122	{
947c8e28	123	// jump check:
	124	if (tmplast != maxindex) { // neighbor to maxindex can share peak with maxindex
	125	if (fData[tmplast] >= prevLast) {
	126	lastJump = true;
	127	}
	128	}
	129	prevLast = fData[tmplast];
d655d7dd	130	tmplast ++;
d655d7dd	131	}
f57baa2d	132
947c8e28	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 --;
	137
	138	*first = tmpfirst;
	139	*last = tmplast;
e7acbc37	140	return;
d655d7dd	141	}
	142
	143
	144
	145	Float_t
57839add	146	AliCaloRawAnalyzer::ReverseAndSubtractPed( const AliCaloBunchInfo bunch, const UInt_t /altrocfg1/, const UInt_t /altrocfg2/, double outarray ) const
d655d7dd	147	{
	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();
	152
57839add	153	double ped = EvaluatePedestal( sig , length);
	154
	155	for( int i=0; i < length; i++ )
	156	{
	157	outarray[i] = sig[length -i -1] - ped;
	158	}
	159
	160	return ped;
	161	}
	162
	163
	164
	165	Float_t
	166	AliCaloRawAnalyzer::EvaluatePedestal(const UShort_t * const data, const int /length/ ) const
	167	{
	168	// double ped = 0;
d655d7dd	169	double tmp = 0;
	170
	171	if( fIsZerosupressed == false )
	172	{
f57baa2d	173	for(int i=0; i < fNsamplePed; i++ )
d655d7dd	174	{
57839add	175	tmp += data[i];
d655d7dd	176	}
f57baa2d	177	}
e37e3c84	178
f57baa2d	179	return tmp/fNsamplePed;
f57baa2d	180	}
d655d7dd	181
	182
	183	short
57839add	184	AliCaloRawAnalyzer::Max( const AliCaloBunchInfo const bunch , int const maxindex ) const
d655d7dd	185	{
	186	//comment
	187	short tmpmax = -1;
	188	int tmpindex = -1;
	189	const UShort_t *sig = bunch->GetData();
	190
	191	for(int i=0; i < bunch->GetLength(); i++ )
	192	{
	193	if( sig[i] > tmpmax )
	194	{
	195	tmpmax = sig[i];
	196	tmpindex = i;
	197	}
	198	}
	199
	200	if(maxindex != 0 )
	201	{
48a2e3eb	202	// *maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
48a2e3eb	203	*maxindex = bunch->GetLength() -1 - tmpindex + bunch->GetStartBin();
d655d7dd	204	}
	205
	206	return tmpmax;
	207	}
	208
	209
507751ce	210	bool
	211	AliCaloRawAnalyzer::CheckBunchEdgesForMax( const AliCaloBunchInfo *const bunch ) const
	212	{
	213	// a bunch is considered invalid if the maximum is in the first or last time-bin
	214	short tmpmax = -1;
	215	int tmpindex = -1;
	216	const UShort_t *sig = bunch->GetData();
	217
	218	for(int i=0; i < bunch->GetLength(); i++ )
	219	{
	220	if( sig[i] > tmpmax )
	221	{
	222	tmpmax = sig[i];
	223	tmpindex = i;
	224	}
	225	}
	226
	227	bool bunchOK = true;
	228	if (tmpindex == 0 \|\| tmpindex == (bunch->GetLength()-1) )
	229	{
	230	bunchOK = false;
	231	}
	232
	233	return bunchOK;
	234	}
	235
	236
d655d7dd	237	int
57839add	238	AliCaloRawAnalyzer::SelectBunch( const vector<AliCaloBunchInfo> &bunchvector,short const maxampbin, short const maxamplitude ) const
d655d7dd	239	{
	240	//We select the bunch with the highest amplitude unless any time constraints is set
	241	short max = -1;
	242	short bunchindex = -1;
	243	short maxall = -1;
	244	int indx = -1;
	245
	246	for(unsigned int i=0; i < bunchvector.size(); i++ )
	247	{
48a2e3eb	248	max = Max( &bunchvector.at(i), &indx ); // CRAP PTH, bug fix, trouble if more than one bunches
d655d7dd	249	if( IsInTimeRange( indx) )
	250	{
	251	if( max > maxall )
	252	{
	253	maxall = max;
	254	bunchindex = i;
48a2e3eb	255	*maxampbin = indx;
48a2e3eb	256	*maxamplitude = max;
d655d7dd	257	}
	258	}
	259	}
	260
507751ce	261	if (bunchindex >= 0) {
	262	bool bunchOK = CheckBunchEdgesForMax( &bunchvector.at(bunchindex) );
	263	if (! bunchOK) {
	264	bunchindex = -1;
	265	}
	266	}
	267
d655d7dd	268	return bunchindex;
	269	}
	270
	271
	272	bool
57839add	273	AliCaloRawAnalyzer::IsInTimeRange( const int maxindex ) const
d655d7dd	274	{
	275	// Ckeck if the index of the max ADC vaue is consistent with trigger.
	276	if( ( fMinTimeIndex < 0 && fMaxTimeIndex < 0) \|\|fMaxTimeIndex < 0 )
	277	{
	278	return true;
	279	}
	280	return ( maxindex < fMaxTimeIndex ) && ( maxindex > fMinTimeIndex ) ? true : false;
	281	}
	282
	283
	284	void
57839add	285	AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) const
d655d7dd	286	{
	287	//comment
	288	cout << __FILE__ << __LINE__<< "************* Printing Bunches *****************" << endl;
	289	cout << __FILE__ << __LINE__<< "*** There are " << bvctr.size() << ", bunches" << endl;
	290
	291	for(unsigned int i=0; i < bvctr.size() ; i++ )
	292	{
	293	PrintBunch( bvctr.at(i) );
	294	cout << " bunch = " << i << endl;
	295	}
	296	cout << __FILE__ << __LINE__<< "************* Done ... *****************" << endl;
	297	}
	298
	299
	300	void
57839add	301	AliCaloRawAnalyzer::PrintBunch( const AliCaloBunchInfo &bunch ) const
d655d7dd	302	{
	303	//comment
	304	cout << __FILE__ << ":" << __LINE__ << endl;
	305	cout << __FILE__ << __LINE__ << ", startimebin = " << bunch.GetStartBin() << ", length =" << bunch.GetLength() << endl;
	306	const UShort_t *sig = bunch.GetData();
	307
	308	for ( Int_t j = 0; j < bunch.GetLength(); j++)
	309	{
	310	printf("%d\t", sig[j] );
	311	}
	312	cout << endl;
	313	}
	314
48a2e3eb	315
947c8e28	316	Double_t
	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,
	320	const Double_t tau)
	321	{
	322	// Input:
	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)
	329	// Output:
	330	// chi2 - chi2
	331
	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;
	335	}
	336
	337	int nsamples = last - first + 1;
	338	// printf(" AliCaloRawAnalyzer::CalculateChi2 : first %i last %i : nsamples %i \n", first, last, nsamples);
	339
	340	Int_t x = 0;
	341	Double_t chi2 = 0;
	342	Double_t dy = 0.0, xx = 0.0, f=0.0;
	343
	344	for (Int_t i=0; i<nsamples; i++) {
	345	x = first + i; // timebin
	346	xx = (x - time + tau) / tau; // help variable
	347	f = 0;
	348	if (xx > 0) {
	349	f = amp * xxxx TMath::Exp(2 * (1 - xx )) ;
	350	}
	351	dy = fReversed[x] - f;
	352	chi2 += dy*dy;
	353	// printf(" AliCaloRawAnalyzer::CalculateChi2 : %i : y %f -> f %f : dy %f \n", i, fReversed[first+i], f, dy);
	354	}
	355
	356	if (adcErr>0.0) { // weight chi2
	357	chi2 /= (adcErr*adcErr);
	358	}
	359	return chi2;
	360	}
	361
	362
	363	void
	364	AliCaloRawAnalyzer::CalculateMeanAndRMS(const Int_t first, const Int_t last,
	365	Double_t & mean, Double_t & rms)
	366	{
	367	// Input:
	368	// first, last - sample array indices to be used
	369	// Output:
	370	// mean and RMS of samples
	371	//
	372	// To possibly be used to differentiate good signals from bad before fitting
	373	//
	374	mean = AliCaloFitResults::kDummy;
	375	rms = AliCaloFitResults::kDummy;
	376
	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
	379	return;
380	}
381
382	int nsamples = last - first + 1;
383	// printf(" AliCaloRawAnalyzer::CalculateMeanAndRMS : first %i last %i : nsamples %i \n", first, last, nsamples);
384
385	int x = 0;
386	Double_t sampleSum = 0; // sum of samples
387	Double_t squaredSampleSum = 0; // sum of samples squared
388
389	for (Int_t i=0; i<nsamples; i++) {
390	x = first + i;
391	sampleSum += fReversed[x];
392	squaredSampleSum += (fReversed[x] * fReversed[x]);
393	}
394
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);
399
400	return;
401	}
402
57839add	403	AliCaloFitResults
	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
507751ce	407	return AliCaloFitResults( 0, 0 );
57839add	408	}
57839add	409
48a2e3eb	410
57839add	411	int
f57baa2d	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)
57839add	413	{ // method to do the selection of what should possibly be fitted
57839add	414	int nsamples = 0;
f57baa2d	415	short maxampindex = 0;
57839add	416	index = SelectBunch( bunchvector, &maxampindex, &maxamp ); // select the bunch with the highest amplitude unless any time constraints is set
	417
	418
1afbf4b1	419	if( index >= 0 && maxamp >= fAmpCut) // something valid was found, and non-zero amplitude
57839add	420	{
	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 );
	424
1afbf4b1	425	if ( maxf >= fAmpCut ) // possibly significant signal
57839add	426	{
57839add	427	// select array around max to possibly be used in fit
f57baa2d	428	maxrev = maxampindex - bunchvector.at(index).GetStartBin();
f57baa2d	429	SelectSubarray( fReversed, bunchvector.at(index).GetLength(), maxrev, &first, &last);
57839add	430
	431	// sanity check: maximum should not be in first or last bin
	432	// if we should do a fit
f57baa2d	433	if (first!=maxrev && last!=maxrev) {
57839add	434	// calculate how many samples we have
	435	nsamples = last - first + 1;
	436	}
	437	}
	438	}
	439
	440	return nsamples;
	441	}
	442