Fitting
Goal: to estimate
simultaneously peak shape parameters in spectra with large number of peaks
peaks can be fitted separately, each peak (or multiplets) in a region or
together all peaks in a spectrum. To fit separately each peak one needs to
determine the fitted region. However it can happen that the regions of
neighboring peaks are overlapping. Then the results of fitting are very poor.
On the other hand, when fitting together all peaks found in a spectrum, one
needs to have a method that is stable (converges) and fast enough to carry out
fitting in reasonable time
we have implemented the nonsymmetrical semiempirical peak shape function
[1]
it contains the symmetrical Gaussian as well as nonsymmetrical terms.
where T and S are relative amplitudes
and B is slope.
algorithm without matrix inversion (AWMI) allows fitting tens, hundreds
of peaks simultaneously that represent sometimes thousands of parameters [2],
[5].
Function:
void TSpectrumFit::FitAwmi(float *fSource)
This function fits the source
spectrum using AWMI algorithm. The calling program should fill in input fitting
parameters of the TSpectrumFit class using a set of TSpectrumFit setters. The
fitted parameters are written into the class and the fitted data are written
into source spectrum.
Parameter:
fSource-pointer to
the vector of source spectrum
Member variables of the
TSpectrumFit class:
Int_t fNPeaks; //number of peaks present in fit, input
parameter, it should be > 0
Int_t fNumberIterations; //number of iterations in fitting
procedure, input parameter, it should be > 0
Int_t fXmin; //first fitted channel
Int_t fXmax; //last fitted channel
Int_t fStatisticType; //type of statistics, possible values
kFitOptimChiCounts (chi square statistics with counts as weighting
coefficients), kFitOptimChiFuncValues (chi square statistics with function
values as weighting coefficients),kFitOptimMaxLikelihood
Int_t fAlphaOptim; //optimization of convergence algorithm, possible
values kFitAlphaHalving, kFitAlphaOptimal
Int_t fPower; //possible values kFitPower2,4,6,8,10,12,
for details see references. It applies only for Awmi fitting function.
Int_t fFitTaylor; //order of Taylor expansion, possible
values kFitTaylorOrderFirst, kFitTaylorOrderSecond. It applies only for Awmi
fitting function.
Double_t fAlpha; //convergence coefficient, input
parameter, it should be positive number and <=1, for details see references
Double_t fChi; //here the fitting functions return
resulting chi square
Double_t *fPositionInit; //[fNPeaks] array of initial values of
peaks positions, input parameters
Double_t *fPositionCalc; //[fNPeaks] array of calculated values of
fitted positions, output parameters
Double_t *fPositionErr; //[fNPeaks] array of position errors
Double_t *fAmpInit; //[fNPeaks] array of initial values of
peaks amplitudes, input parameters
Double_t *fAmpCalc; //[fNPeaks] array of calculated values of
fitted amplitudes, output parameters
Double_t *fAmpErr; //[fNPeaks] array of amplitude errors
Double_t *fArea; //[fNPeaks] array of calculated areas of
peaks
Double_t *fAreaErr; //[fNPeaks] array of errors of peak areas
Double_t fSigmaInit; //initial value of sigma parameter
Double_t fSigmaCalc; //calculated value of sigma parameter
Double_t fSigmaErr; //error value of sigma parameter
Double_t fTInit; //initial value of t parameter (relative
amplitude of tail), for details see html manual and references
Double_t fTCalc; //calculated value of t parameter
Double_t fTErr; //error value of t parameter
Double_t fBInit; //initial value of b parameter (slope),
for details see html manual and references
Double_t fBCalc; //calculated value of b parameter
Double_t fBErr; //error value of b parameter
Double_t fSInit; //initial value of s parameter (relative
amplitude of step), for details see html manual and references
Double_t fSCalc; //calculated value of s parameter
Double_t fSErr; //error value of s parameter
Double_t fA0Init; //initial value of background a0
parameter(backgroud is estimated as a0+a1*x+a2*x*x)
Double_t fA0Calc; //calculated value of background a0
parameter
Double_t fA0Err; //error value of background a0 parameter
Double_t fA1Init; //initial value of background a1
parameter(backgroud is estimated as a0+a1*x+a2*x*x)
Double_t fA1Calc; //calculated value of background a1
parameter
Double_t fA1Err; //error value of background a1 parameter
Double_t fA2Init; //initial value of background a2
parameter(backgroud is estimated as a0+a1*x+a2*x*x)
Double_t fA2Calc; //calculated value of background a2
parameter
Double_t fA2Err; //error value of background a2 parameter
Bool_t *fFixPosition; //[fNPeaks] array of logical values which
allow to fix appropriate positions (not fit). However they are present in the
estimated functional
Bool_t *fFixAmp; //[fNPeaks] array of logical values which
allow to fix appropriate amplitudes (not fit). However they are present in the
estimated functional
Bool_t fFixSigma; //logical value of sigma parameter, which
allows to fix the parameter (not to fit).
Bool_t fFixT; //logical value of t parameter, which
allows to fix the parameter (not to fit).
Bool_t fFixB; //logical value of b parameter, which
allows to fix the parameter (not to fit).
Bool_t fFixS; //logical value of s parameter, which
allows to fix the parameter (not to fit).
Bool_t fFixA0; //logical value of a0 parameter, which
allows to fix the parameter (not to fit).
Bool_t fFixA1; //logical value of a1 parameter, which
allows to fix the parameter (not to fit).
Bool_t fFixA2; //logical value of a2 parameter, which
allows to fix the parameter (not to fit).
References:
[1] Phillps G.W., Marlow K.W.,
NIM 137 (1976) 525.
[2] I. A. Slavic: Nonlinear
least-squares fitting without matrix inversion applied to complex Gaussian
spectra analysis. NIM 134 (1976) 285-289.
[3] T. Awaya: A new method for
curve fitting to the data with low statistics not using chi-square method. NIM
165 (1979) 317-323.
[4] T. Hauschild, M. Jentschel:
Comparison of maximum likelihood estimation and chi-square statistics applied
to counting experiments. NIM A 457 (2001) 384-401.
[5] M. Morháč, J.
Kliman, M. Jandel, Ľ. Krupa, V. Matouek: Study of fitting algorithms
applied to simultaneous analysis of large number of peaks in -ray spectra. Applied Spectroscopy, Vol. 57, No. 7, pp. 753-760, 2003
Example script FitAwmi.c:
Fig. 1 Original spectrum
(black line) and fitted spectrum using AWMI algorithm (red line) and number of
iteration steps = 1000. Positions of fitted peaks are denoted by markers
Script:
// Example to illustrate
fitting function using AWMI algorithm.
// To execute this example,
do
// root > .x FitAwmi.C
void FitAwmi() {
Double_t a;
Int_t
i,nfound=0,bin;
Int_t nbins = 256;
Int_t xmin = 0;
Int_t xmax =
nbins;
Float_t * source =
new float[nbins];
Float_t * dest =
new float[nbins];
TH1F *h = new
TH1F("h","Fitting using AWMI algorithm",nbins,xmin,xmax);
TH1F *d = new
TH1F("d","",nbins,xmin,xmax);
TFile *f = new
TFile("TSpectrum.root");
h=(TH1F*)
f->Get("fit;1");
for (i = 0; i <
nbins; i++) source[i]=h->GetBinContent(i + 1);
TCanvas *Fit1 =
gROOT->GetListOfCanvases()->FindObject("Fit1");
if (!Fit1) Fit1 =
new TCanvas("Fit1","Fit1",10,10,1000,700);
h->Draw("L");
TSpectrum *s = new
TSpectrum();
//searching for
candidate peaks positions
nfound =
s->SearchHighRes(source, dest, nbins, 2, 0.1, kFALSE, 10000, kFALSE, 0);
Bool_t *FixPos =
new Bool_t[nfound];
Bool_t *FixAmp =
new Bool_t[nfound];
for(i = 0; i<
nfound ; i++){
FixPos[i] =
kFALSE;
FixAmp[i] =
kFALSE;
}
//filling in the
initial estimates of the input parameters
Float_t *PosX =
new Float_t[nfound];
Float_t *PosY =
new Float_t[nfound];
PosX =
s->GetPositionX();
for (i = 0; i <
nfound; i++) {
a=PosX[i];
bin = 1 +
Int_t(a + 0.5);
PosY[i] =
h->GetBinContent(bin);
}
TSpectrumFit
*pfit=new TSpectrumFit(nfound);
pfit->SetFitParameters(xmin, xmax-1, 1000, 0.1, pfit->kFitOptimChiCounts,
pfit->kFitAlphaHalving, pfit->kFitPower2,
pfit->kFitTaylorOrderFirst);
pfit->SetPeakParameters(2,
kFALSE, PosX, (Bool_t *) FixPos, PosY, (Bool_t *) FixAmp);
pfit->FitAwmi(source);
Double_t
*CalcPositions = new Double_t[nfound];
Double_t
*CalcAmplitudes = new Double_t[nfound];
CalcPositions=pfit->GetPositions();
CalcAmplitudes=pfit->GetAmplitudes();
for (i = 0; i <
nbins; i++) d->SetBinContent(i + 1,source[i]);
d->SetLineColor(kRed);
d->Draw("SAME L");
for (i = 0; i <
nfound; i++) {
a=CalcPositions[i];
bin = 1 +
Int_t(a + 0.5);
PosX[i] =
d->GetBinCenter(bin);
PosY[i] =
d->GetBinContent(bin);
}
TPolyMarker * pm =
(TPolyMarker*)h->GetListOfFunctions()->FindObject("TPolyMarker");
if (pm) {
h->GetListOfFunctions()->Remove(pm);
delete pm;
}
pm = new
TPolyMarker(nfound, PosX, PosY);
h->GetListOfFunctions()->Add(pm);
pm->SetMarkerStyle(23);
pm->SetMarkerColor(kRed);
pm->SetMarkerSize(1);
}
End_Html
Int_t i, j, k, shift =
2 * fNPeaks + 7, peak_vel, rozmer, 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 * (2 * fNPeaks + 7)];
for (i = 0, j = 0; i < fNPeaks; i++) {
working_space[2 * i] = fAmpInit[i]; //vector parameter
if (fFixAmp[i] == false) {
working_space[shift + j] = fAmpInit[i]; //vector xk
j += 1;
}
working_space[2 * i + 1] = fPositionInit[i]; //vector parameter
if (fFixPosition[i] == false) {
working_space[shift + j] = fPositionInit[i]; //vector xk
j += 1;
}
}
peak_vel = 2 * i;
working_space[2 * i] = fSigmaInit; //vector parameter
if (fFixSigma == false) {
working_space[shift + j] = fSigmaInit; //vector xk
j += 1;
}
working_space[2 * i + 1] = fTInit; //vector parameter
if (fFixT == false) {
working_space[shift + j] = fTInit; //vector xk
j += 1;
}
working_space[2 * i + 2] = fBInit; //vector parameter
if (fFixB == false) {
working_space[shift + j] = fBInit; //vector xk
j += 1;
}
working_space[2 * i + 3] = fSInit; //vector parameter
if (fFixS == false) {
working_space[shift + j] = fSInit; //vector xk
j += 1;
}
working_space[2 * i + 4] = fA0Init; //vector parameter
if (fFixA0 == false) {
working_space[shift + j] = fA0Init; //vector xk
j += 1;
}
working_space[2 * i + 5] = fA1Init; //vector parameter
if (fFixA1 == false) {
working_space[shift + j] = fA1Init; //vector xk
j += 1;
}
working_space[2 * i + 6] = fA2Init; //vector parameter
if (fFixA2 == false) {
working_space[shift + j] = fA2Init; //vector xk
j += 1;
}
rozmer = j;
if (rozmer == 0){
delete [] working_space;
Error ("FitAwmi","All parameters are fixed");
return;
}
if (rozmer >= fXmax - fXmin + 1){
delete [] working_space;
Error ("FitAwmi","Number of fitted parameters is larger than # of fitted points");
return;
}
for (iter = 0; iter < fNumberIterations; iter++) {
for (j = 0; j < rozmer; 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 (i = fXmin; i <= fXmax; i++) {
yw = source[i];
ywm = yw;
f = Shape(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel], working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6]);
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.001)
ywm = 0.001;
}
else {
if (ywm == 0)
ywm = 1;
}
//calculation of gradient vector
for (j = 0, k = 0; j < fNPeaks; j++) {
if (fFixAmp[j] == false) {
a = Deramp((Double_t) i, working_space[2 * j + 1],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2]);
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 (fFixPosition[j] == false) {
a = Deri0((Double_t) i, working_space[2 * j],
working_space[2 * j + 1],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2]);
if (fFitTaylor == kFitTaylorOrderSecond)
d = Derderi0((Double_t) i, working_space[2 * j],
working_space[2 * j + 1],
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 (fFixSigma == false) {
a = Dersigma(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2]);
if (fFitTaylor == kFitTaylorOrderSecond)
d = Derdersigma(fNPeaks, (Double_t) i,
working_space, 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 (fFixT == false) {
a = Dert(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel],
working_space[peak_vel + 2]);
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 (fFixB == false) {
a = Derb(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel], working_space[peak_vel + 1],
working_space[peak_vel + 2]);
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 (fFixS == false) {
a = Ders(fNPeaks, (Double_t) i, 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 (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 (fFixA1 == false) {
a = Dera1((Double_t) i);
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 (fFixA2 == false) {
a = Dera2((Double_t) i);
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 < rozmer; 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 < rozmer; 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 < rozmer; 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[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
j += 1;
}
if (fFixPosition[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[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
j += 1;
}
}
if (fFixSigma == 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 (fFixT == false) {
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixB == 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 + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixS == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixA1 == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixA2 == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
chi2 = 0;
for (i = fXmin; i <= fXmax; i++) {
yw = source[i];
ywm = yw;
f = Shape(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6]);
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 < rozmer; 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[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
j += 1;
}
if (fFixPosition[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[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
j += 1;
}
}
if (fFixSigma == 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 (fFixT == false) {
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixB == 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 + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixS == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixA1 == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixA2 == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
chi = chi_min;
}
}
else {
for (j = 0; j < rozmer; 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[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
j += 1;
}
if (fFixPosition[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[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
j += 1;
}
}
if (fFixSigma == 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 (fFixT == false) {
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixB == 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 + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixS == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixA1 == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixA2 == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
chi = 0;
for (i = fXmin; i <= fXmax; i++) {
yw = source[i];
ywm = yw;
f = Shape(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6]);
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 < rozmer; j++) {
working_space[4 * shift + j] = 0; //temp_xk[j]
working_space[2 * shift + j] = 0; //der[j]
}
for (i = fXmin, chi_cel = 0; i <= fXmax; i++) {
yw = source[i];
if (yw == 0)
yw = 1;
f = Shape(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel], working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6]);
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 = Deramp((Double_t) i, working_space[2 * j + 1],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
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 (fFixPosition[j] == false) {
a = Deri0((Double_t) i, working_space[2 * j],
working_space[2 * j + 1],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
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 (fFixSigma == false) {
a = Dersigma(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
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 (fFixT == false) {
a = Dert(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel],
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 (fFixB == false) {
a = Derb(fNPeaks, (Double_t) i, 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 (fFixS == false) {
a = Ders(fNPeaks, (Double_t) i, 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; //tem_xk[k]
}
k += 1;
}
if (fFixA0 == false) {
a = 1.0;
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 (fFixA1 == false) {
a = Dera1((Double_t) i);
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 (fFixA2 == false) {
a = Dera2((Double_t) i);
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 - rozmer;
chi_er = chi_cel / b;
for (i = 0, j = 0; i < fNPeaks; i++) {
fArea[i] =
Area(working_space[2 * i], working_space[peak_vel],
working_space[peak_vel + 1], working_space[peak_vel + 2]);
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]
if (fArea[i] > 0) {
a = Derpa(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;
fAreaErr[i] = TMath::Sqrt(TMath::Abs(a * a * b * chi_er));
}
else
fAreaErr[i] = 0;
j += 1;
}
else {
fAmpCalc[i] = fAmpInit[i];
fAmpErr[i] = 0;
fAreaErr[i] = 0;
}
if (fFixPosition[i] == false) {
fPositionCalc[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fPositionErr[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 {
fPositionCalc[i] = fPositionInit[i];
fPositionErr[i] = 0;
}
}
if (fFixSigma == false) {
fSigmaCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSigmaErr = 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 {
fSigmaCalc = fSigmaInit;
fSigmaErr = 0;
}
if (fFixT == false) {
fTCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fTErr = 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 {
fTCalc = fTInit;
fTErr = 0;
}
if (fFixB == false) {
fBCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fBErr = 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 {
fBCalc = fBInit;
fBErr = 0;
}
if (fFixS == false) {
fSCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSErr = 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 {
fSCalc = fSInit;
fSErr = 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 (fFixA1 == false) {
fA1Calc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fA1Err = 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 {
fA1Calc = fA1Init;
fA1Err = 0;
}
if (fFixA2 == false) {
fA2Calc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fA2Err = 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 {
fA2Calc = fA2Init;
fA2Err = 0;
}
b = fXmax - fXmin + 1 - rozmer;
fChi = chi_cel / b;
for (i = fXmin; i <= fXmax; i++) {
f = Shape(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel], working_space[peak_vel + 1],
working_space[peak_vel + 3], working_space[peak_vel + 2],
working_space[peak_vel + 4], working_space[peak_vel + 5],
working_space[peak_vel + 6]);
source[i] = f;
}
delete[]working_space;
return;
}
/////////////////FITTING FUNCTION WITH MATRIX INVERSION///////////////////////////////////////
//_______________________________________________________________________________
void TSpectrumFit::StiefelInversion(Double_t **a, Int_t size)
{
//////////////////////////////////////////////////////////////////////////////////
// AUXILIARY FUNCTION //
// //
// This function calculates solution of the system of linear equations. //
// The matrix a should have a dimension size*(size+4) //
// The calling function should fill in the matrix, the column size should //
// contain vector y (right side of the system of equations). The result is //
// placed into size+1 column of the matrix. //
// according to sigma of peaks. //
// Function parameters: //
// -a-matrix with dimension size*(size+4) // //
// -size-number of rows of the matrix //
// //
//////////////////////////////////////////////////////////////////////////////////
Int_t i, j, k = 0;
Double_t sk = 0, b, lambdak, normk, normk_old = 0;
do {
normk = 0;
//calculation of rk and norm
for (i = 0; i < size; i++) {
a[i][size + 2] = -a[i][size]; //rk=-C
for (j = 0; j < size; j++) {
a[i][size + 2] += a[i][j] * a[j][size + 1]; //A*xk-C
}
normk += a[i][size + 2] * a[i][size + 2]; //calculation normk
}
//calculation of sk
if (k != 0) {
sk = normk / normk_old;
}
//calculation of uk
for (i = 0; i < size; i++) {
a[i][size + 3] = -a[i][size + 2] + sk * a[i][size + 3]; //uk=-rk+sk*uk-1
}
//calculation of lambdak
lambdak = 0;
for (i = 0; i < size; i++) {
for (j = 0, b = 0; j < size; j++) {
b += a[i][j] * a[j][size + 3]; //A*uk
}
lambdak += b * a[i][size + 3];
}
if (TMath::Abs(lambdak) > 1e-50) //computer zero
lambdak = normk / lambdak;
else
lambdak = 0;
for (i = 0; i < size; i++)
a[i][size + 1] += lambdak * a[i][size + 3]; //xk+1=xk+lambdak*uk
normk_old = normk;
k += 1;
} while (k < size && TMath::Abs(normk) > 1e-50); //computer zero
return;
}
//_______________________________________________________________________________
void TSpectrumFit::FitStiefel(Float_t *source)
{
/////////////////////////////////////////////////////////////////////////////
// ONE-DIMENSIONAL FIT FUNCTION
// ALGORITHM WITH MATRIX INVERSION (STIEFEL-HESTENS METHOD)
// This function fits the source spectrum. The calling program should
// fill in input parameters
// The fitted parameters are written into
// output parameters and fitted data are written into
// source spectrum.
//
// Function parameters:
// source-pointer to the vector of source spectrum
//
/////////////////////////////////////////////////////////////////////////////
//Begin_Html
Stiefel fitting algorithm
Function:
void TSpectrumFit::FitStiefel(float *fSource)
This function fits the source
spectrum using Stiefel-Hestens method [1] (see Awmi function). The calling
program should fill in input fitting parameters of the TSpectrumFit class using
a set of TSpectrumFit setters. The fitted parameters are written into the class
and the fitted data are written into source spectrum. It converges faster than
Awmi method.
Parameter:
fSource-pointer to
the vector of source spectrum
Example script FitStiefel.c:
Fig. 2 Original spectrum
(black line) and fitted spectrum using Stiefel-Hestens method (red line) and
number of iteration steps = 100. Positions of fitted peaks are denoted by
markers
Script:
// Example to illustrate
fitting function using Stiefel-Hestens method.
// To execute this example, do
// root > .x FitStiefel.C
void FitStiefel() {
Double_t a;
Int_t i,nfound=0,bin;
Int_t nbins = 256;
Int_t xmin = 0;
Int_t xmax =
nbins;
Float_t * source = new float[nbins];
Float_t * dest = new
float[nbins];
TH1F *h = new TH1F("h","Fitting
using AWMI algorithm",nbins,xmin,xmax);
TH1F *d = new
TH1F("d","",nbins,xmin,xmax);
TFile *f = new
TFile("TSpectrum.root");
h=(TH1F*)
f->Get("fit;1");
for (i = 0; i < nbins;
i++) source[i]=h->GetBinContent(i + 1);
TCanvas *Fit1 =
gROOT->GetListOfCanvases()->FindObject("Fit1");
if (!Fit1) Fit1 = new
TCanvas("Fit1","Fit1",10,10,1000,700);
h->Draw("L");
TSpectrum *s = new
TSpectrum();
//searching for candidate
peaks positions
nfound =
s->SearchHighRes(source, dest, nbins, 2, 0.1, kFALSE, 10000, kFALSE, 0);
Bool_t *FixPos = new
Bool_t[nfound];
Bool_t *FixAmp = new
Bool_t[nfound];
for(i = 0; i< nfound ;
i++){
FixPos[i] = kFALSE;
FixAmp[i] = kFALSE;
}
//filling in the initial
estimates of the input parameters
Float_t *PosX = new
Float_t[nfound];
Float_t *PosY = new
Float_t[nfound];
PosX =
s->GetPositionX();
for (i = 0; i < nfound;
i++) {
a=PosX[i];
bin = 1 + Int_t(a +
0.5);
PosY[i] =
h->GetBinContent(bin);
}
TSpectrumFit *pfit=new
TSpectrumFit(nfound);
pfit->SetFitParameters(xmin, xmax-1, 1000, 0.1, pfit->kFitOptimChiCounts,
pfit->kFitAlphaHalving, pfit->kFitPower2,
pfit->kFitTaylorOrderFirst);
pfit->SetPeakParameters(2,
kFALSE, PosX, (Bool_t *) FixPos, PosY, (Bool_t *) FixAmp);
pfit->FitStiefel(source);
Double_t *CalcPositions =
new Double_t[nfound];
Double_t *CalcAmplitudes = new
Double_t[nfound];
CalcPositions=pfit->GetPositions();
CalcAmplitudes=pfit->GetAmplitudes();
for (i = 0; i < nbins;
i++) d->SetBinContent(i + 1,source[i]);
d->SetLineColor(kRed);
d->Draw("SAME
L");
for (i = 0; i < nfound;
i++) {
a=CalcPositions[i];
bin = 1 + Int_t(a +
0.5);
PosX[i] =
d->GetBinCenter(bin);
PosY[i] =
d->GetBinContent(bin);
}
TPolyMarker * pm =
(TPolyMarker*)h->GetListOfFunctions()->FindObject("TPolyMarker");
if (pm) {
h->GetListOfFunctions()->Remove(pm);
delete pm;
}
pm = new
TPolyMarker(nfound, PosX, PosY);
h->GetListOfFunctions()->Add(pm);
pm->SetMarkerStyle(23);
pm->SetMarkerColor(kRed);
pm->SetMarkerSize(1);
}
End_Html
Int_t i, j, k, shift =
2 * fNPeaks + 7, peak_vel, rozmer, iter, regul_cycle,
flag;
Double_t a, b, 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 * (2 * fNPeaks + 7)];
for (i = 0, j = 0; i < fNPeaks; i++) {
working_space[2 * i] = fAmpInit[i]; //vector parameter
if (fFixAmp[i] == false) {
working_space[shift + j] = fAmpInit[i]; //vector xk
j += 1;
}
working_space[2 * i + 1] = fPositionInit[i]; //vector parameter
if (fFixPosition[i] == false) {
working_space[shift + j] = fPositionInit[i]; //vector xk
j += 1;
}
}
peak_vel = 2 * i;
working_space[2 * i] = fSigmaInit; //vector parameter
if (fFixSigma == false) {
working_space[shift + j] = fSigmaInit; //vector xk
j += 1;
}
working_space[2 * i + 1] = fTInit; //vector parameter
if (fFixT == false) {
working_space[shift + j] = fTInit; //vector xk
j += 1;
}
working_space[2 * i + 2] = fBInit; //vector parameter
if (fFixB == false) {
working_space[shift + j] = fBInit; //vector xk
j += 1;
}
working_space[2 * i + 3] = fSInit; //vector parameter
if (fFixS == false) {
working_space[shift + j] = fSInit; //vector xk
j += 1;
}
working_space[2 * i + 4] = fA0Init; //vector parameter
if (fFixA0 == false) {
working_space[shift + j] = fA0Init; //vector xk
j += 1;
}
working_space[2 * i + 5] = fA1Init; //vector parameter
if (fFixA1 == false) {
working_space[shift + j] = fA1Init; //vector xk
j += 1;
}
working_space[2 * i + 6] = fA2Init; //vector parameter
if (fFixA2 == false) {
working_space[shift + j] = fA2Init; //vector xk
j += 1;
}
rozmer = j;
if (rozmer == 0){
Error ("FitAwmi","All parameters are fixed");
delete [] working_space;
return;
}
if (rozmer >= fXmax - fXmin + 1){
Error ("FitAwmi","Number of fitted parameters is larger than # of fitted points");
delete [] working_space;
return;
}
Double_t **working_matrix = new Double_t *[rozmer];
for (i = 0; i < rozmer; i++)
working_matrix[i] = new Double_t[rozmer + 4];
for (iter = 0; iter < fNumberIterations; iter++) {
for (j = 0; j < rozmer; j++) {
working_space[3 * shift + j] = 0; //temp
for (k = 0; k <= rozmer; k++) {
working_matrix[j][k] = 0;
}
}
//filling working matrix
alpha = fAlpha;
chi_opt = 0;
for (i = fXmin; i <= fXmax; i++) {
//calculation of gradient vector
for (j = 0, k = 0; j < fNPeaks; j++) {
if (fFixAmp[j] == false) {
working_space[2 * shift + k] =
Deramp((Double_t) i, working_space[2 * j + 1],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2]);
k += 1;
}
if (fFixPosition[j] == false) {
working_space[2 * shift + k] =
Deri0((Double_t) i, working_space[2 * j],
working_space[2 * j + 1], working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2]);
k += 1;
}
} if (fFixSigma == false) {
working_space[2 * shift + k] =
Dersigma(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2]);
k += 1;
}
if (fFixT == false) {
working_space[2 * shift + k] =
Dert(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel], working_space[peak_vel + 2]);
k += 1;
}
if (fFixB == false) {
working_space[2 * shift + k] =
Derb(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel], working_space[peak_vel + 1],
working_space[peak_vel + 2]);
k += 1;
}
if (fFixS == false) {
working_space[2 * shift + k] =
Ders(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel]);
k += 1;
}
if (fFixA0 == false) {
working_space[2 * shift + k] = 1.;
k += 1;
}
if (fFixA1 == false) {
working_space[2 * shift + k] = Dera1((Double_t) i);
k += 1;
}
if (fFixA2 == false) {
working_space[2 * shift + k] = Dera2((Double_t) i);
k += 1;
}
yw = source[i];
ywm = yw;
f = Shape(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel], working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6]);
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 < rozmer; j++) {
for (k = 0; k < rozmer; 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 < rozmer; j++) {
working_matrix[j][rozmer] -= b * working_space[2 * shift + j];
}
}
for (i = 0; i < rozmer; i++) {
working_matrix[i][rozmer + 1] = 0; //xk
}
StiefelInversion(working_matrix, rozmer);
for (i = 0; i < rozmer; i++) {
working_space[2 * shift + i] = working_matrix[i][rozmer + 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 < rozmer; 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 < rozmer; 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[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
j += 1;
}
if (fFixPosition[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[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
j += 1;
}
}
if (fFixSigma == 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 (fFixT == false) {
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixB == 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 + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixS == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixA1 == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixA2 == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
chi2 = 0;
for (i = fXmin; i <= fXmax; i++) {
yw = source[i];
ywm = yw;
f = Shape(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6]);
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 < rozmer; 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[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
j += 1;
}
if (fFixPosition[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[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
j += 1;
}
}
if (fFixSigma == 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 (fFixT == false) {
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixB == 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 + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixS == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixA1 == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixA2 == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
chi = chi_min;
}
}
else {
for (j = 0; j < rozmer; j++) {
working_space[shift + j] = working_space[4 * shift + j] + alpha * working_space[2 * shift + j]; //xk[j]=temp_xk[j]+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[2 * i] = working_space[shift + j]; //parameter[2*i]=xk[j]
j += 1;
}
if (fFixPosition[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[2 * i + 1] = working_space[shift + j]; //parameter[2*i+1]=xk[j]
j += 1;
}
}
if (fFixSigma == 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 (fFixT == false) {
working_space[peak_vel + 1] = working_space[shift + j]; //parameter[peak_vel+1]=xk[j]
j += 1;
}
if (fFixB == 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 + 2] = working_space[shift + j]; //parameter[peak_vel+2]=xk[j]
j += 1;
}
if (fFixS == false) {
working_space[peak_vel + 3] = working_space[shift + j]; //parameter[peak_vel+3]=xk[j]
j += 1;
}
if (fFixA0 == false) {
working_space[peak_vel + 4] = working_space[shift + j]; //parameter[peak_vel+4]=xk[j]
j += 1;
}
if (fFixA1 == false) {
working_space[peak_vel + 5] = working_space[shift + j]; //parameter[peak_vel+5]=xk[j]
j += 1;
}
if (fFixA2 == false) {
working_space[peak_vel + 6] = working_space[shift + j]; //parameter[peak_vel+6]=xk[j]
j += 1;
}
chi = 0;
for (i = fXmin; i <= fXmax; i++) {
yw = source[i];
ywm = yw;
f = Shape(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6]);
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 < rozmer; j++) {
working_space[4 * shift + j] = 0; //temp_xk[j]
working_space[2 * shift + j] = 0; //der[j]
}
for (i = fXmin, chi_cel = 0; i <= fXmax; i++) {
yw = source[i];
if (yw == 0)
yw = 1;
f = Shape(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel], working_space[peak_vel + 1],
working_space[peak_vel + 3],
working_space[peak_vel + 2],
working_space[peak_vel + 4],
working_space[peak_vel + 5],
working_space[peak_vel + 6]);
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 = Deramp((Double_t) i, working_space[2 * j + 1],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
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 (fFixPosition[j] == false) {
a = Deri0((Double_t) i, working_space[2 * j],
working_space[2 * j + 1],
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
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 (fFixSigma == false) {
a = Dersigma(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel],
working_space[peak_vel + 1],
working_space[peak_vel + 3],
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 (fFixT == false) {
a = Dert(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel],
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 (fFixB == false) {
a = Derb(fNPeaks, (Double_t) i, 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 (fFixS == false) {
a = Ders(fNPeaks, (Double_t) i, 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; //tem_xk[k]
}
k += 1;
}
if (fFixA0 == false) {
a = 1.0;
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 (fFixA1 == false) {
a = Dera1((Double_t) i);
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 (fFixA2 == false) {
a = Dera2((Double_t) i);
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 - rozmer;
chi_er = chi_cel / b;
for (i = 0, j = 0; i < fNPeaks; i++) {
fArea[i] =
Area(working_space[2 * i], working_space[peak_vel],
working_space[peak_vel + 1], working_space[peak_vel + 2]);
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]
if (fArea[i] > 0) {
a = Derpa(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;
fAreaErr[i] = TMath::Sqrt(TMath::Abs(a * a * b * chi_er));
}
else
fAreaErr[i] = 0;
j += 1;
}
else {
fAmpCalc[i] = fAmpInit[i];
fAmpErr[i] = 0;
fAreaErr[i] = 0;
}
if (fFixPosition[i] == false) {
fPositionCalc[i] = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fPositionErr[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 {
fPositionCalc[i] = fPositionInit[i];
fPositionErr[i] = 0;
}
}
if (fFixSigma == false) {
fSigmaCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSigmaErr = 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 {
fSigmaCalc = fSigmaInit;
fSigmaErr = 0;
}
if (fFixT == false) {
fTCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fTErr = 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 {
fTCalc = fTInit;
fTErr = 0;
}
if (fFixB == false) {
fBCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fBErr = 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 {
fBCalc = fBInit;
fBErr = 0;
}
if (fFixS == false) {
fSCalc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fSErr = 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 {
fSCalc = fSInit;
fSErr = 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 (fFixA1 == false) {
fA1Calc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fA1Err = 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 {
fA1Calc = fA1Init;
fA1Err = 0;
}
if (fFixA2 == false) {
fA2Calc = working_space[shift + j]; //xk[j]
if (working_space[3 * shift + j] != 0) //temp[j]
fA2Err = 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 {
fA2Calc = fA2Init;
fA2Err = 0;
}
b = fXmax - fXmin + 1 - rozmer;
fChi = chi_cel / b;
for (i = fXmin; i <= fXmax; i++) {
f = Shape(fNPeaks, (Double_t) i, working_space,
working_space[peak_vel], working_space[peak_vel + 1],
working_space[peak_vel + 3], working_space[peak_vel + 2],
working_space[peak_vel + 4], working_space[peak_vel + 5],
working_space[peak_vel + 6]);
source[i] = f;
}
for (i = 0; i < rozmer; i++)
delete [] working_matrix[i];
delete [] working_matrix;
delete [] working_space;
return;
}
//____________________________________________________________________________
void TSpectrumFit::SetFitParameters(Int_t xmin,Int_t xmax, Int_t numberIterations, Double_t alpha, Int_t statisticType, Int_t alphaOptim, Int_t power, Int_t fitTaylor)
{
//////////////////////////////////////////////////////////////////////////////
// SETTER FUNCTION
//
// This function sets the following fitting parameters:
// -xmin, xmax - 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){
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,fNumberIterations=numberIterations,fAlpha=alpha,fStatisticType=statisticType,fAlphaOptim=alphaOptim,fPower=power,fFitTaylor=fitTaylor;
}
//_______________________________________________________________________________
void TSpectrumFit::SetPeakParameters(Double_t sigma, Bool_t fixSigma, const Float_t *positionInit, const Bool_t *fixPosition, const Float_t *ampInit, const Bool_t *fixAmp)
{
//////////////////////////////////////////////////////////////////////////////
// SETTER FUNCTION
//
// This function sets the following fitting parameters of peaks:
// -sigma - initial value of sigma parameter
// -fixSigma - logical value of sigma parameter, which allows to fix the parameter (not to fit)
// -positionInit - aray of initial values of peaks positions
// -fixPosition - array of logical values which allow to fix appropriate positions (not fit). However they are present in the estimated functional.
// -ampInit - aray of initial values of peaks amplitudes
// -fixAmp - aray of logical values which allow to fix appropriate amplitudes (not fit). However they are present in the estimated functional
//////////////////////////////////////////////////////////////////////////////
Int_t i;
if (sigma <= 0){
Error ("SetPeakParameters","Invalid sigma, must be > than 0");
return;
}
for(i=0; i < fNPeaks; i++){
if((Int_t) positionInit[i] < fXmin || (Int_t) positionInit[i] > fXmax){
Error ("SetPeakParameters","Invalid peak position, must be in the range fXmin, fXmax");
return;
}
if(ampInit[i] < 0){
Error ("SetPeakParameters","Invalid peak amplitude, must be > than 0");
return;
}
}
fSigmaInit = sigma,fFixSigma = fixSigma;
for(i=0; i < fNPeaks; i++){
fPositionInit[i] = (Double_t) positionInit[i];
fFixPosition[i] = fixPosition[i];
fAmpInit[i] = (Double_t) ampInit[i];
fFixAmp[i] = fixAmp[i];
}
}
//_______________________________________________________________________________
void TSpectrumFit::SetBackgroundParameters(Double_t a0Init, Bool_t fixA0, Double_t a1Init, Bool_t fixA1, Double_t a2Init, Bool_t fixA2)
{
//////////////////////////////////////////////////////////////////////////////
// SETTER FUNCTION
//
// This function sets the following fitting parameters of background:
// -a0Init - initial value of a0 parameter (backgroud is estimated as a0+a1*x+a2*x*x)
// -fixA0 - logical value of a0 parameter, which allows to fix the parameter (not to fit)
// -a1Init - initial value of a1 parameter
// -fixA1 - logical value of a1 parameter, which allows to fix the parameter (not to fit)
// -a2Init - initial value of a2 parameter
// -fixA2 - logical value of a2 parameter, which allows to fix the parameter (not to fit)
//////////////////////////////////////////////////////////////////////////////
fA0Init = a0Init;
fFixA0 = fixA0;
fA1Init = a1Init;
fFixA1 = fixA1;
fA2Init = a2Init;
fFixA2 = fixA2;
}
//_______________________________________________________________________________
void TSpectrumFit::SetTailParameters(Double_t tInit, Bool_t fixT, Double_t bInit, Bool_t fixB, Double_t sInit, Bool_t fixS)
{
//////////////////////////////////////////////////////////////////////////////
// SETTER FUNCTION
//
// This function sets the following fitting parameters of tails of peaks
// -tInit - initial value of t parameter
// -fixT - logical value of t parameter, which allows to fix the parameter (not to fit)
// -bInit - initial value of b parameter
// -fixB - logical value of b parameter, which allows to fix the parameter (not to fit)
// -sInit - initial value of s parameter
// -fixS - logical value of s parameter, which allows to fix the parameter (not to fit)
//////////////////////////////////////////////////////////////////////////////
fTInit = tInit;
fFixT = fixT;
fBInit = bInit;
fFixB = fixB;
fSInit = sInit;
fFixS = fixS;
}
//_______________________________________________________________________________
void TSpectrumFit::GetSigma(Double_t &sigma, Double_t &sigmaErr)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the sigma parameter and its error
// -sigma - gets the fitted value of sigma parameter
// -sigmaErr - gets error value of sigma parameter
//////////////////////////////////////////////////////////////////////////////
sigma=fSigmaCalc;
sigmaErr=fSigmaErr;
}
//_______________________________________________________________________________
void TSpectrumFit::GetBackgroundParameters(Double_t &a0, Double_t &a0Err, Double_t &a1, Double_t &a1Err, Double_t &a2, Double_t &a2Err)
{
//////////////////////////////////////////////////////////////////////////////
// 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
// -a1 - gets the fitted value of a1 parameter
// -a1Err - gets error value of a1 parameter
// -a2 - gets the fitted value of a2 parameter
// -a2Err - gets error value of a2 parameter
//////////////////////////////////////////////////////////////////////////////
a0 = fA0Calc;
a0Err = fA0Err;
a1 = fA1Calc;
a1Err = fA1Err;
a2 = fA2Calc;
a2Err = fA2Err;
}
//_______________________________________________________________________________
void TSpectrumFit::GetTailParameters(Double_t &t, Double_t &tErr, Double_t &b, Double_t &bErr, Double_t &s, Double_t &sErr)
{
//////////////////////////////////////////////////////////////////////////////
// GETTER FUNCTION
//
// This function gets the tail parameters and their errors
// -t - gets the fitted value of t parameter
// -tErr - gets error value of t parameter
// -b - gets the fitted value of b parameter
// -bErr - gets error value of b parameter
// -s - gets the fitted value of s parameter
// -sErr - gets error value of s parameter
//////////////////////////////////////////////////////////////////////////////
t = fTCalc;
tErr = fTErr;
b = fBCalc;
bErr = fBErr;
s = fSCalc;
sErr = fSErr;
}
