Logo ROOT   6.10/09
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 
15 Class RooIntegralMorph is an implementation of the histogram interpolation
16 technique described by Alex Read in 'NIM A 425 (1999) 357-369 'Linear interpolation of histograms'
17 for continuous functions rather than histograms. The interpolation method, in short,
18 works 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 
28 From a technical point of view class RooIntegralMorph is a p.d.f that takes
29 two input p.d.fs f1(x,p) an f2(x,q) and an interpolation parameter to
30 make a p.d.f fbar(x,p,q,alpha). The shapes f1 and f2 are always taken
31 to be end the end-points of the parameter alpha, regardless of what
32 the those numeric values are.
33 
34 Since the value of fbar(x) cannot be easily calculated for a given value
35 of x, class RooIntegralMorph is an implementation of RooAbsCachedPdf and
36 calculates the shape of the interpolated p.d.f. fbar(x) for all values
37 of 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
39 can be controlled by the binning named "cache" on the RooRealVar representing
40 the observable x. The fbar sampling algorithm is based on a recursive division
41 mechanism with a built-in precision cutoff: First an initial sampling in
42 64 equally spaced bins is made. Then the value of fbar is calculated in
43 the center of each gap. If the calculated value deviates too much from
44 the value obtained by linear interpolation from the edge bins, gap
45 is recursively divided. This strategy makes it possible to define a very
46 fine cache sampling (e.g. 1000 or 10000) bins without incurring a
47 corresponding CPU penalty.
48 
49 Note on numeric stability of the algorithm. Since the algorithm relies
50 on a numeric inversion of cumulative distributions functions, some precision
51 may be lost at the 'edges' of the same (i.e. at regions in x where the
52 c.d.f. value is close to zero or one). The general sampling strategy is
53 to start with 64 equally spaces samples in the range y=(0.01-0.99).
54 Then the y ranges are pushed outward by reducing y (or the distance of y to 1.0)
55 by a factor of sqrt(10) iteratively up to the point where the corresponding
56 x value no longer changes significantly. For p.d.f.s with very flat tails
57 such as Gaussians some part of the tail may be lost due to limitations
58 in numeric precision in the CDF inversion step.
59 
60 An effect related to the above limitation in numeric precision should
61 be anticipated when floating the alpha parameter in a fit. If a p.d.f
62 with such flat tails is fitted, it is likely that the dataset contains
63 events in the flat tail region. If the alpha parameter is varied, the
64 likelihood contribution from such events may exhibit discontinuities
65 in alpha, causing discontinuities in the summed likelihood as well
66 that will cause convergence problems in MINUIT. To mitigate this effect
67 one can use the setCacheAlpha() method to instruct RooIntegralMorph
68 to construct a two-dimensional cache for its output values in both
69 x and alpha. If linear interpolation is requested on the resulting
70 output histogram, the resulting interpolation of the p.d.f in the
71 alpha dimension will smooth out the discontinuities in the tail regions
72 result in a continuous likelihood distribution that can be fitted.
73 An added advantage of the cacheAlpha option is that if parameters
74 p,q of f1,f2 are fixed, the cached values in RooIntegralMorph are
75 valid for the entire fit session and do not need to be recalculated
76 for each change in alpha, which may result an considerable increase
77 in 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 
92 using 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 
102 RooIntegralMorph::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 {
153  RooArgSet* par1 = pdf1.arg().getParameters(RooArgSet()) ;
154  RooArgSet* par2 = pdf2.arg().getParameters(RooArgSet()) ;
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 = (*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 
259  _rf1 = new RooBrentRootFinder(*_cb1) ;
260  _rf2 = new RooBrentRootFinder(*_cb2) ;
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
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(_yatX[igapHigh]<0 && igapHigh<(_yatXmax)) 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 }
virtual Double_t getMin(const char *name=0) const
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
float xmin
Definition: THbookFile.cxx:93
virtual Double_t getMax(const char *name=0) const
virtual Bool_t add(const RooAbsCollection &col, Bool_t silent=kFALSE)
Add a collection of arguments to this collection by calling add() for each element in the source coll...
Definition: RooArgSet.h:86
virtual const RooArgSet * get() const
Definition: RooDataHist.h:77
std::vector< Double_t > _calcX
std::vector< Double_t > _yatX
virtual RooArgSet * actualObservables(const RooArgSet &nset) const
Observable to be cached for given choice of normalization.
float ymin
Definition: THbookFile.cxx:93
virtual Int_t numBins(const char *rangeName=0) const
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
void set(Double_t weight, Double_t wgtErr=-1)
Increment the weight of the bin enclosing the coordinates given by &#39;row&#39; by the specified amount...
virtual RooArgList containedArgs(Action)
Returns all RooAbsArg objects contained in the cache element.
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
MorphCacheElem(RooIntegralMorph &self, const RooArgSet *nset)
Construct of cache element, copy relevant input from RooIntegralMorph, create the cdfs from the input...
int nbins[3]
STL namespace.
Float_t delta
Iterator abstract base class.
Definition: TIterator.h:30
#define coutP(a)
Definition: RooMsgService.h:32
virtual RooArgSet * actualParameters(const RooArgSet &nset) const
Parameters of the cache.
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:3027
double sqrt(double)
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
virtual void removeAll()
Remove all arguments from our set, deleting them if we own them.
#define oocoutE(o, a)
Definition: RooMsgService.h:47
Int_t binX(Double_t x)
Return the bin number enclosing the given x value.
#define oocxcoutD(o, a)
Definition: RooMsgService.h:81
friend class RooArgSet
Definition: RooAbsArg.h:469
Class RooIntegralMorph is an implementation of the histogram interpolation technique described by Ale...
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:36
virtual Double_t offset() const
Definition: RooAbsReal.h:309
virtual void setVal(Double_t value)
Set value of variable to &#39;value&#39;.
Definition: RooRealVar.cxx:205
virtual Bool_t findRoot(Double_t &result, Double_t xlo, Double_t xhi, Double_t value=0) const
Do the root finding using the Brent-Decker method.
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
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.
virtual const char * inputBaseName() const
Return base name component for cache components in this case a string encoding the names of both end ...
void findRange()
Determine which range of y values can be mapped to x values from the numeric inversion of the input c...
Implement the abstract 1-dimensional root finding interface using the Brent-Decker method...
RooAbsArg * absArg() const
Definition: RooArgProxy.h:37
void setTol(Double_t tol)
float xmax
Definition: THbookFile.cxx:93
#define ccoutP(a)
Definition: RooMsgService.h:39
virtual void fillCacheObject(PdfCacheElem &cache) const
Fill the cache with the interpolated shape.
void fillGap(Int_t ixlo, Int_t ixhi, Double_t splitPoint=0.5)
Fill all empty histogram bins between bins ixlo and ixhi.
static const double x1[5]
RooAbsCachedPdf is the abstract base class for p.d.f.s that need or want to cache their evaluate() ou...
#define ClassImp(name)
Definition: Rtypes.h:336
RooAbsArg * find(const char *name) const
Find object with given name in list.
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
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)...
RooArgSet * getParameters(const RooAbsData *data, Bool_t stripDisconnected=kTRUE) const
Create a list of leaf nodes in the arg tree starting with ourself as top node that don&#39;t match any of...
Definition: RooAbsArg.cxx:560
RooRealProxy pdf1
Double_t y[n]
Definition: legend1.C:17
RooRealProxy alpha
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...
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
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...
#define oocoutW(o, a)
Definition: RooMsgService.h:46
virtual RooArgList containedArgs(Action)
Return all RooAbsArg components contained in this cache.
Double_t evaluate() const
Dummy.
virtual Bool_t remove(const RooAbsArg &var, Bool_t silent=kFALSE, Bool_t matchByNameOnly=kFALSE)
Remove the specified argument from our list.
void setUnitNorm(Bool_t flag)
Definition: RooHistPdf.h:67
void interpolateGap(Int_t ixlo, Int_t ixhi)
Fill empty histogram bins between ixlo and ixhi with values obtained from linear interpolation of ixl...
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
virtual TObject * Next()=0
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
const RooAbsReal & arg() const
Definition: RooRealProxy.h:43
friend class MorphCacheElem
virtual PdfCacheElem * createCache(const RooArgSet *nset) const
Create and return a derived MorphCacheElem.
const Bool_t kTRUE
Definition: RtypesCore.h:91
void calculate(TIterator *iter)
Calculate shape of p.d.f for x,alpha values defined by dIter iterator over cache histogram.