class TSpectrumFit: 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 is

 

 

 

 

 

 

 

 

where a represents vector of fitted parameters (positions p(j), amplitudes A(j), sigma, relative amplitudes T, S and slope B).

 

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

        Function parameters:
        source-pointer to the vector 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 [1]

•         it contains the symmetrical Gaussian as well as nonsymmetrical terms.

 

 

 

 

 

where T and S are relative amplitudes and B is slope.

 

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

Function:

void TSpectrumFit::FitAwmi(double *fSource)

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

 

Parameter:

        fSource-pointer to the vector of source spectrum                 

       

Member variables of the TSpectrumFit 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

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

 

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 FitAwmi.c:

Fig. 1 Original spectrum (black line) and fitted spectrum using AWMI algorithm (red line) and number of iteration steps = 1000. Positions of fitted peaks are denoted by markers

Script:

// Example to illustrate fitting function using AWMI algorithm.

// To execute this example, do

// root > .x FitAwmi.C

 

void FitAwmi() {

   Double_t a;

   Int_t i,nfound=0,bin;

   Int_t nbins = 256;

   Int_t xmin  = 0;

   Int_t xmax  = nbins;

   Double_t * source = new Double_t[nbins];

   Double_t * dest = new Double_t[nbins];  

   TH1F *h = new TH1F("h","Fitting using AWMI algorithm",nbins,xmin,xmax);

   TH1F *d = new TH1F("d","",nbins,xmin,xmax);     

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

   h=(TH1F*) f->Get("fit;1");  

   for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);     

   TCanvas *Fit1 = gROOT->GetListOfCanvases()->FindObject("Fit1");

   if (!Fit1) Fit1 = new TCanvas("Fit1","Fit1",10,10,1000,700);

   h->Draw("L");

   TSpectrum *s = new TSpectrum();

   //searching for candidate peaks positions

   nfound = s->SearchHighRes(source, dest, nbins, 2, 0.1, kFALSE, 10000, kFALSE, 0);

   Bool_t *FixPos = new Bool_t[nfound];

   Bool_t *FixAmp = new Bool_t[nfound];     

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

      FixPos[i] = kFALSE;

      FixAmp[i] = kFALSE;   

   }

   //filling in the initial estimates of the input parameters

   Double_t *PosX = new Double_t[nfound];        

   Double_t *PosY = new Double_t[nfound];

   PosX = s->GetPositionX();

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

                                                a=PosX[i];

        bin = 1 + Int_t(a + 0.5);

        PosY[i] = h->GetBinContent(bin);

   }  

   TSpectrumFit *pfit=new TSpectrumFit(nfound);

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

   pfit->SetPeakParameters(2, kFALSE, PosX, (Bool_t *) FixPos, PosY, (Bool_t *) FixAmp);  

   pfit->FitAwmi(source);

   Double_t *CalcPositions = new Double_t[nfound];     

   Double_t *CalcAmplitudes = new Double_t[nfound];        

   CalcPositions=pfit->GetPositions();

   CalcAmplitudes=pfit->GetAmplitudes();  

   for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);

   d->SetLineColor(kRed);  

   d->Draw("SAME L"); 

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

                                                a=CalcPositions[i];

        bin = 1 + Int_t(a + 0.5);               

        PosX[i] = d->GetBinCenter(bin);

        PosY[i] = d->GetBinContent(bin);

   }

   TPolyMarker * pm = (TPolyMarker*)h->GetListOfFunctions()->FindObject("TPolyMarker");

   if (pm) {

      h->GetListOfFunctions()->Remove(pm);

      delete pm;

   }

   pm = new TPolyMarker(nfound, PosX, PosY);

   h->GetListOfFunctions()->Add(pm);

   pm->SetMarkerStyle(23);

   pm->SetMarkerColor(kRed);

   pm->SetMarkerSize(1);  

}

        ONE-DIMENSIONAL FIT FUNCTION
        ALGORITHM WITH MATRIX INVERSION (STIEFEL-HESTENS METHOD)
        This function fits the source spectrum. The calling program should
        fill in input parameters
        The fitted parameters are written into
        output parameters and fitted data are written into
        source spectrum.

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

