@@ -306,29 +306,31 @@ inline Double_t verticalInterpPdfIntegral(double const* coefList, std::size_t nC
306306 return result > 0 . ? result : integralFloorVal;
307307}
308308
309- inline int parametricHistFindBin (const int N_bins, const double * bins, const double x)
310- {
311- if (x < bins[ 0 ] || x >= bins[N_bins]) return -1 ;
312-
309+ inline int parametricHistFindBin (const int N_bins, const double * bins, const double x) {
310+ if (x < bins[ 0 ] || x >= bins[N_bins])
311+ return -1 ;
312+
313313 // Search for the bin
314314 for (int i = 0 ; i < N_bins; ++i) {
315- if (x >= bins[i] && x < bins[i+1 ]) return i;
315+ if (x >= bins[i] && x < bins[i + 1 ])
316+ return i;
316317 }
317318 return -1 ;
318319}
319320
320- inline int parametricHistFindBin (const int N_bins, std::vector<double > const & bins, const double x)
321- {
321+ inline int parametricHistFindBin (const int N_bins, std::vector<double > const & bins, const double x) {
322322 return parametricHistFindBin (N_bins, bins.data (), x);
323323}
324324
325- inline Double_t parametricHistMorphScale (const double parVal, const int nMorphs,
325+ inline Double_t parametricHistMorphScale (const double parVal,
326+ const int nMorphs,
326327 const double * morphCoeffs,
327- const double * morphDiffs, const double * morphSums,
328- double smoothRegion)
329- {
328+ const double * morphDiffs,
329+ const double * morphSums,
330+ double smoothRegion) {
330331 double morphScale = 1.0 ;
331- if (!morphDiffs || !morphSums) return morphScale;
332+ if (!morphDiffs || !morphSums)
333+ return morphScale;
332334 for (int i = 0 ; i < nMorphs; ++i) {
333335 double coeff = morphCoeffs[i];
334336 double a = 0.5 * coeff;
@@ -338,114 +340,131 @@ inline Double_t parametricHistMorphScale(const double parVal, const int nMorphs,
338340 return morphScale;
339341}
340342
341- inline Double_t parametricHistEvaluate (const int bin_i, const double * parVals, const double * bins,
342- const int N_bins, const double * morphCoeffs, const int nMorphs,
343- const double * morphDiffs, const double * morphSums,
344- const double * widths, const double smoothRegion)
345- {
346- if (bin_i < 0 ) return 0.0 ;
347- // Morphing case
348- if (morphCoeffs != nullptr && nMorphs > 0 ) {
349- // morphDiffs and morphSums are flattened arrays of size N_bins * nMorphs
350- const double * binMorphDiffs = nullptr ;
351- const double * binMorphSums = nullptr ;
352- if (morphDiffs) {
353- binMorphDiffs = morphDiffs + bin_i * nMorphs;
354- }
355- if (morphSums) {
356- binMorphSums = morphSums + bin_i * nMorphs;
357- }
358- double parVal = parVals[bin_i];
359- double scale = parametricHistMorphScale (parVal,
360- nMorphs,
361- morphCoeffs,
362- binMorphDiffs,
363- binMorphSums,
364- smoothRegion);
365- return (parVal * scale) / widths[bin_i];
366- }
367- // No morphing case
368- return parVals[bin_i] / widths[bin_i];
343+ inline Double_t parametricHistEvaluate (const int bin_i,
344+ const double * parVals,
345+ const double * bins,
346+ const int N_bins,
347+ const double * morphCoeffs,
348+ const int nMorphs,
349+ const double * morphDiffs,
350+ const double * morphSums,
351+ const double * widths,
352+ const double smoothRegion) {
353+ if (bin_i < 0 )
354+ return 0.0 ;
355+ // Morphing case
356+ if (morphCoeffs != nullptr && nMorphs > 0 ) {
357+ // morphDiffs and morphSums are flattened arrays of size N_bins * nMorphs
358+ const double * binMorphDiffs = nullptr ;
359+ const double * binMorphSums = nullptr ;
360+ if (morphDiffs) {
361+ binMorphDiffs = morphDiffs + bin_i * nMorphs;
362+ }
363+ if (morphSums) {
364+ binMorphSums = morphSums + bin_i * nMorphs;
365+ }
366+ double parVal = parVals[bin_i];
367+ double scale = parametricHistMorphScale (parVal, nMorphs, morphCoeffs, binMorphDiffs, binMorphSums, smoothRegion);
368+ return (parVal * scale) / widths[bin_i];
369+ }
370+ // No morphing case
371+ return parVals[bin_i] / widths[bin_i];
369372}
370373
371- inline Double_t parametricMorphFunction (const int j, const double parVal, const bool hasMorphs, const int nMorphs,
372- const double * morphCoeffs, const double * morphDiffs,
373- const double * morphSums, double smoothRegion)
374- {
375- double morphScale = 1.0 ;
376- if (!hasMorphs) return morphScale;
377-
378- int ndim = nMorphs;
379- // apply all morphs one by one to the bin
380- // almost certaintly a faster way to do this in a vectorized way ....
381- for (int i = 0 ; i < ndim; ++i) {
382- double x = morphCoeffs[i];
383- double a = 0.5 *x, b = smoothStepFunc (x, smoothRegion);
384- morphScale *= 1 + (1 ./parVal) * a*(morphDiffs[j*nMorphs + i] + b*morphSums[j*nMorphs + i]);
385- }
386- return morphScale;
374+ inline Double_t parametricMorphFunction (const int j,
375+ const double parVal,
376+ const bool hasMorphs,
377+ const int nMorphs,
378+ const double * morphCoeffs,
379+ const double * morphDiffs,
380+ const double * morphSums,
381+ double smoothRegion) {
382+ double morphScale = 1.0 ;
383+ if (!hasMorphs)
384+ return morphScale;
385+
386+ int ndim = nMorphs;
387+ // apply all morphs one by one to the bin
388+ // almost certaintly a faster way to do this in a vectorized way ....
389+ for (int i = 0 ; i < ndim; ++i) {
390+ double x = morphCoeffs[i];
391+ double a = 0.5 * x, b = smoothStepFunc (x, smoothRegion);
392+ morphScale *= 1 + (1 . / parVal) * a * (morphDiffs[j * nMorphs + i] + b * morphSums[j * nMorphs + i]);
393+ }
394+ return morphScale;
387395}
388396
389- inline Double_t parametricHistFullSum (const double * parVals, const int nBins, const bool hasMorphs, const int nMorphs,
390- const double * morphCoeffs, const double * morphDiffs,
391- const double * morphSums, double smoothRegion)
392- {
393- double sum=0 ;
394- for (int i = 0 ; i < nBins; ++i) {
395- double thisVal = parVals[i];
396- if (hasMorphs) {
397- // Apply morphing to this bin, just like in RooParametricHist::evaluate
398- thisVal *= parametricMorphFunction (i, thisVal, hasMorphs, nMorphs,
399- morphCoeffs, morphDiffs, morphSums, smoothRegion);
400- }
401- sum+=thisVal;
402- }
403- return sum;
397+ inline Double_t parametricHistFullSum (const double * parVals,
398+ const int nBins,
399+ const bool hasMorphs,
400+ const int nMorphs,
401+ const double * morphCoeffs,
402+ const double * morphDiffs,
403+ const double * morphSums,
404+ double smoothRegion) {
405+ double sum = 0 ;
406+ for (int i = 0 ; i < nBins; ++i) {
407+ double thisVal = parVals[i];
408+ if (hasMorphs) {
409+ // Apply morphing to this bin, just like in RooParametricHist::evaluate
410+ thisVal *=
411+ parametricMorphFunction (i, thisVal, hasMorphs, nMorphs, morphCoeffs, morphDiffs, morphSums, smoothRegion);
412+ }
413+ sum += thisVal;
414+ }
415+ return sum;
404416}
405417
406- inline Double_t parametricHistIntegral (const double * parVals, const double * bins, const int N_bins, const double * morphCoeffs,
407- const int nMorphs, const double * morphDiffs, const double * morphSums,
408- const double * widths, const double smoothRegion, const char * rangeName,
409- const double xmin, const double xmax)
410- {
411- // No ranges
412- if (!rangeName) {
413- return parametricHistFullSum (parVals, N_bins, morphCoeffs != nullptr , nMorphs,
414- morphCoeffs, morphDiffs, morphSums, smoothRegion);
415- }
416-
417- // Case with ranges, calculate integral explicitly
418- double sum=0 ;
419- int i;
420- for (i=1 ; i<=N_bins; i++) {
421- // Get maybe-morphed bin value
422- double binVal = parVals[i-1 ] / widths[i-1 ];
423- if (morphCoeffs != nullptr ) {
424- binVal *= parametricMorphFunction (i - 1 , parVals[i - 1 ], true , nMorphs,
425- morphCoeffs, morphDiffs, morphSums, smoothRegion);
426- }
418+ inline Double_t parametricHistIntegral (const double * parVals,
419+ const double * bins,
420+ const int N_bins,
421+ const double * morphCoeffs,
422+ const int nMorphs,
423+ const double * morphDiffs,
424+ const double * morphSums,
425+ const double * widths,
426+ const double smoothRegion,
427+ const char * rangeName,
428+ const double xmin,
429+ const double xmax) {
430+ // No ranges
431+ if (!rangeName) {
432+ return parametricHistFullSum (
433+ parVals, N_bins, morphCoeffs != nullptr , nMorphs, morphCoeffs, morphDiffs, morphSums, smoothRegion);
434+ }
427435
428- if (bins[i-1 ] >= xmin && bins[i] <= xmax) {
429- // Bin fully in integration domain
430- sum += (bins[i] - bins[i-1 ]) * binVal;
431- } else if (bins[i-1 ] < xmin && bins[i] > xmax) {
432- // Domain is fully contained in this bin
433- sum += (xmax - xmin) * binVal;
434- // Exit here, this is the last bin to be processed by construction
435- double fullSum = parametricHistFullSum (parVals, N_bins, morphCoeffs != nullptr , nMorphs,
436- morphCoeffs, morphDiffs, morphSums, smoothRegion);
437- return sum / fullSum;
438- } else if (bins[i-1 ] < xmin && bins[i] <= xmax && bins[i] > xmin) {
439- // Lower domain boundary is in bin
440- sum += (bins[i] - xmin) * binVal;
441- } else if (bins[i-1 ] >= xmin && bins[i] > xmax && bins[i-1 ] < xmax) {
442- // Upper domain boundary is in bin
443- // Exit here, this is the last bin to be processed by construction
444- sum += (xmax - bins[i-1 ]) * binVal;
445- return sum;
446- }
447- }
448- return sum;
436+ // Case with ranges, calculate integral explicitly
437+ double sum = 0 ;
438+ int i;
439+ for (i = 1 ; i <= N_bins; i++) {
440+ // Get maybe-morphed bin value
441+ double binVal = parVals[i - 1 ] / widths[i - 1 ];
442+ if (morphCoeffs != nullptr ) {
443+ binVal *= parametricMorphFunction (
444+ i - 1 , parVals[i - 1 ], true , nMorphs, morphCoeffs, morphDiffs, morphSums, smoothRegion);
445+ }
446+
447+ if (bins[i - 1 ] >= xmin && bins[i] <= xmax) {
448+ // Bin fully in integration domain
449+ sum += (bins[i] - bins[i - 1 ]) * binVal;
450+ } else if (bins[i - 1 ] < xmin && bins[i] > xmax) {
451+ // Domain is fully contained in this bin
452+ sum += (xmax - xmin) * binVal;
453+ // Exit here, this is the last bin to be processed by construction
454+ double fullSum = parametricHistFullSum (
455+ parVals, N_bins, morphCoeffs != nullptr , nMorphs, morphCoeffs, morphDiffs, morphSums, smoothRegion);
456+ return sum / fullSum;
457+ } else if (bins[i - 1 ] < xmin && bins[i] <= xmax && bins[i] > xmin) {
458+ // Lower domain boundary is in bin
459+ sum += (bins[i] - xmin) * binVal;
460+ } else if (bins[i - 1 ] >= xmin && bins[i] > xmax && bins[i - 1 ] < xmax) {
461+ // Upper domain boundary is in bin
462+ // Exit here, this is the last bin to be processed by construction
463+ sum += (xmax - bins[i - 1 ]) * binVal;
464+ return sum;
465+ }
466+ }
467+ return sum;
449468}
450469
451470} // namespace MathFuncs
0 commit comments