89 _varItr = _varList.createIterator() ;
91 TIterator* varItr = varList.createIterator() ;
94 if (!dynamic_cast<RooAbsReal*>(var)) {
96 <<
" is not of type RooAbsReal" << endl ;
100 _varName.push_back(var->
GetName());
114 _varList(
"varList",
"List of variables",this),
119 _weights(&_weights0),
127 if (!dynamic_cast<RooAbsReal*>(var)) {
129 <<
" is not of type RooAbsReal" << endl ;
139 coutE(
InputArguments) <<
"ERROR: RhhNDKeysPdf::RhhNDKeysPdf() : The vector-size of rho is different from that of varList."
140 <<
"Unable to create the PDF." << endl;
162 _varList(
"varList",
"List of variables",this),
167 _weights(&_weights0),
177 coutW(
InputArguments) <<
"RooNDKeysPdf::RooNDKeysPdf() : Warning : asymmetric mirror(s) no longer supported." << endl;
193 _varList(
"varList",
"List of variables",this),
198 _weights(&_weights0),
217 _varList(
"varList",this,other._varList),
219 _options(other._options),
220 _widthFactor(other._widthFactor),
221 _nSigma(other._nSigma),
222 _weights(&_weights0),
223 _rotate(other._rotate)
379 <<
"\n\tdebug = " <<
_debug
385 <<
"Calculated normalization could be too large."
409 coutE(
InputArguments) <<
"ERROR: RooNDKeysPdf::initialize() : The observable list is empty. "
410 <<
"Unable to begin generating the PDF." << endl;
415 coutE(
InputArguments) <<
"ERROR: RooNDKeysPdf::initialize() : The input data set is empty. "
416 <<
"Unable to begin generating the PDF." << endl;
491 vector<RooRealVar*> dVars(
_nDim);
502 vector<Double_t>& point =
_dataPts[i];
511 mat(j,k) += dVars[j]->getVal() * dVars[k]->getVal() * myweight;
515 point[j] = pointV[j] = dVars[j]->getVal();
517 _x0[j] += 1. * myweight;
518 _x1[j] += point[j] * myweight ;
519 _x2[j] += point[j] * point[j] * myweight ;
520 if (
_x2[j]!=
_x2[j]) exit(3);
543 covMatRho(j,k) = (mat(j,k)/
_x0[j] - _mean[j]*_mean[k]) *
_rho[j] *
_rho[k];
585 for (
Int_t j=0; j<
_nDim; j++) { sigmaRraw[j] =
sqrt(sigmaRraw[j]); }
598 coutI(
Contents) <<
"RooNDKeysPdf::loadDataSet(" <<
this <<
")"
599 <<
"\n Number of events in dataset: " << _nEvents
600 <<
"\n Weighted number of events in dataset: " <<
_nEventsW
624 vector<Double_t>
dummy(_nDim,0.);
631 vector<vector<Double_t> > mpoints(size,dummy);
632 vector<vector<Int_t> > mjdcs(size);
637 vector<Int_t>& mjdxK = mjdcs[0];
638 vector<Double_t>& mpointK = mpoints[0];
642 mpointK[j] = 2.*
_xDatLo[j]-x[j];
645 mpointK[j] = 2.*
_xDatHi[j]-x[j];
650 vector<Int_t>& mjdx0 = mjdcs[0];
652 if (size==1 && mjdx0.size()==0)
continue;
656 vector<Int_t>& mjdx = mjdcs[0];
657 vector<Double_t>& mpoint = mpoints[0];
660 Int_t eMir = 1 << mjdx.size();
661 vector<vector<Double_t> > epoints(eMir,x);
668 epoints[
l] = epoints[
l-size1];
670 vector<Double_t>& epoint = epoints[
l];
671 epoint[mjdx[
Int_t(mjdx.size()-1)-
m]] = mpoint[mjdx[
Int_t(mjdx.size()-1)-
m]];
677 epoints.erase(epoints.begin());
686 for (
Int_t j=0; j<
_nDim; j++) { pointR[j] = (epoints[
m])[j]; }
715 coutI(
Contents) <<
"RooNDKeysPdf::loadWeightSet(" <<
this <<
") : Number of weighted events : " <<
_wMap.size() << endl;
743 map<Int_t,Double_t>::iterator wMapItr =
_wMap.begin();
746 for (; wMapItr!=
_wMap.end(); ++wMapItr) {
747 Int_t i = (*wMapItr).first;
756 inVarRange = inVarRange &&
kTRUE;
757 }
else { inVarRange = inVarRange &&
kFALSE; }
760 inVarRangePlusShell = inVarRangePlusShell &&
kTRUE;
761 }
else { inVarRangePlusShell = inVarRangePlusShell &&
kFALSE; }
766 bi->
bIdcs.push_back(i);
770 if (inVarRangePlusShell) {
780 if (inShell) bi->
sIdcs.push_back(i);
788 <<
"\n Events in shell " << bi->
sIdcs.size()
789 <<
"\n Events in box " << bi->
bIdcs.size()
790 <<
"\n Events in box and shell " << bi->
bpsIdcs.size()
810 cxcoutD(
Eval) <<
"RooNDKeysPdf::calculatePreNorm() : "
824 vector<TVectorD>::iterator dpRItr =
_dataPtsR.begin();
829 itrVecR.push_back(
itPair(i,dpRItr));
830 }
else itrVecR.push_back(
itPair(i,dpRItr));
840 cxcoutD(
Eval) <<
"RooNDKeysPdf::sortDataIndices() : Number of sorted events : " <<
_sortTVIdcs[j].size() << endl;
850 cxcoutD(
Eval) <<
"RooNDKeysPdf::calculateBandWidth()" << endl;
856 cxcoutD(
Eval) <<
"RooNDKeysPdf::calculateBandWidth() Using static bandwidth." << endl;
862 weight[j] =
_rho[j] *
_n * (*_sigmaR)[j] ;
870 cxcoutD(
Eval) <<
"RooNDKeysPdf::calculateBandWidth() Using adaptive bandwidth." << endl;
872 double sqrt12=
sqrt(12.);
887 weight[j] = norm * f / sqrt12 ;
904 map<Int_t,Bool_t> ibMap;
910 map<Int_t,Bool_t>::iterator ibMapItr = ibMap.begin();
912 for (; ibMapItr!=ibMap.end(); ++ibMapItr) {
913 Int_t i = (*ibMapItr).first;
919 const vector<Double_t>& point =
_dataPts[i];
920 const vector<Double_t>& weight = weights[
_idx[i]];
923 (*_dx)[j] = x[j]-point[j];
932 Double_t c = 1./(2.*weight[j]*weight[j]);
937 z += (g*
_wMap[_idx[i]]);
952 for (
Int_t j=0; j<
_nDim; j++) { xRm[j] = xRp[j] = x[j]; }
963 vector<TVectorD> xvecRm(1,xRm);
964 vector<TVectorD> xvecRp(1,xRp);
966 map<Int_t,Bool_t> ibMapRT;
974 itVec::iterator it=lo;
976 if (_nDim==1) {
for (it=lo; it!=
hi; ++it) ibMap[(*it).first] =
kTRUE; }
977 else {
for (it=lo; it!=
hi; ++it) ibMapRT[(*it).first] =
kTRUE; }
981 for (it=lo; it!=
hi; ++it)
982 if (ibMapRT.find((*it).first)!=ibMapRT.end()) { ibMap[(*it).first] =
kTRUE; }
985 if (j!=_nDim-1) { ibMapRT = ibMap; }
1056 if (rangeName)
return 0 ;
1071 cxcoutD(
Eval) <<
"Calling RooNDKeysPdf::analyticalIntegral(" <<
GetName() <<
") with code " << code
1072 <<
" and rangeName " << (rangeName?rangeName:
"<none>") << endl;
1084 string rangeNameStr(rangeName) ;
1106 cxcoutD(
Eval) <<
"RooNDKeysPdf::analyticalIntegral() : Found new boundaries ... " << (rangeName?rangeName:
"<none>") << endl;
1111 if (!bi->
filled || newBounds) {
1124 cxcoutD(
Eval) <<
"RooNDKeysPdf::analyticalIntegral() : Using mirrored normalization : " << bi->
nEventsBW << endl;
1131 if (norm<0.) norm=0.;
1136 const vector<Double_t>& weight = (*_weights)[
_idx[bi->
sIdcs[i]]];
1138 vector<Double_t> chi(
_nDim,100.);
1141 if(!doInt[j])
continue;
1144 chi[j] = (x[j]-bi->
xVarLo[j])/weight[j];
1146 chi[j] = (bi->
xVarHi[j]-x[j])/weight[j];
1157 cxcoutD(
Eval) <<
"RooNDKeysPdf::analyticalIntegral() : Final normalization : " << norm <<
" " << bi->
nEventsBW << endl;
const TVectorD & GetEigenValues() const
std::vector< itPair > itVec
const RooArgSet * nset() const
std::vector< Double_t > xVarHi
std::vector< Double_t > _xVarLo
virtual TMatrixTBase< Element > & Zero()
Set matrix elements to zero.
void variables(TString fin="TMVA.root", TString dirName="InputVariables_Id", TString title="TMVA Input Variables", Bool_t isRegression=kFALSE, Bool_t useTMVAStyle=kTRUE)
std::vector< std::string > _varName
std::vector< Int_t > bIdcs
void boxInfoInit(BoxInfo *bi, const char *rangeName, Int_t code) const
std::vector< Int_t > _bIdcs
virtual Double_t getMin(const char *name=0) const
std::vector< Double_t > _xDatLo
TMatrixT< Element > & T()
void ToLower()
Change string to lower-case.
Double_t evaluate() const
do not persist
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
std::vector< Double_t > _x1
Iterator abstract base class.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Int_t getAnalyticalIntegral(RooArgSet &allVars, RooArgSet &analVars, const char *rangeName=0) const
Interface function getAnalyticalIntergral advertises the analytical integrals that are supported...
std::vector< Double_t > _xDatHi
std::vector< Double_t > _xVarLoM3s
std::vector< std::vector< Double_t > > _weights1
std::vector< Double_t > _xVarHiM3s
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::map< Int_t, Double_t > _wMap
std::map< std::string, std::string >::const_iterator iter
void calculateBandWidth() const
TIterator * createIterator(Bool_t dir=kIterForward) const
std::map< std::pair< std::string, int >, BoxInfo * > _rangeBoxInfo
Double_t getVal(const RooArgSet *set=0) const
std::vector< Double_t > _xDatLo3s
TVectorT< Double_t > TVectorD
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
std::vector< std::vector< Double_t > > * _weights
TMatrixT< Double_t > TMatrixD
void sortDataIndices(BoxInfo *bi=0) const
sort entries, as needed for loopRange()
void Print(Option_t *name="") const
Print the matrix as a table of elements.
void Print(Option_t *option="") const
Print the vector as a list of elements.
std::vector< Double_t > xVarHiP3s
void calculatePreNorm(BoxInfo *bi) const
bi->nEventsBMSW=0.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
std::vector< Double_t > _x
RooAbsArg * find(const char *name) const
Find object with given name in list.
virtual Double_t weight() const
Return event weight of current event.
std::vector< itVec > _sortTVIdcs
Double_t Erf(Double_t x)
Computation of the error function erf(x).
void setOptions() const
set the configuration
virtual Int_t numEntries() const
void createPdf(Bool_t firstCall=kTRUE) const
evaluation order of constructor.
TVectorT< Element > & Zero()
Set vector elements to zero.
virtual const char * GetName() const
Returns name of object.
std::vector< Double_t > xVarHiM3s
RooNDKeysPdf(const char *name, const char *title, const RooArgList &varList, RooDataSet &data, TString options="a", Double_t rho=1, Double_t nSigma=3, Bool_t rotate=kTRUE)
std::vector< Double_t > _x0
std::vector< Double_t > _mean
std::vector< Double_t > _xVarHiP3s
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
const TMatrixD & GetEigenVectors() const
std::map< Int_t, Bool_t > _bpsIdcs
std::vector< Int_t > sIdcs
void mirrorDataSet() const
determine mirror dataset.
TMatrixTSym< Double_t > TMatrixDSym
static RooMathCoreReg dummy
std::vector< Double_t > xVarLoP3s
std::vector< Double_t > _xVarHi
std::vector< Int_t > _idx
std::pair< Int_t, VecTVecDouble::iterator > itPair
std::vector< Double_t > _xDatHi3s
std::vector< Int_t > _sIdcs
std::map< Int_t, Bool_t > bpsIdcs
virtual Double_t getMax(const char *name=0) const
std::vector< Int_t > bmsIdcs
std::vector< Double_t > _rho
void initialize() const
initialization
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Double_t analyticalIntegral(Int_t code, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral. ...
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
virtual TObject * Next()=0
void calculateShell(BoxInfo *bi) const
determine points in +/- nSigma shell around the box determined by the variable ranges.
std::vector< TVectorD > _dataPtsR
std::vector< std::vector< Double_t > > _dataPts
std::vector< std::vector< Double_t > > _weights0
void loadWeightSet() const
float type_of_call hi(const int &, const int &)
void loadDataSet(Bool_t firstCall) const
copy the dataset and calculate some useful variables
Bool_t matchArgs(const RooArgSet &allDeps, RooArgSet &numDeps, const RooArgProxy &a) const
Utility function for use in getAnalyticalIntegral().
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
std::vector< Double_t > _xVarLoP3s
double norm(double *x, double *p)
std::vector< Double_t > _sigma
virtual const RooArgSet * get(Int_t index) const
Return RooArgSet with coordinates of event 'index'.
std::vector< Double_t > xVarLo
std::vector< Int_t > _bmsIdcs
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()
std::vector< Double_t > xVarLoM3s
std::vector< Double_t > _x2
ClassImp(RooNDKeysPdf) RooNDKeysPdf
Construct N-dimensional kernel estimation p.d.f.