73 :
RooAbsPdf(
name, title), _varList(
"varList",
"List of variables", this),
74 _rhoList(
"rhoList",
"List of rho parameters", this), _data(&data), _options(options), _widthFactor(rho),
75 _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
79 for (
const auto var : varList) {
81 coutE(InputArguments) <<
"RooNDKeysPdf::ctor(" <<
GetName() <<
") ERROR: variable " << var->
GetName()
82 <<
" is not of type RooAbsReal" << endl ;
97 :
RooAbsPdf(
name, title), _varList(
"varList",
"List of variables", this),
98 _rhoList(
"rhoList",
"List of rho parameters", this), _ownedData(createDatasetFromHist(varList, hist)), _data(_ownedData.get()),
99 _options(options), _widthFactor(rho), _nSigma(nSigma), _rotate(rotate),
100 _sortInput(sortInput), _nAdpt(1)
102 for (
const auto var : varList) {
104 coutE(InputArguments) <<
"RooNDKeysPdf::ctor(" <<
GetName() <<
") ERROR: variable " << var->
GetName()
105 <<
" is not of type RooAbsReal" << endl;
120 :
RooAbsPdf(
name, title), _varList(
"varList",
"List of variables", this),
121 _rhoList(
"rhoList",
"List of rho parameters", this), _data(&data), _options(options), _widthFactor(-1.0),
122 _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
124 for (
const auto var : varList) {
126 coutE(InputArguments) <<
"RooNDKeysPdf::ctor(" <<
GetName() <<
") ERROR: variable " << var->
GetName()
127 <<
" is not of type RooAbsReal" << endl;
136 coutE(InputArguments)
137 <<
"ERROR: RooNDKeysPdf::RooNDKeysPdf() : The vector-size of rho is different from that of varList."
138 <<
"Unable to create the PDF." << endl;
159 :
RooAbsPdf(
name, title), _varList(
"varList",
"List of variables", this),
160 _rhoList(
"rhoList",
"List of rho parameters", this), _data(&data), _options(options), _widthFactor(-1.0),
161 _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
163 for (
const auto var : varList) {
165 coutE(InputArguments) <<
"RooNDKeysPdf::ctor(" <<
GetName() <<
") ERROR: variable " << var->
GetName()
166 <<
" is not of type RooAbsReal" << endl;
175 for (
unsigned int i=0; i < rhoList.
size(); ++i) {
176 const auto rho = rhoList.
at(i);
178 coutE(InputArguments) <<
"RooNDKeysPdf::ctor(" <<
GetName() <<
") ERROR: parameter " << rho->
GetName()
179 <<
" is not of type RooRealVar" << endl;
188 coutE(InputArguments) <<
"ERROR: RooNDKeysPdf::RooNDKeysPdf() : The size of rhoList is different from varList."
189 <<
"Unable to create the PDF." << endl;
205 :
RooAbsPdf(
name, title), _varList(
"varList",
"List of variables", this),
206 _rhoList(
"rhoList",
"List of rho parameters", this), _ownedData(createDatasetFromHist(varList, hist)), _data(_ownedData.get()),
207 _options(options), _widthFactor(-1), _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput),
210 for (
const auto var : varList) {
212 coutE(InputArguments) <<
"RooNDKeysPdf::ctor(" <<
GetName() <<
") ERROR: variable " << var->
GetName()
213 <<
" is not of type RooAbsReal" << endl;
223 for (
unsigned int i=0; i < rhoList.
size(); ++i) {
224 const auto rho = rhoList.
at(i);
226 coutE(InputArguments) <<
"RooNDKeysPdf::ctor(" <<
GetName() <<
") ERROR: parameter " << rho->
GetName()
227 <<
" is not of type RooRealVar" << endl;
235 coutE(InputArguments) <<
"ERROR: RooNDKeysPdf::RooNDKeysPdf() : The size of rhoList is different from varList."
236 <<
"Unable to create the PDF." << endl;
252 :
RooAbsPdf(
name, title), _varList(
"varList",
"List of variables", this),
253 _rhoList(
"rhoList",
"List of rho parameters", this), _data(&data), _options(
"a"), _widthFactor(rho),
254 _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
261 coutW(InputArguments) <<
"RooNDKeysPdf::RooNDKeysPdf() : Warning : asymmetric mirror(s) no longer supported."
275 :
RooAbsPdf(
name, title), _varList(
"varList",
"List of variables", this),
276 _rhoList(
"rhoList",
"List of rho parameters", this), _data(&data), _options(options), _widthFactor(rho),
277 _nSigma(nSigma), _rotate(rotate), _sortInput(sortInput), _nAdpt(1)
290 :
RooAbsPdf(other,
name), _varList(
"varList", this, other._varList), _rhoList(
"rhoList", this, other._rhoList),
291 _ownedData(other._ownedData ? new
RooDataSet(*other._ownedData) : nullptr),
292 _data(_ownedData ? _ownedData.get() : other._data), _options(other._options), _widthFactor(other._widthFactor),
293 _nSigma(other._nSigma), _rotate(other._rotate), _sortInput(other._sortInput),
431 cxcoutD(InputArguments) <<
"RooNDKeysPdf::setOptions() options = " <<
_options
434 <<
"\n\tdebug = " <<
_debug
439 coutW(InputArguments) <<
"RooNDKeysPdf::setOptions() : Warning : nSigma = " <<
_nSigma <<
" < 2.0. "
440 <<
"Calculated normalization could be too large."
467 coutE(InputArguments) <<
"ERROR: RooNDKeysPdf::initialize() : The observable list is empty. "
468 <<
"Unable to begin generating the PDF." << endl;
473 coutE(InputArguments) <<
"ERROR: RooNDKeysPdf::initialize() : The input data set is empty. "
474 <<
"Unable to begin generating the PDF." << endl;
480 vector<Double_t> dummy(
_nDim,0.);
543 vector<RooRealVar*> dVars(
_nDim);
554 vector<Double_t>& point =
_dataPts[i];
563 mat(j,k) += dVars[j]->getVal() * dVars[k]->getVal() * myweight;
567 point[j] = pointV[j] = dVars[j]->getVal();
569 _x0[j] += 1. * myweight;
570 _x1[j] += point[j] * myweight ;
571 _x2[j] += point[j] * point[j] * myweight ;
572 if (
_x2[j]!=
_x2[j]) exit(3);
605 for (
Int_t j=0; j<
_nDim; j++) { sigmaRraw[j] = sqrt(sigmaRraw[j]); }
644 coutI(Contents) <<
"RooNDKeysPdf::loadDataSet(" <<
this <<
")"
645 <<
"\n Number of events in dataset: " <<
_nEvents
646 <<
"\n Weighted number of events in dataset: " <<
_nEventsW << endl;
666 vector<Double_t> dummy(
_nDim,0.);
673 vector<vector<Double_t> > mpoints(
size,dummy);
674 vector<vector<Int_t> > mjdcs(
size);
679 vector<Int_t>& mjdxK = mjdcs[0];
680 vector<Double_t>& mpointK = mpoints[0];
692 vector<Int_t>& mjdx0 = mjdcs[0];
694 if (
size==1 && mjdx0.size()==0)
continue;
698 vector<Int_t>& mjdx = mjdcs[0];
699 vector<Double_t>& mpoint = mpoints[0];
702 Int_t eMir = 1 << mjdx.size();
703 vector<vector<Double_t> > epoints(eMir,
x);
710 epoints[
l] = epoints[
l-size1];
712 vector<Double_t>& epoint = epoints[
l];
713 epoint[mjdx[
Int_t(mjdx.size()-1)-
m]] = mpoint[mjdx[
Int_t(mjdx.size()-1)-
m]];
719 epoints.erase(epoints.begin());
728 for (
Int_t j=0; j<
_nDim; j++) { pointR[j] = (epoints[
m])[j]; }
757 coutI(Contents) <<
"RooNDKeysPdf::loadWeightSet(" <<
this <<
") : Number of weighted events : " <<
_wMap.size() << endl;
784 for (
const auto& wMapItr :
_wMap) {
785 Int_t i = wMapItr.first;
794 inVarRange = inVarRange &&
kTRUE;
795 }
else { inVarRange = inVarRange &&
kFALSE; }
798 inVarRangePlusShell = inVarRangePlusShell &&
kTRUE;
799 }
else { inVarRangePlusShell = inVarRangePlusShell &&
kFALSE; }
804 bi->
bIdcs.push_back(i);
808 if (inVarRangePlusShell) {
818 if (inShell) bi->
sIdcs.push_back(i);
825 coutI(Contents) <<
"RooNDKeysPdf::calculateShell() : "
826 <<
"\n Events in shell " << bi->
sIdcs.size()
827 <<
"\n Events in box " << bi->
bIdcs.size()
828 <<
"\n Events in box and shell " << bi->
bpsIdcs.size()
846 cxcoutD(Eval) <<
"RooNDKeysPdf::calculatePreNorm() : "
860 for (
unsigned int i = 0; i <
_dataPtsR.size(); ++i) {
867 vector<TVectorD>::iterator dpRItr =
_dataPtsR.begin();
872 itrVecR.push_back(
itPair(i, dpRItr));
874 itrVecR.push_back(
itPair(i, dpRItr));
880 sort(itrVecR.begin(), itrVecR.end(), [=](
const itPair&
a,
const itPair&
b) {
881 return (*a.second)[j] < (*b.second)[j];
887 cxcoutD(Eval) <<
"RooNDKeysPdf::sortDataIndices() : Number of sorted events : " <<
_sortTVIdcs[j].size() << endl;
895 cxcoutD(Eval) <<
"RooNDKeysPdf::calculateBandWidth()" << endl;
906 cxcoutD(Eval) <<
"RooNDKeysPdf::calculateBandWidth() Using static bandwidth." << endl;
913 weight[j] =
_n * (*_sigmaR)[j];
920 cxcoutD(Eval) <<
"RooNDKeysPdf::calculateBandWidth() Using adaptive bandwidth." << endl;
922 double sqrt12 = sqrt(12.);
925 vector<Double_t> dummy(
_nDim, 0.);
928 std::vector<std::vector<Double_t>> *weights_prev(0);
929 std::vector<std::vector<Double_t>> *weights_new(0);
950 vector<Double_t> &weight = (*weights_new)[i];
952 Double_t norm = (
_n * (*_sigmaR)[j]) / sqrtSigmaAvgR;
953 weight[j] = norm *
f / sqrt12;
972 map<Int_t,Bool_t> ibMap;
983 Int_t i = ibMapItr.first;
991 const vector<Double_t> &point =
_dataPts[i];
992 const vector<Double_t> &weight = weights[
_idx[i]];
995 (*_dx)[j] =
x[j] - point[j];
1004 Double_t c = 1. / (2. * weight[j] * weight[j]);
1008 g *= exp(-
c *
r *
r);
1009 g *= 1. / (sqrt2pi * weight[j]);
1026 xRm[j] = xRp[j] =
x[j];
1034 xRm[j] -=
_nSigma * (
_n * (*_sigmaR)[j]);
1035 xRp[j] +=
_nSigma * (
_n * (*_sigmaR)[j]);
1040 vector<TVectorD> xvecRm(1,xRm);
1041 vector<TVectorD> xvecRp(1,xRp);
1043 map<Int_t,Bool_t> ibMapRT;
1048 return (*
a.second)[j] < (*
b.second)[j];
1056 if (
_nDim==1) {
for (it=lo; it!=
hi; ++it) ibMap[(*it).first] =
kTRUE; }
1057 else {
for (it=lo; it!=
hi; ++it) ibMapRT[(*it).first] =
kTRUE; }
1061 for (it=lo; it!=
hi; ++it)
1062 if (ibMapRT.find((*it).first)!=ibMapRT.end()) { ibMap[(*it).first] =
kTRUE; }
1065 if (j!=
_nDim-1) { ibMapRT = ibMap; }
1096 bi->
xVarLo[j] = var->getMin(rangeName);
1097 bi->
xVarHi[j] = var->getMax(rangeName);
1099 bi->
xVarLo[j] = var->getVal() ;
1100 bi->
xVarHi[j] = var->getVal() ;
1134 if (rangeName)
return 0 ;
1149 cxcoutD(Eval) <<
"Calling RooNDKeysPdf::analyticalIntegral(" <<
GetName() <<
") with code " << code
1150 <<
" and rangeName " << (rangeName?rangeName:
"<none>") << endl;
1162 string rangeNameStr(rangeName) ;
1175 if ((var->getMin(rangeName)-bi->
xVarLo[j]!=0) ||
1176 (var->getMax(rangeName)-bi->
xVarHi[j]!=0)) {
1183 cxcoutD(Eval) <<
"RooNDKeysPdf::analyticalIntegral() : Found new boundaries ... " << (rangeName?rangeName:
"<none>") << endl;
1188 if (!bi->
filled || newBounds) {
1201 cxcoutD(Eval) <<
"RooNDKeysPdf::analyticalIntegral() : Using mirrored normalization : " << bi->
nEventsBW << endl;
1208 if (norm<0.) norm=0.;
1213 const vector<Double_t>& weight = (*_weights)[
_idx[bi->
sIdcs[i]]];
1215 vector<Double_t> chi(
_nDim,100.);
1218 if(!doInt[j])
continue;
1221 chi[j] = (
x[j]-bi->
xVarLo[j])/weight[j];
1223 chi[j] = (bi->
xVarHi[j]-
x[j])/weight[j];
1226 prob *= (0.5 +
TMath::Erf(fabs(chi[j])/sqrt(2.))/2.);
1228 prob *= (0.5 -
TMath::Erf(fabs(chi[j])/sqrt(2.))/2.);
1234 cxcoutD(Eval) <<
"RooNDKeysPdf::analyticalIntegral() : Final normalization : " << norm <<
" " << bi->
nEventsBW << endl;
1244 std::vector<RooRealVar *> varVec;
1247 for (
const auto var : varList) {
1249 coutE(InputArguments) <<
"RooNDKeysPdf::createDatasetFromHist(" <<
GetName() <<
") WARNING: variable "
1250 << var->
GetName() <<
" is not of type RooRealVar. Skip." << endl;
1253 varsAndWeightSet.
add(*var);
1254 varVec.push_back(
static_cast<RooRealVar *
>(var));
1258 RooRealVar weight(
"weight",
"event weight", 0);
1259 varsAndWeightSet.
add(weight);
1262 unsigned int histndim(0);
1263 std::string classname = hist.
ClassName();
1264 if (classname.find(
"TH1") == 0) {
1266 }
else if (classname.find(
"TH2") == 0) {
1268 }
else if (classname.find(
"TH3") == 0) {
1271 assert(histndim == varVec.size());
1273 if (histndim > 3 || histndim <= 0) {
1274 coutE(InputArguments) <<
"RooNDKeysPdf::createDatasetFromHist(" <<
GetName()
1275 <<
") ERROR: input histogram dimension not between [1-3]: " << histndim << endl;
1283 for (
int i = 1; i <= hist.
GetXaxis()->GetNbins(); ++i) {
1287 varVec[0]->setVal(xval);
1289 if (varVec.size() == 1) {
1292 dataFromHist->
add(varsAndWeightSet, fval);
1295 for (
int j = 1; j <= hist.
GetYaxis()->GetNbins(); ++j) {
1297 varVec[1]->setVal(yval);
1299 if (varVec.size() == 2) {
1302 dataFromHist->
add(varsAndWeightSet, fval);
1305 for (
int k = 1; k <= hist.
GetZaxis()->GetNbins(); ++k) {
1307 varVec[2]->setVal(zval);
1311 dataFromHist->
add(varsAndWeightSet, fval);
1318 return dataFromHist;
1330 cxcoutD(Eval) <<
"RooNDKeysPdf::getWeights() Return evaluated weights." << endl;
1333 const vector<Double_t>&
x =
_dataPts[i];
1338 const vector<Double_t>& weight = (*_weights)[i];
1339 mref(i,
_nDim) = weight[k];
1357 covMatRho(j, k) = (*_covMat)(j, k) *
_rho[j] *
_rho[k];
1364 (*_sigmaR)[j] = sqrt((*
_sigmaR)[j]);
typedef void(GLAPIENTRYP _GLUfuncptr)(void)
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
std::vector< itPair > itVec
std::pair< Int_t, VecTVecDouble::iterator > itPair
TMatrixTSym< Double_t > TMatrixDSym
TMatrixT< Double_t > TMatrixD
TVectorT< Double_t > TVectorD
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Storage_t::size_type size() const
const char * GetName() const
Returns name of object.
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual Int_t numEntries() const
Return number of entries in dataset, i.e., count unweighted entries.
const RooArgSet * nset() const
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
Double_t 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.
RooChangeTracker is a meta object that tracks value changes in a given set of RooAbsArgs by registeri...
Bool_t hasChanged(Bool_t clearState)
Returns true if state has changed since last call with clearState=kTRUE.
RooDataSet is a container class to hold unbinned data.
virtual const RooArgSet * get(Int_t index) const override
Return RooArgSet with coordinates of event 'index'.
virtual void add(const RooArgSet &row, Double_t weight=1.0, Double_t weightError=0) override
Add a data point, with its coordinates specified in the 'data' argset, to the data set.
virtual Double_t weight() const override
Return event weight of current event.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Reimplementation of standard RooArgList::add()
Generic N-dimensional implementation of a kernel estimation p.d.f.
std::vector< Int_t > _sIdcs
void initialize()
initialization
std::vector< Double_t > _xVarLo
void calculatePreNorm(BoxInfo *bi) const
bi->nEventsBMSW=0.; bi->nEventsBW=0.;
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
std::vector< Double_t > _xDatLo3s
std::vector< Double_t > _x0
std::vector< Double_t > _xVarLoP3s
std::vector< Double_t > _xDatHi3s
std::vector< Int_t > _bIdcs
std::vector< TVectorD > _dataPtsR
void loopRange(std::vector< Double_t > &x, std::map< Int_t, Bool_t > &ibMap) const
determine closest points to x, to loop over in evaluate()
std::vector< Double_t > _xVarHi
void createPdf(Bool_t firstCall=kTRUE)
evaluation order of constructor.
void calculateShell(BoxInfo *bi) const
determine points in +/- nSigma shell around the box determined by the variable ranges.
void calculateBandWidth()
std::vector< Double_t > _x1
std::vector< Double_t > _mean
std::vector< Int_t > _bmsIdcs
std::vector< Double_t > _rho
std::vector< Double_t > _xDatHi
std::vector< Double_t > _x
std::vector< itVec > _sortTVIdcs
Weights to be used. Points either to _weights0 or _weights1.
void boxInfoInit(BoxInfo *bi, const char *rangeName, Int_t code) const
std::map< std::pair< std::string, int >, BoxInfo * > _rangeBoxInfo
std::vector< Double_t > _xVarLoM3s
std::vector< Int_t > _idx
std::vector< Double_t > _xDatLo
void sortDataIndices(BoxInfo *bi=0)
sort entries, as needed for loopRange()
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported.
void loadDataSet(Bool_t firstCall)
copy the dataset and calculate some useful variables
TString _options
do not persist
std::vector< std::vector< Double_t > > * _weights
TMatrixD getWeights(const int &k) const
Return evaluated weights.
std::vector< Double_t > _xVarHiM3s
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
std::map< Int_t, Bool_t > _ibNoSort
std::vector< std::vector< Double_t > > _weights0
void mirrorDataSet()
determine mirror dataset.
void setOptions()
set the configuration
std::vector< std::vector< Double_t > > _dataPts
RooDataSet * createDatasetFromHist(const RooArgList &varList, const TH1 &hist) const
std::vector< Double_t > _sigma
std::vector< std::string > _varName
Double_t gauss(std::vector< Double_t > &x, std::vector< std::vector< Double_t > > &weights) const
loop over all closest point to x, as determined by loopRange()
void checkInitWeights() const
std::map< Int_t, Double_t > _wMap
std::vector< std::vector< Double_t > > _weights1
std::vector< Double_t > _x2
std::vector< Double_t > _xVarHiP3s
std::map< Int_t, Bool_t > _bpsIdcs
RooChangeTracker * _tracker
RooRealVar represents a variable that can be changed from the outside.
virtual void setVal(Double_t value)
Set value of variable to 'value'.
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
TH1 is the base class of all histogram classes in ROOT.
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
const TVectorD & GetEigenValues() const
const TMatrixD & GetEigenVectors() const
void Print(Option_t *name="") const
Print the matrix as a table of elements.
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
TMatrixT< Element > & T()
virtual const char * GetName() const
Returns name of object.
virtual const char * ClassName() const
Returns name of class to which the object belongs.
void ToLower()
Change string to lower-case.
const char * Data() const
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
TVectorT< Element > & Zero()
Set vector elements to zero.
void Print(Option_t *option="") const
Print the vector as a list of elements.
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Double_t Erf(Double_t x)
Computation of the error function erf(x).
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
constexpr Double_t TwoPi()
std::vector< Double_t > xVarHiP3s
std::vector< Double_t > xVarHi
std::map< Int_t, Bool_t > bpsIdcs
std::vector< Int_t > bIdcs
std::vector< Double_t > xVarLoM3s
std::vector< Double_t > xVarHiM3s
std::vector< Double_t > xVarLoP3s
std::vector< Double_t > xVarLo
std::vector< Int_t > sIdcs
std::vector< Int_t > bmsIdcs