Hello everyone,
I’m trying to adapt CombinedFit.C (from ROOT examples) to solve my problem of optimisation. I want to fit two data sets (as TH2) to two different functions that have two parameters, both shared between the two functions, by minimisation of chi squared. I found that the optimisation problem that CombinedFit.C solves is like mine, so I started there. However, the fit does not result in reasonable values, and I could not identify if the data from the TH2 isn’t being loaded properly or something else. I was able to take a profile using ProfileX from my data and do the fit properly by loading the data as TGraph, but the profile limits the data that is fitted.
I wanted to know how can this type of fit using Fit::Fitter and data from a TH2 for this project and future ones. I’m adding my code as text here, since the Forum is not letting me upload neither the code nor the auxiliary files for running my code. I hope I can upload the files soon.
Thanks in advance.
ROOT Version: 6.28/06
Platform: Ubuntu 22.04.5 LTS
Compiler: gcc 11.4.0
// Include standard stuff
#include <cstdlib>
#include <cwchar>
#include <fstream>
#include <iostream>
#include <sstream>
#include <stdio.h>
#include <stdlib.h>
#include <string>
#include <sys/stat.h>
// Include ROOT stuff
#include "Math/IFunction.h"
#include "TCanvas.h"
#include "TClonesArray.h"
#include "TCutG.h"
#include "TDirectory.h"
#include "TF1.h"
#include "TF2.h"
#include "TFile.h"
#include "TGraph.h"
#include "TH1.h"
#include "TH2.h"
#include "TLatex.h"
#include "TObject.h"
#include "TPDF.h"
#include "TProfile.h"
#include "TRatioPlot.h"
#include "TTree.h"
#include <Fit/BinData.h>
#include <Fit/Chi2FCN.h>
#include <Fit/Fitter.h>
#include <HFitInterface.h>
#include <Math/WrappedMultiTF1.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TStyle.h>
using namespace std;
TPDF pdff("result.pdf");
TFile *outfile = new TFile("data.root", "READONLY");
TGraph *simul1;
Double_t InterpolatedGraph(Double_t *x, Double_t *par) {
return par[0] + par[1] * (simul1->Eval(x[0]));
};
TGraph *simul2;
Double_t InterpolatedGraph2(Double_t *x, Double_t *par) {
return par[0] + par[1] * (simul2->Eval(x[0]));
};
int ipar[2] = {
0, // intercept common parameter
1 // slope common parameter
};
// 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 p[2];
for (int i = 0; i < 2; ++i)
p[i] = par[ipar[i]];
return (*fChi2_1)(p) + (*fChi2_2)(p);
}
const ROOT::Math::IMultiGenFunction *fChi2_1;
const ROOT::Math::IMultiGenFunction *fChi2_2;
};
char filename[4];
int combinedFitTH2() {
TH2F *htemp1 = (TH2F *)outfile->Get("htemp1");
TH2F *htemp2 = (TH2F *)outfile->Get("htemp2");
TCanvas *canvas = new TCanvas("canvas");
canvas->Print("result.pdf[");
canvas->SetLogz();
simul1 = new TGraph("f1.txt", "%lg %lg");
TF1 *f1 = new TF1("f1", InterpolatedGraph, 0, 200, 2);
f1->SetParameters(0, 1);
f1->SetRange(0, 600, 70, 1400);
f1->SetParName(0, "intercept");
f1->SetParName(1, "slope");
simul2 = new TGraph("f2.txt", "%lg %lg");
TF1 *f2 = new TF1("f2", InterpolatedGraph2, 0, 20, 2);
f2->SetParameters(0, 1);
f2->SetRange(0, 0, 20, 200);
f2->SetParName(0, "intercept");
f2->SetParName(1, "slope");
ROOT::Math::WrappedMultiTF1 wf1(*f1, 0);
ROOT::Math::WrappedMultiTF1 wf2(*f2, 0);
ROOT::Fit::DataOptions opt;
// set the data range
ROOT::Fit::DataRange range1;
range1.SetRange(0, 70, 600, 1400);
ROOT::Fit::BinData data1(opt, range1);
ROOT::Fit::FillData(data1, htemp1);
ROOT::Fit::DataRange range2;
range2.SetRange(0, 15, 0, 200);
ROOT::Fit::BinData data2(opt, range2);
ROOT::Fit::FillData(data2, htemp2);
ROOT::Fit::Chi2Function chi2_f1(data1, wf1);
ROOT::Fit::Chi2Function chi2_f2(data2, wf2);
GlobalChi2 globalChi2(chi2_f1, chi2_f2);
ROOT::Fit::Fitter fitter;
const int Npar = 2;
double par0[Npar] = {0, 1};
// create before the parameter settings in order to fix or set range on them
fitter.Config().SetParamsSettings(Npar, par0);
// fitter.Config().ParSettings(0).SetLimits(0,200);
// fitter.Config().ParSettings(1).SetLimits(0,100);
fitter.Config().MinimizerOptions().SetPrintLevel(0);
fitter.Config().SetMinimizer("Minuit2", "Migrad");
// fit FCN function directly
// (specify optionally data size and flag to indicate that is a chi2 fit)
fitter.FitFCN(Npar, globalChi2, par0, data1.Size() + data2.Size(), true);
ROOT::Fit::FitResult result = fitter.Result();
result.Print(std::cout);
std::cout << "Combined fit Chi2 = " << result.Chi2() << std::endl;
// f1->SetFitResult(result, ipar);
f1->SetLineColor(kBlue);
f2->SetLineColor(kRed);
// hSi->GetListOfFunctions()->Add(f1);
// hSi->Draw();
// f1->Draw("same");
htemp1->Draw("colz");
htemp2->Draw("samecolz");
f1->Draw("same");
f2->Draw("same");
TLatex *tex3 = new TLatex(0, 713.64,
Form("Blue: a*f(x) + b, a = %.2f b = %.2f",
f1->GetParameter(1), f1->GetParameter(0)));
tex3->SetTextSize(0.03);
tex3->SetLineWidth(1);
tex3->Draw();
TLatex *tex1 = new TLatex(0, 753.64,
Form("Red: a*f(x) + b, a = %.2f b = %.2f",
f2->GetParameter(1), f2->GetParameter(0)));
tex1->SetTextSize(0.03);
tex1->SetLineWidth(1);
tex1->Draw();
canvas->Update();
canvas->Print("result.pdf", "");
canvas->Print("result.pdf]");
return 0;
}
