Hello, everybody!
I am a new Root user. My professor asked my to fit some spectras with comon peak width for all spectras. I saw beautiful example CombinedFit, written by Lorenzo Moneta and rewrote that script with my functions.
First fitting function has 13 parameters (numbers 0-12), Second function has also 13 parameters (numbers 13-25)
For example peak width corresponds to 0th parameter in first function and 13th in second function.
Then I need to change 13th parameter with 0th parameter in second function. Is it correct?
In any case I have huge errors in all parameters.
Here is rewritten script by me
‘’’ #include “Fit/Fitter.h”
#include “Fit/BinData.h”
#include “Fit/Chi2FCN.h”
#include “TH1.h”
#include “TList.h”
#include “Math/WrappedMultiTF1.h”
#include “HFitInterface.h”
#include “TCanvas.h”
#include “TStyle.h”
//объявляем параметры спеткров
int First_spec_pars[13]={
0,//Gcp1-ширина центрального пика первого спектра
1,//x01 -положение центрального пика первого спектра
2,//Acp1-Амплитуда центрального пика первого спектра
3,//G11- ширина первого фононного пика первого спектра
4,//E11 -положение первого фононного пика первого спектра
5,//S11 -амплитуда первого фононного пика первого спектра
6,//G21- ширина второго фононного пика первого спектра
7,//E21 -положение второго фононного пика первого спектра
8,//S21 -амплитуда второго фононного пика первого спектра
9, //G31- ширина третьего фононного пика первого спектра
10,//E31 -положение третьего фононного пика первого спектра
11,//S31 -амплитуда третьего фононного пика первого спектра
12,//T1 -температура первого спектра в Кельвинах
};
int Second_spec_pars[13]={
13,//Gcp2-ширина центрального пика первого спектра
14,//x02 -положение центрального пика первого спектра
15,//Acp2 - Амплитуда центрального пика первого спектра
16,//G12- ширина первого фононного пика первого спектра
17,//E12 -положение первого фононного пика первого спектра
18,//S12 -амплитуда первого фононного пика первого спектра
19,//G22- ширина второго фононного пика первого спектра
20,//E22 -положение второго фононного пика первого спектра
21,//S22 -амплитуда второго фононного пика первого спектра
22,//G32- ширина третьего фононного пика первого спектра
23,//E32 -положение третьего фононного пика первого спектра
24,//S32 -амплитуда третьего фононного пика первого спектра
25,//T2 -температура второго спектра в Кельвинах
};
// Create the GlobalCHi2 structure
struct GlobalChi2 {
GlobalChi2( ROOT::Math::IMultiGenFunction & f1,
ROOT::Math::IMultiGenFunction & f2) :
fChi2_1(&f1), fChi2_2(&f2) {}
double operator() (const double *par) const {
double p1[13];
for (int i = 0; i < 13; ++i) p1[i] = par[First_spec_pars[i] ];
double p2[13];
for (int i = 0; i < 13; ++i) p2[i] = par[Second_spec_pars[i] ];
return (*fChi2_1)(p1) + (*fChi2_2)(p2);
}
const ROOT::Math::IMultiGenFunction * fChi2_1;
const ROOT::Math::IMultiGenFunction * fChi2_2;
};
void combinedFit_Stas_Stas()
{
spectra1=new TGraph("/home/stanislav/Minuit_fit/PZT_2p4_2s_42_440K_a06_minuit");
//spectra1=new TGraph("/home/stanislav/Minuit_fit/PZT_2p4_2s_19_540K_a06");
spectra2=new TGraph("/home/stanislav/Minuit_fit/PZT_2p4_2s_42_440K_a06_minuit");
TF1 * f1 = new TF1(“f”,"[2]1/3.141592 [0] / ( (x-[1])(x-[1]) + [0][0] )+ [5]x/(1-exp( -9.1x/[12] ))[3]/((xx-[4][4])(xx-[4][4])+xx[3])+ [8]x/(1-exp( -9.1x/[12] ))[6]/((xx-[7][7])(xx-[7][7])+xx[6])+ [11]x/(1-exp( -9.1x/[12] ))[9]/((xx-[10][10])(xx-[10][10])+xx[9])",-20,20);
TF1 * f2 = new TF1(“f”,"[15]1/3.141592 [0] / ( (x-[14])(x-[14]) + [0][0] )+ [18]x/(1-exp( -9.1x/[25] ))[16]/((xx-[17][17])(xx-[17][17])+xx[16])+ [21]x/(1-exp( -9.1x/[25] ))[19]/((xx-[20][20])(xx-[20][20])+xx[19])+ [24]x/(1-exp( -9.1x/[25] ))[22]/((xx-[23][23])(xx-[23][23])+xx[22])",-20,20);
//First spectra
// perform now global fit
ROOT::Math::WrappedMultiTF1 wf1(*f1,1);
ROOT::Math::WrappedMultiTF1 wf2(*f2,1);
ROOT::Fit::DataOptions opt;
ROOT::Fit::DataRange rangeB;
// set the data range
rangeB.SetRange(-20,20);
ROOT::Fit::BinData dataB(opt,rangeB);
ROOT::Fit::FillData(dataB, spectra1);
ROOT::Fit::DataRange rangeSB;
rangeSB.SetRange(-20,20);
ROOT::Fit::BinData dataSB(opt,rangeSB);
ROOT::Fit::FillData(dataSB, spectra2);
ROOT::Fit::Chi2Function chi2_B(dataB, wf1);
ROOT::Fit::Chi2Function chi2_SB(dataSB, wf2);
GlobalChi2 globalChi2(chi2_B, chi2_SB);
ROOT::Fit::Fitter fitter;
const int Npar = 27;
// G x0 A G1 E1 S1 G2 E2 S2 G3 E3 S3 T
double par0[Npar] = {4, 0,180 ,35,4,1, 100,14,1, 1,20,1, 440, 4, 0,180 ,35,4,1, 100,14,1, 1,20,1, 440};
// create before the parameter settings in order to fix or set range on them
fitter.Config().SetParamsSettings(26,par0);
// fix 5-th parameter
// first spectra
fitter.Config().ParSettings(0).SetLimits(0,4);//Gcp1 limits
fitter.Config().ParSettings(1).SetLimits(-0.5,0.5);//x01 limits
fitter.Config().ParSettings(2).SetLimits(0,200);//Acp1 limits
fitter.Config().ParSettings(3).SetLimits(0,40); //G11 limits
fitter.Config().ParSettings(4).SetLimits(2,5); //E11 limits
fitter.Config().ParSettings(5).SetLimits(0,4); //S11 limits
fitter.Config().ParSettings(6).SetLimits(0,150); //G21 limits
fitter.Config().ParSettings(7).SetLimits(5,15); //E21 limits
fitter.Config().ParSettings(8).SetLimits(0,4); //S21 limits
fitter.Config().ParSettings(9).SetLimits(0,10); //G31 limits
fitter.Config().ParSettings(10).SetLimits(10,20); //E31 limits
fitter.Config().ParSettings(11).SetLimits(0,10); //S31 limits
fitter.Config().ParSettings(12).Fix();//Fix temperature in first spectra
// second spectra
fitter.Config().ParSettings(13).SetLimits(0,4);//Gcp1 limits
fitter.Config().ParSettings(14).SetLimits(-0.5,0.5);//x01 limits
fitter.Config().ParSettings(15).SetLimits(0,200);//Acp1 limits
fitter.Config().ParSettings(16).SetLimits(0,40); //G12 limits
fitter.Config().ParSettings(17).SetLimits(2,5); //E12 limits
fitter.Config().ParSettings(18).SetLimits(0,4); //S12 limits
fitter.Config().ParSettings(19).SetLimits(0,10); //G22 limits
fitter.Config().ParSettings(20).SetLimits(10,20); //E22 limits
fitter.Config().ParSettings(21).SetLimits(0,10); //S22 limits
fitter.Config().ParSettings(22).SetLimits(0,10); //G31 limits
fitter.Config().ParSettings(23).SetLimits(10,20); //E31 limits
fitter.Config().ParSettings(24).SetLimits(0,10); //S31 limits
fitter.Config().ParSettings(25).Fix();//Fix temperature in second spectra
//fitter.Config().ParSettings(3).SetLimits(0,10000);
// fitter.Config().ParSettings(3).SetStepSize(5);
fitter.Config().MinimizerOptions().SetPrintLevel(0);
fitter.Config().SetMinimizer(“SIMPLEX”);
// fit FCN function directly
// (specify optionally data size and flag to indicate that is a chi2 fit)
fitter.FitFCN(26,globalChi2,0,dataB.Size()+dataSB.Size(),true);
ROOT::Fit::FitResult result = fitter.Result();
result.Print(std::cout);
TCanvas * c1 = new TCanvas(“Simfit”,“Simultaneous fit of two spectras”,
10,10,700,700);
c1->Divide(1,2);
c1->cd(1);
gStyle->SetOptFit(1111);
f1->SetFitResult( result, First_spec_pars);
f1->SetRange(rangeB().first, rangeB().second);
f1->SetLineColor(kBlue);
spectra1->GetListOfFunctions()->Add(f1);
spectra1->SetMarkerStyle(22);
spectra1->Draw(“AP”);
c1->cd(2);
f2->SetFitResult( result, Second_spec_pars);
f2->SetRange(rangeSB().first, rangeSB().second);
f2->SetLineColor(kRed);
spectra2->GetListOfFunctions()->Add(f2);
spectra2->SetMarkerStyle(22);
spectra2->Draw(“AP”);
}
in terminal I have folowing:
root [0] .x combinedFit_Stas_Stas.C
warning: expression result unused [-Wunused-value]
warning: expression result unused [-Wunused-value]
Info in ROOT::Math::ParameterSettings>: lower/upper bounds outside current parameter value. The value will be set to (low+up)/2
Warning in ROOT::Math::FitConfig::CreateMinimizer: Could not create the SIMPLEX minimizer. Try using the minimizer Minuit
Minimizer is Minuit / Migrad
Chi2 = 2.23316e+306
NDf = 298
Edm = 0
NCalls = 2831
Par_0 = 3.75329 +/- 2.72263 (limited)
Par_1 = 0 +/- 0.746941 (limited)
Par_2 = 180 +/- 163.395 (limited)
Par_3 = 35 +/- 32.8639 (limited)
Par_4 = 4 +/- 2.40947 (limited)
Par_5 = 1 +/- 2.4334 (limited)
Par_6 = 100 +/- 120.474 (limited)
Par_7 = 14 +/- 8.16976 (limited)
Par_8 = 1 +/- 2.4334 (limited)
Par_9 = 1 +/- 5.20646 (limited)
Par_10 = 20 +/- 7.11496 (limited)
Par_11 = 1 +/- 5.20646 (limited)
Par_12 = 440 (fixed)
Par_13 = 4 +/- 2.84598 (limited)
Par_14 = 0 +/- 0.746941 (limited)
Par_15 = 180 +/- 163.395 (limited)
Par_16 = 35 +/- 32.8639 (limited)
Par_17 = 4 +/- 2.40947 (limited)
Par_18 = 1 +/- 2.4334 (limited)
Par_19 = 5 +/- 7.46941 (limited)
Par_20 = 14 +/- 6.99749 (limited)
Par_21 = 1 +/- 5.20646 (limited)
Par_22 = 1 +/- 5.20646 (limited)
Par_23 = 20 +/- 7.11496 (limited)
Par_24 = 1 +/- 5.20646 (limited)
Par_25 = 440 (fixed)
What I do wrong?
Thank you