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
30#include "TMath.h"
31
33
34using namespace std;
35
37
38using namespace RooStats;
39using namespace HistFactory;
40
41////////////////////////////////////////////////////////////////////////////////
42/// Default constructor
43
44FlexibleInterpVar::FlexibleInterpVar()
45{
47 _nominal = 0;
51}
52
53
54////////////////////////////////////////////////////////////////////////////////
55
56FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
57 const RooArgList& paramList,
58 Double_t nominal, vector<double> low, vector<double> high) :
59 RooAbsReal(name, title),
60 _paramList("paramList","List of paramficients",this),
61 _nominal(nominal), _low(low), _high(high), _interpBoundary(1.)
62{
65
66
67 TIterator* paramIter = paramList.createIterator() ;
68 RooAbsArg* param ;
69 while((param = (RooAbsArg*)paramIter->Next())) {
70 if (!dynamic_cast<RooAbsReal*>(param)) {
71 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
72 << " is not of type RooAbsReal" << endl ;
73 R__ASSERT(0) ;
74 }
75 _paramList.add(*param) ;
76 _interpCode.push_back(0); // default code
77 }
78 if (int(_low.size() ) != _paramList.getSize() || _low.size() != _high.size()) {
79 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") invalid input low/high vectors " << endl;
80 R__ASSERT(int(_low.size() ) == _paramList.getSize());
81 R__ASSERT(_low.size() == _high.size());
82 }
83
84 delete paramIter ;
86
87}
88
89////////////////////////////////////////////////////////////////////////////////
90
91FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
92 const RooArgList& paramList,
93 double nominal, const RooArgList& low, const RooArgList& high) :
94 RooAbsReal(name, title),
95 _paramList("paramList","List of paramficients",this),
96 _nominal(nominal), _interpBoundary(1.)
97{
98 RooFIter lowIter = low.fwdIterator() ;
99 RooAbsReal* val ;
100 while ((val = (RooAbsReal*) lowIter.next())) {
101 _low.push_back(val->getVal()) ;
102 }
103
104 RooFIter highIter = high.fwdIterator() ;
105 while ((val = (RooAbsReal*) highIter.next())) {
106 _high.push_back(val->getVal()) ;
107 }
108
109
110 _logInit = kFALSE ;
112
113
114 TIterator* paramIter = paramList.createIterator() ;
115 RooAbsArg* param ;
116 while((param = (RooAbsArg*)paramIter->Next())) {
117 if (!dynamic_cast<RooAbsReal*>(param)) {
118 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
119 << " is not of type RooAbsReal" << endl ;
120 R__ASSERT(0) ;
121 }
122 _paramList.add(*param) ;
123 _interpCode.push_back(0); // default code
124 }
125 if (int(_low.size() ) != _paramList.getSize() || _low.size() != _high.size()) {
126 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") invalid input low/high lists " << endl;
127 R__ASSERT(int(_low.size() ) == _paramList.getSize());
128 R__ASSERT(_low.size() == _high.size());
129 }
130
131 delete paramIter ;
133
134}
135
136
137
138
139////////////////////////////////////////////////////////////////////////////////
140
141FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title,
142 const RooArgList& paramList,
143 double nominal, vector<double> low, vector<double> high,
144 vector<int> code) :
145 RooAbsReal(name, title),
146 _paramList("paramList","List of paramficients",this),
147 _nominal(nominal), _low(low), _high(high), _interpCode(code), _interpBoundary(1.)
148{
149 _logInit = kFALSE ;
151
152
153 TIterator* paramIter = paramList.createIterator() ;
154 RooAbsArg* param ;
155 while((param = (RooAbsArg*)paramIter->Next())) {
156 if (!dynamic_cast<RooAbsReal*>(param)) {
157 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName()
158 << " is not of type RooAbsReal" << endl ;
159 // use R__ASSERT which remains also in release mode
160 R__ASSERT(0) ;
161 }
162 _paramList.add(*param) ;
163 }
164 if (int(_low.size() ) != _paramList.getSize() || _low.size() != _high.size() || _low.size() != _interpCode.size()) {
165 coutE(InputArguments) << "FlexibleInterpVar::ctor(" << GetName() << ") invalid input vectors " << endl;
166 R__ASSERT(int(_low.size() ) == _paramList.getSize());
167 R__ASSERT(_low.size() == _high.size());
168 R__ASSERT(_low.size() == _interpCode.size());
169 }
170 delete paramIter ;
172
173}
174
175////////////////////////////////////////////////////////////////////////////////
176/// Constructor of flat polynomial function
177
178FlexibleInterpVar::FlexibleInterpVar(const char* name, const char* title) :
179 RooAbsReal(name, title),
180 _paramList("paramList","List of coefficients",this),
181 _nominal(0), _interpBoundary(1.)
182{
183 _logInit = kFALSE ;
186}
187
188////////////////////////////////////////////////////////////////////////////////
189
191 RooAbsReal(other, name),
192 _paramList("paramList",this,other._paramList),
193 _nominal(other._nominal), _low(other._low), _high(other._high), _interpCode(other._interpCode), _interpBoundary(other._interpBoundary)
194
195{
196 // Copy constructor
197 _logInit = kFALSE ;
200
201}
202
203
204////////////////////////////////////////////////////////////////////////////////
205/// Destructor
206
208{
209 delete _paramIter ;
211}
212
213
214////////////////////////////////////////////////////////////////////////////////
215
217 int index = _paramList.index(&param);
218 if(index<0){
219 coutE(InputArguments) << "FlexibleInterpVar::setInterpCode ERROR: " << param.GetName()
220 << " is not in list" << endl ;
221 } else {
222 coutW(InputArguments) << "FlexibleInterpVar::setInterpCode : " << param.GetName()
223 << " is now " << code << endl ;
224 _interpCode.at(index) = code;
225 }
226 // GHL: Adding suggestion by Swagato:
227 _logInit = kFALSE ;
229}
230
231////////////////////////////////////////////////////////////////////////////////
232
234 for(unsigned int i=0; i<_interpCode.size(); ++i){
235 _interpCode.at(i) = code;
236 }
237 // GHL: Adding suggestion by Swagato:
238 _logInit = kFALSE ;
240
241}
242
243////////////////////////////////////////////////////////////////////////////////
244
246 coutW(InputArguments) << "FlexibleInterpVar::setNominal : nominal is now " << newNominal << endl ;
247 _nominal = newNominal;
248
249 _logInit = kFALSE ;
250
252}
253
254////////////////////////////////////////////////////////////////////////////////
255
257 int index = _paramList.index(&param);
258 if(index<0){
259 coutE(InputArguments) << "FlexibleInterpVar::setLow ERROR: " << param.GetName()
260 << " is not in list" << endl ;
261 } else {
262 coutW(InputArguments) << "FlexibleInterpVar::setLow : " << param.GetName()
263 << " is now " << newLow << endl ;
264 _low.at(index) = newLow;
265 }
266
267 _logInit = kFALSE ;
268
270}
271
272////////////////////////////////////////////////////////////////////////////////
273
275 int index = _paramList.index(&param);
276 if(index<0){
277 coutE(InputArguments) << "FlexibleInterpVar::setHigh ERROR: " << param.GetName()
278 << " is not in list" << endl ;
279 } else {
280 coutW(InputArguments) << "FlexibleInterpVar::setHigh : " << param.GetName()
281 << " is now " << newHigh << endl ;
282 _high.at(index) = newHigh;
283 }
284
285 _logInit = kFALSE ;
287}
288
289////////////////////////////////////////////////////////////////////////////////
290
292 for(unsigned int i=0; i<_interpCode.size(); ++i){
293 coutI(InputArguments) <<"interp code for " << _paramList.at(i)->GetName() << " = " << _interpCode.at(i) <<endl;
294 // GHL: Adding suggestion by Swagato:
295 if( _low.at(i) <= 0.001 ) coutE(InputArguments) << GetName() << ", " << _paramList.at(i)->GetName() << ": low value = " << _low.at(i) << endl;
296 if( _high.at(i) <= 0.001 ) coutE(InputArguments) << GetName() << ", " << _paramList.at(i)->GetName() << ": high value = " << _high.at(i) << endl;
297 }
298
299}
300
301////////////////////////////////////////////////////////////////////////////////
302
303double FlexibleInterpVar::PolyInterpValue(int i, double x) const {
304 // code for polynomial interpolation used when interpCode=4
305
306 double boundary = _interpBoundary;
307
308 double x0 = boundary;
309
310
311 // cache the polynomial coefficient values
312 // which do not depend on x but on the boundaries values
313 if (!_logInit) {
314
316
317 unsigned int n = _low.size();
318 assert(n == _high.size() );
319
320 _polCoeff.resize(n*6) ;
321
322 for (unsigned int j = 0; j < n ; j++) {
323
324 // location of the 6 coefficient for the j-th variable
325 double * coeff = &_polCoeff[j * 6];
326
327 // GHL: Swagato's suggestions
328 double pow_up = std::pow(_high[j]/_nominal, x0);
329 double pow_down = std::pow(_low[j]/_nominal, x0);
330 double logHi = std::log(_high[j]) ;
331 double logLo = std::log(_low[j] );
332 double pow_up_log = _high[j] <= 0.0 ? 0.0 : pow_up * logHi;
333 double pow_down_log = _low[j] <= 0.0 ? 0.0 : -pow_down * logLo;
334 double pow_up_log2 = _high[j] <= 0.0 ? 0.0 : pow_up_log * logHi;
335 double pow_down_log2= _low[j] <= 0.0 ? 0.0 : -pow_down_log* logLo;
336
337 double S0 = (pow_up+pow_down)/2;
338 double A0 = (pow_up-pow_down)/2;
339 double S1 = (pow_up_log+pow_down_log)/2;
340 double A1 = (pow_up_log-pow_down_log)/2;
341 double S2 = (pow_up_log2+pow_down_log2)/2;
342 double A2 = (pow_up_log2-pow_down_log2)/2;
343
344 //fcns+der+2nd_der are eq at bd
345
346 // cache coefficient of the polynomial
347 coeff[0] = 1./(8*x0) *( 15*A0 - 7*x0*S1 + x0*x0*A2);
348 coeff[1] = 1./(8*x0*x0) *(-24 + 24*S0 - 9*x0*A1 + x0*x0*S2);
349 coeff[2] = 1./(4*pow(x0, 3))*( - 5*A0 + 5*x0*S1 - x0*x0*A2);
350 coeff[3] = 1./(4*pow(x0, 4))*( 12 - 12*S0 + 7*x0*A1 - x0*x0*S2);
351 coeff[4] = 1./(8*pow(x0, 5))*( + 3*A0 - 3*x0*S1 + x0*x0*A2);
352 coeff[5] = 1./(8*pow(x0, 6))*( -8 + 8*S0 - 5*x0*A1 + x0*x0*S2);
353
354 }
355
356 }
357
358 // GHL: Swagato's suggestions
359 // if( _low[i] == 0 ) _low[i] = 0.0001;
360 // if( _high[i] == 0 ) _high[i] = 0.0001;
361
362 // get pointer to location of coefficients in the vector
363
364 assert(int(_polCoeff.size()) > i );
365 const double * coefficients = &_polCoeff.front() + 6*i;
366
367 double a = coefficients[0];
368 double b = coefficients[1];
369 double c = coefficients[2];
370 double d = coefficients[3];
371 double e = coefficients[4];
372 double f = coefficients[5];
373
374
375 // evaluate the 6-th degree polynomial using Horner's method
376 double value = 1. + x * (a + x * ( b + x * ( c + x * ( d + x * ( e + x * f ) ) ) ) );
377 return value;
378}
379
380////////////////////////////////////////////////////////////////////////////////
381/// Calculate and return value of polynomial
382
384{
386 _paramIter->Reset() ;
387
388 RooAbsReal* param ;
389 //const RooArgSet* nset = _paramList.nset() ;
390 int i=0;
391
392 // TString name = GetName();
393 // if (name == TString("ZHWW_ll12_vzll_epsilon") )
394 // // std::cout << "evaluate flexible interp var - init flag is " << _logInit << std::endl;
395
396 while((param=(RooAbsReal*)_paramIter->Next())) {
397 // param->Print("v");
398
399
400 Int_t icode = _interpCode[i] ;
401
402 switch(icode) {
403
404 case 0: {
405 // piece-wise linear
406 if(param->getVal()>0)
407 total += param->getVal()*(_high[i] - _nominal );
408 else
409 total += param->getVal()*(_nominal - _low[i]);
410 break ;
411 }
412 case 1: {
413 // pice-wise log
414 if(param->getVal()>=0)
415 total *= pow(_high[i]/_nominal, +param->getVal());
416 else
417 total *= pow(_low[i]/_nominal, -param->getVal());
418 break ;
419 }
420 case 2: {
421 // parabolic with linear
422 double a = 0.5*(_high[i]+_low[i])-_nominal;
423 double b = 0.5*(_high[i]-_low[i]);
424 double c = 0;
425 if(param->getVal()>1 ){
426 total += (2*a+b)*(param->getVal()-1)+_high[i]-_nominal;
427 } else if(param->getVal()<-1 ) {
428 total += -1*(2*a-b)*(param->getVal()+1)+_low[i]-_nominal;
429 } else {
430 total += a*pow(param->getVal(),2) + b*param->getVal()+c;
431 }
432 break ;
433 }
434 case 3: {
435 //parabolic version of log-normal
436 double a = 0.5*(_high[i]+_low[i])-_nominal;
437 double b = 0.5*(_high[i]-_low[i]);
438 double c = 0;
439 if(param->getVal()>1 ){
440 total += (2*a+b)*(param->getVal()-1)+_high[i]-_nominal;
441 } else if(param->getVal()<-1 ) {
442 total += -1*(2*a-b)*(param->getVal()+1)+_low[i]-_nominal;
443 } else {
444 total += a*pow(param->getVal(),2) + b*param->getVal()+c;
445 }
446 break ;
447 }
448
449 case 4: {
450 double boundary = _interpBoundary;
451 double x = param->getVal();
452 //std::cout << icode << " param " << param->GetName() << " " << param->getVal() << " boundary " << boundary << std::endl;
453
454 if(x >= boundary)
455 {
456 total *= std::pow(_high[i]/_nominal, +param->getVal());
457 }
458 else if (x <= -boundary)
459 {
460 total *= std::pow(_low[i]/_nominal, -param->getVal());
461 }
462 else if (x != 0)
463 {
464 total *= PolyInterpValue(i, x);
465 }
466 break ;
467 }
468 default: {
469 coutE(InputArguments) << "FlexibleInterpVar::evaluate ERROR: " << param->GetName()
470 << " with unknown interpolation code" << endl ;
471 }
472 }
473 ++i;
474 }
475
476 if(total<=0) {
478 }
479
480 return total;
481}
482
483void FlexibleInterpVar::printMultiline(ostream& os, Int_t contents,
485{
487 os << indent << "--- FlexibleInterpVar ---" << endl;
489}
490
492{
493 _paramIter->Reset();
494 for (int i=0;i<(int)_low.size();i++) {
496 os << setw(36) << param->GetName()<<": "<<setw(7) << _low[i]<<" "<<setw(7) << _high[i]
497 <<endl;
498 }
499}
500
501
#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:31
#define coutW(a)
Definition: RooMsgService.h:33
#define coutE(a)
Definition: RooMsgService.h:34
#define TRACE_DESTROY
Definition: RooTrace.h:23
#define TRACE_CREATE
Definition: RooTrace.h:22
int Int_t
Definition: RtypesCore.h:41
const Bool_t kFALSE
Definition: RtypesCore.h:88
bool Bool_t
Definition: RtypesCore.h:59
double Double_t
Definition: RtypesCore.h:55
const Bool_t kTRUE
Definition: RtypesCore.h:87
#define ClassImp(name)
Definition: Rtypes.h:365
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:71
void setValueDirty()
Mark the element dirty. This forces a re-evaluation when a value is requested.
Definition: RooAbsArg.h:466
RooFIter fwdIterator() const R__SUGGEST_ALTERNATIVE("begin()
One-time forward iterator.
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.
TIterator * createIterator(Bool_t dir=kIterForward) const R__SUGGEST_ALTERNATIVE("begin()
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:59
virtual void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const
Structure printing.
Definition: RooAbsReal.cxx:477
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:87
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
@ InputArguments
Definition: RooGlobalFunc.h:68
Namespace for the RooStats classes.
Definition: Asimov.h:20
static T Min()
Returns maximum representation for type T.
Definition: TMath.h:911
auto * a
Definition: textangle.C:12