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)
84 coutE(
InputArguments) <<
"PiecewiseInterpolation::ctor(" <<
GetName() <<
") ERROR: input lists should be of equal length" << endl ;
88 for (
auto *comp : lowSet) {
91 <<
" in first list is not of type RooAbsReal" << endl ;
101 for (
auto *comp : highSet) {
104 <<
" in first list is not of type RooAbsReal" << endl ;
114 for (
auto *comp : paramSet) {
117 <<
" in first list is not of type RooAbsReal" << endl ;
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)
172 double sum(nominal) ;
183 if(param->getVal()>0)
184 sum += param->getVal()*(high->getVal() - nominal );
186 sum += param->getVal()*(nominal - low->getVal());
191 if(param->getVal()>=0)
192 sum *=
pow(high->getVal()/nominal, +param->getVal());
194 sum *=
pow(low->getVal()/nominal, -param->getVal());
199 double a = 0.5*(high->getVal()+low->getVal())-nominal;
200 double b = 0.5*(high->getVal()-low->getVal());
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;
207 sum +=
a*
pow(param->getVal(),2) +
b*param->getVal()+
c;
213 double a = 0.5*(high->getVal()+low->getVal())-nominal;
214 double b = 0.5*(high->getVal()-low->getVal());
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;
221 sum +=
a*
pow(param->getVal(),2) +
b*param->getVal()+
c;
232 double x = param->getVal();
234 sum +=
x*(high->getVal() - nominal );
236 sum +=
x*(nominal - low->getVal());
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);
245 double val = nominal +
x * (
S +
x *
A * ( 15 +
x *
x * (-10 +
x *
x * 3 ) ) );
248 if (val < 0) val = 0;
258 double x = param->getVal();
260 if (
x > x0 ||
x < -x0)
263 sum +=
x*(high->getVal() - nominal );
265 sum +=
x*(nominal - low->getVal());
267 else if (nominal != 0)
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;
276 double b = 3*
A/(2*x0);
278 double d = -
A/(2*x0*x0*x0);
280 double val = nominal +
a*
x +
b*
pow(
x, 2) + 0 +
d*
pow(
x, 4);
281 if (val < 0) val = 0;
291 <<
" with unknown interpolation code" << icode << endl ;
304 cxcoutD(
Tracing) <<
"PiecewiseInterpolation::evaluate - sum < 0, not forcing positive definite"<<endl;
317 for(
unsigned int j=0; j < nominal.size(); ++j) {
330 for (
unsigned int j=0; j < nominal.size(); ++j) {
332 sum[j] += param * (high[j] - nominal[j]);
334 sum[j] += param * (nominal[j] - low[j] );
340 for (
unsigned int j=0; j < nominal.size(); ++j) {
342 sum[j] *=
pow(high[j]/ nominal[j], +param);
344 sum[j] *=
pow(low[j] / nominal[j], -param);
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]);
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];
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]);
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];
380 for (
unsigned int j=0; j < nominal.size(); ++j) {
381 const double x = param;
383 sum[j] +=
x * (high[j] - nominal[j]);
384 }
else if (
x < -1.) {
385 sum[j] +=
x * (nominal[j] - low[j]);
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);
392 double val = nominal[j] +
x * (
S +
x *
A * ( 15. +
x *
x * (-10. +
x *
x * 3. ) ) );
394 if (val < 0.) val = 0.;
395 sum[j] += val - nominal[j];
400 for (
unsigned int j=0; j < nominal.size(); ++j) {
401 if (param > 1. || param < -1.) {
403 sum[j] += param * (high[j] - nominal[j]);
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;
414 const double b = 3*
A/(2*1.);
416 const double d = -
A/(2*1.*1.*1.);
418 double val = nominal[j] +
a * param +
b *
pow(param, 2) +
d *
pow(param, 4);
419 if (val < 0.) val = 0.;
421 sum[j] += val - nominal[j];
427 <<
" with unknown interpolation code" << icode << std::endl;
428 throw std::invalid_argument(
"PiecewiseInterpolation::evaluateSpan() got invalid interpolation code " + std::to_string(icode));
434 for(
unsigned int j=0; j < nominal.size(); ++j) {
452 cout <<
"Currently BinIntegrator only knows how to deal with 1-d "<<endl;
462 const RooArgSet* normSet,
const char* )
const
477 if (allVars.
empty())
return 0 ;
490 cout <<
"can't factorize integral" << endl;
496 analVars.
add(allVars) ;
499 Int_t sterileIdx(-1) ;
611 if( cache==
nullptr ) {
612 std::cout <<
"Error: Cache Element is nullptr" << std::endl;
613 throw std::exception();
624 for (
auto funcInt : static_range_cast<RooAbsReal*>(cache->
_funcIntList)) {
625 value += funcInt->getVal() ;
629 if(i==0 || i>1) { cout <<
"problem, wrong number of nominal functions"<<endl; }
635 for (
auto const *param : static_range_cast<RooAbsReal *>(
_paramSet)) {
639 if(param->getVal() > 0) {
640 value += param->getVal()*(high->
getVal() - nominal);
642 value += param->getVal()*(nominal - low->
getVal());
722 <<
" is not in list" << endl ;
726 <<
" is now " << code << endl ;
797void PiecewiseInterpolation::printMetaArgs(ostream& os) const
806 RooAbsArg* arg1, *arg2 ;
807 if (_highSet.getSize()!=0) {
809 while((arg1=(RooAbsArg*)_lowIter->Next())) {
815 arg2=(RooAbsArg*)_highIter->Next() ;
816 os << arg1->GetName() << " * " << arg2->GetName() ;
821 while((arg1=(RooAbsArg*)_lowIter->Next())) {
827 os << arg1->GetName() ;
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
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 ¶m, int code, bool silent=false)
~PiecewiseInterpolation() override
Destructor.
void Streamer(TBuffer &) override
Stream an object of class PiecewiseInterpolation.
void setAllInterpCodes(int code)
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 printAllInterpCodes()
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...
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...
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.
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...
bool _forceNumInt
Force numerical integration if flag set.
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.
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.
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
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)
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
RooRealVar represents a variable that can be changed from the outside.
const T & arg() const
Return reference to object held in proxy.
Buffer base class used for serializing objects.
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=nullptr)=0
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
const char * GetName() const override
Returns name of object.
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
RooArgSet S(Args_t &&... args)
static uint64_t sum(uint64_t i)