114 fDataHist->get(
bin1);
115 double n1 = fDataHist->weight();
116 fDataHist->get(
bin2);
117 double n2 = fDataHist->weight();
126 double n1 = fSparseHist->GetBinContent(
bin1);
127 double n2 = fSparseHist->GetBinContent(
bin2);
135 fChain(chain), fParam(param) {}
137 double xi = fChain->Get(i)->getRealValue(fParam->GetName());
138 double xj = fChain->Get(
j)->getRealValue(fParam->GetName());
176 x[i] =
fAxes[i]->getVal();
205 <<
"Interval type not set. Returning false." << std::endl;
246 "number of variables in axes (" <<
size <<
247 ") doesn't match number of parameters ("
264 <<
"parameters have not been set." << std::endl;
270 "MCMCInterval::CreateKeysPdf: creation of Keys PDF failed: " <<
271 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
272 "in Markov chain." << std::endl;
287 fCutoffVar = std::make_unique<RooRealVar>(
"cutoff",
"cutoff", 0);
289 fProduct = std::make_unique<RooProduct>(
"product",
"Keys PDF & Heaviside Product",
298 coutE(
Eval) <<
"* Error in MCMCInterval::CreateHist(): " <<
299 "Crucial data member was nullptr." << std::endl;
300 coutE(
Eval) <<
"Make sure to fully construct/initialize." << std::endl;
307 "MCMCInterval::CreateHist: creation of histogram failed: " <<
308 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
309 "in Markov chain." << std::endl;
314 fHist = std::make_unique<TH1F>(
"posterior",
"MCMC Posterior Histogram",
318 fHist = std::make_unique<TH2F>(
"posterior",
"MCMC Posterior Histogram",
323 fHist = std::make_unique<TH3F>(
"posterior",
"MCMC Posterior Histogram",
329 coutE(
Eval) <<
"* Error in MCMCInterval::CreateHist() : " <<
330 "TH1* couldn't handle dimension: " <<
fDimension << std::endl;
340 (
static_cast<TH1F&
>(*fHist)).Fill(
entry->getRealValue(
fAxes[0]->GetName()),
343 (
static_cast<TH2F&
>(*fHist)).Fill(
entry->getRealValue(
fAxes[0]->GetName()),
347 (
static_cast<TH3F &
>(*fHist))
367 <<
"Crucial data member was nullptr." << std::endl;
376 min[i] =
fAxes[i]->getMin();
377 max[i] =
fAxes[i]->getMax();
378 bins[i] =
fAxes[i]->numBins();
380 fSparseHist = std::make_unique<THnSparseF>(
"posterior",
"MCMC Posterior Histogram",
381 fDimension, bins.data(), min.data(), max.data());
390 "MCMCInterval::CreateSparseHist: creation of histogram failed: " <<
391 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
392 "in Markov chain." << std::endl;
412 coutE(
Eval) <<
"* Error in MCMCInterval::CreateDataHist(): " <<
413 "Crucial data member was nullptr or empty." << std::endl;
414 coutE(
Eval) <<
"Make sure to fully construct/initialize." << std::endl;
420 "MCMCInterval::CreateDataHist: creation of histogram failed: " <<
421 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
422 "in Markov chain." << std::endl;
440 "Crucial data member (Markov chain) was nullptr." << std::endl;
448 "MCMCInterval::CreateVector: creation of vector failed: " <<
449 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
450 "in Markov chain." << std::endl;
458 for (i = 0; i <
size; i++) {
478 if (
dynamic_cast<RooRealVar *
>(obj) !=
nullptr) {
481 coutE(
Eval) <<
"* Error in MCMCInterval::SetParameters: " << obj->
GetName() <<
" not a RooRealVar*"
501 "Error: Interval type not set" << std::endl;
523 <<
"Fraction must be in the range [0, 1]. "
530 <<
"Error: Can only find a tail-fraction interval for 1-D intervals"
537 <<
"Crucial data member was nullptr." << std::endl;
574 double ll = param->
getMin();
643 coutW(
Eval) <<
"Warning: Integral of Keys PDF came out to " << full
644 <<
" instead of expected value 1. Will continue using this "
645 <<
"factor to normalize further integrals of this PDF." << std::endl;
655 volume *= (var->getMax() - var->getMin());
760 std::vector<Long_t> bins(numBins);
770 for (i = numBins - 1; i >= 0; i--) {
788 for ( ; i >= 0; i--) {
798 for ( ; i < numBins; i++) {
805 if (i == numBins - 1) {
832 std::vector<Int_t> bins(numBins);
841 for (i = numBins - 1; i >= 0; i--) {
860 for ( ; i >= 0; i--) {
871 for ( ; i < numBins; i++) {
879 if (i == numBins - 1) {
904 <<
"not implemented for this type of interval. Returning 0." << std::endl;
920 "Error: Interval type not set" << std::endl;
936 "Error: Interval type not set" << std::endl;
1017 <<
"Sorry, will not compute lower limit unless dimension == 1" << std::endl;
1025 coutE(
Eval) <<
"In MCMCInterval::LowerLimitBySparseHist: "
1026 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1027 <<
"steps in the Markov chain. Returning param.getMin()." << std::endl;
1037 for (
Long_t i = 0; i < numBins; i++) {
1061 coutE(
Eval) <<
"In MCMCInterval::LowerLimitByDataHist: "
1062 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1063 <<
"steps in the Markov chain. Returning param.getMin()." << std::endl;
1072 for (
Int_t i = 0; i < numBins; i++) {
1094 <<
"Sorry, will not compute upper limit unless dimension == 1" << std::endl;
1102 coutE(
Eval) <<
"In MCMCInterval::UpperLimitBySparseHist: "
1103 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1104 <<
"steps in the Markov chain. Returning param.getMax()." << std::endl;
1114 for (
Long_t i = 0; i < numBins; i++) {
1138 coutE(
Eval) <<
"In MCMCInterval::UpperLimitByDataHist: "
1139 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1140 <<
"steps in the Markov chain. Returning param.getMax()." << std::endl;
1149 for (
Int_t i = 0; i < numBins; i++) {
1177 coutE(
Eval) <<
"in MCMCInterval::LowerLimitByKeys(): "
1178 <<
"couldn't find lower limit, check that the number of burn in "
1179 <<
"steps < number of total steps in the Markov chain. Returning "
1180 <<
"param.getMin()" << std::endl;
1189 for (
Int_t i = 0; i < numBins; i++) {
1217 coutE(
Eval) <<
"in MCMCInterval::UpperLimitByKeys(): "
1218 <<
"couldn't find upper limit, check that the number of burn in "
1219 <<
"steps < number of total steps in the Markov chain. Returning "
1220 <<
"param.getMax()" << std::endl;
1229 for (
Int_t i = 0; i < numBins; i++) {
1256 coutE(
Eval) <<
"in MCMCInterval::KeysMax(): "
1257 <<
"couldn't find Keys max value, check that the number of burn in "
1258 <<
"steps < number of total steps in the Markov chain. Returning 0"
1266 for (
Int_t i = 0; i < numBins; i++) {
1317 <<
"confidence level not set " << std::endl;
1319 if (
fHist ==
nullptr)
1322 if (
fHist ==
nullptr) {
1327 return static_cast<TH1*
>(
fHist->Clone(
"MCMCposterior_hist"));
1336 <<
"confidence level not set " << std::endl;
1355 <<
"confidence level not set " << std::endl;
1389 return (std::abs(
a -
b) < std::abs(
fDelta * (
a +
b)/2));
1425 if (!
fAxes[i]->getBinning(
nullptr,
false,
false).isUniform()) {
1448 fKeysDataHist = std::make_unique<RooDataHist>(
"_productDataHist",
1449 "Keys PDF & Heaviside Product Data Hist",
fParameters);
1469 coutE(
Eval) <<
"MCMCInterval: size is wrong, parameters don't match" << std::endl;
1473 coutE(
Eval) <<
"MCMCInterval: size is ok, but parameters don't match" << std::endl;
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
TRObject operator()(const T1 &t1) const
double getRealValue(const char *name, double defVal=0.0, bool verbose=false) const
Get value of a RooAbsReal stored in set with given name.
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
RooArgList is a container object that can hold multiple RooAbsArg objects.
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Container class to hold N-dimensional binned data.
Container class to hold unbinned data.
RooFit::OwningPtr< RooDataHist > binnedClone(const char *newName=nullptr, const char *newTitle=nullptr) const
Return binned clone of this dataset.
Generic N-dimensional implementation of a kernel estimation p.d.f.
Provides numeric constants used in RooFit.
static constexpr double infinity()
Return internal infinity representation.
Represents the product of a given set of RooAbsReal objects.
Variable that can be changed from the outside.
const RooAbsBinning & getBinning(const char *name=nullptr, bool verbose=true, bool createOnTheFly=false) const override
Return binning definition with name.
void setBins(Int_t nBins, const char *name=nullptr)
Create a uniform binning under name 'name' for this variable.
ConfInterval is an interface class for a generic interval in the RooStats framework.
virtual void CreateDataHist()
virtual void DetermineByDataHist()
virtual void DetermineShortestInterval()
virtual double GetActualConfidenceLevel()
virtual double GetKeysPdfCutoff() { return fKeysCutoff; }
double fKeysConfLevel
the actual conf level determined by keys
virtual void CreateVector(RooRealVar *param)
std::unique_ptr< Heaviside > fHeaviside
the Heaviside function
virtual double LowerLimitByDataHist(RooRealVar ¶m)
determine lower limit using histogram
double fDelta
topCutoff (a) considered == bottomCutoff (b) iff (std::abs(a - b) < std::abs(fDelta * (a + b)/2)); Th...
double fConfidenceLevel
Requested confidence level (eg. 0.95 for 95% CL)
virtual double UpperLimitBySparseHist(RooRealVar ¶m)
determine upper limit using histogram
std::unique_ptr< MarkovChain > fChain
the markov chain
enum IntervalType fIntervalType
virtual double UpperLimit(RooRealVar ¶m)
get the highest value of param that is within the confidence interval
std::unique_ptr< RooDataHist > fDataHist
the binned Markov Chain data
virtual void SetAxes(RooArgList &axes)
Set which parameters go on which axis.
virtual void DetermineByKeys()
double fTFLower
lower limit of the tail-fraction interval
bool fUseSparseHist
whether to use sparse hist (vs. RooDataHist)
double fVecWeight
sum of weights of all entries in fVector
double fLeftSideTF
left side tail-fraction for interval
bool AcceptableConfLevel(double confLevel)
virtual void DetermineBySparseHist()
double GetKeysMax()
Determine the approximate maximum value of the Keys PDF.
virtual double LowerLimitBySparseHist(RooRealVar ¶m)
determine lower limit using histogram
bool WithinDeltaFraction(double a, double b)
double fKeysCutoff
cutoff keys pdf value to be in interval
virtual double LowerLimitShortest(RooRealVar ¶m)
get the lower limit of param in the shortest confidence interval Note that this works better for some...
virtual double LowerLimitTailFraction(RooRealVar ¶m)
determine lower limit of the lower confidence interval
virtual TH1 * GetPosteriorHist()
set the number of bins to use (same for all axes, for now) virtual void SetNumBins(Int_t numBins);
std::unique_ptr< RooRealVar > fCutoffVar
cutoff variable to use for integrating keys pdf
Int_t fDimension
number of variables
void SetConfidenceLevel(double cl) override
set the desired confidence level (see GetActualConfidenceLevel()) Note: calling this function trigger...
MCMCInterval(const char *name=nullptr)
default constructor
std::unique_ptr< RooProduct > fProduct
the (keysPdf * heaviside) product
std::vector< RooRealVar * > fAxes
array of pointers to RooRealVars representing
virtual void DetermineTailFractionInterval()
double fTFConfLevel
the actual conf level of tail-fraction interval
double fHistConfLevel
the actual conf level determined by hist
virtual void CreateSparseHist()
double fFull
Value of intergral of fProduct.
std::unique_ptr< RooDataHist > fKeysDataHist
data hist representing product
virtual double UpperLimitByDataHist(RooRealVar ¶m)
determine upper limit using histogram
virtual void DetermineInterval()
std::unique_ptr< TH1 > fHist
the binned Markov Chain data
virtual double CalcConfLevel(double cutoff, double full)
virtual double LowerLimitByHist(RooRealVar ¶m)
determine lower limit using histogram
virtual double LowerLimit(RooRealVar ¶m)
get the lowest value of param that is within the confidence interval
virtual double UpperLimitShortest(RooRealVar ¶m)
get the upper limit of param in the confidence interval Note that this works better for some distribu...
bool CheckParameters(const RooArgSet &point) const override
check if parameters are correct. (dummy implementation to start)
std::vector< Int_t > fVector
vector containing the Markov chain data
virtual double LowerLimitByKeys(RooRealVar ¶m)
determine lower limit in the shortest interval by using keys pdf
virtual RooProduct * GetPosteriorKeysProduct()
Get a clone of the (keyspdf * heaviside) product of the posterior.
double fTFUpper
upper limit of the tail-fraction interval
double fHistCutoff
cutoff bin size to be in interval
RooArgSet * GetParameters() const override
return a set containing the parameters of this interval the caller owns the returned RooArgSet*
std::unique_ptr< RooNDKeysPdf > fKeysPdf
the kernel estimation pdf
virtual double UpperLimitTailFraction(RooRealVar ¶m)
determine upper limit of the lower confidence interval
bool fIsHistStrict
whether the specified confidence level is a floor for the actual confidence level (strict),...
virtual RooNDKeysPdf * GetPosteriorKeysPdf()
Get a clone of the keys pdf of the posterior.
virtual double UpperLimitByHist(RooRealVar ¶m)
determine upper limit using histogram
virtual void CreateKeysPdf()
virtual double UpperLimitByKeys(RooRealVar ¶m)
determine upper limit in the shortest interval by using keys pdf
Int_t fNumBurnInSteps
number of steps to discard as burn in, starting from the first
std::unique_ptr< THnSparse > fSparseHist
the binned Markov Chain data
virtual void CreateHist()
bool fUseKeys
whether to use kernel estimation
RooArgSet fParameters
parameters of interest for this interval
bool IsInInterval(const RooArgSet &point) const override
determine whether this point is in the confidence interval
virtual void DetermineByHist()
virtual void SetParameters(const RooArgSet ¶meters)
Set the parameters of interest for this interval and change other internal data members accordingly.
virtual void CreateKeysDataHist()
double fEpsilon
acceptable error for Keys interval determination
virtual double GetKeysPdfCutoff()
get the cutoff RooNDKeysPdf value for being considered in the confidence interval
virtual double GetHistCutoff()
get the cutoff bin height for being considered in the confidence interval
Stores the steps in a Markov Chain of points.
1-D histogram with a float per channel (see TH1 documentation)
TH1 is the base class of all histogram classes in ROOT.
2-D histogram with a float per channel (see TH1 documentation)
3-D histogram with a float per channel (see TH1 documentation)
Efficient multidimensional histogram.
const char * GetName() const override
Returns name of object.
RooCmdArg SelectVars(const RooArgSet &vars)
RooCmdArg EventRange(Int_t nStart, Int_t nStop)
RooCmdArg NormSet(Args_t &&... argsOrArgSet)
The namespace RooFit contains mostly switches that change the behaviour of functions of PDFs (or othe...
Namespace for the RooStats classes.
void SetParameters(const RooArgSet *desiredVals, RooArgSet *paramsToChange)
CompareDataHistBins(RooDataHist *hist)
CompareSparseHistBins(THnSparse *hist)
CompareVectorIndices(MarkovChain *chain, RooRealVar *param)
static uint64_t sum(uint64_t i)