diff --git a/tutorials/fit/NumericalMinimization.C b/tutorials/fit/NumericalMinimization.C index 252a6464381ee81062ee69f2e5b226ea61e4184f..7ec522367f9c675eb3dafc988f52313741614e3b 100644 --- a/tutorials/fit/NumericalMinimization.C +++ b/tutorials/fit/NumericalMinimization.C @@ -1,10 +1,12 @@ -// tutorial showing how to use the Minimizer class in ROOT -// and all the possible minimizer +// Example on how to use the new Minimizer class in ROOT +// Show usage with all the possible minimizers. // Minimize the Rosenbrock function (a 2D -function) // This example is described also in // http://root.cern.ch/drupal/content/numerical-minimization#multidim_minim // input : minimizer name + algorithm name // randomSeed: = <0 : fixed value: 0 random with seed 0; >0 random with given seed +// +//Author: L. Moneta Dec 2010 #include "Math/Minimizer.h" #include "Math/Factory.h" @@ -22,19 +24,24 @@ double RosenBrock(const double *xx ) return 100*tmp1*tmp1+tmp2*tmp2; } -int NumericalMinimization(const char * minName = "Minuit2",const char *algoName = "" , int randomSeed = -1) +int NumericalMinimization(const char * minName = "Minuit2", + const char *algoName = "" , + int randomSeed = -1) { - // create minimizer giving a name and a name (optionally) for the specific algorithm + // create minimizer giving a name and a name (optionally) for the specific + // algorithm // possible choices are: // minName algoName // Minuit /Minuit2 Migrad, Simplex,Combined,Scan (default is Migrad) // Minuit2 Fumili2 // Fumili - // GSLMultiMin ConjugateFR, ConjugatePR, BFGS, BFGS2, SteepestDescent + // GSLMultiMin ConjugateFR, ConjugatePR, BFGS, + // BFGS2, SteepestDescent // GSLMultiFit // GSLSimAn // Genetic - ROOT::Math::Minimizer* min = ROOT::Math::Factory::CreateMinimizer(minName, algoName); + ROOT::Math::Minimizer* min = + ROOT::Math::Factory::CreateMinimizer(minName, algoName); // set tolerance , etc... min->SetMaxFunctionCalls(1000000); // for Minuit/Minuit2 @@ -70,9 +77,11 @@ int NumericalMinimization(const char * minName = "Minuit2",const char *algoName // expected minimum is 0 if ( min->MinValue() < 1.E-4 && f(xs) < 1.E-4) - std::cout << "Minimizer " << minName << " - " << algoName << " converged to the right minimum" << std::endl; + std::cout << "Minimizer " << minName << " - " << algoName + << " converged to the right minimum" << std::endl; else { - std::cout << "Minimizer " << minName << " - " << algoName << " failed to converge !!!" << std::endl; + std::cout << "Minimizer " << minName << " - " << algoName + << " failed to converge !!!" << std::endl; Error("NumericalMinimization","fail to converge"); } diff --git a/tutorials/fit/combinedFit.C b/tutorials/fit/combinedFit.C new file mode 100644 index 0000000000000000000000000000000000000000..d8570950d0b93b61d0dbb50f5370879631f8a80b --- /dev/null +++ b/tutorials/fit/combinedFit.C @@ -0,0 +1,139 @@ +//+ Combined (simultaneous) fit of two histogram with separate functions +// and some common parameters +// +// See http://root.cern.ch/phpBB3//viewtopic.php?f=3&t=11740#p50908 +// for a modified version working with Fumili or GSLMultiFit +// +// N.B. this macro must be compiled with ACliC +// +//Author: L. Moneta - Dec 2010 + +#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" + + +// 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 +}; + +struct GlobalChi2 { + GlobalChi2( ROOT::Math::IMultiGenFunction & f1, + ROOT::Math::IMultiGenFunction & f2) : + fChi2_1(&f1), fChi2_2(&f2) {} + + // parameter vector is first background (in common 1 and 2) + // and then is signal (only in 2) + 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; +}; + +void combinedFit() { + + TH1D * hB = new TH1D("hB","histo B",100,0,100); + TH1D * hSB = new TH1D("hSB","histo S+B",100, 0,100); + + TF1 * fB = new TF1("fB","expo",0,100); + 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); + + GlobalChi2 globalChi2(chi2_B, chi2_SB); + + ROOT::Fit::Fitter fitter; + + const int Npar = 6; + double par0[Npar] = { 5,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(4).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(0); + fitter.Config().SetMinimizer("Minuit2","Migrad"); + + // fit FCN function directly + // (specify optionally data size and flag to indicate that is a chi2 fit) + fitter.FitFCN(6,globalChi2,0,dataB.Size()+dataSB.Size(),true); + ROOT::Fit::FitResult result = fitter.Result(); + result.Print(std::cout); + + TCanvas * c1 = new TCanvas("Simfit","Simultaneous fit of two histograms", + 10,10,700,700); + c1->Divide(1,2); + c1->cd(1); + gStyle->SetOptFit(1111); + + fB->SetFitResult( result, iparB); + fB->SetRange(rangeB().first, rangeB().second); + fB->SetLineColor(kBlue); + hB->GetListOfFunctions()->Add(fB); + hB->Draw(); + + c1->cd(2); + fSB->SetFitResult( result, iparSB); + fSB->SetRange(rangeSB().first, rangeSB().second); + fSB->SetLineColor(kRed); + hSB->GetListOfFunctions()->Add(fSB); + hSB->Draw(); + + +}