43using std::cout, std::endl, std::string, std::vector;
58 _cacheMgr(this, 10, true, true),
59 _parList(
"parList",
"List of morph parameters", this),
60 _obsList(
"obsList",
"List of observables", this),
61 _referenceGrid(referenceGrid),
62 _pdfList(
"pdfList",
"List of pdfs", this),
84 _cacheMgr(this, 10, true, true),
85 _parList(
"parList",
"List of morph parameters", this),
86 _obsList(
"obsList",
"List of observables", this),
87 _pdfList(
"pdfList",
"List of pdfs", this),
95 for (
int i = 0; i < mrefpoints.
GetNrows(); ++i) {
97 if (mrefpoints[i] == grid.
array()[j]) {
124 _cacheMgr(this, 10, true, true),
125 _parList(
"parList",
"List of morph parameters", this),
126 _obsList(
"obsList",
"List of observables", this),
127 _pdfList(
"pdfList",
"List of pdfs", this),
134 for (
auto *mref : mrefList) {
136 coutE(InputArguments) <<
"RooMomentMorphFuncND::ctor(" <<
GetName() <<
") ERROR: mref " << mref->GetName()
137 <<
" is not of type RooAbsReal" << endl;
138 throw string(
"RooMomentMorphFuncND::ctor() ERROR mref is not of type RooAbsReal");
141 coutW(InputArguments) <<
"RooMomentMorphFuncND::ctor(" <<
GetName() <<
") WARNING mref point " << i
142 <<
" is not a constant, taking a snapshot of its value" << endl;
150 for (i = 0; i < mrefpoints.
GetNrows(); ++i) {
152 if (mrefpoints[i] == grid.
array()[j]) {
178 _cacheMgr(other._cacheMgr, this),
179 _parList(
"parList", this, other._parList),
180 _obsList(
"obsList", this, other._obsList),
181 _referenceGrid(other._referenceGrid),
182 _pdfList(
"pdfList", this, other._pdfList),
183 _setting(other._setting),
184 _useHorizMorph(other._useHorizMorph)
213 coutE(InputArguments) <<
"RooMomentMorphFuncND::initialize(" <<
GetName() <<
") ERROR: nPar != nDim"
214 <<
": " << nPar <<
" !=" << nDim << endl;
219 coutE(InputArguments) <<
"RooMomentMorphFuncND::initialize(" <<
GetName() <<
") ERROR: nPdf != nRef"
220 <<
": " << nPdf <<
" !=" << nRef << endl;
225 _M = std::make_unique<TMatrixD>(nPdf, nPdf);
226 _MSqr = std::make_unique<TMatrixD>(depth, depth);
230 vector<vector<double>> dm(nPdf);
231 for (
int k = 0; k < nPdf; ++k) {
233 for (
int idim = 0; idim < nPar; idim++) {
235 dm2.push_back(delta);
240 vector<vector<int>> powers;
241 for (
int idim = 0; idim < nPar; idim++) {
247 powers.push_back(xtmp);
250 vector<vector<int>>
output;
252 int nCombs =
output.size();
254 for (
int k = 0; k < nPdf; ++k) {
256 for (
int i = 0; i < nCombs; i++) {
258 for (
int ix = 0; ix < nPar; ix++) {
259 double delta = dm[k][ix];
278 : _pdfList(other._pdfList), _pdfMap(other._pdfMap), _nref(other._nref)
280 for (
unsigned int i = 0; i < other.
_grid.size(); i++) {
296 vector<int> thisBoundaries;
297 vector<double> thisBoundaryCoordinates;
298 thisBoundaries.push_back(bin_x);
299 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
302 _nref.push_back(thisBoundaryCoordinates);
308 vector<int> thisBoundaries;
309 vector<double> thisBoundaryCoordinates;
310 thisBoundaries.push_back(bin_x);
311 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
312 thisBoundaries.push_back(bin_y);
313 thisBoundaryCoordinates.push_back(_grid[1]->array()[bin_y]);
316 _nref.push_back(thisBoundaryCoordinates);
322 vector<int> thisBoundaries;
323 vector<double> thisBoundaryCoordinates;
324 thisBoundaries.push_back(bin_x);
325 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
326 thisBoundaries.push_back(bin_y);
327 thisBoundaryCoordinates.push_back(_grid[1]->array()[bin_y]);
328 thisBoundaries.push_back(bin_z);
329 thisBoundaryCoordinates.push_back(_grid[2]->array()[bin_z]);
332 _nref.push_back(thisBoundaryCoordinates);
338 vector<double> thisBoundaryCoordinates;
339 int nBins = bins.size();
340 thisBoundaryCoordinates.reserve(nBins);
341 for (
int i = 0; i < nBins; i++) {
342 thisBoundaryCoordinates.push_back(_grid[i]->array()[bins[i]]);
346 _nref.push_back(thisBoundaryCoordinates);
361 vector<RooAbsReal *> meanrv(nPdf * nObs, null);
362 vector<RooAbsReal *> sigmarv(nPdf * nObs, null);
363 vector<RooAbsReal *> myrms(nObs, null);
364 vector<RooAbsReal *> mypos(nObs, null);
365 vector<RooAbsReal *> slope(nPdf * nObs, null);
366 vector<RooAbsReal *> offsets(nPdf * nObs, null);
367 vector<RooAbsReal *> transVar(nPdf * nObs, null);
368 vector<RooAbsReal *> transPdf(nPdf, null);
378 for (
int i = 0; i < 3 * nPdf; ++i) {
379 string fracName =
Form(
"frac_%d", i);
386 }
else if (i < 2 * nPdf) {
394 std::unique_ptr<RooAbsReal> theSum;
401 for (
int i = 0; i < nPdf; ++i) {
402 for (
int j = 0; j < nObs; ++j) {
409 sigmarv[
sij(i, j)] = mom;
410 meanrv[
sij(i, j)] = mom->
mean();
412 ownedComps.
add(*sigmarv[
sij(i, j)]);
417 for (
int j = 0; j < nObs; ++j) {
420 for (
int i = 0; i < nPdf; ++i) {
421 meanList.
add(*meanrv[
sij(i, j)]);
422 rmsList.
add(*sigmarv[
sij(i, j)]);
424 string myrmsName =
Form(
"%s_rms_%d",
GetName(), j);
425 string myposName =
Form(
"%s_pos_%d",
GetName(), j);
426 mypos[j] =
new RooAddition(myposName.c_str(), myposName.c_str(), meanList, coefList2);
427 myrms[j] =
new RooAddition(myrmsName.c_str(), myrmsName.c_str(), rmsList, coefList3);
434 for (
auto const *pdf : static_range_cast<Base_t *>(
_pdfList)) {
436 string pdfName =
Form(
"pdf_%d", i);
440 for (
auto *var : static_range_cast<RooRealVar *>(obsList)) {
442 string slopeName =
Form(
"%s_slope_%d_%d",
GetName(), i, j);
443 string offsetName =
Form(
"%s_offset_%d_%d",
GetName(), i, j);
452 string transVarName =
Form(
"%s_transVar_%d_%d",
GetName(), i, j);
453 transVar[
sij(i, j)] =
new RooLinearVar(transVarName.c_str(), transVarName.c_str(), *var, *slope[
sij(i, j)],
454 *offsets[
sij(i, j)]);
460 ownedComps.
add(*transVar[
sij(i, j)]);
465 transPdfList.
add(*transPdf[i]);
466 ownedComps.
add(*transPdf[i]);
474 theSum = std::make_unique<RooAddPdf>(sumName.c_str(), sumName.c_str(), pdfList, coefList);
476 theSum = std::make_unique<RooRealSumFunc>(sumName.c_str(), sumName.c_str(), pdfList, coefList);
482 theSum->addOwnedComponents(ownedComps);
485 std::string trackerName = std::string(
GetName()) +
"_frac_tracker";
489 std::make_unique<RooChangeTracker>(trackerName.c_str(), trackerName.c_str(),
_parList,
true),
497 std::unique_ptr<RooChangeTracker> &&tracker,
const RooArgList &flist)
498 : _sum(std::move(
sumFunc)), _tracker(std::move(tracker))
525 if (cache->
_tracker->hasChanged(
true)) {
528 return cache->
_sum.get();
536 if (cache->
_tracker->hasChanged(
true)) {
548 return static_cast<RooRealVar *
>(_frac.at(i));
554 return static_cast<RooRealVar *
>(_frac.at(i));
563 double fracLinear(1.);
564 double fracNonLinear(1.);
569 for (
int idim = 0; idim < nPar; idim++) {
571 dm2.push_back(delta);
574 vector<vector<int>> powers;
575 for (
int idim = 0; idim < nPar; idim++) {
581 powers.push_back(xtmp);
584 vector<vector<int>>
output;
586 int nCombs =
output.size();
588 vector<double> deltavec(nPdf, 1.0);
591 for (
int i = 0; i < nCombs; i++) {
593 for (
int ix = 0; ix < nPar; ix++) {
594 double delta = dm2[ix];
597 deltavec[nperm] = tmpDm;
601 double sumposfrac = 0.0;
602 for (
int i = 0; i < nPdf; ++i) {
605 for (
int j = 0; j < nPdf; ++j) {
606 ffrac += (*self.
_M)(j, i) * deltavec[j] * fracNonLinear;
615 const_cast<RooRealVar *
>(frac(i))->setVal(ffrac);
619 const_cast<RooRealVar *
>(frac(nPdf + i))->setVal(ffrac);
620 const_cast<RooRealVar *
>(frac(2 * nPdf + i))->setVal(ffrac);
623 cout <<
"NonLinear fraction " << ffrac << endl;
625 frac(nPdf + i)->Print();
626 frac(2 * nPdf + i)->Print();
631 for (
int i = 0; i < nPdf; ++i) {
632 if (frac(i)->
getVal() < 0)
633 const_cast<RooRealVar *
>(frac(i))->setVal(0.);
634 const_cast<RooRealVar *
>(frac(i))->setVal(frac(i)->getVal() / sumposfrac);
642 for (
int i = 0; i < nPdf; ++i) {
644 const_cast<RooRealVar *
>(frac(i))->setVal(initval);
645 const_cast<RooRealVar *
>(frac(nPdf + i))->setVal(initval);
646 const_cast<RooRealVar *
>(frac(2 * nPdf + i))->setVal(initval);
649 std::vector<double> mtmp;
652 for (
auto *
m : static_range_cast<RooRealVar *>(self.
_parList)) {
653 mtmp.push_back(
m->getVal());
659 vector<double> deltavec(depth, 1.0);
665 for (
int ix = 0; ix < nPar; ix++) {
669 for (
int iperm = 1; iperm <= nPar; ++iperm) {
671 double dtmp = mtmp[xtmp[0]] - self.
_squareVec[0][xtmp[0]];
672 for (
int itmp = 1; itmp < iperm; ++itmp) {
673 dtmp *= mtmp[xtmp[itmp]] - self.
_squareVec[0][xtmp[itmp]];
675 deltavec[nperm + 1] = dtmp;
680 double origFrac1(0.);
681 double origFrac2(0.);
682 for (
int i = 0; i < depth; ++i) {
684 for (
int j = 0; j < depth; ++j) {
685 ffrac += (*self.
_MSqr)(j, i) * deltavec[j] * fracLinear;
689 origFrac1 = frac(self.
_squareIdx[i])->getVal();
701 cout <<
"Linear fraction " << ffrac << endl;
729 vector<vector<double>> boundaries(nPar);
730 for (
int idim = 0; idim < nPar; idim++) {
734 boundaries[idim].push_back(lo);
735 boundaries[idim].push_back(
hi);
738 vector<vector<double>>
output;
742 for (
int isq = 0; isq < depth; isq++) {
743 for (
int iref = 0; iref < nRef; iref++) {
767 for (
int ix = 0; ix < nPar; ix++) {
771 for (
int k = 0; k < depth; ++k) {
777 for (
int iperm = 1; iperm <= nPar; ++iperm) {
779 double dtmp =
_squareVec[k][xtmp[0]] - squareBase[xtmp[0]];
780 for (
int itmp = 1; itmp < iperm; ++itmp) {
781 dtmp *=
_squareVec[k][xtmp[itmp]] - squareBase[xtmp[itmp]];
783 M(k, nperm + 1) = dtmp;
796 if (allVars.
size() == 1) {
803 cout <<
"Currently BinIntegrator only knows how to deal with 1-d " << endl;
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
void setLocalNoDirtyInhibit(bool flag) const
Abstract base class for RooRealVar binning definitions.
Abstract container object that can hold multiple RooAbsArg objects.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * first() const
bool addTyped(const RooAbsCollection &list, bool silent=false)
Adds elements of a given RooAbsCollection to the container if they match the specified type.
bool setRealValue(const char *name, double newVal=0.0, bool verbose=false)
Set value of a RooAbsRealLValue stored in set with given name to newVal No error messages are printed...
const RooArgSet * nset() const
Abstract base class for objects that represent a real value and implements functionality common to al...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
virtual double getValV(const RooArgSet *normalisationSet=nullptr) const
Return value of object.
RooNumIntConfig * specialIntegratorConfig() const
Returns the specialized integrator configuration for this RooAbsReal.
RooAbsMoment * sigma(RooRealVar &obs)
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.
Implements a RooAbsBinning in terms of an array of boundary values, posing no constraints on the choi...
double * array() const override
Return array of boundary values.
Int_t numBoundaries() const override
Return the number boundaries.
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=nullptr)
Setter function without integration set.
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=nullptr, const TNamed *isetRangeName=nullptr)
Getter function without integration set.
bool setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
Represents 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 occurrence of arg 'orig' with arg 'subst'.
RooAbsArg * build(const char *masterCatState, bool verbose=false)
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...
void calculateFractions(const RooMomentMorphFuncND &self, bool verbose=true) const
std::unique_ptr< RooChangeTracker > _tracker
std::unique_ptr< RooAbsReal > _sum
RooArgList containedArgs(Action) override
CacheElem(std::unique_ptr< RooAbsReal > &&sumFunc, std::unique_ptr< RooChangeTracker > &&tracker, const RooArgList &flist)
std::vector< int > _nnuis
void addBinning(const RooAbsBinning &binning)
std::vector< RooAbsBinning * > _grid
std::vector< std::vector< double > > _nref
void addPdf(const RooAbsReal &func, int bin_x)
RooObjCacheManager _cacheMgr
! Transient cache manager
std::unique_ptr< TMatrixD > _MSqr
void findShape(const std::vector< double > &x) const
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooAbsReal * sumFunc(const RooArgSet *nset)
CacheElem * getCache(const RooArgSet *nset) const
std::unique_ptr< TMatrixD > _M
RooArgSet * _curNormSet
! Transient cache manager
bool setBinIntegrator(RooArgSet &allVars)
double getValV(const RooArgSet *set=nullptr) const override
Return value of object.
std::vector< std::vector< double > > _squareVec
int sij(const int &i, const int &j) const
~RooMomentMorphFuncND() override
std::vector< int > _squareIdx
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
Variable that can be changed from the outside.
TMatrixT< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant.
const char * GetName() const override
Returns name of object.
Element * GetMatrixArray()
void cartesianProduct(std::vector< std::vector< T > > &out, std::vector< std::vector< T > > &in)
bool nextCombination(const Iterator first, Iterator k, const Iterator last)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.