129 fDataHist->get(bin1);
130 double n1 = fDataHist->weight();
131 fDataHist->get(bin2);
132 double n2 = fDataHist->weight();
141 double n1 = fSparseHist->GetBinContent(bin1);
142 double n2 = fSparseHist->GetBinContent(bin2);
150 fChain(chain), fParam(param) {}
152 double xi = fChain->Get(i)->getRealValue(fParam->GetName());
153 double xj = fChain->Get(j)->getRealValue(fParam->GetName());
191 x[i] =
fAxes[i]->getVal();
220 <<
"Interval type not set. Returning false." << endl;
261 "number of variables in axes (" <<
size <<
262 ") doesn't match number of parameters ("
279 <<
"parameters have not been set." << endl;
285 "MCMCInterval::CreateKeysPdf: creation of Keys PDF failed: " <<
286 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
287 "in Markov chain." << endl;
317 coutE(
Eval) <<
"* Error in MCMCInterval::CreateHist(): " <<
318 "Crucial data member was nullptr." << endl;
319 coutE(
Eval) <<
"Make sure to fully construct/initialize." << endl;
322 if (
fHist !=
nullptr) {
329 "MCMCInterval::CreateHist: creation of histogram failed: " <<
330 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
331 "in Markov chain." << endl;
336 fHist =
new TH1F(
"posterior",
"MCMC Posterior Histogram",
340 fHist =
new TH2F(
"posterior",
"MCMC Posterior Histogram",
345 fHist =
new TH3F(
"posterior",
"MCMC Posterior Histogram",
351 coutE(
Eval) <<
"* Error in MCMCInterval::CreateHist() : " <<
352 "TH1* couldn't handle dimension: " <<
fDimension << endl;
389 <<
"Crucial data member was nullptr." << endl;
406 fDimension, bins.data(), min.data(), max.data());
415 "MCMCInterval::CreateSparseHist: creation of histogram failed: " <<
416 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
417 "in Markov chain." << endl;
437 coutE(
Eval) <<
"* Error in MCMCInterval::CreateDataHist(): " <<
438 "Crucial data member was nullptr or empty." << endl;
439 coutE(
Eval) <<
"Make sure to fully construct/initialize." << endl;
445 "MCMCInterval::CreateDataHist: creation of histogram failed: " <<
446 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
447 "in Markov chain." << endl;
465 "Crucial data member (Markov chain) was nullptr." << endl;
473 "MCMCInterval::CreateVector: creation of vector failed: " <<
474 "Number of burn-in steps (num steps to ignore) >= number of steps " <<
475 "in Markov chain." << endl;
483 for (i = 0; i <
size; i++) {
500 if (
fAxes !=
nullptr)
505 if (
dynamic_cast<RooRealVar *
>(obj) !=
nullptr) {
508 coutE(
Eval) <<
"* Error in MCMCInterval::SetParameters: " << obj->
GetName() <<
" not a RooRealVar*"
528 "Error: Interval type not set" << endl;
548 if (fLeftSideTF < 0 || fLeftSideTF > 1) {
550 <<
"Fraction must be in the range [0, 1]. "
557 <<
"Error: Can only find a tail-fraction interval for 1-D intervals"
562 if (
fAxes ==
nullptr) {
564 <<
"Crucial data member was nullptr." << endl;
597 double leftTailSum = 0;
598 double rightTailSum = 0;
601 double ll = param->
getMin();
602 double ul = param->
getMax();
616 if (std::abs(leftTailSum +
w - leftTailCutoff) <
617 std::abs(leftTailSum - leftTailCutoff)) {
631 if (std::abs(rightTailSum +
w - rightTailCutoff) <
632 std::abs(rightTailSum - rightTailCutoff)) {
670 coutW(
Eval) <<
"Warning: Integral of Keys PDF came out to " << full
671 <<
" instead of expected value 1. Will continue using this "
672 <<
"factor to normalize further integrals of this PDF." << endl;
681 for (
auto *var : static_range_cast<RooRealVar*>(
fParameters))
682 volume *= (var->getMax() - var->getMin());
684 double topCutoff = full / volume;
685 double bottomCutoff = topCutoff;
692 bool changed =
false;
705 bottomCutoff = topCutoff / 2.0;
726 topCutoff = bottomCutoff * 2.0;
730 coutI(
Eval) <<
"range set: [" << bottomCutoff <<
", " << topCutoff <<
"]"
733 cutoff = (topCutoff + bottomCutoff) / 2.0;
745 bottomCutoff = cutoff;
749 cutoff = (topCutoff + bottomCutoff) / 2.0;
750 coutI(
Eval) <<
"cutoff range: [" << bottomCutoff <<
", "
751 << topCutoff <<
"]" << endl;
787 std::vector<Long_t> bins(numBins);
788 for (
Int_t ibin = 0; ibin < numBins; ibin++)
789 bins[ibin] = (
Long_t)ibin;
797 for (i = numBins - 1; i >= 0; i--) {
815 for ( ; i >= 0; i--) {
825 for ( ; i < numBins; i++) {
832 if (i == numBins - 1) {
859 std::vector<Int_t> bins(numBins);
860 for (
Int_t ibin = 0; ibin < numBins; ibin++)
868 for (i = numBins - 1; i >= 0; i--) {
887 for ( ; i >= 0; i--) {
898 for ( ; i < numBins; i++) {
906 if (i == numBins - 1) {
931 <<
"not implemented for this type of interval. Returning 0." << endl;
947 "Error: Interval type not set" << endl;
963 "Error: Interval type not set" << endl;
1044 <<
"Sorry, will not compute lower limit unless dimension == 1" << endl;
1052 coutE(
Eval) <<
"In MCMCInterval::LowerLimitBySparseHist: "
1053 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1054 <<
"steps in the Markov chain. Returning param.getMin()." << endl;
1062 double lowerLimit = param.
getMax();
1064 for (
Long_t i = 0; i < numBins; i++) {
1067 if (val < lowerLimit)
1088 coutE(
Eval) <<
"In MCMCInterval::LowerLimitByDataHist: "
1089 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1090 <<
"steps in the Markov chain. Returning param.getMin()." << endl;
1097 double lowerLimit = param.
getMax();
1099 for (
Int_t i = 0; i < numBins; i++) {
1103 if (val < lowerLimit)
1121 <<
"Sorry, will not compute upper limit unless dimension == 1" << endl;
1129 coutE(
Eval) <<
"In MCMCInterval::UpperLimitBySparseHist: "
1130 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1131 <<
"steps in the Markov chain. Returning param.getMax()." << endl;
1139 double upperLimit = param.
getMin();
1141 for (
Long_t i = 0; i < numBins; i++) {
1144 if (val > upperLimit)
1165 coutE(
Eval) <<
"In MCMCInterval::UpperLimitByDataHist: "
1166 <<
"couldn't determine cutoff. Check that num burn in steps < num "
1167 <<
"steps in the Markov chain. Returning param.getMax()." << endl;
1174 double upperLimit = param.
getMin();
1176 for (
Int_t i = 0; i < numBins; i++) {
1180 if (val > upperLimit)
1204 coutE(
Eval) <<
"in MCMCInterval::LowerLimitByKeys(): "
1205 <<
"couldn't find lower limit, check that the number of burn in "
1206 <<
"steps < number of total steps in the Markov chain. Returning "
1207 <<
"param.getMin()" << endl;
1214 double lowerLimit = param.
getMax();
1216 for (
Int_t i = 0; i < numBins; i++) {
1220 if (val < lowerLimit)
1244 coutE(
Eval) <<
"in MCMCInterval::UpperLimitByKeys(): "
1245 <<
"couldn't find upper limit, check that the number of burn in "
1246 <<
"steps < number of total steps in the Markov chain. Returning "
1247 <<
"param.getMax()" << endl;
1254 double upperLimit = param.
getMin();
1256 for (
Int_t i = 0; i < numBins; i++) {
1260 if (val > upperLimit)
1283 coutE(
Eval) <<
"in MCMCInterval::KeysMax(): "
1284 <<
"couldn't find Keys max value, check that the number of burn in "
1285 <<
"steps < number of total steps in the Markov chain. Returning 0"
1293 for (
Int_t i = 0; i < numBins; i++) {
1333 double confLevel = integral->getVal(
fParameters) / full;
1334 coutI(
Eval) <<
"cutoff = " << cutoff <<
", conf = " << confLevel << endl;
1344 <<
"confidence level not set " << endl;
1346 if (
fHist ==
nullptr)
1349 if (
fHist ==
nullptr) {
1354 return static_cast<TH1*
>(
fHist->
Clone(
"MCMCposterior_hist"));
1363 <<
"confidence level not set " << endl;
1382 <<
"confidence level not set " << endl;
1416 return (std::abs(
a -
b) < std::abs(
fDelta * (
a +
b)/2));
1423 if (
fAxes ==
nullptr)
1450 bool tempChangeBinning =
true;
1452 if (!
fAxes[i]->getBinning(
nullptr,
false,
false).isUniform()) {
1453 tempChangeBinning =
false;
1462 tempChangeBinning =
false;
1464 if (tempChangeBinning) {
1476 "Keys PDF & Heaviside Product Data Hist",
fParameters);
1479 if (tempChangeBinning) {
1496 coutE(
Eval) <<
"MCMCInterval: size is wrong, parameters don't match" << std::endl;
1500 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
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
THnSparseT< TArrayF > THnSparseF
TRObject operator()(const T1 &t1) const
TObject * Clone(const char *newname=nullptr) const override
Make a clone of an object using the Streamer facility.
Int_t numBins() const
Return number of bins.
bool equals(const RooAbsCollection &otherColl) const
Check if this and other collection have identically-named contents.
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
RooFit::OwningPtr< RooAbsData > reduce(const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create a reduced copy of this dataset.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
Int_t numBins(const char *rangeName=nullptr) const override
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.
RooDataHist * fillDataHist(RooDataHist *hist, const RooArgSet *nset, double scaleFactor, bool correctForBinVolume=false, bool showProgress=false) const
Fill a RooDataHist with values sampled from this function at the bin centers.
RooFit::OwningPtr< RooAbsReal > createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2={}, const RooCmdArg &arg3={}, const RooCmdArg &arg4={}, const RooCmdArg &arg5={}, const RooCmdArg &arg6={}, const RooCmdArg &arg7={}, const RooCmdArg &arg8={}) const
Create an object that represents the integral of the function over one or more observables listed in ...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
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.
double sum(bool correctForBinSize, bool inverseCorr=false) const
Return the sum of the weights of all bins in the histogram.
Int_t getIndex(const RooAbsCollection &coord, bool fast=false) const
Calculate bin number of the given coordinates.
double weight(std::size_t i) const
Return weight of i-th bin.
const RooArgSet * get() const override
Get bin centre of current bin.
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.
void setVal(double value) override
Set value of variable to 'value'.
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.
Represents the Heaviside function.
MCMCInterval is a concrete implementation of the RooStats::ConfInterval interface.
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)
RooDataHist * fDataHist
the binned Markov Chain data
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)
TH1 * fHist
the binned Markov Chain data
virtual double UpperLimitBySparseHist(RooRealVar ¶m)
determine upper limit using histogram
enum IntervalType fIntervalType
virtual double UpperLimit(RooRealVar ¶m)
get the highest value of param that is within the confidence interval
RooRealVar * fCutoffVar
cutoff variable to use for integrating keys pdf
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.
RooProduct * fProduct
the (keysPdf * heaviside) product
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);
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
RooRealVar ** fAxes
array of pointers to RooRealVars representing the axes of the histogram fAxes[0] represents x-axis,...
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()
Heaviside * fHeaviside
the Heaviside function
double fFull
Value of intergral of fProduct.
virtual double UpperLimitByDataHist(RooRealVar ¶m)
determine upper limit using histogram
virtual void DetermineInterval()
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...
THnSparse * fSparseHist
the binned Markov Chain data
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
MarkovChain * fChain
the markov chain
RooArgSet * GetParameters() const override
return a set containing the parameters of this interval the caller owns the returned RooArgSet*
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
virtual void CreateHist()
bool fUseKeys
whether to use kernel estimation
RooDataHist * fKeysDataHist
data hist representing product
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()
RooNDKeysPdf * fKeysPdf
the kernel estimation pdf
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.
virtual const RooArgSet * Get(Int_t i) const
get the entry at position i
virtual double Weight() const
get the weight of the current (last indexed) entry
virtual const RooDataSet * GetAsConstDataSet() const
virtual Int_t Size() const
get the number of steps in the chain
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
1-D histogram with a float per channel (see TH1 documentation)
TH1 is the base class of all histogram classes in ROOT.
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
2-D histogram with a float per channel (see TH1 documentation)
3-D histogram with a float per channel (see TH1 documentation)
Long64_t Fill(const Double_t *x, Double_t w=1.)
TAxis * GetAxis(Int_t dim) const
Efficient multidimensional histogram.
Double_t GetBinContent(const Int_t *idx) const
Forwards to THnBase::GetBinContent() overload.
Long64_t GetBin(const Int_t *idx) const override
void Sumw2() override
Enable calculation of errors.
Long64_t GetNbins() const override
virtual void SetTitle(const char *title="")
Set the title of the TNamed.
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)