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

// -*- 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 <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),
  fOverflowCut(950),
  fNsamplePed(3),
  fIsZerosupressed( false ),
  fVerbose( false ),
  fAlgo(Algo::kNONE), 
  fL1Phase(0),
  fAmp(0),
  fTof(0),
  fTau( EMCAL::TAU ),
  fFixTau( true )
{
  // Ctor
  
  snprintf(fName,     256, "%s", name);
  snprintf(fNameShort,256, "%s", nameshort);
  
  for(int i=0; i < ALTROMAXSAMPLES; i++ )
    {
      fReversed[i] = 0;
    }
}


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 > ALTROMAXSAMPLES  || max > ALTROMAXSAMPLES  )
    {
      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
{
  // Get maximum of array
  
  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 *data, int length, short maxindex,
                                    int * first, int * last, 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
  //Until the ADC value is less that fFitArrayCut, or derivative changes sign (data jump)
  int tmpfirst =  maxindex;
  int tmplast  =  maxindex;
  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) && (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 ++;
    }

  // 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,
                                           UInt_t /*altrocfg1*/, 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 * data, int /*length*/ ) const
{
  //  Pedestal evaluation if not zero suppressed
  
  double tmp = 0;

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

  return  tmp/fNsamplePed;
 }


short  
AliCaloRawAnalyzer::Max( const AliCaloBunchInfo * bunch , int * maxindex ) const
{
  // Get maximum in bunch array
  
  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 * maxampbin, short * maxamplitude )
{
  //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, fMaxTimeIndex, fMinTimeIndex) )
	{
	  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 int maxtindx, const int mintindx ) const
{
  // Ckeck if the index of the max ADC vaue is consistent with trigger.
  if( ( mintindx  < 0 && maxtindx   < 0) ||maxtindx  < 0 )
    {
      return true; 
    }
  return ( maxindex < maxtindx ) && ( maxindex > mintindx  ) ? true : false;
}


void 
AliCaloRawAnalyzer::PrintBunches( const vector<AliCaloBunchInfo> &bvctr ) 
{
  // Print bunch vector info
  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 )
{
  // Print bunch infor
  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( Double_t amp,    Double_t time,
                                   Int_t    first,  Int_t    last,
                                   Double_t adcErr, 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(Int_t first, 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,
                                          UInt_t altrocfg1, UInt_t altrocfg2, Int_t & index,
                                          Float_t & maxf, short & maxamp,
                                          short & maxrev, Float_t & ped,
                                          int & first, int & last, 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 >= 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 >= acut  ) // 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, acut );

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