Skip to content
Snippets Groups Projects
Commit e28186c2 authored by Lorenzo Moneta's avatar Lorenzo Moneta
Browse files

remove the rs50* tutorials 

git-svn-id: http://root.cern.ch/svn/root/trunk@37545 27541ba8-7e3a-0410-8455-c3a389f83636
parent a9467e1a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 0 additions and 751 deletions
tutorials/roostats/rs500a_PrepareWorkspace_Poisson.C deleted 100644 → 0
+ 0
− 80
View file @ a9467e1a
//////////////////////////////////////////////////////////////////////////
//
// RooStats tutorial macro #500a
// 2009/08 - Nils Ruthmann, Gregory Schott
//
// Prepare a workspace (stored in a ROOT file) containing a models,
// data and other objects needed to run statistical classes in
// RooStats.
//
// In this macro a PDF model is built for a counting analysis. A
// certain number of events are observed (this can be enforced or left
// free) while a number of background events is expected. In this
// macro, no systematic uncertainty is considered (see
// rs500b_PrepareWorkspace_Poisson_withSystematics.C) The parameter of
// interest is the signal yield and we assume for it a flat prior.
// All needed objects are stored in a ROOT file (within a RooWorkspace
// container); this ROOT file can then be fed as input to various
// statistical methods.
//
// root -q -x -l 'rs500a_PrepareWorkspace_Poisson.C()'
//
//////////////////////////////////////////////////////////////////////////
#include "RooAbsPdf.h"
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "RooPlot.h"
#include "RooWorkspace.h"
#include "TFile.h"
using namespace RooFit;
// prepare the workspace
// type = 0 : binned data with fixed number of total events
// type = 1 : binned data with total events fluctuationg with Poisson statistics
// type = 2 : binned data without any bin-by bin fluctuations (Asimov data)
void rs500a_PrepareWorkspace_Poisson( TString fileName = "WS_Poisson.root", int type = 1 )
{
// use a RooWorkspace to store the pdf models, prior informations, list of parameters,...
RooWorkspace myWS("myWS");
// Observable
myWS.factory("x[0,0,1]") ;
// Pdf in observable,
myWS.factory("Uniform::sigPdf(x)") ;
myWS.factory("Uniform::bkgPdf(x)") ;
myWS.factory("SUM::model(S[0,0,60]*sigPdf,B[10]*bkgPdf)") ;
// Background only pdf
myWS.factory("ExtendPdf::modelBkg(bkgPdf,B)") ;
// Prior
myWS.factory("Uniform::priorPOI(S)") ;
// Definition of observables and parameters of interest
myWS.defineSet("observables","x");
myWS.defineSet("POI","S");
// Generate data
RooAbsData* data = 0;
// binned data with fixed number of events
if (type ==0) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),myWS.var("S")->getVal(),Name("data"));
// binned data with Poisson fluctuations
if (type ==1) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),Extended(),Name("data"));
// Asimov data: binned data without any fluctuations (average case)
if (type ==2) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),Name("data"),ExpectedData());
myWS.import(*data) ;
// control plot of the generated data
RooPlot* plot = myWS.var("x")->frame();
data->plotOn(plot);
plot->DrawClone();
myWS.writeToFile(fileName);
std::cout << "\nRooFit model initialized and stored in " << fileName << std::endl;
}
tutorials/roostats/rs500b_PrepareWorkspace_Poisson_withSystematics.C deleted 100644 → 0
+ 0
− 84
View file @ a9467e1a
//////////////////////////////////////////////////////////////////////////
//
// RooStats tutorial macro #500b
// 2009/08 - Nils Ruthmann, Gregory Schott
//
// Prepare a workspace (stored in a ROOT file) containing a models,
// data and other objects needed to run statistical classes in
// RooStats.
//
// In this macro a PDF model is built for a counting analysis. A
// certain number of events are observed (this can be enforced or left
// free) while a number of background events is expected. It is also
// assumed there is a systematic uncertainty on the number of expected
// background events. The parameter of interest is the signal yield
// and we assume for it a flat prior. All needed objects are stored
// in a ROOT file (within a RooWorkspace container); this ROOT file
// can then be fed as input to various statistical methods.
//
// root -q -x -l 'rs500b_PrepareWorkspace_Poisson_withSystematics.C()'
//
//////////////////////////////////////////////////////////////////////////
#include "RooAbsPdf.h"
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "RooPlot.h"
#include "RooWorkspace.h"
#include "TFile.h"
using namespace RooFit;
// prepare the workspace
// type = 0 : binned data with fixed number of total events
// type = 1 : binned data with N with POisson fluctuations
// type = 2 : binned data without any bin-by bin fluctuations (Asimov data)
void rs500b_PrepareWorkspace_Poisson_withSystematics( TString fileName = "WS_Poisson_withSystematics.root", int type = 1 )
{
// use a RooWorkspace to store the pdf models, prior informations, list of parameters,...
RooWorkspace myWS("myWS");
// Observable
myWS.factory("x[0,0,1]") ;
// Pdf in observable,
myWS.factory("Uniform::sigPdf(x)") ;
myWS.factory("Uniform::bkgPdf(x)") ;
myWS.factory("SUM::model(S[100,0,1500]*sigPdf,B[1000,0,3000]*bkgPdf)") ;
// Background only pdf
myWS.factory("ExtendPdf::modelBkg(bkgPdf,B)") ;
// Priors
myWS.factory("Gaussian::priorNuisance(B,1000,200)") ;
myWS.factory("Uniform::priorPOI(S)") ;
// Definition of observables and parameters of interest
myWS.defineSet("observables","x");
myWS.defineSet("POI","S");
myWS.defineSet("parameters","B") ;
// Generate data
RooAbsData* data = 0;
// binned data with fixed number of events
if (type ==0) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),myWS.var("S")->getVal(),Name("data"));
// binned data with Poisson fluctuations
if (type ==1) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),Extended(),Name("data"));
// Asimov data: binned data without any fluctuations (average case)
if (type == 2) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),Name("data"),ExpectedData());
myWS.import(*data) ;
myWS.writeToFile(fileName);
std::cout << "\nRooFit model initialized and stored in " << fileName << std::endl;
// control plot of the generated data
RooPlot* plot = myWS.var("x")->frame();
data->plotOn(plot);
plot->DrawClone();
}
tutorials/roostats/rs500c_PrepareWorkspace_GaussOverFlat.C deleted 100644 → 0
+ 0
− 81
View file @ a9467e1a
//////////////////////////////////////////////////////////////////////////
//
// RooStats tutorial macro #500c
// 2009/08 - Nils Ruthmann, Gregory Schott
//
// Prepare a workspace (stored in a ROOT file) containing a models,
// data and other objects needed to run statistical classes in
// RooStats.
//
// In this macro a PDF model is built assuming signal has a Gaussian
// PDF and the background a flat PDF. The parameter of interest is
// the signal yield and we assume for it a flat prior. In this macro,
// no systematic uncertainty is considered (see
// rs500d_PrepareWorkspace_GaussOverFlat_withSystematics.C). All needed
// objects are stored in a ROOT file (within a RooWorkspace
// container); this ROOT file can then be fed as input to various
// statistical methods.
//
// root -q -x -l 'rs500c_PrepareWorkspace_GaussOverFlat.C()'
//
//////////////////////////////////////////////////////////////////////////
#include "RooAbsPdf.h"
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "RooDataSet.h"
#include "RooPlot.h"
#include "RooWorkspace.h"
#include "TFile.h"
using namespace RooFit;
// prepare the workspace
// type = 0 : unbinned data with total number of events fluctuating using Poisson statistics
// type = 1 : binned data with total number of events fluctuating using Poisson statistics
// type = 2 : binned data without any bin-bin fuctuation (Asimov data)
void rs500c_PrepareWorkspace_GaussOverFlat( TString fileName = "WS_GaussOverFlat.root", int type= 1 )
{
// use a RooWorkspace to store the pdf models, prior informations, list of parameters,...
RooWorkspace myWS("myWS");
// Observable
myWS.factory("mass[0,500]") ;
// Pdf in observable,
myWS.factory("Gaussian::sigPdf(mass,200,50)") ;
myWS.factory("Uniform::bkgPdf(mass)") ;
myWS.factory("SUM::model(S[20,0,60]*sigPdf,B[10]*bkgPdf)") ;
// Background only pdf
myWS.factory("ExtendPdf::modelBkg(bkgPdf,B)") ;
// Priors
myWS.factory("Uniform::priorPOI(S)") ;
// Definition of observables and parameters of interest
myWS.defineSet("observables","mass");
myWS.defineSet("POI","S");
// Generate data
RooAbsData* data = 0;
// unbinned data with Poisson fluctuations for total number of events
if (type == 0) data = myWS.pdf("model")->generate(*myWS.set("observables"),Extended(),Name("data"));
// binned data with Poisson fluctuations for total number of events
if (type == 1) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),Extended(),Name("data"));
// binned without any fluctuations (average case)
if (type == 2) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),Name("data"),ExpectedData());
myWS.import(*data) ;
myWS.writeToFile(fileName);
std::cout << "\nRooFit model initialized and stored in " << fileName << std::endl;
// control plot of the generated data
RooPlot* plot = myWS.var("mass")->frame();
data->plotOn(plot);
plot->DrawClone();
}
tutorials/roostats/rs500d_PrepareWorkspace_GaussOverFlat_withSystematics.C deleted 100644 → 0
+ 0
− 86
View file @ a9467e1a
//////////////////////////////////////////////////////////////////////////
//
// RooStats tutorial macro #500d
// 2009/08 - Nils Ruthmann, Gregory Schott
//
// Prepare a workspace (stored in a ROOT file) containing a models,
// data and other objects needed to run statistical classes in
// RooStats.
//
// In this macro a PDF model is built assuming signal has a Gaussian
// PDF and the background a flat PDF. The parameter of interest is
// the signal yield and we assume for it a flat prior. It is shown
// how two types of systematics uncertainties can be expressed; those
// are a sytematic uncertainty on the background yield and another on
// one of the parameters (sigma) of the signal shape. All needed
// objects are stored in a ROOT file (within a RooWorkspace
// container); this ROOT file can then be fed as input to various
// statistical methods.
//
// root -q -x -l 'rs500d_PrepareWorkspace_GaussOverFlat_withSystematics.C()'
//
//////////////////////////////////////////////////////////////////////////
#include "RooAbsPdf.h"
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "RooDataSet.h"
#include "RooPlot.h"
#include "RooWorkspace.h"
#include "TFile.h"
using namespace RooFit;
// prepare the workspace
// type = 0 : unbinned data with total number of events fluctuating using Poisson statistics
// type = 1 : binned data with total number of events fluctuating using Poisson statistics
// type = 2 : binned data without any bin-bin fuctuation (Asimov data)
void rs500d_PrepareWorkspace_GaussOverFlat_withSystematics( TString fileName = "WS_GaussOverFlat_withSystematics.root", int type = 1 )
{
// use a RooWorkspace to store the pdf models, prior informations, list of parameters,...
RooWorkspace myWS("myWS");
// Observable
myWS.factory("mass[0,500]") ;
// Pdf in observable
myWS.factory("Gaussian::sigPdf(mass,200,sigSigma[0,100])") ;
myWS.factory("Uniform::bkgPdf(mass)") ;
myWS.factory("SUM::model(S[1000,0,100000]*sigPdf,B[1000,0,100000]*bkgPdf)") ;
// Background only pdf
myWS.factory("ExtendPdf::modelBkg(bkgPdf,B)") ;
// Priors
myWS.factory("Gaussian::prior_sigSigma(sigSigma,50,5)") ;
myWS.factory("Gaussian::prior_B(B,1000,200)") ;
myWS.factory("PROD::priorNuisance(prior_sigSigma,prior_B)") ;
myWS.factory("Uniform::priorPOI(S)") ;
// Definition of observables and parameters of interest
myWS.defineSet("observables","mass");
myWS.defineSet("parameters","B,sigSigma");
myWS.defineSet("POI","S");
// Generate data
RooAbsData* data = 0;
// unbinned data with Poisson fluctuations for total number of events
if (type == 0) data = myWS.pdf("model")->generate(*myWS.set("observables"),Extended(),Name("data"));
// binned data with Poisson fluctuations for total number of events
if (type == 1) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),Extended(),Name("data"));
// binned without any fluctuations (average case)
if (type == 2) data = myWS.pdf("model")->generateBinned(*myWS.set("observables"),Name("data"),ExpectedData());
myWS.import(*data) ;
myWS.writeToFile(fileName);
std::cout << "\nRooFit model initialized and stored in " << fileName << std::endl;
// control plot of the generated data
RooPlot* plot = myWS.var("mass")->frame();
data->plotOn(plot);
plot->DrawClone();
}
tutorials/roostats/rs501_ProfileLikelihoodCalculator_limit.C deleted 100644 → 0
+ 0
− 66
View file @ a9467e1a
//////////////////////////////////////////////////////////////////////////
//
// RooStats tutorial macro #501
// 2009/08 - Nils Ruthmann, Gregory Schott
//
// Show how to run the RooStats classes to perform specific tasks. The
// ROOT file containing a workspace holding the models, data and other
// objects needed to run can be prepared with any of the rs500*.C
// tutorial macros.
//
// Compute with ProfileLikelihoodCalculator a 95% CL upper limit on
// the parameter of interest for the given data.
//
//////////////////////////////////////////////////////////////////////////
#include "RooRealVar.h"
#include "RooProdPdf.h"
#include "RooWorkspace.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/LikelihoodIntervalPlot.h"
#include "TFile.h"
using namespace RooFit;
using namespace RooStats;
void rs501_ProfileLikelihoodCalculator_limit( const char* fileName="WS_GaussOverFlat.root" )
{
// Open the ROOT file and import from the workspace the objects needed for this tutorial
TFile* file = new TFile(fileName);
RooWorkspace* myWS = (RooWorkspace*) file->Get("myWS");
RooAbsPdf* modelTmp = myWS->pdf("model");
RooAbsData* data = myWS->data("data");
RooAbsPdf* priorNuisance = myWS->pdf("priorNuisance");
const RooArgSet* POI = myWS->set("POI");
RooRealVar* parameterOfInterest = dynamic_cast<RooRealVar*>(POI->first());
assert(parameterOfInterest);
// If there are nuisance parameters, multiply their prior distribution to the full model
RooAbsPdf* model = modelTmp;
if( priorNuisance!=0 ) model = new RooProdPdf("constrainedModel","Model with nuisance parameters",*modelTmp,*priorNuisance);
// Set up the ProfileLikelihoodCalculator
ProfileLikelihoodCalculator plc(*data, *model, *POI);
// ProfileLikelihoodCalculator usually make intervals: the 95% CL one-sided upper-limit is the same as the two-sided upper-limit of a 90% CL interval
plc.SetTestSize(0.10);
// Pointer to the confidence interval
model->fitTo(*data,SumW2Error(kFALSE)); // <-- problem
LikelihoodInterval* interval = plc.GetInterval();
// Compute the upper limit: a fit is needed first in order to locate the minimum of the -log(likelihood) and ease the upper limit computation
model->fitTo(*data,SumW2Error(kFALSE)); // <-- problem
const double upperLimit = interval->UpperLimit(*parameterOfInterest); // <-- to simplify
file->Close();
// Make a plot of the profile-likelihood and confidence interval
LikelihoodIntervalPlot plot(interval);
plot.Draw();
std::cout << "One sided upper limit at 95% CL: "<< upperLimit << std::endl;
delete model;
}
tutorials/roostats/rs502_ProfileLikelihoodCalculator_significance.C deleted 100644 → 0
+ 0
− 102
View file @ a9467e1a
//////////////////////////////////////////////////////////////////////////
//
// RooStats tutorial macro #502
// 2009/08 - Nils Ruthmann, Gregory Schott
//
// Show how to run the RooStats classes to perform specific tasks. The
// ROOT file containing a workspace holding the models, data and other
// objects needed to run can be prepared with any of the rs500*.C
// tutorial macros.
//
// Compute with ProfileLikelihoodCalculator a 68% CL interval and the
// signal significance for the given data.
//
//////////////////////////////////////////////////////////////////////////
#include "RooRealVar.h"
#include "RooProdPdf.h"
#include "RooWorkspace.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/LikelihoodIntervalPlot.h"
#include "RooStats/HypoTestResult.h"
#include "TFile.h"
using namespace RooFit;
using namespace RooStats;
void rs502_ProfileLikelihoodCalculator_significance( const char* fileName="WS_GaussOverFlat.root" )
{
// Open the ROOT file and import from the workspace the objects needed for this tutorial
TFile* file = new TFile(fileName);
RooWorkspace* myWS = (RooWorkspace*) file->Get("myWS");
RooAbsPdf* modelTmp = myWS->pdf("model");
RooAbsData* data = myWS->data("data");
RooAbsPdf* priorNuisance = myWS->pdf("priorNuisance");
const RooArgSet* POI = myWS->set("POI");
RooRealVar* parameterOfInterest = dynamic_cast<RooRealVar*>(POI->first());
assert(parameterOfInterest);
// If there are nuisance parameters, multiply their prior distribution to the full model
RooAbsPdf* model = modelTmp;
if( priorNuisance!=0 ) model = new RooProdPdf("constrainedModel","Model with nuisance parameters",*modelTmp,*priorNuisance);
//myWS->var("B")->setConstant();
// Set up the ProfileLikelihoodCalculator
ProfileLikelihoodCalculator plc(*data, *model, *POI);
// The 68% CL ProfileLikelihoodCalculator interval correspond to test size of 0.32
plc.SetTestSize(0.32);
// Compute the confidence interval: a fit is needed first in order to locate the minimum of the -log(likelihood) and ease the upper limit computation
model->fitTo(*data,SumW2Error(kFALSE));
// Pointer to the confidence interval
LikelihoodInterval* interval = plc.GetInterval();
const double MLE = parameterOfInterest->getVal();
const double lowerLimit = interval->LowerLimit(*parameterOfInterest);
parameterOfInterest->setVal(MLE);
const double upperLimit = interval->UpperLimit(*parameterOfInterest);
parameterOfInterest->setVal(MLE);
// Make a plot of the profile-likelihood and confidence interval
LikelihoodIntervalPlot plot(interval);
plot.Draw();
// Get the significance using the GetHypoTest function: (a plot is not possible)
// create a copy of the POI parameters to set the values to zero
RooArgSet nullparams;
nullparams.addClone(*parameterOfInterest);
( (RooRealVar *) (nullparams.first()))->setVal(0);
plc.SetNullParameters(nullparams);
HypoTestResult* testresult=plc.GetHypoTest();
const double significance = testresult->Significance();
// // Another way is to use directly
// // the profile-log-likelihood function
// parameterOfInterest->setVal(MLE);
// RooAbsReal* profile = interval->GetLikelihoodRatio();
// profile->getVal();
// //Go to the Bkg Hypothesis
// RooRealVar* myarg = (RooRealVar *) profile->getVariables()->find(parameterOfInterest->GetName()); // <-- cloned!
// myarg->setVal(0);
// Double_t delta_nll = profile->getVal();
// double significance = sqrt(fabs(2*delta_nll));
// if (delta_nll<0) significance *= -1;
std::cout << "68% CL interval: [ " << lowerLimit << " ; " << upperLimit << " ]\n";
std::cout << "significance estimation: " << significance << std::endl;
delete model;
// closing the file will delete the workspace
//file->Close();
}
tutorials/roostats/rs504_ProfileLikelihoodCalculator_averageSignificance.C deleted 100644 → 0
+ 0
− 112
View file @ a9467e1a
//////////////////////////////////////////////////////////////////////////
//
// RooStats tutorial macro #504
// 2009/08 - Nils Ruthmann, Gregory Schott
//
///////////////////////////////////////////////////////////////////////
#include "RooRealVar.h"
#include "RooProdPdf.h"
#include "RooWorkspace.h"
#include "RooRandom.h"
#include "RooMCStudy.h"
#include "RooDataSet.h"
#include "RooPlot.h"
#include "RooDLLSignificanceMCSModule.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/LikelihoodIntervalPlot.h"
#include "RooStats/HypoTestResult.h"
#include "RooStats/UpperLimitMCSModule.h"
#include "TFile.h"
#include "TStopwatch.h"
#include "TCanvas.h"
void rs504_ProfileLikelihoodCalculator_averageSignificance(const char* fname="WS_GaussOverFlat_withSystematics.root",int ntoys=100,const char* outputplot="pll_avSign.ps"){
using namespace RooFit;
using namespace RooStats;
TStopwatch t;
t.Start();
TFile* file =new TFile(fname);
RooWorkspace* my_WS = (RooWorkspace*) file->Get("myWS");
//Import the objects needed
RooAbsPdf* model_naked=my_WS->pdf("model");
RooAbsPdf* priorNuisance=my_WS->pdf("priorNuisance");
const RooArgSet* paramInterestSet=my_WS->set("POI");
RooRealVar* paramInterest= (RooRealVar*) paramInterestSet->first();
const RooArgSet* observable=my_WS->set("observables");
const RooArgSet* nuisanceParam=my_WS->set("parameters");
//If there are nuisance parameters present, multiply their prior to the model
RooAbsPdf* model=model_naked;
if(priorNuisance!=0) {
model=new RooProdPdf("constrainedModel","Model with nuisance parameters",*model_naked,*priorNuisance);
//From now work with the product of both
}
//Save the default values of the parameters:
RooArgSet* parameters=model->getVariables();
RooArgSet* default_parameters=new RooArgSet("default_parameters");
TIterator* it=parameters->createIterator();
RooRealVar* currentparam=(RooRealVar*) it->Next();
do {
default_parameters->addClone(*currentparam,false);
currentparam=(RooRealVar*) it->Next();
}while(currentparam!=0);
if(priorNuisance!=0)
RooFormulaVar nll_nuisance("nllSyst","-TMath::Log(@0)",RooArgList(*priorNuisance));
else
RooFormulaVar nll_nuisance("nllSyst","0",RooArgList(*priorNuisance));
RooRandom::randomGenerator()->SetSeed(110);
//--------------------------------------------------------------------
//ROOMCSTUDY
//For simplicity use RooMCStudy.
RooMCStudy* mcs;
if (nuisanceParam)
mcs = new RooMCStudy(*model,*observable,Extended(kTRUE),
FitOptions(Extended(kTRUE),PrintEvalErrors(-1),Minos(kFALSE)),Constrain(*nuisanceParam)) ;
else
mcs = new RooMCStudy(*model,*observable,Extended(kTRUE),
FitOptions(Extended(kTRUE),Minos(kFALSE),PrintEvalErrors(-1))) ;
//Adding a module which allows to compute the significance in each toy experiment
RooDLLSignificanceMCSModule sigModule(*paramInterest,0) ;
//If there are nuisance parameters present, they should be generated according to their pdf for every new toy experiment.
//this is done using a MCSModule
mcs->addModule(sigModule);
mcs->generateAndFit(ntoys);
TString signstring("significance_nullhypo_");
TH1* mcssign_histo=(TH1F*)mcs->fitParDataSet().createHistogram(signstring+paramInterest->GetName());
TCanvas* c2 =new TCanvas();
c2->Divide(2,2);
c2->cd(1);
mcssign_histo->Draw();
c2->cd(2);
mcs->plotPull(*paramInterest)->Draw();
c2->cd(3);
if (my_WS->var("B")) { mcs->plotParam(*(my_WS->var("B")))->Draw();
c2->cd(4);
mcs->plotError(*(my_WS->var("B")))->Draw();
}
c2->Print(outputplot);
std::cout<<"The average significance after "<<ntoys<<" toys is: "<<mcssign_histo->GetMean()<<std::endl;
//file->Close();
t.Stop();
t.Print();
}
tutorials/roostats/rs505_HybridCalculator_significance.C deleted 100644 → 0
+ 0
− 72
View file @ a9467e1a
#include "RooRandom.h"
#include "RooRealVar.h"
#include "RooProdPdf.h"
#include "RooWorkspace.h"
#include "RooStats/HybridCalculatorOriginal.h"
#include "RooStats/HybridResult.h"
#include "RooStats/HybridPlot.h"
#include "TFile.h"
#include "TStopwatch.h"
#include "TCanvas.h"
using namespace RooFit;
using namespace RooStats;
void rs505_HybridCalculator_significance(const char* fname="WS_GaussOverFlat_withSystematics.root",int ntoys=5000,const char* outputplot="hc_sign_shape_nosyst.pdf"){
RooRandom::randomGenerator()->SetSeed(100);
TStopwatch t;
t.Start();
TFile* file =new TFile(fname);
RooWorkspace* my_WS = (RooWorkspace*) file->Get("myWS");
if (!my_WS) return;
//Import the objects needed
RooAbsPdf* model=my_WS->pdf("model");
RooAbsData* data=my_WS->data("data");
RooAbsPdf* nuisanceTerm=my_WS->pdf("nuisanceTerm");
//const RooArgSet* paramInterestSet=my_WS->set("POI");
//RooRealVar* paramInterest= (RooRealVar*) paramInterestSet->first();
RooAbsPdf* modelBkg=my_WS->pdf("modelBkg");
//const RooArgSet* observable=my_WS->set("observables");
const RooArgSet* nuisanceParam=my_WS->set("parameters");
HybridCalculatorOriginal * hc=new HybridCalculatorOriginal(*data,*model,*modelBkg);
hc->SetNumberOfToys(ntoys);
bool useNuisance=false;
if(nuisanceTerm!=0){
hc->UseNuisance(kTRUE);
hc->SetNuisancePdf(*nuisanceTerm);
useNuisance=true;
cout<<"The following nuisance parameters are taken into account:"<<endl;
nuisanceParam->Print();
hc->SetNuisanceParameters(*nuisanceParam);
}else{
hc->UseNuisance(kFALSE);
}
HybridResult* hcresult=hc->GetHypoTest();
double clsb_data = hcresult->CLsplusb();
double clb_data = hcresult->CLb();
double cls_data = hcresult->CLs();
double data_significance = hcresult->Significance();
cout<<"CL_b:"<<clb_data<<endl;
cout<<"CL_s:"<<cls_data<<endl;
cout<<"CL_sb:"<<clsb_data<<endl;
cout<<"significance:" <<data_significance<<endl;
HybridPlot* hcPlot=hcresult->GetPlot("hcPlot","p Values Plot",100);
TCanvas*c1=new TCanvas();
c1->cd();
hcPlot->Draw();
c1->Print(outputplot);
file->Close();
t.Stop();
t.Print();
}
tutorials/roostats/rs506_HybridCalculator_averageSignificance.C deleted 100644 → 0
+ 0
− 68
View file @ a9467e1a
#include "RooRandom.h"
#include "RooRealVar.h"
#include "RooProdPdf.h"
#include "RooWorkspace.h"
#include "RooStats/HybridCalculatorOriginal.h"
#include "RooStats/HybridResult.h"
#include "RooStats/HybridPlot.h"
#include "TFile.h"
#include "TStopwatch.h"
#include "TCanvas.h"
using namespace RooFit;
using namespace RooStats;
void rs506_HybridCalculator_averageSignificance(const char* fname="WS_GaussOverFlat_withSystematics.root",int ntoys=2500,const char* outputplot="hc_averagesign_counting_syst.ps"){
RooRandom::randomGenerator()->SetSeed(100);
TStopwatch t;
t.Start();
TFile* file =new TFile(fname);
RooWorkspace* my_WS = (RooWorkspace*) file->Get("myWS");
if (!my_WS) return;
//Import the objects needed
RooAbsPdf* model=my_WS->pdf("model");
RooAbsPdf* nuisanceTerm=my_WS->pdf("nuisanceTerm");
//const RooArgSet* paramInterestSet=my_WS->set("paramInterest");
//RooRealVar* paramInterest=(RooRealVar*) paramInterestSet->first();
RooAbsPdf* modelBkg=my_WS->pdf("modelBkg");
const RooArgSet* nuisanceParam=my_WS->set("parameters");
RooArgList observable(*(my_WS->set("observables") ) );
HybridCalculatorOriginal * hc=new HybridCalculatorOriginal(*model,*modelBkg,observable);
hc->SetNumberOfToys(ntoys);
hc->SetTestStatistic(1);
bool useNuisance=false;
if(nuisanceTerm!=0){
hc->UseNuisance(kTRUE);
hc->SetNuisancePdf(*nuisanceTerm);
useNuisance=true;
nuisanceParam->Print();
hc->SetNuisanceParameters(*nuisanceParam);
}else{
hc->UseNuisance(kFALSE);
}
RooRandom::randomGenerator()->SetSeed(0);
HybridResult* hcresult=hc->Calculate(ntoys,useNuisance);
HybridPlot* hcplot = hcresult->GetPlot("HybridPlot","Toy MC Q ",100);
double mean_sb_toys_test_stat = hcplot->GetSBmean();
hcresult->SetDataTestStatistics(mean_sb_toys_test_stat);
double mean_significance = hcresult->Significance();
cout<<"significance:" <<mean_significance<<endl;
hcplot = hcresult->GetPlot("HybridPlot","Toy MC -2ln Q ",100);
TCanvas*c1=new TCanvas();
c1->cd();
hcplot->Draw(outputplot);
c1->Print(outputplot);
file->Close();
t.Stop();
t.Print();
}
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment