Fit multiple histograms with the same function

This kind of syntax:

works only with pre-defined functions. With interpreted function you can use the [], as you are doing. It should work fine.
However in this case, I think you should define the global chi2 function differently. Your both TF1 will have 6 parameters. So you comute the global chi2 just as

without copying the parameter vectors.

The other error you are having is that you need to compile the macro with ACliC.

Lorenzo

Hi,

after some time I am back on this topic. I would have one question concerning the class BinData: since the values I would like to fit are stored in a TGraphErrors, when I use the FillData() command, are they stored and used for the fit both data and errors?

Thank you,

best

Hi,

Yes, when using FillData() passing a TGraphErrors, both the error on x and y will be included by default in the BinData class and sued for fitting. You can control this, by using the possible options in the ROOT::Fit::DataOptions class.
By setting DataOptions::fErrors1 to true, the fit will not use any error (it will be like a TGraph fitting), or by setting DataOptions::fCoordErrors to false, you will not include in the fit the errors on the coordinate (x).

Best Regards

Lorenzo

Thank you Lorenzo for your answer. I would still have one question. Is it possible to have in the output the value of the fitted global chi2? I would like to know this value since the fit procedure seems to work now, and I also get realistic values for my parameter, but when I ask for the probability of the chi2 double prob = result.Prob(); I always get 1.0, which sounds strange…

Thank you again

Hi,

When fitting, you are minimizing the combined chi square. You can evaluate the chi square function that you pass to the fitter at the fitted parameter values.
If you want to have this chi2 value stored in the Fit Result you can simply do, by passing a true flag in Fitter::FitFCN

fitter.FitFCN(6,globalChi2,par0,dataone.Size()+datatwo.Size(),true);
ROOT::Fit::FitResult result = fitter.Result();
double chi2 = result.Chi2();

Lorenzo

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

Hi,

The code looks fine to me, but I would need to have a file with your histograms to be sure is correct.
Are you sure your calculation of the chi2 is correct ? Also, remember, if you have empty bins, they will not by default be considered in the chi2. If you include empty bins with an error of 1 you might get a chi2/ndf that is lower. You can include empty bins by setting ROOT::Fit::DataOptions::fUseEmpty to true.

Best Regards

Lorenzo

Hi Lorenzo,

I think it’s ok in my case to use the default fUseEmpty = false. If you think it’s needed, perhaps I could send you the macro by private message. Just let me know if it would be ok for you.

Best

Hi Lorenzo,

first of all I wish you a good start with the New Year! Back to the problem, I didn’t get any reply to my last post, do you think it would be ok for you to have a quick look at my script?

Thanks a lot!

Best

Hi,

As I said before, the code you have posted looks fine to me. If I have understood well, the only problem you have is a chi2 value that is different from what you expect. This is strange if you are also excluding the empty bins in your manual chi2 computation. Are you sure your manual computation is correct ?
If you want to me to look further in your problem, please post (or send by private message or by email,
moneta at cern . ch ) the macro and the root file with the equired histograms.

Best regards and Happy New Year !

Lorenzo

Hi,

I realise that this is an old thread, but I am approaching this problem in my fit at the moment. I am basing my fitting code on the combinedFit.C tutorial and I have 24 functions to combine. Therefore, I would like to use arrays of ROOT::Math::WrappedMultiTF1, ROOT::Fit::BinData, ROOT::Fit::Chi2Function. I understand that this isn’t possible, but can be done with “a vector of object pointers”. Could you possibly post an example of how to initialise and use such an array please?

Secondly, is it possible to give an array of chi2 values to the global chi2 function? If so, could you possibly post an example of a function like this?

Best regards,
Liam

I am working on a similar problem fitting a variable number of histograms to functions that share some parameters.

One question I have is where is a specification for all the valid arguments for fitter.Config().MinimizerOptions().SetMinimizerType()
I cannot find it.

I am trying to make a vector of ROOT::Fit::Chi2Function one for each histogram, but I run into the problem when I try to push_back an individual ROOT::Fit::Chi2Function into the vector as I loop over the histograms.

The error is

/usr/include/c++/4.7/bits/stl_vector.h:893:4: required from ‘void std::vector<_Tp, _Alloc>::push_back(const value_type&) [with _Tp = ROOT::Fit::Chi2FCNROOT::Math::IBaseFunctionMultiDim; _Alloc = std::allocator<ROOT::Fit::Chi2FCNROOT::Math::IBaseFunctionMultiDim >; std::vector<_Tp, _Alloc>::value_type = ROOT::Fit::Chi2FCNROOT::Math::IBaseFunctionMultiDim]’
/home/wam/src/tone/./mslrfit-v4.cpp:1066:29: required from here
/home/wam/src/root/include/Fit/Chi2FCN.h:64:7: error: non-static reference member ‘const ROOT::Fit::BinData& ROOT::Fit::Chi2FCNROOT::Math::IBaseFunctionMultiDim::fData’, can’t use default assignment operator
/home/wam/src/root/include/Fit/Chi2FCN.h:64:7: error: non-static reference member ‘const IModelFunction& ROOT::Fit::Chi2FCNROOT::Math::IBaseFunctionMultiDim::fFunc’, can’t use default assignment operator

Thanks for any advice!

Andrew

Hi,

Currently you cannot use the assignment operator for the ROOT::Fit::Chi2FCN class, since it contains reference members.
You should make a vector of pointers, a
std::vector< ROOT::Fit::Chi2FCNROOT::Math::IBaseFunctionMultiDim *>

Cheers

Lorenzo

Thanks!

I am still struggling with how to pass the pointers to the individual histogram chi2 functions into the globalChi2 class, so they are available for it to use, but I think I sorted this out.

Andrew

I have now implemented fitting multiple histograms with the same function. In doing this I made a vector<TF1*> for the fit functions. In a loop I set up the function parameters, including fixing some of the parameters, before pushing back the TF1 into the vector. This works fine EXCEPT the fixed parameters are no longer fixed when I later use the vector entries for fitting.

Hi,
I try to run this macro
combinedFit.C (3.5 KB)
with ROOT 5.34 but, it gave some errors:


I don’t know how to solve this problem. Could you help me?

If I try to compile this macro with ACliC even then some problems occur: