| 1 |
#include <iostream>
|
| 2 |
#include <iomanip>
|
| 3 |
#include <sstream>
|
| 4 |
#include <cstdlib>
|
| 5 |
#include <cmath>
|
| 6 |
#include <algorithm>
|
| 7 |
|
| 8 |
|
| 9 |
#ifndef __CINT__
|
| 10 |
#include "RooGlobalFunc.h"
|
| 11 |
#endif
|
| 12 |
|
| 13 |
#include "RooProfileLL.h"
|
| 14 |
#include "RooAbsPdf.h"
|
| 15 |
#include "RooStats/HypoTestResult.h"
|
| 16 |
#include "RooRealVar.h"
|
| 17 |
#include "RooPlot.h"
|
| 18 |
#include "RooDataSet.h"
|
| 19 |
#include "RooTreeDataStore.h"
|
| 20 |
#include "TTree.h"
|
| 21 |
#include "TCanvas.h"
|
| 22 |
#include "TLine.h"
|
| 23 |
#include "TStopwatch.h"
|
| 24 |
|
| 25 |
#include "RooStats/ProfileLikelihoodCalculator.h"
|
| 26 |
#include "RooStats/MCMCCalculator.h"
|
| 27 |
#include "RooStats/UniformProposal.h"
|
| 28 |
#include "RooStats/FeldmanCousins.h"
|
| 29 |
#include "RooStats/NumberCountingPdfFactory.h"
|
| 30 |
#include "RooStats/ConfInterval.h"
|
| 31 |
#include "RooStats/PointSetInterval.h"
|
| 32 |
#include "RooStats/LikelihoodInterval.h"
|
| 33 |
#include "RooStats/LikelihoodIntervalPlot.h"
|
| 34 |
#include "RooStats/RooStatsUtils.h"
|
| 35 |
#include "RooStats/ModelConfig.h"
|
| 36 |
#include "RooStats/MCMCInterval.h"
|
| 37 |
#include "RooStats/MCMCIntervalPlot.h"
|
| 38 |
#include "RooStats/ProposalFunction.h"
|
| 39 |
#include "RooStats/ProposalHelper.h"
|
| 40 |
#include "RooFitResult.h"
|
| 41 |
|
| 42 |
#include <iostream>
|
| 43 |
#include <iomanip>
|
| 44 |
#include <cmath>
|
| 45 |
#include <fstream>
|
| 46 |
#include <time.h>
|
| 47 |
#include "ConfigFile.h"
|
| 48 |
|
| 49 |
// use this order for safety on library loading
|
| 50 |
using namespace RooFit ;
|
| 51 |
using namespace RooStats ;
|
| 52 |
|
| 53 |
|
| 54 |
|
| 55 |
double simpleProfile2(ConfigFile * config, string Type,
|
| 56 |
double dat, double bkg, double bkg_uncertainty, double sig, double
|
| 57 |
sig_uncertainty, double xsec, double ExpNsigLimit ) //absolute uncertainties!
|
| 58 |
{
|
| 59 |
TStopwatch t;
|
| 60 |
t.Start();
|
| 61 |
|
| 62 |
/////////////////////////////////////////
|
| 63 |
// The Model building stage
|
| 64 |
/////////////////////////////////////////
|
| 65 |
char * name = new char[1024];
|
| 66 |
RooWorkspace* wspace = new RooWorkspace();
|
| 67 |
//sprintf(name,"Poisson::countingModel(obs[50,0,200],sum(s[50,0,100]*ratioSigEff[1.,0,2.],b[50,0,100]*ratioBkgEff[1.,0.,2.]))");
|
| 68 |
double d=(dat>sig+bkg+bkg_uncertainty+sig_uncertainty?dat:sig+bkg+bkg_uncertainty+sig_uncertainty);
|
| 69 |
sprintf(name,"Poisson::countingModel(obs[%f,0,%f],sum(s[%f,%f,%f]*ratioSigEff[1.,0,2.],b[%f,0,%f]*ratioBkgEff[1.,0.,2.]))",
|
| 70 |
dat,d*2.0, ExpNsigLimit, 0.3*ExpNsigLimit, ExpNsigLimit*2., bkg, bkg*2.);
|
| 71 |
wspace->factory(name); // counting model
|
| 72 |
//wspace->factory("Gaussian::sigConstraint(1,ratioSigEff,0.257)"); // 5% signal efficiency uncertainty
|
| 73 |
//wspace->factory("Gaussian::bkgConstraint(1,ratioBkgEff,0.136)"); // 10% background efficiency uncertainty
|
| 74 |
sprintf(name,"Gaussian::sigConstraint(1,ratioSigEff,%f)",sig_uncertainty/sig);
|
| 75 |
wspace->factory(name); // 5% signal efficiency uncertainty
|
| 76 |
sprintf(name,"Gaussian::bkgConstraint(1,ratioBkgEff,%f)",bkg_uncertainty/bkg);
|
| 77 |
wspace->factory(name); // 10% background efficiency uncertainty
|
| 78 |
wspace->factory("PROD::modelWithConstraints(countingModel,sigConstraint,bkgConstraint)"); // product of terms
|
| 79 |
wspace->Print();
|
| 80 |
|
| 81 |
RooAbsPdf* modelWithConstraints = wspace->pdf("modelWithConstraints"); // get the model
|
| 82 |
RooRealVar* obs = wspace->var("obs"); // get the observable
|
| 83 |
RooRealVar* s = wspace->var("s"); // get the signal we care about
|
| 84 |
RooRealVar* b = wspace->var("b"); // get the background and set it to a constant. Uncertainty included in ratioBkgEff
|
| 85 |
b->setConstant();
|
| 86 |
RooRealVar* ratioSigEff = wspace->var("ratioSigEff"); // get uncertaint parameter to constrain
|
| 87 |
RooRealVar* ratioBkgEff = wspace->var("ratioBkgEff"); // get uncertaint parameter to constrain
|
| 88 |
RooArgSet constrainedParams(*ratioSigEff, *ratioBkgEff); // need to constrain these in the fit (should change default behavior)
|
| 89 |
|
| 90 |
// Create an example dataset with 160 observed events
|
| 91 |
obs->setVal(dat);
|
| 92 |
b->setVal(bkg);
|
| 93 |
s->setVal(sig);
|
| 94 |
RooDataSet* data = new RooDataSet("exampleData", "exampleData", RooArgSet(*obs));
|
| 95 |
data->add(*obs);
|
| 96 |
|
| 97 |
RooArgSet all(*s, *ratioBkgEff, *ratioSigEff);
|
| 98 |
|
| 99 |
// not necessary
|
| 100 |
modelWithConstraints->fitTo(*data, RooFit::Constrain(RooArgSet(*ratioSigEff, *ratioBkgEff)));
|
| 101 |
|
| 102 |
// Now let's make some confidence intervals for s, our parameter of interest
|
| 103 |
RooArgSet paramOfInterest(*s);
|
| 104 |
|
| 105 |
ModelConfig modelConfig(new RooWorkspace());
|
| 106 |
modelConfig.SetPdf(*modelWithConstraints);
|
| 107 |
modelConfig.SetParametersOfInterest(paramOfInterest);
|
| 108 |
|
| 109 |
|
| 110 |
// First, let's use a Calculator based on the Profile Likelihood Ratio
|
| 111 |
//ProfileLikelihoodCalculator plc(*data, *modelWithConstraints, paramOfInterest);
|
| 112 |
ProfileLikelihoodCalculator plc(*data, modelConfig);
|
| 113 |
plc.SetTestSize(.10);
|
| 114 |
ConfInterval* lrint = plc.GetInterval(); // that was easy.
|
| 115 |
|
| 116 |
// Let's make a plot
|
| 117 |
TCanvas* dataCanvas = new TCanvas("dataCanvas");
|
| 118 |
dataCanvas->Divide(2,2);
|
| 119 |
|
| 120 |
dataCanvas->cd(3);
|
| 121 |
// --- Generate a toyMC sample from composite PDF ---
|
| 122 |
//RooDataSet *pseudodata = modelWithConstraints->generate(*b,2000) ;
|
| 123 |
// --- Perform extended ML fit of composite PDF to toy data ---
|
| 124 |
//modelWithConstraints.fitTo(*pseudodata,Extended());
|
| 125 |
|
| 126 |
// don't skip drawing, doesnt work without!!!
|
| 127 |
RooPlot* mesframe = b->frame();
|
| 128 |
modelWithConstraints->plotOn(mesframe);
|
| 129 |
mesframe->Draw();
|
| 130 |
|
| 131 |
dataCanvas->cd(1);
|
| 132 |
LikelihoodIntervalPlot plotInt((LikelihoodInterval*)lrint);
|
| 133 |
plotInt.SetTitle("Profile Likelihood Ratio and Posterior for S");
|
| 134 |
plotInt.Draw();
|
| 135 |
|
| 136 |
// Get Lower and Upper limits from Profile Calculator
|
| 137 |
double up = ((LikelihoodInterval*) lrint)->UpperLimit(*s);
|
| 138 |
cout << Type << ": d:"<<dat<<", b:"<<bkg<<"+-"<<bkg_uncertainty
|
| 139 |
<<"; s:"<<sig<<"+-"<<sig_uncertainty <<std::endl;
|
| 140 |
cout << "Profile lower limit on s = " << ((LikelihoodInterval*) lrint)->LowerLimit(*s) << endl;
|
| 141 |
cout << "Profile upper limit on s = " << up << endl;
|
| 142 |
config->add("RooSimpleProfile.signal."+Type+"UpperLimit", up);
|
| 143 |
config->add("RooSimpleProfile.xsec."+Type+"UpperLimit", up/sig * xsec);
|
| 144 |
|
| 145 |
/*
|
| 146 |
|
| 147 |
// Second, use a Calculator based on the Feldman Cousins technique
|
| 148 |
FeldmanCousins fc(*data, modelConfig);
|
| 149 |
fc.UseAdaptiveSampling(true);
|
| 150 |
fc.FluctuateNumDataEntries(false); // number counting analysis: dataset always has 1 entry with N events observed
|
| 151 |
fc.SetNBins(50); // number of points to test per parameter
|
| 152 |
fc.SetTestSize(.10); //95% single sided
|
| 153 |
fc.AdditionalNToysFactor(2);
|
| 154 |
// fc.SaveBeltToFile(true); // optional
|
| 155 |
ConfInterval* fcint = NULL;
|
| 156 |
fcint = fc.GetInterval(); // that was easy.
|
| 157 |
|
| 158 |
|
| 159 |
// Get Lower and Upper limits from FeldmanCousins with profile construction
|
| 160 |
if (fcint != NULL) {
|
| 161 |
double fcul = ((PointSetInterval*) fcint)->UpperLimit(*s);
|
| 162 |
double fcll = ((PointSetInterval*) fcint)->LowerLimit(*s);
|
| 163 |
cout << "FC lower limit on s = " << fcll << endl;
|
| 164 |
cout << "FC upper limit on s = " << fcul << endl;
|
| 165 |
config->add("RooFC.signal."+Type+"UpperLimit", fcul);
|
| 166 |
config->add("RooFC.xsec."+Type+"UpperLimit", fcul/sig * xsec);
|
| 167 |
//TLine* fcllLine = new TLine(fcll, 0, fcll, 1);
|
| 168 |
//TLine* fculLine = new TLine(fcul, 0, fcul, 1);
|
| 169 |
//fcllLine->SetLineColor(kRed);
|
| 170 |
//fculLine->SetLineColor(kRed);
|
| 171 |
//fcllLine->Draw("same");
|
| 172 |
//fculLine->Draw("same");
|
| 173 |
dataCanvas->Update();
|
| 174 |
}
|
| 175 |
*/
|
| 176 |
RooFitResult* fit = modelWithConstraints->fitTo(*data, Save(true));
|
| 177 |
|
| 178 |
// Third, use a Calculator based on Markov Chain monte carlo
|
| 179 |
// Before configuring the calculator, let's make a ProposalFunction
|
| 180 |
// that will achieve a high acceptance rate
|
| 181 |
ProposalHelper ph;
|
| 182 |
ph.SetVariables((RooArgSet&)fit->floatParsFinal());
|
| 183 |
ph.SetCovMatrix(fit->covarianceMatrix());
|
| 184 |
ph.SetUpdateProposalParameters(true);
|
| 185 |
ph.SetCacheSize(100);
|
| 186 |
ProposalFunction* pdfProp = ph.GetProposalFunction(); // that was easy
|
| 187 |
|
| 188 |
MCMCCalculator mc(*data, modelConfig);
|
| 189 |
mc.SetNumIters(100000); // steps to propose in the chain
|
| 190 |
mc.SetTestSize(.10); // 95% CL single-sided
|
| 191 |
mc.SetNumBurnInSteps(100); // ignore first N steps in chain as "burn in"
|
| 192 |
mc.SetProposalFunction(*pdfProp);
|
| 193 |
mc.SetLeftSideTailFraction(0.5); // find a "central" interval
|
| 194 |
MCMCInterval* mcInt = (MCMCInterval*)mc.GetInterval(); // that was easy
|
| 195 |
|
| 196 |
|
| 197 |
/*
|
| 198 |
|
| 199 |
// Plot MCMC interval and print some statistics
|
| 200 |
MCMCIntervalPlot mcPlot(*mcInt);
|
| 201 |
mcPlot.SetLineColor(kMagenta);
|
| 202 |
mcPlot.SetLineWidth(2);
|
| 203 |
mcPlot.Draw("same");
|
| 204 |
*/
|
| 205 |
double mcul = mcInt->UpperLimit(*s);
|
| 206 |
double mcll = mcInt->LowerLimit(*s);
|
| 207 |
cout << "MCMC lower limit on s = " << mcll << endl;
|
| 208 |
cout << "MCMC upper limit on s = " << mcul << endl;
|
| 209 |
cout << "MCMC Actual confidence level: "
|
| 210 |
<< mcInt->GetActualConfidenceLevel() << endl;
|
| 211 |
cout << "CLs Exp upper limit on s = "<<ExpNsigLimit << endl;
|
| 212 |
config->add("RooMCMC.signal."+Type+"UpperLimit", mcul);
|
| 213 |
config->add("RooMCMC.xsec."+Type+"UpperLimit", mcul/sig * xsec);
|
| 214 |
/*
|
| 215 |
// 3-d plot of the parameter points
|
| 216 |
dataCanvas->cd(2);
|
| 217 |
// also plot the points in the markov chain
|
| 218 |
TTree& chain = ((RooTreeDataStore*) mcInt->GetChainAsDataSet()->store())->tree();
|
| 219 |
chain.SetMarkerStyle(6);
|
| 220 |
chain.SetMarkerColor(kRed);
|
| 221 |
chain.Draw("s:ratioSigEff:ratioBkgEff","weight_MarkovChain_local_","box"); // 3-d box proporional to posterior
|
| 222 |
|
| 223 |
// the points used in the profile construction
|
| 224 |
TTree& parameterScan = ((RooTreeDataStore*) fc.GetPointsToScan()->store())->tree();
|
| 225 |
parameterScan.SetMarkerStyle(24);
|
| 226 |
parameterScan.Draw("s:ratioSigEff:ratioBkgEff","","same");
|
| 227 |
|
| 228 |
delete wspace;
|
| 229 |
delete lrint;
|
| 230 |
delete mcInt;
|
| 231 |
delete fcint;
|
| 232 |
delete data;
|
| 233 |
*/
|
| 234 |
/// print timing info
|
| 235 |
t.Stop();
|
| 236 |
t.Print();
|
| 237 |
|
| 238 |
return up;
|
| 239 |
}
|
| 240 |
|
| 241 |
void rooStat(int argc, char *argv[])
|
| 242 |
{
|
| 243 |
//myPDF();
|
| 244 |
//rooStatOrig();
|
| 245 |
if (argc<1){
|
| 246 |
std::cerr<<"Error: Expected '"<<argv[0]<<" <limit config>'!"<<std::endl;
|
| 247 |
exit(1);
|
| 248 |
}
|
| 249 |
|
| 250 |
for (int file=1; file<argc; ++file){
|
| 251 |
std::cout << "opening: "<<argv[file]<<std::endl;
|
| 252 |
ConfigFile config(argv[file]);
|
| 253 |
int data = config.read<int>("data");
|
| 254 |
double bkg = config.read<double>("background");
|
| 255 |
double bkgUncert = config.read<double>("background.uncertainty");
|
| 256 |
double sig = config.read<double>("signal.LO");
|
| 257 |
double sigUncert = config.read<double>("signal.LO.uncertainty");
|
| 258 |
double sigNLO = config.read<double>("signal.NLO");
|
| 259 |
double sigUncertNLO = config.read<double>("signal.NLO.uncertainty");
|
| 260 |
double xsec = config.read<double>("Xsection");
|
| 261 |
double kfactor = config.read<double>("signal.kFactor");
|
| 262 |
double ExpNsigLimit = config.read<double>("RooMCMC.signal.LO.ExpUpperLimit",18.0);
|
| 263 |
double ExpNsigLimitNLO = config.read<double>("RooMCMC.signal.NLO.ExpUpperLimit",18.0);
|
| 264 |
double ObsNsigLimit = config.read<double>("RooMCMC.signal.LO.ObsUpperLimit",18.0);
|
| 265 |
double ObsNsigLimitNLO = config.read<double>("RooMCMC.signal.NLO.ObsUpperLimit",18.0);
|
| 266 |
|
| 267 |
|
| 268 |
double sigSigTau = config.read<double>("signal.LO.signalregion.Tau");
|
| 269 |
double sigSignal = config.read<double>("signal.LO.signalregion.IsoMuon");
|
| 270 |
double sigSigTauNLO = config.read<double>("signal.NLO.signalregion.Tau");
|
| 271 |
double sigSignalNLO = config.read<double>("signal.NLO.signalregion.IsoMuon");
|
| 272 |
double bkgControl = config.read<double>("background.controlregion.IsoMuon");
|
| 273 |
double bkgSignal = config.read<double>("background.signalregion.IsoMuon");
|
| 274 |
|
| 275 |
double b_reduc = (bkg-sigSignal-sigSigTau<0?0:bkg-sigSignal-sigSigTau);
|
| 276 |
|
| 277 |
simpleProfile2(&config, "LO.Obs", data, b_reduc, bkgUncert, sig, sigUncert, xsec, ObsNsigLimit);
|
| 278 |
simpleProfile2(&config, "LO.Exp", b_reduc, b_reduc, bkgUncert, sig, sigUncert, xsec, ExpNsigLimit);
|
| 279 |
|
| 280 |
b_reduc = (bkg-sigSignalNLO-sigSigTauNLO<0?0:bkg-sigSignalNLO-sigSigTauNLO);
|
| 281 |
simpleProfile2(&config, "NLO.Obs", data, b_reduc, bkgUncert, sigNLO, sigUncertNLO, xsec*kfactor, ObsNsigLimitNLO);
|
| 282 |
simpleProfile2(&config, "NLO.Exp", b_reduc, b_reduc, bkgUncert, sigNLO, sigUncertNLO, xsec*kfactor, ExpNsigLimitNLO);
|
| 283 |
|
| 284 |
//write stuff:
|
| 285 |
time_t rawtime;
|
| 286 |
struct tm * timeinfo;
|
| 287 |
time ( &rawtime );
|
| 288 |
timeinfo = localtime ( &rawtime );
|
| 289 |
ofstream ofile;
|
| 290 |
ofile.open (argv[file]);
|
| 291 |
if (ofile.good())
|
| 292 |
std::cout << "writing '"<<argv[file]<<"'"<<std::endl;
|
| 293 |
else if ( (ofile.rdstate() & ifstream::failbit ) != 0 )
|
| 294 |
std::cerr << "ERROR opening '"<<argv[file]<<"'! Does the directory exist?"<<std::endl;
|
| 295 |
ofile << config.getComment() << asctime (timeinfo)
|
| 296 |
<< config.getComment()<< "\n"
|
| 297 |
<< config;
|
| 298 |
}
|
| 299 |
}
|
| 300 |
|
| 301 |
|
| 302 |
int main(int argc, char *argv[])
|
| 303 |
{
|
| 304 |
try{
|
| 305 |
rooStat(argc,argv);
|
| 306 |
}
|
| 307 |
catch(exception& e){
|
| 308 |
cout << "Exception catched: " << e.what();
|
| 309 |
}
|
| 310 |
}
|
