// @(#)root/hist:$Name:  $:$Id: TSpectrum2.h,v 1.3 2003/07/11 09:42:22 brun Exp $
// Author: Miroslav Morhac   27/05/99

/*************************************************************************
 * Copyright (C) 1995-2000, 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_TSpectrum2
#define ROOT_TSpectrum2

/////////////////////////////////////////////////////////////////////////////
//	THIS FILE CONTAINS HEADERS FOR ADVANCED				   //
//	SPECTRA PROCESSING FUNCTIONS.	   				   //
//									   //
//	ONE-DIMENSIONAL BACKGROUND ESTIMATION FUNCTION			   //
//	TWO-DIMENSIONAL BACKGROUND ESTIMATION FUNCTION			   //
//	ONE-DIMENSIONAL DECONVOLUTION FUNCTION				   //
//	TWO-DIMENSIONAL DECONVOLUTION FUNCTION				   //
//	ONE-DIMENSIONAL PEAK SEARCH FUNCTION				   //
//	TWO-DIMENSIONAL PEAK SEARCH FUNCTION				   //
//									   //
//	Miroslav Morhac							   //
//	Institute of Physics						   //
//	Slovak Academy of Sciences					   //
//	Dubravska cesta 9, 842 28 BRATISLAVA				   //
//	SLOVAKIA							   //
//									   //
//	email:fyzimiro@savba.sk,	 fax:+421 7 54772479		   //
//									   //
/////////////////////////////////////////////////////////////////////////////


#ifndef ROOT_TH2
#include "TH2.h"
#endif

const int BACK2_ORDER2 =0;
const int BACK2_ORDER4 =1;
const int BACK2_ORDER6 =2;
const int BACK2_ORDER8 =3;
const int BACK2_INCREASING_WINDOW =0;
const int BACK2_DECREASING_WINDOW =1;
const int BACK2_SUCCESSIVE_FILTERING =0;
const int BACK2_ONE_STEP_FILTERING =1;

const int SMOOTH2_3POINTS =3;
const int SMOOTH2_5POINTS =5;
const int SMOOTH2_7POINTS =7;
const int SMOOTH2_9POINTS =9;
const int SMOOTH2_11POINTS =11;
const int SMOOTH2_13POINTS =13;
const int SMOOTH2_15POINTS =15;

const bool SEARCH2_EXCLUDE_MARKOV =false;
const bool SEARCH2_INCLUDE_MARKOV =true;

const int FIT2_OPTIM_CHI_COUNTS =0;
const int FIT2_OPTIM_CHI_FUNC_VALUES =1;
const int FIT2_OPTIM_MAX_LIKELIHOOD =2;
const int FIT2_ALPHA_HALVING =0;
const int FIT2_ALPHA_OPTIMAL =1;
const int FIT2_FIT_POWER2 =2;
const int FIT2_FIT_POWER4 =4;
const int FIT2_FIT_POWER6 =6;
const int FIT2_FIT_POWER8 =8;
const int FIT2_FIT_POWER10 =10;
const int FIT2_FIT_POWER12 =12;
const int FIT2_TAYLOR_ORDER_FIRST =0;
const int FIT2_TAYLOR_ORDER_SECOND =1;
const int FIT2_NUM_OF_REGUL_CYCLES =100;

const int TRANSFORM2_HAAR =0;
const int TRANSFORM2_WALSH =1;
const int TRANSFORM2_COS =2;
const int TRANSFORM2_SIN =3;
const int TRANSFORM2_FOURIER =4;
const int TRANSFORM2_HARTLEY =5;
const int TRANSFORM2_FOURIER_WALSH =6;
const int TRANSFORM2_FOURIER_HAAR =7;
const int TRANSFORM2_WALSH_HAAR =8;
const int TRANSFORM2_COS_WALSH =9;
const int TRANSFORM2_COS_HAAR =10;
const int TRANSFORM2_SIN_WALSH =11;
const int TRANSFORM2_SIN_HAAR =12;
const int TRANSFORM2_FORWARD =0;
const int TRANSFORM2_INVERSE =1;

const int MAX_NUMBER_OF_PEAKS2= 1000;

	class TSpectrumTwoDimFit{
        public:
		int number_of_peaks;//input parameter, should be>0
		int number_of_iterations;//input parameter, should be >0
                int xmin;//first fitted channel in x direction
                int xmax;//last fitted channel in x direction
                int ymin;//first fitted channel in y direction
                int ymax;//last fitted channel in y direction
        	double alpha;//convergence coefficient, input parameter, it should be positive number and <=1
        	double chi; //here the function returns resulting chi square
                int statistic_type; //type of statistics, possible values FIT2_OPTIM_CHI_COUNTS (chi square statistics with counts as weighting coefficients), FIT2_OPTIM_CHI_FUNC_VALUES (chi square statistics with function values as weighting coefficients),FIT2_OPTIM_MAX_LIKELIHOOD
                int alpha_optim;//optimization of convergence coefficients, possible values FIT2_ALPHA_HALVING, FIT2_ALPHA_OPTIMAL, see manual
                int power;//possible values FIT21_FIT_POWER2,4,6,8,10,12, see manual
                int fit_taylor;//order of Taylor expansion, possible values FIT2_TAYLOR_ORDER_FIRST, FIT2_TAYLOR_ORDER_SECOND
		double position_init_x[MAX_NUMBER_OF_PEAKS2];//initial values of x positions of 2D peaks, input parameters
		double position_calc_x[MAX_NUMBER_OF_PEAKS2];//calculated values of fitted x positions of 2D peaks, output parameters
		double position_err_x[MAX_NUMBER_OF_PEAKS2];//x position errors of 2D peaks
                bool fix_position_x[MAX_NUMBER_OF_PEAKS2];//logical vector which allows to fix appropriate x positions of 2D peaks (not fit). However they are present in the estimated functional
		double position_init_y[MAX_NUMBER_OF_PEAKS2];//initial values of y positions of 2D peaks, input parameters
		double position_calc_y[MAX_NUMBER_OF_PEAKS2];//calculated values of fitted y positions of 2D peaks, output parameters
		double position_err_y[MAX_NUMBER_OF_PEAKS2];//y position errors of 2D peaks
                bool fix_position_y[MAX_NUMBER_OF_PEAKS2];//logical vector which allows to fix appropriate y positions of 2D peaks (not fit). However they are present in the estimated functional
		double position_init_x1[MAX_NUMBER_OF_PEAKS2];//initial values of x positions of 1D ridges, input parameters
		double position_calc_x1[MAX_NUMBER_OF_PEAKS2];//calculated values of fitted x positions of 1D ridges, output parameters
		double position_err_x1[MAX_NUMBER_OF_PEAKS2];//x position errors of 1D ridges
                bool fix_position_x1[MAX_NUMBER_OF_PEAKS2];//logical vector which allows to fix appropriate x positions of 1D ridges (not fit). However they are present in the estimated functional
		double position_init_y1[MAX_NUMBER_OF_PEAKS2];//initial values of y positions of 1D ridges, input parameters
		double position_calc_y1[MAX_NUMBER_OF_PEAKS2];//calculated values of fitted y positions of 1D ridges, output parameters
		double position_err_y1[MAX_NUMBER_OF_PEAKS2];//y position errors of 1D ridges
                bool fix_position_y1[MAX_NUMBER_OF_PEAKS2];//logical vector which allows to fix appropriate y positions of 1D ridges (not fit). However they are present in the estimated functional
        	double amp_init[MAX_NUMBER_OF_PEAKS2];//initial values of 2D peaks amplitudes, input parameters
	        double amp_calc[MAX_NUMBER_OF_PEAKS2];//calculated values of fitted amplitudes of 2D peaks, output parameters
        	double amp_err[MAX_NUMBER_OF_PEAKS2];//amplitude errors of 2D peaks
	        bool fix_amp[MAX_NUMBER_OF_PEAKS2];//logical vector which allows to fix appropriate amplitudes of 2D peaks (not fit). However they are present in the estimated functional
        	double amp_init_x1[MAX_NUMBER_OF_PEAKS2];//initial values of 1D ridges amplitudes, input parameters
	        double amp_calc_x1[MAX_NUMBER_OF_PEAKS2];//calculated values of fitted amplitudes of 1D ridges, output parameters
        	double amp_err_x1[MAX_NUMBER_OF_PEAKS2];//amplitude errors of 1D ridges
	        bool fix_amp_x1[MAX_NUMBER_OF_PEAKS2];//logical vector which allows to fix appropriate amplitudes of 1D ridges (not fit). However they are present in the estimated functional
        	double amp_init_y1[MAX_NUMBER_OF_PEAKS2];//initial values of 1D ridges amplitudes, input parameters
	        double amp_calc_y1[MAX_NUMBER_OF_PEAKS2];//calculated values of fitted amplitudes of 1D ridges, output parameters
        	double amp_err_y1[MAX_NUMBER_OF_PEAKS2];//amplitude errors of 1D ridges
	        bool fix_amp_y1[MAX_NUMBER_OF_PEAKS2];//logical vector which allows to fix appropriate amplitudes of 1D ridges (not fit). However they are present in the estimated functional
        	double volume[MAX_NUMBER_OF_PEAKS2];//calculated areas of peaks
	        double volume_err[MAX_NUMBER_OF_PEAKS2];//errors of peak areas
        	double sigma_init_x;//sigma x parameter, meaning analogical to the above given parameters, see manual
	        double sigma_calc_x;
        	double sigma_err_x;
	        bool fix_sigma_x;
        	double sigma_init_y;//sigma y parameter, meaning analogical to the above given parameters, see manual
	        double sigma_calc_y;
        	double sigma_err_y;
	        bool fix_sigma_y;
        	double ro_init;// correlation coefficient - see manual
	        double ro_calc;
        	double ro_err;
	        bool fix_ro;
        	double txy_init;//t parameter for 2D peaks, meaning analogical to the above given parameters, see manual
	        double txy_calc;
        	double txy_err;
	        bool fix_txy;
        	double sxy_init;//s parameter for 2D peaks, meaning analogical to the above given parameters, see manual
	        double sxy_calc;
        	double sxy_err;
	        bool fix_sxy;
        	double tx_init;//t parameter for 1D ridges (x direction), meaning analogical to the above given parameters, see manual
	        double tx_calc;
        	double tx_err;
	        bool fix_tx;
        	double ty_init;//t parameter for 1D ridges (y direction), meaning analogical to the above given parameters, see manual
	        double ty_calc;
        	double ty_err;
	        bool fix_ty;
        	double sx_init;//s parameter for 1D ridges (x direction), meaning analogical to the above given parameters, see manual
	        double sx_calc;
        	double sx_err;
	        bool fix_sx;
        	double sy_init;//s parameter for 1D ridges (y direction), meaning analogical to the above given parameters, see manual
	        double sy_calc;
        	double sy_err;
	        bool fix_sy;
        	double bx_init;//b parameter for 1D ridges (x direction), meaning analogical to the above given parameters, see manual
	        double bx_calc;
        	double bx_err;
	        bool fix_bx;
        	double by_init;//b parameter for 1D ridges (y direction), meaning analogical to the above given parameters, see manual
	        double by_calc;
        	double by_err;
	        bool fix_by;
        	double a0_init;//backgroud is estimated as a0+ax*x+ay*y
	        double a0_calc;
        	double a0_err;
	        bool fix_a0;
        	double ax_init;
	        double ax_calc;
        	double ax_err;
	        bool fix_ax;
        	double ay_init;
	        double ay_calc;
        	double ay_err;
	        bool fix_ay;
	};

class TSpectrum2 : public TNamed {
protected:
   Int_t         fMaxPeaks;       //Maximum number of peaks to be found
   Int_t         fNPeaks;         //number of peaks found
   Float_t      *fPosition;       //!array of current peak positions
   Float_t      *fPositionX;      //!X position of peaks
   Float_t      *fPositionY;      //!Y position of peaks
   Float_t       fResolution;     //resolution of the neighboring peaks
   TH1          *fHistogram;      //resulting histogram

public:
   TSpectrum2();
   TSpectrum2(Int_t maxpositions, Float_t resolution=1);
   virtual ~TSpectrum2();
   virtual const char *Background(TH1 *hist,int niter, Option_t *option="goff");
   const char *Background2(float **spectrum,int sizex,int sizey,int niter);
   const char *Deconvolution2(float **source,const float **resp,int sizex,int sizey,int niter);
   TH1        *GetHistogram() const {return fHistogram;}
   Int_t       GetNPeaks() const {return fNPeaks;}
   Float_t    *GetPositionX() const {return fPositionX;}
   Float_t    *GetPositionY() const {return fPositionY;}
   int         PeakEvaluate(const double *temp,int size,int xmax,double xmin,bool markov);
   virtual Int_t Search(TH1 *hist, Double_t sigma, Option_t *option="goff", Double_t threshold=0.05);
   void        SetResolution(Float_t resolution=1);

   //new functions April 2003
   const char *Background2RectangularRidges(float **spectrum,int sizex,int sizey,int number_of_iterations_x,int number_of_iterations_y,int direction,int filter_order,int filter_type);
   const char *Background2RectangularRidgesX(float **spectrum,int sizex,int sizey,int number_of_iterations,int direction,int filter_order);
   const char *Background2RectangularRidgesY(float **spectrum,int sizex,int sizey,int number_of_iterations,int direction,int filter_order);
   const char *Background2SkewRidges(float **spectrum,int sizex,int sizey,int number_of_iterations_x,int number_of_iterations_y,int direction,int filter_order);
   const char *Background2NonlinearRidges(float **spectrum,int sizex,int sizey,int number_of_iterations_x,int number_of_iterations_y,int direction,int filter_order);
   const char *Smooth2Markov(float **source, int sizex, int sizey, int aver_window);
   void        DecFourier2(double *working_space,int num,int iter,int inv);
   const char *Deconvolution2HighResolution(float** source,const float** resp,int sizex,int sizey,int number_of_iterations,int number_of_repetitions,double boost);
   Int_t       Search2HighRes(float **source,float **dest, int sizex, int sizey,
                              double sigma, double threshold,
                              bool background_remove,int decon_iterations,
                              bool markov, int aver_window);
   double      Lls(double a);

      //auxiliary functions for 2. parameter fit functions
   double Erfc(double x);
   double Derfc(double x);
   double Ourpowl(double a,int pw);
   void   StiefelInversion(double **a,int rozmer);

   double Shape2(int num_of_fitted_peaks,double x,double y,const double *parameter,double sigmax,double sigmay,double ro,double a0,double ax,double ay,double txy,double sxy,double tx,double ty,double sx,double sy,double bx,double by);
   double Deramp2(double x,double y,double x0,double y0,double sigmax,double sigmay,double ro,double txy,double sxy,double bx,double by);
   double Derampx(double x,double x0,double sigmax,double tx,double sx,double bx);
   double Deri02(double x,double y,double a,double x0,double y0,double sigmax,double sigmay,double ro,double txy,double sxy,double bx,double by);
   double Derderi02(double x,double y,double a,double x0,double y0,double sigmax,double sigmay,double ro);
   double Derj02(double x,double y,double a,double x0,double y0,double sigmax,double sigmay,double ro,double txy,double sxy,double bx,double by);
   double Derderj02(double x,double y,double a,double x0,double y0,double sigmax,double sigmay,double ro);
   double Deri01(double x,double ax,double x0,double sigmax,double tx,double sx,double bx);
   double Derderi01(double x,double ax,double x0,double sigmax);
   double Dersigmax(int num_of_fitted_peaks,double x,double y,const double *parameter,double sigmax,double sigmay,double ro,double txy,double sxy,double tx,double sx,double bx,double by);
   double Derdersigmax(int num_of_fitted_peaks,double x,double y,const double *parameter,double sigmax,double sigmay,double ro);
   double Dersigmay(int num_of_fitted_peaks,double x,double y,const double *parameter,double sigmax,double sigmay,double ro,double txy,double sxy,double ty,double sy,double bx,double by);
   double Derdersigmay(int num_of_fitted_peaks,double x,double y,const double *parameter,double sigmax,double sigmay,double ro);
   double Derro(int num_of_fitted_peaks,double x,double y,const double *parameter,double sx,double sy,double r);
   double Dertxy(int num_of_fitted_peaks,double x,double y,const double *parameter,double sigmax,double sigmay,double bx,double by);
   double Dersxy(int num_of_fitted_peaks,double x,double y,const double *parameter,double sigmax,double sigmay);
   double Dertx(int num_of_fitted_peaks,double x,const double *parameter,double sigmax,double bx);
   double Derty(int num_of_fitted_peaks,double x,const double *parameter,double sigmax,double bx);
   double Dersx(int num_of_fitted_peaks,double x,const double *parameter,double sigmax);
   double Dersy(int num_of_fitted_peaks,double x,const double *parameter,double sigmax);
   double Derbx(int num_of_fitted_peaks,double x,double y,const double *parameter,double sigmax,double sigmay,double txy,double tx,double bx,double by);
   double Derby(int num_of_fitted_peaks,double x,double y,const double *parameter,double sigmax,double sigmay,double txy,double ty,double bx,double by);
   double Derpa2(double sx,double sy,double ro);
   double Derpsigmax(double a,double sy,double ro);
   double Derpsigmay(double a,double sx,double ro);
   double Derpro(double a,double sx,double sy,double ro);
   double Volume(double a,double sx,double sy,double ro);

   const char* Fit2Awmi(float **source, TSpectrumTwoDimFit* p,int sizex,int sizey);
   const char* Fit2Stiefel(float **source, TSpectrumTwoDimFit* p,int sizex,int sizey);

//////////AUXILIARY FUNCTIONS FOR 2. DIMENSIONAL TRANSFORM, FILTER AND ENHANCE FUNCTIONS ////////////////////////
   void  Haar(float *working_space,int num,int direction);
   void  Walsh(float *working_space,int num);
   void  BitReverse(float *working_space,int num);
   void  Fourier(float *working_space,int num,int hartley,int direction,int zt_clear);
   void  BitReverseHaar(float *working_space,int shift,int num,int start);
   int   GeneralExe(float *working_space,int zt_clear,int num,int degree,int type);
   int   GeneralInv(float *working_space,int num,int degree,int type);

   void  HaarWalsh2(float **working_matrix,float *working_vector,int numx,int numy,int direction,int type);
   void  FourCos2(float **working_matrix,float *working_vector,int numx,int numy,int direction,int type);
   void  General2(float **working_matrix,float *working_vector,int numx,int numy,int direction,int type,int degree);

   const char *Transform2(const float **source,float **dest,int sizex,int sizey,int type,int direction,int degree);
   const char *Filter2Zonal(const float **source,float **dest,int sizex,int sizey,int type,int degree,int xmin, int xmax,int ymin,int ymax,float filter_coeff);
   const char *Enhance2(const float **source,float **dest,int sizex,int sizey,int type,int degree,int xmin, int xmax,int ymin,int ymax,float enhance_coeff);

   ClassDef(TSpectrum2,1)  //Peak Finder, background estimator, Deconvolution for 2-D histograms
};

#endif

