Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooIntegralMorph.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2 * Project: RooFit *
3 * *
4 * Copyright (c) 2000-2005, Regents of the University of California *
5 * and Stanford University. All rights reserved. *
6 * *
7 * Redistribution and use in source and binary forms, *
8 * with or without modification, are permitted according to the terms *
9 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
10 *****************************************************************************/
11
12/** \class RooIntegralMorph
13 \ingroup Roofit
14
15Class RooIntegralMorph is an implementation of the histogram interpolation
16technique described by Alex Read in 'NIM A 425 (1999) 357-369 'Linear interpolation of histograms'
17for continuous functions rather than histograms. The interpolation method, in short,
18works as follows.
19
20 - Given a p.d.f f1(x) with c.d.f F1(x) and p.d.f f2(x) with c.d.f F2(x)
21
22 - One finds takes a value 'y' of both c.d.fs and determines the corresponding x
23 values x(1,2) at which F(1,2)(x)==y.
24
25 - The value of the interpolated p.d.f fbar(x) is then calculated as
26 fbar(alpha*x1+(1-alpha)*x2) = f1(x1)*f2(x2) / ( alpha*f2(x2) + (1-alpha)*f1(x1) ) ;
27
28From a technical point of view class RooIntegralMorph is a p.d.f that takes
29two input p.d.fs f1(x,p) an f2(x,q) and an interpolation parameter to
30make a p.d.f fbar(x,p,q,alpha). The shapes f1 and f2 are always taken
31to be end the end-points of the parameter alpha, regardless of what
32the those numeric values are.
33
34Since the value of fbar(x) cannot be easily calculated for a given value
35of x, class RooIntegralMorph is an implementation of RooAbsCachedPdf and
36calculates the shape of the interpolated p.d.f. fbar(x) for all values
37of x for a given value of alpha,p,q and caches these values in a histogram
38(as implemented by RooAbsCachedPdf). The binning granularity of the cache
39can be controlled by the binning named "cache" on the RooRealVar representing
40the observable x. The fbar sampling algorithm is based on a recursive division
41mechanism with a built-in precision cutoff: First an initial sampling in
4264 equally spaced bins is made. Then the value of fbar is calculated in
43the center of each gap. If the calculated value deviates too much from
44the value obtained by linear interpolation from the edge bins, gap
45is recursively divided. This strategy makes it possible to define a very
46fine cache sampling (e.g. 1000 or 10000) bins without incurring a
47corresponding CPU penalty.
48
49Note on numeric stability of the algorithm. Since the algorithm relies
50on a numeric inversion of cumulative distributions functions, some precision
51may be lost at the 'edges' of the same (i.e. at regions in x where the
52c.d.f. value is close to zero or one). The general sampling strategy is
53to start with 64 equally spaces samples in the range y=(0.01-0.99).
54Then the y ranges are pushed outward by reducing y (or the distance of y to 1.0)
55by a factor of sqrt(10) iteratively up to the point where the corresponding
56x value no longer changes significantly. For p.d.f.s with very flat tails
57such as Gaussians some part of the tail may be lost due to limitations
58in numeric precision in the CDF inversion step.
59
60An effect related to the above limitation in numeric precision should
61be anticipated when floating the alpha parameter in a fit. If a p.d.f
62with such flat tails is fitted, it is likely that the dataset contains
63events in the flat tail region. If the alpha parameter is varied, the
64likelihood contribution from such events may exhibit discontinuities
65in alpha, causing discontinuities in the summed likelihood as well
66that will cause convergence problems in MINUIT. To mitigate this effect
67one can use the setCacheAlpha() method to instruct RooIntegralMorph
68to construct a two-dimensional cache for its output values in both
69x and alpha. If linear interpolation is requested on the resulting
70output histogram, the resulting interpolation of the p.d.f in the
71alpha dimension will smooth out the discontinuities in the tail regions
72result in a continuous likelihood distribution that can be fitted.
73An added advantage of the cacheAlpha option is that if parameters
74p,q of f1,f2 are fixed, the cached values in RooIntegralMorph are
75valid for the entire fit session and do not need to be recalculated
76for each change in alpha, which may result an considerable increase
77in calculation speed.
78
79**/
80
81#include "Riostream.h"
82
83#include "RooIntegralMorph.h"
84#include "RooAbsCategory.h"
85#include "RooBrentRootFinder.h"
86#include "RooAbsFunc.h"
87#include "RooRealVar.h"
88#include "RooDataHist.h"
89#include "TH1.h"
90
91 using std::flush, std::endl;
92
93
94////////////////////////////////////////////////////////////////////////////////
95/// Constructor with observables x, pdf shapes pdf1 and pdf2 which represent
96/// the shapes at the end points of the interpolation parameter alpha
97/// If doCacheAlpha is true, a two-dimensional cache is constructed in
98/// both alpha and x
99
100RooIntegralMorph::RooIntegralMorph(const char *name, const char *title,
101 RooAbsReal& _pdf1,
102 RooAbsReal& _pdf2,
103 RooAbsReal& _x,
104 RooAbsReal& _alpha,
105 bool doCacheAlpha) :
106 RooAbsCachedPdf(name,title,2),
107 pdf1("pdf1","pdf1",this,_pdf1),
108 pdf2("pdf2","pdf2",this,_pdf2),
109 x("x","x",this,_x),
110 alpha("alpha","alpha",this,_alpha),
111 _cacheAlpha(doCacheAlpha)
112{
113}
114
115////////////////////////////////////////////////////////////////////////////////
116/// Copy constructor
117
120 pdf1("pdf1",this,other.pdf1),
121 pdf2("pdf2",this,other.pdf2),
122 x("x",this,other.x),
123 alpha("alpha",this,other.alpha),
124 _cacheAlpha(other._cacheAlpha)
125{
126}
127
128////////////////////////////////////////////////////////////////////////////////
129/// Observable to be cached for given choice of normalization.
130/// Returns the 'x' observable unless doCacheAlpha is set in which
131/// case a set with both x and alpha
132
134{
135 RooArgSet* obs = new RooArgSet ;
136 if (_cacheAlpha) {
137 obs->add(alpha.arg()) ;
138 }
139 obs->add(x.arg()) ;
141}
142
143////////////////////////////////////////////////////////////////////////////////
144/// Parameters of the cache. Returns parameters of both pdf1 and pdf2
145/// and parameter cache, in case doCacheAlpha is not set.
146
148{
149 std::unique_ptr<RooArgSet> par1{pdf1->getParameters(static_cast<RooArgSet*>(nullptr))};
151 pdf2->getParameters(nullptr, par2);
152 par1->add(par2,true) ;
153 par1->remove(x.arg(),true,true) ;
154 if (!_cacheAlpha) {
155 par1->add(alpha.arg()) ;
156 }
157 return RooFit::makeOwningPtr(std::move(par1));
158}
159
160////////////////////////////////////////////////////////////////////////////////
161/// Return base name component for cache components in this case
162/// a string encoding the names of both end point p.d.f.s
163
165{
166 static TString name ;
167
168 name = pdf1.arg().GetName() ;
169 name.Append("_MORPH_") ;
170 name.Append(pdf2.arg().GetName()) ;
171 return name.Data() ;
172}
173
174////////////////////////////////////////////////////////////////////////////////
175/// Fill the cache with the interpolated shape.
176
178{
179 MorphCacheElem& mcache = static_cast<MorphCacheElem&>(cache) ;
180
181 // If cacheAlpha is true employ slice iterator here to fill all slices
182
183 if (!_cacheAlpha) {
184
185 std::unique_ptr<TIterator> dIter{cache.hist()->sliceIterator(const_cast<RooAbsReal&>(x.arg()),RooArgSet())};
186 mcache.calculate(dIter.get());
187
188 } else {
189 std::unique_ptr<TIterator> slIter{cache.hist()->sliceIterator(const_cast<RooAbsReal&>(alpha.arg()),RooArgSet())};
190
191 double alphaSave = alpha ;
193 coutP(Eval) << "RooIntegralMorph::fillCacheObject(" << GetName() << ") filling multi-dimensional cache" ;
194 while(slIter->Next()) {
195 alphaSet.assign(*cache.hist()->get()) ;
196 std::unique_ptr<TIterator> dIter{cache.hist()->sliceIterator(const_cast<RooAbsReal&>(x.arg()),RooArgSet(alpha.arg()))};
197 mcache.calculate(dIter.get());
198 ccoutP(Eval) << "." << flush;
199 }
200 ccoutP(Eval) << std::endl;
201
202 const_cast<RooIntegralMorph*>(this)->alpha = alphaSave ;
203 }
204}
205
206////////////////////////////////////////////////////////////////////////////////
207/// Create and return a derived MorphCacheElem.
208
210{
211 return new MorphCacheElem(const_cast<RooIntegralMorph&>(*this),nset) ;
212}
213
214////////////////////////////////////////////////////////////////////////////////
215/// Return all RooAbsArg components contained in this cache
216
218{
220 ret.add(PdfCacheElem::containedArgs(action)) ;
221 ret.add(*_self) ;
222 ret.add(*_pdf1) ;
223 ret.add(*_pdf2) ;
224 ret.add(*_x ) ;
225 ret.add(*_alpha) ;
226 ret.add(*_c1) ;
227 ret.add(*_c2) ;
228
229 return ret ;
230}
231
232////////////////////////////////////////////////////////////////////////////////
233/// Construct of cache element, copy relevant input from RooIntegralMorph,
234/// create the cdfs from the input p.d.fs and instantiate the root finders
235/// on the cdfs to perform the inversion.
236
239 _self(&self),
240 _pdf1(static_cast<RooAbsPdf *>(self.pdf1.absArg())),
241 _pdf2(static_cast<RooAbsPdf *>(self.pdf2.absArg())),
242 _x(static_cast<RooRealVar *>(self.x.absArg())),
243 _alpha(static_cast<RooAbsReal *>(self.alpha.absArg())),
244 _yatXmin(0),
245 _yatXmax(0),
246 _ccounter(0),
247 _ycutoff(1e-7)
248{
249 // Mark in base class that normalization of cached pdf is invariant under pdf parameters
250
251 _nset = std::make_unique<RooArgSet>(*_x);
252
253 _c1 = std::unique_ptr<RooAbsReal>{_pdf1->createCdf(*_x)};
254 _c2 = std::unique_ptr<RooAbsReal>{_pdf2->createCdf(*_x)};
255 _cb1 = std::unique_ptr<RooAbsFunc>{_c1->bindVars(*_x,_nset.get())};
256 _cb2 = std::unique_ptr<RooAbsFunc>{_c2->bindVars(*_x,_nset.get())};
257
258 _rf1 = std::make_unique<RooBrentRootFinder>(*_cb1);
259 _rf2 = std::make_unique<RooBrentRootFinder>(*_cb2);
260
261 _rf1->setTol(1e-12) ;
262 _rf2->setTol(1e-12) ;
263
264 // _yatX = 0 ;
265 // _calcX = 0 ;
266
267 // Must do this here too: fillCache() may not be called if cache contents is retrieved from EOcache
268 pdf()->setUnitNorm(true) ;
269}
270
271////////////////////////////////////////////////////////////////////////////////
272/// Destructor
273
277
278////////////////////////////////////////////////////////////////////////////////
279/// Calculate the x value of the output p.d.f at the given cdf value y.
280/// The ok boolean is filled with the success status of the operation.
281
283{
284 if (y<0 || y>1) {
285 oocoutW(_self,Eval) << "RooIntegralMorph::MorphCacheElem::calcX() WARNING: requested root finding for unphysical CDF value " << y << std::endl ;
286 }
287 double x1;
288 double x2;
289
290 double xmax = _x->getMax("cache") ;
291 double xmin = _x->getMin("cache") ;
292
293 ok=true ;
294 ok &= _rf1->findRoot(x1,xmin,xmax,y) ;
295 ok &= _rf2->findRoot(x2,xmin,xmax,y) ;
296 if (!ok) return 0 ;
297 _ccounter++ ;
298
299 return _alpha->getVal()*x1 + (1-_alpha->getVal())*x2 ;
300}
301
302////////////////////////////////////////////////////////////////////////////////
303/// Return the bin number enclosing the given x value
304
306{
307 double xmax = _x->getMax("cache") ;
308 double xmin = _x->getMin("cache") ;
309 return (Int_t)(_x->numBins("cache")*(X-xmin)/(xmax-xmin)) ;
310}
311
312////////////////////////////////////////////////////////////////////////////////
313/// Calculate shape of p.d.f for x,alpha values
314/// defined by dIter iterator over cache histogram
315
317{
318 double xsave = _self->x ;
319
320 // if (!_yatX) {
321 // _yatX = new double[_x->numBins("cache")+1] ;
322 // _calcX = new double[_x->numBins("cache")+1] ;
323 // }
324
325 _yatX.resize(_x->numBins("cache")+1);
326 _calcX.resize(_x->numBins("cache")+1);
327
328 _ccounter = 0 ;
329
330 // Get number of bins from PdfCacheElem histogram
331 Int_t nbins = _x->numBins("cache") ;
332
333 // Initialize yatX array to 'un-calculated values (-1)'
334 for (int i=0 ; i<nbins ; i++) {
335 _yatX[i] = -1 ;
336 _calcX[i] = 0 ;
337 }
338
339 // Find low and high point
340 findRange() ;
341
342 // Perform initial scan of 100 points
343 for (int i=0 ; i<10 ; i++) {
344
345 // Take a point in y
346 double offset = _yatX[_yatXmin] ;
347 double delta = (_yatX[_yatXmax] - _yatX[_yatXmin])/10 ;
348 double y = offset + i*delta ;
349
350 // Calculate corresponding X
351 bool ok ;
352 double X = calcX(y,ok) ;
353 if (ok) {
354 Int_t iX = binX(X) ;
355 _yatX[iX] = y ;
356 _calcX[iX] = X ;
357 }
358 }
359
360 // Now take an iteration filling the 'gaps'
361 Int_t igapLow = _yatXmin+1 ;
362 while(true) {
363 // Find next gap
365 while(igapHigh<(_yatXmax) && _yatX[igapHigh]<0) igapHigh++ ;
366
367 // Fill the gap (iteratively and/or using interpolation)
368 fillGap(igapLow-1,igapHigh) ;
369
370 // Terminate after processing of last gap
371 if (igapHigh>=_yatXmax-1) break ;
372 igapLow = igapHigh+1 ;
373 }
374
375 // Make one more iteration to recalculate Y value at bin centers
376 double xmax = _x->getMax("cache") ;
377 double xmin = _x->getMin("cache") ;
378 double binw = (xmax-xmin)/_x->numBins("cache") ;
379 for (int i=_yatXmin+1 ; i<_yatXmax-1 ; i++) {
380
381 // Calculate additional offset to apply if bin ixlo does not have X value calculated at bin center
382 double xBinC = xmin + (i+0.5)*binw ;
383 double xOffset = xBinC-_calcX[i] ;
384 if (std::abs(xOffset/binw)>1e-3) {
385 double slope = (_yatX[i+1]-_yatX[i-1])/(_calcX[i+1]-_calcX[i-1]) ;
386 double newY = _yatX[i] + slope*xOffset ;
387 //cout << "bin " << i << " needs to be re-centered " << xOffset/binw << " slope = " << slope << " origY = " << _yatX[i] << " newY = " << newY << std::endl ;
388 _yatX[i] = newY ;
389 }
390 }
391
392 // Zero output histogram below lowest calculable X value
393 for (int i=0; i<_yatXmin ; i++) {
394 dIter->Next() ;
395 //_hist->get(i) ;
396 hist()->set(0) ;
397 }
398
399 double x1 = _x->getMin("cache");
400 double x2 = _x->getMin("cache");
401
402 double xMax = _x->getMax("cache");
403
404 // Transfer calculated values to histogram
405 for (int i=_yatXmin ; i<_yatXmax ; i++) {
406
407 double y = _yatX[i] ;
408
409 // Little optimization here exploiting the fact that th cumulative
410 // distribution functions increase monotonically, so we already know that
411 // the next x-value must be higher than the last one as y is increasing. So
412 // we can use the previous x values as lower bounds.
413 _rf1->findRoot(x1,x1,xMax,y) ;
414 _rf2->findRoot(x2,x2,xMax,y) ;
415
416 _x->setVal(x1);
417 double f1x1 = _pdf1->getVal(_nset.get());
418 _x->setVal(x2);
419 double f2x2 = _pdf2->getVal(_nset.get());
420 double fbarX = f1x1*f2x2 / ( _alpha->getVal()*f2x2 + (1-_alpha->getVal())*f1x1 ) ;
421
422 dIter->Next() ;
423 //_hist->get(i) ;
424 hist()->set(fbarX) ;
425 }
426 // Zero output histogram above highest calculable X value
427 for (int i=_yatXmax+1 ; i<nbins ; i++) {
428 dIter->Next() ;
429 //_hist->get(i) ;
430 hist()->set(0) ;
431 }
432
433 pdf()->setUnitNorm(true) ;
434 _self->x = xsave ;
435
436 oocxcoutD(_self,Eval) << "RooIntegralMorph::MorphCacheElem::calculate(" << _self->GetName() << ") calculation required " << _ccounter << " samplings of cdfs" << std::endl ;
437}
438
439////////////////////////////////////////////////////////////////////////////////
440/// Fill all empty histogram bins between bins ixlo and ixhi. The value of 'splitPoint'
441/// defines the split point for the recursive division strategy to fill the gaps
442/// If the midpoint value of y is very close to the midpoint in x, use interpolation
443/// to fill the gaps, otherwise the intervals again.
444
446{
447 // CONVENTION: _yatX[ixlo] is filled, _yatX[ixhi] is filled, elements in between are empty
448 // std::cout << "fillGap: gap from _yatX[" << ixlo << "]=" << _yatX[ixlo] << " to _yatX[" << ixhi << "]=" << _yatX[ixhi] << ", size = " << ixhi-ixlo << std::endl ;
449
450 if (_yatX[ixlo]<0) {
451 oocoutE(_self,Eval) << "RooIntegralMorph::MorphCacheElme::fillGap(" << _self->GetName() << "): ERROR in fillgap " << ixlo << " = " << ixhi
452 << " splitPoint= " << splitPoint << " _yatX[ixlo] = " << _yatX[ixlo] << std::endl ;
453 }
454 if (_yatX[ixhi]<0) {
455 oocoutE(_self,Eval) << "RooIntegralMorph::MorphCacheElme::fillGap(" << _self->GetName() << "): ERROR in fillgap " << ixlo << " = " << ixhi
456 << " splitPoint " << splitPoint << " _yatX[ixhi] = " << _yatX[ixhi] << std::endl ;
457 }
458
459 // Determine where half-way Y value lands
460 double ymid = _yatX[ixlo]*splitPoint + _yatX[ixhi]*(1-splitPoint) ;
461 bool ok ;
462 double Xmid = calcX(ymid,ok) ;
463 if (!ok) {
464 oocoutW(_self,Eval) << "RooIntegralMorph::MorphCacheElem::fillGap(" << _self->GetName() << ") unable to calculate midpoint in gap ["
465 << ixlo << "," << ixhi << "], resorting to interpolation" << std::endl ;
466 interpolateGap(ixlo,ixhi) ;
467 }
468
469 Int_t iX = binX(Xmid) ;
470 double cq = (Xmid-_calcX[ixlo])/(_calcX[ixhi]-_calcX[ixlo])-0.5 ;
471
472 // Store midway point
473 _yatX[iX] = ymid ;
474 _calcX[iX] = Xmid ;
475
476
477 // Policy: If centration quality is better than 1% OR better than 1/10 of a bin, fill interval with linear interpolation
478 if (std::abs(cq)<0.01 || std::abs(cq*(ixhi-ixlo))<0.1 || ymid<_ycutoff ) {
479
480 // Fill remaining gaps on either side with linear interpolation
481 if (iX-ixlo>1) {
482 interpolateGap(ixlo,iX) ;
483 }
484 if (ixhi-iX>1) {
485 interpolateGap(iX,ixhi) ;
486 }
487
488 } else {
489
490 if (iX==ixlo) {
491
492 if (splitPoint<0.95) {
493 // Midway value lands on lowest bin, retry split with higher split point
494 double newSplit = splitPoint + 0.5*(1-splitPoint) ;
495 fillGap(ixlo,ixhi,newSplit) ;
496 } else {
497 // Give up and resort to interpolation
498 interpolateGap(ixlo,ixhi) ;
499 }
500
501 } else if (iX==ixhi) {
502
503 // Midway value lands on highest bin, retry split with lower split point
504 if (splitPoint>0.05) {
505 double newSplit = splitPoint/2 ;
506 fillGap(ixlo,ixhi,newSplit) ;
507 } else {
508 // Give up and resort to interpolation
509 interpolateGap(ixlo,ixhi) ;
510 }
511
512 } else {
513
514 // Midway point reasonable, iterate on interval on both sides
515 if (iX-ixlo>1) {
516 fillGap(ixlo,iX) ;
517 }
518 if (ixhi-iX>1) {
519 fillGap(iX,ixhi) ;
520 }
521 }
522 }
523}
524
525////////////////////////////////////////////////////////////////////////////////
526/// Fill empty histogram bins between ixlo and ixhi with values obtained
527/// from linear interpolation of ixlo,ixhi elements.
528
530{
531 //cout << "filling gap with linear interpolation ixlo=" << ixlo << " ixhi=" << ixhi << std::endl ;
532
533 double xmax = _x->getMax("cache") ;
534 double xmin = _x->getMin("cache") ;
535 double binw = (xmax-xmin)/_x->numBins("cache") ;
536
537 // Calculate deltaY in terms of actual X difference calculate, not based on nominal bin width
538 double deltaY = (_yatX[ixhi]-_yatX[ixlo])/((_calcX[ixhi]-_calcX[ixlo])/binw) ;
539
540 // Calculate additional offset to apply if bin ixlo does not have X value calculated at bin center
541 double xBinC = xmin + (ixlo+0.5)*binw ;
542 double xOffset = xBinC-_calcX[ixlo] ;
543
544 for (int j=ixlo+1 ; j<ixhi ; j++) {
545 _yatX[j] = _yatX[ixlo]+(xOffset+(j-ixlo))*deltaY ;
546 _calcX[j] = xmin + (j+0.5)*binw ;
547 }
548
549}
550
551////////////////////////////////////////////////////////////////////////////////
552/// Determine which range of y values can be mapped to x values
553/// from the numeric inversion of the input c.d.fs.
554/// Start with a y range of [0.1-0.9] and push boundaries
555/// outward with a factor of 1/sqrt(10). Stop iteration if
556/// inverted x values no longer change
557
559{
560 double xmin = _x->getMin("cache") ;
561 double xmax = _x->getMax("cache") ;
562 Int_t nbins = _x->numBins("cache") ;
563
564 double x1;
565 double x2;
566 bool ok = true ;
567 double ymin = 0.1;
568 double yminSave(-1);
569 double Xsave(-1);
570 double Xlast = xmax;
571
572 // Find lowest Y value that can be measured
573 // Start at 0.1 and iteratively lower limit by sqrt(10)
574 while(true) {
575 ok &= _rf1->findRoot(x1,xmin,xmax,ymin) ;
576 ok &= _rf2->findRoot(x2,xmin,xmax,ymin) ;
577 oocxcoutD(_self,Eval) << "RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName() << ") findMin: x1 = " << x1 << " x2 = " << x2 << " ok = " << (ok?"T":"F") << std::endl ;
578
579 // Terminate in case of non-convergence
580 if (!ok) break ;
581
582 // Terminate if value of X no longer moves by >0.1 bin size
583 double X = _alpha->getVal()*x1 + (1-_alpha->getVal())*x2 ;
584 if (std::abs(X-Xlast)/(xmax-xmin)<0.0001) {
585 break ;
586 }
587 Xlast=X ;
588
589 // Store new Y value
590 _yatXmin = (Int_t)(nbins*(X-xmin)/(xmax-xmin)) ;
591 _yatX[_yatXmin] = ymin ;
592 _calcX[_yatXmin] = X ;
593 yminSave = ymin ;
594 Xsave=X ;
595
596 // Reduce ymin by half an order of magnitude
597 ymin /=sqrt(10.) ;
598
599 // Emergency break
600 if (ymin<_ycutoff) break ;
601 }
602 _yatX[_yatXmin] = yminSave ;
603 _calcX[_yatXmin] = Xsave ;
604
605 // Find highest Y value that can be measured
606 // Start at 1 - 0.1 and iteratively lower delta by sqrt(10)
607 ok = true ;
608 double deltaymax = 0.1;
609 double deltaymaxSave(-1);
610 Xlast=xmin ;
611 while(true) {
612 ok &= _rf1->findRoot(x1,xmin,xmax,1-deltaymax) ;
613 ok &= _rf2->findRoot(x2,xmin,xmax,1-deltaymax) ;
614
615 oocxcoutD(_self,Eval) << "RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName() << ") findMax: x1 = " << x1 << " x2 = " << x2 << " ok = " << (ok?"T":"F") << std::endl ;
616
617 // Terminate in case of non-convergence
618 if (!ok) break ;
619
620 // Terminate if value of X no longer moves by >0.1 bin size
621 double X = _alpha->getVal()*x1 + (1-_alpha->getVal())*x2 ;
622 if (std::abs(X-Xlast)/(xmax-xmin)<0.0001) {
623 break ;
624 }
625 Xlast=X ;
626
627 // Store new Y value
628 _yatXmax = (Int_t)(nbins*(X-xmin)/(xmax-xmin)) ;
629 _yatX[_yatXmax] = 1-deltaymax ;
630 _calcX[_yatXmax] = X ;
632 Xsave=X ;
633
634 // Reduce ymin by half an order of magnitude
635 deltaymax /=sqrt(10.) ;
636
637 // Emergency break
638 if (deltaymax<_ycutoff) break ;
639 }
640
641 _yatX[_yatXmax] = 1-deltaymaxSave ;
642 _calcX[_yatXmax] = Xsave ;
643
644
645 // Initialize values out of range to 'out-of-range' (-2)
646 for (int i=0 ; i<_yatXmin ; i++) _yatX[i] = -2 ;
647 for (int i=_yatXmax+1 ; i<nbins; i++) _yatX[i] = -2 ;
648 oocxcoutD(_self,Eval) << "RooIntegralMorph::findRange(" << _self->GetName() << "): ymin = " << _yatX[_yatXmin] << " ymax = " << _yatX[_yatXmax] << std::endl;
649 oocxcoutD(_self,Eval) << "RooIntegralMorph::findRange(" << _self->GetName() << "): xmin = " << _calcX[_yatXmin] << " xmax = " << _calcX[_yatXmax] << std::endl;
650}
651
652////////////////////////////////////////////////////////////////////////////////
653/// Dummy
654
656{
657 return 0 ;
658}
659
660////////////////////////////////////////////////////////////////////////////////
661/// Indicate to the RooAbsCachedPdf base class that for the filling of the
662/// cache the traversal of the x should be in the innermost loop, to minimize
663/// recalculation of the one-dimensional internal cache for a fixed value of alpha
664
666{
667 // Put x last to minimize cache faulting
668 orderedObs.removeAll() ;
669
670 orderedObs.add(obs) ;
671 RooAbsArg* obsX = obs.find(x.arg().GetName()) ;
672 if (obsX) {
673 orderedObs.remove(*obsX) ;
674 orderedObs.add(*obsX) ;
675 }
676}
#define e(i)
Definition RSha256.hxx:103
#define coutP(a)
#define oocoutW(o, a)
#define oocxcoutD(o, a)
#define ccoutP(a)
#define oocoutE(o, a)
int Int_t
Definition RtypesCore.h:45
#define X(type, name)
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Option_t Option_t TPoint TPoint const char x2
Option_t Option_t TPoint TPoint const char x1
char name[80]
Definition TGX11.cxx:110
float xmin
float ymin
float xmax
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:77
RooFit::OwningPtr< 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...
Abstract base class for p.d.f.s that need or want to cache their evaluate() output in a RooHistPdf de...
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
RooAbsArg * find(const char *name) const
Find object with given name in list.
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
RooFit::OwningPtr< 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.
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
virtual double offset() const
Definition RooAbsReal.h:373
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
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.
Definition RooDataHist.h:82
void setUnitNorm(bool flag)
Definition RooHistPdf.h:78
std::unique_ptr< RooBrentRootFinder > _rf1
std::unique_ptr< RooAbsReal > _c1
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...
std::unique_ptr< RooAbsFunc > _cb2
void fillGap(Int_t ixlo, Int_t ixhi, double splitPoint=0.5)
Fill all empty histogram bins between bins ixlo and ixhi.
std::unique_ptr< RooArgSet > _nset
void findRange()
Determine which range of y values can be mapped to x values from the numeric inversion of the input c...
std::unique_ptr< RooAbsFunc > _cb1
std::unique_ptr< RooBrentRootFinder > _rf2
RooArgList containedArgs(Action) override
Return all RooAbsArg components contained in this cache.
std::unique_ptr< RooAbsReal > _c2
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...
RooIntegralMorph()=default
RooFit::OwningPtr< RooArgSet > actualObservables(const RooArgSet &nset) const override
Observable to be cached for given choice of normalization.
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.
double evaluate() const override
Dummy.
RooFit::OwningPtr< RooArgSet > actualParameters(const RooArgSet &nset) const override
Parameters of the cache.
Variable that can be changed from the outside.
Definition RooRealVar.h:37
const T & arg() const
Return reference to object held in proxy.
Iterator abstract base class.
Definition TIterator.h:30
const char * GetName() const override
Returns name of object.
Definition TNamed.h:49
Basic string class.
Definition TString.h:139
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
T * OwningPtr
An alias for raw pointers for indicating that the return type of a RooFit function is an owning point...
Definition Config.h:35
OwningPtr< T > makeOwningPtr(std::unique_ptr< T > &&ptr)
Internal helper to turn a std::unique_ptr<T> into an OwningPtr.
Definition Config.h:40