class TSpectrum2Fit: public TNamed

numberPeaks: number of fitted peaks (must be greater than zero)
the constructor allocates arrays for all fitted parameters (peak positions, amplitudes etc) and sets the member
variables to their default values. One can change these variables by member functions (setters) of TSpectrumFit class.

Shape function of the fitted peaks contains the two-dimensional symmetrical Gaussian two one-dimensional symmetrical Gaussian ridges as well as nonsymmetrical terms and background.

 TWO-DIMENSIONAL FIT FUNCTION
  ALGORITHM WITHOUT MATRIX INVERSION
  This function fits the source spectrum. The calling program should
  fill in input parameters of the TSpectrum2Fit class.
  The fitted parameters are written into
  TSpectrum2Fit class output parameters and fitted data are written into
  source spectrum.

        Function parameters:
        source-pointer to the matrix of source spectrum

Fitting

 

Goal: to estimate simultaneously peak shape parameters in spectra with large number of peaks

 

•         peaks can be fitted separately, each peak (or multiplets) in a region or together all peaks in a spectrum. To fit separately each peak one needs to determine the fitted region. However it can happen that the regions of neighboring peaks are overlapping. Then the results of fitting are very poor. On the other hand, when fitting together all peaks found in a  spectrum, one needs to have a method that is  stable (converges) and fast enough to carry out fitting in reasonable time

•         we have implemented the nonsymmetrical semiempirical peak shape function

•         it contains the two-dimensional symmetrical Gaussian two one-dimensional symmetrical Gaussian ridges as well as nonsymmetrical terms and background.

where Txy, Tx, Ty, Sxy, Sx, Sy are relative amplitudes and Bx, By are slopes.

 

•         algorithm without matrix inversion (AWMI) allows fitting tens, hundreds of peaks simultaneously that represent sometimes thousands of parameters [2], [5].

 

 

Function:

void TSpectrumFit2::FitAwmi(float **fSource)

 

This function fits the source spectrum using AWMI algorithm. The calling program should fill in input parameters of the TSpectrumFit2 class using a set of TSpectrumFit2 setters. The fitted parameters are written into the class and fitted data are written into source spectrum.

 

 

Parameter:

        fSource-pointer to the matrix of source spectrum                 

 

 

Member variables of  TSpectrumFit2 class:

   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).

 

References:

[1] Phillps G.W., Marlow K.W., NIM 137 (1976) 525.

[2] I. A. Slavic: Nonlinear least-squares fitting without matrix inversion applied to complex Gaussian spectra analysis. NIM 134 (1976) 285-289.

[3] T. Awaya: A new method for curve fitting to the data with low statistics not using chi-square method. NIM 165 (1979) 317-323.

[4] T. Hauschild, M. Jentschel: Comparison of maximum likelihood estimation and chi-square statistics applied to counting experiments. NIM A 457 (2001) 384-401.

 [5]  M. Morháč,  J. Kliman,  M. Jandel,  Ľ. Krupa, V. Matoušek: Study of fitting algorithms applied to simultaneous analysis of large number of peaks in -ray spectra. Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003

 

Example  – script FitAwmi2.c:

Fig. 1 Original two-dimensional spectrum with found peaks (using TSpectrum2 peak searching function). The positions of peaks were used as initial estimates in fitting procedure.

Fig. 2 Fitted (generated from fitted parameters) spectrum of the data from Fig. 1 using Algorithm Without Matrix Inversion. Each peak was represented by 7 parameters, which together with Sigmax, Sigmay and a0 resulted in 38 fitted parameters. The chi-square after 1000 iterations was 0.642342.

 

Script:

 

// Example to illustrate fitting function, algorithm without matrix inversion (AWMI) (class TSpectrumFit2).

// To execute this example, do

// root > .x FitAwmi2.C

 

