// combined fit of two histogram with separate functions but a common parameter // Version2: In this case the FCN is impelmented using the ROOT::Mat::FitMethodFunciton interface // This makes the fit possible also with GSLMultiFit, Fumili and Fumili2 // // N.B. this macro must be compiled with ACliC #include "Fit/Fitter.h" #include "Fit/BinData.h" #include "Fit/Chi2FCN.h" #include "TH1.h" #include "TList.h" #include "Math/WrappedMultiTF1.h" #include "Math/FitMethodFunction.h" #include "HFitInterface.h" #include "TCanvas.h" #include "TStyle.h" // definition of shared parameter // background function int iparB[2] = { 0, // exp amplitude in B histo 2 // exp common parameter }; // signal + background function int iparSB[5] = { 1, // exp amplitude in S+B histo 2, // exp common parameter 3, // gaussian amplitude 4, // gaussian mean 5 // gaussian sigma }; class GlobalChi2 : public ROOT::Math::FitMethodFunction { public: GlobalChi2( int dim, int npoints, ROOT::Math::FitMethodFunction & f1, ROOT::Math::FitMethodFunction & f2) : ROOT::Math::FitMethodFunction(dim,npoints), fChi2_1(&f1), fChi2_2(&f2) {} ROOT::Math::IMultiGenFunction * Clone() const { // copy using default copy-ctor // i.e. function pointer will be copied (and not the functions) return new GlobalChi2(*this); } double DataElement(const double *par, unsigned int ipoint, double *g = 0) const { // implement evaluation of single chi2 element 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] ]; double g1[2]; double g2[5]; double value = 0; if (g != 0) { for (int i = 0; i < 6; ++i) g[i] = 0; } if (ipoint < fChi2_1->NPoints() ) { if (g != 0) { value = fChi2_1->DataElement(p1, ipoint, g1); // update gradient values for (int i= 0; i < 2; ++i) g[iparB[i]] = g1[i]; } else // no need to compute gradient in this case value = fChi2_1->DataElement(p1, ipoint); } else { unsigned int jpoint = ipoint- fChi2_1->NPoints(); assert (jpoint < fChi2_2->NPoints() ); if ( g != 0) { value = fChi2_2->DataElement(p2, jpoint, g2); // update gradient values for (int i= 0; i < 5; ++i) g[iparSB[i]] = g2[i]; } else // no need to compute gradient in this case value = fChi2_2->DataElement(p2, jpoint); } return value; } // needed if want to use Fumili or Fumili2 virtual Type_t Type() const { return ROOT::Math::FitMethodFunction::kLeastSquare; } private: // parameter vector is first background (in common 1 and 2) and then is signal (only in 2) virtual double DoEval (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::FitMethodFunction * fChi2_1; const ROOT::Math::FitMethodFunction * fChi2_2; }; Double_t old_par=286.26; Double_t myexp(Double_t *x, Double_t *par) { if(old_par!=par[0]) { std::cout << "Change " << old_par << " " << par[0] << std::endl; old_par=par[0]; } return par[0]*exp(x[0]*par[1]); } void combinedFit2() { int nbins = 10; TH1D * hB = new TH1D("hB","histo B",nbins,0,100); TH1D * hSB = new TH1D("hSB","histo S+B",nbins, 0,100); // TF1 * fB = new TF1("fB","expo",0,100); TF1 * fB = new TF1("fB",myexp,0,100, 2); fB->SetParameters(286.26,-0.05); // fB->SetParameters(1,-0.05); hB->FillRandom("fB"); TF1 * fS = new TF1("fS","gaus",0,100); fS->SetParameters(1,30,5); hSB->FillRandom("fB",2000); hSB->FillRandom("fS",1000); // perform now global fit TF1 * fSB = new TF1("fSB","expo + gaus(2)",0,100); ROOT::Math::WrappedMultiTF1 wfB(*fB,1); ROOT::Math::WrappedMultiTF1 wfSB(*fSB,1); ROOT::Fit::DataOptions opt; ROOT::Fit::DataRange rangeB; // set the data range rangeB.SetRange(10,90); ROOT::Fit::BinData dataB(opt,rangeB); ROOT::Fit::FillData(dataB, hB); ROOT::Fit::DataRange rangeSB; rangeSB.SetRange(10,50); ROOT::Fit::BinData dataSB(opt,rangeSB); ROOT::Fit::FillData(dataSB, hSB); ROOT::Fit::Chi2Function chi2_B(dataB, wfB); ROOT::Fit::Chi2Function chi2_SB(dataSB, wfSB); ROOT::Fit::Fitter fitter; const int Npar = 6; double par0[Npar] = { 286.26,5,-0.1,100, 30,10}; // create before the parameter settings in order to fix or set range on them fitter.Config().SetParamsSettings(6,par0); // fix 5-th parameter fitter.Config().ParSettings(0).Fix(); // set limits on the third and 4-th parameter fitter.Config().ParSettings(2).SetLimits(-10,-1.E-4); fitter.Config().ParSettings(3).SetLimits(0,10000); fitter.Config().ParSettings(3).SetStepSize(5); fitter.Config().MinimizerOptions().SetPrintLevel(3); fitter.Config().MinimizerOptions().SetMinimizerType("GSLMultiFit"); // fitter.Config().MinimizerOptions().SetMinimizerType("Fumili"); // fit FCN function directly (specify optionally data size and flag to indicate that is a chi2 fit) GlobalChi2 globalChi2(Npar,dataB.Size()+dataSB.Size(),chi2_B, chi2_SB); fitter.FitFCN(globalChi2,0,dataB.Size()+dataSB.Size(),true); ROOT::Fit::FitResult result = fitter.Result(); result.Print(std::cout); TCanvas * c1 = new TCanvas("Simultaneous fit"); c1->Divide(1,2); c1->cd(1); gStyle->SetOptFit(1111); fB->SetFitResult( result, iparB); fB->SetRange(rangeB().first, rangeB().second); hB->GetListOfFunctions()->Add(fB); hB->Draw(); c1->cd(2); fSB->SetFitResult( result, iparSB); fSB->SetRange(rangeSB().first, rangeSB().second); hSB->GetListOfFunctions()->Add(fSB); hSB->Draw(); } int main() { combinedFit2(); }