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