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