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