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 _cache(0)
114{
115}
116
117////////////////////////////////////////////////////////////////////////////////
118/// Copy constructor
119
121 RooAbsCachedPdf(other,name),
122 pdf1("pdf1",this,other.pdf1),
123 pdf2("pdf2",this,other.pdf2),
124 x("x",this,other.x),
125 alpha("alpha",this,other.alpha),
126 _cacheAlpha(other._cacheAlpha),
127 _cache(0)
128{
129}
130
131////////////////////////////////////////////////////////////////////////////////
132/// Observable to be cached for given choice of normalization.
133/// Returns the 'x' observable unless doCacheAlpha is set in which
134/// case a set with both x and alpha
135
137{
138 RooArgSet* obs = new RooArgSet ;
139 if (_cacheAlpha) {
140 obs->add(alpha.arg()) ;
141 }
142 obs->add(x.arg()) ;
144}
145
146////////////////////////////////////////////////////////////////////////////////
147/// Parameters of the cache. Returns parameters of both pdf1 and pdf2
148/// and parameter cache, in case doCacheAlpha is not set.
149
151{
152 auto par1 = pdf1->getParameters(static_cast<RooArgSet*>(nullptr));
153 RooArgSet par2;
154 pdf2->getParameters(nullptr, par2);
155 par1->add(par2,true) ;
156 par1->remove(x.arg(),true,true) ;
157 if (!_cacheAlpha) {
158 par1->add(alpha.arg()) ;
159 }
160 return RooFit::OwningPtr<RooArgSet>{std::move(par1)};
161}
162
163////////////////////////////////////////////////////////////////////////////////
164/// Return base name component for cache components in this case
165/// a string encoding the names of both end point p.d.f.s
166
168{
169 static TString name ;
170
171 name = pdf1.arg().GetName() ;
172 name.Append("_MORPH_") ;
173 name.Append(pdf2.arg().GetName()) ;
174 return name.Data() ;
175}
176
177////////////////////////////////////////////////////////////////////////////////
178/// Fill the cache with the interpolated shape.
179
181{
182 MorphCacheElem& mcache = static_cast<MorphCacheElem&>(cache) ;
183
184 // If cacheAlpha is true employ slice iterator here to fill all slices
185
186 if (!_cacheAlpha) {
187
188 std::unique_ptr<TIterator> dIter{cache.hist()->sliceIterator(const_cast<RooAbsReal&>(x.arg()),RooArgSet())};
189 mcache.calculate(dIter.get());
190
191 } else {
192 std::unique_ptr<TIterator> slIter{cache.hist()->sliceIterator(const_cast<RooAbsReal&>(alpha.arg()),RooArgSet())};
193
194 double alphaSave = alpha ;
195 RooArgSet alphaSet(alpha.arg()) ;
196 coutP(Eval) << "RooIntegralMorph::fillCacheObject(" << GetName() << ") filling multi-dimensional cache" ;
197 while(slIter->Next()) {
198 alphaSet.assign(*cache.hist()->get()) ;
199 std::unique_ptr<TIterator> dIter{cache.hist()->sliceIterator(const_cast<RooAbsReal&>(x.arg()),RooArgSet(alpha.arg()))};
200 mcache.calculate(dIter.get());
201 ccoutP(Eval) << "." << flush;
202 }
203 ccoutP(Eval) << std::endl;
204
205 const_cast<RooIntegralMorph*>(this)->alpha = alphaSave ;
206 }
207}
208
209////////////////////////////////////////////////////////////////////////////////
210/// Create and return a derived MorphCacheElem.
211
213{
214 return new MorphCacheElem(const_cast<RooIntegralMorph&>(*this),nset) ;
215}
216
217////////////////////////////////////////////////////////////////////////////////
218/// Return all RooAbsArg components contained in this cache
219
221{
222 RooArgList ret ;
223 ret.add(PdfCacheElem::containedArgs(action)) ;
224 ret.add(*_self) ;
225 ret.add(*_pdf1) ;
226 ret.add(*_pdf2) ;
227 ret.add(*_x ) ;
228 ret.add(*_alpha) ;
229 ret.add(*_c1) ;
230 ret.add(*_c2) ;
231
232 return ret ;
233}
234
235////////////////////////////////////////////////////////////////////////////////
236/// Construct of cache element, copy relevant input from RooIntegralMorph,
237/// create the cdfs from the input p.d.fs and instantiate the root finders
238/// on the cdfs to perform the inversion.
239
241{
242 // Mark in base class that normalization of cached pdf is invariant under pdf parameters
243 _x = (RooRealVar*)self.x.absArg() ;
244 _nset = std::make_unique<RooArgSet>(*_x);
245
246 _alpha = (RooAbsReal*)self.alpha.absArg() ;
247 _pdf1 = (RooAbsPdf*)(self.pdf1.absArg()) ;
248 _pdf2 = (RooAbsPdf*)(self.pdf2.absArg()) ;
249 _c1 = std::unique_ptr<RooAbsReal>{_pdf1->createCdf(*_x)};
250 _c2 = std::unique_ptr<RooAbsReal>{_pdf2->createCdf(*_x)};
251 _cb1 = _c1->bindVars(*_x,_nset.get());
252 _cb2 = _c2->bindVars(*_x,_nset.get());
253 _self = &self ;
254
255 _rf1 = std::make_unique<RooBrentRootFinder>(*_cb1);
256 _rf2 = std::make_unique<RooBrentRootFinder>(*_cb2);
257 _ccounter = 0 ;
258
259 _rf1->setTol(1e-12) ;
260 _rf2->setTol(1e-12) ;
261 _ycutoff = 1e-7 ;
262
263 // _yatX = 0 ;
264 // _calcX = 0 ;
265
266 // Must do this here too: fillCache() may not be called if cache contents is retrieved from EOcache
267 pdf()->setUnitNorm(true) ;
268
269 _yatXmax = 0 ;
270 _yatXmin = 0 ;
271}
272
273////////////////////////////////////////////////////////////////////////////////
274/// Destructor
275
277{
278}
279
280////////////////////////////////////////////////////////////////////////////////
281/// Calculate the x value of the output p.d.f at the given cdf value y.
282/// The ok boolean is filled with the success status of the operation.
283
285{
286 if (y<0 || y>1) {
287 oocoutW(_self,Eval) << "RooIntegralMorph::MorphCacheElem::calcX() WARNING: requested root finding for unphysical CDF value " << y << endl ;
288 }
289 double x1,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,x2 ;
566 bool ok = true ;
567 double ymin=0.1,yminSave(-1) ;
568 double Xsave(-1),Xlast=xmax ;
569
570 // Find lowest Y value that can be measured
571 // Start at 0.1 and iteratively lower limit by sqrt(10)
572 while(true) {
573 ok &= _rf1->findRoot(x1,xmin,xmax,ymin) ;
574 ok &= _rf2->findRoot(x2,xmin,xmax,ymin) ;
575 oocxcoutD(_self,Eval) << "RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName() << ") findMin: x1 = " << x1 << " x2 = " << x2 << " ok = " << (ok?"T":"F") << endl ;
576
577 // Terminate in case of non-convergence
578 if (!ok) break ;
579
580 // Terminate if value of X no longer moves by >0.1 bin size
581 double X = _alpha->getVal()*x1 + (1-_alpha->getVal())*x2 ;
582 if (std::abs(X-Xlast)/(xmax-xmin)<0.0001) {
583 break ;
584 }
585 Xlast=X ;
586
587 // Store new Y value
588 _yatXmin = (Int_t)(nbins*(X-xmin)/(xmax-xmin)) ;
589 _yatX[_yatXmin] = ymin ;
590 _calcX[_yatXmin] = X ;
591 yminSave = ymin ;
592 Xsave=X ;
593
594 // Reduce ymin by half an order of magnitude
595 ymin /=sqrt(10.) ;
596
597 // Emergency break
598 if (ymin<_ycutoff) break ;
599 }
600 _yatX[_yatXmin] = yminSave ;
601 _calcX[_yatXmin] = Xsave ;
602
603 // Find highest Y value that can be measured
604 // Start at 1 - 0.1 and iteratively lower delta by sqrt(10)
605 ok = true ;
606 double deltaymax=0.1, deltaymaxSave(-1) ;
607 Xlast=xmin ;
608 while(true) {
609 ok &= _rf1->findRoot(x1,xmin,xmax,1-deltaymax) ;
610 ok &= _rf2->findRoot(x2,xmin,xmax,1-deltaymax) ;
611
612 oocxcoutD(_self,Eval) << "RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName() << ") findMax: x1 = " << x1 << " x2 = " << x2 << " ok = " << (ok?"T":"F") << endl ;
613
614 // Terminate in case of non-convergence
615 if (!ok) break ;
616
617 // Terminate if value of X no longer moves by >0.1 bin size
618 double X = _alpha->getVal()*x1 + (1-_alpha->getVal())*x2 ;
619 if (std::abs(X-Xlast)/(xmax-xmin)<0.0001) {
620 break ;
621 }
622 Xlast=X ;
623
624 // Store new Y value
625 _yatXmax = (Int_t)(nbins*(X-xmin)/(xmax-xmin)) ;
626 _yatX[_yatXmax] = 1-deltaymax ;
627 _calcX[_yatXmax] = X ;
628 deltaymaxSave = deltaymax ;
629 Xsave=X ;
630
631 // Reduce ymin by half an order of magnitude
632 deltaymax /=sqrt(10.) ;
633
634 // Emergency break
635 if (deltaymax<_ycutoff) break ;
636 }
637
638 _yatX[_yatXmax] = 1-deltaymaxSave ;
639 _calcX[_yatXmax] = Xsave ;
640
641
642 // Initialize values out of range to 'out-of-range' (-2)
643 for (int i=0 ; i<_yatXmin ; i++) _yatX[i] = -2 ;
644 for (int i=_yatXmax+1 ; i<nbins; i++) _yatX[i] = -2 ;
645 oocxcoutD(_self,Eval) << "RooIntegralMorph::findRange(" << _self->GetName() << "): ymin = " << _yatX[_yatXmin] << " ymax = " << _yatX[_yatXmax] << endl;
646 oocxcoutD(_self,Eval) << "RooIntegralMorph::findRange(" << _self->GetName() << "): xmin = " << _calcX[_yatXmin] << " xmax = " << _calcX[_yatXmax] << endl;
647}
648
649////////////////////////////////////////////////////////////////////////////////
650/// Dummy
651
653{
654 return 0 ;
655}
656
657////////////////////////////////////////////////////////////////////////////////
658/// Indicate to the RooAbsCachedPdf base class that for the filling of the
659/// cache the traversal of the x should be in the innermost loop, to minimize
660/// recalculation of the one-dimensional internal cache for a fixed value of alpha
661
663{
664 // Put x last to minimize cache faulting
665 orderedObs.removeAll() ;
666
667 orderedObs.add(obs) ;
668 RooAbsArg* obsX = obs.find(x.arg().GetName()) ;
669 if (obsX) {
670 orderedObs.remove(*obsX) ;
671 orderedObs.add(*obsX) ;
672 }
673}
#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
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition RooAbsArg.h:78
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...
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.
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.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition RooAbsReal.h:61
virtual double offset() const
Definition RooAbsReal.h:384
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooAbsArg * absArg() const
Return pointer to contained argument.
Definition RooArgProxy.h:47
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:79
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...
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< 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...
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.
RooRealVar represents a variable that can be changed from the outside.
Definition RooRealVar.h:40
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
T * OwningPtr
An alias for raw pointers for indicating that the return type of a RooFit function is an owning point...
Definition Config.h:43