// @(#)root/spectrum:$Name:  $:$Id: TSpectrumFit.h,v 1.2 2006/10/02 10:58:42 brun Exp $
// 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_TSpectrumFit
#define ROOT_TSpectrumFit

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// TSpectrumFit                                                         //
//                                                                      //
// Class for fitting 1D spectra using AWMI (algorithm without matrix    //
// inversion) and conjugate gradient algorithms for symmetrical         //
// matrices (Stiefel-Hestens method). AWMI method allows to fit         //
// simulaneously 100s up to 1000s peaks. Stiefel method is very stable, //
// it converges faster, but is more time consuming                      //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

#ifndef ROOT_TH1
#include "TH1.h"
#endif

class TSpectrumFit : 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
   Int_t     fXmax;                      //last fitted channel
   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 *fPositionInit;              //[fNPeaks] array of initial values of peaks positions, input parameters
   Double_t *fPositionCalc;              //[fNPeaks] array of calculated values of fitted positions, output parameters
   Double_t *fPositionErr;               //[fNPeaks] array of position errors
   Double_t *fAmpInit;                   //[fNPeaks] array of initial values of peaks amplitudes, input parameters
   Double_t *fAmpCalc;                   //[fNPeaks] array of calculated values of fitted amplitudes, output parameters
   Double_t *fAmpErr;                    //[fNPeaks] array of amplitude errors
   Double_t *fArea;                      //[fNPeaks] array of calculated areas of peaks
   Double_t *fAreaErr;                   //[fNPeaks] array of errors of peak areas
   Double_t  fSigmaInit;                 //initial value of sigma parameter
   Double_t  fSigmaCalc;                 //calculated value of sigma parameter
   Double_t  fSigmaErr;                  //error value of sigma parameter
   Double_t  fTInit;                     //initial value of t parameter (relative amplitude of tail), for details see html manual and references
   Double_t  fTCalc;                     //calculated value of t parameter
   Double_t  fTErr;                      //error value of t parameter
   Double_t  fBInit;                     //initial value of b parameter (slope), for details see html manual and references
   Double_t  fBCalc;                     //calculated value of b parameter
   Double_t  fBErr;                      //error value of b parameter
   Double_t  fSInit;                     //initial value of s parameter (relative amplitude of step), for details see html manual and references
   Double_t  fSCalc;                     //calculated value of s parameter
   Double_t  fSErr;                      //error value of s parameter
   Double_t  fA0Init;                    //initial value of background a0 parameter(backgroud is estimated as a0+a1*x+a2*x*x)
   Double_t  fA0Calc;                    //calculated value of background a0 parameter
   Double_t  fA0Err;                     //error value of background a0 parameter
   Double_t  fA1Init;                    //initial value of background a1 parameter(backgroud is estimated as a0+a1*x+a2*x*x)
   Double_t  fA1Calc;                    //calculated value of background a1 parameter
   Double_t  fA1Err;                     //error value of background a1 parameter
   Double_t  fA2Init;                    //initial value of background a2 parameter(backgroud is estimated as a0+a1*x+a2*x*x)
   Double_t  fA2Calc;                    //calculated value of background a2 parameter
   Double_t  fA2Err;                     //error value of background a2 parameter
   Bool_t   *fFixPosition;               //[fNPeaks] array of logical values which allow to fix appropriate positions (not fit). However they are present in the estimated functional   
   Bool_t   *fFixAmp;                    //[fNPeaks] array of logical values which allow to fix appropriate amplitudes (not fit). However they are present in the estimated functional      
   Bool_t    fFixSigma;                  //logical value of sigma parameter, which allows to fix the parameter (not to fit).   
   Bool_t    fFixT;                      //logical value of t parameter, which allows to fix the parameter (not to fit).      
   Bool_t    fFixB;                      //logical value of b parameter, which allows to fix the parameter (not to fit).   
   Bool_t    fFixS;                      //logical value of s parameter, 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    fFixA1;                     //logical value of a1 parameter, which allows to fix the parameter (not to fit).   
   Bool_t    fFixA2;                     //logical value of a2 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
   };
   TSpectrumFit(void); //default constructor     
   TSpectrumFit(Int_t numberPeaks);
   virtual ~TSpectrumFit();

   //auxiliary functions for 1. parameter fit functions