Stiefel fitting algorithm

 

Function:

void TSpectrumFit::FitStiefel(double *fSource)

 

This function fits the source spectrum using Stiefel-Hestens method [1] (see Awmi function).  The calling program should fill in input fitting parameters of the TSpectrumFit class using a set of TSpectrumFit 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 vector of source spectrum                 

       

Example – script FitStiefel.c:

Fig. 2 Original spectrum (black line) and fitted spectrum using Stiefel-Hestens method (red line) and number of iteration steps = 100. Positions of fitted peaks are denoted by markers

 

Script:

// Example to illustrate fitting function using Stiefel-Hestens method.

// To execute this example, do

// root > .x FitStiefel.C

 

void FitStiefel() {

   Double_t a;

   Int_t i,nfound=0,bin;

   Int_t nbins = 256;

   Int_t xmin  = 0;

   Int_t xmax  = nbins;

   Double_t * source = new Double_t[nbins];

   Double_t * dest = new Double_t[nbins];  

   TH1F *h = new TH1F("h","Fitting using AWMI algorithm",nbins,xmin,xmax);

   TH1F *d = new TH1F("d","",nbins,xmin,xmax);     

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

   h=(TH1F*) f->Get("fit;1");  

   for (i = 0; i < nbins; i++) source[i]=h->GetBinContent(i + 1);     

   TCanvas *Fit1 = gROOT->GetListOfCanvases()->FindObject("Fit1");

   if (!Fit1) Fit1 = new TCanvas("Fit1","Fit1",10,10,1000,700);

   h->Draw("L");

   TSpectrum *s = new TSpectrum();

   //searching for candidate peaks positions

   nfound = s->SearchHighRes(source, dest, nbins, 2, 0.1, kFALSE, 10000, kFALSE, 0);

   Bool_t *FixPos = new Bool_t[nfound];

   Bool_t *FixAmp = new Bool_t[nfound];     

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

      FixPos[i] = kFALSE;

      FixAmp[i] = kFALSE;   

   }

   //filling in the initial estimates of the input parameters

   Double_t *PosX = new Double_t[nfound];        

   Double_t *PosY = new Double_t[nfound];

   PosX = s->GetPositionX();

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

                                                a=PosX[i];

        bin = 1 + Int_t(a + 0.5);

        PosY[i] = h->GetBinContent(bin);

   }  

   TSpectrumFit *pfit=new TSpectrumFit(nfound);

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

   pfit->SetPeakParameters(2, kFALSE, PosX, (Bool_t *) FixPos, PosY, (Bool_t *) FixAmp);  

   pfit->FitStiefel(source);

   Double_t *CalcPositions = new Double_t[nfound];     

   Double_t *CalcAmplitudes = new Double_t[nfound];        

   CalcPositions=pfit->GetPositions();

   CalcAmplitudes=pfit->GetAmplitudes();  

   for (i = 0; i < nbins; i++) d->SetBinContent(i + 1,source[i]);

   d->SetLineColor(kRed);  

   d->Draw("SAME L"); 

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

                                                a=CalcPositions[i];

        bin = 1 + Int_t(a + 0.5);               

        PosX[i] = d->GetBinCenter(bin);

        PosY[i] = d->GetBinContent(bin);

   }

   TPolyMarker * pm = (TPolyMarker*)h->GetListOfFunctions()->FindObject("TPolyMarker");

   if (pm) {

      h->GetListOfFunctions()->Remove(pm);

      delete pm;

   }

   pm = new TPolyMarker(nfound, PosX, PosY);

   h->GetListOfFunctions()->Add(pm);

   pm->SetMarkerStyle(23);

   pm->SetMarkerColor(kRed);

   pm->SetMarkerSize(1);  

}