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");
}
}