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 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,true) ;
156 par1->remove(x.arg(),true,true) ;
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 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(true) ;
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 x1,x2 ;
298
299 double xmax = _x->getMax("cache") ;
300 double xmin = _x->getMin("cache") ;
301
302 ok=true ;
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 xmax = _x->getMax("cache") ;
317 double 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 xsave = _self->x ;
328
329 // if (!_yatX) {
330 // _yatX = new double[_x->numBins("cache")+1] ;
331 // _calcX = new double[_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 offset = _yatX[_yatXmin] ;
357 double delta = (_yatX[_yatXmax] - _yatX[_yatXmin])/10 ;
358 double y = offset + i*delta ;
359
360 // Calculate corresponding X
361 bool ok ;
362 double 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 xmax = _x->getMax("cache") ;
387 double xmin = _x->getMin("cache") ;
388 double binw = (xmax-xmin)/_x->numBins("cache") ;
389 for (int i=_yatXmin+1 ; i<_yatXmax-1 ; i++) {
390
391 // Calculate additional offset to apply if bin ixlo does not have X value calculated at bin center
392 double xBinC = xmin + (i+0.5)*binw ;
393 double xOffset = xBinC-_calcX[i] ;
394 if (fabs(xOffset/binw)>1e-3) {
395 double slope = (_yatX[i+1]-_yatX[i-1])/(_calcX[i+1]-_calcX[i-1]) ;
396 double 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
409 double x1 = _x->getMin("cache");
410 double x2 = _x->getMin("cache");
411
412 double xMax = _x->getMax("cache");
413
414 // Transfer calculated values to histogram
415 for (int i=_yatXmin ; i<_yatXmax ; i++) {
416
417 double y = _yatX[i] ;
418
419 // Little optimization here exploiting the fact that th cumulative
420 // distribution functions increase monotonically, so we already know that
421 // the next x-value must be higher than the last one as y is increasing. So
422 // we can use the previous x values as lower bounds.
423 _rf1->findRoot(x1,x1,xMax,y) ;
424 _rf2->findRoot(x2,x2,xMax,y) ;
425
426 _x->setVal(x1) ; double f1x1 = _pdf1->getVal(&nsetTmp) ;
427 _x->setVal(x2) ; double f2x2 = _pdf2->getVal(&nsetTmp) ;
428 double fbarX = f1x1*f2x2 / ( _alpha->getVal()*f2x2 + (1-_alpha->getVal())*f1x1 ) ;
429
430 dIter->Next() ;
431 //_hist->get(i) ;
432 hist()->set(fbarX) ;
433 }
434 // Zero output histogram above highest calculable X value
435 for (int i=_yatXmax+1 ; i<nbins ; i++) {
436 dIter->Next() ;
437 //_hist->get(i) ;
438 hist()->set(0) ;
439 }
440
441 pdf()->setUnitNorm(true) ;
442 _self->x = xsave ;
443
444 oocxcoutD(_self,Eval) << "RooIntegralMorph::MorphCacheElem::calculate(" << _self->GetName() << ") calculation required " << _ccounter << " samplings of cdfs" << endl ;
445}
446
447////////////////////////////////////////////////////////////////////////////////
448/// Fill all empty histogram bins between bins ixlo and ixhi. The value of 'splitPoint'
449/// defines the split point for the recursive division strategy to fill the gaps
450/// If the midpoint value of y is very close to the midpoint in x, use interpolation
451/// to fill the gaps, otherwise the intervals again.
452
453void RooIntegralMorph::MorphCacheElem::fillGap(Int_t ixlo, Int_t ixhi, double splitPoint)
454{
455 // CONVENTION: _yatX[ixlo] is filled, _yatX[ixhi] is filled, elements in between are empty
456 // cout << "fillGap: gap from _yatX[" << ixlo << "]=" << _yatX[ixlo] << " to _yatX[" << ixhi << "]=" << _yatX[ixhi] << ", size = " << ixhi-ixlo << endl ;
457
458 if (_yatX[ixlo]<0) {
459 oocoutE(_self,Eval) << "RooIntegralMorph::MorphCacheElme::fillGap(" << _self->GetName() << "): ERROR in fillgap " << ixlo << " = " << ixhi
460 << " splitPoint= " << splitPoint << " _yatX[ixlo] = " << _yatX[ixlo] << endl ;
461 }
462 if (_yatX[ixhi]<0) {
463 oocoutE(_self,Eval) << "RooIntegralMorph::MorphCacheElme::fillGap(" << _self->GetName() << "): ERROR in fillgap " << ixlo << " = " << ixhi
464 << " splitPoint " << splitPoint << " _yatX[ixhi] = " << _yatX[ixhi] << endl ;
465 }
466
467 // Determine where half-way Y value lands
468 double ymid = _yatX[ixlo]*splitPoint + _yatX[ixhi]*(1-splitPoint) ;
469 bool ok ;
470 double Xmid = calcX(ymid,ok) ;
471 if (!ok) {
472 oocoutW(_self,Eval) << "RooIntegralMorph::MorphCacheElem::fillGap(" << _self->GetName() << ") unable to calculate midpoint in gap ["
473 << ixlo << "," << ixhi << "], resorting to interpolation" << endl ;
474 interpolateGap(ixlo,ixhi) ;
475 }
476
477 Int_t iX = binX(Xmid) ;
478 double cq = (Xmid-_calcX[ixlo])/(_calcX[ixhi]-_calcX[ixlo])-0.5 ;
479
480 // Store midway point
481 _yatX[iX] = ymid ;
482 _calcX[iX] = Xmid ;
483
484
485 // Policy: If centration quality is better than 1% OR better than 1/10 of a bin, fill interval with linear interpolation
486 if (fabs(cq)<0.01 || fabs(cq*(ixhi-ixlo))<0.1 || ymid<_ycutoff ) {
487
488 // Fill remaining gaps on either side with linear interpolation
489 if (iX-ixlo>1) {
490 interpolateGap(ixlo,iX) ;
491 }
492 if (ixhi-iX>1) {
493 interpolateGap(iX,ixhi) ;
494 }
495
496 } else {
497
498 if (iX==ixlo) {
499
500 if (splitPoint<0.95) {
501 // Midway value lands on lowest bin, retry split with higher split point
502 double newSplit = splitPoint + 0.5*(1-splitPoint) ;
503 fillGap(ixlo,ixhi,newSplit) ;
504 } else {
505 // Give up and resort to interpolation
506 interpolateGap(ixlo,ixhi) ;
507 }
508
509 } else if (iX==ixhi) {
510
511 // Midway value lands on highest bin, retry split with lower split point
512 if (splitPoint>0.05) {
513 double newSplit = splitPoint/2 ;
514 fillGap(ixlo,ixhi,newSplit) ;
515 } else {
516 // Give up and resort to interpolation
517 interpolateGap(ixlo,ixhi) ;
518 }
519
520 } else {
521
522 // Midway point reasonable, iterate on interval on both sides
523 if (iX-ixlo>1) {
524 fillGap(ixlo,iX) ;
525 }
526 if (ixhi-iX>1) {
527 fillGap(iX,ixhi) ;
528 }
529 }
530 }
531}
532
533////////////////////////////////////////////////////////////////////////////////
534/// Fill empty histogram bins between ixlo and ixhi with values obtained
535/// from linear interpolation of ixlo,ixhi elements.
536
538{
539 //cout << "filling gap with linear interpolation ixlo=" << ixlo << " ixhi=" << ixhi << endl ;
540
541 double xmax = _x->getMax("cache") ;
542 double xmin = _x->getMin("cache") ;
543 double binw = (xmax-xmin)/_x->numBins("cache") ;
544
545 // Calculate deltaY in terms of actual X difference calculate, not based on nominal bin width
546 double deltaY = (_yatX[ixhi]-_yatX[ixlo])/((_calcX[ixhi]-_calcX[ixlo])/binw) ;
547
548 // Calculate additional offset to apply if bin ixlo does not have X value calculated at bin center
549 double xBinC = xmin + (ixlo+0.5)*binw ;
550 double xOffset = xBinC-_calcX[ixlo] ;
551
552 for (int j=ixlo+1 ; j<ixhi ; j++) {
553 _yatX[j] = _yatX[ixlo]+(xOffset+(j-ixlo))*deltaY ;
554 _calcX[j] = xmin + (j+0.5)*binw ;
555 }
556
557}
558
559////////////////////////////////////////////////////////////////////////////////
560/// Determine which range of y values can be mapped to x values
561/// from the numeric inversion of the input c.d.fs.
562/// Start with a y range of [0.1-0.9] and push boundaries
563/// outward with a factor of 1/sqrt(10). Stop iteration if
564/// inverted x values no longer change
565
567{
568 double xmin = _x->getMin("cache") ;
569 double xmax = _x->getMax("cache") ;
570 Int_t nbins = _x->numBins("cache") ;
571
572 double x1,x2 ;
573 bool ok = true ;
574 double ymin=0.1,yminSave(-1) ;
575 double Xsave(-1),Xlast=xmax ;
576
577 // Find lowest Y value that can be measured
578 // Start at 0.1 and iteratively lower limit by sqrt(10)
579 while(true) {
580 ok &= _rf1->findRoot(x1,xmin,xmax,ymin) ;
581 ok &= _rf2->findRoot(x2,xmin,xmax,ymin) ;
582 oocxcoutD(_self,Eval) << "RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName() << ") findMin: x1 = " << x1 << " x2 = " << x2 << " ok = " << (ok?"T":"F") << endl ;
583
584 // Terminate in case of non-convergence
585 if (!ok) break ;
586
587 // Terminate if value of X no longer moves by >0.1 bin size
588 double X = _alpha->getVal()*x1 + (1-_alpha->getVal())*x2 ;
589 if (fabs(X-Xlast)/(xmax-xmin)<0.0001) {
590 break ;
591 }
592 Xlast=X ;
593
594 // Store new Y value
595 _yatXmin = (Int_t)(nbins*(X-xmin)/(xmax-xmin)) ;
596 _yatX[_yatXmin] = ymin ;
597 _calcX[_yatXmin] = X ;
598 yminSave = ymin ;
599 Xsave=X ;
600
601 // Reduce ymin by half an order of magnitude
602 ymin /=sqrt(10.) ;
603
604 // Emergency break
605 if (ymin<_ycutoff) break ;
606 }
607 _yatX[_yatXmin] = yminSave ;
608 _calcX[_yatXmin] = Xsave ;
609
610 // Find highest Y value that can be measured
611 // Start at 1 - 0.1 and iteratively lower delta by sqrt(10)
612 ok = true ;
613 double deltaymax=0.1, deltaymaxSave(-1) ;
614 Xlast=xmin ;
615 while(true) {
616 ok &= _rf1->findRoot(x1,xmin,xmax,1-deltaymax) ;
617 ok &= _rf2->findRoot(x2,xmin,xmax,1-deltaymax) ;
618
619 oocxcoutD(_self,Eval) << "RooIntegralMorph::MorphCacheElem::findRange(" << _self->GetName() << ") findMax: x1 = " << x1 << " x2 = " << x2 << " ok = " << (ok?"T":"F") << endl ;
620
621 // Terminate in case of non-convergence
622 if (!ok) break ;
623
624 // Terminate if value of X no longer moves by >0.1 bin size
625 double X = _alpha->getVal()*x1 + (1-_alpha->getVal())*x2 ;
626 if (fabs(X-Xlast)/(xmax-xmin)<0.0001) {
627 break ;
628 }
629 Xlast=X ;
630
631 // Store new Y value
632 _yatXmax = (Int_t)(nbins*(X-xmin)/(xmax-xmin)) ;
633 _yatX[_yatXmax] = 1-deltaymax ;
634 _calcX[_yatXmax] = X ;
635 deltaymaxSave = deltaymax ;
636 Xsave=X ;
637
638 // Reduce ymin by half an order of magnitude
639 deltaymax /=sqrt(10.) ;
640
641 // Emergency break
642 if (deltaymax<_ycutoff) break ;
643 }
644
645 _yatX[_yatXmax] = 1-deltaymaxSave ;
646 _calcX[_yatXmax] = Xsave ;
647
648
649 // Initialize values out of range to 'out-of-range' (-2)
650 for (int i=0 ; i<_yatXmin ; i++) _yatX[i] = -2 ;
651 for (int i=_yatXmax+1 ; i<nbins; i++) _yatX[i] = -2 ;
652 oocxcoutD(_self,Eval) << "RooIntegralMorph::findRange(" << _self->GetName() << "): ymin = " << _yatX[_yatXmin] << " ymax = " << _yatX[_yatXmax] << endl;
653 oocxcoutD(_self,Eval) << "RooIntegralMorph::findRange(" << _self->GetName() << "): xmin = " << _calcX[_yatXmin] << " xmax = " << _calcX[_yatXmax] << endl;
654}
655
656////////////////////////////////////////////////////////////////////////////////
657/// Dummy
658
660{
661 return 0 ;
662}
663
664////////////////////////////////////////////////////////////////////////////////
665/// Indicate to the RooAbsCachedPdf base class that for the filling of the
666/// cache the traversal of the x should be in the innermost loop, to minimize
667/// recalculation of the one-dimensional internal cache for a fixed value of alpha
668
670{
671 // Put x last to minimize cache faulting
672 orderedObs.removeAll() ;
673
674 orderedObs.add(obs) ;
675 RooAbsArg* obsX = obs.find(x.arg().GetName()) ;
676 if (obsX) {
677 orderedObs.remove(*obsX) ;
678 orderedObs.add(*obsX) ;
679 }
680}
#define e(i)
Definition: RSha256.hxx:103
#define coutP(a)
Definition: RooMsgService.h:35
#define oocoutW(o, a)
Definition: RooMsgService.h:51
#define oocxcoutD(o, a)
Definition: RooMsgService.h:87
#define ccoutP(a)
Definition: RooMsgService.h:43
#define oocoutE(o, a)
Definition: RooMsgService.h:52
int Int_t
Definition: RtypesCore.h:45
#define ClassImp(name)
Definition: Rtypes.h:375
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
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:72
friend class RooArgSet
Definition: RooAbsArg.h:645
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:565
RooArgList containedArgs(Action) override
Returns all RooAbsArg objects contained in the cache element.
RooAbsCachedPdf is the abstract base class for p.d.f.s that need or want to cache their evaluate() ou...
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
virtual bool remove(const RooAbsArg &var, bool silent=false, bool matchByNameOnly=false)
Remove the specified argument from our list.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
void assign(const RooAbsCollection &other) const
Sets the value, cache and constant attribute of any argument in our set that also appears in the othe...
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsReal * createCdf(const RooArgSet &iset, const RooArgSet &nset=RooArgSet())
Create a cumulative distribution function of this p.d.f in terms of the observables listed in iset.
Definition: RooAbsPdf.cxx:3300
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:64
RooAbsFunc * bindVars(const RooArgSet &vars, const RooArgSet *nset=0, bool clipInvalid=false) const
Create an interface adaptor f(vars) that binds us to the specified variables (in arbitrary order).
virtual double offset() const
Definition: RooAbsReal.h:373
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:40
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:57
Implement the abstract 1-dimensional root finding interface using the Brent-Decker method.
void setTol(double tol)
Set convergence tolerance parameter.
TIterator * sliceIterator(RooAbsArg &sliceArg, const RooArgSet &otherArgs)
Create an iterator over all bins in a slice defined by the subset of observables listed in sliceArg.
const RooArgSet * get() const override
Get bin centre of current bin.
Definition: RooDataHist.h:82
void setUnitNorm(bool flag)
Definition: RooHistPdf.h:82
~MorphCacheElem() override
Destructor.
void calculate(TIterator *iter)
Calculate shape of p.d.f for x,alpha values defined by dIter iterator over cache histogram.
void interpolateGap(Int_t ixlo, Int_t ixhi)
Fill empty histogram bins between ixlo and ixhi with values obtained from linear interpolation of ixl...
MorphCacheElem(RooIntegralMorph &self, const RooArgSet *nset)
Construct of cache element, copy relevant input from RooIntegralMorph, create the cdfs from the input...
void fillGap(Int_t ixlo, Int_t ixhi, double splitPoint=0.5)
Fill all empty histogram bins between bins ixlo and ixhi.
void findRange()
Determine which range of y values can be mapped to x values from the numeric inversion of the input c...
RooArgList containedArgs(Action) override
Return all RooAbsArg components contained in this cache.
double calcX(double y, bool &ok)
Calculate the x value of the output p.d.f at the given cdf value y.
Int_t binX(double x)
Return the bin number enclosing the given x value.
Class RooIntegralMorph is an implementation of the histogram interpolation technique described by Ale...
friend class MorphCacheElem
PdfCacheElem * createCache(const RooArgSet *nset) const override
Create and return a derived MorphCacheElem.
const char * inputBaseName() const override
Return base name component for cache components in this case a string encoding the names of both end ...
void preferredObservableScanOrder(const RooArgSet &obs, RooArgSet &orderedObs) const override
Indicate to the RooAbsCachedPdf base class that for the filling of the cache the traversal of the x s...
void fillCacheObject(PdfCacheElem &cache) const override
Fill the cache with the interpolated shape.
RooArgSet * actualParameters(const RooArgSet &nset) const override
Parameters of the cache.
RooRealProxy alpha
double evaluate() const override
Dummy.
RooArgSet * actualObservables(const RooArgSet &nset) const override
Observable to be cached for given choice of normalization.
RooRealVar represents a variable that can be changed from the outside.
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: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)