protected:   
   Double_t            Area(Double_t a,Double_t sigma,Double_t t,Double_t b);
   Double_t            Dera1(Double_t i);
   Double_t            Dera2(Double_t i);
   Double_t            Deramp(Double_t i,Double_t i0,Double_t sigma,Double_t t,Double_t s,Double_t b);
   Double_t            Derb(Int_t num_of_fitted_peaks,Double_t i,const Double_t* parameter,Double_t sigma,Double_t t,Double_t b);
   Double_t            Derderi0(Double_t i,Double_t amp,Double_t i0,Double_t sigma);
   Double_t            Derdersigma(Int_t num_of_fitted_peaks,Double_t i,const Double_t* parameter,Double_t sigma);
   Double_t            Derfc(Double_t x);
   Double_t            Deri0(Double_t i,Double_t amp,Double_t i0,Double_t sigma,Double_t t,Double_t s,Double_t b);
   Double_t            Derpa(Double_t sigma,Double_t t,Double_t b);
   Double_t            Derpb(Double_t a,Double_t sigma,Double_t t,Double_t b);
   Double_t            Derpsigma(Double_t a,Double_t t,Double_t b);
   Double_t            Derpt(Double_t a,Double_t sigma,Double_t b);
   Double_t            Ders(Int_t num_of_fitted_peaks,Double_t i,const Double_t* parameter,Double_t sigma);
   Double_t            Dersigma(Int_t num_of_fitted_peaks,Double_t i,const Double_t* parameter,Double_t sigma,Double_t t,Double_t s,Double_t b);
   Double_t            Dert(Int_t num_of_fitted_peaks,Double_t i,const Double_t* parameter,Double_t sigma,Double_t b);
   Double_t            Erfc(Double_t x);
   Double_t            Ourpowl(Double_t a,Int_t pw);
   Double_t            Shape(Int_t num_of_fitted_peaks,Double_t i,const Double_t *parameter,Double_t sigma,Double_t t,Double_t s,Double_t b,Double_t a0,Double_t a1,Double_t a2);
   void                StiefelInversion(Double_t **a,Int_t rozmer);

public:
   void                FitAwmi(float *source);
   void                FitStiefel(float *source);
   Double_t           *GetAmplitudes() const {return fAmpCalc;}
   Double_t           *GetAmplitudesErrors() const {return fAmpErr;}
   Double_t           *GetAreas() const {return fArea;}
   Double_t           *GetAreasErrors() const {return fAreaErr;}
   void                GetBackgroundParameters(Double_t &a0, Double_t &a0Err, Double_t &a1, Double_t &a1Err, Double_t &a2, Double_t &a2Err);
   Double_t            GetChi() const {return fChi;}
   Double_t           *GetPositions() const {return fPositionCalc;}
   Double_t           *GetPositionsErrors() const {return fPositionErr;}
   void                GetSigma(Double_t &sigma, Double_t &sigmaErr);
   void                GetTailParameters(Double_t &t, Double_t &tErr, Double_t &b, Double_t &bErr, Double_t &s, Double_t &sErr);
   void                SetBackgroundParameters(Double_t a0Init, Bool_t fixA0, Double_t a1Init, Bool_t fixA1, Double_t a2Init, Bool_t fixA2);
   void                SetFitParameters(Int_t xmin,Int_t xmax, Int_t numberIterations, Double_t alpha, Int_t statisticType, Int_t alphaOptim, Int_t power, Int_t fitTaylor);
   void                SetPeakParameters(Double_t sigma, Bool_t fixSigma, const Float_t *positionInit, const Bool_t *fixPosition, const Float_t *ampInit, const Bool_t *fixAmp);
   void                SetTailParameters(Double_t tInit, Bool_t fixT, Double_t bInit, Bool_t fixB, Double_t sInit, Bool_t fixS);

   ClassDef(TSpectrumFit,1)  //Spectrum Fitter using algorithm without matrix inversion and conjugate gradient method for symmetrical matrices (Stiefel-Hestens method)
};

#endif