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 *
6  * *
7  * Redistribution and use in source and binary forms, *
8  * with or without modification, are permitted according to the terms *
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) {
142  }
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()) ;
156  par1->remove(x.arg(),kTRUE,kTRUE) ;
157  if (!_cacheAlpha) {
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 ;
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
669  RooAbsArg* obsX = obs.find(x.arg().GetName()) ;
670  if (obsX) {
671  orderedObs.remove(*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.