38 : _curNormSet(0), _M(0), _MSqr(0), _setting(
RooMomentMorphND::Linear), _useHorizMorph(true)
40 _parItr = _parList.createIterator();
41 _obsItr = _obsList.createIterator();
49 :
RooAbsPdf(
name, title), _cacheMgr(this, 10,
kTRUE,
kTRUE), _parList(
"parList",
"List of morph parameters", this),
50 _obsList(
"obsList",
"List of observables", this), _referenceGrid(referenceGrid),
51 _pdfList(
"pdfList",
"List of pdfs", this), _setting(setting), _useHorizMorph(true)
70 :
RooAbsPdf(
name, title), _cacheMgr(this, 10,
kTRUE,
kTRUE), _parList(
"parList",
"List of morph parameters", this),
71 _obsList(
"obsList",
"List of observables", this), _pdfList(
"pdfList",
"List of pdfs", this), _setting(setting),
78 for (
int i = 0; i < mrefpoints.
GetNrows(); ++i) {
80 if (mrefpoints[i] == grid.
array()[j]) {
106 :
RooAbsPdf(
name, title), _cacheMgr(this, 10,
kTRUE,
kTRUE), _parList(
"parList",
"List of morph parameters", this),
107 _obsList(
"obsList",
"List of observables", this), _pdfList(
"pdfList",
"List of pdfs", this), _setting(setting),
114 for (
int i = 0; (mref =
dynamic_cast<RooAbsReal *
>(mrefItr->
Next())); ++i) {
116 coutE(InputArguments) <<
"RooMomentMorphND::ctor(" <<
GetName() <<
") ERROR: mref " << mref->
GetName()
117 <<
" is not of type RooAbsReal" << endl;
118 throw string(
"RooMomentMorphND::ctor() ERROR mref is not of type RooAbsReal");
121 coutW(InputArguments) <<
"RooMomentMorphND::ctor(" <<
GetName() <<
") WARNING mref point " << i
122 <<
" is not a constant, taking a snapshot of its value" << endl;
124 mrefpoints[i] = mref->
getVal();
131 for (
int i = 0; i < mrefpoints.
GetNrows(); ++i) {
133 if (mrefpoints[i] == grid.
array()[j]) {
158 :
RooAbsPdf(other,
name), _cacheMgr(other._cacheMgr, this), _curNormSet(0),
159 _parList(
"parList", this, other._parList), _obsList(
"obsList", this, other._obsList),
160 _referenceGrid(other._referenceGrid), _pdfList(
"pdfList", this, other._pdfList), _M(0), _MSqr(0),
161 _setting(other._setting), _useHorizMorph(other._useHorizMorph)
194 coutE(InputArguments) <<
"RooMomentMorphND::ctor(" <<
GetName() <<
") ERROR: parameter " << par->
GetName()
195 <<
" is not of type RooAbsReal" << endl;
196 throw string(
"RooMomentMorphND::initializeParameters() ERROR parameter is not of type RooAbsReal");
212 coutE(InputArguments) <<
"RooMomentMorphND::ctor(" <<
GetName() <<
") ERROR: variable " << var->
GetName()
213 <<
" is not of type RooAbsReal" << endl;
214 throw string(
"RooMomentMorphND::initializeObservables() ERROR variable is not of type RooAbsReal");
227 typename vector<T>::const_iterator begin;
228 typename vector<T>::const_iterator end;
229 typename vector<T>::const_iterator me;
235 vector<Digits<T>> vd;
237 for (
typename vector<vector<T>>::const_iterator it = in.begin(); it != in.end(); ++it) {
238 Digits<T>
d = {(*it).begin(), (*it).end(), (*it).begin()};
244 for (
typename vector<Digits<T>>::const_iterator it = vd.begin(); it != vd.end(); ++it) {
245 result.push_back(*(it->me));
247 out.push_back(result);
249 for (
typename vector<Digits<T>>::iterator it = vd.begin();;) {
251 if (it->me == it->end) {
252 if (it + 1 == vd.end()) {
292 coutE(InputArguments) <<
"RooMomentMorphND::initialize(" <<
GetName() <<
") ERROR: nPar != nDim"
293 <<
": " << nPar <<
" !=" << nDim << endl;
298 coutE(InputArguments) <<
"RooMomentMorphND::initialize(" <<
GetName() <<
") ERROR: nPdf != nRef"
299 <<
": " << nPdf <<
" !=" << nRef << endl;
309 vector<vector<double>> dm(nPdf);
310 for (
int k = 0; k < nPdf; ++k) {
312 for (
int idim = 0; idim < nPar; idim++) {
314 dm2.push_back(delta);
319 vector<vector<int>> powers;
320 for (
int idim = 0; idim < nPar; idim++) {
325 powers.push_back(xtmp);
328 vector<vector<int>>
output;
330 int nCombs =
output.size();
332 for (
int k = 0; k < nPdf; ++k) {
334 for (
int i = 0; i < nCombs; i++) {
336 for (
int ix = 0; ix < nPar; ix++) {
356 : _grid(other._grid), _pdfList(other._pdfList), _pdfMap(other._pdfMap), _nref(other._nref)
368 vector<int> thisBoundaries;
369 vector<double> thisBoundaryCoordinates;
370 thisBoundaries.push_back(bin_x);
371 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
374 _nref.push_back(thisBoundaryCoordinates);
380 vector<int> thisBoundaries;
381 vector<double> thisBoundaryCoordinates;
382 thisBoundaries.push_back(bin_x);
383 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
384 thisBoundaries.push_back(bin_y);
385 thisBoundaryCoordinates.push_back(_grid[1]->array()[bin_y]);
388 _nref.push_back(thisBoundaryCoordinates);
394 vector<int> thisBoundaries;
395 vector<double> thisBoundaryCoordinates;
396 thisBoundaries.push_back(bin_x);
397 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
398 thisBoundaries.push_back(bin_y);
399 thisBoundaryCoordinates.push_back(_grid[1]->array()[bin_y]);
400 thisBoundaries.push_back(bin_z);
401 thisBoundaryCoordinates.push_back(_grid[2]->array()[bin_z]);
404 _nref.push_back(thisBoundaryCoordinates);
410 vector<double> thisBoundaryCoordinates;
411 int nBins = bins.size();
412 for (
int i = 0; i < nBins; i++) {
413 thisBoundaryCoordinates.push_back(_grid[i]->array()[bins[i]]);
417 _nref.push_back(thisBoundaryCoordinates);
434 vector<RooAbsReal *> meanrv(nPdf * nObs, null);
435 vector<RooAbsReal *> sigmarv(nPdf * nObs, null);
436 vector<RooAbsReal *> myrms(nObs, null);
437 vector<RooAbsReal *> mypos(nObs, null);
438 vector<RooAbsReal *> slope(nPdf * nObs, null);
439 vector<RooAbsReal *> offsetrv(nPdf * nObs, null);
440 vector<RooAbsReal *> transVar(nPdf * nObs, null);
441 vector<RooAbsReal *> transPdf(nPdf, null);
451 for (
int i = 0; i < 3 * nPdf; ++i) {
452 string fracName =
Form(
"frac_%d", i);
458 else if (i < 2 * nPdf)
471 for (
int i = 0; i < nPdf; ++i) {
472 for (
int j = 0; j < nObs; ++j) {
479 sigmarv[
sij(i, j)] = mom;
480 meanrv[
sij(i, j)] = mom->
mean();
482 ownedComps.
add(*sigmarv[
sij(i, j)]);
487 for (
int j = 0; j < nObs; ++j) {
490 for (
int i = 0; i < nPdf; ++i) {
491 meanList.
add(*meanrv[
sij(i, j)]);
492 rmsList.
add(*sigmarv[
sij(i, j)]);
494 string myrmsName =
Form(
"%s_rms_%d",
GetName(), j);
495 string myposName =
Form(
"%s_pos_%d",
GetName(), j);
496 mypos[j] =
new RooAddition(myposName.c_str(), myposName.c_str(), meanList, coefList2);
497 myrms[j] =
new RooAddition(myrmsName.c_str(), myrmsName.c_str(), rmsList, coefList3);
506 for (
int i = 0; i < nPdf; ++i) {
511 string pdfName =
Form(
"pdf_%d", i);
514 for (
int j = 0; j < nObs; ++j) {
516 string slopeName =
Form(
"%s_slope_%d_%d",
GetName(), i, j);
517 string offsetName =
Form(
"%s_offset_%d_%d",
GetName(), i, j);
527 string transVarName =
Form(
"%s_transVar_%d_%d",
GetName(), i, j);
528 transVar[
sij(i, j)] =
new RooLinearVar(transVarName.c_str(), transVarName.c_str(), *var, *slope[
sij(i, j)],
529 *offsetrv[
sij(i, j)]);
535 ownedComps.
add(*transVar[
sij(i, j)]);
539 transPdfList.
add(*transPdf[i]);
540 ownedComps.
add(*transPdf[i]);
544 theSumPdf =
new RooAddPdf(sumpdfName.c_str(), sumpdfName.c_str(), transPdfList, coefList);
546 theSumPdf =
new RooAddPdf(sumpdfName.c_str(), sumpdfName.c_str(),
_pdfList, coefList);
555 string trackerName =
Form(
"%s_frac_tracker",
GetName());
559 cache =
new CacheElem(*theSumPdf, *tracker, fracl);
626template <
typename Iterator>
629 if ((
first == last) || (
first == k) || (last == k)) {
632 Iterator itr1 =
first;
633 Iterator itr2 = last;
642 while (
first != itr1) {
643 if (*--itr1 < *itr2) {
645 while (!(*itr1 < *j)) ++j;
650 rotate(itr1, j, last);
655 rotate(k, itr2, last);
659 rotate(
first, k, last);
675 for (
int idim = 0; idim < nPar; idim++) {
677 dm2.push_back(delta);
680 vector<vector<int>> powers;
681 for (
int idim = 0; idim < nPar; idim++) {
686 powers.push_back(xtmp);
689 vector<vector<int>>
output;
691 int nCombs =
output.size();
693 vector<double> deltavec(nPdf, 1.0);
696 for (
int i = 0; i < nCombs; i++) {
698 for (
int ix = 0; ix < nPar; ix++) {
702 deltavec[nperm] = tmpDm;
706 double sumposfrac = 0.0;
707 for (
int i = 0; i < nPdf; ++i) {
710 for (
int j = 0; j < nPdf; ++j) {
711 ffrac += (*self.
_M)(j, i) * deltavec[j] * fracNonLinear;
724 ((
RooRealVar *)frac(nPdf + i))->setVal(ffrac);
725 ((
RooRealVar *)frac(2 * nPdf + i))->setVal(ffrac);
728 cout <<
"NonLinear fraction " << ffrac << endl;
730 frac(nPdf + i)->Print();
731 frac(2 * nPdf + i)->Print();
736 for (
int i = 0; i < nPdf; ++i) {
750 for (
int i = 0; i < nPdf; ++i) {
753 ((
RooRealVar *)frac(nPdf + i))->setVal(initval);
754 ((
RooRealVar *)frac(2 * nPdf + i))->setVal(initval);
759 for (
int j = 0; j < nPar; j++) {
761 mtmp.push_back(
m->getVal());
767 vector<double> deltavec(depth, 1.0);
772 for (
int ix = 0; ix < nPar; ix++) {
776 for (
int iperm = 1; iperm <= nPar; ++iperm) {
778 double dtmp = mtmp[xtmp[0]] - self.
_squareVec[0][xtmp[0]];
779 for (
int itmp = 1; itmp < iperm; ++itmp) {
780 dtmp *= mtmp[xtmp[itmp]] - self.
_squareVec[0][xtmp[itmp]];
782 deltavec[nperm + 1] = dtmp;
787 Double_t origFrac1(0.), origFrac2(0.);
788 for (
int i = 0; i < depth; ++i) {
790 for (
int j = 0; j < depth; ++j) {
791 ffrac += (*self.
_MSqr)(j, i) * deltavec[j] * fracLinear;
807 cout <<
"Linear fraction " << ffrac << endl;
835 vector<vector<double>> boundaries(nPar);
836 for (
int idim = 0; idim < nPar; idim++) {
840 boundaries[idim].push_back(lo);
841 boundaries[idim].push_back(
hi);
844 vector<vector<double>>
output;
848 for (
int isq = 0; isq < depth; isq++) {
849 for (
int iref = 0; iref < nRef; iref++) {
872 for (
int ix = 0; ix < nPar; ix++) {
876 for (
int k = 0; k < depth; ++k) {
882 for (
int iperm = 1; iperm <= nPar; ++iperm) {
884 double dtmp =
_squareVec[k][xtmp[0]] - squareBase[xtmp[0]];
885 for (
int itmp = 1; itmp < iperm; ++itmp) {
886 dtmp *=
_squareVec[k][xtmp[itmp]] - squareBase[xtmp[itmp]];
888 M(k, nperm + 1) = dtmp;
908 cout <<
"Currently BinIntegrator only knows how to deal with 1-d " << endl;
bool next_combination(const Iterator first, Iterator k, const Iterator last)
void cartesian_product(vector< vector< T > > &out, vector< vector< T > > &in)
TMatrixT< Double_t > TMatrixD
char * Form(const char *fmt,...)
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Bool_t addOwnedComponents(const RooArgSet &comps)
Take ownership of the contents of 'comps'.
void addServerList(RooAbsCollection &serverList, Bool_t valueProp=kTRUE, Bool_t shapeProp=kFALSE)
Register a list of RooAbsArg as servers to us by calling addServer() for each arg in the list.
void setLocalNoDirtyInhibit(Bool_t flag) const
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Bool_t setRealValue(const char *name, Double_t newVal=0, Bool_t verbose=kFALSE)
Set value of a RooAbsRealLValye stored in set with given name to newVal No error messages are printed...
RooAbsArg * first() const
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsMoment represents the first, second, or third order derivative of any RooAbsReal as calculated ...
const RooArgSet * nset() const
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
RooNumIntConfig * specialIntegratorConfig() const
Returns the specialized integrator configuration for this RooAbsReal.
RooAddPdf is an efficient implementation of a sum of PDFs of the form.
RooAddition calculates the sum of a set of RooAbsReal terms, or when constructed with two sets,...
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.
Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Add element to non-owning set.
Class RooBinning is an implements RooAbsBinning in terms of an array of boundary values,...
virtual Int_t numBoundaries() const
virtual Double_t * array() const
Return array of boundary values.
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=0, const TNamed *isetRangeName=0)
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=0)
virtual Bool_t setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
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.
RooConstVar represent a constant real-valued object.
RooCustomizer is a factory class to produce clones of a prototype composite PDF object with the same ...
void replaceArg(const RooAbsArg &orig, const RooAbsArg &subst)
Replace any occurence of arg 'orig' with arg 'subst'.
RooAbsArg * build(const char *masterCatState, Bool_t verbose=kFALSE)
Build a clone of the prototype executing all registered 'replace' rules and 'split' rules for the mas...
RooLinearVar is the most general form of a derived real-valued object that can be used by RooRealInte...
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Reimplementation of standard RooArgList::add()
RooChangeTracker * _tracker
virtual RooArgList containedArgs(Action)
void calculateFractions(const RooMomentMorphND &self, Bool_t verbose=kTRUE) const
void addBinning(const RooAbsBinning &binning)
std::vector< std::vector< double > > _nref
std::vector< RooAbsBinning * > _grid
void addPdf(const RooAbsPdf &pdf, int bin_x)
std::vector< int > _nnuis
virtual Double_t getVal(const RooArgSet *set=0) const
std::vector< int > _squareIdx
Grid _referenceGrid
Do not persist.
CacheElem * getCache(const RooArgSet *nset) const
RooObjCacheManager _cacheMgr
Bool_t setBinIntegrator(RooArgSet &allVars)
std::vector< std::vector< double > > _squareVec
void initializeObservables(const RooArgList &obsList)
int sij(const int &i, const int &j) const
void findShape(const std::vector< double > &x) const
Double_t evaluate() const
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
TIterator * _obsItr
Do not persist.
virtual ~RooMomentMorphND()
RooAbsPdf * sumPdf(const RooArgSet *nset)
void initializeParameters(const RooArgList &parList)
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
RooRealVar represents a variable that can be changed from the outside.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Overloaded RooArgSet::add() method inserts 'var' into set and registers 'var' as server to owner with...
Iterator abstract base class.
virtual TObject * Next()=0
TMatrixT< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant.
virtual const char * GetName() const
Returns name of object.
Element * GetMatrixArray()
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
static void output(int code)