Home | News | Documentation | Download

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)?
Thank you in advance.

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);
  fitter.Config().SetMinimizer("Minuit2","Migrad");
 
   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];
instead of

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.