void FitAwmi2() {

   Int_t i, j, nfound;

   Int_t nbinsx = 64;

   Int_t nbinsy = 64;  

   Int_t xmin  = 0;

   Int_t xmax  = nbinsx;

   Int_t ymin  = 0;

   Int_t ymax  = nbinsy;

   Float_t ** source = new float *[nbinsx];  

   Float_t ** dest = new float *[nbinsx];     

   for (i=0;i<nbinsx;i++)

                                                source[i]=new float[nbinsy];

   for (i=0;i<nbinsx;i++)

                                                dest[i]=new float[nbinsy];

   TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);

   TFile *f = new TFile("TSpectrum2.root");

   search=(TH2F*) f->Get("search4;1");

   TCanvas *Searching = new TCanvas("Searching","Two-dimensional fitting using Algorithm Without Matrix Inversion",10,10,1000,700);

   TSpectrum2 *s = new TSpectrum2();

   for (i = 0; i < nbinsx; i++){

     for (j = 0; j < nbinsy; j++){

                    source[i][j] = search->GetBinContent(i + 1,j + 1);

                 }

   }

   //searching for candidate peaks positions    

   nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 5, kTRUE, 3, kFALSE, 3);  

   Bool_t *FixPosX = new Bool_t[nfound];

   Bool_t *FixPosY = new Bool_t[nfound];  

   Bool_t *FixAmp = new Bool_t[nfound];     

   Float_t *PosX = new Float_t[nfound];        

   Float_t *PosY = new Float_t[nfound];

   Float_t *Amp = new Float_t[nfound];     

   Float_t *AmpXY = new Float_t[nfound];        

   PosX = s->GetPositionX();

   PosY = s->GetPositionY();      

   printf("Found %d candidate peaks\n",nfound);  

   for(i = 0; i< nfound ; i++){

      FixPosX[i] = kFALSE;

      FixPosY[i] = kFALSE;     

      FixAmp[i] = kFALSE;   

      Amp[i] = source[(int)(PosX[i]+0.5)][(int)(PosY[i]+0.5)];      //initial values of peaks amplitudes, input parameters         

      AmpXY[i] = 0;

   }

   //filling in the initial estimates of the input parameters

   TSpectrumFit2 *pfit=new TSpectrumFit2(nfound);

   pfit->SetFitParameters(xmin, xmax-1, ymin, ymax-1, 1000, 0.1, pfit->kFitOptimChiCounts, pfit->kFitAlphaHalving, pfit->kFitPower2, pfit->kFitTaylorOrderFirst);  

   pfit->SetPeakParameters(2, kFALSE, 2, kFALSE, 0, kTRUE, PosX, (Bool_t *) FixPosX, PosY, (Bool_t *) FixPosY, PosX, (Bool_t *) FixPosX, PosY, (Bool_t *) FixPosY, Amp, (Bool_t *) FixAmp, AmpXY, (Bool_t *) FixAmp, AmpXY, (Bool_t *) FixAmp);     

   pfit->SetBackgroundParameters(0, kFALSE, 0, kTRUE, 0, kTRUE);  

   pfit->FitAwmi(source);

    for (i = 0; i < nbinsx; i++){

     for (j = 0; j < nbinsy; j++){

                  search->SetBinContent(i + 1, j + 1,source[i][j]);

                 }

   }  

   search->Draw("SURF");

}

 

Example 2 – script FitA2.c:

Fig. 3 Original two-dimensional gamma-gamma-ray spectrum with found peaks (using TSpectrum2 peak searching function).

Fig. 4 Fitted (generated from fitted parameters) spectrum of the data from Fig. 3 using Algorithm Without Matrix Inversion. 152 peaks were identified. Each peak was represented by 7 parameters, which together with Sigmax, Sigmay and a0 resulted in 1067 fitted parameters. The chi-square after 1000 iterations was 0.728675. One can observe good correspondence with the original data.

 

Script:

// Example to illustrate fitting function, algorithm without matrix inversion (AWMI) (class TSpectrumFit2).

// To execute this example, do

// root > .x FitA2.C

