// @(#)root/spectrum:$Id$
// Author: Miroslav Morhac 25/09/06
/*************************************************************************
* Copyright (C) 1995-2006, Rene Brun and Fons Rademakers. *
* All rights reserved. *
* *
* For the licensing terms see $ROOTSYS/LICENSE. *
* For the list of contributors see $ROOTSYS/README/CREDITS. *
*************************************************************************/
#ifndef ROOT_TSpectrum2Fit
#define ROOT_TSpectrum2Fit
#include "TNamed.h"
class TSpectrum2Fit : public TNamed {
protected:
Int_t fNPeaks; ///< number of peaks present in fit, input parameter, it should be > 0
Int_t fNumberIterations; ///< number of iterations in fitting procedure, input parameter, it should be > 0
Int_t fXmin; ///< first fitted channel in x direction
Int_t fXmax; ///< last fitted channel in x direction
Int_t fYmin; ///< first fitted channel in y direction
Int_t fYmax; ///< last fitted channel in y direction
Int_t fStatisticType; ///< type of statistics, possible values kFitOptimChiCounts (chi square statistics with counts as weighting coefficients), kFitOptimChiFuncValues (chi square statistics with function values as weighting coefficients),kFitOptimMaxLikelihood
Int_t fAlphaOptim; ///< optimization of convergence algorithm, possible values kFitAlphaHalving, kFitAlphaOptimal
Int_t fPower; ///< possible values kFitPower2,4,6,8,10,12, for details see references. It applies only for Awmi fitting function.
Int_t fFitTaylor; ///< order of Taylor expansion, possible values kFitTaylorOrderFirst, kFitTaylorOrderSecond. It applies only for Awmi fitting function.
Double_t fAlpha; ///< convergence coefficient, input parameter, it should be positive number and <=1, for details see references
Double_t fChi; ///< here the fitting functions return resulting chi square
Double_t *fPositionInitX; ///< [fNPeaks] array of initial values of x positions of 2D peaks, input parameters
Double_t *fPositionCalcX; ///< [fNPeaks] array of calculated values of x positions of 2D peaks, output parameters
Double_t *fPositionErrX; ///< [fNPeaks] array of error values of x positions of 2D peaks, output parameters
Double_t *fPositionInitY; ///< [fNPeaks] array of initial values of y positions of 2D peaks, input parameters
Double_t *fPositionCalcY; ///< [fNPeaks] array of calculated values of y positions of 2D peaks, output parameters
Double_t *fPositionErrY; ///< [fNPeaks] array of error values of y positions of 2D peaks, output parameters
Double_t *fPositionInitX1; ///< [fNPeaks] array of initial x positions of 1D ridges, input parameters
Double_t *fPositionCalcX1; ///< [fNPeaks] array of calculated x positions of 1D ridges, output parameters
Double_t *fPositionErrX1; ///< [fNPeaks] array of x positions errors of 1D ridges, output parameters
Double_t *fPositionInitY1; ///< [fNPeaks] array of initial y positions of 1D ridges, input parameters
Double_t *fPositionCalcY1; ///< [fNPeaks] array of calculated y positions of 1D ridges, output parameters
Double_t *fPositionErrY1; ///< [fNPeaks] array of y positions errors of 1D ridges, output parameters
Double_t *fAmpInit; ///< [fNPeaks] array of initial values of amplitudes of 2D peaks, input parameters
Double_t *fAmpCalc; ///< [fNPeaks] array of calculated values of amplitudes of 2D peaks, output parameters
Double_t *fAmpErr; ///< [fNPeaks] array of amplitudes errors of 2D peaks, output parameters
Double_t *fAmpInitX1; ///< [fNPeaks] array of initial values of amplitudes of 1D ridges in x direction, input parameters
Double_t *fAmpCalcX1; ///< [fNPeaks] array of calculated values of amplitudes of 1D ridges in x direction, output parameters
Double_t *fAmpErrX1; ///< [fNPeaks] array of amplitudes errors of 1D ridges in x direction, output parameters
Double_t *fAmpInitY1; ///< [fNPeaks] array of initial values of amplitudes of 1D ridges in y direction, input parameters
Double_t *fAmpCalcY1; ///< [fNPeaks] array of calculated values of amplitudes of 1D ridges in y direction, output parameters
Double_t *fAmpErrY1; ///< [fNPeaks] array of amplitudes errors of 1D ridges in y direction, output parameters
Double_t *fVolume; ///< [fNPeaks] array of calculated volumes of 2D peaks, output parameters
Double_t *fVolumeErr; ///< [fNPeaks] array of volumes errors of 2D peaks, output parameters
Double_t fSigmaInitX; ///< initial value of sigma x parameter
Double_t fSigmaCalcX; ///< calculated value of sigma x parameter
Double_t fSigmaErrX; ///< error value of sigma x parameter
Double_t fSigmaInitY; ///< initial value of sigma y parameter
Double_t fSigmaCalcY; ///< calculated value of sigma y parameter
Double_t fSigmaErrY; ///< error value of sigma y parameter
Double_t fRoInit; ///< initial value of correlation coefficient
Double_t fRoCalc; ///< calculated value of correlation coefficient
Double_t fRoErr; ///< error value of correlation coefficient
Double_t fTxyInit; ///< initial value of t parameter for 2D peaks (relative amplitude of tail), for details see html manual and references
Double_t fTxyCalc; ///< calculated value of t parameter for 2D peaks
Double_t fTxyErr; ///< error value of t parameter for 2D peaks
Double_t fSxyInit; ///< initial value of s parameter for 2D peaks (relative amplitude of step), for details see html manual and references
Double_t fSxyCalc; ///< calculated value of s parameter for 2D peaks
Double_t fSxyErr; ///< error value of s parameter for 2D peaks
Double_t fTxInit; ///< initial value of t parameter for 1D ridges in x direction (relative amplitude of tail), for details see html manual and references
Double_t fTxCalc; ///< calculated value of t parameter for 1D ridges in x direction
Double_t fTxErr; ///< error value of t parameter for 1D ridges in x direction
Double_t fTyInit; ///< initial value of t parameter for 1D ridges in y direction (relative amplitude of tail), for details see html manual and references
Double_t fTyCalc; ///< calculated value of t parameter for 1D ridges in y direction
Double_t fTyErr; ///< error value of t parameter for 1D ridges in y direction
Double_t fSxInit; ///< initial value of s parameter for 1D ridges in x direction (relative amplitude of step), for details see html manual and references
Double_t fSxCalc; ///< calculated value of s parameter for 1D ridges in x direction
Double_t fSxErr; ///< error value of s parameter for 1D ridges in x direction
Double_t fSyInit; ///< initial value of s parameter for 1D ridges in y direction (relative amplitude of step), for details see html manual and references
Double_t fSyCalc; ///< calculated value of s parameter for 1D ridges in y direction
Double_t fSyErr; ///< error value of s parameter for 1D ridges in y direction
Double_t fBxInit; ///< initial value of b parameter for 1D ridges in x direction (slope), for details see html manual and references
Double_t fBxCalc; ///< calculated value of b parameter for 1D ridges in x direction
Double_t fBxErr; ///< error value of b parameter for 1D ridges in x direction
Double_t fByInit; ///< initial value of b parameter for 1D ridges in y direction (slope), for details see html manual and references
Double_t fByCalc; ///< calculated value of b parameter for 1D ridges in y direction
Double_t fByErr; ///< error value of b parameter for 1D ridges in y direction
Double_t fA0Init; ///< initial value of background a0 parameter(backgroud is estimated as a0+ax*x+ay*y)
Double_t fA0Calc; ///< calculated value of background a0 parameter
Double_t fA0Err; ///< error value of background a0 parameter
Double_t fAxInit; ///< initial value of background ax parameter(backgroud is estimated as a0+ax*x+ay*y)
Double_t fAxCalc; ///< calculated value of background ax parameter
Double_t fAxErr; ///< error value of background ax parameter
Double_t fAyInit; ///< initial value of background ay parameter(backgroud is estimated as a0+ax*x+ay*y)
Double_t fAyCalc; ///< calculated value of background ay parameter
Double_t fAyErr; ///< error value of background ay parameter
Bool_t *fFixPositionX; ///< [fNPeaks] array of logical values which allow to fix appropriate x positions of 2D peaks (not fit). However they are present in the estimated functional
Bool_t *fFixPositionY; ///< [fNPeaks] array of logical values which allow to fix appropriate y positions of 2D peaks (not fit). However they are present in the estimated functional
Bool_t *fFixPositionX1; ///< [fNPeaks] array of logical values which allow to fix appropriate x positions of 1D ridges (not fit). However they are present in the estimated functional
Bool_t *fFixPositionY1; ///< [fNPeaks] array of logical values which allow to fix appropriate y positions of 1D ridges (not fit). However they are present in the estimated functional
Bool_t *fFixAmp; ///< [fNPeaks] array of logical values which allow to fix appropriate amplitudes of 2D peaks (not fit). However they are present in the estimated functional
Bool_t *fFixAmpX1; ///< [fNPeaks] array of logical values which allow to fix appropriate amplitudes of 1D ridges in x direction (not fit). However they are present in the estimated functional
Bool_t *fFixAmpY1; ///< [fNPeaks] array of logical values which allow to fix appropriate amplitudes of 1D ridges in y direction (not fit). However they are present in the estimated functional
Bool_t fFixSigmaX; ///< logical value of sigma x parameter, which allows to fix the parameter (not to fit).
Bool_t fFixSigmaY; ///< logical value of sigma y parameter, which allows to fix the parameter (not to fit).
Bool_t fFixRo; ///< logical value of correlation coefficient, which allows to fix the parameter (not to fit).
Bool_t fFixTxy; ///< logical value of t parameter for 2D peaks, which allows to fix the parameter (not to fit).
Bool_t fFixSxy; ///< logical value of s parameter for 2D peaks, which allows to fix the parameter (not to fit).
Bool_t fFixTx; ///< logical value of t parameter for 1D ridges in x direction, which allows to fix the parameter (not to fit).
Bool_t fFixTy; ///< logical value of t parameter for 1D ridges in y direction, which allows to fix the parameter (not to fit).
Bool_t fFixSx; ///< logical value of s parameter for 1D ridges in x direction, which allows to fix the parameter (not to fit).
Bool_t fFixSy; ///< logical value of s parameter for 1D ridges in y direction, which allows to fix the parameter (not to fit).
Bool_t fFixBx; ///< logical value of b parameter for 1D ridges in x direction, which allows to fix the parameter (not to fit).
Bool_t fFixBy; ///< logical value of b parameter for 1D ridges in y direction, which allows to fix the parameter (not to fit).
Bool_t fFixA0; ///< logical value of a0 parameter, which allows to fix the parameter (not to fit).
Bool_t fFixAx; ///< logical value of ax parameter, which allows to fix the parameter (not to fit).
Bool_t fFixAy; ///< logical value of ay parameter, which allows to fix the parameter (not to fit).
public:
enum {
kFitOptimChiCounts =0,
kFitOptimChiFuncValues =1,
kFitOptimMaxLikelihood =2,
kFitAlphaHalving =0,
kFitAlphaOptimal =1,
kFitPower2 =2,
kFitPower4 =4,
kFitPower6 =6,
kFitPower8 =8,
kFitPower10 =10,
kFitPower12 =12,
kFitTaylorOrderFirst =0,
kFitTaylorOrderSecond =1,
kFitNumRegulCycles =100
};
TSpectrum2Fit(void); //default constructor
TSpectrum2Fit(Int_t numberPeaks);
virtual ~TSpectrum2Fit();
//auxiliary functions for 2. parameter fit functions
protected:
Double_t Deramp2(Double_t x,Double_t y,Double_t x0,Double_t y0,Double_t sigmax,Double_t sigmay,Double_t ro,Double_t txy,Double_t sxy,Double_t bx,Double_t by);
Double_t Derampx(Double_t x,Double_t x0,Double_t sigmax,Double_t tx,Double_t sx,Double_t bx);
Double_t Derbx(Int_t numOfFittedPeaks,Double_t x,Double_t y,const Double_t *parameter,Double_t sigmax,Double_t sigmay,Double_t txy,Double_t tx,Double_t bx,Double_t by);
Double_t Derby(Int_t numOfFittedPeaks,Double_t x,Double_t y,const Double_t *parameter,Double_t sigmax,Double_t sigmay,Double_t txy,Double_t ty,Double_t bx,Double_t by);
Double_t Derderi01(Double_t x,Double_t ax,Double_t x0,Double_t sigmax);
Double_t Derderi02(Double_t x,Double_t y,Double_t a,Double_t x0,Double_t y0,Double_t sigmax,Double_t sigmay,Double_t ro);
Double_t Derderj02(Double_t x,Double_t y,Double_t a,Double_t x0,Double_t y0,Double_t sigmax,Double_t sigmay,Double_t ro);
Double_t Derdersigmax(Int_t numOfFittedPeaks,Double_t x,Double_t y,const Double_t *parameter,Double_t sigmax,Double_t sigmay,Double_t ro);
Double_t Derdersigmay(Int_t numOfFittedPeaks,Double_t x,Double_t y,const Double_t *parameter,Double_t sigmax,Double_t sigmay,Double_t ro);
Double_t Derfc(Double_t x);
Double_t Deri01(Double_t x,Double_t ax,Double_t x0,Double_t sigmax,Double_t tx,Double_t sx,Double_t bx);
Double_t Deri02(Double_t x,Double_t y,Double_t a,Double_t x0,Double_t y0,Double_t sigmax,Double_t sigmay,Double_t ro,Double_t txy,Double_t sxy,Double_t bx,Double_t by);
Double_t Derj02(Double_t x,Double_t y,Double_t a,Double_t x0,Double_t y0,Double_t sigmax,Double_t sigmay,Double_t ro,Double_t txy,Double_t sxy,Double_t bx,Double_t by);
Double_t Derpa2(Double_t sx,Double_t sy,Double_t ro);
Double_t Derpro(Double_t a,Double_t sx,Double_t sy,Double_t ro);
Double_t Derpsigmax(Double_t a,Double_t sy,Double_t ro);
Double_t Derpsigmay(Double_t a,Double_t sx,Double_t ro);
Double_t Derro(Int_t numOfFittedPeaks,Double_t x,Double_t y,const Double_t *parameter,Double_t sx,Double_t sy,Double_t r);
Double_t Dersigmax(Int_t numOfFittedPeaks,Double_t x,Double_t y,const Double_t *parameter,Double_t sigmax,Double_t sigmay,Double_t ro,Double_t txy,Double_t sxy,Double_t tx,Double_t sx,Double_t bx,Double_t by);
Double_t Dersigmay(Int_t numOfFittedPeaks,Double_t x,Double_t y,const Double_t *parameter,Double_t sigmax,Double_t sigmay,Double_t ro,Double_t txy,Double_t sxy,Double_t ty,Double_t sy,Double_t bx,Double_t by);
Double_t Dersx(Int_t numOfFittedPeaks,Double_t x,const Double_t *parameter,Double_t sigmax);
Double_t Dersxy(Int_t numOfFittedPeaks,Double_t x,Double_t y,const Double_t *parameter,Double_t sigmax,Double_t sigmay);
Double_t Dersy(Int_t numOfFittedPeaks,Double_t x,const Double_t *parameter,Double_t sigmax);
Double_t Dertx(Int_t numOfFittedPeaks,Double_t x,const Double_t *parameter,Double_t sigmax,Double_t bx);
Double_t Dertxy(Int_t numOfFittedPeaks,Double_t x,Double_t y,const Double_t *parameter,Double_t sigmax,Double_t sigmay,Double_t bx,Double_t by);
Double_t Derty(Int_t numOfFittedPeaks,Double_t x,const Double_t *parameter,Double_t sigmax,Double_t bx);
Double_t Erfc(Double_t x);
Double_t Ourpowl(Double_t a,Int_t pw);
Double_t Shape2(Int_t numOfFittedPeaks,Double_t x,Double_t y,const Double_t *parameter,Double_t sigmax,Double_t sigmay,Double_t ro,Double_t a0,Double_t ax,Double_t ay,Double_t txy,Double_t sxy,Double_t tx,Double_t ty,Double_t sx,Double_t sy,Double_t bx,Double_t by);
void StiefelInversion(Double_t **a,Int_t size);
Double_t Volume(Double_t a,Double_t sx,Double_t sy,Double_t ro);
public:
void FitAwmi(Double_t **source);
void FitStiefel(Double_t **source);
void GetAmplitudes(Double_t *amplitudes, Double_t *amplitudesX1, Double_t *amplitudesY1);
void GetAmplitudeErrors(Double_t *amplitudeErrors, Double_t *amplitudeErrorsX1, Double_t *amplitudeErrorsY1);
void GetBackgroundParameters(Double_t &a0, Double_t &a0Err, Double_t &ax, Double_t &axErr, Double_t &ay, Double_t &ayErr);
Double_t GetChi() const {return fChi;}
void GetPositions(Double_t *positionsX, Double_t *positionsY, Double_t *positionsX1, Double_t *positionsY1);
void GetPositionErrors(Double_t *positionErrorsX, Double_t *positionErrorsY, Double_t *positionErrorsX1, Double_t *positionErrorsY1);
void GetRo(Double_t &ro, Double_t &roErr);
void GetSigmaX(Double_t &sigmaX, Double_t &sigmaErrX);
void GetSigmaY(Double_t &sigmaY, Double_t &sigmaErrY);
void GetTailParameters(Double_t &txy, Double_t &txyErr, Double_t &tx, Double_t &txErr, Double_t &ty, Double_t &tyErr, Double_t &bx, Double_t &bxErr, Double_t &by, Double_t &byErr, Double_t &sxy, Double_t &sxyErr, Double_t &sx, Double_t &sxErr, Double_t &sy, Double_t &syErr);
void GetVolumes(Double_t *volumes);
void GetVolumeErrors(Double_t *volumeErrors);
void SetBackgroundParameters(Double_t a0Init, Bool_t fixA0, Double_t axInit, Bool_t fixAx, Double_t ayInit, Bool_t fixAy);
void SetFitParameters(Int_t xmin,Int_t xmax,Int_t ymin,Int_t ymax, Int_t numberIterations, Double_t alpha, Int_t statisticType, Int_t alphaOptim, Int_t power, Int_t fitTaylor);
void SetPeakParameters(Double_t sigmaX, Bool_t fixSigmaX, Double_t sigmaY, Bool_t fixSigmaY, Double_t ro, Bool_t fixRo, const Double_t *positionInitX, const Bool_t *fixPositionX, const Double_t *positionInitY, const Bool_t *fixPositionY, const Double_t *positionInitX1, const Bool_t *fixPositionX1, const Double_t *positionInitY1, const Bool_t *fixPositionY1, const Double_t *ampInit, const Bool_t *fixAmp, const Double_t *ampInitX1, const Bool_t *fixAmpX1, const Double_t *ampInitY1, const Bool_t *fixAmpY1);
void SetTailParameters(Double_t tInitXY, Bool_t fixTxy, Double_t tInitX, Bool_t fixTx, Double_t tInitY, Bool_t fixTy, Double_t bInitX, Bool_t fixBx, Double_t bInitY, Bool_t fixBy, Double_t sInitXY, Bool_t fixSxy, Double_t sInitX, Bool_t fixSx, Double_t sInitY, Bool_t fixSy);
ClassDef(TSpectrum2Fit,1) //Spectrum2 Fitter using algorithm without matrix inversion and conjugate gradient method for symmetrical matrices (Stiefel-Hestens method)
};
#endif