Logo ROOT  
Reference Guide
FlexibleInterpVar.cxx
Go to the documentation of this file.
1// @(#)root/roostats:$Id: cranmer $
2// Author: Kyle Cranmer, Akira Shibata
3// Author: Giovanni Petrucciani (UCSD) (log-interpolation)
4/*************************************************************************
5 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. *
6 * All rights reserved. *
7 * *
8 * For the licensing terms see $ROOTSYS/LICENSE. *
9 * For the list of contributors see $ROOTSYS/README/CREDITS. *
10 *************************************************************************/
11
12////////////////////////////////////////////////////////////////////////////////
13
14/** \class RooStats::HistFactory::FlexibleInterpVar
15 * \ingroup HistFactory
16 */
17
18#include "RooFit.h"
19
20#include "Riostream.h"
21#include <math.h>
22#include "TMath.h"
23
24#include "RooAbsReal.h"
25#include "RooRealVar.h"
26#include "RooArgList.h"
27#include "RooMsgService.h"
28#include "RooTrace.h"
29
31
32using namespace std;
33
35
36using namespace RooStats;
37using namespace HistFactory;
38
39////////////////////////////////////////////////////////////////////////////////
40/// Default constructor
41
42FlexibleInterpVar::FlexibleInterpVar()
43{
45 _nominal = 0;
49}
50
51
52////////////////////////////////////////////////////////////////////////////////
53
54FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
55 const RooArgList& paramList,
56 Double_t nominal, vector<double> low, vector<double> high) :
57 RooAbsReal(name, title),
58 _paramList("paramList","List of paramficients",this),
59 _nominal(nominal), _low(low), _high(high), _interpBoundary(1.)
60{
63
64
65 TIterator* paramIter = paramList.createIterator() ;
66 RooAbsArg* param ;
67 while((param = (RooAbsArg*)paramIter->Next())) {
68 if (!dynamic_cast<RooAbsReal*>(param)) {
69 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
70 << " is not of type RooAbsReal" << endl ;
71 R__ASSERT(0) ;
72 }
73 _paramList.add(*param) ;
74 _interpCode.push_back(0); // default code
75 }
76 if (int(_low.size() ) != _paramList.getSize() || _low.size() != _high.size()) {
77 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") invalid input low/high vectors " << endl;
78 R__ASSERT(int(_low.size() ) == _paramList.getSize());
79 R__ASSERT(_low.size() == _high.size());
80 }
81
82 delete paramIter ;
84
85}
86
87////////////////////////////////////////////////////////////////////////////////
88
89FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
90 const RooArgList& paramList,
91 double nominal, const RooArgList& low, const RooArgList& high) :
92 RooAbsReal(name, title),
93 _paramList("paramList","List of paramficients",this),
94 _nominal(nominal), _interpBoundary(1.)
95{
96 RooFIter lowIter = low.fwdIterator() ;
97 RooAbsReal* val ;
98 while ((val = (RooAbsReal*) lowIter.next())) {
99 _low.push_back(val->getVal()) ;
100 }
101
102 RooFIter highIter = high.fwdIterator() ;
103 while ((val = (RooAbsReal*) highIter.next())) {
104 _high.push_back(val->getVal()) ;
105 }
106
107
108 _logInit = kFALSE ;
110
111
112 TIterator* paramIter = paramList.createIterator() ;
113 RooAbsArg* param ;
114 while((param = (RooAbsArg*)paramIter->Next())) {
115 if (!dynamic_cast<RooAbsReal*>(param)) {
116 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
117 << " is not of type RooAbsReal" << endl ;
118 R__ASSERT(0) ;
119 }
120 _paramList.add(*param) ;
121 _interpCode.push_back(0); // default code
122 }
123 if (int(_low.size() ) != _paramList.getSize() || _low.size() != _high.size()) {
124 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") invalid input low/high lists " << endl;
125 R__ASSERT(int(_low.size() ) == _paramList.getSize());
126 R__ASSERT(_low.size() == _high.size());
127 }
128
129 delete paramIter ;
131
132}
133
134
135
136
137////////////////////////////////////////////////////////////////////////////////
138
139FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
140 const RooArgList& paramList,
141 double nominal, vector<double> low, vector<double> high,
142 vector<int> code) :
143 RooAbsReal(name, title),
144 _paramList("paramList","List of paramficients",this),
145 _nominal(nominal), _low(low), _high(high), _interpCode(code), _interpBoundary(1.)
146{
147 _logInit = kFALSE ;
149
150
151 TIterator* paramIter = paramList.createIterator() ;
152 RooAbsArg* param ;
153 while((param = (RooAbsArg*)paramIter->Next())) {
154 if (!dynamic_cast<RooAbsReal*>(param)) {
155 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
156 << " is not of type RooAbsReal" << endl ;
157 // use R__ASSERT which remains also in release mode
158 R__ASSERT(0) ;
159 }
160 _paramList.add(*param) ;
161 }
162 if (int(_low.size() ) != _paramList.getSize() || _low.size() != _high.size() || _low.size() != _interpCode.size()) {
163 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") invalid input vectors " << endl;
164 R__ASSERT(int(_low.size() ) == _paramList.getSize());
165 R__ASSERT(_low.size() == _high.size());
166 R__ASSERT(_low.size() == _interpCode.size());
167 }
168 delete paramIter ;
170
171}
172
173////////////////////////////////////////////////////////////////////////////////
174/// Constructor of flat polynomial function
175
176FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title) :
177 RooAbsReal(name, title),
178 _paramList("paramList","List of coefficients",this),
179 _nominal(0), _interpBoundary(1.)
180{
181 _logInit = kFALSE ;
184}
185
186////////////////////////////////////////////////////////////////////////////////
187
189 RooAbsReal(other, name),
190 _paramList("paramList",this,other._paramList),
191 _nominal(other._nominal), _low(other._low), _high(other._high), _interpCode(other._interpCode), _interpBoundary(other._interpBoundary)
192
193{
194 // Copy constructor
195 _logInit = kFALSE ;
198
199}
200
201
202////////////////////////////////////////////////////////////////////////////////
203/// Destructor
204
206{
207 delete _paramIter ;
209}
210
211
212////////////////////////////////////////////////////////////////////////////////
213
215 int index = _paramList.index(&param);
216 if(index<0){
217 coutE(InputArguments) << "FlexibleInterpVar::setInterpCode ERROR: " << param.GetName()
218 << " is not in list" << endl ;
219 } else {
220 coutW(InputArguments) << "FlexibleInterpVar::setInterpCode : " << param.GetName()
221 << " is now " << code << endl ;
222 _interpCode.at(index) = code;
223 }
224 // GHL: Adding suggestion by Swagato:
225 _logInit = kFALSE ;
227}
228
229////////////////////////////////////////////////////////////////////////////////
230
232 for(unsigned int i=0; i<_interpCode.size(); ++i){
233 _interpCode.at(i) = code;
234 }
235 // GHL: Adding suggestion by Swagato:
236 _logInit = kFALSE ;
238
239}
240
241////////////////////////////////////////////////////////////////////////////////
242
244 coutW(InputArguments) << "FlexibleInterpVar::setNominal : nominal is now " << newNominal << endl ;
245 _nominal = newNominal;
246
247 _logInit = kFALSE ;
248
250}
251
252////////////////////////////////////////////////////////////////////////////////
253
255 int index = _paramList.index(&param);
256 if(index<0){
257 coutE(InputArguments) << "FlexibleInterpVar::setLow ERROR: " << param.GetName()
258 << " is not in list" << endl ;
259 } else {
260 coutW(InputArguments) << "FlexibleInterpVar::setLow : " << param.GetName()
261 << " is now " << newLow << endl ;
262 _low.at(index) = newLow;
263 }
264
265 _logInit = kFALSE ;
266
268}
269
270////////////////////////////////////////////////////////////////////////////////
271
273 int index = _paramList.index(&param);
274 if(index<0){
275 coutE(InputArguments) << "FlexibleInterpVar::setHigh ERROR: " << param.GetName()
276 << " is not in list" << endl ;
277 } else {
278 coutW(InputArguments) << "FlexibleInterpVar::setHigh : " << param.GetName()
279 << " is now " << newHigh << endl ;
280 _high.at(index) = newHigh;
281 }
282
283 _logInit = kFALSE ;
285}
286
287////////////////////////////////////////////////////////////////////////////////
288
290 for(unsigned int i=0; i<_interpCode.size(); ++i){
291 coutI(InputArguments) <<"interp code for " << _paramList.at(i)->GetName() << " = " << _interpCode.at(i) <<endl;
292 // GHL: Adding suggestion by Swagato:
293 if( _low.at(i) <= 0.001 ) coutE(InputArguments) << GetName() << ", " << _paramList.at(i)->GetName() << ": low value = " << _low.at(i) << endl;
294 if( _high.at(i) <= 0.001 ) coutE(InputArguments) << GetName() << ", " << _paramList.at(i)->GetName() << ": high value = " << _high.at(i) << endl;
295 }
296
297}
298
299////////////////////////////////////////////////////////////////////////////////
300
301double FlexibleInterpVar::PolyInterpValue(int i, double x) const {
302 // code for polynomial interpolation used when interpCode=4
303
304 double boundary = _interpBoundary;
305
306 double x0 = boundary;
307
308
309 // cache the polynomial coefficient values
310 // which do not depend on x but on the boundaries values
311 if (!_logInit) {
312
314
315 unsigned int n = _low.size();
316 assert(n == _high.size() );
317
318 _polCoeff.resize(n*6) ;
319
320 for (unsigned int j = 0; j < n ; j++) {
321
322 // location of the 6 coefficient for the j-th variable
323 double * coeff = &_polCoeff[j * 6];
324
325 // GHL: Swagato's suggestions
326 double pow_up = std::pow(_high[j]/_nominal, x0);
327 double pow_down = std::pow(_low[j]/_nominal, x0);
328 double logHi = std::log(_high[j]) ;
329 double logLo = std::log(_low[j] );
330 double pow_up_log = _high[j] <= 0.0 ? 0.0 : pow_up * logHi;
331 double pow_down_log = _low[j] <= 0.0 ? 0.0 : -pow_down * logLo;
332 double pow_up_log2 = _high[j] <= 0.0 ? 0.0 : pow_up_log * logHi;
333 double pow_down_log2= _low[j] <= 0.0 ? 0.0 : -pow_down_log* logLo;
334
335 double S0 = (pow_up+pow_down)/2;
336 double A0 = (pow_up-pow_down)/2;
337 double S1 = (pow_up_log+pow_down_log)/2;
338 double A1 = (pow_up_log-pow_down_log)/2;
339 double S2 = (pow_up_log2+pow_down_log2)/2;
340 double A2 = (pow_up_log2-pow_down_log2)/2;
341
342 //fcns+der+2nd_der are eq at bd
343
344 // cache coefficient of the polynomial
345 coeff[0] = 1./(8*x0) *( 15*A0 - 7*x0*S1 + x0*x0*A2);
346 coeff[1] = 1./(8*x0*x0) *(-24 + 24*S0 - 9*x0*A1 + x0*x0*S2);
347 coeff[2] = 1./(4*pow(x0, 3))*( - 5*A0 + 5*x0*S1 - x0*x0*A2);
348 coeff[3] = 1./(4*pow(x0, 4))*( 12 - 12*S0 + 7*x0*A1 - x0*x0*S2);
349 coeff[4] = 1./(8*pow(x0, 5))*( + 3*A0 - 3*x0*S1 + x0*x0*A2);
350 coeff[5] = 1./(8*pow(x0, 6))*( -8 + 8*S0 - 5*x0*A1 + x0*x0*S2);
351
352 }
353
354 }
355
356 // GHL: Swagato's suggestions
357 // if( _low[i] == 0 ) _low[i] = 0.0001;
358 // if( _high[i] == 0 ) _high[i] = 0.0001;
359
360 // get pointer to location of coefficients in the vector
361
362 assert(int(_polCoeff.size()) > i );
363 const double * coefficients = &_polCoeff.front() + 6*i;
364
365 double a = coefficients[0];
366 double b = coefficients[1];
367 double c = coefficients[2];
368 double d = coefficients[3];
369 double e = coefficients[4];
370 double f = coefficients[5];
371
372
373 // evaluate the 6-th degree polynomial using Horner's method
374 double value = 1. + x * (a + x * ( b + x * ( c + x * ( d + x * ( e + x * f ) ) ) ) );
375 return value;
376}
377
378////////////////////////////////////////////////////////////////////////////////
379/// Calculate and return value of polynomial
380
382{
384 _paramIter->Reset() ;
385
386 RooAbsReal* param ;
387 //const RooArgSet* nset = _paramList.nset() ;
388 int i=0;
389
390 // TString name = GetName();
391 // if (name == TString("ZHWW_ll12_vzll_epsilon") )
392 // // std::cout << "evaluate flexible interp var - init flag is " << _logInit << std::endl;
393
394 while((param=(RooAbsReal*)_paramIter->Next())) {
395 // param->Print("v");
396
397
398 Int_t icode = _interpCode[i] ;
399
400 switch(icode) {
401
402 case 0: {
403 // piece-wise linear
404 if(param->getVal()>0)
405 total += param->getVal()*(_high[i] - _nominal );
406 else
407 total += param->getVal()*(_nominal - _low[i]);
408 break ;
409 }
410 case 1: {
411 // pice-wise log
412 if(param->getVal()>=0)
413 total *= pow(_high[i]/_nominal, +param->getVal());
414 else
415 total *= pow(_low[i]/_nominal, -param->getVal());
416 break ;
417 }
418 case 2: {
419 // parabolic with linear
420 double a = 0.5*(_high[i]+_low[i])-_nominal;
421 double b = 0.5*(_high[i]-_low[i]);
422 double c = 0;
423 if(param->getVal()>1 ){
424 total += (2*a+b)*(param->getVal()-1)+_high[i]-_nominal;
425 } else if(param->getVal()<-1 ) {
426 total += -1*(2*a-b)*(param->getVal()+1)+_low[i]-_nominal;
427 } else {
428 total += a*pow(param->getVal(),2) + b*param->getVal()+c;
429 }
430 break ;
431 }
432 case 3: {
433 //parabolic version of log-normal
434 double a = 0.5*(_high[i]+_low[i])-_nominal;
435 double b = 0.5*(_high[i]-_low[i]);
436 double c = 0;
437 if(param->getVal()>1 ){
438 total += (2*a+b)*(param->getVal()-1)+_high[i]-_nominal;
439 } else if(param->getVal()<-1 ) {
440 total += -1*(2*a-b)*(param->getVal()+1)+_low[i]-_nominal;
441 } else {
442 total += a*pow(param->getVal(),2) + b*param->getVal()+c;
443 }
444 break ;
445 }
446
447 case 4: {
448 double boundary = _interpBoundary;
449 double x = param->getVal();
450 //std::cout << icode << " param " << param->GetName() << " " << param->getVal() << " boundary " << boundary << std::endl;
451
452 if(x >= boundary)
453 {
454 total *= std::pow(_high[i]/_nominal, +param->getVal());
455 }
456 else if (x <= -boundary)
457 {
458 total *= std::pow(_low[i]/_nominal, -param->getVal());
459 }
460 else if (x != 0)
461 {
462 total *= PolyInterpValue(i, x);
463 }
464 break ;
465 }
466 default: {
467 coutE(InputArguments) << "FlexibleInterpVar::evaluate ERROR: " << param->GetName()
468 << " with unknown interpolation code" << endl ;
469 }
470 }
471 ++i;
472 }
473
474 if(total<=0) {
476 }
477
478 return total;
479}
480
481void FlexibleInterpVar::printMultiline(ostream& os, Int_t contents,
483{
485 os << indent << "--- FlexibleInterpVar ---" << endl;
487}
488
490{
491 _paramIter->Reset();
492 for (int i=0;i<(int)_low.size();i++) {
494 os << setw(36) << param->GetName()<<": "<<setw(7) << _low[i]<<" "<<setw(7) << _high[i]
495 <<endl;
496 }
497}
498
499
#define d(i)
Definition: RSha256.hxx:102
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define S0(x)
Definition: RSha256.hxx:88
#define S1(x)
Definition: RSha256.hxx:89
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
#define coutI(a)
Definition: RooMsgService.h:30
#define coutW(a)
Definition: RooMsgService.h:32
#define coutE(a)
Definition: RooMsgService.h:33
#define TRACE_DESTROY
Definition: RooTrace.h:23
#define TRACE_CREATE
Definition: RooTrace.h:22
const Bool_t kFALSE
Definition: RtypesCore.h:90
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
static void indent(ostringstream &buf, int indent_level)
#define R__ASSERT(e)
Definition: TError.h:96
static unsigned int total
char name[80]
Definition: TGX11.cxx:109
double pow(double, double)
double log(double)
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:73
void setValueDirty()
Mark the element dirty. This forces a re-evaluation when a value is requested.
Definition: RooAbsArg.h:487
Int_t getSize() const
Int_t index(const RooAbsArg *arg) const
Returns index of given arg, or -1 if arg is not in the collection.
RooFIter fwdIterator() const
One-time forward iterator.
TIterator * createIterator(Bool_t dir=kIterForward) const
TIterator-style iteration over contained elements.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:60
virtual void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const
Structure printing.
Definition: RooAbsReal.cxx:475
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:90
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:21
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:74
A one-time forward iterator working on RooLinkedList or RooAbsCollection.
RooAbsArg * next()
Return next element or nullptr if at end.
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Reimplementation of standard RooArgList::add()
std::vector< double > _polCoeff
flag used for chaching polynomial coefficients
Double_t evaluate() const
cached polynomial coefficients
void setInterpCode(RooAbsReal &param, int code)
virtual void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const
Structure printing.
void setLow(RooAbsReal &param, Double_t newLow)
double PolyInterpValue(int i, double x) const
virtual void printFlexibleInterpVars(std::ostream &os) const
void setHigh(RooAbsReal &param, Double_t newHigh)
Iterator abstract base class.
Definition: TIterator.h:30
virtual void Reset()=0
virtual TObject * Next()=0
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
Basic string class.
Definition: TString.h:131
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
@ HistFactory
Definition: RooGlobalFunc.h:69
@ InputArguments
Definition: RooGlobalFunc.h:68
Namespace for the RooStats classes.
Definition: Asimov.h:19
static T Min()
Returns maximum representation for type T.
Definition: TMath.h:911
auto * a
Definition: textangle.C:12