rf601_intminuit.C

Go to the documentation of this file.

/// \file
/// \ingroup tutorial_roofit
/// \notebook -js
/// Likelihood and minimization: interactive minimization with MINUIT
///
/// \macro_image
/// \macro_code
/// \macro_output
///
/// \date July 2008
/// \author Wouter Verkerke
 
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "RooProdPdf.h"
#include "RooAddPdf.h"
#include "RooMinimizer.h"
#include "RooFitResult.h"
#include "RooPlot.h"
#include "TCanvas.h"
#include "TAxis.h"
#include "TH1.h"
using namespace RooFit;
 
void rf601_intminuit()
{
   // S e t u p   p d f   a n d   l i k e l i h o o d
   // -----------------------------------------------
 
   // Observable
   RooRealVar x("x", "x", -20, 20);
 
   // Model (intentional strong correlations)
   RooRealVar mean("mean", "mean of g1 and g2", 0);
   RooRealVar sigma_g1("sigma_g1", "width of g1", 3);
   RooGaussian g1("g1", "g1", x, mean, sigma_g1);
 
   RooRealVar sigma_g2("sigma_g2", "width of g2", 4, 3.0, 6.0);
   RooGaussian g2("g2", "g2", x, mean, sigma_g2);
 
   RooRealVar frac("frac", "frac", 0.5, 0.0, 1.0);
   RooAddPdf model("model", "model", RooArgList(g1, g2), frac);
 
   // Generate 1000 events
   std::unique_ptr<RooDataSet> data{model.generate(x, 1000)};
 
   // Construct unbinned likelihood of model w.r.t. data
   std::unique_ptr<RooAbsReal> nll{model.createNLL(*data)};
 
   // I n t e r a c t i v e   m i n i m i z a t i o n ,   e r r o r   a n a l y s i s
   // -------------------------------------------------------------------------------
 
   // Create MINUIT interface object
   RooMinimizer m(*nll);
 
   // Activate verbose logging of MINUIT parameter space stepping
   m.setVerbose(true);
 
   // Call MIGRAD to minimize the likelihood
   m.migrad();
 
   // Print values of all parameters, that reflect values (and error estimates)
   // that are back propagated from MINUIT
   std::unique_ptr<RooArgSet>{model.getParameters(x)}->Print("s");
 
   // Disable verbose logging
   m.setVerbose(false);
 
   // Run HESSE to calculate errors from d2L/dp2
   m.hesse();
 
   // Print value (and error) of sigma_g2 parameter, that reflects
   // value and error back propagated from MINUIT
   sigma_g2.Print();
 
   // Run MINOS on sigma_g2 parameter only
   m.minos(sigma_g2);
 
   // Print value (and error) of sigma_g2 parameter, that reflects
   // value and error back propagated from MINUIT
   sigma_g2.Print();
 
   // S a v i n g   r e s u l t s ,   c o n t o u r   p l o t s
   // ---------------------------------------------------------
 
   // Save a snapshot of the fit result. This object contains the initial
   // fit parameters, the final fit parameters, the complete correlation
   // matrix, the EDM, the minimized FCN , the last MINUIT status code and
   // the number of times the RooFit function object has indicated evaluation
   // problems (e.g. zero probabilities during likelihood evaluation)
   std::unique_ptr<RooFitResult> fitResult{m.save()};
 
   // Make contour plot of mx vs sx at 1,2,3 sigma
   RooPlot *frame = m.contour(frac, sigma_g2, 1, 2, 3);
   frame->SetTitle("Minuit contour plot");
 
   // Print the fit result snapshot
   fitResult->Print("v");
 
   // C h a n g e   p a r a m e t e r   v a l u e s ,   f l o a t i n g
   // -----------------------------------------------------------------
 
   // At any moment you can manually change the value of a (constant)
   // parameter
   mean = 0.3;
 
   // Rerun MIGRAD,HESSE
   m.migrad();
   m.hesse();
   frac.Print();
 
   // Now fix sigma_g2
   sigma_g2.setConstant(true);
 
   // Rerun MIGRAD,HESSE
   m.migrad();
   m.hesse();
   frac.Print();
 
   new TCanvas("rf601_intminuit", "rf601_intminuit", 600, 600);
   gPad->SetLeftMargin(0.15);
   frame->GetYaxis()->SetTitleOffset(1.4);
   frame->Draw();
}