Logo ROOT  
Reference Guide
PiecewiseInterpolation.cxx
Go to the documentation of this file.
1 /** \class PiecewiseInterpolation
2  * \ingroup HistFactory
3  * The PiecewiseInterpolation is a class that can morph distributions into each other, which
4  * is useful to estimate systematic uncertainties. Given a nominal distribution and one or
5  * more altered or distorted ones, it computes a new shape depending on the value of the nuisance
6  * parameters \f$ \alpha_i \f$:
7  * \f[
8  * A = \sum_i \mathrm{Interpolate}(\mathrm{low}_i, \mathrm{nominal}, \mathrm{high}_i, \alpha_i).
9  * \f]
10  * If an \f$ \alpha_i \f$ is zero, the distribution is identical to the nominal distribution, at
11  * \f$ \pm 1 \f$ it is identical to the up/down distribution for that specific \f$ i \f$.
12  *
13  * The class supports several interpolation methods, which can be selected for each parameter separately
14  * using setInterpCode(). The default interpolation code is 4. This performs
15  * - \f$ |\alpha | > 1 \f$: Linear extrapolation.
16  * - \f$ |\alpha | < 1 \f$: Polynomial interpolation. A sixth-order polynomial is used. Its coefficients
17  * are chosen such that function, first, and second derivative at \f$ \alpha \pm 1 \f$ match the values
18  * that the extrapolation procedure uses.
19  */
20 
22 
23 #include "RooFit.h"
24 
25 #include "Riostream.h"
26 #include "TBuffer.h"
27 
28 #include "RooAbsReal.h"
29 #include "RooAbsPdf.h"
30 #include "RooErrorHandler.h"
31 #include "RooArgSet.h"
32 #include "RooRealVar.h"
33 #include "RooMsgService.h"
34 #include "RooNumIntConfig.h"
35 #include "RooTrace.h"
36 #include "RunContext.h"
37 
38 #include <exception>
39 #include <math.h>
40 #include <algorithm>
41 
42 using namespace std;
43 
45 ;
46 
47 
48 ////////////////////////////////////////////////////////////////////////////////
49 
51 {
52  _positiveDefinite=false;
54 }
55 
56 
57 
58 ////////////////////////////////////////////////////////////////////////////////
59 /// Construct a new interpolation. The value of the function will be
60 /// \f[
61 /// A = \sum_i \mathrm{Interpolate}(\mathrm{low}_i, \mathrm{nominal}, \mathrm{high}_i).
62 /// \f]
63 /// \param name Name of the object.
64 /// \param title Title (for e.g. plotting)
65 /// \param nominal Nominal value of the function.
66 /// \param lowSet Set of down variations.
67 /// \param highSet Set of up variations.
68 /// \param paramSet Parameters that control the interpolation.
69 /// \param takeOwnership If true, the PiecewiseInterpolation object will take ownership of the arguments in the low, high and parameter sets.
70 PiecewiseInterpolation::PiecewiseInterpolation(const char* name, const char* title, const RooAbsReal& nominal,
71  const RooArgList& lowSet,
72  const RooArgList& highSet,
73  const RooArgList& paramSet,
74  Bool_t takeOwnership) :
75  RooAbsReal(name, title),
76  _nominal("!nominal","nominal value", this, (RooAbsReal&)nominal),
77  _lowSet("!lowSet","low-side variation",this),
78  _highSet("!highSet","high-side variation",this),
79  _paramSet("!paramSet","high-side variation",this),
80  _positiveDefinite(false)
81 
82 {
83  // KC: check both sizes
84  if (lowSet.getSize() != highSet.getSize()) {
85  coutE(InputArguments) << "PiecewiseInterpolation::ctor(" << GetName() << ") ERROR: input lists should be of equal length" << endl ;
87  }
88 
89  RooFIter inputIter1 = lowSet.fwdIterator() ;
90  RooAbsArg* comp ;
91  while((comp = inputIter1.next())) {
92  if (!dynamic_cast<RooAbsReal*>(comp)) {
93  coutE(InputArguments) << "PiecewiseInterpolation::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
94  << " in first list is not of type RooAbsReal" << endl ;
96  }
97  _lowSet.add(*comp) ;
98  if (takeOwnership) {
99  _ownedList.addOwned(*comp) ;
100  }
101  }
102 
103 
104  RooFIter inputIter2 = highSet.fwdIterator() ;
105  while((comp = inputIter2.next())) {
106  if (!dynamic_cast<RooAbsReal*>(comp)) {
107  coutE(InputArguments) << "PiecewiseInterpolation::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
108  << " in first list is not of type RooAbsReal" << endl ;
110  }
111  _highSet.add(*comp) ;
112  if (takeOwnership) {
113  _ownedList.addOwned(*comp) ;
114  }
115  }
116 
117 
118  RooFIter inputIter3 = paramSet.fwdIterator() ;
119  while((comp = inputIter3.next())) {
120  if (!dynamic_cast<RooAbsReal*>(comp)) {
121  coutE(InputArguments) << "PiecewiseInterpolation::ctor(" << GetName() << ") ERROR: component " << comp->GetName()
122  << " in first list is not of type RooAbsReal" << endl ;
124  }
125  _paramSet.add(*comp) ;
126  if (takeOwnership) {
127  _ownedList.addOwned(*comp) ;
128  }
129  _interpCode.push_back(0); // default code: linear interpolation
130  }
131 
132 
133  // Choose special integrator by default
134  specialIntegratorConfig(kTRUE)->method1D().setLabel("RooBinIntegrator") ;
136 }
137 
138 
139 
140 ////////////////////////////////////////////////////////////////////////////////
141 /// Copy constructor
142 
144  RooAbsReal(other, name),
145  _nominal("!nominal",this,other._nominal),
146  _lowSet("!lowSet",this,other._lowSet),
147  _highSet("!highSet",this,other._highSet),
148  _paramSet("!paramSet",this,other._paramSet),
149  _positiveDefinite(other._positiveDefinite),
150  _interpCode(other._interpCode)
151 {
152  // Member _ownedList is intentionally not copy-constructed -- ownership is not transferred
154 }
155 
156 
157 
158 ////////////////////////////////////////////////////////////////////////////////
159 /// Destructor
160 
162 {
164 }
165 
166 
167 
168 
169 ////////////////////////////////////////////////////////////////////////////////
170 /// Calculate and return current value of self
171 
173 {
174  ///////////////////
175  Double_t nominal = _nominal;
176  Double_t sum(nominal) ;
177 
178  for (unsigned int i=0; i < _paramSet.size(); ++i) {
179  auto param = static_cast<RooAbsReal*>(_paramSet.at(i));
180  auto low = static_cast<RooAbsReal*>(_lowSet.at(i));
181  auto high = static_cast<RooAbsReal*>(_highSet.at(i));
182  Int_t icode = _interpCode[i] ;
183 
184  switch(icode) {
185  case 0: {
186  // piece-wise linear
187  if(param->getVal()>0)
188  sum += param->getVal()*(high->getVal() - nominal );
189  else
190  sum += param->getVal()*(nominal - low->getVal());
191  break ;
192  }
193  case 1: {
194  // pice-wise log
195  if(param->getVal()>=0)
196  sum *= pow(high->getVal()/nominal, +param->getVal());
197  else
198  sum *= pow(low->getVal()/nominal, -param->getVal());
199  break ;
200  }
201  case 2: {
202  // parabolic with linear
203  double a = 0.5*(high->getVal()+low->getVal())-nominal;
204  double b = 0.5*(high->getVal()-low->getVal());
205  double c = 0;
206  if(param->getVal()>1 ){
207  sum += (2*a+b)*(param->getVal()-1)+high->getVal()-nominal;
208  } else if(param->getVal()<-1 ) {
209  sum += -1*(2*a-b)*(param->getVal()+1)+low->getVal()-nominal;
210  } else {
211  sum += a*pow(param->getVal(),2) + b*param->getVal()+c;
212  }
213  break ;
214  }
215  case 3: {
216  //parabolic version of log-normal
217  double a = 0.5*(high->getVal()+low->getVal())-nominal;
218  double b = 0.5*(high->getVal()-low->getVal());
219  double c = 0;
220  if(param->getVal()>1 ){
221  sum += (2*a+b)*(param->getVal()-1)+high->getVal()-nominal;
222  } else if(param->getVal()<-1 ) {
223  sum += -1*(2*a-b)*(param->getVal()+1)+low->getVal()-nominal;
224  } else {
225  sum += a*pow(param->getVal(),2) + b*param->getVal()+c;
226  }
227  break ;
228  }
229  case 4: {
230 
231  // WVE ****************************************************************
232  // WVE *** THIS CODE IS CRITICAL TO HISTFACTORY FIT CPU PERFORMANCE ***
233  // WVE *** Do not modify unless you know what you are doing... ***
234  // WVE ****************************************************************
235 
236  double x = param->getVal();
237  if (x>1) {
238  sum += x*(high->getVal() - nominal );
239  } else if (x<-1) {
240  sum += x*(nominal - low->getVal());
241  } else {
242  double eps_plus = high->getVal() - nominal;
243  double eps_minus = nominal - low->getVal();
244  double S = 0.5 * (eps_plus + eps_minus);
245  double A = 0.0625 * (eps_plus - eps_minus);
246 
247  //fcns+der+2nd_der are eq at bd
248 
249  double val = nominal + x * (S + x * A * ( 15 + x * x * (-10 + x * x * 3 ) ) );
250 
251 
252  if (val < 0) val = 0;
253  sum += val-nominal;
254  }
255  break ;
256 
257  // WVE ****************************************************************
258  }
259  case 5: {
260 
261  double x0 = 1.0;//boundary;
262  double x = param->getVal();
263 
264  if (x > x0 || x < -x0)
265  {
266  if(x>0)
267  sum += x*(high->getVal() - nominal );
268  else
269  sum += x*(nominal - low->getVal());
270  }
271  else if (nominal != 0)
272  {
273  double eps_plus = high->getVal() - nominal;
274  double eps_minus = nominal - low->getVal();
275  double S = (eps_plus + eps_minus)/2;
276  double A = (eps_plus - eps_minus)/2;
277 
278  //fcns+der are eq at bd
279  double a = S;
280  double b = 3*A/(2*x0);
281  //double c = 0;
282  double d = -A/(2*x0*x0*x0);
283 
284  double val = nominal + a*x + b*pow(x, 2) + 0/*c*pow(x, 3)*/ + d*pow(x, 4);
285  if (val < 0) val = 0;
286 
287  //cout << "Using interp code 5, val = " << val << endl;
288 
289  sum += val-nominal;
290  }
291  break ;
292  }
293  default: {
294  coutE(InputArguments) << "PiecewiseInterpolation::evaluate ERROR: " << param->GetName()
295  << " with unknown interpolation code" << icode << endl ;
296  break ;
297  }
298  }
299  }
300 
301  if(_positiveDefinite && (sum<0)){
302  sum = 0;
303  // cout <<"sum < 0 forcing positive definite"<<endl;
304  // int code = 1;
305  // RooArgSet* myset = new RooArgSet();
306  // cout << "integral = " << analyticalIntegralWN(code, myset) << endl;
307  } else if(sum<0){
308  cxcoutD(Tracing) <<"PiecewiseInterpolation::evaluate - sum < 0, not forcing positive definite"<<endl;
309  }
310  return sum;
311 
312 }
313 
314 
315 ////////////////////////////////////////////////////////////////////////////////
316 /// Interpolate between input distributions for all values of the observable in `evalData`.
317 /// \param[in/out] evalData Struct holding spans pointing to input data. The results of this function will be stored here.
318 /// \param[in] normSet Arguments to normalise over.
320  auto nominal = _nominal->getValues(evalData, normSet);
321  auto sum = evalData.makeBatch(this, nominal.size());
322  std::copy(nominal.begin(), nominal.end(), sum.begin());
323 
324  for (unsigned int i=0; i < _paramSet.size(); ++i) {
325  const double param = static_cast<RooAbsReal*>(_paramSet.at(i))->getVal();
326  auto low = static_cast<RooAbsReal*>(_lowSet.at(i) )->getValues(evalData, normSet);
327  auto high = static_cast<RooAbsReal*>(_highSet.at(i))->getValues(evalData, normSet);
328  const int icode = _interpCode[i];
329 
330  switch(icode) {
331  case 0: {
332  // piece-wise linear
333  for (unsigned int j=0; j < nominal.size(); ++j) {
334  if(param >0)
335  sum[j] += param * (high[j] - nominal[j]);
336  else
337  sum[j] += param * (nominal[j] - low[j] );
338  }
339  break;
340  }
341  case 1: {
342  // pice-wise log
343  for (unsigned int j=0; j < nominal.size(); ++j) {
344  if(param >=0)
345  sum[j] *= pow(high[j]/ nominal[j], +param);
346  else
347  sum[j] *= pow(low[j] / nominal[j], -param);
348  }
349  break;
350  }
351  case 2:
352  // parabolic with linear
353  for (unsigned int j=0; j < nominal.size(); ++j) {
354  const double a = 0.5*(high[j]+low[j])-nominal[j];
355  const double b = 0.5*(high[j]-low[j]);
356  const double c = 0;
357  if (param > 1.) {
358  sum[j] += (2*a+b)*(param -1)+high[j]-nominal[j];
359  } else if (param < -1.) {
360  sum[j] += -1*(2*a-b)*(param +1)+low[j]-nominal[j];
361  } else {
362  sum[j] += a*pow(param ,2) + b*param +c;
363  }
364  }
365  break;
366  case 3: {
367  //parabolic version of log-normal
368  for (unsigned int j=0; j < nominal.size(); ++j) {
369  const double a = 0.5*(high[j]+low[j])-nominal[j];
370  const double b = 0.5*(high[j]-low[j]);
371  const double c = 0;
372  if (param > 1.) {
373  sum[j] += (2*a+b)*(param -1)+high[j]-nominal[j];
374  } else if (param < -1.) {
375  sum[j] += -1*(2*a-b)*(param +1)+low[j]-nominal[j];
376  } else {
377  sum[j] += a*pow(param ,2) + b*param +c;
378  }
379  }
380  break;
381  }
382  case 4:
383  for (unsigned int j=0; j < nominal.size(); ++j) {
384  const double x = param;
385  if (x > 1.) {
386  sum[j] += x * (high[j] - nominal[j]);
387  } else if (x < -1.) {
388  sum[j] += x * (nominal[j] - low[j]);
389  } else {
390  const double eps_plus = high[j] - nominal[j];
391  const double eps_minus = nominal[j] - low[j];
392  const double S = 0.5 * (eps_plus + eps_minus);
393  const double A = 0.0625 * (eps_plus - eps_minus);
394 
395  double val = nominal[j] + x * (S + x * A * ( 15. + x * x * (-10. + x * x * 3. ) ) );
396 
397  if (val < 0.) val = 0.;
398  sum[j] += val - nominal[j];
399  }
400  }
401  break;
402  case 5:
403  for (unsigned int j=0; j < nominal.size(); ++j) {
404  if (param > 1. || param < -1.) {
405  if(param>0)
406  sum[j] += param * (high[j] - nominal[j]);
407  else
408  sum[j] += param * (nominal[j] - low[j] );
409  } else if (nominal[j] != 0) {
410  const double eps_plus = high[j] - nominal[j];
411  const double eps_minus = nominal[j] - low[j];
412  const double S = (eps_plus + eps_minus)/2;
413  const double A = (eps_plus - eps_minus)/2;
414 
415  //fcns+der are eq at bd
416  const double a = S;
417  const double b = 3*A/(2*1.);
418  //double c = 0;
419  const double d = -A/(2*1.*1.*1.);
420 
421  double val = nominal[j] + a * param + b * pow(param, 2) + d * pow(param, 4);
422  if (val < 0.) val = 0.;
423 
424  sum[j] += val - nominal[j];
425  }
426  }
427  break;
428  default:
429  coutE(InputArguments) << "PiecewiseInterpolation::evaluateSpan(): " << _paramSet[i].GetName()
430  << " with unknown interpolation code" << icode << std::endl;
431  throw std::invalid_argument("PiecewiseInterpolation::evaluateSpan() got invalid interpolation code " + std::to_string(icode));
432  break;
433  }
434  }
435 
436  if (_positiveDefinite) {
437  for (double& val : sum) {
438  if (val < 0.)
439  val = 0.;
440  }
441  }
442 
443  return sum;
444 }
445 
446 ////////////////////////////////////////////////////////////////////////////////
447 
449 {
450  if(allVars.getSize()==1){
451  RooAbsReal* temp = const_cast<PiecewiseInterpolation*>(this);
452  temp->specialIntegratorConfig(kTRUE)->method1D().setLabel("RooBinIntegrator") ;
453  int nbins = ((RooRealVar*) allVars.first())->numBins();
454  temp->specialIntegratorConfig(kTRUE)->getConfigSection("RooBinIntegrator").setRealValue("numBins",nbins);
455  return true;
456  }else{
457  cout << "Currently BinIntegrator only knows how to deal with 1-d "<<endl;
458  return false;
459  }
460  return false;
461 }
462 
463 ////////////////////////////////////////////////////////////////////////////////
464 /// Advertise that all integrals can be handled internally.
465 
467  const RooArgSet* normSet, const char* /*rangeName*/) const
468 {
469  /*
470  cout << "---------------------------\nin PiecewiseInterpolation get analytic integral " <<endl;
471  cout << "all vars = "<<endl;
472  allVars.Print("v");
473  cout << "anal vars = "<<endl;
474  analVars.Print("v");
475  cout << "normset vars = "<<endl;
476  if(normSet2)
477  normSet2->Print("v");
478  */
479 
480 
481  // Handle trivial no-integration scenario
482  if (allVars.getSize()==0) return 0 ;
483  if (_forceNumInt) return 0 ;
484 
485 
486  // Force using numeric integration
487  // use special numeric integrator
488  return 0;
489 
490 
491  // KC: check if interCode=0 for all
492  RooFIter paramIterExtra(_paramSet.fwdIterator()) ;
493  int i=0;
494  while( paramIterExtra.next() ) {
495  if(!_interpCode.empty() && _interpCode[i]!=0){
496  // can't factorize integral
497  cout <<"can't factorize integral"<<endl;
498  return 0;
499  }
500  ++i;
501  }
502 
503  // Select subset of allVars that are actual dependents
504  analVars.add(allVars) ;
505  // RooArgSet* normSet = normSet2 ? getObservables(normSet2) : 0 ;
506  // RooArgSet* normSet = getObservables();
507  // RooArgSet* normSet = 0;
508 
509 
510  // Check if this configuration was created before
511  Int_t sterileIdx(-1) ;
512  CacheElem* cache = (CacheElem*) _normIntMgr.getObj(normSet,&analVars,&sterileIdx) ;
513  if (cache) {
514  return _normIntMgr.lastIndex()+1 ;
515  }
516 
517  // Create new cache element
518  cache = new CacheElem ;
519 
520  // Make list of function projection and normalization integrals
521  RooAbsReal *func ;
522  // const RooArgSet* nset = _paramList.nset() ;
523 
524  // do nominal
525  func = (RooAbsReal*)(&_nominal.arg()) ;
526  RooAbsReal* funcInt = func->createIntegral(analVars) ;
527  cache->_funcIntList.addOwned(*funcInt) ;
528 
529  // do variations
530  RooFIter lowIter(_lowSet.fwdIterator()) ;
531  RooFIter highIter(_highSet.fwdIterator()) ;
532  RooFIter paramIter(_paramSet.fwdIterator()) ;
533 
534  // int i=0;
535  i=0;
536  while(paramIter.next() ) {
537  func = (RooAbsReal*)lowIter.next() ;
538  funcInt = func->createIntegral(analVars) ;
539  cache->_lowIntList.addOwned(*funcInt) ;
540 
541  func = (RooAbsReal*)highIter.next() ;
542  funcInt = func->createIntegral(analVars) ;
543  cache->_highIntList.addOwned(*funcInt) ;
544  ++i;
545  }
546 
547  // Store cache element
548  Int_t code = _normIntMgr.setObj(normSet,&analVars,(RooAbsCacheElement*)cache,0) ;
549 
550  return code+1 ;
551 }
552 
553 
554 
555 
556 ////////////////////////////////////////////////////////////////////////////////
557 /// Implement analytical integrations by doing appropriate weighting from component integrals
558 /// functions to integrators of components
559 
560 Double_t PiecewiseInterpolation::analyticalIntegralWN(Int_t code, const RooArgSet* /*normSet2*/,const char* /*rangeName*/) const
561 {
562  /*
563  cout <<"Enter analytic Integral"<<endl;
564  printDirty(true);
565  // _nominal.arg().setDirtyInhibit(kTRUE) ;
566  _nominal.arg().setShapeDirty() ;
567  RooAbsReal* temp ;
568  RooFIter lowIter(_lowSet.fwdIterator()) ;
569  while((temp=(RooAbsReal*)lowIter.next())) {
570  // temp->setDirtyInhibit(kTRUE) ;
571  temp->setShapeDirty() ;
572  }
573  RooFIter highIter(_highSet.fwdIterator()) ;
574  while((temp=(RooAbsReal*)highIter.next())) {
575  // temp->setDirtyInhibit(kTRUE) ;
576  temp->setShapeDirty() ;
577  }
578  */
579 
580  /*
581  RooAbsArg::setDirtyInhibit(kTRUE);
582  printDirty(true);
583  cout <<"done setting dirty inhibit = true"<<endl;
584 
585  // old integral, only works for linear and not positive definite
586  CacheElem* cache = (CacheElem*) _normIntMgr.getObjByIndex(code-1) ;
587 
588 
589  std::unique_ptr<RooArgSet> vars2( getParameters(RooArgSet()) );
590  std::unique_ptr<RooArgSet> iset( _normIntMgr.nameSet2ByIndex(code-1)->select(*vars2) );
591  cout <<"iset = "<<endl;
592  iset->Print("v");
593 
594  double sum = 0;
595  RooArgSet* vars = getVariables();
596  vars->remove(_paramSet);
597  _paramSet.Print("v");
598  vars->Print("v");
599  if(vars->getSize()==1){
600  RooRealVar* obs = (RooRealVar*) vars->first();
601  for(int i=0; i<obs->numBins(); ++i){
602  obs->setVal( obs->getMin() + (.5+i)*(obs->getMax()-obs->getMin())/obs->numBins());
603  sum+=evaluate()*(obs->getMax()-obs->getMin())/obs->numBins();
604  cout << "obs = " << obs->getVal() << " sum = " << sum << endl;
605  }
606  } else{
607  cout <<"only know how to deal with 1 observable right now"<<endl;
608  }
609  */
610 
611  /*
612  _nominal.arg().setDirtyInhibit(kFALSE) ;
613  RooFIter lowIter2(_lowSet.fwdIterator()) ;
614  while((temp=(RooAbsReal*)lowIter2.next())) {
615  temp->setDirtyInhibit(kFALSE) ;
616  }
617  RooFIter highIter2(_highSet.fwdIterator()) ;
618  while((temp=(RooAbsReal*)highIter2.next())) {
619  temp->setDirtyInhibit(kFALSE) ;
620  }
621  */
622 
623  /*
624  RooAbsArg::setDirtyInhibit(kFALSE);
625  printDirty(true);
626  cout <<"done"<<endl;
627  cout << "sum = " <<sum<<endl;
628  //return sum;
629  */
630 
631  // old integral, only works for linear and not positive definite
632  CacheElem* cache = (CacheElem*) _normIntMgr.getObjByIndex(code-1) ;
633  if( cache==NULL ) {
634  std::cout << "Error: Cache Element is NULL" << std::endl;
635  throw std::exception();
636  }
637 
638  // old integral, only works for linear and not positive definite
639  RooFIter funcIntIter = cache->_funcIntList.fwdIterator() ;
640  RooFIter lowIntIter = cache->_lowIntList.fwdIterator() ;
641  RooFIter highIntIter = cache->_highIntList.fwdIterator() ;
642  RooAbsReal *funcInt(0), *low(0), *high(0), *param(0) ;
643  Double_t value(0) ;
644  Double_t nominal(0);
645 
646  // get nominal
647  int i=0;
648  while(( funcInt = (RooAbsReal*)funcIntIter.next())) {
649  value += funcInt->getVal() ;
650  nominal = value;
651  i++;
652  }
653  if(i==0 || i>1) { cout << "problem, wrong number of nominal functions"<<endl; }
654 
655  // now get low/high variations
656  i = 0;
657  RooFIter paramIter(_paramSet.fwdIterator()) ;
658 
659  // KC: old interp code with new iterator
660  while( (param=(RooAbsReal*)paramIter.next()) ) {
661  low = (RooAbsReal*)lowIntIter.next() ;
662  high = (RooAbsReal*)highIntIter.next() ;
663 
664  if(param->getVal()>0) {
665  value += param->getVal()*(high->getVal() - nominal );
666  } else {
667  value += param->getVal()*(nominal - low->getVal());
668  }
669  ++i;
670  }
671 
672  /* // MB : old bit of interpolation code
673  while( (param=(RooAbsReal*)_paramIter->Next()) ) {
674  low = (RooAbsReal*)lowIntIter->Next() ;
675  high = (RooAbsReal*)highIntIter->Next() ;
676 
677  if(param->getVal()>0) {
678  value += param->getVal()*(high->getVal() - nominal );
679  } else {
680  value += param->getVal()*(nominal - low->getVal());
681  }
682  ++i;
683  }
684  */
685 
686  /* KC: the code below is wrong. Can't pull out a constant change to a non-linear shape deformation.
687  while( (param=(RooAbsReal*)paramIter.next()) ) {
688  low = (RooAbsReal*)lowIntIter.next() ;
689  high = (RooAbsReal*)highIntIter.next() ;
690 
691  if(_interpCode.empty() || _interpCode.at(i)==0){
692  // piece-wise linear
693  if(param->getVal()>0)
694  value += param->getVal()*(high->getVal() - nominal );
695  else
696  value += param->getVal()*(nominal - low->getVal());
697  } else if(_interpCode.at(i)==1){
698  // pice-wise log
699  if(param->getVal()>=0)
700  value *= pow(high->getVal()/nominal, +param->getVal());
701  else
702  value *= pow(low->getVal()/nominal, -param->getVal());
703  } else if(_interpCode.at(i)==2){
704  // parabolic with linear
705  double a = 0.5*(high->getVal()+low->getVal())-nominal;
706  double b = 0.5*(high->getVal()-low->getVal());
707  double c = 0;
708  if(param->getVal()>1 ){
709  value += (2*a+b)*(param->getVal()-1)+high->getVal()-nominal;
710  } else if(param->getVal()<-1 ) {
711  value += -1*(2*a-b)*(param->getVal()+1)+low->getVal()-nominal;
712  } else {
713  value += a*pow(param->getVal(),2) + b*param->getVal()+c;
714  }
715  } else if(_interpCode.at(i)==3){
716  //parabolic version of log-normal
717  double a = 0.5*(high->getVal()+low->getVal())-nominal;
718  double b = 0.5*(high->getVal()-low->getVal());
719  double c = 0;
720  if(param->getVal()>1 ){
721  value += (2*a+b)*(param->getVal()-1)+high->getVal()-nominal;
722  } else if(param->getVal()<-1 ) {
723  value += -1*(2*a-b)*(param->getVal()+1)+low->getVal()-nominal;
724  } else {
725  value += a*pow(param->getVal(),2) + b*param->getVal()+c;
726  }
727 
728  } else {
729  coutE(InputArguments) << "PiecewiseInterpolation::analyticalIntegralWN ERROR: " << param->GetName()
730  << " with unknown interpolation code" << endl ;
731  }
732  ++i;
733  }
734  */
735 
736  // cout << "value = " << value <<endl;
737  return value;
738 }
739 
740 
741 ////////////////////////////////////////////////////////////////////////////////
742 
744  int index = _paramSet.index(&param);
745  if(index<0){
746  coutE(InputArguments) << "PiecewiseInterpolation::setInterpCode ERROR: " << param.GetName()
747  << " is not in list" << endl ;
748  } else {
749  coutW(InputArguments) << "PiecewiseInterpolation::setInterpCode : " << param.GetName()
750  << " is now " << code << endl ;
751  _interpCode.at(index) = code;
752  }
753 }
754 
755 
756 ////////////////////////////////////////////////////////////////////////////////
757 
759  for(unsigned int i=0; i<_interpCode.size(); ++i){
760  _interpCode.at(i) = code;
761  }
762 }
763 
764 
765 ////////////////////////////////////////////////////////////////////////////////
766 
768  for(unsigned int i=0; i<_interpCode.size(); ++i){
769  coutI(InputArguments) <<"interp code for " << _paramSet.at(i)->GetName() << " = " << _interpCode.at(i) <<endl;
770  }
771 }
772 
773 
774 ////////////////////////////////////////////////////////////////////////////////
775 /// WVE note: assumes nominal and alternates have identical structure, must add explicit check
776 
777 std::list<Double_t>* PiecewiseInterpolation::binBoundaries(RooAbsRealLValue& obs, Double_t xlo, Double_t xhi) const
778 {
779  return _nominal.arg().binBoundaries(obs,xlo,xhi) ;
780 }
781 
782 
783 ////////////////////////////////////////////////////////////////////////////////
784 /// WVE note: assumes nominal and alternates have identical structure, must add explicit check
785 
787 {
788  return _nominal.arg().isBinnedDistribution(obs) ;
789 }
790 
791 
792 
793 ////////////////////////////////////////////////////////////////////////////////
794 
796 {
797  return _nominal.arg().plotSamplingHint(obs,xlo,xhi) ;
798 }
799 
800 ////////////////////////////////////////////////////////////////////////////////
801 /// Stream an object of class PiecewiseInterpolation.
802 
803 void PiecewiseInterpolation::Streamer(TBuffer &R__b)
804 {
805  if (R__b.IsReading()) {
807  specialIntegratorConfig(kTRUE)->method1D().setLabel("RooBinIntegrator") ;
808  if (_interpCode.empty()) _interpCode.resize(_paramSet.getSize());
809  } else {
811  }
812 }
813 
814 
815 /*
816 ////////////////////////////////////////////////////////////////////////////////
817 /// Customized printing of arguments of a PiecewiseInterpolation to more intuitively reflect the contents of the
818 /// product operator construction
819 
820 void PiecewiseInterpolation::printMetaArgs(ostream& os) const
821 {
822  _lowIter->Reset() ;
823  if (_highIter) {
824  _highIter->Reset() ;
825  }
826 
827  Bool_t first(kTRUE) ;
828 
829  RooAbsArg* arg1, *arg2 ;
830  if (_highSet.getSize()!=0) {
831 
832  while((arg1=(RooAbsArg*)_lowIter->Next())) {
833  if (!first) {
834  os << " + " ;
835  } else {
836  first = kFALSE ;
837  }
838  arg2=(RooAbsArg*)_highIter->Next() ;
839  os << arg1->GetName() << " * " << arg2->GetName() ;
840  }
841 
842  } else {
843 
844  while((arg1=(RooAbsArg*)_lowIter->Next())) {
845  if (!first) {
846  os << " + " ;
847  } else {
848  first = kFALSE ;
849  }
850  os << arg1->GetName() ;
851  }
852 
853  }
854 
855  os << " " ;
856 }
857 
858 */
RooErrorHandler::softAbort
static void softAbort()
Definition: RooErrorHandler.h:30
c
#define c(i)
Definition: RSha256.hxx:101
RooCacheManager::setObj
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=0)
Definition: RooCacheManager.h:50
RooNumIntConfig::getConfigSection
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
Definition: RooNumIntConfig.cxx:222
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:100
RooAbsReal.h
RooBatchCompute::RunContext::makeBatch
RooSpan< double > makeBatch(const RooAbsReal *owner, std::size_t size)
Create a writable batch.
Definition: RunContext.cxx:87
PiecewiseInterpolation::CacheElem::_highIntList
RooArgList _highIntList
Definition: PiecewiseInterpolation.h:83
RooMsgService.h
RooAbsCollection::first
RooAbsArg * first() const
Definition: RooAbsCollection.h:236
RooAbsReal::plotSamplingHint
virtual std::list< Double_t > * plotSamplingHint(RooAbsRealLValue &obs, Double_t xlo, Double_t xhi) const
Interface for returning an optional hint for initial sampling points when constructing a curve projec...
Definition: RooAbsReal.cxx:3816
RooNumIntConfig.h
PiecewiseInterpolation::analyticalIntegralWN
Double_t analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=0) const
Implement analytical integrations by doing appropriate weighting from component integrals functions t...
Definition: PiecewiseInterpolation.cxx:560
RooFit.h
RooTemplateProxy::arg
const T & arg() const
Return reference to object held in proxy.
Definition: RooTemplateProxy.h:290
RooFit::InputArguments
@ InputArguments
Definition: RooGlobalFunc.h:61
RooCacheManager::getObj
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=0, const TNamed *isetRangeName=0)
Definition: RooCacheManager.h:45
RooArgSet.h
PiecewiseInterpolation::_normIntMgr
RooObjCacheManager _normIntMgr
Definition: PiecewiseInterpolation.h:86
RooAbsReal::isBinnedDistribution
virtual Bool_t isBinnedDistribution(const RooArgSet &) const
Tests if the distribution is binned. Unless overridden by derived classes, this always returns false.
Definition: RooAbsReal.h:341
PiecewiseInterpolation::setInterpCode
void setInterpCode(RooAbsReal &param, int code)
Definition: PiecewiseInterpolation.cxx:743
RooAbsReal::createIntegral
RooAbsReal * createIntegral(const RooArgSet &iset, const RooCmdArg &arg1, const RooCmdArg &arg2=RooCmdArg::none(), const RooCmdArg &arg3=RooCmdArg::none(), const RooCmdArg &arg4=RooCmdArg::none(), const RooCmdArg &arg5=RooCmdArg::none(), const RooCmdArg &arg6=RooCmdArg::none(), const RooCmdArg &arg7=RooCmdArg::none(), const RooCmdArg &arg8=RooCmdArg::none()) const
Create an object that represents the integral of the function over one or more observables listed in ...
Definition: RooAbsReal.cxx:548
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
RooNumIntConfig::method1D
RooCategory & method1D()
Definition: RooNumIntConfig.h:34
PiecewiseInterpolation::CacheElem::_funcIntList
RooArgList _funcIntList
Definition: PiecewiseInterpolation.h:81
TBuffer::ReadClassBuffer
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
RooAbsCollection::fwdIterator
RooFIter fwdIterator() const
One-time forward iterator.
Definition: RooAbsCollection.h:193
sum
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345
coutE
#define coutE(a)
Definition: RooMsgService.h:33
RooAbsReal::binBoundaries
virtual std::list< Double_t > * binBoundaries(RooAbsRealLValue &obs, Double_t xlo, Double_t xhi) const
Retrieve bin boundaries if this distribution is binned in obs.
Definition: RooAbsReal.cxx:3805
coutW
#define coutW(a)
Definition: RooMsgService.h:32
RooArgList
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooAbsReal::getVal
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
PiecewiseInterpolation.h
RooAbsCollection::setRealValue
Bool_t setRealValue(const char *name, Double_t newVal=0, Bool_t verbose=kFALSE)
Set value of a RooAbsRealLValye stored in set with given name to newVal No error messages are printed...
Definition: RooAbsCollection.cxx:874
PiecewiseInterpolation::CacheElem::_lowIntList
RooArgList _lowIntList
Definition: PiecewiseInterpolation.h:82
x
Double_t x[n]
Definition: legend1.C:17
RooAbsReal::_forceNumInt
Bool_t _forceNumInt
Definition: RooAbsReal.h:480
coutI
#define coutI(a)
Definition: RooMsgService.h:30
PiecewiseInterpolation::isBinnedDistribution
virtual Bool_t isBinnedDistribution(const RooArgSet &obs) const
WVE note: assumes nominal and alternates have identical structure, must add explicit check.
Definition: PiecewiseInterpolation.cxx:786
RooFitShortHand::S
RooArgSet S(const RooAbsArg &v1)
Definition: RooGlobalFunc.cxx:385
RooAbsReal
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:61
RooArgList::at
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:72
TBuffer
Buffer base class used for serializing objects.
Definition: TBuffer.h:43
ROOT::Math::Cephes::A
static double A[]
Definition: SpecFuncCephes.cxx:170
RooAbsCacheElement
RooAbsCacheElement is the abstract base class for objects to be stored in RooAbsCache cache manager o...
Definition: RooAbsCacheElement.h:26
PiecewiseInterpolation::_ownedList
RooArgList _ownedList
Definition: PiecewiseInterpolation.h:89
RooAbsCollection::GetName
const char * GetName() const
Returns name of object.
Definition: RooAbsCollection.h:286
PiecewiseInterpolation::setBinIntegrator
Bool_t setBinIntegrator(RooArgSet &allVars)
Definition: PiecewiseInterpolation.cxx:448
PiecewiseInterpolation::_lowSet
RooListProxy _lowSet
Definition: PiecewiseInterpolation.h:90
b
#define b(i)
Definition: RSha256.hxx:100
bool
PiecewiseInterpolation::_positiveDefinite
Bool_t _positiveDefinite
Definition: PiecewiseInterpolation.h:94
PiecewiseInterpolation::CacheElem
Definition: PiecewiseInterpolation.h:71
RooTrace.h
PiecewiseInterpolation::plotSamplingHint
virtual std::list< Double_t > * plotSamplingHint(RooAbsRealLValue &obs, Double_t xlo, Double_t xhi) const
Interface for returning an optional hint for initial sampling points when constructing a curve projec...
Definition: PiecewiseInterpolation.cxx:795
PiecewiseInterpolation::_nominal
RooRealProxy _nominal
Definition: PiecewiseInterpolation.h:88
RooCategory::setLabel
virtual Bool_t setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
Definition: RooCategory.cxx:185
RooCacheManager::getObjByIndex
T * getObjByIndex(Int_t index) const
Retrieve payload object by slot index.
Definition: RooCacheManager.h:310
TBuffer.h
TRACE_DESTROY
#define TRACE_DESTROY
Definition: RooTrace.h:24
RooAbsPdf.h
a
auto * a
Definition: textangle.C:12
PiecewiseInterpolation::PiecewiseInterpolation
PiecewiseInterpolation()
Definition: PiecewiseInterpolation.cxx:50
RooFIter
A one-time forward iterator working on RooLinkedList or RooAbsCollection.
Definition: RooLinkedListIter.h:40
TBuffer::WriteClassBuffer
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
RooFIter::next
RooAbsArg * next()
Return next element or nullptr if at end.
Definition: RooLinkedListIter.h:49
TRACE_CREATE
#define TRACE_CREATE
Definition: RooTrace.h:23
PiecewiseInterpolation::evaluate
Double_t evaluate() const
Calculate and return current value of self.
Definition: PiecewiseInterpolation.cxx:172
RooAbsCollection::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Definition: RooAbsCollection.cxx:455
RooAbsCollection::size
Storage_t::size_type size() const
Definition: RooAbsCollection.h:214
RooRealVar.h
RooAbsCollection::addOwned
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Definition: RooAbsCollection.cxx:403
PiecewiseInterpolation::getAnalyticalIntegralWN
Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &analVars, const RooArgSet *normSet, const char *rangeName=0) const
Advertise that all integrals can be handled internally.
Definition: PiecewiseInterpolation.cxx:466
RooFit::Tracing
@ Tracing
Definition: RooGlobalFunc.h:61
TBuffer::IsReading
Bool_t IsReading() const
Definition: TBuffer.h:86
PiecewiseInterpolation::_paramSet
RooListProxy _paramSet
Definition: PiecewiseInterpolation.h:92
PiecewiseInterpolation::evaluateSpan
RooSpan< double > evaluateSpan(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet) const
Interpolate between input distributions for all values of the observable in evalData.
Definition: PiecewiseInterpolation.cxx:319
Double_t
double Double_t
Definition: RtypesCore.h:59
RooAbsReal::specialIntegratorConfig
RooNumIntConfig * specialIntegratorConfig() const
Returns the specialized integrator configuration for this RooAbsReal.
Definition: RooAbsReal.cxx:3523
name
char name[80]
Definition: TGX11.cxx:110
d
#define d(i)
Definition: RSha256.hxx:102
RooErrorHandler.h
RunContext.h
RooCacheManager::lastIndex
Int_t lastIndex() const
Definition: RooCacheManager.h:66
PiecewiseInterpolation::~PiecewiseInterpolation
virtual ~PiecewiseInterpolation()
Destructor.
Definition: PiecewiseInterpolation.cxx:161
PiecewiseInterpolation::setAllInterpCodes
void setAllInterpCodes(int code)
Definition: PiecewiseInterpolation.cxx:758
RooAbsArg
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:72
RooAbsReal::getValues
virtual RooSpan< const double > getValues(RooBatchCompute::RunContext &evalData, const RooArgSet *normSet=nullptr) const
by this change, please consult the release notes for ROOT 6.24 for guidance on how to make this trans...
Definition: RooAbsReal.cxx:312
RooListProxy::add
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Reimplementation of standard RooArgList::add()
Definition: RooListProxy.cxx:104
PiecewiseInterpolation::binBoundaries
virtual std::list< Double_t > * binBoundaries(RooAbsRealLValue &, Double_t, Double_t) const
WVE note: assumes nominal and alternates have identical structure, must add explicit check.
Definition: PiecewiseInterpolation.cxx:777
TNamed::GetName
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
PiecewiseInterpolation::_highSet
RooListProxy _highSet
Definition: PiecewiseInterpolation.h:91
Class
void Class()
Definition: Class.C:29
RooRealVar
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
Riostream.h
pow
double pow(double, double)
RooAbsRealLValue
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
Definition: RooAbsRealLValue.h:31
PiecewiseInterpolation::_interpCode
std::vector< int > _interpCode
Definition: PiecewiseInterpolation.h:96
RooBatchCompute::RunContext
This struct enables passing computation data around between elements of a computation graph.
Definition: RunContext.h:31
PiecewiseInterpolation
The PiecewiseInterpolation is a class that can morph distributions into each other,...
Definition: PiecewiseInterpolation.h:30
RooAbsCollection::getSize
Int_t getSize() const
Definition: RooAbsCollection.h:231
RooSpan
A simple container to hold a batch of data values.
Definition: RooSpan.h:34
RooArgSet
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:33
RooAbsCollection::index
Int_t index(const RooAbsArg *arg) const
Returns index of given arg, or -1 if arg is not in the collection.
Definition: RooAbsCollection.h:247
cxcoutD
#define cxcoutD(a)
Definition: RooMsgService.h:81
int
PiecewiseInterpolation::printAllInterpCodes
void printAllInterpCodes()
Definition: PiecewiseInterpolation.cxx:767