void FitA2() {

   Int_t i, j, nfound;

   Int_t nbinsx = 256;

   Int_t nbinsy = 256;  

   Int_t xmin  = 0;

   Int_t xmax  = nbinsx;

   Int_t ymin  = 0;

   Int_t ymax  = nbinsy;

   Float_t ** source = new float *[nbinsx];  

   Float_t ** dest = new float *[nbinsx];     

   for (i=0;i<nbinsx;i++)

                                                source[i]=new float[nbinsy];

   for (i=0;i<nbinsx;i++)

                                                dest[i]=new float[nbinsy];  

   TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);

   TFile *f = new TFile("TSpectrum2.root");

   search=(TH2F*) f->Get("fit1;1");

   TCanvas *Searching = new TCanvas("Searching","Two-dimensional fitting using Algorithm Without Matrix Inversion",10,10,1000,700);

   TSpectrum2 *s = new TSpectrum2(1000,1);

   for (i = 0; i < nbinsx; i++){

     for (j = 0; j < nbinsy; j++){

                    source[i][j] = search->GetBinContent(i + 1,j + 1);

                 }

   }  

   nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 2, kTRUE, 100, kFALSE, 3);  

   printf("Found %d candidate peaks\n",nfound);

   Bool_t *FixPosX = new Bool_t[nfound];

   Bool_t *FixPosY = new Bool_t[nfound];  

   Bool_t *FixAmp = new Bool_t[nfound];     

   Float_t *PosX = new Float_t[nfound];        

   Float_t *PosY = new Float_t[nfound];

   Float_t *Amp = new Float_t[nfound];     

   Float_t *AmpXY = new Float_t[nfound];        

   PosX = s->GetPositionX();

   PosY = s->GetPositionY();     

   for(i = 0; i< nfound ; i++){

      FixPosX[i] = kFALSE;

      FixPosY[i] = kFALSE;     

      FixAmp[i] = kFALSE;   

      Amp[i] = source[(int)(PosX[i]+0.5)][(int)(PosY[i]+0.5)];      //initial values of peaks amplitudes, input parameters         

      AmpXY[i] = 0;

   }

   filling in the initial estimates of the input parameters

   TSpectrumFit2 *pfit=new TSpectrumFit2(nfound);

   pfit->SetFitParameters(xmin, xmax-1, ymin, ymax-1, 1000, 0.1, pfit->kFitOptimChiCounts, pfit->kFitAlphaHalving, pfit->kFitPower2, pfit->kFitTaylorOrderFirst);  

   pfit->SetPeakParameters(2, kFALSE, 2, kFALSE, 0, kTRUE, PosX, (Bool_t *) FixPosX, PosY, (Bool_t *) FixPosY, PosX, (Bool_t *) FixPosX, PosY, (Bool_t *) FixPosY, Amp, (Bool_t *) FixAmp, AmpXY, (Bool_t *) FixAmp, AmpXY, (Bool_t *) FixAmp);     

   pfit->SetBackgroundParameters(0, kFALSE, 0, kTRUE, 0, kTRUE);  

   pfit->FitAwmi(source);

   for (i = 0; i < nbinsx; i++){

     for (j = 0; j < nbinsy; j++){

                  search->SetBinContent(i + 1, j + 1,source[i][j]);

                 }

   }  

   search->Draw("SURF");  

}

 

 TWO-DIMENSIONAL FIT FUNCTION USING STIEFEL-HESTENS
  ALGORITHM.
  This function fits the source spectrum. The calling program should
  fill in input parameters of the TSpectrum2Fit class.
  The fitted parameters are written into
  TSpectrum2Fit class output parameters and fitted data are written into
  source spectrum.

        Function parameters:
        source-pointer to the matrix of source spectrum

Stiefel fitting algorithm

 

Function:

void TSpectrumFit2::FitStiefel(float **fSource)

This function fits the source spectrum using Stiefel-Hestens method [1].  The calling program should fill in input fitting parameters of the TSpectrumFit2 class using a set of TSpectrumFit2 setters. The fitted parameters are written into the class and the fitted data are written into source spectrum. It converges faster than Awmi method.

 

Parameter:

        fSource-pointer to the matrix of source spectrum                 

       

Reference:

[1] B. Mihaila: Analysis of complex gamma spectra, Rom. Jorn. Phys., Vol. 39, No. 2, (1994), 139-148.

 

Example 1 – script FitS.c:

