This simple script plots the sampling distribution of the profile likelihood ratio test statistic based on the input Model File. To do this one needs to specify the value of the parameter of interest that will be used for evaluating the test statistic and the value of the parameters used for generating the toy data. In this case, it uses the upper-limit estimated from the ProfileLikleihoodCalculator, which assumes the asymptotic chi-square distribution for -2 log profile likelihood ratio. Thus, the script is handy for checking to see if the asymptotic approximations are valid. To aid, that comparison, the script overlays a chi-square distribution as well. The most common parameter of interest is a parameter proportional to the signal rate, and often that has a lower-limit of 0, which breaks the standard chi-square distribution. Thus the script allows the parameter to be negative so that the overlay chi-square is the correct asymptotic distribution.
Processing /builddir/build/BUILD/root-6.10.00/tutorials/roostats/StandardTestStatDistributionDemo.C...
[1mRooFit v3.60 -- Developed by Wouter Verkerke and David Kirkby[0m
Copyright (C) 2000-2013 NIKHEF, University of California & Stanford University
All rights reserved, please
read http:
=== Using the following for ModelConfig ===
Observables:
RooArgSet:: = (obs_x_channel1,weightVar,channelCat)
Nuisance
Parameters:
RooArgSet:: = (alpha_syst2,alpha_syst3,gamma_stat_channel1_bin_0,gamma_stat_channel1_bin_1)
Global Observables:
RooArgSet:: = (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.174888
[#1]
INFO:Minization -- p.d.f. provides expected number of events, including extended term in likelihood.
[#1]
INFO:Minization -- Including the following contraint terms in minimization: (alpha_syst2Constraint,alpha_syst3Constraint,gamma_stat_channel1_bin_0_constraint,gamma_stat_channel1_bin_1_constraint)
[#0]
PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoGLobalFit - find MLE
[#0]
PROGRESS:Minization -- ProfileLikelihoodCalcultor::DoMinimizeNLL -
using Minuit / Migrad with
strategy 1
[#1]
INFO:Minization -- The following expressions have been identified
as constant and will be precalculated and cached: (signal_channel1_nominal,background1_channel1_nominal,background2_channel1_nominal)
[#1]
INFO:Minization -- The following expressions will be evaluated in cache-and-track mode: (mc_stat_channel1)
RooFitResult: minimized FCN value: -1044.38, estimated distance to minimum: 2.40386e-08
covariance matrix quality: Full, accurate covariance matrix
Status : MINIMIZE=0
Floating Parameter FinalValue +/- Error
-------------------- --------------------------
SigXsecOverSM 1.1153e+00 +/- 5.86e-01
alpha_syst2 -8.9017e-03 +/- 9.82e-01
alpha_syst3 1.7918e-02 +/- 9.48e-01
gamma_stat_channel1_bin_0 9.9955e-01 +/- 4.93e-02
gamma_stat_channel1_bin_1 1.0036e+00 +/- 8.01e-02
Warning: lower value
for SigXsecOverSM is at limit 0
--------------------------------------
Will
generate sampling distribution at SigXsecOverSM = 2.32744
1) 0x20386010 RooRealVar:: SigXsecOverSM = 2.32744 +/- 0.586079
L(-3 - 3)
"SigXsecOverSM" 2) 0x203821c0 RooRealVar:: alpha_syst2 = 0.710797 +/- 0.98218
L(-5 - 5)
"alpha_syst2" 3) 0x2038a3a0 RooRealVar:: alpha_syst3 = 0.260222 +/- 0.947589
L(-5 - 5)
"alpha_syst3" 4) 0x202df1d0 RooRealVar:: gamma_stat_channel1_bin_0 = 1.03668 +/- 0.0493172
L(0 - 1.25)
"gamma_stat_channel1_bin_0" 5) 0x202f6020 RooRealVar:: gamma_stat_channel1_bin_1 = 1.05156 +/- 0.0800649
L(0 - 1.5)
"gamma_stat_channel1_bin_1"
using namespace RooFit;
using namespace RooStats;
bool useProof = false;
int nworkers = 0;
void StandardTestStatDistributionDemo(const char* infile = "",
const char* workspaceName = "combined",
const char* modelConfigName = "ModelConfig",
const char* dataName = "obsData"){
int nToyMC = 1000;
bool allowNegativeMu=true;
const char* filename = "";
if (!strcmp(infile,"")) {
filename = "results/example_combined_GaussExample_model.root";
if (!fileExist) {
#ifdef _WIN32
cout << "HistFactory file cannot be generated on Windows - exit" << endl;
return;
#endif
cout <<"will run standard hist2workspace example"<<endl;
gROOT->ProcessLine(
".! prepareHistFactory .");
gROOT->ProcessLine(
".! hist2workspace config/example.xml");
cout <<"\n\n---------------------"<<endl;
cout <<"Done creating example input"<<endl;
cout <<"---------------------\n\n"<<endl;
}
}
else
filename = infile;
if(!file ){
cout <<"StandardRooStatsDemoMacro: Input file " << filename << " is not found" << endl;
return;
}
if(!w){
cout <<"workspace not found" << endl;
return;
}
ModelConfig* mc = (ModelConfig*) w->
obj(modelConfigName);
if(!data || !mc){
cout << "data or ModelConfig was not found" <<endl;
return;
}
double plcUpperLimit = interval->
UpperLimit(*firstPOI);
delete interval;
cout << "\n\n--------------------------------------"<<endl;
cout <<
"Will generate sampling distribution at " << firstPOI->
GetName() <<
" = " << plcUpperLimit <<endl;
if(nPOI>1){
cout <<"not sure what to do with other parameters of interest, but here are their values"<<endl;
}
if(allowNegativeMu)
cout << "tell it to use 1 event"<<endl;
}
firstPOI->
setVal(plcUpperLimit);
if (useProof) {
}
firstPOI->
setVal(plcUpperLimit);
allParameters.
Print(
"v");
TF1* f =
new TF1(
"f",
Form(
"2*ROOT::Math::chisquared_pdf(2*x,%d,0)",nPOI),min,max);
c1->
SaveAs(
"standard_test_stat_distribution.pdf");
}
ROOT 6.10/00 - Reference Guide Generated on Fri Jun 16 2017 12:57:58 using Doxygen 1.8.5.