hf001_example.C File Reference

Detailed Description

A ROOT script demonstrating an example of writing and fitting a HistFactory model using C++ only.

 
#include <RooStats/HistFactory/MakeModelAndMeasurementsFast.h>
#include <RooStats/HistFactory/Measurement.h>
 
#include <TFile.h>
#include <TROOT.h>
 
void hf001_example()
{
   std::string InputFile = "./data/example.root";
   // in case the file is not found
   bool bfile = gSystem->AccessPathName(InputFile.c_str());
   if (bfile) {
      std::cout << "Input file is not found - run prepareHistFactory script " << std::endl;
      gROOT->ProcessLine(".! prepareHistFactory .");
      bfile = gSystem->AccessPathName(InputFile.c_str());
      if (bfile) {
         std::cout << "Still no " << InputFile << ", giving up.\n";
         exit(1);
      }
   }
 
   // Create the measurement
   RooStats::HistFactory::Measurement meas("meas", "meas");
 
   meas.SetOutputFilePrefix("./results/example_UsingC");
   meas.SetPOI("SigXsecOverSM");
   meas.AddConstantParam("alpha_syst1");
   meas.AddConstantParam("Lumi");
 
   meas.SetLumi(1.0);
   meas.SetLumiRelErr(0.10);
   meas.SetBinHigh(2);
 
   // Create a channel
 
   RooStats::HistFactory::Channel chan("channel1");
   chan.SetData("data", InputFile);
   chan.SetStatErrorConfig(0.05, "Poisson");
 
   // Now, create some samples
 
   // Create the signal sample
   RooStats::HistFactory::Sample signal("signal", "signal", InputFile);
   signal.AddOverallSys("syst1", 0.95, 1.05);
   signal.AddNormFactor("SigXsecOverSM", 1, 0, 3);
   chan.AddSample(signal);
 
   // Background 1
   RooStats::HistFactory::Sample background1("background1", "background1", InputFile);
   background1.ActivateStatError("background1_statUncert", InputFile);
   background1.AddOverallSys("syst2", 0.95, 1.05);
   chan.AddSample(background1);
 
   // Background 1
   RooStats::HistFactory::Sample background2("background2", "background2", InputFile);
   background2.ActivateStatError();
   background2.AddOverallSys("syst3", 0.95, 1.05);
   chan.AddSample(background2);
 
   // Done with this channel
   // Add it to the measurement:
   meas.AddChannel(chan);
 
   // Collect the histograms from their files,
   // print some output,
   meas.CollectHistograms();
   meas.PrintTree();
 
   // One can print XML code to an output directory:
   // meas.PrintXML( "xmlFromCCode", meas.GetOutputFilePrefix() );
 
   // Now, do the measurement
   std::unique_ptr<RooWorkspace> ws{MakeModelAndMeasurementFast(meas)};
 
   RooStats::ModelConfig *modelConfig = static_cast<RooStats::ModelConfig *>(ws->obj("ModelConfig"));
 
   // Get probability density function and parameters list from model
   RooAbsPdf *pdf = modelConfig->GetPdf();
   RooArgSet globalObservables{*modelConfig->GetGlobalObservables()};
 
   // Perform the fit using Minos to get the correct asymmetric uncertainties
   using namespace RooFit;
   std::unique_ptr<RooFitResult> result{
      pdf->fitTo(*ws->data("obsData"), Save(), PrintLevel(-1), GlobalObservables(globalObservables), Minos(true))};
 
   // Getting list of Parameters of Interest and getting first from them
   RooRealVar *poi = static_cast<RooRealVar *>(modelConfig->GetParametersOfInterest()->first());
 
   std::unique_ptr<RooAbsReal> nll{pdf->createNLL(*ws->data("obsData"))};
   std::unique_ptr<RooAbsReal> profile{nll->createProfile(*poi)};
 
   // frame for future plot
   RooPlot *frame = poi->frame();
 
   frame->SetTitle("");
   frame->GetYaxis()->SetTitle("-log likelihood");
   frame->GetXaxis()->SetTitle(poi->GetTitle());
 
   TCanvas *profileLikelihoodCanvas = new TCanvas{"combined", "", 800, 600};
 
   double xmin = poi->getMin();
   double xmax = poi->getMax();
   TLine *line = new TLine(xmin, .5, xmax, .5);
   line->SetLineColor(kGreen);
   TLine *line90 = new TLine(xmin, 2.71 / 2., xmax, 2.71 / 2.);
   line90->SetLineColor(kGreen);
   TLine *line95 = new TLine(xmin, 3.84 / 2., xmax, 3.84 / 2.);
   line95->SetLineColor(kGreen);
   frame->addObject(line);
   frame->addObject(line90);
   frame->addObject(line95);
 
   nll->plotOn(frame, ShiftToZero(), LineColor(kRed), LineStyle(kDashed));
   profile->plotOn(frame);
 
   frame->SetMinimum(0);
   frame->SetMaximum(2.);
 
   frame->Draw();
 
   // Print fit results to console in verbose mode to see asymmetric uncertainties
   result->Print("v");
}

 
[#2] PROGRESS:HistFactory -- Getting histogram ./data/example.root:/data
[#2] INFO:HistFactory -- Opened input file: ./data/example.root: 
[#2] PROGRESS:HistFactory -- Getting histogram ./data/example.root:/signal
[#2] PROGRESS:HistFactory -- Getting histogram ./data/example.root:/background1
[#2] PROGRESS:HistFactory -- Getting histogram ./data/example.root:/background1_statUncert
[#2] PROGRESS:HistFactory -- Getting histogram ./data/example.root:/background2
Measurement Name: meas   OutputFilePrefix: ./results/example_UsingC   POI: SigXsecOverSM   Lumi: 1  LumiRelErr: 0.1   BinLow: 0   BinHigh: 2  ExportOnly: 1
Constant Params:  alpha_syst1 Lumi
Channels:
    Channel Name: channel1  InputFile: 
    Data:
       InputFile: ./data/example.root   HistoName: data   HistoPath:     HistoAddress: 0x56478f09c750
    statErrorConfig:
       RelErrorThreshold: 0.05    ConstraintType: Poisson
    Samples: 
       Name: signal      Channel: channel1    NormalizeByTheory: True    StatErrorActivate: False
                InputFile: ./data/example.root   HistName: signal  HistoPath:     HistoAddress: 0x564790d57770
       Name: background1       Channel: channel1    NormalizeByTheory: True    StatErrorActivate: False
                InputFile: ./data/example.root   HistName: background1   HistoPath:     HistoAddress: 0x56478d4fa610
          StatError Activate: 1   InputFile: ./data/example.root   HistName: background1_statUncert    HistoPath:     HistoAddress: 0x56478c8bfdc0
       Name: background2       Channel: channel1    NormalizeByTheory: True    StatErrorActivate: False
                InputFile: ./data/example.root   HistName: background2   HistoPath:     HistoAddress: 0x564790ee17a0
          StatError Activate: 1   InputFile: ./data/example.root   HistName:   HistoPath:     HistoAddress: 0
    End of Channel channel1
[#2] INFO:HistFactory -- End Measurement: meas
[#2] INFO:HistFactory -- Making Model and Measurements (Fast) for measurement: meas
[#2] INFO:HistFactory -- using lumi = 1 and lumiError = 0.1 including bins between 0 and 2
[#2] INFO:HistFactory -- fixing the following parameters:
   alpha_syst1
   Lumi
[#2] INFO:HistFactory -- Creating the output file: ./results/example_UsingC_meas.root
[#2] INFO:HistFactory -- Creating the HistoToWorkspaceFactoryFast factory
[#2] INFO:HistFactory -- Setting preprocess functions
[#2] PROGRESS:HistFactory -- Starting to process channel: channel1
[#2] PROGRESS:HistFactory -- 
-----------------------------------------
   Starting to process 'channel1' channel with 1 observables
-----------------------------------------
 
[#2] INFO:HistFactory -- making normFactor: SigXsecOverSM
[#2] INFO:HistFactory -- processing hist signal
[#2] INFO:HistFactory -- signal_channel1 has no variation histograms 
[#2] INFO:HistFactory -- processing hist background1
[#2] INFO:HistFactory -- background1_channel1 has no variation histograms 
[#2] INFO:HistFactory -- Sample: background1 to be included in Stat Error for channel channel1
[#2] INFO:HistFactory -- Using external histogram for Stat Errors for   Channel: channel1 Sample: background1  Error Histogram: background1_statUncert
[#2] INFO:HistFactory -- processing hist background2
[#2] INFO:HistFactory -- background2_channel1 has no variation histograms 
[#2] INFO:HistFactory -- Sample: background2 to be included in Stat Error for channel channel1
[#2] INFO:HistFactory -- Making Statistical Uncertainty Hist for  Channel: channel1 Sample: background2
[#2] INFO:HistFactory -- Making Total Uncertainty for bin 1 Error = 5 CentralVal = 100 RelativeError = 0.05
[#2] INFO:HistFactory -- Making Total Uncertainty for bin 2 Error = 10 CentralVal = 100 RelativeError = 0.1
[#2] INFO:HistFactory -- About to create Constraint Terms from: mc_stat_channel1 params: (gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1)
[#2] INFO:HistFactory -- Using Poisson StatErrors in channel: channel1
[#2] INFO:HistFactory -- Creating constraint for: gamma_stat_channel1_bin_0. Type of constraint: 1
[#2] INFO:HistFactory -- Creating constraint for: gamma_stat_channel1_bin_1. Type of constraint: 1
[#2] PROGRESS:HistFactory -- 
-----------------------------------------
   import model into workspace
-----------------------------------------
 
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
RooDataSet::AsimovData[obs_x_channel1,weight:binWeightAsimov] = 2 entries (230 weighted)
 
RooWorkspace(channel1) channel1 workspace contents
 
variables
---------
(Lumi,SigXsecOverSM,alpha_syst1,alpha_syst2,alpha_syst3,gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1,nominalLumi,obs_x_channel1)
 
p.d.f.s
-------
RooGaussian::alpha_syst1Constraint[ x=alpha_syst1 mean=nom_alpha_syst1 sigma=1 ] = 1
RooGaussian::alpha_syst2Constraint[ x=alpha_syst2 mean=nom_alpha_syst2 sigma=1 ] = 1
RooGaussian::alpha_syst3Constraint[ x=alpha_syst3 mean=nom_alpha_syst3 sigma=1 ] = 1
RooRealSumPdf::channel1_model[ signal_channel1_scaleFactors * signal_channel1_shapes + background1_channel1_scaleFactors * background1_channel1_shapes + background2_channel1_scaleFactors * background2_channel1_shapes ] = 220/230
RooPoisson::gamma_stat_channel1_bin_0_constraint[ x=nom_gamma_stat_channel1_bin_0 mean=gamma_stat_channel1_bin_0_poisMean ] = 0.019943
RooPoisson::gamma_stat_channel1_bin_1_constraint[ x=nom_gamma_stat_channel1_bin_1 mean=gamma_stat_channel1_bin_1_poisMean ] = 0.039861
RooGaussian::lumiConstraint[ x=Lumi mean=nominalLumi sigma=0.1 ] = 1
RooProdPdf::model_channel1[ lumiConstraint * alpha_syst1Constraint * alpha_syst2Constraint * alpha_syst3Constraint * gamma_stat_channel1_bin_0_constraint * gamma_stat_channel1_bin_1_constraint * channel1_model(obs_x_channel1) ] = 0.174888
 
functions
--------
RooHistFunc::background1_channel1_Hist_alphanominal[ depList=(obs_x_channel1) ] = 0
RooStats::HistFactory::FlexibleInterpVar::background1_channel1_epsilon[ paramList=(alpha_syst2) ] = 1
RooProduct::background1_channel1_scaleFactors[ background1_channel1_epsilon * Lumi ] = 1
RooProduct::background1_channel1_shapes[ background1_channel1_Hist_alphanominal * mc_stat_channel1 * channel1_model_binWidth ] = 0
RooHistFunc::background2_channel1_Hist_alphanominal[ depList=(obs_x_channel1) ] = 100
RooStats::HistFactory::FlexibleInterpVar::background2_channel1_epsilon[ paramList=(alpha_syst3) ] = 1
RooProduct::background2_channel1_scaleFactors[ background2_channel1_epsilon * Lumi ] = 1
RooProduct::background2_channel1_shapes[ background2_channel1_Hist_alphanominal * mc_stat_channel1 * channel1_model_binWidth ] = 200
RooBinWidthFunction::channel1_model_binWidth[ HistFuncForBinWidth=signal_channel1_Hist_alphanominal ] = 2
RooProduct::gamma_stat_channel1_bin_0_poisMean[ gamma_stat_channel1_bin_0 * gamma_stat_channel1_bin_0_tau ] = 400
RooProduct::gamma_stat_channel1_bin_1_poisMean[ gamma_stat_channel1_bin_1 * gamma_stat_channel1_bin_1_tau ] = 100
ParamHistFunc::mc_stat_channel1[ ] = 1
RooHistFunc::signal_channel1_Hist_alphanominal[ depList=(obs_x_channel1) ] = 10
RooStats::HistFactory::FlexibleInterpVar::signal_channel1_epsilon[ paramList=(alpha_syst1) ] = 1
RooProduct::signal_channel1_scaleFactors[ signal_channel1_epsilon * SigXsecOverSM * Lumi ] = 1
RooProduct::signal_channel1_shapes[ signal_channel1_Hist_alphanominal * channel1_model_binWidth ] = 20
 
datasets
--------
RooDataSet::asimovData(obs_x_channel1)
RooDataSet::obsData(obs_x_channel1)
 
embedded datasets (in pdfs and functions)
-----------------------------------------
RooDataHist::signal_channel1_Hist_alphanominalDHist(obs_x_channel1)
RooDataHist::background1_channel1_Hist_alphanominalDHist(obs_x_channel1)
RooDataHist::background2_channel1_Hist_alphanominalDHist(obs_x_channel1)
 
named sets
----------
ModelConfig_GlobalObservables:(nominalLumi,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
ModelConfig_Observables:(obs_x_channel1)
constraintTerms:(lumiConstraint,alpha_syst1Constraint,alpha_syst2Constraint,alpha_syst3Constraint,gamma_stat_channel1_bin_0_constraint,gamma_stat_channel1_bin_1_constraint)
globalObservables:(nominalLumi,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
likelihoodTerms:(channel1_model)
observables:(obs_x_channel1)
observablesSet:(obs_x_channel1)
 
generic objects
---------------
RooStats::ModelConfig::ModelConfig
 
[#2] INFO:HistFactory -- Setting Parameter(s) of Interest as: SigXsecOverSM 
 
=== Using the following for ModelConfig ===
Observables:             RooArgSet:: = (obs_x_channel1)
Parameters of Interest:  RooArgSet:: = (SigXsecOverSM)
Nuisance Parameters:     RooArgSet:: = (alpha_syst2,alpha_syst3,gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1)
Global Observables:      RooArgSet:: = (nominalLumi,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
PDF:                     RooProdPdf::model_channel1[ lumiConstraint * alpha_syst1Constraint * alpha_syst2Constraint * alpha_syst3Constraint * gamma_stat_channel1_bin_0_constraint * gamma_stat_channel1_bin_1_constraint * channel1_model(obs_x_channel1) ] = 0.174888
 
[#2] INFO:HistFactory -- Opening File to hold channel: ./results/example_UsingC_channel1_meas_model.root
[#2] INFO:HistFactory -- About to write channel measurement to file
[#2] PROGRESS:HistFactory -- Writing sample: signal
[#2] PROGRESS:HistFactory -- Writing sample: background1
[#2] PROGRESS:HistFactory -- Writing sample: background2
[#2] PROGRESS:HistFactory -- Saved all histograms
[#2] PROGRESS:HistFactory -- Saved Measurement
[#2] PROGRESS:HistFactory -- Successfully wrote channel to file
[#2] INFO:HistFactory -- full list of observables:
(obs_x_channel1)
[#2] PROGRESS:HistFactory -- 
-----------------------------------------
   Entering combination
-----------------------------------------
 
 
RooWorkspace(combined) combined contents
 
variables
---------
(channelCat,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1,nominalLumi,obs_x_channel1)
 
datasets
--------
RooDataSet::obsData(obs_x_channel1,channelCat)
 
named sets
----------
ModelConfig_GlobalObservables:(nominalLumi,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
ModelConfig_Observables:(obs_x_channel1,channelCat)
globalObservables:(nominalLumi,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
observables:(obs_x_channel1,channelCat)
 
[#2] PROGRESS:HistFactory -- 
-----------------------------------------
   Importing combined model
-----------------------------------------
 
[#2] INFO:HistFactory -- setting alpha_syst1 constant
[#2] INFO:HistFactory -- setting Lumi constant
[#2] PROGRESS:HistFactory -- 
-----------------------------------------
   create toy data
-----------------------------------------
 
[#1] INFO:NumericIntegration -- RooRealIntegral::init(channel1_model_Int[obs_x_channel1]) using numeric integrator RooBinIntegrator to calculate Int(obs_x_channel1)
RooDataSet::AsimovData0[obs_x_channel1,channelCat,weight:binWeightAsimov] = 2 entries (230 weighted)
[#2] INFO:HistFactory -- Setting Parameter(s) of Interest as: SigXsecOverSM 
 
=== Using the following for ModelConfig ===
Observables:             RooArgSet:: = (obs_x_channel1,channelCat)
Parameters of Interest:  RooArgSet:: = (SigXsecOverSM)
Nuisance Parameters:     RooArgSet:: = (alpha_syst2,alpha_syst3,gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1)
Global Observables:      RooArgSet:: = (nominalLumi,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
PDF:                     RooSimultaneous::simPdf[ indexCat=channelCat channel1=model_channel1 ] = 0.190787
 
[#2] PROGRESS:HistFactory -- Writing combined workspace to file: ./results/example_UsingC_combined_meas_model.root
[#2] PROGRESS:HistFactory -- Writing combined measurement to file: ./results/example_UsingC_combined_meas_model.root
[#2] PROGRESS:HistFactory -- Writing sample: signal
[#2] PROGRESS:HistFactory -- Writing sample: background1
[#2] PROGRESS:HistFactory -- Writing sample: background2
[#2] PROGRESS:HistFactory -- Saved all histograms
[#2] PROGRESS:HistFactory -- Saved Measurement
[#1] INFO:Minimization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:Minimization --  Including the following constraint terms in minimization: (lumiConstraint,alpha_syst1Constraint,alpha_syst2Constraint,alpha_syst3Constraint,gamma_stat_channel1_bin_0_constraint,gamma_stat_channel1_bin_1_constraint)
[#1] INFO:Minimization -- The following global observables have been defined and their values are taken from the model: (nominalLumi,nom_alpha_syst1,nom_alpha_syst2,nom_alpha_syst3,nom_gamma_stat_channel1_bin_0,nom_gamma_stat_channel1_bin_1)
[#1] INFO:Fitting -- RooAbsPdf::fitTo(simPdf) fixing normalization set for coefficient determination to observables in data
[#1] INFO:Fitting -- using CPU computation library compiled with -mavx2
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_simPdf_obsData) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: activating const optimization
[#1] INFO:Minimization -- RooAbsMinimizerFcn::setOptimizeConst: deactivating const optimization
[#1] INFO:Minimization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1] INFO:Minimization --  Including the following constraint terms in minimization: (lumiConstraint,alpha_syst1Constraint,alpha_syst2Constraint,alpha_syst3Constraint,gamma_stat_channel1_bin_0_constraint,gamma_stat_channel1_bin_1_constraint)
[#1] INFO:Minimization -- The global observables are not defined , normalize constraints with respect to the parameters (Lumi,SigXsecOverSM,alpha_syst1,alpha_syst2,alpha_syst3,gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1)
[#1] INFO:Fitting -- RooAbsPdf::fitTo(simPdf) fixing normalization set for coefficient determination to observables in data
[#1] INFO:NumericIntegration -- RooRealIntegral::init(gamma_stat_channel1_bin_1_constraint_Int[gamma_stat_channel1_bin_1]) using numeric integrator RooIntegrator1D to calculate Int(gamma_stat_channel1_bin_1)
[#1] INFO:NumericIntegration -- RooRealIntegral::init(gamma_stat_channel1_bin_0_constraint_Int[gamma_stat_channel1_bin_0]) using numeric integrator RooIntegrator1D to calculate Int(gamma_stat_channel1_bin_0)
[#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[SigXsecOverSM]) Creating instance of MINUIT
[#1] INFO:Fitting -- RooAddition::defaultErrorLevel(nll_simPdf_obsData) Summation contains a RooNLLVar, using its error level
[#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[SigXsecOverSM]) determining minimum likelihood for current configurations w.r.t all observable
[#1] INFO:Minimization -- RooProfileLL::evaluate(RooEvaluatorWrapper_Profile[SigXsecOverSM]) minimum found at (SigXsecOverSM=1.11284)
..........................................................................................................................................................................................................
  RooFitResult: minimized FCN value: 15.1129, estimated distance to minimum: 0.000182975
                covariance matrix quality: Full, accurate covariance matrix
                Status : MINIMIZE=0 HESSE=0 MINOS=0 
 
    Constant Parameter    Value     
  --------------------  ------------
                  Lumi    1.0000e+00
           alpha_syst1    0.0000e+00
       nom_alpha_syst1    0.0000e+00
       nom_alpha_syst2    0.0000e+00
       nom_alpha_syst3    0.0000e+00
  nom_gamma_stat_channel1_bin_0    4.0000e+02
  nom_gamma_stat_channel1_bin_1    1.0000e+02
           nominalLumi    1.0000e+00
 
    Floating Parameter  InitialValue    FinalValue (+HiError,-LoError)    GblCorr.
  --------------------  ------------  ----------------------------------  --------
         SigXsecOverSM    1.0000e+00    1.1109e+00 (+6.16e-01,-5.93e-01)  <none>
           alpha_syst2    0.0000e+00   -1.0039e-02 (+9.90e-01,-9.88e-01)  <none>
           alpha_syst3    0.0000e+00    1.5339e-02 (+9.58e-01,-9.50e-01)  <none>
  gamma_stat_channel1_bin_0    1.0000e+00    9.9953e-01 (+5.03e-02,-4.86e-02)  <none>
  gamma_stat_channel1_bin_1    1.0000e+00    1.0050e+00 (+8.08e-02,-7.98e-02)  <none>
 

Author

George Lewis

Definition in file hf001_example.C.