Fig. 1 Original two-dimensional spectrum with found peaks (using TSpectrum2 peak searching function). The positions of peaks were used as initial estimates in fitting procedure.

Fig. 2 Fitted (generated from fitted parameters) spectrum of the data from Fig. 1 using Stiefel-Hestens method. Each peak was represented by 7 parameters, which together with Sigmax, Sigmay and a0 resulted in 38 fitted parameters. The chi-square after 1000 iterations was 0.642157.

 

Script:

// Example to illustrate fitting function, algorithm without matrix inversion (AWMI) (class TSpectrumFit2).

// To execute this example, do

// root > .x FitStiefel2.C

 

void FitStiefel2() {

   Int_t i, j, nfound;

   Int_t nbinsx = 64;

   Int_t nbinsy = 64;  

   Int_t xmin  = 0;

   Int_t xmax  = nbinsx;

   Int_t ymin  = 0;

   Int_t ymax  = nbinsy;

   Float_t ** source = new float *[nbinsx];  

   Float_t ** dest = new float *[nbinsx];     

   for (i=0;i<nbinsx;i++)

                                                source[i]=new float[nbinsy];

   for (i=0;i<nbinsx;i++)

                                                dest[i]=new float[nbinsy];

   TH2F *search = new TH2F("search","High resolution peak searching",nbinsx,xmin,xmax,nbinsy,ymin,ymax);

   TFile *f = new TFile("TSpectrum2.root");

   search=(TH2F*) f->Get("search4;1");

   TCanvas *Searching = new TCanvas("Searching","Two-dimensional fitting using Stiefel-Hestens method",10,10,1000,700);

   TSpectrum2 *s = new TSpectrum2();

   for (i = 0; i < nbinsx; i++){

     for (j = 0; j < nbinsy; j++){

                    source[i][j] = search->GetBinContent(i + 1,j + 1);

                 }

   }  

   nfound = s->SearchHighRes(source, dest, nbinsx, nbinsy, 2, 5, kTRUE, 3, kFALSE, 3);  

   printf("Found %d candidate peaks\n",nfound);

   Bool_t *FixPosX = new Bool_t[nfound];

   Bool_t *FixPosY = new Bool_t[nfound];   

   Bool_t *FixAmp = new Bool_t[nfound];     

   Float_t *PosX = new Float_t[nfound];        

   Float_t *PosY = new Float_t[nfound];

   Float_t *Amp = new Float_t[nfound];     

   Float_t *AmpXY = new Float_t[nfound];        

   PosX = s->GetPositionX();

   PosY = s->GetPositionY();     

   for(i = 0; i< nfound ; i++){

      FixPosX[i] = kFALSE;

      FixPosY[i] = kFALSE;     

      FixAmp[i] = kFALSE;   

      Amp[i] = source[(int)(PosX[i]+0.5)][(int)(PosY[i]+0.5)];      //initial values of peaks amplitudes, input parameters         

      AmpXY[i] = 0;

   }

   //filling in the initial estimates of the input parameters

   TSpectrumFit2 *pfit=new TSpectrumFit2(nfound);

   pfit->SetFitParameters(xmin, xmax-1, ymin, ymax-1, 1000, 0.1, pfit->kFitOptimChiCounts, pfit->kFitAlphaHalving, pfit->kFitPower2, pfit->kFitTaylorOrderFirst);  

   pfit->SetPeakParameters(2, kFALSE, 2, kFALSE, 0, kTRUE, PosX, (Bool_t *) FixPosX, PosY, (Bool_t *) FixPosY, PosX, (Bool_t *) FixPosX, PosY, (Bool_t *) FixPosY, Amp, (Bool_t *) FixAmp, AmpXY, (Bool_t *) FixAmp, AmpXY, (Bool_t *) FixAmp);     

   pfit->SetBackgroundParameters(0, kFALSE, 0, kTRUE, 0, kTRUE);  

   pfit->FitStiefel(source);

    for (i = 0; i < nbinsx; i++){

     for (j = 0; j < nbinsy; j++){

                  search->SetBinContent(i + 1, j + 1,source[i][j]);

                 }

   }  

   search->Draw("SURF");

}