Skip to content

Commit

Permalink
[RF] RooTruthModel analytical integrals for general integration ranges
Browse files Browse the repository at this point in the history
The analytical integral code of the `RooTruthModel` was making the wrong
assumption that if one uses the single-sided bases, the minimum x value
is always at zero (or the maximum value at zero, for the case of the
flipped bases).

This resulted in wrong integral values when integrating over a subrange,
as reported here on the forum:

https://root-forum.cern.ch/t/possible-bug-in-integration-of-roobdecay-and-rooabsanaconvpdf/56968

This commit rewrites the RooTruthModel analytical integral code to also
consider these cases. To avoid that with the additional code branches
the code becomes too verbose, the code was refactored to use a helper
function for evaluating indefinite integrals of symmetric or asymmetric
basis functions.

The refactored code is tested by the integration tests in `stressRooFit`,
and the problem that was reported on the forum is covered by a new unit
test.
  • Loading branch information
guitargeek committed Nov 6, 2023
1 parent 2d17eab commit d12d1e0
Show file tree
Hide file tree
Showing 5 changed files with 242 additions and 187 deletions.
1 change: 0 additions & 1 deletion roofit/roofitcore/inc/RooTruthModel.h
Original file line number Diff line number Diff line change
Expand Up @@ -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 ;

Expand Down
97 changes: 44 additions & 53 deletions roofit/roofitcore/src/RooAbsAnaConvPdf.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -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<RooAbsPdf*>(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<RooAbsPdf *>(_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<RooAbsPdf*>(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<RooAbsPdf *>(_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;
}


Expand Down

0 comments on commit d12d1e0

Please sign in to comment.