-
Notifications
You must be signed in to change notification settings - Fork 393
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
Merge pull request #1017 from guitargeek/combine_math_funcs
Factor out some mathematical details into free functions
- Loading branch information
Showing
7 changed files
with
200 additions
and
99 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,152 @@ | ||
| #ifndef CombineMathFuncs_h | ||
| #define CombineMathFuncs_h | ||
|
|
||
| #include <cmath> | ||
|
|
||
| namespace RooFit { | ||
| namespace Detail { | ||
| namespace MathFuncs { | ||
|
|
||
| inline double smoothStepFunc(double x, double smoothRegion) | ||
| { | ||
| if (std::abs(x) >= smoothRegion) | ||
| return x > 0 ? +1 : -1; | ||
| double xnorm = x / smoothRegion; | ||
| double xnorm2 = xnorm * xnorm; | ||
| return 0.125 * xnorm * (xnorm2 * (3. * xnorm2 - 10.) + 15); | ||
| } | ||
|
|
||
| inline void fastVerticalInterpHistPdf2(int nBins, int nCoefs, double const *coefs, double const *nominal, | ||
| double const *binWidth, double const *morphsSum, double const *morphsDiff, | ||
| double smoothRegion, double *out) | ||
| { | ||
| for (int iBin = 0; iBin < nBins; ++iBin) { | ||
| out[iBin] = nominal[iBin]; | ||
| } | ||
|
|
||
| double normSum = 0.0; | ||
|
|
||
| for (int iBin = 0; iBin < nBins; ++iBin) { | ||
| // apply all morphs one by one | ||
| for (int iCoef = 0; iCoef < nCoefs; ++iCoef) { | ||
| double const *sum = morphsSum + iCoef * nBins; | ||
| double const *diff = morphsDiff + iCoef * nBins; | ||
| double x = coefs[iCoef]; | ||
| double a = 0.5 * x; | ||
| double b = smoothStepFunc(x, smoothRegion); | ||
|
|
||
| out[iBin] += a * (diff[iBin] + b * sum[iBin]); | ||
| } | ||
|
|
||
| out[iBin] = std::max(1e-9, out[iBin]); | ||
| normSum += out[iBin] * binWidth[iBin]; | ||
| } | ||
|
|
||
| if (normSum > 0.0) { | ||
| double normSumInv = 1. / normSum; | ||
| for (int iBin = 0; iBin < nBins; ++iBin) { | ||
| out[iBin] *= normSumInv; | ||
| } | ||
| } | ||
| } | ||
|
|
||
| inline void fastVerticalInterpHistPdf2D2(int nBinsX, int nBinsY, int nCoefs, double const *coefs, | ||
| double const *nominal, double const *binWidth, double const *morphsSum, | ||
| double const *morphsDiff, double smoothRegion, double *out) | ||
| { | ||
| int nBins = nBinsX * nBinsY; | ||
|
|
||
| for (int iBin = 0; iBin < nBins; ++iBin) { | ||
| out[iBin] = nominal[iBin]; | ||
| } | ||
|
|
||
| for (int iBinX = 0; iBinX < nBinsX; ++iBinX) { | ||
|
|
||
| double normSum = 0.0; | ||
|
|
||
| for (int iBinY = 0; iBinY < nBinsY; ++iBinY) { | ||
| int iBin = iBinY + nBinsY * iBinX; | ||
| // apply all morphs one by one | ||
| for (int iCoef = 0; iCoef < nCoefs; ++iCoef) { | ||
| double const *sum = morphsSum + iCoef * nBinsY; | ||
| double const *diff = morphsDiff + iCoef * nBinsY; | ||
| double x = coefs[iCoef]; | ||
| double a = 0.5 * x; | ||
| double b = smoothStepFunc(x, smoothRegion); | ||
|
|
||
| out[iBin] += a * (diff[iBin] + b * sum[iBin]); | ||
| } | ||
|
|
||
| out[iBin] = std::max(1e-9, out[iBin]); | ||
| normSum += out[iBin] * binWidth[iBin]; | ||
| } | ||
|
|
||
| if (normSum > 0.0) { | ||
| double normSumInv = 1. / normSum; | ||
| for (int iBinY = 0; iBinY < nBinsY; ++iBinY) { | ||
| int iBin = iBinY + nBinsY * iBinX; | ||
| out[iBin] *= normSumInv; | ||
| } | ||
| } | ||
| } | ||
| } | ||
|
|
||
| inline double logKappaForX(double theta, double logKappaLow, double logKappaHigh) | ||
| { | ||
| double logKappa = 0.0; | ||
|
|
||
| if (std::abs(theta) >= 0.5) { | ||
| logKappa = theta >= 0 ? logKappaHigh : -logKappaLow; | ||
| } else { | ||
| // interpolate between log(kappaHigh) and -log(kappaLow) | ||
| // logKappa(x) = avg + halfdiff * h(2x) | ||
| // where h(x) is the 3th order polynomial | ||
| // h(x) = (3 x^5 - 10 x^3 + 15 x)/8; | ||
| // chosen so that h(x) satisfies the following: | ||
| // h (+/-1) = +/-1 | ||
| // h'(+/-1) = 0 | ||
| // h"(+/-1) = 0 | ||
| double logKhi = logKappaHigh; | ||
| double logKlo = -logKappaLow; | ||
| double avg = 0.5 * (logKhi + logKlo); | ||
| double halfdiff = 0.5 * (logKhi - logKlo); | ||
| double twox = theta + theta; | ||
| double twox2 = twox * twox; | ||
| double alpha = 0.125 * twox * (twox2 * (3 * twox2 - 10.) + 15.); | ||
| logKappa = avg + alpha * halfdiff; | ||
| } | ||
|
|
||
| return logKappa; | ||
| } | ||
|
|
||
| inline double asymPow(double theta, double kappaLow, double kappaHigh) | ||
| { | ||
| return std::exp(logKappaForX(theta, std::log(kappaLow), std::log(kappaHigh)) * theta); | ||
| } | ||
|
|
||
| inline double processNormalization(double nominalValue, std::size_t nThetas, std::size_t nAsymmThetas, | ||
| std::size_t nOtherFactors, double const *thetas, double const *logKappas, | ||
| double const *asymmThetas, double const *asymmLogKappasLow, | ||
| double const *asymmLogKappasHigh, double const *otherFactors) | ||
| { | ||
| double logVal = 0.0; | ||
| for (std::size_t i = 0; i < nThetas; i++) { | ||
| logVal += thetas[i] * logKappas[i]; | ||
| } | ||
| for (std::size_t i = 0; i < nAsymmThetas; i++) { | ||
| double x = asymmThetas[i]; | ||
| logVal += x * logKappaForX(x, asymmLogKappasLow[i], asymmLogKappasHigh[i]); | ||
| } | ||
| double norm = nominalValue; | ||
| norm *= std::exp(logVal); | ||
| for (std::size_t i = 0; i < nOtherFactors; i++) { | ||
| norm *= otherFactors[i]; | ||
| } | ||
| return norm; | ||
| } | ||
|
|
||
| } // namespace MathFuncs | ||
| } // namespace Detail | ||
| } // namespace RooFit | ||
|
|
||
| #endif |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters