Thank you again. Now the macro seems to work quite nicely. I still have some problems with the chi2 I get from the fit, which is much higher than the chi2 I calculate “manually” ( (y_i - f(a_i))^2 / ey_i^2 ).
I report here the main part of the script, so that perhaps when you have some time you could have a look if there is any mistake. To briefly summarize, I have to fit 8 different sets of data with a function of 4 parameters, two of them are fixed for each set, the other two have to be globally fitted.
#include <stdio.h>
#include <stdarg.h>
#include <math.h>
#include <iostream>
#include <fstream>
#include <ostream>
#include <sstream>
#include "TMath.h"
#include "TROOT.h"
#include "Math/DistFunc.h"
#include "TCanvas.h"
#include "TH1D.h"
#include "TStyle.h"
#include "TF1.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TLegend.h"
#include "TFrame.h"
#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 "Fit/FitResult.h"
struct GlobalChi2 {
GlobalChi2( ROOT::Math::IMultiGenFunction & f1, ROOT::Math::IMultiGenFunction & f2, ROOT::Math::IMultiGenFunction &f3, ROOT::Math::IMultiGenFunction &f4, ROOT::Math::IMultiGenFunction &f5, ROOT::Math::IMultiGenFunction &f6, ROOT::Math::IMultiGenFunction &f7, ROOT::Math::IMultiGenFunction &f8) :
fChi2_1(&f1), fChi2_2(&f2), fChi2_3(&f3), fChi2_4(&f4), fChi2_5(&f5), fChi2_6(&f6), fChi2_7(&f7), fChi2_8(&f8) {}
double operator() (const double *par) const {
double p1[4];
p1[0] = par[0];
p1[1] = par[1];
p1[2] = par[2];
p1[3] = par[3];
double p2[4];
p2[0] = par[4];
p2[1] = par[5];
p2[2] = par[2];
p2[3] = par[3];
double p3[4];
p3[0] = par[6];
p3[1] = par[7];
p3[2] = par[2];
p3[3] = par[3];
double p4[4];
p4[0] = par[8];
p4[1] = par[9];
p4[2] = par[2];
p4[3] = par[3];
double p5[4];
p5[0] = par[10];
p5[1] = par[11];
p5[2] = par[2];
p5[3] = par[3];
double p6[4];
p6[0] = par[12];
p6[1] = par[13];
p6[2] = par[2];
p6[3] = par[3];
double p7[4];
p7[0] = par[14];
p7[1] = par[15];
p7[2] = par[2];
p7[3] = par[3];
double p8[4];
p8[0] = par[16];
p8[1] = par[17];
p8[2] = par[2];
p8[3] = par[3];
return (*fChi2_1)(p1) + (*fChi2_2)(p2) + (*fChi2_3)(p3) + (*fChi2_4)(p4) + (*fChi2_5)(p5) + (*fChi2_6)(p6) + (*fChi2_7)(p7) + (*fChi2_8)(p8);
}
const ROOT::Math::IMultiGenFunction * fChi2_1;
const ROOT::Math::IMultiGenFunction * fChi2_2;
const ROOT::Math::IMultiGenFunction * fChi2_3;
const ROOT::Math::IMultiGenFunction * fChi2_4;
const ROOT::Math::IMultiGenFunction * fChi2_5;
const ROOT::Math::IMultiGenFunction * fChi2_6;
const ROOT::Math::IMultiGenFunction * fChi2_7;
const ROOT::Math::IMultiGenFunction * fChi2_8;
};
int CombFit_mod( ) {
gROOT -> SetStyle("Plain");
gStyle -> SetOptFit(1111111);
gStyle -> SetOptStat(111111);
/*
....input data and TGraphErrors.....
*/
TH1D *h1 = new TH1D("h1","h1",300, 0., 300.);
TH1D *h2 = new TH1D("h2","h2",300, 0., 300.);
TH1D *h3 = new TH1D("h3","h3",300, 0., 300.);
TH1D *h4 = new TH1D("h4","h4",300, 0., 300.);
TH1D *h5 = new TH1D("h5","h5",300, 0., 300.);
TH1D *h6 = new TH1D("h6","h6",300, 0., 300.);
TH1D *h7 = new TH1D("h7","h7",300, 0., 300.);
TH1D *h8 = new TH1D("h8","h8",300, 0., 300.);
TF1 *f1 = new TF1("f1","([0]*exp(-(x)/[2]) + [1]*exp(-(x-1.0)/[3])+[4])",0.,180.);
TF1 *f2 = new TF1("f2","([2]*exp(-(x)/[0]) + [3]*exp(-(x-1.0)/[1])+[4])",0.,180.);
TF1 *f3 = new TF1("f3","([2]*exp(-(x)/[0]) + [3]*exp(-(x-1.0)/[1])+[4])",0.,180.);
TF1 *f4 = new TF1("f4","([2]*exp(-(x)/[0]) + [3]*exp(-(x-1.0)/[1])+[4])",0.,180.);
TF1 *f5 = new TF1("f5","([2]*exp(-(x)/[0]) + [3]*exp(-(x-1.0)/[1])+[4])",0.,180.);
TF1 *f6 = new TF1("f6","([2]*exp(-(x)/[0]) + [3]*exp(-(x-1.0)/[1])+[4])",0.,180.);
TF1 *f7 = new TF1("f7","([2]*exp(-(x)/[0]) + [3]*exp(-(x-1.0)/[1])+[4])",0.,180.);
TF1 *f8 = new TF1("f8","([2]*exp(-(x)/[0]) + [3]*exp(-(x-1.0)/[1])+[4])",0.,180.);
// perform now global fit
ROOT::Math::WrappedMultiTF1 wf1(*f1,1);
ROOT::Math::WrappedMultiTF1 wf2(*f2,1);
ROOT::Math::WrappedMultiTF1 wf3(*f3,1);
ROOT::Math::WrappedMultiTF1 wf4(*f4,1);
ROOT::Math::WrappedMultiTF1 wf5(*f5,1);
ROOT::Math::WrappedMultiTF1 wf6(*f6,1);
ROOT::Math::WrappedMultiTF1 wf7(*f7,1);
ROOT::Math::WrappedMultiTF1 wf8(*f8,1);
ROOT::Fit::DataOptions opt;
ROOT::Fit::DataRange rangeD;
rangeD.SetRange(0.,180.);
ROOT::Fit::BinData data1(opt,rangeD);
ROOT::Fit::BinData data2(opt,rangeD);
ROOT::Fit::BinData data3(opt,rangeD);
ROOT::Fit::BinData data4(opt,rangeD);
ROOT::Fit::BinData data5(opt,rangeD);
ROOT::Fit::BinData data6(opt,rangeD);
ROOT::Fit::BinData data7(opt,rangeD);
ROOT::Fit::BinData data8(opt,rangeD);
ROOT::Fit::FillData(data1, expH_11);
ROOT::Fit::FillData(data2, expH_31);
ROOT::Fit::FillData(data3, expH_13);
ROOT::Fit::FillData(data4, expH_62);
ROOT::Fit::FillData(data5, expHe3_53);
ROOT::Fit::FillData(data6, expHe3_81);
ROOT::Fit::FillData(data7, expHe4_53);
ROOT::Fit::FillData(data8, expHe4_123);
ROOT::Fit::Chi2Function chi2_1(data1, wf1);
ROOT::Fit::Chi2Function chi2_2(data2, wf2);
ROOT::Fit::Chi2Function chi2_3(data3, wf3);
ROOT::Fit::Chi2Function chi2_4(data4, wf4);
ROOT::Fit::Chi2Function chi2_5(data5, wf5);
ROOT::Fit::Chi2Function chi2_6(data6, wf6);
ROOT::Fit::Chi2Function chi2_7(data7, wf7);
ROOT::Fit::Chi2Function chi2_8(data8, wf8);
GlobalChi2 globalChi2(chi2_1, chi2_2, chi2_3, chi2_4, chi2_5, chi2_6, chi2_7, chi2_8);
ROOT::Fit::Fitter fitter;
const int Npar = 18;
double par0[Npar] = {0.60, 0.40, 15.0, 120.0, 0.40, 0.60, 0.58, 0.42, 0.17, 0.83, 0.35, 0.65, 0.19, 0.81, 0.347, 0.653, 0.065, 0.935};
fitter.Config().SetParamsSettings(18,par0);
fitter.Config().ParSettings(0).Fix();
fitter.Config().ParSettings(1).Fix();
fitter.Config().ParSettings(4).Fix();
fitter.Config().ParSettings(5).Fix();
fitter.Config().ParSettings(6).Fix();
fitter.Config().ParSettings(7).Fix();
fitter.Config().ParSettings(8).Fix();
fitter.Config().ParSettings(9).Fix();
fitter.Config().ParSettings(10).Fix();
fitter.Config().ParSettings(11).Fix();
fitter.Config().ParSettings(12).Fix();
fitter.Config().ParSettings(13).Fix();
fitter.Config().ParSettings(14).Fix();
fitter.Config().ParSettings(15).Fix();
fitter.Config().ParSettings(16).Fix();
fitter.Config().ParSettings(17).Fix();
fitter.FitFCN(18,globalChi2,par0,data1.Size()+data2.Size()+data3.Size()+data4.Size()+data5.Size()+data6.Size()+data7.Size()+data8.Size(), true);
fitter.Config().MinimizerOptions().SetPrintLevel(1);
ROOT::Fit::FitResult result = fitter.Result();
result.Print(std::cout);
TCanvas * c1 = new TCanvas("Simultaneous fit");
c1->Divide(4,2);
c1->cd(1);
f1->SetParameters(result.Parameter(0), result.Parameter(1), result.Parameter(2), result.Parameter(3));
f1->SetParError(0, result.Error(0) );
f1->SetParError(1, result.Error(1) );
f1->SetParError(2, result.Error(2) );
f1->SetParError(3, result.Error(3) );
f1->SetChisquare(result.MinFcnValue() );
f1->SetNDF(result.Ndf() );
h1->GetListOfFunctions()->Add(f1);
f1->SetRange(rangeD().first, rangeD().second);
f1->DrawClone("l same");
c1->cd(2);
f2->SetParameters(result.Parameter(2), result.Parameter(3), result.Parameter(4), result.Parameter(5));
f2->SetParError(2, result.Error(2) );
f2->SetParError(3, result.Error(3) );
f2->SetParError(4, result.Error(4) );
f2->SetParError(5, result.Error(5) );
f2->SetChisquare(result.MinFcnValue() );
f2->SetNDF(result.Ndf() );
h2->GetListOfFunctions()->Add(f2);
f2->SetRange(rangeD().first, rangeD().second);
//h2->Draw("l");
f2->DrawClone("l same");
c1->cd(3);
f3->SetParameters(result.Parameter(2), result.Parameter(3), result.Parameter(6), result.Parameter(7));
f3->SetParError(2, result.Error(2) );
f3->SetParError(3, result.Error(3) );
f3->SetParError(6, result.Error(6) );
f3->SetParError(7, result.Error(7) );
f3->SetChisquare(result.MinFcnValue() );
f3->SetNDF(result.Ndf() );
h3->GetListOfFunctions()->Add(f3);
f3->SetRange(rangeD().first, rangeD().second);
f3->DrawClone("l same");
c1->cd(4);
f4->SetParameters(result.Parameter(2), result.Parameter(3), result.Parameter(8), result.Parameter(9));
f4->SetParError(2, result.Error(2) );
f4->SetParError(3, result.Error(3) );
f4->SetParError(8, result.Error(8) );
f4->SetParError(9, result.Error(9) );
f4->SetChisquare(result.MinFcnValue() );
f4->SetNDF(result.Ndf() );
h4->GetListOfFunctions()->Add(f4);
f4->SetRange(rangeD().first, rangeD().second);
f4->DrawClone("l same");
c1->cd(5);
f5->SetParameters(result.Parameter(2), result.Parameter(3), result.Parameter(10), result.Parameter(11));
f5->SetParError(2, result.Error(2) );
f5->SetParError(3, result.Error(3) );
f5->SetParError(10, result.Error(10) );
f5->SetParError(11, result.Error(11) );
f5->SetChisquare(result.MinFcnValue() );
f5->SetNDF(result.Ndf() );
h5->GetListOfFunctions()->Add(f5);
f5->SetRange(rangeD().first, rangeD().second);
f5->DrawClone("l same");
c1->cd(6);
f6->SetParameters(result.Parameter(2), result.Parameter(3), result.Parameter(12), result.Parameter(13));
f6->SetParError(2, result.Error(2) );
f6->SetParError(3, result.Error(3) );
f6->SetParError(12, result.Error(12) );
f6->SetParError(13, result.Error(13) );
f6->SetChisquare(result.MinFcnValue() );
f6->SetNDF(result.Ndf() );
h6->GetListOfFunctions()->Add(f6);
f6->SetRange(rangeD().first, rangeD().second);
f6->DrawClone("l same");
c1->cd(7);
f7->SetParameters(result.Parameter(2), result.Parameter(3), result.Parameter(14), result.Parameter(15));
f7->SetParError(2, result.Error(2) );
f7->SetParError(3, result.Error(3) );
f7->SetParError(14, result.Error(14) );
f7->SetParError(15, result.Error(15) );
f7->SetChisquare(result.MinFcnValue() );
f7->SetNDF(result.Ndf() );
h7->GetListOfFunctions()->Add(f7);
f7->SetRange(rangeD().first, rangeD().second);
f7->DrawClone("l same");
c1->cd(8);
f8->SetParameters(result.Parameter(2), result.Parameter(3), result.Parameter(16), result.Parameter(17));
f8->SetParError(2, result.Error(2) );
f8->SetParError(3, result.Error(3) );
f8->SetParError(16, result.Error(16) );
f8->SetParError(17, result.Error(17) );
f8->SetChisquare(result.MinFcnValue() );
f8->SetNDF(result.Ndf() );
h8->GetListOfFunctions()->Add(f7);
f8->SetRange(rangeD().first, rangeD().second);
f8->DrawClone("l same");
return 0;
}
Thanks a lot for your support.
Best regards