diff --git a/interface/AtlasPdfs.h b/interface/AtlasPdfs.h index 11964a8ba66..4dcf1239d43 100644 --- a/interface/AtlasPdfs.h +++ b/interface/AtlasPdfs.h @@ -387,4 +387,130 @@ class RooStarMomentMorph : public RooAbsPdf { #endif +#ifndef QuickStats_Shapes_RooTwoSidedCBShape +#define QuickStats_Shapes_RooTwoSidedCBShape + +#include "RooAbsPdf.h" +#include "RooRealProxy.h" +#include "RooAbsReal.h" +#include "RooArgSet.h" +#include "RooRealVar.h" + +class RooRealVar; + +class RooTwoSidedCBShape : public RooAbsPdf { +public: + RooTwoSidedCBShape() {} + RooTwoSidedCBShape(const char *name, const char *title, RooAbsReal& _m, + RooAbsReal& _m0, RooAbsReal& _sigma, + RooAbsReal& _alphaLo, RooAbsReal& _nLo, + RooAbsReal& _alphaHi, RooAbsReal& _nHi); + + RooTwoSidedCBShape(const RooTwoSidedCBShape& other, const char* name = 0); + virtual TObject* clone(const char* newname) const { return new RooTwoSidedCBShape(*this,newname); } + + inline virtual ~RooTwoSidedCBShape() { } + + Int_t getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* rangeName=0) const; + Double_t analyticalIntegral(Int_t code, const char* rangeName=0) const; + +protected: + + Double_t gaussianIntegral(Double_t tmin, Double_t tmax) const; + Double_t powerLawIntegral(Double_t tmin, Double_t tmax, Double_t alpha, Double_t n) const; + Double_t evaluate() const; + + RooRealProxy m; + RooRealProxy m0; + RooRealProxy sigma; + RooRealProxy alphaLo; + RooRealProxy nLo; + RooRealProxy alphaHi; + RooRealProxy nHi; + + +private: + + ClassDef(RooTwoSidedCBShape,1) +}; + + +#endif + +// @(#)root/roostats:$Id: cranmer $ +// Author: Kyle Cranmer, Akira Shibata +// Modified by Hongtao Yang for xmlAnaWSBuilder: we need to have both upper and lower uncertainties as functions, not just numbers +/************************************************************************* + * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +#ifndef ROOSTATS_ResponseFunction +#define ROOSTATS_ResponseFunction + +#include +#include +#include +#include +#include + +class RooRealVar; +class RooArgList; + +class ResponseFunction : public RooAbsReal +{ +public: + ResponseFunction(); + ResponseFunction(const char *name, const char *title, + const RooArgList &_paramList, + std::vector nominal, const RooArgList &lowList, const RooArgList &highList, std::vector code); + + ResponseFunction(const ResponseFunction &, const char *); + + void setInterpCode(RooAbsReal ¶m, int code); + void setAllInterpCodes(int code); + void setGlobalBoundary(double boundary) { _interpBoundary = boundary; } + + virtual TObject *clone(const char *newname) const { return new ResponseFunction(*this, newname); } + virtual ~ResponseFunction(); + + virtual void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose = kFALSE, TString indent = "") const; + virtual void printResponseFunctions(std::ostream &os) const; + +private: + double PolyInterpValue(int i, double x) const; + inline double logn(double nominal, double delta, double theta) const { return pow(1 + delta / nominal, theta); }; + inline void cacheCoef(unsigned j) const; + +protected: + RooListProxy _paramList; + std::vector _nominal; + RooListProxy _lowList; + RooListProxy _highList; + std::vector _interpCode; + Double_t _interpBoundary; + + mutable Bool_t _logInit; //! flag used for chaching polynomial coefficients + mutable std::vector _polCoeff; //! cached polynomial coefficients + + Double_t evaluate() const; + +public: + double interpolationBoundary() const { return _interpBoundary; } + const std::vector& interpolationCodes() const { return _interpCode; } + const std::vector& nominal() const { return _nominal; } +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,26,0) + const RooAbsCollection::Storage_t& parameters() const { return _paramList.get(); } + const RooAbsCollection::Storage_t& low() const { return _lowList.get(); } + const RooAbsCollection::Storage_t& high() const { return _highList.get(); } +#endif + + + ClassDef(ResponseFunction, 2) // flexible interpolation +}; + +#endif diff --git a/src/AtlasPdfs.cxx b/src/AtlasPdfs.cxx index 34e758a3ad2..212ef9aa50a 100644 --- a/src/AtlasPdfs.cxx +++ b/src/AtlasPdfs.cxx @@ -1924,3 +1924,546 @@ Bool_t RooStarMomentMorph::setBinIntegrator(RooArgSet& allVars) } return false; } + +#include + +#include "TMath.h" +#include "Math/ProbFuncMathCore.h" + +#include "RooFit.h" +#include "RooAbsReal.h" +#include "RooRealVar.h" +#include "RooMath.h" + + +ClassImp(RooTwoSidedCBShape) + +//_____________________________________________________________________________ +RooTwoSidedCBShape::RooTwoSidedCBShape(const char *name, const char *title, + RooAbsReal& _m, RooAbsReal& _m0, RooAbsReal& _sigma, + RooAbsReal& _alphaLo, RooAbsReal& _nLo, + RooAbsReal& _alphaHi, RooAbsReal& _nHi) : + RooAbsPdf(name, title), + m("m", "Dependent", this, _m), + m0("m0", "M0", this, _m0), + sigma("sigma", "Sigma", this, _sigma), + alphaLo("alphaLo", "Low-side Alpha", this, _alphaLo), + nLo("nLo", "Low-side Order", this, _nLo), + alphaHi("alphaHi", "High-side Alpha", this, _alphaHi), + nHi("nHi", "Hig-side Order", this, _nHi) +{ +} + + +//_____________________________________________________________________________ +RooTwoSidedCBShape::RooTwoSidedCBShape(const RooTwoSidedCBShape& other, const char* name) : + RooAbsPdf(other, name), m("m", this, other.m), m0("m0", this, other.m0), + sigma("sigma", this, other.sigma), + alphaLo("alphaLo", this, other.alphaLo), nLo("nLo", this, other.nLo), + alphaHi("alphaHi", this, other.alphaHi), nHi("nHi", this, other.nHi) +{ +} + + +//_____________________________________________________________________________ +Double_t RooTwoSidedCBShape::evaluate() const { + + Double_t t = (m-m0)/sigma; + + if (t < -alphaLo) { + Double_t a = exp(-0.5*alphaLo*alphaLo); + Double_t b = nLo/alphaLo - alphaLo; + return a/TMath::Power(alphaLo/nLo*(b - t), nLo); + } + else if (t > alphaHi) { + Double_t a = exp(-0.5*alphaHi*alphaHi); + Double_t b = nHi/alphaHi - alphaHi; + return a/TMath::Power(alphaHi/nHi*(b + t), nHi); + } + return exp(-0.5*t*t); +} + + +//_____________________________________________________________________________ +Int_t RooTwoSidedCBShape::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const +{ + if (matchArgs(allVars,analVars,m)) return 1; + return 0; +} + + +//_____________________________________________________________________________ +Double_t RooTwoSidedCBShape::analyticalIntegral(Int_t code, const char* rangeName) const +{ + assert(code==1); + Double_t result = 0; + + Double_t sig = fabs((Double_t)sigma); + Double_t tmin = (m.min(rangeName)-m0)/sig; + Double_t tmax = (m.max(rangeName)-m0)/sig; + + if (tmin < -alphaLo) + result += powerLawIntegral(tmin, TMath::Min(tmax, -alphaLo), alphaLo, nLo); + if (tmin < alphaHi && tmax > -alphaLo) + result += gaussianIntegral(TMath::Max(tmin, -alphaLo), TMath::Min(tmax, alphaHi)); + if (tmax > alphaHi) + result += powerLawIntegral(-tmax, TMath::Min(-tmin, -alphaHi), alphaHi, nHi); + + return sig*result; +} + + +//_____________________________________________________________________________ +Double_t RooTwoSidedCBShape::gaussianIntegral(Double_t tmin, Double_t tmax) const +{ + return sqrt(TMath::TwoPi())*(ROOT::Math::gaussian_cdf(tmax) - ROOT::Math::gaussian_cdf(tmin)); +} + + +//_____________________________________________________________________________ +Double_t RooTwoSidedCBShape::powerLawIntegral(Double_t tmin, Double_t tmax, Double_t alpha, Double_t n) const +{ + // Implement protection for n = 1 from RooCBShape.cxx + bool useLog = false; + if( fabs(n-1.0) < 1.0e-05 ) useLog = true; + + Double_t a = exp(-0.5*alpha*alpha); + Double_t b = n/alpha - alpha; + + if( useLog ) return a * TMath::Power(n/alpha, n) * ( log( b-tmin ) - log( b-tmax ) ); + return a/(1 - n)*( (b - tmin)/(TMath::Power(alpha/n*(b - tmin), n)) - (b - tmax)/(TMath::Power(alpha/n*(b - tmax), n)) ); +} + +// @(#)root/roostats:$Id: cranmer $ +// Author: Kyle Cranmer, Akira Shibata +// Author: Giovanni Petrucciani (UCSD) (log-interpolation) +// Modified by Hongtao Yang for xmlAnaWSBuilder +/************************************************************************* + * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers. * + * All rights reserved. * + * * + * For the licensing terms see $ROOTSYS/LICENSE. * + * For the list of contributors see $ROOTSYS/README/CREDITS. * + *************************************************************************/ + +//_________________________________________________ +/* +BEGIN_HTML +

+

+END_HTML +*/ +// + +#include "RooFit.h" + +#include "Riostream.h" +#include +#include "TMath.h" + +#include "RooAbsReal.h" +#include "RooRealVar.h" +#include "RooArgList.h" +#include "RooMsgService.h" +#include "RooTrace.h" + +#include "TMath.h" + + +ClassImp(ResponseFunction) + +//_____________________________________________________________________________ +ResponseFunction::ResponseFunction() +{ + // Default constructor + _interpBoundary = 1.; + _logInit = kFALSE; + TRACE_CREATE +} + +//_____________________________________________________________________________ +ResponseFunction::ResponseFunction(const char *name, const char *title, + const RooArgList ¶mList, + std::vector nominal, const RooArgList &lowList, const RooArgList &highList, std::vector code) + : RooAbsReal(name, title), + _paramList("paramList", "List of paramficients", this), + _nominal(nominal), + _lowList("lowList", "List of low variation parameters", this), + _highList("highList", "List of high variation parameters", this), + _interpCode(code), + _interpBoundary(1.) +{ + + _logInit = kFALSE; + + std::unique_ptr paramIter(paramList.createIterator()); + RooAbsArg *param; + while ((param = (RooAbsArg *)paramIter->Next())) + { + if (!dynamic_cast(param)) + { + coutE(InputArguments) << "ResponseFunction::ctor(" << GetName() << ") ERROR: paramficient " << param->GetName() + << " is not of type RooAbsReal" << std::endl; + assert(0); + } + _paramList.add(*param); + } + + + std::unique_ptr lowIter(lowList.createIterator()); + RooAbsArg *low; + while ((low = (RooAbsArg *)lowIter->Next())) + { + if (!dynamic_cast(low)) + { + coutE(InputArguments) << "ResponseFunction::ctor(" << GetName() << ") ERROR: paramficient " << low->GetName() + << " is not of type RooAbsReal" << std::endl; + assert(0); + } + _lowList.add(*low); + } + + + std::unique_ptr highIter(highList.createIterator()); + RooAbsArg *high; + while ((high = (RooAbsArg *)highIter->Next())) + { + if (!dynamic_cast(high)) + { + coutE(InputArguments) << "ResponseFunction::ctor(" << GetName() << ") ERROR: paramficient " << high->GetName() + << " is not of type RooAbsReal" << std::endl; + assert(0); + } + _highList.add(*high); + } + + TRACE_CREATE +} + +//_____________________________________________________________________________ +ResponseFunction::ResponseFunction(const ResponseFunction &other, const char *name) + : RooAbsReal(other, name), + _paramList("paramList", this, other._paramList), + _nominal(other._nominal), + _lowList("lowList", this, other._lowList), + _highList("highList", this, other._highList), + _interpCode(other._interpCode), _interpBoundary(other._interpBoundary) + +{ + // Copy constructor + _logInit = kFALSE; + TRACE_CREATE +} + +//_____________________________________________________________________________ +ResponseFunction::~ResponseFunction() +{ + // Destructor + + TRACE_DESTROY +} + +//_____________________________________________________________________________ +void ResponseFunction::setInterpCode(RooAbsReal ¶m, int code) +{ + + int index = _paramList.index(¶m); + if (index < 0) + { + coutE(InputArguments) << "ResponseFunction::setInterpCode ERROR: " << param.GetName() + << " is not in list" << std::endl; + } + else + { + coutW(InputArguments) << "ResponseFunction::setInterpCode : " << param.GetName() + << " is now " << code << std::endl; + _interpCode.at(index) = code; + } + // GHL: Adding suggestion by Swagato: + _logInit = kFALSE; + setValueDirty(); +} + +//_____________________________________________________________________________ +void ResponseFunction::setAllInterpCodes(int code) +{ + + for (unsigned int i = 0; i < _interpCode.size(); ++i) + { + _interpCode.at(i) = code; + } + // GHL: Adding suggestion by Swagato: + _logInit = kFALSE; + setValueDirty(); +} + +//_____________________________________________________________________________ +double ResponseFunction::PolyInterpValue(int i, double x) const +{ + // If this is running for first time + if (!_logInit) + { + // Create polynomial array + const unsigned int n = _lowList.getSize(); + assert(n == (unsigned int)_highList.getSize()); + _polCoeff.resize(n * 6); + // cache the polynomial coefficient values + // which do not dpened on x but on the boundaries values + for (unsigned j = 0; j < n; j++) + cacheCoef(j); + _logInit = kTRUE; + } + else + { + // Check upper and lower variations of i to see whether it is constant + // If not, need to update the coefficient everytime + if (!_highList[i].isConstant() || !_lowList[i].isConstant()) + cacheCoef(i); + } + + // GHL: Swagato's suggestions + // if( _low[i] == 0 ) _low[i] = 0.0001; + // if( _high[i] == 0 ) _high[i] = 0.0001; + + // get pointer to location of coefficients in the vector + const double *coeff = &_polCoeff.front() + 6 * i; + + double a = coeff[0]; + double b = coeff[1]; + double c = coeff[2]; + double d = coeff[3]; + double e = coeff[4]; + double f = coeff[5]; + + // evaluate the 6-th degree polynomial using Horner's method + double value = 1. + x * (a + x * (b + x * (c + x * (d + x * (e + x * f))))); + return value; +} + +//_____________________________________________________________________________ +Double_t ResponseFunction::evaluate() const +{ + // Calculate and return value of polynomial + + Double_t total = 1; + + RooAbsReal *param; + int i = 0; + + for (Int_t i=0; i < _paramList.getSize(); i++) { + const auto& param = static_cast(_paramList[i]); + const auto& low = static_cast(_lowList[i]); + const auto& high = static_cast(_highList[i]); + + Int_t icode = _interpCode[i]; + + double highVar = fabs(high.getVal()); + double nominal = _nominal[i]; + double x = (high.getVal() < 0) ? -param.getVal() : +param.getVal(); + + switch (icode) + { + + case 0: + { + // piece-wise linear: only consider upper uncertainty as it should always be symmetric in XML workspace builder + total *= nominal + x * highVar; + break; + } + case 1: + { + // pice-wise log: only consider upper uncertainty as it should always be symmetric in XML workspace builder + total *= logn(nominal, highVar, x); + break; + } + case 4: + { + double boundary = _interpBoundary; + double lowVar = fabs(low.getVal()); + + if (x >= boundary) + { + total *= logn(nominal, highVar, x); + } + else if (x <= -boundary) + { + total *= logn(nominal, lowVar, x); + } + else if (x != 0) + { + total *= PolyInterpValue(i, x); + } + break; + } + default: + { + coutE(InputArguments) << "ResponseFunction::evaluate ERROR: " << param.GetName() + << " with unknown interpolation code" << std::endl; + } + } + } + + return total; +} + +void ResponseFunction::printMultiline(std::ostream &os, Int_t contents, + Bool_t verbose, TString indent) const +{ + RooAbsReal::printMultiline(os, contents, verbose, indent); + os << indent << "--- ResponseFunction ---" << std::endl; + printResponseFunctions(os); +} + +void ResponseFunction::printResponseFunctions(std::ostream &os) const +{ + for (Int_t i = 0; i < _paramList.getSize(); i++) + { + const auto& param = static_cast(_paramList[i]); + const auto& low = static_cast(_lowList[i]); + const auto& high = static_cast(_highList[i]); + os << std::setw(36) << param.GetName() << ": " + << std::setw(7) << _nominal.at(i) << " " + << std::setw(7) << low.getVal() << " " + << std::setw(7) << high.getVal() + << std::setw(7) << _interpCode.at(i) + << std::endl; + } +} + +void ResponseFunction::cacheCoef(unsigned j) const +{ + // code for polynomial interpolation used when interpCode=4 + double boundary = _interpBoundary; + double x0 = boundary; + + // location of the 6 coefficient for the j-th variable + double *coeff = &_polCoeff[j * 6]; + double highVar = _nominal[j] + fabs((static_cast(_highList[j])).getVal()); + double lowVar = (_nominal[j] * _nominal[j]) / (_nominal[j] + fabs((static_cast(_lowList[j])).getVal())); + + // GHL: Swagato's suggestions + double pow_up = std::pow(highVar / _nominal[j], x0); + double pow_down = std::pow(lowVar / _nominal[j], x0); + double logHi = std::log(highVar); + double logLo = std::log(lowVar); + double pow_up_log = highVar <= 0.0 ? 0.0 : pow_up * logHi; + double pow_down_log = lowVar <= 0.0 ? 0.0 : -pow_down * logLo; + double pow_up_log2 = highVar <= 0.0 ? 0.0 : pow_up_log * logHi; + double pow_down_log2 = lowVar <= 0.0 ? 0.0 : -pow_down_log * logLo; + + double S0 = (pow_up + pow_down) / 2; + double A0 = (pow_up - pow_down) / 2; + double S1 = (pow_up_log + pow_down_log) / 2; + double A1 = (pow_up_log - pow_down_log) / 2; + double S2 = (pow_up_log2 + pow_down_log2) / 2; + double A2 = (pow_up_log2 - pow_down_log2) / 2; + + // fcns+der+2nd_der are eq at bd + + // cache coefficient of the polynomial + coeff[0] = 1. / (8 * x0) * (15 * A0 - 7 * x0 * S1 + x0 * x0 * A2); + coeff[1] = 1. / (8 * x0 * x0) * (-24 + 24 * S0 - 9 * x0 * A1 + x0 * x0 * S2); + coeff[2] = 1. / (4 * pow(x0, 3)) * (-5 * A0 + 5 * x0 * S1 - x0 * x0 * A2); + coeff[3] = 1. / (4 * pow(x0, 4)) * (12 - 12 * S0 + 7 * x0 * A1 - x0 * x0 * S2); + coeff[4] = 1. / (8 * pow(x0, 5)) * (+3 * A0 - 3 * x0 * S1 + x0 * x0 * A2); + coeff[5] = 1. / (8 * pow(x0, 6)) * (-8 + 8 * S0 - 5 * x0 * A1 + x0 * x0 * S2); +} + +#if ROOT_VERSION_CODE >= ROOT_VERSION(6,29,0) + +#include +#include +#include +#include + +using namespace RooFit::Detail; + +namespace { + class ResponseFunctionStreamer : public RooFit::JSONIO::Exporter { +public: + virtual const std::string& key() const override { + static const std::string _key = "resonse"; + return _key; + } + + virtual bool exportObject(RooJSONFactoryWSTool *, const RooAbsArg *func, JSONNode &elem) const override + { + const ResponseFunction *pip = static_cast(func); + elem["type"] << key(); + elem["interpolationCodes"].fill_seq(pip->interpolationCodes()); + elem["vars"].fill_seq(pip->parameters(), [](auto const &f) { return f->GetName(); }); + elem["nom"].fill_seq(pip->nominal()); + elem["high"].fill_seq(pip->high(), [](auto const &f) { return f->GetName(); }); + elem["low"].fill_seq(pip->low(), [](auto const &f) { return f->GetName(); }); + return true; + } +}; +} // namespace + + +namespace { +class ResponseFunctionFactory : public RooFit::JSONIO::Importer { +public: + virtual bool importFunction(RooJSONFactoryWSTool *tool, const JSONNode &p) const override + { + std::string name(RooJSONFactoryWSTool::name(p)); + if (!p.has_child("vars")) { + RooJSONFactoryWSTool::error("no vars of '" + name + "'"); + } + if (!p.has_child("high")) { + RooJSONFactoryWSTool::error("no high variations of '" + name + "'"); + } + if (!p.has_child("low")) { + RooJSONFactoryWSTool::error("no low variations of '" + name + "'"); + } + if (!p.has_child("nom")) { + RooJSONFactoryWSTool::error("no nominal variation of '" + name + "'"); + } + + std::vector nom; + nom << p["nom"]; + + std::vector codes; + codes << p["interpolationCodes"]; + + RooArgList vars; + for (const auto &d : p["vars"].children()) { + std::string objname(RooJSONFactoryWSTool::name(d)); + RooRealVar *obj = tool->request(objname, name); + vars.add(*obj); + } + + RooArgList high; + for (const auto &d : p["high"].children()) { + std::string objname(RooJSONFactoryWSTool::name(d)); + RooAbsReal *obj = tool->request(objname, name); + high.add(*obj); + } + + RooArgList low; + for (const auto &d : p["low"].children()) { + std::string objname(RooJSONFactoryWSTool::name(d)); + RooAbsReal *obj = tool->request(objname, name); + low.add(*obj); + } + + ResponseFunction pip(name.c_str(), name.c_str(), vars, nom, low, high, codes); + + tool->workspace()->import(pip, RooFit::RecycleConflictNodes(true), RooFit::Silence(true)); + return true; + } +}; +} // namespace + + +namespace { + int register_serializations(){ + RooFit::JSONIO::registerImporter("response", std::make_unique()); + RooFit::JSONIO::registerExporter(ResponseFunction::Class(), std::make_unique()); + return 1; + } + int _dummy = register_serializations(); +} + +#endif + diff --git a/src/classes_def.xml b/src/classes_def.xml index 338f4a01a03..c8890186f59 100644 --- a/src/classes_def.xml +++ b/src/classes_def.xml @@ -173,6 +173,8 @@ + +