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