|
|
//+ Combined (simultaneous) fit of two histogram with separate functions
// and some common parameters
//
// See http://root.cern.ch/phpBB3//viewtopic.php?f=3&t=11740#p50908
// for a modified version working with Fumili or GSLMultiFit
//
// N.B. this macro must be compiled with ACliC
//
//Author: L. Moneta - Dec 2010
#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"
#include "TGraphErrors.h"
#include "TLegend.h"
#include "TMultiGraph.h"
// definition of shared parameter
// background function
int iparB[3] = { 0, // exp amplitude in B histo
1, // T_0
2 // common parameter
};
// signal + background function
int iparSB2[3] = { 0, // amplitude in S+B histo
3, // T_0
2, // common parameter
};
int iparSB3[3] = { 0, // amplitude
4, // T_0
2, // common parameter
};
int iparSB4[3] = { 0, // amplitude
5, // T_0
2, // common parameter
};
struct GlobalChi2 { GlobalChi2( ROOT::Math::IMultiGenFunction & f1, ROOT::Math::IMultiGenFunction & f2, ROOT::Math::IMultiGenFunction & f3, ROOT::Math::IMultiGenFunction & f4) : fChi2_1(&f1), fChi2_2(&f2), fChi2_3(&f3), fChi2_4(&f4) {}
// 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[3]; for (int i = 0; i < 3; ++i) p1[i] = par[iparB[i] ];
double p2[3]; for (int i = 0; i < 3; ++i) p2[i] = par[iparSB2[i] ];
double p3[3]; for (int i = 0; i < 3; ++i) p3[i] = par[iparSB3[i] ];
double p4[3]; for (int i = 0; i < 3; ++i) p4[i] = par[iparSB4[i] ];
return (*fChi2_1)(p1) + (*fChi2_2)(p2) + (*fChi2_3)(p3) + (*fChi2_4)(p4); }
const ROOT::Math::IMultiGenFunction * fChi2_1; const ROOT::Math::IMultiGenFunction * fChi2_2; const ROOT::Math::IMultiGenFunction * fChi2_3; const ROOT::Math::IMultiGenFunction * fChi2_4;
};
void combinedFit_hit2() {
#if defined(__CINT__) && !defined(__MAKECINT__)
cout << "ERROR: This tutorial can run only using ACliC, you must run it by doing: " << endl; cout << "\t .x $ROOTSYS/tutorials/fit/combinedFit_hit.C+" << endl; return;#endif
TF1 * fB = new TF1("fB","[0]*( (1- 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x )) / (1+[2]* (1- 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x ))*1.0*1.65901*TMath::Power(x,-1.7218)) + 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x ) ) ",0,1);//energy units in keV
//TF1 * fB = new TF1("fB","[0]*( (1-[1]) / (1+[2]* (1- [1])*1.0*1.65901*TMath::Power(x,-1.7218)) +[1] ) ",0,1);
TF1 * fSB2 = new TF1("fSB2","[0]*( (1- 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x )) / (1+[2]* (1- 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x ))*4.0*1.65901*TMath::Power(x,-1.7218)) + 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x ) ) ",0,1);//energy units in keV
//TF1 * fSB2 = new TF1("fSB2","[0]*( (1-[1]) / (1+[2]* (1- [1])*4.0*1.65901*TMath::Power(x,-1.7218)) +[1] ) ",0,1);
TF1 * fSB3 = new TF1("fSB3","[0]*( (1- 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x )) / (1+[2]* (1- 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x ))*36.0*1.65901*TMath::Power(x,-1.7218)) + 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x ) ) ",0,1);//energy units in keV
//TF1 * fSB3 = new TF1("fSB3","[0]*( (1-[1]) / (1+[2]* (1- [1])*36.0*1.65901*TMath::Power(x,-1.7218)) +[1] ) ",0,1);
TF1 * fSB4 = new TF1("fSB4","[0]*( (1- 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x )) / (1+[2]* (1- 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x ))*64.0*1.65901*TMath::Power(x,-1.7218)) + 0.5*(TMath::Log(2*0.511E6*x*x/(1-x*x)/[1]) -x*x )/(TMath::Log(2*0.511E6*x*x/(1-x*x)/68.7) -x*x ) ) ",0,1); //energy units in keV
// TF1 * fSB4 = new TF1("fSB4","[0]*( (1-[1]) / (1+[2]* (1- [1])*64.0*1.65901*TMath::Power(x,-1.7218)) +[1] ) ",0,1);
TGraphErrors * hB = new TGraphErrors("energylist_p_bpmbeta1.txt","%lg %lg %lg"); TGraphErrors * hSB2 = new TGraphErrors("energylist_he_bpmbeta1.txt","%lg %lg %lg"); TGraphErrors * hSB3 = new TGraphErrors("energylist_c_bpmbeta1.txt","%lg %lg %lg"); TGraphErrors * hSB4 = new TGraphErrors("energylist_o_bpmbeta1.txt","%lg %lg %lg");
TGraphErrors * hB_2 = new TGraphErrors("energylist_p_bpmbeta1.txt","%lg %lg %lg"); TGraphErrors * hSB2_2 = new TGraphErrors("energylist_he_bpmbeta1.txt","%lg %lg %lg"); TGraphErrors * hSB3_2 = new TGraphErrors("energylist_c_bpmbeta1.txt","%lg %lg %lg"); TGraphErrors * hSB4_2 = new TGraphErrors("energylist_o_bpmbeta1.txt","%lg %lg %lg"); // perform now global fit
// TF1 * fSB2 = new TF1("fSB2","expo + gaus(2)",0,100);
ROOT::Math::WrappedMultiTF1 wfB(*fB,1); ROOT::Math::WrappedMultiTF1 wfSB2(*fSB2,1); ROOT::Math::WrappedMultiTF1 wfSB3(*fSB3,1); ROOT::Math::WrappedMultiTF1 wfSB4(*fSB4,1);
ROOT::Fit::DataOptions opt; ROOT::Fit::DataRange rangeB; // set the data range
rangeB.SetRange(0.3,0.59); ROOT::Fit::BinData dataB(opt,rangeB); ROOT::Fit::FillData(dataB, hB);
ROOT::Fit::DataRange rangeSB2; rangeSB2.SetRange(0.34,0.59); ROOT::Fit::BinData dataSB2(opt,rangeSB2); ROOT::Fit::FillData(dataSB2, hSB2);
ROOT::Fit::DataRange rangeSB3; rangeSB3.SetRange(0.4,0.73); ROOT::Fit::BinData dataSB3(opt,rangeSB3); ROOT::Fit::FillData(dataSB3, hSB3);
ROOT::Fit::DataRange rangeSB4; rangeSB4.SetRange(0.43,0.73); ROOT::Fit::BinData dataSB4(opt,rangeSB4); ROOT::Fit::FillData(dataSB4, hSB4);
ROOT::Fit::Chi2Function chi2_B(dataB, wfB); ROOT::Fit::Chi2Function chi2_SB2(dataSB2, wfSB2); ROOT::Fit::Chi2Function chi2_SB3(dataSB3, wfSB3); ROOT::Fit::Chi2Function chi2_SB4(dataSB4, wfSB4);
GlobalChi2 globalChi2(chi2_B, chi2_SB2, chi2_SB3, chi2_SB4);
ROOT::Fit::Fitter fitter;
const int Npar = 6; double par0[Npar] = {1.4,1E5, 0.005, 1e4,1e4,1e4}; // double par0[Npar] = {1.4,0.000007,0.15,1.4, 0.00250,2, 0.005,2, 0.007};
// create before the parameter settings in order to fix or set range on them
fitter.Config().SetParamsSettings(Npar,par0); // fix 5-th parameter
// fitter.Config().ParSettings(1).Fix();
// set limits on the third and 4-th parameter
fitter.Config().ParSettings(0).SetLimits(0,5); fitter.Config().ParSettings(2).SetLimits(0,1.0); fitter.Config().ParSettings(1).SetLimits(0,2E9); fitter.Config().ParSettings(3).SetLimits(0,2E6);//1
fitter.Config().ParSettings(4).SetLimits(0,2E6);//10
fitter.Config().ParSettings(5).SetLimits(0,2E6);//10
// fitter.Config().ParSettings(4).SetStepSize(2);
// fitter.Config().ParSettings(1).SetStepSize(0.1);
// fitter.Config().ParSettings(6).SetStepSize(10);
// fitter.Config().ParSettings(8).SetStepSize(10);
fitter.Config().MinimizerOptions().SetPrintLevel(0); fitter.Config().SetMinimizer("Minuit2","Migrad"); //Minuit2
//fitter.Config().SetNormErrors(true);
// fit FCN function directly
// (specify optionally data size and flag to indicate that is a chi2 fit)
fitter.FitFCN(Npar,globalChi2,0,dataB.Size()+dataSB2.Size()+dataSB3.Size()+dataSB4.Size(),true); ROOT::Fit::FitResult result = fitter.Result(); result.Print(std::cout, true); //result.PrintCovMatrix(std::cout);
cout << "FitResult.Status() = " << result.Status() << endl;
TCanvas * c1 = new TCanvas("Simfit","Simultaneous fit"); // c1->Divide(2,2);
// c1->cd(1);
// gStyle->SetOptFit(1111);
fB->SetFitResult( result, iparB); fB->SetRange(rangeB().first, rangeB().second); fB->SetLineColor(kRed); hB->GetListOfFunctions()->Add(fB);
//hB->Draw("AP");
// c1->cd(2);
fSB2->SetFitResult( result, iparSB2); fSB2->SetRange(rangeSB2().first, rangeSB2().second); fSB2->SetLineColor(kRed); hSB2->GetListOfFunctions()->Add(fSB2); // hSB2->Draw("AP");
// c1->cd(3);
fSB3->SetFitResult( result, iparSB3); fSB3->SetRange(rangeSB3().first, rangeSB3().second); fSB3->SetLineColor(kRed); hSB3->GetListOfFunctions()->Add(fSB3); // hSB3->Draw("AP");
// c1->cd(4);
fSB4->SetFitResult( result, iparSB4); fSB4->SetRange(rangeSB4().first, rangeSB4().second); fSB4->SetLineColor(kRed); hSB4->GetListOfFunctions()->Add(fSB4); // hSB4->Draw("AP");
TMultiGraph * mg_e9 = new TMultiGraph(); mg_e9->Add(hB,"p"); hB->SetMarkerStyle(20); mg_e9->Add(hSB2,"p"); hSB2->SetMarkerStyle(21); mg_e9->Add(hSB3,"p"); hSB3->SetMarkerStyle(22); mg_e9->Add(hSB4,"p"); hSB4->SetMarkerStyle(23); mg_e9->Draw("a"); //must draw first before labeling axes
mg_e9->SetTitle(" "); mg_e9->GetXaxis()->SetTitle("Lorentz #beta"); mg_e9->GetYaxis()->SetTitle("dA/dE / (10^{-6} a.u./(Mev#upointion#upoints^{-1})"); mg_e9->SetMinimum(0.5); mg_e9->SetMaximum(1.4); TF1 * tf1_birk1 = new TF1("tf1_birk1","[0]*(1/(1+[1]*1.65901*TMath::Power(x,-1.7218)))",0.3,0.59); TF1 * tf1_birk2 = new TF1("tf1_birk2","[0]*(1/(1+[1]*4*1.65901*TMath::Power(x,-1.7218)))",0.34,0.59); TF1 * tf1_birk3 = new TF1("tf1_birk3","[0]*(1/(1+[1]*36*1.65901*TMath::Power(x,-1.7218)))",0.4,0.73); TF1 * tf1_birk4 = new TF1("tf1_birk4","[0]*(1/(1+[1]*64*1.65901*TMath::Power(x,-1.7218)))",0.43,0.73); tf1_birk1->SetParameters(1.46123,0.0177651);tf1_birk1->SetLineStyle(2);tf1_birk1->SetLineColor(kBlack); tf1_birk2->SetParameters(1.46123,0.00682331);tf1_birk2->SetLineStyle(2);tf1_birk2->SetLineColor(kBlack); tf1_birk3->SetParameters(1.46123,0.00389243);tf1_birk3->SetLineStyle(2);tf1_birk3->SetLineColor(kBlack); tf1_birk4->SetParameters(1.46123,0.00310792);tf1_birk4->SetLineStyle(2);tf1_birk4->SetLineColor(kBlack); // hB_2->Fit(tf1_birk1);
// hSB2_2->Fit(tf1_birk2);
//hSB3_2->Fit(tf1_birk3);
//hSB4_2->Fit(tf1_birk4);
tf1_birk1->Draw("same"); tf1_birk2->Draw("same"); tf1_birk3->Draw("same"); tf1_birk4->Draw("same"); TLegend * mylegende9 = new TLegend(0.65,0.65,0.9,0.9); mylegende9->SetFillColor(0); // white background
mylegende9->SetTextFont(22); mylegende9->SetBorderSize(0); // get rid of the box
mylegende9->SetTextSize(0.035); // set text size
mylegende9->AddEntry(hB,"Protons","p"); // options: p,l,f
mylegende9->AddEntry(hSB2,"Helium","p"); // options: p,l,f
mylegende9->AddEntry(hSB3,"Carbon","p"); // options: p,l,f
mylegende9->AddEntry(hSB4,"Oxygen","p"); // options: p,l,f
mylegende9->AddEntry(fB,"BVT Fit","l"); // options: p,l,f
mylegende9->AddEntry(tf1_birk1,"Birks Fit","l"); // options: p,l,f
// mylegende9->AddEntry(tf1_pow1,"[1]#upointpow(#beta,[2])}","l"); // options: p,l,f
mylegende9->Draw(); gPad->Modified(); c1->SaveAs("figs/bvt_beta_combinedfit.pdf"); c1->SaveAs("figs/bvt_beta_combinedfit.png"); c1->SaveAs("figs/bvt_beta_combinedfit.C");
}
|
