diff --git a/roofit/roofitcore/inc/RooTruthModel.h b/roofit/roofitcore/inc/RooTruthModel.h index ccc16f34626e1..05a0b33334b9d 100644 --- a/roofit/roofitcore/inc/RooTruthModel.h +++ b/roofit/roofitcore/inc/RooTruthModel.h @@ -39,7 +39,6 @@ class RooTruthModel : public RooResolutionModel { RooTruthModel(const char *name, const char *title, RooAbsRealLValue& x) ; RooTruthModel(const RooTruthModel& other, const char* name=nullptr); TObject* clone(const char* newname) const override { return new RooTruthModel(*this,newname) ; } - ~RooTruthModel() override; Int_t basisCode(const char* name) const override ; diff --git a/roofit/roofitcore/src/RooAbsAnaConvPdf.cxx b/roofit/roofitcore/src/RooAbsAnaConvPdf.cxx index d269383767a18..8972cbe187169 100644 --- a/roofit/roofitcore/src/RooAbsAnaConvPdf.cxx +++ b/roofit/roofitcore/src/RooAbsAnaConvPdf.cxx @@ -467,66 +467,57 @@ Int_t RooAbsAnaConvPdf::getAnalyticalIntegralWN(RooArgSet& allVars, /// Set \f$ x \f$ must be contained in \f$ v \f$ and set \f$ y \f$ must be contained in \f$ w \f$. /// -double RooAbsAnaConvPdf::analyticalIntegralWN(Int_t code, const RooArgSet* normSet, const char* rangeName) const +double RooAbsAnaConvPdf::analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName) const { - // WVE needs adaptation to handle new rangeName feature - - // Handle trivial passthrough scenario - if (code==0) return getVal(normSet) ; - - // Unpack master code - RooArgSet *intCoefSet, *intConvSet, *normCoefSet, *normConvSet ; - _codeReg.retrieve(code-1,intCoefSet,intConvSet,normCoefSet,normConvSet) ; - - Int_t index(0) ; - double answer(0) ; - - if (normCoefSet==nullptr&&normConvSet==nullptr) { - - // Integral over unnormalized function - double integral(0) ; - const TNamed *_rangeName = RooNameReg::ptr(rangeName); - for (auto convArg : _convSet) { - auto conv = static_cast(convArg); - double coef = getCoefNorm(index++,intCoefSet,_rangeName) ; - //cout << "coefInt[" << index << "] = " << coef << " " ; intCoefSet->Print("1") ; - if (coef!=0) { - integral += coef* conv->getNormObj(nullptr,intConvSet,_rangeName)->getVal(); - cxcoutD(Eval) << "RooAbsAnaConv::aiWN(" << GetName() << ") [" << index-1 << "] integral += " << conv->getNorm(intConvSet) << endl ; + // WVE needs adaptation to handle new rangeName feature + + // Handle trivial passthrough scenario + if (code == 0) + return getVal(normSet); + + // Unpack master code + RooArgSet *intCoefSet, *intConvSet, *normCoefSet, *normConvSet; + _codeReg.retrieve(code - 1, intCoefSet, intConvSet, normCoefSet, normConvSet); + + Int_t index(0); + + if (normCoefSet == nullptr && normConvSet == nullptr) { + // Integral over unnormalized function + double integral(0); + const TNamed *rangeNamePtr = RooNameReg::ptr(rangeName); + for (auto *conv : static_range_cast(_convSet)) { + double coef = getCoefNorm(index++, intCoefSet, rangeNamePtr); + if (coef != 0) { + const double term = coef * conv->getNormObj(nullptr, intConvSet, rangeNamePtr)->getVal(); + integral += term; + cxcoutD(Eval) << "RooAbsAnaConv::aiWN(" << GetName() << ") [" << index - 1 << "] integral += " << term + << std::endl; + } } + return integral; + } - } - answer = integral ; - - } else { - - // Integral over normalized function - double integral(0) ; - double norm(0) ; - const TNamed *_rangeName = RooNameReg::ptr(rangeName); - for (auto convArg : _convSet) { - auto conv = static_cast(convArg); - - double coefInt = getCoefNorm(index,intCoefSet,_rangeName) ; - //cout << "coefInt[" << index << "] = " << coefInt << "*" << term << " " << (intCoefSet?*intCoefSet:RooArgSet()) << endl ; - if (coefInt!=0) { - double term = conv->getNormObj(nullptr,intConvSet,_rangeName)->getVal(); - integral += coefInt*term ; - } + // Integral over normalized function + double integral(0); + double norm(0); + const TNamed *rangeNamePtr = RooNameReg::ptr(rangeName); + for (auto *conv : static_range_cast(_convSet)) { - double coefNorm = getCoefNorm(index,normCoefSet) ; - //cout << "coefNorm[" << index << "] = " << coefNorm << "*" << term << " " << (normCoefSet?*normCoefSet:RooArgSet()) << endl ; - if (coefNorm!=0) { - double term = conv->getNormObj(nullptr,normConvSet)->getVal(); - norm += coefNorm*term ; + double coefInt = getCoefNorm(index, intCoefSet, rangeNamePtr); + if (coefInt != 0) { + double term = conv->getNormObj(nullptr, intConvSet, rangeNamePtr)->getVal(); + integral += coefInt * term; } - index++ ; - } - answer = integral/norm ; - } + double coefNorm = getCoefNorm(index, normCoefSet); + if (coefNorm != 0) { + double term = conv->getNormObj(nullptr, normConvSet)->getVal(); + norm += coefNorm * term; + } - return answer ; + index++; + } + return integral / norm; } diff --git a/roofit/roofitcore/src/RooTruthModel.cxx b/roofit/roofitcore/src/RooTruthModel.cxx index 78abcd4a51949..851ed46aa5031 100644 --- a/roofit/roofitcore/src/RooTruthModel.cxx +++ b/roofit/roofitcore/src/RooTruthModel.cxx @@ -25,20 +25,22 @@ as a RooFormulaVar. The 6 basis functions used in B mixing and decay and 2 basi functions used in D mixing have been hand coded for increased execution speed. **/ -#include "Riostream.h" -#include "RooBatchCompute.h" -#include "RooTruthModel.h" -#include "RooGenContext.h" -#include "RooAbsAnaConvPdf.h" +#include -#include "TError.h" +#include +#include +#include + +#include + +#include #include +#include + using namespace std ; ClassImp(RooTruthModel); -; - //////////////////////////////////////////////////////////////////////////////// @@ -61,15 +63,6 @@ RooTruthModel::RooTruthModel(const RooTruthModel& other, const char* name) : -//////////////////////////////////////////////////////////////////////////////// -/// Destructor - -RooTruthModel::~RooTruthModel() -{ -} - - - //////////////////////////////////////////////////////////////////////////////// /// Return basis code for given basis definition string. Return special /// codes for 'known' bases for which compiled definition exists. Return @@ -327,127 +320,159 @@ Int_t RooTruthModel::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVa } +namespace { + +// From asking WolframAlpha: integrate exp(-x/tau) over x. +inline double indefiniteIntegralExpBasisPlus(double x, double tau, double /*dm*/) +{ + // Restrict to positive x + x = std::max(x, 0.0); + return -tau * std::exp(-x / tau); +} + +// From asking WolframAlpha: integrate exp(-x/tau)* x / tau over x. +inline double indefiniteIntegralLinBasisPlus(double x, double tau, double /*dm*/) +{ + // Restrict to positive x + x = std::max(x, 0.0); + return -(tau + x) * std::exp(-x / tau); +} + +// From asking WolframAlpha: integrate exp(-x/tau) * (x / tau)^2 over x. +inline double indefiniteIntegralQuadBasisPlus(double x, double tau, double /*dm*/) +{ + // Restrict to positive x + x = std::max(x, 0.0); + return -(std::exp(-x / tau) * (2 * tau * tau + x * x + 2 * tau * x)) / tau; +} + +// A common factor that appears in the integrals of the trigonometric +// function bases (sin and cos). +inline double commonFactorPlus(double x, double tau, double dm) +{ + const double num = tau * std::exp(-x / tau); + const double den = dm * dm * tau * tau + 1.0; + return num / den; +} + +// A common factor that appears in the integrals of the hyperbolic +// trigonometric function bases (sinh and cosh). +inline double commonFactorHyperbolicPlus(double x, double tau, double dm) +{ + const double num = 2 * tau * std::exp(-x / tau); + const double den = dm * dm * tau * tau - 4.0; + return num / den; +} + +// From asking WolframAlpha: integrate exp(-x/tau)*sin(x*m) over x. +inline double indefiniteIntegralSinBasisPlus(double x, double tau, double dm) +{ + // Restrict to positive x + x = std::max(x, 0.0); + const double fac = commonFactorPlus(x, tau, dm); + // Only multiply with the sine term if the coefficient is non zero, + // i.e. if x was not infinity. Otherwise, we are evaluating the + // sine of infinity, which is NAN! + return fac != 0.0 ? fac * (-tau * dm * std::cos(dm * x) - std::sin(dm * x)) : 0.0; +} + +// From asking WolframAlpha: integrate exp(-x/tau)*cos(x*m) over x. +inline double indefiniteIntegralCosBasisPlus(double x, double tau, double dm) +{ + // Restrict to positive x + x = std::max(x, 0.0); + const double fac = commonFactorPlus(x, tau, dm); + return fac != 0.0 ? fac * (tau * dm * std::sin(dm * x) - std::cos(dm * x)) : 0.0; +} + +// From asking WolframAlpha: integrate exp(-x/tau)*sinh(x*m/2) over x. +inline double indefiniteIntegralSinhBasisPlus(double x, double tau, double dm) +{ + // Restrict to positive x + x = std::max(x, 0.0); + const double fac = commonFactorHyperbolicPlus(x, tau, dm); + const double arg = 0.5 * dm * x; + return fac != 0.0 ? fac * (tau * dm * std::cosh(arg) - 2. * std::sinh(arg)) : 0.0; +} + +// From asking WolframAlpha: integrate exp(-x/tau)*cosh(x*m/2) over x. +inline double indefiniteIntegralCoshBasisPlus(double x, double tau, double dm) +{ + // Restrict to positive x + x = std::max(x, 0.0); + const double fac = commonFactorHyperbolicPlus(x, tau, dm); + const double arg = 0.5 * dm * x; + return fac != 0.0 ? fac * (tau * dm * std::sinh(arg) + 2. * std::cosh(arg)) : 0.0; +} + +// Integrate one of the basis functions. Takes a function that represents the +// indefinite integral, some parameters, and a flag that indicats whether the +// basis function is symmetric or antisymmetric. This information is used to +// evaluate the integrals for the "Minus" and "Sum" cases. +template +double definiteIntegral(Function indefiniteIntegral, double xmin, double xmax, double tau, double dm, + RooTruthModel::BasisSign basisSign, bool isSymmetric) +{ + // Note: isSymmetric == false implies antisymmetric + if (tau == 0.0) + return isSymmetric ? 1.0 : 0.0; + double result = 0.0; + if (basisSign != RooTruthModel::Minus) { + result += indefiniteIntegral(xmax, tau, dm) - indefiniteIntegral(xmin, tau, dm); + } + if (basisSign != RooTruthModel::Plus) { + const double resultMinus = indefiniteIntegral(-xmax, tau, dm) - indefiniteIntegral(-xmin, tau, dm); + result += isSymmetric ? -resultMinus : resultMinus; + } + return result; +} + +} // namespace //////////////////////////////////////////////////////////////////////////////// /// Implement analytical integrals when used as p.d.f and for compiled /// basis functions. -double RooTruthModel::analyticalIntegral(Int_t code, const char* rangeName) const +double RooTruthModel::analyticalIntegral(Int_t code, const char *rangeName) const { + // Code must be 1 + R__ASSERT(code == 1); - // Code must be 1 - R__ASSERT(code==1) ; + // Unconvoluted PDF + if (_basisCode == noBasis) + return 1; - // Unconvoluted PDF - if (_basisCode==noBasis) return 1 ; + // Precompiled basis functions + BasisType basisType = (BasisType)((_basisCode == 0) ? 0 : (_basisCode / 10) + 1); + BasisSign basisSign = (BasisSign)(_basisCode - 10 * (basisType - 1) - 2); - // Precompiled basis functions - BasisType basisType = (BasisType)( (_basisCode == 0) ? 0 : (_basisCode/10) + 1 ); - BasisSign basisSign = (BasisSign)( _basisCode - 10*(basisType-1) - 2 ) ; - //cout << " calling RooTruthModel::analyticalIntegral with basisType " << basisType << endl; + const bool needsDm = + basisType == sinBasis || basisType == cosBasis || basisType == sinhBasis || basisType == coshBasis; - double tau = ((RooAbsReal*)basis().getParameter(1))->getVal() ; + const double tau = ((RooAbsReal *)basis().getParameter(1))->getVal(); + const double dm = + needsDm ? ((RooAbsReal *)basis().getParameter(2))->getVal() : std::numeric_limits::quiet_NaN(); - const double xmin = x.min(rangeName); - const double xmax = x.max(rangeName); - - switch (basisType) { - case expBasis: - { - // WVE fixed for ranges - double result(0) ; - if (tau==0) return 1 ; - if ((basisSign != Minus) && (xmax>0)) { - result += tau*(-exp(-xmax/tau) - -exp(-max(0.,xmin)/tau) ) ; // plus and both - } - if ((basisSign != Plus) && (xmin<0)) { - result -= tau*(-exp(-max(0.,xmin)/tau)) - -tau*exp(-xmax/tau) ; // minus and both - } + const double xmin = x.min(rangeName); + const double xmax = x.max(rangeName); - return result ; - } - case sinBasis: - { - double result(0) ; - if (tau==0) return 0 ; - double dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ; - if (basisSign != Minus) { - double term = exp(-xmax/tau); - // We only multiply with the sine term if the coefficient is non zero, - // i.e. if xmax was not infinity. Otherwise, we are evaluating the - // sine of infinity, which is NAN! Same applies to the other terms - // below. - if(term > 0.0) term *= -1/tau*sin(dm*xmax) - dm*cos(dm*xmax); - term += dm; - result += term; - } - if (basisSign != Plus) { - double term = exp(xmin/tau); - if (term > 0.0) term *= -1/tau*sin(dm*(-xmin)) - dm*cos(dm*(-xmin)); - term += dm; - result -= term; - } - return result / (1/(tau*tau) + dm*dm) ; - } - case cosBasis: - { - double result(0) ; - if (tau==0) return 1 ; - double dm = ((RooAbsReal*)basis().getParameter(2))->getVal() ; - if (basisSign != Minus) { - double term = exp(-xmax/tau); - if(term > 0.0) term *= -1/tau*cos(dm*xmax) + dm*sin(dm*xmax); - term += 1/tau; - result += term; - } - if (basisSign != Plus) { - double term = exp( xmin/tau); - if(term > 0.0) term *= -1/tau*cos(dm*(-xmin)) + dm*sin(dm*(-xmin)); - term += 1/tau; - result += term; - } - return result / (1/(tau*tau) + dm*dm) ; - } - case linBasis: - { - if (tau==0) return 0 ; - double t_max = xmax/tau ; - return tau*( 1 - (1 + t_max)*exp(-t_max) ) ; - } - case quadBasis: - { - if (tau==0) return 0 ; - double t_max = xmax/tau ; - return tau*( 2 - (2 + (2 + t_max)*t_max)*exp(-t_max) ) ; - } - case sinhBasis: - { - double result(0) ; - if (tau==0) return 0 ; - double dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ; - double taup = 2*tau/(2-tau*dg); - double taum = 2*tau/(2+tau*dg); - if (basisSign != Minus) result += 0.5*( taup*(1-exp(-xmax/taup)) - taum*(1-exp(-xmax/taum)) ) ; - if (basisSign != Plus) result -= 0.5*( taup*(1-exp( xmin/taup)) - taum*(1-exp( xmin/taum)) ) ; - return result ; - } - case coshBasis: - { - double result(0) ; - if (tau==0) return 1 ; - double dg = ((RooAbsReal*)basis().getParameter(2))->getVal() ; - double taup = 2*tau/(2-tau*dg); - double taum = 2*tau/(2+tau*dg); - if (basisSign != Minus) result += 0.5*( taup*(1-exp(-xmax/taup)) + taum*(1-exp(-xmax/taum)) ) ; - if (basisSign != Plus) result += 0.5*( taup*(1-exp( xmin/taup)) + taum*(1-exp( xmin/taum)) ) ; - return result ; - } - default: - R__ASSERT(0) ; - } + auto integrate = [&](auto indefiniteIntegral, bool isSymmetric) { + return definiteIntegral(indefiniteIntegral, xmin, xmax, tau, dm, basisSign, isSymmetric); + }; - R__ASSERT(0) ; - return 0 ; + switch (basisType) { + case expBasis: return integrate(indefiniteIntegralExpBasisPlus, /*isSymmetric=*/true); + case sinBasis: return integrate(indefiniteIntegralSinBasisPlus, /*isSymmetric=*/false); + case cosBasis: return integrate(indefiniteIntegralCosBasisPlus, /*isSymmetric=*/true); + case linBasis: return integrate(indefiniteIntegralLinBasisPlus, /*isSymmetric=*/false); + case quadBasis: return integrate(indefiniteIntegralQuadBasisPlus, /*isSymmetric=*/true); + case sinhBasis: return integrate(indefiniteIntegralSinhBasisPlus, /*isSymmetric=*/false); + case coshBasis: return integrate(indefiniteIntegralCoshBasisPlus, /*isSymmetric=*/true); + default: R__ASSERT(0); + } + + R__ASSERT(0); + return 0; } diff --git a/roofit/roofitcore/test/CMakeLists.txt b/roofit/roofitcore/test/CMakeLists.txt index ab2b4bcdd4704..3d5564fafa22c 100644 --- a/roofit/roofitcore/test/CMakeLists.txt +++ b/roofit/roofitcore/test/CMakeLists.txt @@ -52,19 +52,20 @@ if(NOT MSVC OR win_broken_tests) endif() ROOT_ADD_GTEST(testTestStatistics testTestStatistics.cxx LIBRARIES RooFitCore) endif() -ROOT_ADD_GTEST(testNaNPacker testNaNPacker.cxx LIBRARIES RooFitCore) -ROOT_ADD_GTEST(testRooSimultaneous testRooSimultaneous.cxx LIBRARIES RooFitCore) -ROOT_ADD_GTEST(testRooSTLRefCountList testRooSTLRefCountList.cxx LIBRARIES RooFitCore) +ROOT_ADD_GTEST(testGlobalObservables testGlobalObservables.cxx LIBRARIES RooFitCore) +ROOT_ADD_GTEST(testInterface TestStatistics/testInterface.cxx LIBRARIES RooFitCore) ROOT_ADD_GTEST(testLikelihoodSerial TestStatistics/testLikelihoodSerial.cxx LIBRARIES RooFitCore) +ROOT_ADD_GTEST(testNaNPacker testNaNPacker.cxx LIBRARIES RooFitCore) ROOT_ADD_GTEST(testRooAbsL TestStatistics/testRooAbsL.cxx LIBRARIES RooFitCore) -ROOT_ADD_GTEST(testRooRealL TestStatistics/testRooRealL.cxx LIBRARIES RooFitCore) -ROOT_ADD_GTEST(testInterface TestStatistics/testInterface.cxx LIBRARIES RooFitCore) -ROOT_ADD_GTEST(testGlobalObservables testGlobalObservables.cxx LIBRARIES RooFitCore) -ROOT_ADD_GTEST(testRooPolyFunc testRooPolyFunc.cxx LIBRARIES Gpad RooFitCore) -ROOT_ADD_GTEST(testSumW2Error testSumW2Error.cxx LIBRARIES Gpad RooFitCore) ROOT_ADD_GTEST(testRooHist testRooHist.cxx LIBRARIES RooFitCore) ROOT_ADD_GTEST(testRooHistPdf testRooHistPdf.cxx LIBRARIES RooFitCore) +ROOT_ADD_GTEST(testRooPolyFunc testRooPolyFunc.cxx LIBRARIES Gpad RooFitCore) +ROOT_ADD_GTEST(testRooRealL TestStatistics/testRooRealL.cxx LIBRARIES RooFitCore) ROOT_ADD_GTEST(testRooRombergIntegrator testRooRombergIntegrator.cxx LIBRARIES MathCore RooFitCore) +ROOT_ADD_GTEST(testRooSTLRefCountList testRooSTLRefCountList.cxx LIBRARIES RooFitCore) +ROOT_ADD_GTEST(testRooSimultaneous testRooSimultaneous.cxx LIBRARIES RooFitCore) +ROOT_ADD_GTEST(testRooTruthModel testRooTruthModel.cxx LIBRARIES RooFitCore RooFit) +ROOT_ADD_GTEST(testSumW2Error testSumW2Error.cxx LIBRARIES Gpad RooFitCore) if (roofit_multiprocess) ROOT_ADD_GTEST(testTestStatisticsPlot TestStatistics/testPlot.cxx LIBRARIES RooFitMultiProcess RooFitCore RooFit COPY_TO_BUILDDIR ${CMAKE_CURRENT_SOURCE_DIR}/TestStatistics/TestStatistics_ref.root) diff --git a/roofit/roofitcore/test/testRooTruthModel.cxx b/roofit/roofitcore/test/testRooTruthModel.cxx new file mode 100644 index 0000000000000..a97fcb7c4e434 --- /dev/null +++ b/roofit/roofitcore/test/testRooTruthModel.cxx @@ -0,0 +1,39 @@ +// Tests for the RooTruthModel +// Authors: Jonas Rembser, CERN 11/2023 + +#include +#include +#include +#include + +#include + +#include + +/// Check that the integration over a subrange works when using an analytical +/// convolution with the RooTruthModel. +TEST(RooTruthModel, IntegrateSubrange) +{ + using namespace RooFit; + + RooRealVar dt{"dt", "dt", 0, 10}; + + RooTruthModel truthModel{"tm", "truth model", dt}; + + RooBDecay bcpg{"bcpg0", + "bcpg0", + dt, + RooConst(1.547), + RooConst(0.323206), + RooConst(1), + RooConst(1), + RooConst(1.547), + RooConst(1.547), + RooConst(0.472), + truthModel, + RooBDecay::DecayType::SingleSided}; + dt.setRange("integral", 2, 2); + + std::unique_ptr integ{bcpg.createIntegral({dt}, "integral")}; + EXPECT_FLOAT_EQ(integ->getVal(), 0.0); +}