# GlobalChi2 function for many histograms in ROOT::Fit::Fitter::FitFCN

Hi,

How do I generalize the GlobalChi2 struct from the combinedFit example to fit a flexible amount of histograms??
This struct in then used in the function ROOT::Fit::Fitter::FitFCN.

I want to convert these lines in the code in the example:

`````` struct GlobalChi2 {
GlobalChi2(  ROOT::Math::IMultiGenFunction & f1,
ROOT::Math::IMultiGenFunction & f2) :
fChi2_1(&f1), fChi2_2(&f2) {}
double operator() (const double *par) const {
double p1[2];
for (int i = 0; i < 2; ++i) p1[i] = par[iparB[i] ];

double p2[5];
for (int i = 0; i < 5; ++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;
};
``````

Into something like:

``````struct GlobalChi2 {
GlobalChi2( std::vector<ROOT::Fit::Chi2Function> &f)
{
for(int k=0; k<N; k++) fChi2.push_back(&f[k]);
}
double operator() (const double *par) const
{
double p[N][10];
for(int k=0; k<N; k++) for (int i = 0; i < 10; ++i) p[k][i] = par[ipar[k][i] ];

return for(int k=0; k<N; k++) (*fChi2[k])(p[k]); //this line is wrong
}
std::vector< ROOT::Math::IMultiGenFunction *> fChi2;
};
``````

Thank you in advance and best regards,
J

Hi Jasot. I am inviting @StephanH to help with this.

Cheers,
J.

Thank you @jalopezg and @StephanH, for looking at this. I guess it is more a c++ question than a ROOT question, but I’m sure this is a common use case…

Hi,

What you are doing seems correct. You need only to sum the chi2 contributions when implementing `GlobalChi2::operator()`, i.e. doing something like:

``````double operator() (const double *par) const
{
// use std::vector in case N is not known at compiled time
std::vector<std::vector<double>> p(N);
for(int k=0; k<N; k++) {
p[k].resize(10);
for (int i = 0; i < 10; ++i)
p[k][i] = par[ipar[k][i] ];
}
double sum = 0;
for(int k=0; k<N; k++)
sum += (*fChi2[k])(p[k].data());
return sum;
}
``````

Thanks, that works. I don’t know why I thought it was more complicated.
I just had to change this:
`sum += (*fChi2[k])(p[k].data());`
for:
`sum += (*fChi2[k])(p[k]);`

Thank you!

Dear Jasot,
I am also interested in this generalization.
Could you please provide a working example for two or
more histograms (e.g. the modified version of combinedFit.C)?

Hi @smgrig,
Sorry for the late response, I didn’t see your message. I copy a light version of my macro. It won’t run directly, because you will need to plug the histograms and re-define the fitting function accordingly.
It won’t run in the last version of root, because there is some issues with `std::vector<ROOT::Fit::DataRange>` now. I run it in 6.04 without problems.
Let me know if it helps you.
Cheers,
J.

``````
#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"

const int N=2; //number of histograms to fit
const double Par1[N]={6.0,6.0};
const double Par2[N]={2.5,4.5};
float MaxRange[N]={static_cast<float>(4.5e-6),static_cast<float>(4.5e-6)};
const int r[N]={3542,3639};
const string filename="Histogram_%i.root";

const int NCommonPar=10;
const int NNotCommonPar=5;
const int NPar=15;
const int NTotPar=NCommonPar+NNotCommonPar*N;

struct FuncStruct {
double myPar1;
double myPar2;
int k;
double ExpoConvSimple(double *x,double t0, double sigma, double tau)
{
double t = x[0]; // Variable
double c = 4.;
c/=(2.*tau);
double t1 = (sigma *sigma) / (2 * tau * tau);
double t2 = t - t0;
t2 /= tau;
double t3 = (sigma * sigma) - (tau * (t-t0));
t3 /= (TMath::Sqrt(2) * sigma * tau);
return c * TMath::Exp( t1 - t2 ) * ( 1 - TMath::Erf(t3) );
}
double myFunction(double *x, double *par)
{
double A=par[0];//common
double B=par[1];//common
double C=par[2];//common
double D=par[3];//common
double E=par[4];//common
double F=par[5];//common
double t0=par[6];//common
double sigma=par[7];//common
double T1=par[8];//common
double T2=par[9];//common
double P1=par[10];
double P2=par[11];
double P3=par[12];
double P4=par[13];
double Baseline=par[14];
double T3=1.0/(A*myPar1 + B*myPar2 + C);
double T4=1.0/(D*myPar1 + E*myPar2 + F);
return Baseline + p1*ExpoConvSimple(x, t0, sigma, T1) + P2*ExpoConvSimple(x, t0, sigma, T2) + P3*ExpoConvSimple(x, t0, sigma, T3) + P4*ExpoConvSimple(x, t0, sigma, T4);
}
};

int ipar[N][NPar];

struct GlobalChi2 {
GlobalChi2( std::vector<ROOT::Fit::Chi2Function> &f)
{
for(int k=0; k<f.size(); k++) fChi2.push_back(&f[k]);
Fillipar();
}

void Fillipar()
{
for(int k=0; k<N; k++)for(int i=0; i<NCommonPar; i++){ipar[k][i] = i; }
for(int k=0; k<N; k++) for(int i=NCommonPar; i<NPar; i++){ipar[k][i] = i+k*5;}
}
double operator() (const double *par) const
{
double p[N][15];
for(int k=0; k<N; k++) for (int i = 0; i < NPar; ++i) p[k][i] = par[ipar[k][i] ];
double sum=0;
for(int k=0; k<fChi2.size(); k++)
sum += (*fChi2[k])(p[k%N]);
return sum;
}
std::vector< ROOT::Math::IMultiGenFunction *> fChi2;
};

void GlobalFit_Example() {

TH1D * h[N];
for(int i=0; i<N; i++) h[i]=(TH1D*)GetFromFile(Form(filename.c_str(),r[i]));

TF1 *f[N];
FuncStruct F1[N];
for(int k=0;k<N;k++) F1[k]=FuncStruct({Par1[k],Par2[k],k});

for(int k=0;k<N;k++) f[k]= new TF1(Form("f%i",k),&F1[k],&FuncStruct::myFunction,MinRange,MaxRange[k],NPar);

cout << "Wrapping TF1: " << endl;
std::vector<ROOT::Math::WrappedMultiTF1> wf;
for(int k=0;k<N;k++)wf.push_back(ROOT::Math::WrappedMultiTF1(*f[k],1));

ROOT::Fit::DataOptions opt;
const int NRanges=3;
cout << "Setting DataRanges: " << endl;
std::vector<ROOT::Fit::DataRange> range1;
std::vector<ROOT::Fit::DataRange> range2;
std::vector<ROOT::Fit::DataRange> range3;
for(int k=0;k<N;k++) range1.push_back(ROOT::Fit::DataRange());
for(int k=0;k<N;k++) range2.push_back(ROOT::Fit::DataRange());
for(int k=0;k<N;k++) range3.push_back(ROOT::Fit::DataRange());

double myExcludeMin=160e-9;
double myExcludeMax=300e-9;
double myExcludeMin2=400e-9;
double myExcludeMax2=550e-9;
for(int k=0;k<N;k++)range1[k].SetRange(MinRange,myExcludeMin);
for(int k=0;k<N;k++)range2[k].SetRange(myExcludeMax,myExcludeMin2);
for(int k=0;k<N;k++)range3[k].SetRange(myExcludeMax2,MaxRange[k]);

cout << "Setting BinData: " << endl;
std::vector<ROOT::Fit::BinData> data;cout << 1 << endl;

for(int k=0;k<N;k++){ data.push_back(ROOT::Fit::BinData(opt,range1[k])); cout << k << endl; }
for(int k=0;k<N;k++){ data.push_back(ROOT::Fit::BinData(opt,range2[k])); cout << k << endl; }
for(int k=0;k<N;k++){ data.push_back(ROOT::Fit::BinData(opt,range3[k])); cout << k << endl; }
for(int k=0;k<NRanges*N;k++) {ROOT::Fit::FillData(data[k], h[k%N]);cout << k << " " << k%N << endl;}

cout << "Setting Chi2Function: " << endl;
std::vector<ROOT::Fit::Chi2Function> chi2;
for(int k=0;k<NRanges*N;k++) {chi2.push_back(ROOT::Fit::Chi2Function(data[k], wf[k%N])); cout << k << " " << k%N << endl;}

GlobalChi2 globalChi2(chi2);

ROOT::Fit::Fitter fitter;

std::vector<double> par0 = SetInitialParameters(h);
for(int k=0;k<N;k++) for (int i=0; i<NPar; i++) {f[k]->SetParameter(i,par0[ipar[k][i]]); cout << k << " " << i << " " << ipar[k][i] << " " << par0[ipar[k][i]] << endl;} lets_pause();

fitter.Config().SetParamsSettings(NTotPar,&(par0[0]));
std::vector<std::pair<double,double>> parlimits; parlimits.resize(NTotPar);

for (int i=0; i<NCommonPar; i++) parlimits[i]=std::pair<double,double>({0.3*par0[i],1.9*par0[i]});
for (int i=NCommonPar; i<NTotPar; i++) parlimits[i]=std::pair<double,double>({par0[i]-TMath::Abs(0.9*par0[i]),par0[i]+TMath::Abs(1.2*par0[i])});
for(int k=0;k<N;k++){parlimits[10+5*k]=std::pair<double,double>({par0[10+5*k]-TMath::Abs(1.0*par0[10+5*k]),par0[10+5*k]+TMath::Abs(0.9*par0[10+5*k])}); cout << "limit " << par0[10+5*k] << endl;}
for(int k=0;k<N;k++){parlimits[11+5*k]=std::pair<double,double>({par0[11+5*k]-TMath::Abs(1.0*par0[11+5*k]),par0[11+5*k]+TMath::Abs(3*par0[11+5*k])}); cout << "limit " << par0[11+5*k] << endl;}
for(int k=0;k<N;k++){ fitter.Config().ParSettings(10+5*k).SetLimits(0.0, parlimits[10+5*k].second);}
for(int k=0;k<N;k++){ fitter.Config().ParSettings(11+5*k).SetLimits(0.0, parlimits[11+5*k].second);}
for (int i=0; i<NTotPar; i++) fitter.Config().ParSettings(i).SetStepSize(TMath::Abs(0.01*par0[i]));

fitter.Config().ParSettings(2).Fix();
fitter.Config().ParSettings(8).Fix();

fitter.Config().MinimizerOptions().SetPrintLevel(3);

int datasize=0;
for(int k=0;k<NRanges*N;k++) datasize+=data[k].Size();
fitter.FitFCN(NTotPar,globalChi2,0,datasize,true);
ROOT::Fit::FitResult result = fitter.Result();
result.Print(std::cout);

for (int i=0; i<NTotPar; i++)
{
if(TMath::Abs(1-result.Parameters()[i]/parlimits[i].first)<0.01) cout << "WARNING, PARAMETER " << i << " is on the low imit" << endl;
if(TMath::Abs(1-result.Parameters()[i]/parlimits[i].second)<0.01) cout << "WARNING, PARAMETER " << i << " is on the up limit" << endl;
}

TCanvas * c2 = new TCanvas("c2","Simultaneous fit of two histograms",
10,10,700,700);
c2->DivideSquare(N);
for(int k=0;k<N;k++)
{
h[k]->SetTitle(Form("%.1f par1, %.1f par2", Par1[k],Par2[k]));
f[k]->SetFitResult( result, ipar[k]);
c2->cd(k+1);
h[k]->Draw("HIST"); f[k]->Draw("SAME");
}

}
``````

Hi,
thanks a lot for the help.
I think your macro will run in any root version by using
ROOT::Fit::DataRange range[N];
ROOT::Fit::BinData data[3*N];

``````std::vector<ROOT::Fit::DataRange> range;
std::vector<ROOT::Fit::BinData> data;
``````

Hi,
Which error are you having when using `std::vector<ROOT::Fit::DataRange> ` with the latest ROOT version ? It should work both case, if using an std::vector or a C array as shown in the post above.

Lorenzo

Hi,
I get this:

``````1
root.exe: /builddir/build/BUILD/root-6.22.02/math/mathcore/src/BinData.cxx:235: ROOT::Fit::BinData& ROOT::Fit::BinData::operator=(const ROOT::Fit::BinData&): Assertion `Opt().fIntegral' failed.
``````

When I do:

``````   std::vector<ROOT::Fit::BinData> data;cout << 1 << endl;

for(int k=0;k<N;k++){ data.push_back(ROOT::Fit::BinData(opt,range1[k])); cout << k << endl; }
for(int k=0;k<N;k++){ data.push_back(ROOT::Fit::BinData(opt,range2[k])); cout << k << endl; }
for(int k=0;k<N;k++){ data.push_back(ROOT::Fit::BinData(opt,range3[k])); cout << k << endl; }
``````

I think the same error was reported some time ago and fixed, but it seems that its there again…
Cheers,

Hi,
attached is a modified version of the combinedFit.C tutorial which works with C array but fails with std::vector on lxplus.combinedFitMod3.C (4.2 KB)

``````Processing combinedFitMod3.C+...
Info in <TUnixSystem::ACLiC>: creating shared library /afs/cern.ch/user/s/smbat/newFits/./combinedFitMod3_C.so

*** Break *** segmentation violation

===========================================================
There was a crash.
This is the entire stack trace of all threads:
===========================================================
#0  0x00007f3e7c85246c in waitpid () from /lib64/libc.so.6
#1  0x00007f3e7c7cff62 in do_system () from /lib64/libc.so.6
#2  0x00007f3e7d84ab24 in TUnixSystem::StackTrace() () from /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.22.02/x86_64-centos7-gcc48-opt/lib/libCore.so.6.22
#3  0x00007f3e7d84c7ba in TUnixSystem::DispatchSignals(ESignals) () from /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.22.02/x86_64-centos7-gcc48-opt/lib/libCore.so.6.22
#4  <signal handler called>
#5  0x00007f3e6c8ee5aa in ROOT::Fit::FitData::operator=(ROOT::Fit::FitData const&) () from /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.22.02/x86_64-centos7-gcc48-opt/lib/libMathCore.so.6.22.02
#6  0x00007f3e6c8ddba9 in ROOT::Fit::BinData::operator=(ROOT::Fit::BinData const&) () from /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.22.02/x86_64-centos7-gcc48-opt/lib/libMathCore.so.6.22.02
#7  0x00007f3e6a3ca0b2 in combinedFitMod3() () from /afs/cern.ch/user/s/smbat/newFits/combinedFitMod3_C.so
#8  0x00007f3e7e063066 in ?? ()
#9  0x0000000101923008 in ?? ()
#10 0x000000000072ae30 in ?? ()
#11 0x00007f3e785ef980 in ?? () from /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.22.02/x86_64-centos7-gcc48-opt/lib/libCling.so
#12 0x00007ffcdda5c920 in ?? ()
#13 0x0000000000000000 in ?? ()
===========================================================

The lines below might hint at the cause of the crash.
You may get help by asking at the ROOT forum http://root.cern.ch/forum
Only if you are really convinced it is a bug in ROOT then please submit a
report at http://root.cern.ch/bugs Please post the ENTIRE stack trace
from above as an attachment in addition to anything else
that might help us fixing this issue.
===========================================================
#5  0x00007f3e6c8ee5aa in ROOT::Fit::FitData::operator=(ROOT::Fit::FitData const&) () from /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.22.02/x86_64-centos7-gcc48-opt/lib/libMathCore.so.6.22.02
#6  0x00007f3e6c8ddba9 in ROOT::Fit::BinData::operator=(ROOT::Fit::BinData const&) () from /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.22.02/x86_64-centos7-gcc48-opt/lib/libMathCore.so.6.22.02
#7  0x00007f3e6a3ca0b2 in combinedFitMod3() () from /afs/cern.ch/user/s/smbat/newFits/combinedFitMod3_C.so
#8  0x00007f3e7e063066 in ?? ()
#9  0x0000000101923008 in ?? ()
#10 0x000000000072ae30 in ?? ()
#11 0x00007f3e785ef980 in ?? () from /cvmfs/sft.cern.ch/lcg/app/releases/ROOT/6.22.02/x86_64-centos7-gcc48-opt/lib/libCling.so
#12 0x00007ffcdda5c920 in ?? ()
#13 0x0000000000000000 in ?? ()
``````

Thank you for sharing the code ! I will look into this!

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.