[u/mrichter/AliRoot.git] / HLT / PHOS / AliPHOSFitter.h

#ifndef ALIPHOSFITTER_H
#define ALIPHOSFITTER_H
/* Copyright(c) 2006, ALICE Experiment at CERN, All rights reserved. *
 * See cxx source for full Copyright notice                               */


/**
 * The PhosFitter class is the class for extracting the basic signal parameters
 * "timing" and "energy" from the PHOS raw data. Physical data will for a given readout channel be
 * a sequense of ADC digitized 10 bit integer values, however for performance reasons all values used in
 * calculation is of type double.
 **/
class AliPHOSFitter
      {
      public:
	/**
	 * Default constructor
	 **/
	AliPHOSFitter();

	/**
	 * Main constructor
	 * @param dataPtr Data array for wich a subarray will be taken to perform the fit
	 * @param fs the sampling frequency in entities of MHz. Needed in order to calculate physical time
	 **/
	AliPHOSFitter(double *dataPtr, double fs);

	/**
	 * Destructor
	 **/
	~AliPHOSFitter();


	/**
	 * Attemps to level the basline to zero.
	 * The baseline will be calculated from the pretrigger samples and subtracted from the 
	 * data array. 
	 * If pretrigger samples are not present then the basline correction will be incorrect. 
	 * @param dataPtr array for wich to correct the basline 
	 * @param N the number of pretrigger samples used to calculate the baseline.
	 **/
	void BaselineCorrection(double *dataPtr, int N);

	/**
	 * Shifts the basline with the amount given by baselineValue
	 * If pretrigger samples are not present then the basline correction will be incorrect. 
	 * @param dataPtr array for wich to correct the basline 
	 * @param BaslineValue the basline value to subtract..
	 **/
	void BaselineCorrection(double *dataPtr, double baselineValue);

	/**
	 * Extraction of timing and energy using a least square fit assuming uncorrelated measurement errors.
	 * This is also called a "Chi square fit". The default method is the Levenberg Marquardt nonlinear fit. 
	 * If the errors are correlated (which is typically the case)
	 * the timing and energy will not be mimimum variance lower bound. Correlated errors might also give a
	 * systematic error. The parameters "start" and "length" defines the subarray of the data array set in the
	 * constructor that will be used in the fit. "start + length" cannot exeed the length of the data array.
	 * The baseline must be subtracted before performing the fit othervise the result can be biased.
	 * the initial guess parameters is found by the method "MakeInitialGuess".	
	 * @param start the start index of the subarray of the data array. 
	 * @param length the number of samples to use starting from index 
 	 * @param start the start index of the subarray of the data array. 
	 * @param length the number of samples to use starting from index 
	 **/
	void FitChiSquare(int start, int lenght);

	/** 
	 * Extraction of timing and energy using a least square fit assuming uncorrelated measurement errors.
	 * This is also called a "Chi square fit". The default method is the Levenberg Marquardt nonlinear fit. 
	 * If the errors are correlated (which is typically the case)
	 * the timing and energy will not be mimimum variance lower bound. Correlated errors might also give a
	 * systematic error. The parameters "start" and "length" defines the subarray of the data array set in the
	 * constructor that will be used in the fit. "start + length" cannot exeed the length of the data array.
	 * The baseline must be subtracted before performing the fit othervise the result can be biased.
	 * A good initial guess greatly enchanes performance. A poor initial guess might give a solution that is a local 
	 * minima instead of a global minima. If the startindex of the pulse is the first sample above a low (2-3 ADC levels)
	 * above the baseline, then the initial guess of t0 an be set to zero. A good initial guess for amplitude will
	 * for a gamma 2 function be for instance (Max sample valu - baseline)*exp(2). If a Gamma N function is used
	 * multiply with exp(N).
 	 * @param start the start index of the subarray of the data array. 
	 * @param length the number of samples to use starting from index 
	 * @param tGuess initial guess for  timing (in entities of samples)
	 * @param aGuess initial guess for energy in entities of (ADC channels)*exp(2) (for a gamma 2 fuction)
	 **/
	void FitChiSquare(int start, int lenght, double tGuess, double aGuess);

	/**
	 * Extraction of timing an energy using the
	 * K-level method.
	 * @param kLevel the K-level
	 * "start" and "length" defines the subarray of the data array (set in the constructor)
	 * @param start index of samples to use relative to the the data array
	 * @param length the number of samples to use (starting at "start")
	 **/
	void FitKLevel(double kLevel, int start, int lenght);
	
	/** 
	 * This method gives the statistically best possible unbiased estimate of t0 and amplitude provided that a good
	 * estimate of the covariance matrix exists.
	 * Extraction of timing and energy using a least square fit assuming correlated errors.
	 * The user must supply the autocovariance matrix. For N number of samples this matrix
	 * will be a NxN matrix. 
	 * The default method is the Levenberg Marquardt nonlinear fit. 
	 * The parameters "start" and "length" defines the subarray of the data array set in the
	 * constructor that will be used in the fit. "start + length" cannot exeed the length of the data array.
	 * The baseline must be subtracted before performing the fit othervise the result wil biased.
	 * If the startindex of the pulse is the first sample above a low (2-3 ADC levels)
	 * above the baseline, then the initial guess of t0 an be set to zero. A good initial guess for amplitude will
	 * for a gamma 2 function be for instance (Max sample valu - baseline)*exp(2). If a Gamma N function is used
	 * multiply with exp(N).
	 * The correlation matrix for the parameters is written to the matrix pointe to by pCovar 
	 * @param start the start index of the subarray of the data array. 
	 * @param length the number of samples to use starting from "start" 
	 * @param kmCovarPtrPtr the measurement correlation matrix 	
	 * @param pCovarPtrPtr the correlation matrix of the estimated parameters
	 * @param tGuess initial guess of t0 in entities of sampling intervals
	 * @param aGuess initial guess of the amplitude in entities of ADC levels*exp(2) 
	**/
	void FitLeastMeanSquare(int start, int lenght, const double **kmCovarPtrPtr, double **pCovarPtrPtr);

	/** 
	 * This method gives the statistically best possible unbiased estimate of t0 and amplitude provided that a good
	 * estimate of the covariance matrix exists.
	 * For N number of samples this covariance matrix must be an NxN matrix.
	 * The default method is the Levenberg Marquardt nonlinear fit. 
	 * The parameters "start" and "length" defines the subarray of the data array set in the
	 * constructor that will be used in the fit. "start + length" cannot exeed the length of the data array.
	 * The baseline must be subtracted before performing the fit othervise the result wil biased.
	 * A good initial guess will greatly enchance performance.
	 * a poor initial guess will slow down the algorithm and possibly lead to convergens to local minima
	 * of the cost function instead of the global minima. A very bad intial guess might give divergence of the solution.
	 * If the startindex of the pulse is the first sample above a low (2-3 ADC levels)
	 * above the baseline, then the initial guess of t0 can be set to zero. A good initial guess for amplitude will
	 * for a gamma 2 function be for instance (Max sample value - baseline)*exp(2). If a Gamma N function is used
	 * multiply with exp(N).
	 * @param start the start index of the subarray of the data array. 
	 * @param length the number of samples to use starting from index 
	 * @param mCovar the measurement correlation matrix 	
	 * @param pCovar the correlation matrix of the estimated parameters
	 * @param tGuess initial guess of t0 in entities of sampling intervals
	 * @param aGuess initial guess of the amplitude in entities of ADC levels*exp(2)
	**/
	void FitLeastMeanSquare(int start, int lenght, const double **kmCovarPtrPtr, double **pCovarPtrPtr, double tGuess, double aGuess);

	/**
	 * Extraction of timing and energy using the Peakfinde Algorithm.
	 * The. The parameters "start" and "length" defines a sub array  of the data array
	 * that will be used for the the fit. If start+length must not exeed the total length
	 * of the Data array. "start" must be chosen as close as possible to t0.
	 * The baseline must also be subtracted.
	 * The length of "tVector" and "aVector" mus be equal to length.
	 * "index + length" must not exeed the length of the data array set in the constructor.
	 * @param start the start index of the subarray of the data array. 
	 * @param length the number of samples to use starting from index 
	 * @param tVector the peakfinder vector for timing
	 * @param aVector the peakfinder vector for amplitude (energy)
	 **/
	void FitPeakFinder(int start, int lenght, double *tVector, double *aVector);

	/**
	 * This method finds the start index of the pulse (relative to the data array) by serching for
	 * three or more continious samples above trheshold.
	 * @param treshold treshold to use when searchin for the start of the pulse
	**/
	int FindStartIndex(double treshold);
	
	/**
	 * Gives the timing in entities of sample indexes
	 * Physical time is found by multiplying  with the sampling intervall (Ts).
	 **/	
	float GetTiming();

	/**
	 * Gives the time in entities of ADC channels (quantization levels).  
	 * Absolute enrgy is found by multiplying with offline calibration constants.
	**/
	float GetEnergy();

	//	/**
	//	 * Set data array. Overrrides the data array set in the constructor.
	//	 **/
	//	void SetData(int *data);

	/**
	 * Set data array. Overrides data data array set in the constructor.
	 **/
	void SetData(double *data);


      private:

	/**
	 * This function applies only to the Chi and Least mean square fit. An initial guess is made
	 * based on the average of the first 5 samples and the first value exeeding this value.
	 **/
	void MakeInitialGuess();

	/**
	 * This function applies only to the Chi and Least mean square fit. An initial guess is made
	 * based on the average of the first 5 samples and the first value exeeding threshold + this value.
	 * @param treshold The index of the first value above treshold is ntaken to be the first value.
	 **/
	void MakeInitialGuess(int treshold);
	
	double    *fFloatDataPtr;    /**<Float representation of data that should be fitted */
	double     fSampleFrequency; /**<The ADC sample frequency in MHz used under data taking */
	double     fTau;	     /**<The risetime in micro seconds*/		 
	double     fDTof;            /**<Time of flight in entities of sample intervals */
	double     fDAmpl;           /**<Amplitude in entities of ADC levels*/
	double     fDTofGuess;       /**<Initial guess for t0*/
	double     fDAmplGuess;      /**<Initial guess for amplitude*/
	double   **kfMCovarPtrPtr;   /**<Covariance matrix of the measurements*/
	double   **fPCovarPtrPtr;    /**<Covariance matrix of the estimated parameters*/
      };

#endif
Commit	Line	Data
9b37718a	1	#ifndef ALIPHOSFITTER_H
	2	#define ALIPHOSFITTER_H
	3	/* Copyright(c) 2006, ALICE Experiment at CERN, All rights reserved. *
	4	* See cxx source for full Copyright notice */
	5
	6
	7
	8	/**
	9	* The PhosFitter class is the class for extracting the basic signal parameters
	10	* "timing" and "energy" from the PHOS raw data. Physical data will for a given readout channel be
	11	* a sequense of ADC digitized 10 bit integer values, however for performance reasons all values used in
	12	* calculation is of type double.
	13	**/
	14	class AliPHOSFitter
	15	{
	16	public:
	17	/**
	18	* Default constructor
	19	**/
	20	AliPHOSFitter();
	21
	22	/**
	23	* Main constructor
	24	* @param dataPtr Data array for wich a subarray will be taken to perform the fit
	25	* @param fs the sampling frequency in entities of MHz. Needed in order to calculate physical time
	26	**/
	27	AliPHOSFitter(double *dataPtr, double fs);
	28
	29	/**
	30	* Destructor
	31	**/
	32	~AliPHOSFitter();
	33
	34
	35	/**
	36	* Attemps to level the basline to zero.
	37	* The baseline will be calculated from the pretrigger samples and subtracted from the
	38	* data array.
	39	* If pretrigger samples are not present then the basline correction will be incorrect.
	40	* @param dataPtr array for wich to correct the basline
	41	* @param N the number of pretrigger samples used to calculate the baseline.
	42	**/
	43	void BaselineCorrection(double *dataPtr, int N);
	44
	45	/**
	46	* Shifts the basline with the amount given by baselineValue
	47	* If pretrigger samples are not present then the basline correction will be incorrect.
	48	* @param dataPtr array for wich to correct the basline
	49	* @param BaslineValue the basline value to subtract..
	50	**/
	51	void BaselineCorrection(double *dataPtr, double baselineValue);
	52
	53	/**
	54	* Extraction of timing and energy using a least square fit assuming uncorrelated measurement errors.
	55	* This is also called a "Chi square fit". The default method is the Levenberg Marquardt nonlinear fit.
	56	* If the errors are correlated (which is typically the case)
	57	* the timing and energy will not be mimimum variance lower bound. Correlated errors might also give a
	58	* systematic error. The parameters "start" and "length" defines the subarray of the data array set in the
	59	* constructor that will be used in the fit. "start + length" cannot exeed the length of the data array.
	60	* The baseline must be subtracted before performing the fit othervise the result can be biased.
	61	* the initial guess parameters is found by the method "MakeInitialGuess".
	62	* @param start the start index of the subarray of the data array.
	63	* @param length the number of samples to use starting from index
	64	* @param start the start index of the subarray of the data array.
65	* @param length the number of samples to use starting from index
66	**/
67	void FitChiSquare(int start, int lenght);
68
69	/**
70	* Extraction of timing and energy using a least square fit assuming uncorrelated measurement errors.
71	* This is also called a "Chi square fit". The default method is the Levenberg Marquardt nonlinear fit.
72	* If the errors are correlated (which is typically the case)
73	* the timing and energy will not be mimimum variance lower bound. Correlated errors might also give a
74	* systematic error. The parameters "start" and "length" defines the subarray of the data array set in the
75	* constructor that will be used in the fit. "start + length" cannot exeed the length of the data array.
76	* The baseline must be subtracted before performing the fit othervise the result can be biased.
77	* A good initial guess greatly enchanes performance. A poor initial guess might give a solution that is a local
78	* minima instead of a global minima. If the startindex of the pulse is the first sample above a low (2-3 ADC levels)
79	* above the baseline, then the initial guess of t0 an be set to zero. A good initial guess for amplitude will
80	* for a gamma 2 function be for instance (Max sample valu - baseline)*exp(2). If a Gamma N function is used
81	* multiply with exp(N).
82	* @param start the start index of the subarray of the data array.
83	* @param length the number of samples to use starting from index
84	* @param tGuess initial guess for timing (in entities of samples)
85	* @param aGuess initial guess for energy in entities of (ADC channels)*exp(2) (for a gamma 2 fuction)
86	**/
87	void FitChiSquare(int start, int lenght, double tGuess, double aGuess);
88
89	/**
90	* Extraction of timing an energy using the
91	* K-level method.
92	* @param kLevel the K-level
93	* "start" and "length" defines the subarray of the data array (set in the constructor)
94	* @param start index of samples to use relative to the the data array
95	* @param length the number of samples to use (starting at "start")
96	**/
97	void FitKLevel(double kLevel, int start, int lenght);
98
99	/**
100	* This method gives the statistically best possible unbiased estimate of t0 and amplitude provided that a good
101	* estimate of the covariance matrix exists.
102	* Extraction of timing and energy using a least square fit assuming correlated errors.
103	* The user must supply the autocovariance matrix. For N number of samples this matrix
104	* will be a NxN matrix.
105	* The default method is the Levenberg Marquardt nonlinear fit.
106	* The parameters "start" and "length" defines the subarray of the data array set in the
107	* constructor that will be used in the fit. "start + length" cannot exeed the length of the data array.
108	* The baseline must be subtracted before performing the fit othervise the result wil biased.
109	* If the startindex of the pulse is the first sample above a low (2-3 ADC levels)
110	* above the baseline, then the initial guess of t0 an be set to zero. A good initial guess for amplitude will
111	* for a gamma 2 function be for instance (Max sample valu - baseline)*exp(2). If a Gamma N function is used
112	* multiply with exp(N).
113	* The correlation matrix for the parameters is written to the matrix pointe to by pCovar
114	* @param start the start index of the subarray of the data array.
115	* @param length the number of samples to use starting from "start"
116	* @param kmCovarPtrPtr the measurement correlation matrix
117	* @param pCovarPtrPtr the correlation matrix of the estimated parameters
118	* @param tGuess initial guess of t0 in entities of sampling intervals
119	* @param aGuess initial guess of the amplitude in entities of ADC levels*exp(2)
120	**/
121	void FitLeastMeanSquare(int start, int lenght, const double kmCovarPtrPtr, double pCovarPtrPtr);
122
123	/**
124	* This method gives the statistically best possible unbiased estimate of t0 and amplitude provided that a good
125	* estimate of the covariance matrix exists.
126	* For N number of samples this covariance matrix must be an NxN matrix.
127	* The default method is the Levenberg Marquardt nonlinear fit.
128	* The parameters "start" and "length" defines the subarray of the data array set in the
129	* constructor that will be used in the fit. "start + length" cannot exeed the length of the data array.
130	* The baseline must be subtracted before performing the fit othervise the result wil biased.
131	* A good initial guess will greatly enchance performance.
132	* a poor initial guess will slow down the algorithm and possibly lead to convergens to local minima
133	* of the cost function instead of the global minima. A very bad intial guess might give divergence of the solution.
134	* If the startindex of the pulse is the first sample above a low (2-3 ADC levels)
135	* above the baseline, then the initial guess of t0 can be set to zero. A good initial guess for amplitude will
136	* for a gamma 2 function be for instance (Max sample value - baseline)*exp(2). If a Gamma N function is used
137	* multiply with exp(N).
138	* @param start the start index of the subarray of the data array.
139	* @param length the number of samples to use starting from index
140	* @param mCovar the measurement correlation matrix
141	* @param pCovar the correlation matrix of the estimated parameters
142	* @param tGuess initial guess of t0 in entities of sampling intervals
143	* @param aGuess initial guess of the amplitude in entities of ADC levels*exp(2)
144	**/
145	void FitLeastMeanSquare(int start, int lenght, const double kmCovarPtrPtr, double pCovarPtrPtr, double tGuess, double aGuess);
146
147	/**
148	* Extraction of timing and energy using the Peakfinde Algorithm.
149	* The. The parameters "start" and "length" defines a sub array of the data array
150	* that will be used for the the fit. If start+length must not exeed the total length
151	* of the Data array. "start" must be chosen as close as possible to t0.
152	* The baseline must also be subtracted.
153	* The length of "tVector" and "aVector" mus be equal to length.
154	* "index + length" must not exeed the length of the data array set in the constructor.
155	* @param start the start index of the subarray of the data array.
156	* @param length the number of samples to use starting from index
157	* @param tVector the peakfinder vector for timing
158	* @param aVector the peakfinder vector for amplitude (energy)
159	**/
160	void FitPeakFinder(int start, int lenght, double tVector, double aVector);
161
162	/**
163	* This method finds the start index of the pulse (relative to the data array) by serching for
164	* three or more continious samples above trheshold.
165	* @param treshold treshold to use when searchin for the start of the pulse
166	**/
167	int FindStartIndex(double treshold);
168
169	/**
170	* Gives the timing in entities of sample indexes
171	* Physical time is found by multiplying with the sampling intervall (Ts).
172	**/
173	float GetTiming();
174
175	/**
176	* Gives the time in entities of ADC channels (quantization levels).
177	* Absolute enrgy is found by multiplying with offline calibration constants.
178	**/
179	float GetEnergy();
180
181	// /**
182	// * Set data array. Overrrides the data array set in the constructor.
183	// **/
184	// void SetData(int *data);
185
186	/**
187	* Set data array. Overrides data data array set in the constructor.
188	**/
189	void SetData(double *data);
190
191
192	private:
193
194	/**
195	* This function applies only to the Chi and Least mean square fit. An initial guess is made
196	* based on the average of the first 5 samples and the first value exeeding this value.
197	**/
198	void MakeInitialGuess();
199
200	/**
201	* This function applies only to the Chi and Least mean square fit. An initial guess is made
202	* based on the average of the first 5 samples and the first value exeeding threshold + this value.
203	* @param treshold The index of the first value above treshold is ntaken to be the first value.
204	**/
205	void MakeInitialGuess(int treshold);
206
207	double fFloatDataPtr; /<Float representation of data that should be fitted /
208	double fSampleFrequency; /*<The ADC sample frequency in MHz used under data taking /
209	double fTau; /*<The risetime in micro seconds/
210	double fDTof; /*<Time of flight in entities of sample intervals /
211	double fDAmpl; /*<Amplitude in entities of ADC levels/
212	double fDTofGuess; /*<Initial guess for t0/
213	double fDAmplGuess; /*<Initial guess for amplitude/
214	double kfMCovarPtrPtr; /<Covariance matrix of the measurements*/
215	double fPCovarPtrPtr; /<Covariance matrix of the estimated parameters*/
216	};
217
218	#endif