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
91using namespace std;
92
94
95////////////////////////////////////////////////////////////////////////////////
96/// Constructor with observables x, pdf shapes pdf1 and pdf2 which represent
97/// the shapes at the end points of the interpolation parameter alpha
98/// If doCacheAlpha is true, a two-dimensional cache is constructed in
99/// both alpha and x
100
101RooIntegralMorph::RooIntegralMorph(const char *name, const char *title,
102 RooAbsReal& _pdf1,
103 RooAbsReal& _pdf2,
104 RooAbsReal& _x,
105 RooAbsReal& _alpha,
106 bool doCacheAlpha) :
107 RooAbsCachedPdf(name,title,2),
108 pdf1("pdf1","pdf1",this,_pdf1),
109 pdf2("pdf2","pdf2",this,_pdf2),
110 x("x","x",this,_x),
111 alpha("alpha","alpha",this,_alpha),
112 _cacheAlpha(doCacheAlpha)
113{
114}
115
116////////////////////////////////////////////////////////////////////////////////
117/// Copy constructor
118
120 RooAbsCachedPdf(other,name),
121 pdf1("pdf1",this,other.pdf1),
122 pdf2("pdf2",this,other.pdf2),
123 x("x",this,other.x),
124 alpha("alpha",this,other.alpha),
125 _cacheAlpha(other._cacheAlpha)
126{
127}
128
129////////////////////////////////////////////////////////////////////////////////
130/// Observable to be cached for given choice of normalization.
131/// Returns the 'x' observable unless doCacheAlpha is set in which
132/// case a set with both x and alpha
133
135{
136 RooArgSet* obs = new RooArgSet ;
137 if (_cacheAlpha) {
138 obs->add(alpha.arg()) ;
139 }
140 obs->add(x.arg()) ;
142}
143
144////////////////////////////////////////////////////////////////////////////////
145/// Parameters of the cache. Returns parameters of both pdf1 and pdf2
146/// and parameter cache, in case doCacheAlpha is not set.
147
149{
150 std::unique_ptr<RooArgSet> par1{pdf1->getParameters(static_cast<RooArgSet*>(nullptr))};
151 RooArgSet par2;
152 pdf2->getParameters(nullptr, par2);
153 par1->add(par2,true) ;
154 par1->remove(x.arg(),true,true) ;
155 if (!_cacheAlpha) {
156 par1->add(alpha.arg()) ;
157 }
158 return RooFit::makeOwningPtr(std::move(par1));
159}
160
161////////////////////////////////////////////////////////////////////////////////
162/// Return base name component for cache components in this case
163/// a string encoding the names of both end point p.d.f.s
164
166{
167 static TString name ;
168
169 name = pdf1.arg().GetName() ;
170 name.Append("_MORPH_") ;
171 name.Append(pdf2.arg().GetName()) ;
172 return name.Data() ;
173}
174
175////////////////////////////////////////////////////////////////////////////////
176/// Fill the cache with the interpolated shape.
177
179{
180 MorphCacheElem& mcache = static_cast<MorphCacheElem&>(cache) ;
181
182 // If cacheAlpha is true employ slice iterator here to fill all slices
183
184 if (!_cacheAlpha) {
185
186 std::unique_ptr<TIterator> dIter{cache.hist()->sliceIterator(const_cast<RooAbsReal&>(x.arg()),RooArgSet())};
187 mcache.calculate(dIter.get());
188
189 } else {
190 std::unique_ptr<TIterator> slIter{cache.hist()->sliceIterator(const_cast<RooAbsReal&>(alpha.arg()),RooArgSet())};
191
192 double alphaSave = alpha ;
193 RooArgSet alphaSet(alpha.arg()) ;
194 coutP(Eval) << "RooIntegralMorph::fillCacheObject(" << GetName() << ") filling multi-dimensional cache" ;
195 while(slIter->Next()) {
196 alphaSet.assign(*cache.hist()->get()) ;
197 std::unique_ptr<TIterator> dIter{cache.hist()->sliceIterator(const_cast<RooAbsReal&>(x.arg()),RooArgSet(alpha.arg()))};
198 mcache.calculate(dIter.get());
199 ccoutP(Eval) << "." << flush;
200 }
201 ccoutP(Eval) << std::endl;
202
203 const_cast<RooIntegralMorph*>(this)->alpha = alphaSave ;
204 }
205}
206
207////////////////////////////////////////////////////////////////////////////////
208/// Create and return a derived MorphCacheElem.
209
211{
212 return new MorphCacheElem(const_cast<RooIntegralMorph&>(*this),nset) ;
213}
214
215////////////////////////////////////////////////////////////////////////////////
216/// Return all RooAbsArg components contained in this cache
217
219{
220 RooArgList ret ;
221 ret.add(PdfCacheElem::containedArgs(action)) ;
222 ret.add(*_self) ;
223 ret.add(*_pdf1) ;
224 ret.add(*_pdf2) ;
225 ret.add(*_x ) ;
226 ret.add(*_alpha) ;
227 ret.add(*_c1) ;
228 ret.add(*_c2) ;
229
230 return ret ;
231}
232
233////////////////////////////////////////////////////////////////////////////////
234/// Construct of cache element, copy relevant input from RooIntegralMorph,
235/// create the cdfs from the input p.d.fs and instantiate the root finders
236/// on the cdfs to perform the inversion.
237
239 : PdfCacheElem(self, nsetIn),
240 _self(&self),
241 _pdf1(static_cast<RooAbsPdf *>(self.pdf1.absArg())),
242 _pdf2(static_cast<RooAbsPdf *>(self.pdf2.absArg())),
243 _x(static_cast<RooRealVar *>(self.x.absArg())),
244 _alpha(static_cast<RooAbsReal *>(self.alpha.absArg())),
245 _yatXmin(0),
246 _yatXmax(0),
247 _ccounter(0),
248 _ycutoff(1e-7)
249{
250 // Mark in base class that normalization of cached pdf is invariant under pdf parameters
251
252 _nset = std::make_unique<RooArgSet>(*_x);
253
254 _c1 = std::unique_ptr<RooAbsReal>{_pdf1->createCdf(*_x)};
255 _c2 = std::unique_ptr<RooAbsReal>{_pdf2->createCdf(*_x)};
256 _cb1 = std::unique_ptr<RooAbsFunc>{_c1->bindVars(*_x,_nset.get())};
257 _cb2 = std::unique_ptr<RooAbsFunc>{_c2->bindVars(*_x,_nset.get())};
258
259 _rf1 = std::make_unique<RooBrentRootFinder>(*_cb1);
260 _rf2 = std::make_unique<RooBrentRootFinder>(*_cb2);
261
262 _rf1->setTol(1e-12) ;
263 _rf2->setTol(1e-12) ;
264
265 // _yatX = 0 ;
266 // _calcX = 0 ;
267
268 // Must do this here too: fillCache() may not be called if cache contents is retrieved from EOcache
269 pdf()->setUnitNorm(true) ;
270}
271
272////////////////////////////////////////////////////////////////////////////////
273/// Destructor
274
276{
277}
278
279////////////////////////////////////////////////////////////////////////////////
280/// Calculate the x value of the output p.d.f at the given cdf value y.
281/// The ok boolean is filled with the success status of the operation.
282
284{
285 if (y<0 || y>1) {
286 oocoutW(_self,Eval) << "RooIntegralMorph::MorphCacheElem::calcX() WARNING: requested root finding for unphysical CDF value " << y << endl ;
287 }
288 double x1;
289 double x2;
290
291 double xmax = _x->getMax("cache") ;
292 double xmin = _x->getMin("cache") ;
293
294 ok=true ;
295 ok &= _rf1->findRoot(x1,xmin,xmax,y) ;
296 ok &= _rf2->findRoot(x2,xmin,xmax,y) ;
297 if (!ok) return 0 ;
298 _ccounter++ ;
299
300 return _alpha->getVal()*x1 + (1-_alpha->getVal())*x2 ;
301}
302
303////////////////////////////////////////////////////////////////////////////////
304/// Return the bin number enclosing the given x value
305
307{
308 double xmax = _x->getMax("cache") ;
309 double xmin = _x->getMin("cache") ;
310 return (Int_t)(_x->numBins("cache")*(X-xmin)/(xmax-xmin)) ;
311}
312
313////////////////////////////////////////////////////////////////////////////////
314/// Calculate shape of p.d.f for x,alpha values
315/// defined by dIter iterator over cache histogram
316
318{
319 double xsave = _self->x ;
320
321 // if (!_yatX) {
322 // _yatX = new double[_x->numBins("cache")+1] ;
323 // _calcX = new double[_x->numBins("cache")+1] ;
324 // }
325
326 _yatX.resize(_x->numBins("cache")+1);
327 _calcX.resize(_x->numBins("cache")+1);
328
329 _ccounter = 0 ;
330
331 // Get number of bins from PdfCacheElem histogram
332 Int_t nbins = _x->numBins("cache") ;
333
334 // Initialize yatX array to 'un-calculated values (-1)'
335 for (int i=0 ; i<nbins ; i++) {
336 _yatX[i] = -1 ;
337 _calcX[i] = 0 ;
338 }
339
340 // Find low and high point
341 findRange() ;
342
343 // Perform initial scan of 100 points
344 for (int i=0 ; i<10 ; i++) {
345
346 // Take a point in y
347 double offset = _yatX[_yatXmin] ;
348 double delta = (_yatX[_yatXmax] - _yatX[_yatXmin])/10 ;
349 double y = offset + i*delta ;
350
351 // Calculate corresponding X
352 bool ok ;
353 double X = calcX(y,ok) ;
354 if (ok) {
355 Int_t iX = binX(X) ;
356 _yatX[iX] = y ;
357 _calcX[iX] = X ;
358 }
359 }
360
361 // Now take an iteration filling the 'gaps'
362 Int_t igapLow = _yatXmin+1 ;
363 while(true) {
364 // Find next gap
365 Int_t igapHigh = igapLow+1 ;
366 while(igapHigh<(_yatXmax) && _yatX[igapHigh]<0) igapHigh++ ;
367
368 // Fill the gap (iteratively and/or using interpolation)
369 fillGap(igapLow-1,igapHigh) ;
370
371 // Terminate after processing of last gap
372 if (igapHigh>=_yatXmax-1) break ;
373 igapLow = igapHigh+1 ;
374 }
375
376 // Make one more iteration to recalculate Y value at bin centers
377 double xmax = _x->getMax("cache") ;
378 double xmin = _x->getMin("cache") ;
379 double binw = (xmax-xmin)/_x->numBins("cache") ;
380 for (int i=_yatXmin+1 ; i<_yatXmax-1 ; i++) {
381
382 // Calculate additional offset to apply if bin ixlo does not have X value calculated at bin center
383 double xBinC = xmin + (i+0.5)*binw ;
384 double xOffset = xBinC-_calcX[i] ;
385 if (std::abs(xOffset/binw)>1e-3) {
386 double slope = (_yatX[i+1]-_yatX[i-1])/(_calcX[i+1]-_calcX[i-1]) ;
387 double newY = _yatX[i] + slope*xOffset ;
388 //cout << "bin " << i << " needs to be re-centered " << xOffset/binw << " slope = " << slope << " origY = " << _yatX[i] << " newY = " << newY << endl ;
389 _yatX[i] = newY ;
390 }
391 }
392
393 // Zero output histogram below lowest calculable X value
394 for (int i=0; i<_yatXmin ; i++) {
395 dIter->Next() ;
396 //_hist->get(i) ;
397 hist()->set(0) ;
398 }
399
400 double x1 = _x->getMin("cache");
401 double x2 = _x->getMin("cache");
402
403 double xMax = _x->getMax("cache");
404
405 // Transfer calculated values to histogram
406 for (int i=_yatXmin ; i<_yatXmax ; i++) {
407
408 double y = _yatX[i] ;
409
410 // Little optimization here exploiting the fact that th cumulative
411 // distribution functions increase monotonically, so we already know that
412 // the next x-value must be higher than the last one as y is increasing. So
413 // we can use the previous x values as lower bounds.
414 _rf1->findRoot(x1,x1,xMax,y) ;
415 _rf2->findRoot(x2,x2,xMax,y) ;
416
417 _x->setVal(x1);
418 double f1x1 = _pdf1->getVal(_nset.get());
419 _x->setVal(x2);
420 double f2x2 = _pdf2->getVal(_nset.get());
421 double fbarX = f1x1*f2x2 / ( _alpha->getVal()*f2x2 + (1-_alpha->getVal())*f1x1 ) ;
422
423 dIter->Next() ;
424 //_hist->get(i) ;
425 hist()->set(fbarX) ;
426 }
427 // Zero output histogram above highest calculable X value
428 for (int i=_yatXmax+1 ; i<nbins ; i++) {
429 dIter->Next() ;
430 //_hist->get(i) ;
431 hist()->set(0) ;
432 }
433
434 pdf()->setUnitNorm(true) ;
435 _self->x = xsave ;
436
437 oocxcoutD(_self,Eval) << "RooIntegralMorph::MorphCacheElem::calculate(" << _self->GetName() << ") calculation required " << _ccounter << " samplings of cdfs" << endl ;
438}
439
440////////////////////////////////////////////////////////////////////////////////
441/// Fill all empty histogram bins between bins ixlo and ixhi. The value of 'splitPoint'
442/// defines the split point for the recursive division strategy to fill the gaps
443/// If the midpoint value of y is very close to the midpoint in x, use interpolation
444/// to fill the gaps, otherwise the intervals again.
445
446void RooIntegralMorph::MorphCacheElem::fillGap(Int_t ixlo, Int_t ixhi, double splitPoint)
447{
448 // CONVENTION: _yatX[ixlo] is filled, _yatX[ixhi] is filled, elements in between are empty
449 // cout << "fillGap: gap from _yatX[" << ixlo << "]=" << _yatX[ixlo] << " to _yatX[" << ixhi << "]=" << _yatX[ixhi] << ", size = " << ixhi-ixlo << endl ;
450
451 if (_yatX[ixlo]<0) {
452 oocoutE(_self,Eval) << "RooIntegralMorph::MorphCacheElme::fillGap(" << _self->GetName() << "): ERROR in fillgap " << ixlo << " = " << ixhi
453 << " splitPoint= " << splitPoint << " _yatX[ixlo] = " << _yatX[ixlo] << endl ;
454 }
455 if (_yatX[ixhi]<0) {
456 oocoutE(_self,Eval) << "RooIntegralMorph::MorphCacheElme::fillGap(" << _self->GetName() << "): ERROR in fillgap " << ixlo << " = " << ixhi
457 << " splitPoint " << splitPoint << " _yatX[ixhi] = " << _yatX[ixhi] << endl ;
458 }
459
460 // Determine where half-way Y value lands
461 double ymid = _yatX[ixlo]*splitPoint + _yatX[ixhi]*(1-splitPoint) ;
462 bool ok ;
463 double Xmid = calcX(ymid,ok) ;
464 if (!ok) {
465 oocoutW(_self,Eval) << "RooIntegralMorph::MorphCacheElem::fillGap(" << _self->GetName() << ") unable to calculate midpoint in gap ["
466 << ixlo << "," << ixhi << "], resorting to interpolation" << endl ;
467 interpolateGap(ixlo,ixhi) ;
468 }
469
470 Int_t iX = binX(Xmid) ;
471 double cq = (Xmid-_calcX[ixlo])/(_calcX[ixhi]-_calcX[ixlo])-0.5 ;
472
473 // Store midway point
474 _yatX[iX] = ymid ;
475 _calcX[iX] = Xmid ;
476
477
478 // Policy: If centration quality is better than 1% OR better than 1/10 of a bin, fill interval with linear interpolation
479 if (std::abs(cq)<0.01 || std::abs(cq*(ixhi-ixlo))<0.1 || ymid<_ycutoff ) {
480
481 // Fill remaining gaps on either side with linear interpolation
482 if (iX-ixlo>1) {
483 interpolateGap(ixlo,iX) ;
484 }
485 if (ixhi-iX>1) {
486 interpolateGap(iX,ixhi) ;
487 }
488
489 } else {
490
491 if (iX==ixlo) {
492
493 if (splitPoint<0.95) {
494 // Midway value lands on lowest bin, retry split with higher split point
495 double newSplit = splitPoint + 0.5*(1-splitPoint) ;
496 fillGap(ixlo,ixhi,newSplit) ;
497 } else {
498 // Give up and resort to interpolation
499 interpolateGap(ixlo,ixhi) ;
500 }
501
502 } else if (iX==ixhi) {
503
504 // Midway value lands on highest bin, retry split with lower split point
505 if (splitPoint>0.05) {
506 double newSplit = splitPoint/2 ;
507 fillGap(ixlo,ixhi,newSplit) ;
508 } else {
509 // Give up and resort to interpolation
510 interpolateGap(ixlo,ixhi) ;
511 }
512
513 } else {
514
515 // Midway point reasonable, iterate on interval on both sides
516 if (iX-ixlo>1) {
517 fillGap(ixlo,iX) ;
518 }
519 if (ixhi-iX>1) {
520 fillGap(iX,ixhi) ;
521 }
522 }
523 }
524}
525
526////////////////////////////////////////////////////////////////////////////////
527/// Fill empty histogram bins between ixlo and ixhi with values obtained
528/// from linear interpolation of ixlo,ixhi elements.
529
531{
532 //cout << "filling gap with linear interpolation ixlo=" << ixlo << " ixhi=" << ixhi << endl ;
533
534 double xmax = _x->getMax("cache") ;
535 double xmin = _x->getMin("cache") ;
536 double binw = (xmax-xmin)/_x->numBins("cache") ;
537
538 // Calculate deltaY in terms of actual X difference calculate, not based on nominal bin width
539 double deltaY = (_yatX[ixhi]-_yatX[ixlo])/((_calcX[ixhi]-_calcX[ixlo])/binw) ;
540
541 // Calculate additional offset to apply if bin ixlo does not have X value calculated at bin center
542 double xBinC = xmin + (ixlo+0.5)*binw ;
543 double xOffset = xBinC-_calcX[ixlo] ;
544
545 for (int j=ixlo+1 ; j<ixhi ; j++) {
546 _yatX[j] = _yatX[ixlo]+(xOffset+(j-ixlo))*deltaY ;
547 _calcX[j] = xmin + (j+0.5)*binw ;
548 }
549
550}
551
552////////////////////////////////////////////////////////////////////////////////
553/// Determine which range of y values can be mapped to x values
554/// from the numeric inversion of the input c.d.fs.
555/// Start with a y range of [0.1-0.9] and push boundaries
556/// outward with a factor of 1/sqrt(10). Stop iteration if
557/// inverted x values no longer change
558
560{
561 double xmin = _x->getMin("cache") ;
562 double xmax = _x->getMax("cache") ;
563 Int_t nbins = _x->numBins("cache") ;
564
565 double x1;
566 double x2;
567 bool ok = true ;
568 double ymin = 0.1;
569 double yminSave(-1);
570 double Xsave(-1);
571 double Xlast = xmax;
572
573 // Find lowest Y value that can be measured
574 // Start at 0.1 and iteratively lower limit by sqrt(10)
575 while(true) {
576 ok &= _rf1->findRoot(x1,xmin,xmax,ymin) ;
577 ok &= _rf2->findRoot(x2,xmin,xmax,ymin) ;
578 oocxcoutD(_self,Eval) << "RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName() << ") findMin: x1 = " << x1 << " x2 = " << x2 << " ok = " << (ok?"T":"F") << endl ;
579
580 // Terminate in case of non-convergence
581 if (!ok) break ;
582
583 // Terminate if value of X no longer moves by >0.1 bin size
584 double X = _alpha->getVal()*x1 + (1-_alpha->getVal())*x2 ;
585 if (std::abs(X-Xlast)/(xmax-xmin)<0.0001) {
586 break ;
587 }
588 Xlast=X ;
589
590 // Store new Y value
591 _yatXmin = (Int_t)(nbins*(X-xmin)/(xmax-xmin)) ;
592 _yatX[_yatXmin] = ymin ;
593 _calcX[_yatXmin] = X ;
594 yminSave = ymin ;
595 Xsave=X ;
596
597 // Reduce ymin by half an order of magnitude
598 ymin /=sqrt(10.) ;
599
600 // Emergency break
601 if (ymin<_ycutoff) break ;
602 }
603 _yatX[_yatXmin] = yminSave ;
604 _calcX[_yatXmin] = Xsave ;
605
606 // Find highest Y value that can be measured
607 // Start at 1 - 0.1 and iteratively lower delta by sqrt(10)
608 ok = true ;
609 double deltaymax = 0.1;
610 double deltaymaxSave(-1);
611 Xlast=xmin ;
612 while(true) {
613 ok &= _rf1->findRoot(x1,xmin,xmax,1-deltaymax) ;
614 ok &= _rf2->findRoot(x2,xmin,xmax,1-deltaymax) ;
615
616 oocxcoutD(_self,Eval) << "RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName() << ") findMax: x1 = " << x1 << " x2 = " << x2 << " ok = " << (ok?"T":"F") << endl ;
617
618 // Terminate in case of non-convergence
619 if (!ok) break ;
620
621 // Terminate if value of X no longer moves by >0.1 bin size
622 double X = _alpha->getVal()*x1 + (1-_alpha->getVal())*x2 ;
623 if (std::abs(X-Xlast)/(xmax-xmin)<0.0001) {
624 break ;
625 }
626 Xlast=X ;
627
628 // Store new Y value
629 _yatXmax = (Int_t)(nbins*(X-xmin)/(xmax-xmin)) ;
630 _yatX[_yatXmax] = 1-deltaymax ;
631 _calcX[_yatXmax] = X ;
632 deltaymaxSave = deltaymax ;
633 Xsave=X ;
634
635 // Reduce ymin by half an order of magnitude
636 deltaymax /=sqrt(10.) ;
637
638 // Emergency break
639 if (deltaymax<_ycutoff) break ;
640 }
641
642 _yatX[_yatXmax] = 1-deltaymaxSave ;
643 _calcX[_yatXmax] = Xsave ;
644
645
646 // Initialize values out of range to 'out-of-range' (-2)
647 for (int i=0 ; i<_yatXmin ; i++) _yatX[i] = -2 ;
648 for (int i=_yatXmax+1 ; i<nbins; i++) _yatX[i] = -2 ;
649 oocxcoutD(_self,Eval) << "RooIntegralMorph::findRange(" << _self->GetName() << "): ymin = " << _yatX[_yatXmin] << " ymax = " << _yatX[_yatXmax] << endl;
650 oocxcoutD(_self,Eval) << "RooIntegralMorph::findRange(" << _self->GetName() << "): xmin = " << _calcX[_yatXmin] << " xmax = " << _calcX[_yatXmax] << endl;
651}
652
653////////////////////////////////////////////////////////////////////////////////
654/// Dummy
655
657{
658 return 0 ;
659}
660
661////////////////////////////////////////////////////////////////////////////////
662/// Indicate to the RooAbsCachedPdf base class that for the filling of the
663/// cache the traversal of the x should be in the innermost loop, to minimize
664/// recalculation of the one-dimensional internal cache for a fixed value of alpha
665
667{
668 // Put x last to minimize cache faulting
669 orderedObs.removeAll() ;
670
671 orderedObs.add(obs) ;
672 RooAbsArg* obsX = obs.find(x.arg().GetName()) ;
673 if (obsX) {
674 orderedObs.remove(*obsX) ;
675 orderedObs.add(*obsX) ;
676 }
677}
#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 ClassImp(name)
Definition Rtypes.h:377
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 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.
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:378
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:55
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:76
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.
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
virtual TObject * Next()=0
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
Basic string class.
Definition TString.h:139
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
OwningPtr< T > makeOwningPtr(std::unique_ptr< T > &&ptr)
Internal helper to turn a std::unique_ptr<T> into an OwningPtr.
Definition Config.h:40
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