// @(#)root/spectrum:$Id$
// Author: Miroslav Morhac 25/09/2006
/////////////////////////////////////////////////////////////////////////////
// THIS CLASS CONTAINS ADVANCED TWO-DIMENSIONAL SPECTRA //
// FITTING FUNCTIONS //
// //
// These functions were written by: //
// Miroslav Morhac //
// Institute of Physics //
// Slovak Academy of Sciences //
// Dubravska cesta 9, 842 28 BRATISLAVA //
// SLOVAKIA //
// //
// email:fyzimiro@savba.sk, fax:+421 7 54772479 //
// //
// The original code in C has been repackaged as a C++ class by R.Brun //
// //
// The algorithms in this class have been published in the following //
// references: //
// [1] M. Morhac et al.: Efficient fitting algorithms applied to //
// analysis of coincidence gamma-ray spectra. Computer Physics //
// Communications, Vol 172/1 (2005) pp. 19-41. //
// //
// [2] M. Morhac et al.: Study of fitting algorithms applied to //
// simultaneous analysis of large number of peaks in gamma-ray spectra. //
// Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003. //
// //
// //
//____________________________________________________________________________
#include "TSpectrum2Fit.h"
#include "TMath.h"
ClassImp(TSpectrum2Fit)
//______________________________________________________________________________
TSpectrum2Fit::TSpectrum2Fit() :TNamed("Spectrum2Fit", "Miroslav Morhac peak fitter")
{
//default constructor
fNPeaks = 0;
fNumberIterations = 1;
fXmin = 0;
fXmax = 100;
fYmin = 0;
fYmax = 100;
fStatisticType = kFitOptimChiCounts;
fAlphaOptim = kFitAlphaHalving;
fPower = kFitPower2;
fFitTaylor = kFitTaylorOrderFirst;
fAlpha = 1;
fChi = 0;
fPositionInitX = 0;
fPositionCalcX = 0;
fPositionErrX = 0;
fPositionInitY = 0;
fPositionCalcY = 0;
fPositionErrY = 0;
fPositionInitX1 = 0;
fPositionCalcX1 = 0;
fPositionErrX1 = 0;
fPositionInitY1 = 0;
fPositionCalcY1 = 0;
fPositionErrY1 = 0;
fAmpInit = 0;
fAmpCalc = 0;
fAmpErr = 0;
fAmpInitX1 = 0;
fAmpCalcX1 = 0;
fAmpErrX1 = 0;
fAmpInitY1 = 0;
fAmpCalcY1 = 0;
fAmpErrY1 = 0;
fVolume = 0;
fVolumeErr = 0;
fSigmaInitX = 2;
fSigmaCalcX = 0;
fSigmaErrX = 0;
fSigmaInitY = 2;
fSigmaCalcY = 0;
fSigmaErrY = 0;
fRoInit = 0;
fRoCalc = 0;
fRoErr = 0;
fTxyInit = 0;
fTxyCalc = 0;
fTxyErr = 0;
fTxInit = 0;
fTxCalc = 0;
fTxErr = 0;
fTyInit = 0;
fTyCalc = 0;
fTyErr = 0;
fBxInit = 1;
fBxCalc = 0;
fBxErr = 0;
fByInit = 1;
fByCalc = 0;
fByErr = 0;
fSxyInit = 0;
fSxyCalc = 0;
fSxyErr = 0;
fSxInit = 0;
fSxCalc = 0;
fSxErr = 0;
fSyInit = 0;
fSyCalc = 0;
fSyErr = 0;
fA0Init = 0;
fA0Calc = 0;
fA0Err = 0;
fAxInit = 0;
fAxCalc = 0;
fAxErr = 0;
fAyInit = 0;
fAyCalc = 0;
fAyErr = 0;
fFixPositionX = 0;
fFixPositionY = 0;
fFixPositionX1 = 0;
fFixPositionY1 = 0;
fFixAmp = 0;
fFixAmpX1 = 0;
fFixAmpY1 = 0;
fFixSigmaX = false;
fFixSigmaY = false;
fFixRo = true;
fFixTxy = true;
fFixTx = true;
fFixTy = true;
fFixBx = true;
fFixBy = true;
fFixSxy = true;
fFixSx = true;
fFixSy = true;
fFixA0 = true;
fFixAx = true;
fFixAy = true;
}
//______________________________________________________________________________
TSpectrum2Fit::TSpectrum2Fit(Int_t numberPeaks) :TNamed("Spectrum2Fit", "Miroslav Morhac peak fitter")
{
//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.
//Begin_Html
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");
}
End_Html
Int_t i, i1, i2, j, k, shift =
7 * fNPeaks + 14, peak_vel, size, iter, pw,
regul_cycle, flag;
Double_t a, b, c, d = 0, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi =
0, pi, pmin = 0, chi_cel = 0, chi_er;
Double_t *working_space = new Double_t[5 * (7 * fNPeaks + 14)];
for (i = 0, j = 0; i < fNPeaks; i++) {
working_space[7 * i] = fAmpInit[i]; //vector parameter
if (fFixAmp[i] == false) {
working_space[shift + j] = fAmpInit[i]; //vector xk
j += 1;
}
working_space[7 * i + 1] = fPositionInitX[i]; //vector parameter
if (fFixPositionX[i] == false) {
working_space[shift + j] = fPositionInitX[i]; //vector xk
j += 1;
}
working_space[7 * i + 2] = fPositionInitY[i]; //vector parameter
if (fFixPositionY[i] == false) {
working_space[shift + j] = fPositionInitY[i]; //vector xk
j += 1;
}
working_space[7 * i + 3] = fAmpInitX1[i]; //vector parameter
if (fFixAmpX1[i] == false) {
working_space[shift + j] = fAmpInitX1[i]; //vector xk
j += 1;
}
working_space[7 * i + 4] = fAmpInitY1[i]; //vector parameter
if (fFixAmpY1[i] == false) {
working_space[shift + j] = fAmpInitY1[i]; //vector xk
j += 1;
}
working_space[7 * i + 5] = fPositionInitX1[i]; //vector parameter
if (fFixPositionX1[i] == false) {
working_space[shift + j] = fPositionInitX1[i]; //vector xk
j += 1;
}
working_space[7 * i + 6] = fPositionInitY1[i]; //vector parameter
if (fFixPositionY1[i] == false) {
working_space[shift + j] = fPositionInitY1[i]; //vector xk
j += 1;
}
}
peak_vel = 7 * i;
working_space[7 * i] = fSigmaInitX; //vector parameter
if (fFixSigmaX == false) {
working_space[shift + j] = fSigmaInitX; //vector xk
j += 1;
}
working_space[7 * i + 1] = fSigmaInitY; //vector parameter
if (fFixSigmaY == false) {
working_space[shift + j] = fSigmaInitY; //vector xk
j += 1;
}
working_space[7 * i + 2] = fRoInit; //vector parameter
if (fFixRo == false) {
working_space[shift + j] = fRoInit; //vector xk
j += 1;
}
working_space[7 * i + 3] = fA0Init; //vector parameter
if (fFixA0 == false) {
working_space[shift + j] = fA0Init; //vector xk
j += 1;
}
working_space[7 * i + 4] = fAxInit; //vector parameter
if (fFixAx == false) {
working_space[shift + j] = fAxInit; //vector xk
j += 1;
}
working_space[7 * i + 5] = fAyInit; //vector parameter
if (fFixAy == false) {
working_space[shift + j] = fAyInit; //vector xk
j += 1;
}
working_space[7 * i + 6] = fTxyInit; //vector parameter
if (fFixTxy == false) {
working_space[shift + j] = fTxyInit; //vector xk
j += 1;
}
working_space[7 * i + 7] = fSxyInit; //vector parameter
if (fFixSxy == false) {
working_space[shift + j] = fSxyInit; //vector xk
j += 1;
}
working_space[7 * i + 8] = fTxInit; //vector parameter
if (fFixTx == false) {
working_space[shift + j] = fTxInit; //vector xk
j += 1;
}
working_space[7 * i + 9] = fTyInit; //vector parameter
if (fFixTy == false) {
working_space[shift + j] = fTyInit; //vector xk
j += 1;
}
working_space[7 * i + 10] = fSxyInit; //vector parameter
if (fFixSx == false) {
working_space[shift + j] = fSxInit; //vector xk
j += 1;
}
working_space[7 * i + 11] = fSyInit; //vector parameter
if (fFixSy == false) {
working_space[shift + j] = fSyInit; //vector xk
j += 1;
}
working_space[7 * i + 12] = fBxInit; //vector parameter
if (fFixBx == false) {
working_space[shift + j] = fBxInit; //vector xk
j += 1;
}
working_space[7 * i + 13] = fByInit; //vector parameter
if (fFixBy == false) {
working_space[shift + j] = fByInit; //vector xk
j += 1;
}
size = j;
for (iter = 0; iter < fNumberIterations; iter++) {
for (j = 0; j < size; j++) {
working_space[2 * shift + j] = 0, working_space[3 * shift + j] = 0; //der,temp
}
//filling vectors
alpha = fAlpha;
chi_opt = 0, pw = fPower - 2;
for (i1 = fXmin; i1 <= fXmax; i1++) {
for (i2 = fYmin; i2 <= fYmax; i2++) {
yw = source[i1][i2];
ywm = yw;
f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 3],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 9],
working_space[peak_vel + 10],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (fStatisticType == kFitOptimMaxLikelihood) {
if (f > 0.00001)
chi_opt += yw * TMath::Log(f) - f;
}
else {
if (ywm != 0)
chi_opt += (yw - f) * (yw - f) / ywm;
}
if (fStatisticType == kFitOptimChiFuncValues) {
ywm = f;
if (f < 0.00001)
ywm = 0.00001;
}
else if (fStatisticType == kFitOptimMaxLikelihood) {
ywm = f;
if (f < 0.00001)
ywm = 0.00001;
}
else {
if (ywm == 0)
ywm = 1;
}
//calculation of gradient vector
for (j = 0, k = 0; j < fNPeaks; j++) {
if (fFixAmp[j] == false) {
a = Deramp2((Double_t) i1, (Double_t) i2,
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixPositionX[j] == false) {
a = Deri02((Double_t) i1, (Double_t) i2,
working_space[7 * j],
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (fFitTaylor == kFitTaylorOrderSecond)
d = Derderi02((Double_t) i1, (Double_t) i2,
working_space[7 * j],
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (TMath::Abs(a) > 0.00000001
&& fFitTaylor == kFitTaylorOrderSecond) {
d = d * TMath::Abs(yw - f) / (2 * a * ywm);
if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
&& a <= 0))
d = 0;
}
else
d = 0;
a = a + d;
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixPositionY[j] == false) {
a = Derj02((Double_t) i1, (Double_t) i2,
working_space[7 * j],
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (fFitTaylor == kFitTaylorOrderSecond)
d = Derderj02((Double_t) i1, (Double_t) i2,
working_space[7 * j],
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (TMath::Abs(a) > 0.00000001
&& fFitTaylor == kFitTaylorOrderSecond) {
d = d * TMath::Abs(yw - f) / (2 * a * ywm);
if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
&& a <= 0))
d = 0;
}
else
d = 0;
a = a + d;
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixAmpX1[j] == false) {
a = Derampx((Double_t) i1, working_space[7 * j + 5],
working_space[peak_vel],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixAmpY1[j] == false) {
a = Derampx((Double_t) i2, working_space[7 * j + 6],
working_space[peak_vel + 1],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 13]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixPositionX1[j] == false) {
a = Deri01((Double_t) i1, working_space[7 * j + 3],
working_space[7 * j + 5],
working_space[peak_vel],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12]);
if (fFitTaylor == kFitTaylorOrderSecond)
d = Derderi01((Double_t) i1, working_space[7 * j + 3],
working_space[7 * j + 5],
working_space[peak_vel]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (TMath::Abs(a) > 0.00000001
&& fFitTaylor == kFitTaylorOrderSecond) {
d = d * TMath::Abs(yw - f) / (2 * a * ywm);
if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
&& a <= 0))
d = 0;
}
else
d = 0;
a = a + d;
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //Der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //Der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixPositionY1[j] == false) {
a = Deri01((Double_t) i2, working_space[7 * j + 4],
working_space[7 * j + 6],
working_space[peak_vel + 1],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 13]);
if (fFitTaylor == kFitTaylorOrderSecond)
d = Derderi01((Double_t) i2, working_space[7 * j + 4],
working_space[7 * j + 6],
working_space[peak_vel + 1]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (TMath::Abs(a) > 0.00000001
&& fFitTaylor == kFitTaylorOrderSecond) {
d = d * TMath::Abs(yw - f) / (2 * a * ywm);
if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0
&& a <= 0))
d = 0;
}
else
d = 0;
a = a + d;
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
}
if (fFixSigmaX == false) {
a = Dersigmax(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (fFitTaylor == kFitTaylorOrderSecond)
d = Derdersigmax(fNPeaks, (Double_t) i1,
(Double_t) i2, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (TMath::Abs(a) > 0.00000001
&& fFitTaylor == kFitTaylorOrderSecond) {
d = d * TMath::Abs(yw - f) / (2 * a * ywm);
if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
d = 0;
}
else
d = 0;
a = a + d;
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixSigmaY == false) {
a = Dersigmay(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (fFitTaylor == kFitTaylorOrderSecond)
d = Derdersigmay(fNPeaks, (Double_t) i1,
(Double_t) i2, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (TMath::Abs(a) > 0.00000001
&& fFitTaylor == kFitTaylorOrderSecond) {
d = d * TMath::Abs(yw - f) / (2 * a * ywm);
if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
d = 0;
}
else
d = 0;
a = a + d;
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixRo == false) {
a = Derro(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (TMath::Abs(a) > 0.00000001
&& fFitTaylor == kFitTaylorOrderSecond) {
d = d * TMath::Abs(yw - f) / (2 * a * ywm);
if (((a + d) <= 0 && a >= 0) || ((a + d) >= 0 && a <= 0))
d = 0;
}
else
d = 0;
a = a + d;
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixA0 == false) {
a = 1.;
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixAx == false) {
a = i1;
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixAy == false) {
a = i2;
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixTxy == false) {
a = Dertxy(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixSxy == false) {
a = Dersxy(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixTx == false) {
a = Dertx(fNPeaks, (Double_t) i1, working_space,
working_space[peak_vel],
working_space[peak_vel + 12]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixTy == false) {
a = Derty(fNPeaks, (Double_t) i2, working_space,
working_space[peak_vel + 1],
working_space[peak_vel + 13]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixSx == false) {
a = Dersx(fNPeaks, (Double_t) i1, working_space,
working_space[peak_vel]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixSy == false) {
a = Dersy(fNPeaks, (Double_t) i2, working_space,
working_space[peak_vel + 1]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixBx == false) {
a = Derbx(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 6],
working_space[peak_vel + 8],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
if (fFixBy == false) {
a = Derby(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 6],
working_space[peak_vel + 8],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (ywm != 0) {
c = Ourpowl(a, pw);
if (fStatisticType == kFitOptimChiFuncValues) {
b = a * (yw * yw - f * f) / (ywm * ywm);
working_space[2 * shift + k] += b * c; //der
b = a * a * (4 * yw - 2 * f) / (ywm * ywm);
working_space[3 * shift + k] += b * c; //temp
}
else {
b = a * (yw - f) / ywm;
working_space[2 * shift + k] += b * c; //der
b = a * a / ywm;
working_space[3 * shift + k] += b * c; //temp
}
}
k += 1;
}
}
}
for (j = 0; j < size; j++) {
if (TMath::Abs(working_space[3 * shift + j]) > 0.000001)
working_space[2 * shift + j] = working_space[2 * shift + j] / TMath::Abs(working_space[3 * shift + j]); //der[j]=der[j]/temp[j]
else
working_space[2 * shift + j] = 0; //der[j]
}
//calculate chi_opt
chi2 = chi_opt;
chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
//calculate new parameters
regul_cycle = 0;
for (j = 0; j < size; j++) {
working_space[4 * shift + j] = working_space[shift + j]; //temp_xk[j]=xk[j]
}
do {
if (fAlphaOptim == kFitAlphaOptimal) {
if (fStatisticType != kFitOptimMaxLikelihood)
chi_min = 10000 * chi2;
else
chi_min = 0.1 * chi2;
flag = 0;
for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
for (j = 0; j < size; j++) {
working_space[shift + j] = working_space[4 * shift + j] + pi * alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pi*alpha*der[j]
}
for (i = 0, j = 0; i < fNPeaks; i++) {
if (fFixAmp[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
j += 1;
}
if (fFixPositionX[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
j += 1;
}
if (fFixPositionY[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
j += 1;
}
if (fFixAmpX1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
j += 1;
}
if (fFixAmpY1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
j += 1;
}
if (fFixPositionX1[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
j += 1;
}
if (fFixPositionY1[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
j += 1;
}
}
if (fFixSigmaX == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
j += 1;
}
if (fFixSigmaY == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixRo == false) {
if (working_space[shift + j] < -1) { //xk[j]
working_space[shift + j] = -1; //xk[j]
}
if (working_space[shift + j] > 1) { //xk[j]
working_space[shift + j] = 1; //xk[j]
}
working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixAx == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixAy == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixTxy == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
if (fFixSxy == false) {
working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
j += 1;
}
if (fFixTx == false) {
working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
j += 1;
}
if (fFixTy == false) {
working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
j += 1;
}
if (fFixSx == false) {
working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
j += 1;
}
if (fFixSy == false) {
working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
j += 1;
}
if (fFixBx == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
j += 1;
}
if (fFixBy == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
j += 1;
}
chi2 = 0;
for (i1 = fXmin; i1 <= fXmax; i1++) {
for (i2 = fYmin; i2 <= fYmax; i2++) {
yw = source[i1][i2];
ywm = yw;
f = Shape2(fNPeaks, (Double_t) i1,
(Double_t) i2, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 3],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 9],
working_space[peak_vel + 10],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (fStatisticType == kFitOptimChiFuncValues) {
ywm = f;
if (f < 0.00001)
ywm = 0.00001;
}
if (fStatisticType == kFitOptimMaxLikelihood) {
if (f > 0.00001)
chi2 += yw * TMath::Log(f) - f;
}
else {
if (ywm != 0)
chi2 += (yw - f) * (yw - f) / ywm;
}
}
}
if ((chi2 < chi_min
&& fStatisticType != kFitOptimMaxLikelihood)
|| (chi2 > chi_min
&& fStatisticType == kFitOptimMaxLikelihood)) {
pmin = pi, chi_min = chi2;
}
else
flag = 1;
if (pi == 0.1)
chi_min = chi2;
chi = chi_min;
}
if (pmin != 0.1) {
for (j = 0; j < size; j++) {
working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pmin*alpha*der[j]
}
for (i = 0, j = 0; i < fNPeaks; i++) {
if (fFixAmp[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
j += 1;
}
if (fFixPositionX[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
j += 1;
}
if (fFixPositionY[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
j += 1;
}
if (fFixAmpX1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
j += 1;
}
if (fFixAmpY1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
j += 1;
}
if (fFixPositionX1[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
j += 1;
}
if (fFixPositionY1[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
j += 1;
}
}
if (fFixSigmaX == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
j += 1;
}
if (fFixSigmaY == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixRo == false) {
if (working_space[shift + j] < -1) { //xk[j]
working_space[shift + j] = -1; //xk[j]
}
if (working_space[shift + j] > 1) { //xk[j]
working_space[shift + j] = 1; //xk[j]
}
working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixAx == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixAy == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixTxy == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
if (fFixSxy == false) {
working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
j += 1;
}
if (fFixTx == false) {
working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
j += 1;
}
if (fFixTy == false) {
working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
j += 1;
}
if (fFixSx == false) {
working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
j += 1;
}
if (fFixSy == false) {
working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
j += 1;
}
if (fFixBx == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
j += 1;
}
if (fFixBy == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
j += 1;
}
chi = chi_min;
}
}
else {
for (j = 0; j < size; j++) {
working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pi*alpha*der[j]
}
for (i = 0, j = 0; i < fNPeaks; i++) {
if (fFixAmp[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
j += 1;
}
if (fFixPositionX[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
j += 1;
}
if (fFixPositionY[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
j += 1;
}
if (fFixAmpX1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
j += 1;
}
if (fFixAmpY1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
j += 1;
}
if (fFixPositionX1[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
j += 1;
}
if (fFixPositionY1[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
j += 1;
}
}
if (fFixSigmaX == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
j += 1;
}
if (fFixSigmaY == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixRo == false) {
if (working_space[shift + j] < -1) { //xk[j]
working_space[shift + j] = -1; //xk[j]
}
if (working_space[shift + j] > 1) { //xk[j]
working_space[shift + j] = 1; //xk[j]
}
working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixAx == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixAy == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixTxy == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
if (fFixSxy == false) {
working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
j += 1;
}
if (fFixTx == false) {
working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
j += 1;
}
if (fFixTy == false) {
working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
j += 1;
}
if (fFixSx == false) {
working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
j += 1;
}
if (fFixSy == false) {
working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
j += 1;
}
if (fFixBx == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
j += 1;
}
if (fFixBy == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
j += 1;
}
chi = 0;
for (i1 = fXmin; i1 <= fXmax; i1++) {
for (i2 = fYmin; i2 <= fYmax; i2++) {
yw = source[i1][i2];
ywm = yw;
f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 3],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 9],
working_space[peak_vel + 10],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (fStatisticType == kFitOptimChiFuncValues) {
ywm = f;
if (f < 0.00001)
ywm = 0.00001;
}
if (fStatisticType == kFitOptimMaxLikelihood) {
if (f > 0.00001)
chi += yw * TMath::Log(f) - f;
}
else {
if (ywm != 0)
chi += (yw - f) * (yw - f) / ywm;
}
}
}
}
chi2 = chi;
chi = TMath::Sqrt(TMath::Abs(chi));
if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
alpha = alpha * chi_opt / (2 * chi);
else if (fAlphaOptim == kFitAlphaOptimal)
alpha = alpha / 10.0;
iter += 1;
regul_cycle += 1;
} while (((chi > chi_opt
&& fStatisticType != kFitOptimMaxLikelihood)
|| (chi < chi_opt
&& fStatisticType == kFitOptimMaxLikelihood))
&& regul_cycle < kFitNumRegulCycles);
for (j = 0; j < size; j++) {
working_space[4 * shift + j] = 0; //temp_xk[j]
working_space[2 * shift + j] = 0; //der[j]
}
for (i1 = fXmin, chi_cel = 0; i1 <= fXmax; i1++) {
for (i2 = fYmin; i2 <= fYmax; i2++) {
yw = source[i1][i2];
if (yw == 0)
yw = 1;
f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 3],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 9],
working_space[peak_vel + 10],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
chi_opt = (yw - f) * (yw - f) / yw;
chi_cel += (yw - f) * (yw - f) / yw;
//calculate gradient vector
for (j = 0, k = 0; j < fNPeaks; j++) {
if (fFixAmp[j] == false) {
a = Deramp2((Double_t) i1, (Double_t) i2,
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixPositionX[j] == false) {
a = Deri02((Double_t) i1, (Double_t) i2,
working_space[7 * j],
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixPositionY[j] == false) {
a = Derj02((Double_t) i1, (Double_t) i2,
working_space[7 * j],
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixAmpX1[j] == false) {
a = Derampx((Double_t) i1, working_space[7 * j + 5],
working_space[peak_vel],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixAmpY1[j] == false) {
a = Derampx((Double_t) i2, working_space[7 * j + 6],
working_space[peak_vel + 1],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 13]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixPositionX1[j] == false) {
a = Deri01((Double_t) i1, working_space[7 * j + 3],
working_space[7 * j + 5],
working_space[peak_vel],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixPositionY1[j] == false) {
a = Deri01((Double_t) i2, working_space[7 * j + 4],
working_space[7 * j + 6],
working_space[peak_vel + 1],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 13]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
}
if (fFixSigmaX == false) {
a = Dersigmax(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixSigmaY == false) {
a = Dersigmay(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixRo == false) {
a = Derro(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixA0 == false) {
a = 1.;
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixAx == false) {
a = i1;
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixAy == false) {
a = i2;
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixTxy == false) {
a = Dertxy(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixSxy == false) {
a = Dersxy(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixTx == false) {
a = Dertx(fNPeaks, (Double_t) i1, working_space,
working_space[peak_vel],
working_space[peak_vel + 12]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixTy == false) {
a = Derty(fNPeaks, (Double_t) i2, working_space,
working_space[peak_vel + 1],
working_space[peak_vel + 13]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixSx == false) {
a = Dersx(fNPeaks, (Double_t) i1, working_space,
working_space[peak_vel]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixSy == false) {
a = Dersy(fNPeaks, (Double_t) i2, working_space,
working_space[peak_vel + 1]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixBx == false) {
a = Derbx(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 6],
working_space[peak_vel + 8],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
if (fFixBy == false) {
a = Derby(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 6],
working_space[peak_vel + 8],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
c = Ourpowl(a, pw);
working_space[2 * shift + k] += chi_opt * c; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b * c; //temp_xk[k]
}
k += 1;
}
}
}
}
b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
chi_er = chi_cel / b;
for (i = 0, j = 0; i < fNPeaks; i++) {
fVolume[i] =
Volume(working_space[7 * i], working_space[peak_vel],
working_space[peak_vel + 1], working_space[peak_vel + 2]);
if (fVolume[i] > 0) {
c = 0;
if (fFixAmp[i] == false) {
a = Derpa2(working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
b = working_space[4 * shift + j]; //temp_xk[j]
if (b == 0)
b = 1;
else
b = 1 / b;
c = c + a * a * b;
}
if (fFixSigmaX == false) {
a = Derpsigmax(working_space[shift + j],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
b = working_space[4 * shift + peak_vel]; //temp_xk[j]
if (b == 0)
b = 1;
else
b = 1 / b;
c = c + a * a * b;
}
if (fFixSigmaY == false) {
a = Derpsigmay(working_space[shift + j],
working_space[peak_vel],
working_space[peak_vel + 2]);
b = working_space[4 * shift + peak_vel + 1]; //temp_xk[j]
if (b == 0)
b = 1;
else
b = 1 / b;
c = c + a * a * b;
}
if (fFixRo == false) {
a = Derpro(working_space[shift + j], working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
b = working_space[4 * shift + peak_vel + 2]; //temp_xk[j]
if (b == 0)
b = 1;
else
b = 1 / b;
c = c + a * a * b;
}
fVolumeErr[i] = TMath::Sqrt(TMath::Abs(chi_er * c));
}
else {
fVolumeErr[i] = 0;
}
if (fFixAmp[i] == false) {
fAmpCalc[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fAmpCalc[i] = fAmpInit[i];
fAmpErr[i] = 0;
}
if (fFixPositionX[i] == false) {
fPositionCalcX[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fPositionErrX[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fPositionCalcX[i] = fPositionInitX[i];
fPositionErrX[i] = 0;
}
if (fFixPositionY[i] == false) {
fPositionCalcY[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fPositionErrY[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fPositionCalcY[i] = fPositionInitY[i];
fPositionErrY[i] = 0;
}
if (fFixAmpX1[i] == false) {
fAmpCalcX1[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fAmpErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fAmpCalcX1[i] = fAmpInitX1[i];
fAmpErrX1[i] = 0;
}
if (fFixAmpY1[i] == false) {
fAmpCalcY1[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fAmpErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fAmpCalcY1[i] = fAmpInitY1[i];
fAmpErrY1[i] = 0;
}
if (fFixPositionX1[i] == false) {
fPositionCalcX1[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fPositionErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fPositionCalcX1[i] = fPositionInitX1[i];
fPositionErrX1[i] = 0;
}
if (fFixPositionY1[i] == false) {
fPositionCalcY1[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fPositionErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fPositionCalcY1[i] = fPositionInitY1[i];
fPositionErrY1[i] = 0;
}
}
if (fFixSigmaX == false) {
fSigmaCalcX = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSigmaErrX = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fSigmaCalcX = fSigmaInitX;
fSigmaErrX = 0;
}
if (fFixSigmaY == false) {
fSigmaCalcY = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSigmaErrY = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fSigmaCalcY = fSigmaInitY;
fSigmaErrY = 0;
}
if (fFixRo == false) {
fRoCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fRoErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fRoCalc = fRoInit;
fRoErr = 0;
}
if (fFixA0 == false) {
fA0Calc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fA0Calc = fA0Init;
fA0Err = 0;
}
if (fFixAx == false) {
fAxCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fAxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fAxCalc = fAxInit;
fAxErr = 0;
}
if (fFixAy == false) {
fAyCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fAyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fAyCalc = fAyInit;
fAyErr = 0;
}
if (fFixTxy == false) {
fTxyCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fTxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fTxyCalc = fTxyInit;
fTxyErr = 0;
}
if (fFixSxy == false) {
fSxyCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fSxyCalc = fSxyInit;
fSxyErr = 0;
}
if (fFixTx == false) {
fTxCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fTxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fTxCalc = fTxInit;
fTxErr = 0;
}
if (fFixTy == false) {
fTyCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fTyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fTyCalc = fTyInit;
fTyErr = 0;
}
if (fFixSx == false) {
fSxCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fSxCalc = fSxInit;
fSxErr = 0;
}
if (fFixSy == false) {
fSyCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fSyCalc = fSyInit;
fSyErr = 0;
}
if (fFixBx == false) {
fBxCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fBxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fBxCalc = fBxInit;
fBxErr = 0;
}
if (fFixBy == false) {
fByCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fByErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fByCalc = fByInit;
fByErr = 0;
}
b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
fChi = chi_cel / b;
for (i1 = fXmin; i1 <= fXmax; i1++) {
for (i2 = fYmin; i2 <= fYmax; i2++) {
f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 3],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 9],
working_space[peak_vel + 10],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
source[i1][i2] = f;
}
}
delete [] working_space;
return;
}
// ____________________________________________________________________________________________________________________________
void TSpectrum2Fit::FitStiefel(Float_t **source)
{
/////////////////////////////////////////////////////////////////////////////
// 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
//
/////////////////////////////////////////////////////////////////////////////
//
//Begin_Html
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");
}
End_Html
Int_t i, i1, i2, j, k, shift =
7 * fNPeaks + 14, peak_vel, size, iter, regul_cycle,
flag;
Double_t a, b, c, alpha, chi_opt, yw, ywm, f, chi2, chi_min, chi = 0
, pi, pmin = 0, chi_cel = 0, chi_er;
Double_t *working_space = new Double_t[5 * (7 * fNPeaks + 14)];
for (i = 0, j = 0; i < fNPeaks; i++) {
working_space[7 * i] = fAmpInit[i]; //vector parameter
if (fFixAmp[i] == false) {
working_space[shift + j] = fAmpInit[i]; //vector xk
j += 1;
}
working_space[7 * i + 1] = fPositionInitX[i]; //vector parameter
if (fFixPositionX[i] == false) {
working_space[shift + j] = fPositionInitX[i]; //vector xk
j += 1;
}
working_space[7 * i + 2] = fPositionInitY[i]; //vector parameter
if (fFixPositionY[i] == false) {
working_space[shift + j] = fPositionInitY[i]; //vector xk
j += 1;
}
working_space[7 * i + 3] = fAmpInitX1[i]; //vector parameter
if (fFixAmpX1[i] == false) {
working_space[shift + j] = fAmpInitX1[i]; //vector xk
j += 1;
}
working_space[7 * i + 4] = fAmpInitY1[i]; //vector parameter
if (fFixAmpY1[i] == false) {
working_space[shift + j] = fAmpInitY1[i]; //vector xk
j += 1;
}
working_space[7 * i + 5] = fPositionInitX1[i]; //vector parameter
if (fFixPositionX1[i] == false) {
working_space[shift + j] = fPositionInitX1[i]; //vector xk
j += 1;
}
working_space[7 * i + 6] = fPositionInitY1[i]; //vector parameter
if (fFixPositionY1[i] == false) {
working_space[shift + j] = fPositionInitY1[i]; //vector xk
j += 1;
}
}
peak_vel = 7 * i;
working_space[7 * i] = fSigmaInitX; //vector parameter
if (fFixSigmaX == false) {
working_space[shift + j] = fSigmaInitX; //vector xk
j += 1;
}
working_space[7 * i + 1] = fSigmaInitY; //vector parameter
if (fFixSigmaY == false) {
working_space[shift + j] = fSigmaInitY; //vector xk
j += 1;
}
working_space[7 * i + 2] = fRoInit; //vector parameter
if (fFixRo == false) {
working_space[shift + j] = fRoInit; //vector xk
j += 1;
}
working_space[7 * i + 3] = fA0Init; //vector parameter
if (fFixA0 == false) {
working_space[shift + j] = fA0Init; //vector xk
j += 1;
}
working_space[7 * i + 4] = fAxInit; //vector parameter
if (fFixAx == false) {
working_space[shift + j] = fAxInit; //vector xk
j += 1;
}
working_space[7 * i + 5] = fAyInit; //vector parameter
if (fFixAy == false) {
working_space[shift + j] = fAyInit; //vector xk
j += 1;
}
working_space[7 * i + 6] = fTxyInit; //vector parameter
if (fFixTxy == false) {
working_space[shift + j] = fTxyInit; //vector xk
j += 1;
}
working_space[7 * i + 7] = fSxyInit; //vector parameter
if (fFixSxy == false) {
working_space[shift + j] = fSxyInit; //vector xk
j += 1;
}
working_space[7 * i + 8] = fTxInit; //vector parameter
if (fFixTx == false) {
working_space[shift + j] = fTxInit; //vector xk
j += 1;
}
working_space[7 * i + 9] = fTyInit; //vector parameter
if (fFixTy == false) {
working_space[shift + j] = fTyInit; //vector xk
j += 1;
}
working_space[7 * i + 10] = fSxyInit; //vector parameter
if (fFixSx == false) {
working_space[shift + j] = fSxInit; //vector xk
j += 1;
}
working_space[7 * i + 11] = fSyInit; //vector parameter
if (fFixSy == false) {
working_space[shift + j] = fSyInit; //vector xk
j += 1;
}
working_space[7 * i + 12] = fBxInit; //vector parameter
if (fFixBx == false) {
working_space[shift + j] = fBxInit; //vector xk
j += 1;
}
working_space[7 * i + 13] = fByInit; //vector parameter
if (fFixBy == false) {
working_space[shift + j] = fByInit; //vector xk
j += 1;
}
size = j;
Double_t **working_matrix = new Double_t *[size];
for (i = 0; i < size; i++)
working_matrix[i] = new Double_t[size + 4];
for (iter = 0; iter < fNumberIterations; iter++) {
for (j = 0; j < size; j++) {
working_space[3 * shift + j] = 0; //temp
for (k = 0; k <= size; k++) {
working_matrix[j][k] = 0;
}
}
//filling working matrix
alpha = fAlpha;
chi_opt = 0;
for (i1 = fXmin; i1 <= fXmax; i1++) {
for (i2 = fYmin; i2 <= fYmax; i2++) {
//calculation of gradient vector
for (j = 0, k = 0; j < fNPeaks; j++) {
if (fFixAmp[j] == false) {
working_space[2 * shift + k] =
Deramp2((Double_t) i1, (Double_t) i2,
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
k += 1;
}
if (fFixPositionX[j] == false) {
working_space[2 * shift + k] =
Deri02((Double_t) i1, (Double_t) i2,
working_space[7 * j],
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
k += 1;
}
if (fFixPositionY[j] == false) {
working_space[2 * shift + k] =
Derj02((Double_t) i1, (Double_t) i2,
working_space[7 * j],
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
k += 1;
}
if (fFixAmpX1[j] == false) {
working_space[2 * shift + k] =
Derampx((Double_t) i1, working_space[7 * j + 5],
working_space[peak_vel],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12]);
k += 1;
}
if (fFixAmpY1[j] == false) {
working_space[2 * shift + k] =
Derampx((Double_t) i2, working_space[7 * j + 6],
working_space[peak_vel + 1],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 13]);
k += 1;
}
if (fFixPositionX1[j] == false) {
working_space[2 * shift + k] =
Deri01((Double_t) i1, working_space[7 * j + 3],
working_space[7 * j + 5],
working_space[peak_vel],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12]);
k += 1;
}
if (fFixPositionY1[j] == false) {
working_space[2 * shift + k] =
Deri01((Double_t) i2, working_space[7 * j + 4],
working_space[7 * j + 6],
working_space[peak_vel + 1],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 13]);
k += 1;
}
} if (fFixSigmaX == false) {
working_space[2 * shift + k] =
Dersigmax(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
k += 1;
}
if (fFixSigmaY == false) {
working_space[2 * shift + k] =
Dersigmay(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
k += 1;
}
if (fFixRo == false) {
working_space[2 * shift + k] =
Derro(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
k += 1;
}
if (fFixA0 == false) {
working_space[2 * shift + k] = 1.;
k += 1;
}
if (fFixAx == false) {
working_space[2 * shift + k] = i1;
k += 1;
}
if (fFixAy == false) {
working_space[2 * shift + k] = i2;
k += 1;
}
if (fFixTxy == false) {
working_space[2 * shift + k] =
Dertxy(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
k += 1;
}
if (fFixSxy == false) {
working_space[2 * shift + k] =
Dersxy(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1]);
k += 1;
}
if (fFixTx == false) {
working_space[2 * shift + k] =
Dertx(fNPeaks, (Double_t) i1, working_space,
working_space[peak_vel],
working_space[peak_vel + 12]);
k += 1;
}
if (fFixTy == false) {
working_space[2 * shift + k] =
Derty(fNPeaks, (Double_t) i2, working_space,
working_space[peak_vel + 1],
working_space[peak_vel + 13]);
k += 1;
}
if (fFixSx == false) {
working_space[2 * shift + k] =
Dersx(fNPeaks, (Double_t) i1, working_space,
working_space[peak_vel]);
k += 1;
}
if (fFixSy == false) {
working_space[2 * shift + k] =
Dersy(fNPeaks, (Double_t) i2, working_space,
working_space[peak_vel + 1]);
k += 1;
}
if (fFixBx == false) {
working_space[2 * shift + k] =
Derbx(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 6],
working_space[peak_vel + 8],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
k += 1;
}
if (fFixBy == false) {
working_space[2 * shift + k] =
Derby(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 6],
working_space[peak_vel + 8],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
k += 1;
}
yw = source[i1][i2];
ywm = yw;
f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 3],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 9],
working_space[peak_vel + 10],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (fStatisticType == kFitOptimMaxLikelihood) {
if (f > 0.00001)
chi_opt += yw * TMath::Log(f) - f;
}
else {
if (ywm != 0)
chi_opt += (yw - f) * (yw - f) / ywm;
}
if (fStatisticType == kFitOptimChiFuncValues) {
ywm = f;
if (f < 0.00001)
ywm = 0.00001;
}
else if (fStatisticType == kFitOptimMaxLikelihood) {
ywm = f;
if (f < 0.00001)
ywm = 0.00001;
}
else {
if (ywm == 0)
ywm = 1;
}
for (j = 0; j < size; j++) {
for (k = 0; k < size; k++) {
b = working_space[2 * shift +
j] * working_space[2 * shift +
k] / ywm;
if (fStatisticType == kFitOptimChiFuncValues)
b = b * (4 * yw - 2 * f) / ywm;
working_matrix[j][k] += b;
if (j == k)
working_space[3 * shift + j] += b;
}
}
if (fStatisticType == kFitOptimChiFuncValues)
b = (f * f - yw * yw) / (ywm * ywm);
else
b = (f - yw) / ywm;
for (j = 0; j < size; j++) {
working_matrix[j][size] -=
b * working_space[2 * shift + j];
}
}
}
for (i = 0; i < size; i++) {
working_matrix[i][size + 1] = 0; //xk
}
StiefelInversion(working_matrix, size);
for (i = 0; i < size; i++) {
working_space[2 * shift + i] = working_matrix[i][size + 1]; //der
}
//calculate chi_opt
chi2 = chi_opt;
chi_opt = TMath::Sqrt(TMath::Abs(chi_opt));
//calculate new parameters
regul_cycle = 0;
for (j = 0; j < size; j++) {
working_space[4 * shift + j] = working_space[shift + j]; //temp_xk[j]=xk[j]
}
do {
if (fAlphaOptim == kFitAlphaOptimal) {
if (fStatisticType != kFitOptimMaxLikelihood)
chi_min = 10000 * chi2;
else
chi_min = 0.1 * chi2;
flag = 0;
for (pi = 0.1; flag == 0 && pi <= 100; pi += 0.1) {
for (j = 0; j < size; j++) {
working_space[shift + j] = working_space[4 * shift + j] + pi * alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pi*alpha*der[j]
}
for (i = 0, j = 0; i < fNPeaks; i++) {
if (fFixAmp[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
j += 1;
}
if (fFixPositionX[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
j += 1;
}
if (fFixPositionY[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
j += 1;
}
if (fFixAmpX1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
j += 1;
}
if (fFixAmpY1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
j += 1;
}
if (fFixPositionX1[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
j += 1;
}
if (fFixPositionY1[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
j += 1;
}
}
if (fFixSigmaX == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
j += 1;
}
if (fFixSigmaY == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixRo == false) {
if (working_space[shift + j] < -1) { //xk[j]
working_space[shift + j] = -1; //xk[j]
}
if (working_space[shift + j] > 1) { //xk[j]
working_space[shift + j] = 1; //xk[j]
}
working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixAx == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixAy == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixTxy == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
if (fFixSxy == false) {
working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
j += 1;
}
if (fFixTx == false) {
working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
j += 1;
}
if (fFixTy == false) {
working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
j += 1;
}
if (fFixSx == false) {
working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
j += 1;
}
if (fFixSy == false) {
working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
j += 1;
}
if (fFixBx == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
j += 1;
}
if (fFixBy == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
j += 1;
}
chi2 = 0;
for (i1 = fXmin; i1 <= fXmax; i1++) {
for (i2 = fYmin; i2 <= fYmax; i2++) {
yw = source[i1][i2];
ywm = yw;
f = Shape2(fNPeaks, (Double_t) i1,
(Double_t) i2, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 3],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 9],
working_space[peak_vel + 10],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (fStatisticType == kFitOptimChiFuncValues) {
ywm = f;
if (f < 0.00001)
ywm = 0.00001;
}
if (fStatisticType == kFitOptimMaxLikelihood) {
if (f > 0.00001)
chi2 += yw * TMath::Log(f) - f;
}
else {
if (ywm != 0)
chi2 += (yw - f) * (yw - f) / ywm;
}
}
}
if ((chi2 < chi_min
&& fStatisticType != kFitOptimMaxLikelihood)
|| (chi2 > chi_min
&& fStatisticType == kFitOptimMaxLikelihood)) {
pmin = pi, chi_min = chi2;
}
else
flag = 1;
if (pi == 0.1)
chi_min = chi2;
chi = chi_min;
}
if (pmin != 0.1) {
for (j = 0; j < size; j++) {
working_space[shift + j] = working_space[4 * shift + j] + pmin * alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pmin*alpha*der[j]
}
for (i = 0, j = 0; i < fNPeaks; i++) {
if (fFixAmp[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
j += 1;
}
if (fFixPositionX[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
j += 1;
}
if (fFixPositionY[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
j += 1;
}
if (fFixAmpX1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
j += 1;
}
if (fFixAmpY1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
j += 1;
}
if (fFixPositionX1[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
j += 1;
}
if (fFixPositionY1[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
j += 1;
}
}
if (fFixSigmaX == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
j += 1;
}
if (fFixSigmaY == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixRo == false) {
if (working_space[shift + j] < -1) { //xk[j]
working_space[shift + j] = -1; //xk[j]
}
if (working_space[shift + j] > 1) { //xk[j]
working_space[shift + j] = 1; //xk[j]
}
working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixAx == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixAy == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixTxy == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
if (fFixSxy == false) {
working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
j += 1;
}
if (fFixTx == false) {
working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
j += 1;
}
if (fFixTy == false) {
working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
j += 1;
}
if (fFixSx == false) {
working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
j += 1;
}
if (fFixSy == false) {
working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
j += 1;
}
if (fFixBx == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
j += 1;
}
if (fFixBy == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
j += 1;
}
chi = chi_min;
}
}
else {
for (j = 0; j < size; j++) {
working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+pi*alpha*der[j]
}
for (i = 0, j = 0; i < fNPeaks; i++) {
if (fFixAmp[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i] = working_space[shift + j]; //parameter[7*i]=xk[j]
j += 1;
}
if (fFixPositionX[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 1] = working_space[shift + j]; //parameter[7*i+1]=xk[j]
j += 1;
}
if (fFixPositionY[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 2] = working_space[shift + j]; //parameter[7*i+2]=xk[j]
j += 1;
}
if (fFixAmpX1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 3] = working_space[shift + j]; //parameter[7*i+3]=xk[j]
j += 1;
}
if (fFixAmpY1[i] == false) {
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = 0; //xk[j]
working_space[7 * i + 4] = working_space[shift + j]; //parameter[7*i+4]=xk[j]
j += 1;
}
if (fFixPositionX1[i] == false) {
if (working_space[shift + j] < fXmin) //xk[j]
working_space[shift + j] = fXmin; //xk[j]
if (working_space[shift + j] > fXmax) //xk[j]
working_space[shift + j] = fXmax; //xk[j]
working_space[7 * i + 5] = working_space[shift + j]; //parameter[7*i+5]=xk[j]
j += 1;
}
if (fFixPositionY1[i] == false) {
if (working_space[shift + j] < fYmin) //xk[j]
working_space[shift + j] = fYmin; //xk[j]
if (working_space[shift + j] > fYmax) //xk[j]
working_space[shift + j] = fYmax; //xk[j]
working_space[7 * i + 6] = working_space[shift + j]; //parameter[7*i+6]=xk[j]
j += 1;
}
}
if (fFixSigmaX == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel] = working_space[shift + j]; //parameter[peak_vel]=xk[j]
j += 1;
}
if (fFixSigmaY == false) {
if (working_space[shift + j] < 0.001) { //xk[j]
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixRo == false) {
if (working_space[shift + j] < -1) { //xk[j]
working_space[shift + j] = -1; //xk[j]
}
if (working_space[shift + j] > 1) { //xk[j]
working_space[shift + j] = 1; //xk[j]
}
working_space[peak_vel + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixAx == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixAy == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixTxy == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
if (fFixSxy == false) {
working_space[peak_vel + 7] = working_space[shift + j]; //parameter[peak_vel+7]=xk[j]
j += 1;
}
if (fFixTx == false) {
working_space[peak_vel + 8] = working_space[shift + j]; //parameter[peak_vel+8]=xk[j]
j += 1;
}
if (fFixTy == false) {
working_space[peak_vel + 9] = working_space[shift + j]; //parameter[peak_vel+9]=xk[j]
j += 1;
}
if (fFixSx == false) {
working_space[peak_vel + 10] = working_space[shift + j]; //parameter[peak_vel+10]=xk[j]
j += 1;
}
if (fFixSy == false) {
working_space[peak_vel + 11] = working_space[shift + j]; //parameter[peak_vel+11]=xk[j]
j += 1;
}
if (fFixBx == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 12] = working_space[shift + j]; //parameter[peak_vel+12]=xk[j]
j += 1;
}
if (fFixBy == false) {
if (TMath::Abs(working_space[shift + j]) < 0.001) { //xk[j]
if (working_space[shift + j] < 0) //xk[j]
working_space[shift + j] = -0.001; //xk[j]
else
working_space[shift + j] = 0.001; //xk[j]
}
working_space[peak_vel + 13] = working_space[shift + j]; //parameter[peak_vel+13]=xk[j]
j += 1;
}
chi = 0;
for (i1 = fXmin; i1 <= fXmax; i1++) {
for (i2 = fYmin; i2 <= fYmax; i2++) {
yw = source[i1][i2];
ywm = yw;
f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 3],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 9],
working_space[peak_vel + 10],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (fStatisticType == kFitOptimChiFuncValues) {
ywm = f;
if (f < 0.00001)
ywm = 0.00001;
}
if (fStatisticType == kFitOptimMaxLikelihood) {
if (f > 0.00001)
chi += yw * TMath::Log(f) - f;
}
else {
if (ywm != 0)
chi += (yw - f) * (yw - f) / ywm;
}
}
}
}
chi2 = chi;
chi = TMath::Sqrt(TMath::Abs(chi));
if (fAlphaOptim == kFitAlphaHalving && chi > 1E-6)
alpha = alpha * chi_opt / (2 * chi);
else if (fAlphaOptim == kFitAlphaOptimal)
alpha = alpha / 10.0;
iter += 1;
regul_cycle += 1;
} while (((chi > chi_opt
&& fStatisticType != kFitOptimMaxLikelihood)
|| (chi < chi_opt
&& fStatisticType == kFitOptimMaxLikelihood))
&& regul_cycle < kFitNumRegulCycles);
for (j = 0; j < size; j++) {
working_space[4 * shift + j] = 0; //temp_xk[j]
working_space[2 * shift + j] = 0; //der[j]
}
for (i1 = fXmin, chi_cel = 0; i1 <= fXmax; i1++) {
for (i2 = fYmin; i2 <= fYmax; i2++) {
yw = source[i1][i2];
if (yw == 0)
yw = 1;
f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 3],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 9],
working_space[peak_vel + 10],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
chi_opt = (yw - f) * (yw - f) / yw;
chi_cel += (yw - f) * (yw - f) / yw;
//calculate gradient vector
for (j = 0, k = 0; j < fNPeaks; j++) {
if (fFixAmp[j] == false) {
a = Deramp2((Double_t) i1, (Double_t) i2,
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixPositionX[j] == false) {
a = Deri02((Double_t) i1, (Double_t) i2,
working_space[7 * j],
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixPositionY[j] == false) {
a = Derj02((Double_t) i1, (Double_t) i2,
working_space[7 * j],
working_space[7 * j + 1],
working_space[7 * j + 2],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixAmpX1[j] == false) {
a = Derampx((Double_t) i1, working_space[7 * j + 5],
working_space[peak_vel],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixAmpY1[j] == false) {
a = Derampx((Double_t) i2, working_space[7 * j + 6],
working_space[peak_vel + 1],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 13]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixPositionX1[j] == false) {
a = Deri01((Double_t) i1, working_space[7 * j + 3],
working_space[7 * j + 5],
working_space[peak_vel],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixPositionY1[j] == false) {
a = Deri01((Double_t) i2, working_space[7 * j + 4],
working_space[7 * j + 6],
working_space[peak_vel + 1],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 13]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
}
if (fFixSigmaX == false) {
a = Dersigmax(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 10],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixSigmaY == false) {
a = Dersigmay(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 9],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixRo == false) {
a = Derro(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixA0 == false) {
a = 1.;
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixAx == false) {
a = i1;
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixAy == false) {
a = i2;
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixTxy == false) {
a = Dertxy(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixSxy == false) {
a = Dersxy(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixTx == false) {
a = Dertx(fNPeaks, (Double_t) i1, working_space,
working_space[peak_vel],
working_space[peak_vel + 12]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixTy == false) {
a = Derty(fNPeaks, (Double_t) i2, working_space,
working_space[peak_vel + 1],
working_space[peak_vel + 13]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixSx == false) {
a = Dersx(fNPeaks, (Double_t) i1, working_space,
working_space[peak_vel]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixSy == false) {
a = Dersy(fNPeaks, (Double_t) i2, working_space,
working_space[peak_vel + 1]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixBx == false) {
a = Derbx(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 6],
working_space[peak_vel + 8],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
if (fFixBy == false) {
a = Derby(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 6],
working_space[peak_vel + 8],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
if (yw != 0) {
working_space[2 * shift + k] += chi_opt; //der[k]
b = a * a / yw;
working_space[4 * shift + k] += b; //temp_xk[k]
}
k += 1;
}
}
}
}
b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
chi_er = chi_cel / b;
for (i = 0, j = 0; i < fNPeaks; i++) {
fVolume[i] =
Volume(working_space[7 * i], working_space[peak_vel],
working_space[peak_vel + 1], working_space[peak_vel + 2]);
if (fVolume[i] > 0) {
c = 0;
if (fFixAmp[i] == false) {
a = Derpa2(working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
b = working_space[4 * shift + j]; //temp_xk[j]
if (b == 0)
b = 1;
else
b = 1 / b;
c = c + a * a * b;
}
if (fFixSigmaX == false) {
a = Derpsigmax(working_space[shift + j],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
b = working_space[4 * shift + peak_vel]; //temp_xk[j]
if (b == 0)
b = 1;
else
b = 1 / b;
c = c + a * a * b;
}
if (fFixSigmaY == false) {
a = Derpsigmay(working_space[shift + j],
working_space[peak_vel],
working_space[peak_vel + 2]);
b = working_space[4 * shift + peak_vel + 1]; //temp_xk[j]
if (b == 0)
b = 1;
else
b = 1 / b;
c = c + a * a * b;
}
if (fFixRo == false) {
a = Derpro(working_space[shift + j], working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2]);
b = working_space[4 * shift + peak_vel + 2]; //temp_xk[j]
if (b == 0)
b = 1;
else
b = 1 / b;
c = c + a * a * b;
}
fVolumeErr[i] = TMath::Sqrt(TMath::Abs(chi_er * c));
}
else {
fVolumeErr[i] = 0;
}
if (fFixAmp[i] == false) {
fAmpCalc[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fAmpErr[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fAmpCalc[i] = fAmpInit[i];
fAmpErr[i] = 0;
}
if (fFixPositionX[i] == false) {
fPositionCalcX[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fPositionErrX[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fPositionCalcX[i] = fPositionInitX[i];
fPositionErrX[i] = 0;
}
if (fFixPositionY[i] == false) {
fPositionCalcY[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fPositionErrY[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fPositionCalcY[i] = fPositionInitY[i];
fPositionErrY[i] = 0;
}
if (fFixAmpX1[i] == false) {
fAmpCalcX1[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fAmpErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fAmpCalcX1[i] = fAmpInitX1[i];
fAmpErrX1[i] = 0;
}
if (fFixAmpY1[i] == false) {
fAmpCalcY1[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fAmpErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fAmpCalcY1[i] = fAmpInitY1[i];
fAmpErrY1[i] = 0;
}
if (fFixPositionX1[i] == false) {
fPositionCalcX1[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fPositionErrX1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fPositionCalcX1[i] = fPositionInitX1[i];
fPositionErrX1[i] = 0;
}
if (fFixPositionY1[i] == false) {
fPositionCalcY1[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0)
fPositionErrY1[i] = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fPositionCalcY1[i] = fPositionInitY1[i];
fPositionErrY1[i] = 0;
}
}
if (fFixSigmaX == false) {
fSigmaCalcX = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSigmaErrX = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fSigmaCalcX = fSigmaInitX;
fSigmaErrX = 0;
}
if (fFixSigmaY == false) {
fSigmaCalcY = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSigmaErrY = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fSigmaCalcY = fSigmaInitY;
fSigmaErrY = 0;
}
if (fFixRo == false) {
fRoCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fRoErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fRoCalc = fRoInit;
fRoErr = 0;
}
if (fFixA0 == false) {
fA0Calc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fA0Err = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fA0Calc = fA0Init;
fA0Err = 0;
}
if (fFixAx == false) {
fAxCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fAxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fAxCalc = fAxInit;
fAxErr = 0;
}
if (fFixAy == false) {
fAyCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fAyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fAyCalc = fAyInit;
fAyErr = 0;
}
if (fFixTxy == false) {
fTxyCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fTxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fTxyCalc = fTxyInit;
fTxyErr = 0;
}
if (fFixSxy == false) {
fSxyCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSxyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fSxyCalc = fSxyInit;
fSxyErr = 0;
}
if (fFixTx == false) {
fTxCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fTxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fTxCalc = fTxInit;
fTxErr = 0;
}
if (fFixTy == false) {
fTyCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fTyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fTyCalc = fTyInit;
fTyErr = 0;
}
if (fFixSx == false) {
fSxCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fSxCalc = fSxInit;
fSxErr = 0;
}
if (fFixSy == false) {
fSyCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSyErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fSyCalc = fSyInit;
fSyErr = 0;
}
if (fFixBx == false) {
fBxCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fBxErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fBxCalc = fBxInit;
fBxErr = 0;
}
if (fFixBy == false) {
fByCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fByErr = TMath::Sqrt(TMath::Abs(working_space[2 * shift + j])) / TMath::Sqrt(TMath::Abs(working_space[3 * shift + j])); //der[j]/temp[j]
j += 1;
}
else {
fByCalc = fByInit;
fByErr = 0;
}
b = (fXmax - fXmin + 1) * (fYmax - fYmin + 1) - size;
fChi = chi_cel / b;
for (i1 = fXmin; i1 <= fXmax; i1++) {
for (i2 = fYmin; i2 <= fYmax; i2++) {
f = Shape2(fNPeaks, (Double_t) i1, (Double_t) i2,
working_space, working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 2],
working_space[peak_vel + 3],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6],
working_space[peak_vel + 7],
working_space[peak_vel + 8],
working_space[peak_vel + 9],
working_space[peak_vel + 10],
working_space[peak_vel + 11],
working_space[peak_vel + 12],
working_space[peak_vel + 13]);
source[i1][i2] = f;
}
}
for (i = 0; i < size; i++) delete [] working_matrix[i];
delete [] working_matrix;
delete [] working_space;
return;
}
void TSpectrum2Fit::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)
{
//////////////////////////////////////////////////////////////////////////////
// SETTER FUNCTION
//
// This function sets the following fitting parameters:
// -xmin, xmax, ymin, ymax - fitting region
// -numberIterations - # of desired iterations in the fit
// -alpha - convergence coefficient, it should be positive number and <=1, for details see references
// -statisticType - 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
// -alphaOptim - optimization of convergence algorithm, possible values kFitAlphaHalving, kFitAlphaOptimal
// -power - possible values kFitPower2,4,6,8,10,12, for details see references. It applies only for Awmi fitting function.
// -fitTaylor - order of Taylor expansion, possible values kFitTaylorOrderFirst, kFitTaylorOrderSecond. It applies only for Awmi fitting function.
//////////////////////////////////////////////////////////////////////////////
if(xmin<0 || xmax <= xmin || ymin<0 || ymax <= ymin){
Error("SetFitParameters", "Wrong range");
return;
}
if (numberIterations <= 0){
Error("SetFitParameters","Invalid number of iterations, must be positive");
return;
}
if (alpha <= 0 || alpha > 1){
Error ("SetFitParameters","Invalid step coefficient alpha, must be > than 0 and <=1");
return;
}
if (statisticType != kFitOptimChiCounts
&& statisticType != kFitOptimChiFuncValues
&& statisticType != kFitOptimMaxLikelihood){
Error("SetFitParameters","Wrong type of statistic");
return;
}
if (alphaOptim != kFitAlphaHalving
&& alphaOptim != kFitAlphaOptimal){
Error("SetFitParameters","Wrong optimization algorithm");
return;
}
if (power != kFitPower2 && power != kFitPower4
&& power != kFitPower6 && power != kFitPower8
&& power != kFitPower10 && power != kFitPower12){
Error("SetFitParameters","Wrong power");
return;
}
if (fitTaylor != kFitTaylorOrderFirst
&& fitTaylor != kFitTaylorOrderSecond){
Error("SetFitParameters","Wrong order of Taylor development");
return;
}
fXmin=xmin,fXmax=xmax,fYmin=ymin,fYmax=ymax,fNumberIterations=numberIterations,fAlpha=alpha,fStatisticType=statisticType,fAlphaOptim=alphaOptim,fPower=power,fFitTaylor=fitTaylor;
}
//______________________________________________________________________________
void TSpectrum2Fit::SetPeakParameters(Double_t sigmaX, Bool_t fixSigmaX, Double_t sigmaY, Bool_t fixSigmaY, Double_t ro, Bool_t fixRo, const Float_t *positionInitX, const Bool_t *fixPositionX, const Float_t *positionInitY, const Bool_t *fixPositionY, const Float_t *positionInitX1, const Bool_t *fixPositionX1, const Float_t *positionInitY1, const Bool_t *fixPositionY1, const Float_t *ampInit, const Bool_t *fixAmp, const Float_t *ampInitX1, const Bool_t *fixAmpX1, const Float_t *ampInitY1, const Bool_t *fixAmpY1)
{
//////////////////////////////////////////////////////////////////////////////
// SETTER FUNCTION
//
// This function sets the following fitting parameters of peaks:
// -sigmaX - initial value of sigma x parameter
// -fixSigmaX - logical value of sigma x parameter, which allows to fix the parameter (not to fit)
// -sigmaY - initial value of sigma y parameter
// -fixSigmaY - logical value of sigma y parameter, which allows to fix the parameter (not to fit)
// -ro - initial value of ro parameter (correlation coefficient)
// -fixRo - logical value of ro parameter, which allows to fix the parameter (not to fit)
// -positionInitX - aray of initial values of peaks x positions
// -fixPositionX - array of logical values which allow to fix appropriate x positions (not fit). However they are present in the estimated functional.
// -positionInitY - aray of initial values of peaks y positions
// -fixPositionY - array of logical values which allow to fix appropriate y positions (not fit). However they are present in the estimated functional.
// -ampInit - aray of initial values of 2D peaks amplitudes
// -fixAmp - aray of logical values which allow to fix appropriate amplitudes of 2D peaks (not fit). However they are present in the estimated functional
// -ampInitX1 - aray of initial values of amplitudes of 1D ridges in x direction
// -fixAmpX1 - aray 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
// -ampInitY1 - aray of initial values of amplitudes of 1D ridges in y direction
// -fixAmpY1 - aray 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
//////////////////////////////////////////////////////////////////////////////
if (sigmaX <= 0 || sigmaY <= 0){
Error ("SetPeakParameters","Invalid sigma, must be > than 0");
return;
}
if (ro < -1 || ro > 1){
Error ("SetPeakParameters","Invalid ro, must be from region <-1,1>");
return;
}
Int_t i;
for(i=0; i < fNPeaks; i++){
if(positionInitX[i] < fXmin || positionInitX[i] > fXmax){
Error ("SetPeakParameters","Invalid peak position, must be in the range fXmin, fXmax");
return;
}
if(positionInitY[i] < fYmin || positionInitY[i] > fYmax){
Error ("SetPeakParameters","Invalid peak position, must be in the range fYmin, fYmax");
return;
}
if(positionInitX1[i] < fXmin || positionInitX1[i] > fXmax){
Error ("SetPeakParameters","Invalid ridge position, must be in the range fXmin, fXmax");
return;
}
if(positionInitY1[i] < fYmin || positionInitY1[i] > fYmax){
Error ("SetPeakParameters","Invalid ridge position, must be in the range fYmin, fYmax");
return;
}
if(ampInit[i] < 0){
Error ("SetPeakParameters","Invalid peak amplitude, must be > than 0");
return;
}
if(ampInitX1[i] < 0){
Error ("SetPeakParameters","Invalid x ridge amplitude, must be > than 0");
return;
}
if(ampInitY1[i] < 0){
Error ("SetPeakParameters","Invalid y ridge amplitude, must be > than 0");
return;
}
}
fSigmaInitX = sigmaX, fFixSigmaX = fixSigmaX, fSigmaInitY = sigmaY, fFixSigmaY = fixSigmaY, fRoInit = ro, fFixRo = fixRo;
for(i=0; i < fNPeaks; i++){
fPositionInitX[i] = (Double_t) positionInitX[i];
fFixPositionX[i] = fixPositionX[i];
fPositionInitY[i] = (Double_t) positionInitY[i];
fFixPositionY[i] = fixPositionY[i];
fPositionInitX1[i] = (Double_t) positionInitX1[i];
fFixPositionX1[i] = fixPositionX1[i];
fPositionInitY1[i] = (Double_t) positionInitY1[i];
fFixPositionY1[i] = fixPositionY1[i];
fAmpInit[i] = (Double_t) ampInit[i];
fFixAmp[i] = fixAmp[i];
fAmpInitX1[i] = (Double_t) ampInitX1[i];
fFixAmpX1[i] = fixAmpX1[i];
fAmpInitY1[i] = (Double_t) ampInitY1[i];
fFixAmpY1[i] = fixAmpY1[i];
}
}
//______________________________________________________________________________
void TSpectrum2Fit::SetBackgroundParameters(Double_t a0Init, Bool_t fixA0, Double_t axInit, Bool_t fixAx, Double_t ayInit, Bool_t fixAy)
{
//////////////////////////////////////////////////////////////////////////////
// SETTER FUNCTION
//
// This function sets the following fitting parameters of background:
// -a0Init - initial value of a0 parameter (backgroud is estimated as a0+ax*x+ay*y)
// -fixA0 - logical value of a0 parameter, which allows to fix the parameter (not to fit)
// -axInit - initial value of ax parameter
// -fixAx - logical value of ax parameter, which allows to fix the parameter (not to fit)
// -ayInit - initial value of ay parameter
// -fixAy - logical value of ay parameter, which allows to fix the parameter (not to fit)
//////////////////////////////////////////////////////////////////////////////
fA0Init = a0Init;
fFixA0 = fixA0;
fAxInit = axInit;
fFixAx = fixAx;
fAyInit = ayInit;
fFixAy = fixAy;
}
//______________________________________________________________________________
void TSpectrum2Fit::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)
{
//////////////////////////////////////////////////////////////////////////////
// SETTER FUNCTION
//
// This function sets the following fitting parameters of tails of peaks
// -tInitXY - initial value of txy parameter
// -fixTxy - logical value of txy parameter, which allows to fix the parameter (not to fit)
// -tInitX - initial value of tx parameter
// -fixTx - logical value of tx parameter, which allows to fix the parameter (not to fit)
// -tInitY - initial value of ty parameter
// -fixTy - logical value of ty parameter, which allows to fix the parameter (not to fit)
// -bInitX - initial value of bx parameter
// -fixBx - logical value of bx parameter, which allows to fix the parameter (not to fit)
// -bInitY - initial value of by parameter
// -fixBy - logical value of by parameter, which allows to fix the parameter (not to fit)
// -sInitXY - initial value of sxy parameter
// -fixSxy - logical value of sxy parameter, which allows to fix the parameter (not to fit)
// -sInitX - initial value of sx parameter
// -fixSx - logical value of sx parameter, which allows to fix the parameter (not to fit)
// -sInitY - initial value of sy parameter
// -fixSy - logical value of sy parameter, which allows to fix the parameter (not to fit)
//////////////////////////////////////////////////////////////////////////////
fTxyInit = tInitXY;
fFixTxy = fixTxy;
fTxInit = tInitX;
fFixTx = fixTx;
fTyInit = tInitY;
fFixTy = fixTy;
fBxInit = bInitX;
fFixBx = fixBx;
fByInit = bInitY;
fFixBy = fixBy;
fSxyInit = sInitXY;
fFixSxy = fixSxy;
fSxInit = sInitX;
fFixSx = fixSx;
fSyInit = sInitY;
fFixSy = fixSy;
}
//______________________________________________________________________________
void TSpectrum2Fit::GetPositions(Float_t *positionsX, Float_t *positionsY, Float_t *positionsX1, Float_t *positionsY1)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the positions of fitted 2D peaks and 1D ridges
// -positionX - gets vector of x positions of 2D peaks
// -positionY - gets vector of y positions of 2D peaks
// -positionX1 - gets vector of x positions of 1D ridges
// -positionY1 - gets vector of y positions of 1D ridges
//////////////////////////////////////////////////////////////////////////////
for( Int_t i=0; i < fNPeaks; i++){
positionsX[i] = (Float_t) fPositionCalcX[i];
positionsY[i] = (Float_t) fPositionCalcY[i];
positionsX1[i] = (Float_t) fPositionCalcX1[i];
positionsY1[i] = (Float_t) fPositionCalcY1[i];
}
}
//______________________________________________________________________________
void TSpectrum2Fit::GetPositionErrors(Float_t *positionErrorsX, Float_t *positionErrorsY, Float_t *positionErrorsX1, Float_t *positionErrorsY1)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the errors of positions of fitted 2D peaks and 1D ridges
// -positionErrorsX - gets vector of errors of x positions of 2D peaks
// -positionErrorsY - gets vector of errors of y positions of 2D peaks
// -positionErrorsX1 - gets vector of errors of x positions of 1D ridges
// -positionErrorsY1 - gets vector of errors of y positions of 1D ridges
//////////////////////////////////////////////////////////////////////////////
for( Int_t i=0; i < fNPeaks; i++){
positionErrorsX[i] = (Float_t) fPositionErrX[i];
positionErrorsY[i] = (Float_t) fPositionErrY[i];
positionErrorsX1[i] = (Float_t) fPositionErrX1[i];
positionErrorsY1[i] = (Float_t) fPositionErrY1[i];
}
}
//______________________________________________________________________________
void TSpectrum2Fit::GetAmplitudes(Float_t *amplitudes, Float_t *amplitudesX1, Float_t *amplitudesY1)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the amplitudes of fitted 2D peaks and 1D ridges
// -amplitudes - gets vector of amplitudes of 2D peaks
// -amplitudesX1 - gets vector of amplitudes of 1D ridges in x direction
// -amplitudesY1 - gets vector of amplitudes of 1D ridges in y direction
//////////////////////////////////////////////////////////////////////////////
for( Int_t i=0; i < fNPeaks; i++){
amplitudes[i] = (Float_t) fAmpCalc[i];
amplitudesX1[i] = (Float_t) fAmpCalcX1[i];
amplitudesY1[i] = (Float_t) fAmpCalcY1[i];
}
}
//______________________________________________________________________________
void TSpectrum2Fit::GetAmplitudeErrors(Float_t *amplitudeErrors, Float_t *amplitudeErrorsX1, Float_t *amplitudeErrorsY1)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the amplitudes of fitted 2D peaks and 1D ridges
// -amplitudeErrors - gets vector of amplitudes errors of 2D peaks
// -amplitudeErrorsX1 - gets vector of amplitudes errors of 1D ridges in x direction
// -amplitudesErrorY1 - gets vector of amplitudes errors of 1D ridges in y direction
//////////////////////////////////////////////////////////////////////////////
for( Int_t i=0; i < fNPeaks; i++){
amplitudeErrors[i] = (Float_t) fAmpErr[i];
amplitudeErrorsX1[i] = (Float_t) fAmpErrX1[i];
amplitudeErrorsY1[i] = (Float_t) fAmpErrY1[i];
}
}
//______________________________________________________________________________
void TSpectrum2Fit::GetVolumes(Float_t *volumes)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the volumes of fitted 2D peaks
// -volumes - gets vector of volumes of 2D peaks
//////////////////////////////////////////////////////////////////////////////
for( Int_t i=0; i < fNPeaks; i++){
volumes[i] = (Float_t) fVolume[i];
}
}
//______________________________________________________________________________
void TSpectrum2Fit::GetVolumeErrors(Float_t *volumeErrors)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets errors of the volumes of fitted 2D peaks
// -volumeErrors - gets vector of volumes errors of 2D peaks
//////////////////////////////////////////////////////////////////////////////
for( Int_t i=0; i < fNPeaks; i++){
volumeErrors[i] = (Float_t) fVolumeErr[i];
}
}
//______________________________________________________________________________
void TSpectrum2Fit::GetSigmaX(Double_t &sigmaX, Double_t &sigmaErrX)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the sigma x parameter and its error
// -sigmaX - gets the fitted value of sigma x parameter
// -sigmaErrX - gets error value of sigma x parameter
//////////////////////////////////////////////////////////////////////////////
sigmaX=fSigmaCalcX;
sigmaErrX=fSigmaErrX;
}
//______________________________________________________________________________
void TSpectrum2Fit::GetSigmaY(Double_t &sigmaY, Double_t &sigmaErrY)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the sigma y parameter and its error
// -sigmaY - gets the fitted value of sigma y parameter
// -sigmaErrY - gets error value of sigma y parameter
//////////////////////////////////////////////////////////////////////////////
sigmaY=fSigmaCalcY;
sigmaErrY=fSigmaErrY;
}
//______________________________________________________________________________
void TSpectrum2Fit::GetRo(Double_t &ro, Double_t &roErr)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the ro parameter and its error
// -ro - gets the fitted value of ro parameter
// -roErr - gets error value of ro parameter
//////////////////////////////////////////////////////////////////////////////
ro=fRoCalc;
roErr=fRoErr;
}
//______________________________________________________________________________
void TSpectrum2Fit::GetBackgroundParameters(Double_t &a0, Double_t &a0Err, Double_t &ax, Double_t &axErr, Double_t &ay, Double_t &ayErr)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the background parameters and their errors
// -a0 - gets the fitted value of a0 parameter
// -a0Err - gets error value of a0 parameter
// -ax - gets the fitted value of ax parameter
// -axErr - gets error value of ax parameter
// -ay - gets the fitted value of ay parameter
// -ayErr - gets error value of ay parameter
//////////////////////////////////////////////////////////////////////////////
a0 = fA0Calc;
a0Err = fA0Err;
ax = fAxCalc;
axErr = fAxErr;
ay = fAyCalc;
ayErr = fAyErr;
}
//______________________________________________________________________________
void TSpectrum2Fit::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)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the tail parameters and their errors
// -txy - gets the fitted value of txy parameter
// -txyErr - gets error value of txy parameter
// -tx - gets the fitted value of tx parameter
// -txErr - gets error value of tx parameter
// -ty - gets the fitted value of ty parameter
// -tyErr - gets error value of ty parameter
// -bx - gets the fitted value of bx parameter
// -bxErr - gets error value of bx parameter
// -by - gets the fitted value of by parameter
// -byErr - gets error value of by parameter
// -sxy - gets the fitted value of sxy parameter
// -sxyErr - gets error value of sxy parameter
// -sx - gets the fitted value of sx parameter
// -sxErr - gets error value of sx parameter
// -sy - gets the fitted value of sy parameter
// -syErr - gets error value of sy parameter
//////////////////////////////////////////////////////////////////////////////
txy = fTxyCalc;
txyErr = fTxyErr;
tx = fTxCalc;
txErr = fTxErr;
ty = fTyCalc;
tyErr = fTyErr;
bx = fBxCalc;
bxErr = fBxErr;
by = fByCalc;
byErr = fByErr;
sxy = fSxyCalc;
sxyErr = fSxyErr;
sx = fSxCalc;
sxErr = fSxErr;
sy = fSyCalc;
syErr = fSyErr;
}