diff --git a/hist/hist/src/TFractionFitter.cxx b/hist/hist/src/TFractionFitter.cxx index 6a395a154deff..0d9f23324f08d 100644 --- a/hist/hist/src/TFractionFitter.cxx +++ b/hist/hist/src/TFractionFitter.cxx @@ -3,8 +3,10 @@ // with additions by Bram Wijngaarden /** \class TFractionFitter -Fits MC fractions to data histogram. A la HMCMLL, see R. Barlow and C. Beeston, -Comp. Phys. Comm. 77 (1993) 219-228, and http://www.hep.man.ac.uk/~roger/hfrac.f +Fits MC fractions to data histogram. À la [HMCMLL](https://cds.cern.ch/record/2296378/files/hbook.pdf), +see R. Barlow and C. Beeston, Comp. Phys. Comm. 77 (1993) 219-228 +\see https://indico.in2p3.fr/event/2635/contributions/25070/ +\note An alternative interface is described in https://root.cern.ch/doc/master/rf709__BarlowBeeston_8C_source.html The virtue of this fit is that it takes into account both data and Monte Carlo statistical uncertainties. The way in which this is done is through a standard diff --git a/hist/hist/test/CMakeLists.txt b/hist/hist/test/CMakeLists.txt index c5c9199b56d8e..25c651500f39d 100644 --- a/hist/hist/test/CMakeLists.txt +++ b/hist/hist/test/CMakeLists.txt @@ -5,6 +5,7 @@ # For the list of contributors see $ROOTSYS/README/CREDITS. ROOT_ADD_GTEST(testTProfile2Poly test_tprofile2poly.cxx LIBRARIES Hist Matrix MathCore RIO) +ROOT_ADD_GTEST(testTFractionFitter test_TFractionFitter.cxx LIBRARIES Hist MathCore) ROOT_ADD_GTEST(testTH2PolyBinError test_TH2Poly_BinError.cxx LIBRARIES Hist Matrix MathCore RIO) ROOT_ADD_GTEST(testTH2PolyAdd test_TH2Poly_Add.cxx LIBRARIES Hist Matrix MathCore RIO) ROOT_ADD_GTEST(testTH2PolyGetNumberOfBins test_TH2Poly_GetNumberOfBins.cxx LIBRARIES Hist Matrix MathCore RIO) diff --git a/hist/hist/test/test_TFractionFitter.cxx b/hist/hist/test/test_TFractionFitter.cxx new file mode 100644 index 0000000000000..ecf28e6896d1b --- /dev/null +++ b/hist/hist/test/test_TFractionFitter.cxx @@ -0,0 +1,143 @@ +#include "gtest/gtest.h" +#include "ROOT/TestSupport.hxx" + +#include "TH1F.h" +#include "TF1.h" +#include "TFractionFitter.h" +#include "TRandom.h" + +// https://its.cern.ch/jira/browse/ROOT-9330 +TEST(TFractionFitter, FitExample) +{ + // pointers to the data + TH1F *data; // data histogram + TH1F *mc0; // first MC histogram + TH1F *mc1; // second MC histogram + TH1F *mc2; // third MC histogram + + // parameters and functions to generate the data + Int_t Ndata = 1000; + Int_t N0 = 1000; + Int_t N1 = 1000; + Int_t N2 = 1000; + + Int_t nBins = 40; + + Double_t trueP0 = .01; + Double_t trueP1 = .3; + // Double_t trueP2 = 1. - trueP0 - trueP1; + + // contribution 0 + TF1 *f0 = new TF1("f0", "[0]*(1-cos(x))/TMath::Pi()", 0., TMath::Pi()); + f0->SetParameter(0, 1.); + f0->SetLineColor(2); + // Double_t int0 = f0->Integral(0., TMath::Pi()); + + // contribution 1 + TF1 *f1 = new TF1("f1", "[0]*(1-cos(x)*cos(x))*2./TMath::Pi()", 0., TMath::Pi()); + f1->SetParameter(0, 1.); + f1->SetLineColor(3); + // Double_t int1 = f1->Integral(0., TMath::Pi()); + + // contribution 2 + TF1 *f2 = new TF1("f2", "[0]*(1+cos(x))/TMath::Pi()", 0., TMath::Pi()); + f2->SetParameter(0, 1.); + f2->SetLineColor(4); + // Double_t int2 = f2->Integral(0., TMath::Pi()); + + // generate data + data = new TH1F("data", "Data angle distribution", nBins, 0, TMath::Pi()); + data->SetXTitle("x"); + data->SetMarkerStyle(20); + data->SetMarkerSize(.7); + data->SetMinimum(0); + TH1F *htruemc0 = new TH1F(*data); + htruemc0->SetLineColor(2); + TH1F *htruemc1 = new TH1F(*data); + htruemc1->SetLineColor(3); + TH1F *htruemc2 = new TH1F(*data); + htruemc2->SetLineColor(4); + Double_t p, x; + for (Int_t i = 0; i < Ndata; i++) { + p = gRandom->Uniform(); + if (p < trueP0) { + x = f0->GetRandom(); + htruemc0->Fill(x); + } + else if (p < trueP0 + trueP1) { + x = f1->GetRandom(); + htruemc1->Fill(x); + } + else { + x = f2->GetRandom(); + htruemc2->Fill(x); + } + data->Fill(x); + } + + // generate MC samples + mc0 = new TH1F("mc0", "MC sample 0 angle distribution", nBins, 0, TMath::Pi()); + mc0->SetXTitle("x"); + mc0->SetLineColor(2); + mc0->SetMarkerColor(2); + mc0->SetMarkerStyle(24); + mc0->SetMarkerSize(.7); + for (Int_t i = 0; i < N0; i++) { + mc0->Fill(f0->GetRandom()); + } + + mc1 = new TH1F("mc1", "MC sample 1 angle distribution", nBins, 0, TMath::Pi()); + mc1->SetXTitle("x"); + mc1->SetLineColor(3); + mc1->SetMarkerColor(3); + mc1->SetMarkerStyle(24); + mc1->SetMarkerSize(.7); + for (Int_t i = 0; i < N1; i++) { + mc1->Fill(f1->GetRandom()); + } + + mc2 = new TH1F("mc2", "MC sample 2 angle distribution", nBins, 0, TMath::Pi()); + mc2->SetXTitle("x"); + mc2->SetLineColor(4); + mc2->SetMarkerColor(4); + mc2->SetMarkerStyle(24); + mc2->SetMarkerSize(.7); + for (Int_t i = 0; i < N2; i++) { + mc2->Fill(f2->GetRandom()); + } + + // FractionFitter + TObjArray *mc = new TObjArray(3); // MC histograms are put in this array + mc->Add(mc0); + mc->Add(mc1); + mc->Add(mc2); + TFractionFitter *fit = new TFractionFitter(data, mc); // initialise + fit->Constrain(0, 0.0, 1.0); // constrain fraction 1 to be between 0 and 1 + fit->Constrain(1, 0.0, 1.0); // constrain fraction 1 to be between 0 and 1 + fit->Constrain(2, 0.0, 1.0); // constrain fraction 1 to be between 0 and 1 + // fit->SetRangeX(1,15); // use only the first 15 bins in the fit + ROOT::TestSupport::CheckDiagsRAII diags; + diags.requiredDiag(kWarning, "Minuit2", "DavidonErrorUpdator", false); + diags.requiredDiag(kWarning, "Minuit2", "VariableMetricBuilder Matrix", false); + diags.requiredDiag(kWarning, "Minuit2", "MnPosDef non-positive diagonal", false); + diags.requiredDiag(kWarning, "Minuit2", "MnPosDef Added", false); + diags.requiredDiag(kWarning, "Minuit2", "VariableMetricBuilder gdel", false); + diags.requiredDiag(kWarning, "Minuit2", "VariableMetricBuilder No", false); + diags.requiredDiag(kWarning, "Minuit2", "VariableMetricBuilder Iterations", false); + Int_t status = fit->Fit(); // perform the fit + EXPECT_EQ(status, 0); + + // Cleanup + delete fit; + delete data; + delete mc; + delete mc0; + delete mc1; + delete mc2; + delete htruemc0; + delete htruemc1; + delete htruemc2; + delete f0; + delete f1; + delete f2; +} diff --git a/tutorials/math/fit/fitFraction.C b/tutorials/math/fit/fitFraction.C new file mode 100644 index 0000000000000..d905fed4d1502 --- /dev/null +++ b/tutorials/math/fit/fitFraction.C @@ -0,0 +1,220 @@ +/// \file +/// \ingroup tutorial_fit +/// \notebook +/// FractionFitter example À la [HMCMLL](https://cds.cern.ch/record/2296378/files/hbook.pdf), +/// see R. Barlow and C. Beeston, Comp. Phys. Comm. 77 (1993) 219-228, +/// \see https://indico.in2p3.fr/event/2635/contributions/25070/ +/// \note An alternative interface is described in https://root.cern.ch/doc/master/rf709__BarlowBeeston_8C_source.html +/// +/// \macro_image +/// \macro_output +/// \macro_code +/// +/// \author + +#include "TH1F.h" +#include "TF1.h" +#include "TFractionFitter.h" +#include "TRandom.h" +#include "TCanvas.h" +#include "TStyle.h" +#include "TLatex.h" +#include + +void fitFraction() +{ + // pointers to the data + TH1F *data; // data histogram + TH1F *mc0; // first MC histogram + TH1F *mc1; // second MC histogram + TH1F *mc2; // third MC histogram + + // parameters and functions to generate the data + Int_t Ndata = 1000; + Int_t N0 = 1000; + Int_t N1 = 1000; + Int_t N2 = 1000; + + Int_t nBins = 40; + + Double_t trueP0 = .01; + Double_t trueP1 = .3; + Double_t trueP2 = 1. - trueP0 - trueP1; + + // contribution 0 + TF1 *f0 = new TF1("f0", "[0]*(1-cos(x))/TMath::Pi()", 0., TMath::Pi()); + f0->SetParameter(0, 1.); + f0->SetLineColor(2); + Double_t int0 = f0->Integral(0., TMath::Pi()); + + // contribution 1 + TF1 *f1 = new TF1("f1", "[0]*(1-cos(x)*cos(x))*2./TMath::Pi()", 0., TMath::Pi()); + f1->SetParameter(0, 1.); + f1->SetLineColor(3); + Double_t int1 = f1->Integral(0., TMath::Pi()); + + // contribution 2 + TF1 *f2 = new TF1("f2", "[0]*(1+cos(x))/TMath::Pi()", 0., TMath::Pi()); + f2->SetParameter(0, 1.); + f2->SetLineColor(4); + Double_t int2 = f2->Integral(0., TMath::Pi()); + + // generate data + data = new TH1F("data", "Data angle distribution", nBins, 0, TMath::Pi()); + data->SetXTitle("x"); + data->SetMarkerStyle(20); + data->SetMarkerSize(.7); + data->SetMinimum(0); + TH1F *htruemc0 = new TH1F(*data); + htruemc0->SetLineColor(2); + TH1F *htruemc1 = new TH1F(*data); + htruemc1->SetLineColor(3); + TH1F *htruemc2 = new TH1F(*data); + htruemc2->SetLineColor(4); + Double_t p, x; + for (Int_t i = 0; i < Ndata; i++) { + p = gRandom->Uniform(); + if (p < trueP0) { + x = f0->GetRandom(); + htruemc0->Fill(x); + } + else if (p < trueP0 + trueP1) { + x = f1->GetRandom(); + htruemc1->Fill(x); + } + else { + x = f2->GetRandom(); + htruemc2->Fill(x); + } + data->Fill(x); + } + + // generate MC samples + mc0 = new TH1F("mc0", "MC sample 0 angle distribution", nBins, 0, TMath::Pi()); + mc0->SetXTitle("x"); + mc0->SetLineColor(2); + mc0->SetMarkerColor(2); + mc0->SetMarkerStyle(24); + mc0->SetMarkerSize(.7); + for (Int_t i = 0; i < N0; i++) { + mc0->Fill(f0->GetRandom()); + } + + mc1 = new TH1F("mc1", "MC sample 1 angle distribution", nBins, 0, TMath::Pi()); + mc1->SetXTitle("x"); + mc1->SetLineColor(3); + mc1->SetMarkerColor(3); + mc1->SetMarkerStyle(24); + mc1->SetMarkerSize(.7); + for (Int_t i = 0; i < N1; i++) { + mc1->Fill(f1->GetRandom()); + } + + mc2 = new TH1F("mc2", "MC sample 2 angle distribution", nBins, 0, TMath::Pi()); + mc2->SetXTitle("x"); + mc2->SetLineColor(4); + mc2->SetMarkerColor(4); + mc2->SetMarkerStyle(24); + mc2->SetMarkerSize(.7); + for (Int_t i = 0; i < N2; i++) { + mc2->Fill(f2->GetRandom()); + } + + // FractionFitter + TObjArray *mc = new TObjArray(3); // MC histograms are put in this array + mc->Add(mc0); + mc->Add(mc1); + mc->Add(mc2); + TFractionFitter *fit = new TFractionFitter(data, mc); // initialise + fit->Constrain(0, 0.0, 1.0); // constrain fraction 1 to be between 0 and 1 + fit->Constrain(1, 0.0, 1.0); // constrain fraction 1 to be between 0 and 1 + fit->Constrain(2, 0.0, 1.0); // constrain fraction 1 to be between 0 and 1 + // fit->SetRangeX(1,15); // use only the first 15 bins in the fit + Int_t status = fit->Fit(); // perform the fit + std::cout << "Status: " << status << std::endl; + + // Display + gStyle->SetOptStat(0); + TCanvas c("c", "FractionFitter example", 700, 700); + c.Divide(2,2); + + c.cd(1); + f0->DrawClone(); + f0->GetHistogram()->SetTitle("Original MC distributions"); + f1->DrawClone("same"); + f2->DrawClone("same"); + + c.cd(2); + data->SetTitle("Data distribution with true contributions"); + data->DrawClone("EP"); + htruemc0->Draw("same"); + htruemc1->Draw("same"); + htruemc2->Draw("same"); + + c.cd(3); + mc0->SetTitle("MC generated samples with fit predictions"); + mc0->Draw("PE"); + mc1->Draw("PEsame"); + mc2->Draw("PEsame"); + if (status == 0) { // check on fit status + auto mcp0 = (TH1F*)fit->GetMCPrediction(0); + mcp0->SetLineColor(2); + mcp0->Draw("same"); + auto mcp1 = (TH1F*)fit->GetMCPrediction(1); + mcp1->SetLineColor(3); + mcp1->Draw("same"); + auto mcp2 = (TH1F*)fit->GetMCPrediction(2); + mcp2->SetLineColor(4); + mcp2->Draw("same"); + } + + c.cd(4); + Double_t p0, p1, p2, errP0, errP1, errP2; + TH1F *mcp0, *mcp1, *mcp2; + TLatex l; + l.SetTextSize(.035); + Char_t texte[200]; + if (status == 0) { // check on fit status + TH1F* result = (TH1F*) fit->GetPlot(); + fit->GetResult( 0, p0, errP0); + printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 0, trueP0, p0, errP0); + fit->GetResult( 1, p1, errP1); + printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 1, trueP1, p1, errP1); + fit->GetResult( 2, p2, errP2); + printf(" Parameter %d: true %.3f, estim. %.3f +/- %.3f\n", 2, trueP2, p2, errP2); + data->SetTitle("Data distribution with fitted contributions"); + data->DrawClone("Ep"); + result->Draw("same"); + f0->SetParameter(0,Ndata*p0/int0*data->GetBinWidth(1)); + f0->SetLineStyle(2); + f0->DrawClone("same"); + f1->SetParameter(0,Ndata*p1/int1*data->GetBinWidth(1)); + f1->SetLineStyle(2); + f1->DrawClone("same"); + f2->SetParameter(0,Ndata*p2/int2*data->GetBinWidth(1)); + f2->SetLineStyle(2); + f2->DrawClone("same"); + sprintf( texte, "%d: true %.2f, estimated %.2f +/- %.2f\n", 0, trueP0, p0, errP0); + l.DrawTextNDC( .45, .30, texte); + sprintf( texte, "%d: true %.2f, estimated %.2f +/- %.2f\n", 1, trueP1, p1, errP1); + l.DrawTextNDC( .45, .25, texte); + sprintf( texte, "%d: true %.2f, estimated %.2f +/- %.2f\n", 2, trueP2, p2, errP2); + l.DrawTextNDC( .45, .20, texte); + } + + auto cnew = c.DrawClone(); + // Cleanup + delete fit; + delete data; + delete mc; + delete mc0; + delete mc1; + delete mc2; + delete htruemc0; + delete htruemc1; + delete htruemc2; + delete f0; + delete f1; + delete f2; + // delete cnew; +}