Hello, everyone
I tried to adapt combinedFit.C script for fiiting to datasets of [X1,Y1], [X2,Y2] with 2 functions with some common parameters.
I have changed original script this way
#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 iparB[10] = { 0,1,2,3,4,5,6,7,8,9};
// signal + background function
int iparSB[8] = { 10,11,12,13,14,15,16,17};
// Create the GlobalCHi2 structure
struct GlobalChi2 {
GlobalChi2( ROOT::Math::IMultiGenFunction & f1,
ROOT::Math::IMultiGenFunction & f2) :
fChi2_1(&f1), fChi2_2(&f2) {}
// parameter vector is first background (in common 1 and 2)
// and then is signal (only in 2)
double operator() (const double *par) const {
double p1[10];
for (int i = 0; i < 10; ++i) p1[i] = par[iparB[i] ];
double p2[8];
for (int i = 0; i < 8; ++i) p2[i] = par[iparSB[i] ];
return (*fChi2_1)(p1) + (*fChi2_2)(p2);
}
const ROOT::Math::IMultiGenFunction * fChi2_1;
const ROOT::Math::IMultiGenFunction * fChi2_2;
};
void combinedFit_Stas() {
Sp1=new TGraph("/home/stanislav/Minuit_fit/PZT_2p4_2s_28_500K_a06_minuit");
Sp2=new TGraph("/home/stanislav/Minuit_fit/PZT_2p4_2s_41_440K_a06_minuit");
//TH1D * hB = new TH1D("hB","histo B",100,0,100);
//TH1D * hSB = new TH1D("hSB","histo S+B",100, 0,100);
//First specrtra
//0-Gcp1
//1-x01
//2-Acp1
//3-G11
//4-E11
//5-S11
//6-G21
//7-E21
//8-S21
//9-T1
//Second spectra
//10-Gcp2
//11-x02
//12-Acp2
//13-G12
//14-E12
//5-S12
//15-G22
//16-E22
//17-S22
//18-T2
//TF1 * fB = new TF1("fB","expo",0,100);
TF1 * fB = new TF1("fB","[2]*1/3.141592* [0] / ( (x-[1])*(x-[1]) + [0]*[0])+ [5]*x/(1-exp( -9.1*x/[9] ))*[3]/((x*x-[4]*[4])*(x*x-[4]*[4])+x*x*[3]+ [8]*x/(1-exp( -9.1*x/[9] ))*[6]/((x*x-[7]*[7])*(x*x-[7]*[7])+x*x*[6]))",-20,20);
fB->SetParameters(0.75,0,70, 10,4,1, 30,9,5, 500);
//hB->FillRandom("fB");
//TF1 * fS = new TF1("fS","gaus",0,100);
TF1 * fSB = new TF1("fSB","[12]*1/3.141592* [10] / ( (x-[11])*(x-[11]) + [10]*[10])+ [5]*x/(1-exp( -9.1*x/[17] ))*[13]/((x*x-[14]*[14])*(x*x-[14]*[14])+x*x*[13]+ [8]*x/(1-exp( -9.1*x/[17] ))*[15]/((x*x-[16]*[16])*(x*x-[16]*[16])+x*x*[15]))",-20,20);
fSB->SetParameters(0.75,0,70, 10,4,1, 30,9,5, 440);
//hSB->FillRandom("fB",2000);
//hSB->FillRandom("fS",1000);
// perform now global fit
//hB->Sp1
//hSB->Sp2
//TF1 * fSB = new TF1("fSB","expo + gaus(2)",0,100);
ROOT::Math::WrappedMultiTF1 wfB(*fB,1);
ROOT::Math::WrappedMultiTF1 wfSB(*fSB,1);
ROOT::Fit::DataOptions opt;
ROOT::Fit::DataRange rangeB;
// set the data range
rangeB.SetRange(-20,20);
ROOT::Fit::UnBinData dataB(opt,rangeB);
ROOT::Fit::FillData(dataB, Sp1);
ROOT::Fit::DataRange rangeSB;
rangeSB.SetRange(-20,20);
ROOT::Fit::UnBinData dataSB(opt,rangeSB);
ROOT::Fit::FillData(dataSB, Sp2);
ROOT::Fit::Chi2Function chi2_B(dataB, wfB);
ROOT::Fit::Chi2Function chi2_SB(dataSB, wfSB);
GlobalChi2 globalChi2(chi2_B, chi2_SB);
ROOT::Fit::Fitter fitter;
const int Npar = 18;
double par0[Npar] = {0.75,0,70, 10,4,1, 30,9,5, 500, 0.75,0,1000, 10,4, 30,9, 440};
// create before the parameter settings in order to fix or set range on them
fitter.Config().SetParamsSettings(18,par0);
// set limits for parameters for first spectra
//Cp
fitter.Config().ParSettings(0).Fix();//Gcp1
fitter.Config().ParSettings(1).SetLimits(-0.1,0.1);//x01
fitter.Config().ParSettings(2).SetLimits(0,500);//Acp1
//Phonon1
fitter.Config().ParSettings(3).Fix();//G11
fitter.Config().ParSettings(4).SetLimits(0,7);//E11
fitter.Config().ParSettings(5).SetLimits(0,500);//S11 //common parameter
//Phonon2
fitter.Config().ParSettings(6).Fix();//G21
fitter.Config().ParSettings(7).SetLimits(7,20);//E21
fitter.Config().ParSettings(8).SetLimits(0,500);//S21 //common parameter
//Temperature
fitter.Config().ParSettings(9).Fix();
// set limits for parameters for second spectra
//Cp
fitter.Config().ParSettings(10).Fix();//Gcp2
fitter.Config().ParSettings(11).SetLimits(-0.1,0.1);//x02
fitter.Config().ParSettings(12).SetLimits(0,1000);//Acp2
//Phonon1
fitter.Config().ParSettings(13).Fix();//G12
fitter.Config().ParSettings(14).SetLimits(0,7);//E12
fitter.Config().ParSettings(5).SetLimits(0,500);//S12 //common parameter
//Phonon2
fitter.Config().ParSettings(15).Fix();//G22
fitter.Config().ParSettings(16).SetLimits(7,20);//E22
fitter.Config().ParSettings(8).SetLimits(0,500);//S22 //common parameter
//Temperature
fitter.Config().ParSettings(17).Fix();
//fitter.Config().ParSettings(3).SetStepSize(5);
fitter.Config().MinimizerOptions().SetPrintLevel(5);
fitter.Config().SetMinimizer("GSL");
// fit FCN function directly
// (specify optionally data size and flag to indicate that is a chi2 fit)
fitter.FitFCN(18,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);
fB->SetFitResult( result, iparB);
fB->SetRange(rangeB().first, rangeB().second);
fB->SetLineColor(kBlue);
Sp1->GetListOfFunctions()->Add(fB);
Sp1->Draw("scat");
c1->cd(2);
fSB->SetFitResult( result, iparSB);
fSB->SetRange(rangeSB().first, rangeSB().second);
fSB->SetLineColor(kRed);
Sp2->GetListOfFunctions()->Add(fSB);
Sp2->Draw("scat");
}
Common parameters 5 and 8
Dataset looks like txt file
-20.75 0.411825
-20.5 0.274575
-20.25 0.159675
-20 0.119325
-19.75 0.479625
-19.5 0.556875
-19.25 0.314925
-19 0.867
-18.75 0.1764
-18.5 0.62835
-18.25 0.39165
-18 0.29595
-17.75 0.471825
-17.5 0.415125
-17.25 0.552075
-17 0.87
-16.75 0.649875
-16.5 1.18575
-16.25 0.888
-16 0.5502
-15.75 0.78825
-15.5 0.532425
-15.25 0.554325
-15 1.002
-14.75 0.54765
-14.5 0.6498
-14.25 0.792
-14 1.0215
-13.75 1.00725
-13.5 0.945
-10 2.4405
-9.75 2.9715
-9.5 2.2905
-9.25 3.3
-9 2.316
-8.75 2.5545
-8.5 3.399
-8.25 2.154
-8 2.97525
-7.75 2.63325
-7.5 2.586
-7.25 2.88
-7 3.7845
-6.75 2.90475
-6.5 4.017
-6.25 3.921
-6 3.62925
-5.75 4.02825
-5.5 4.524
-5.25 3.43575
-5 4.503
-4.75 5.71125
-4.5 6.1905
-4.25 6.22425
-4 6.2295
-3.75 7.695
-3.5 9.2025
-3.25 10.86
-3 13.395
-2.75 13.56
-2.5 12.93
-2.25 16.86
-2 16.0425
-1.75 18.615
-1.5 19.8075
-1.25 18.9075
-1 17.865
-0.75 18.51
-0.5 20.3325
-0.25 20.2275
0 23.0025
0.25 23.745
0.5 21.4425
0.75 21.93
1 18.3975
1.25 20.49
1.5 18.96
1.75 19.92
2 19.4475
2.25 17.655
2.5 14.4225
2.75 15.525
3 13.245
3.25 13.785
3.5 11.49
3.75 10.3875
4 8.8575
4.25 7.98
4.5 6.61725
4.75 6.0975
5 5.72775
5.25 6.04725
5.5 4.83225
5.75 4.41375
6 3.5895
6.25 3.885
6.5 3.95325
6.75 3.10125
7 3.33825
7.25 3.6165
7.5 2.727
7.75 3.31275
8 3.41175
8.25 3.4215
8.5 3.108
8.75 3.50475
9 2.82525
9.25 2.10975
9.5 3.1695
9.75 3.6555
10 2.328
10.25 2.9385
10.5 2.27325
10.75 1.866
11 1.866
11.25 1.9245
11.5 2.25825
11.75 1.7025
12 2.16525
12.25 1.944
12.5 1.3215
12.75 2.45175
13 2.01525
13.25 1.395
13.5 1.6335
13.75 1.85475
14 1.77675
14.25 1.54275
14.5 1.392
14.75 1.70325
15 1.4685
15.25 0.77175
15.5 1.23075
15.75 0.84375
16 0.9855
16.25 1.2285
16.5 0.846
16.75 0.693075
17 0.6174
17.25 0.92025
17.5 0.9225
17.75 0.1539
18 0.76875
18.25 0.9285
18.5 0.69315
18.75 0.63345
19 0.386025
19.25 0.5202
19.5 0.6177
19.75 0.508125
20 0.6897
And as a result I got this
Help me please!