109 pdf1(
"pdf1",
"pdf1",this,_pdf1),
110 pdf2(
"pdf2",
"pdf2",this,_pdf2),
112 alpha(
"alpha",
"alpha",this,_alpha),
113 _cacheAlpha(doCacheAlpha),
123 pdf1(
"pdf1",this,other.pdf1),
124 pdf2(
"pdf2",this,other.pdf2),
126 alpha(
"alpha",this,other.alpha),
127 _cacheAlpha(other._cacheAlpha),
155 par1->
add(*par2,
true) ;
173 name.Append(
"_MORPH_") ;
196 double alphaSave =
alpha ;
198 coutP(
Eval) <<
"RooIntegralMorph::fillCacheObject(" <<
GetName() <<
") filling multi-dimensional cache" ;
199 while(slIter->
Next()) {
295 oocoutW(_self,
Eval) <<
"RooIntegralMorph::MorphCacheElem::calcX() WARNING: requested root finding for unphysical CDF value " <<
y << endl ;
299 double xmax = _x->getMax(
"cache") ;
300 double xmin = _x->getMin(
"cache") ;
308 return _alpha->getVal()*
x1 + (1-_alpha->getVal())*
x2 ;
316 double xmax = _x->getMax(
"cache") ;
317 double xmin = _x->getMin(
"cache") ;
327 double xsave = _self->x ;
334 _yatX.resize(_x->numBins(
"cache")+1);
335 _calcX.resize(_x->numBins(
"cache")+1);
341 Int_t nbins = _x->numBins(
"cache") ;
344 for (
int i=0 ; i<nbins ; i++) {
353 for (
int i=0 ; i<10 ; i++) {
356 double offset = _yatX[_yatXmin] ;
357 double delta = (_yatX[_yatXmax] - _yatX[_yatXmin])/10 ;
362 double X = calcX(
y,ok) ;
371 Int_t igapLow = _yatXmin+1 ;
374 Int_t igapHigh = igapLow+1 ;
375 while(igapHigh<(_yatXmax) && _yatX[igapHigh]<0) igapHigh++ ;
378 fillGap(igapLow-1,igapHigh) ;
381 if (igapHigh>=_yatXmax-1) break ;
382 igapLow = igapHigh+1 ;
386 double xmax = _x->getMax(
"cache") ;
387 double xmin = _x->getMin(
"cache") ;
388 double binw = (
xmax-
xmin)/_x->numBins(
"cache") ;
389 for (
int i=_yatXmin+1 ; i<_yatXmax-1 ; i++) {
392 double xBinC =
xmin + (i+0.5)*binw ;
393 double xOffset = xBinC-_calcX[i] ;
394 if (
fabs(xOffset/binw)>1
e-3) {
395 double slope = (_yatX[i+1]-_yatX[i-1])/(_calcX[i+1]-_calcX[i-1]) ;
396 double newY = _yatX[i] + slope*xOffset ;
403 for (
int i=0; i<_yatXmin ; i++) {
409 double x1 = _x->getMin(
"cache");
410 double x2 = _x->getMin(
"cache");
412 double xMax = _x->getMax(
"cache");
415 for (
int i=_yatXmin ; i<_yatXmax ; i++) {
417 double y = _yatX[i] ;
423 _rf1->findRoot(
x1,
x1,xMax,
y) ;
424 _rf2->findRoot(
x2,
x2,xMax,
y) ;
426 _x->setVal(
x1) ;
double f1x1 = _pdf1->getVal(&nsetTmp) ;
427 _x->setVal(
x2) ;
double f2x2 = _pdf2->getVal(&nsetTmp) ;
428 double fbarX = f1x1*f2x2 / ( _alpha->getVal()*f2x2 + (1-_alpha->getVal())*f1x1 ) ;
435 for (
int i=_yatXmax+1 ; i<nbins ; i++) {
441 pdf()->setUnitNorm(
true) ;
444 oocxcoutD(_self,
Eval) <<
"RooIntegralMorph::MorphCacheElem::calculate(" << _self->GetName() <<
") calculation required " << _ccounter <<
" samplings of cdfs" << endl ;
459 oocoutE(_self,
Eval) <<
"RooIntegralMorph::MorphCacheElme::fillGap(" << _self->GetName() <<
"): ERROR in fillgap " << ixlo <<
" = " << ixhi
460 <<
" splitPoint= " << splitPoint <<
" _yatX[ixlo] = " << _yatX[ixlo] << endl ;
463 oocoutE(_self,
Eval) <<
"RooIntegralMorph::MorphCacheElme::fillGap(" << _self->GetName() <<
"): ERROR in fillgap " << ixlo <<
" = " << ixhi
464 <<
" splitPoint " << splitPoint <<
" _yatX[ixhi] = " << _yatX[ixhi] << endl ;
468 double ymid = _yatX[ixlo]*splitPoint + _yatX[ixhi]*(1-splitPoint) ;
470 double Xmid = calcX(ymid,ok) ;
472 oocoutW(_self,
Eval) <<
"RooIntegralMorph::MorphCacheElem::fillGap(" << _self->GetName() <<
") unable to calculate midpoint in gap ["
473 << ixlo <<
"," << ixhi <<
"], resorting to interpolation" << endl ;
474 interpolateGap(ixlo,ixhi) ;
477 Int_t iX = binX(Xmid) ;
478 double cq = (Xmid-_calcX[ixlo])/(_calcX[ixhi]-_calcX[ixlo])-0.5 ;
486 if (
fabs(cq)<0.01 ||
fabs(cq*(ixhi-ixlo))<0.1 || ymid<_ycutoff ) {
490 interpolateGap(ixlo,iX) ;
493 interpolateGap(iX,ixhi) ;
500 if (splitPoint<0.95) {
502 double newSplit = splitPoint + 0.5*(1-splitPoint) ;
503 fillGap(ixlo,ixhi,newSplit) ;
506 interpolateGap(ixlo,ixhi) ;
509 }
else if (iX==ixhi) {
512 if (splitPoint>0.05) {
513 double newSplit = splitPoint/2 ;
514 fillGap(ixlo,ixhi,newSplit) ;
517 interpolateGap(ixlo,ixhi) ;
541 double xmax = _x->getMax(
"cache") ;
542 double xmin = _x->getMin(
"cache") ;
543 double binw = (
xmax-
xmin)/_x->numBins(
"cache") ;
546 double deltaY = (_yatX[ixhi]-_yatX[ixlo])/((_calcX[ixhi]-_calcX[ixlo])/binw) ;
549 double xBinC =
xmin + (ixlo+0.5)*binw ;
550 double xOffset = xBinC-_calcX[ixlo] ;
552 for (
int j=ixlo+1 ; j<ixhi ; j++) {
553 _yatX[j] = _yatX[ixlo]+(xOffset+(j-ixlo))*deltaY ;
554 _calcX[j] =
xmin + (j+0.5)*binw ;
568 double xmin = _x->getMin(
"cache") ;
569 double xmax = _x->getMax(
"cache") ;
570 Int_t nbins = _x->numBins(
"cache") ;
574 double ymin=0.1,yminSave(-1) ;
575 double Xsave(-1),Xlast=
xmax ;
582 oocxcoutD(_self,
Eval) <<
"RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName() <<
") findMin: x1 = " <<
x1 <<
" x2 = " <<
x2 <<
" ok = " << (ok?
"T":
"F") << endl ;
588 double X = _alpha->getVal()*
x1 + (1-_alpha->getVal())*
x2 ;
596 _yatX[_yatXmin] =
ymin ;
597 _calcX[_yatXmin] = X ;
605 if (
ymin<_ycutoff) break ;
607 _yatX[_yatXmin] = yminSave ;
608 _calcX[_yatXmin] = Xsave ;
613 double deltaymax=0.1, deltaymaxSave(-1) ;
616 ok &= _rf1->findRoot(
x1,
xmin,
xmax,1-deltaymax) ;
617 ok &= _rf2->findRoot(
x2,
xmin,
xmax,1-deltaymax) ;
619 oocxcoutD(_self,
Eval) <<
"RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName() <<
") findMax: x1 = " <<
x1 <<
" x2 = " <<
x2 <<
" ok = " << (ok?
"T":
"F") << endl ;
625 double X = _alpha->getVal()*
x1 + (1-_alpha->getVal())*
x2 ;
633 _yatX[_yatXmax] = 1-deltaymax ;
634 _calcX[_yatXmax] = X ;
635 deltaymaxSave = deltaymax ;
639 deltaymax /=
sqrt(10.) ;
642 if (deltaymax<_ycutoff) break ;
645 _yatX[_yatXmax] = 1-deltaymaxSave ;
646 _calcX[_yatXmax] = Xsave ;
650 for (
int i=0 ; i<_yatXmin ; i++) _yatX[i] = -2 ;
651 for (
int i=_yatXmax+1 ; i<nbins; i++) _yatX[i] = -2 ;
652 oocxcoutD(_self,
Eval) <<
"RooIntegralMorph::findRange(" << _self->GetName() <<
"): ymin = " << _yatX[_yatXmin] <<
" ymax = " << _yatX[_yatXmax] << endl;
653 oocxcoutD(_self,
Eval) <<
"RooIntegralMorph::findRange(" << _self->GetName() <<
"): xmin = " << _calcX[_yatXmin] <<
" xmax = " << _calcX[_yatXmax] << endl;
674 orderedObs.
add(obs) ;
677 orderedObs.
remove(*obsX) ;
678 orderedObs.
add(*obsX) ;
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
RooArgSet * getParameters(const RooAbsData *data, bool stripDisconnected=true) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don't match any of...
RooArgList containedArgs(Action) override
Returns all RooAbsArg objects contained in the cache element.
RooAbsCachedPdf is the abstract base class for p.d.f.s that need or want to cache their evaluate() ou...
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsReal * createCdf(const RooArgSet &iset, const RooArgSet &nset=RooArgSet())
Create a cumulative distribution function of this p.d.f in terms of the observables listed in iset.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
RooAbsFunc * bindVars(const RooArgSet &vars, const RooArgSet *nset=0, bool clipInvalid=false) const
Create an interface adaptor f(vars) that binds us to the specified variables (in arbitrary order).
virtual double offset() const
RooArgList is a container object that can hold multiple RooAbsArg objects.
RooAbsArg * absArg() const
Return pointer to contained argument.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Implement the abstract 1-dimensional root finding interface using the Brent-Decker method.
void setTol(double tol)
Set convergence tolerance parameter.
TIterator * sliceIterator(RooAbsArg &sliceArg, const RooArgSet &otherArgs)
Create an iterator over all bins in a slice defined by the subset of observables listed in sliceArg.
const RooArgSet * get() const override
Get bin centre of current bin.
void setUnitNorm(bool flag)
~MorphCacheElem() override
Destructor.
void calculate(TIterator *iter)
Calculate shape of p.d.f for x,alpha values defined by dIter iterator over cache histogram.
void interpolateGap(Int_t ixlo, Int_t ixhi)
Fill empty histogram bins between ixlo and ixhi with values obtained from linear interpolation of ixl...
MorphCacheElem(RooIntegralMorph &self, const RooArgSet *nset)
Construct of cache element, copy relevant input from RooIntegralMorph, create the cdfs from the input...
void fillGap(Int_t ixlo, Int_t ixhi, double splitPoint=0.5)
Fill all empty histogram bins between bins ixlo and ixhi.
RooBrentRootFinder * _rf2
RooBrentRootFinder * _rf1
void findRange()
Determine which range of y values can be mapped to x values from the numeric inversion of the input c...
RooArgList containedArgs(Action) override
Return all RooAbsArg components contained in this cache.
double calcX(double y, bool &ok)
Calculate the x value of the output p.d.f at the given cdf value y.
Int_t binX(double x)
Return the bin number enclosing the given x value.
Class RooIntegralMorph is an implementation of the histogram interpolation technique described by Ale...
friend class MorphCacheElem
PdfCacheElem * createCache(const RooArgSet *nset) const override
Create and return a derived MorphCacheElem.
const char * inputBaseName() const override
Return base name component for cache components in this case a string encoding the names of both end ...
void preferredObservableScanOrder(const RooArgSet &obs, RooArgSet &orderedObs) const override
Indicate to the RooAbsCachedPdf base class that for the filling of the cache the traversal of the x s...
void fillCacheObject(PdfCacheElem &cache) const override
Fill the cache with the interpolated shape.
RooArgSet * actualParameters(const RooArgSet &nset) const override
Parameters of the cache.
double evaluate() const override
Dummy.
RooArgSet * actualObservables(const RooArgSet &nset) const override
Observable to be cached for given choice of normalization.
RooRealVar represents a variable that can be changed from the outside.
const T & arg() const
Return reference to object held in proxy.
Iterator abstract base class.
virtual TObject * Next()=0
const char * GetName() const override
Returns name of object.
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)