1#ifndef RIVET_RIVETYODA_HH
2#define RIVET_RIVETYODA_HH
4#include "Rivet/Config/RivetCommon.hh"
5#include "Rivet/Tools/TypeTraits.hh"
6#include "YODA/AnalysisObject.h"
7#include "YODA/Counter.h"
9#include "YODA/Profile.h"
10#include "YODA/Estimate0D.h"
11#include "YODA/BinnedEstimate.h"
12#include "YODA/Scatter.h"
20#include <unordered_map>
26 template<
size_t DbnN,
typename ... AxisT>
27 using BinnedDbnPtr = std::shared_ptr<YODA::BinnedDbn<DbnN, AxisT...>>;
29 template<
typename ... AxisT>
30 using BinnedHistoPtr = BinnedDbnPtr<
sizeof...(AxisT), AxisT...>;
32 template<
typename ... AxisT>
33 using BinnedProfilePtr = BinnedDbnPtr<
sizeof...(AxisT)+1, AxisT...>;
35 template<
typename ... AxisT>
36 using BinnedEstimatePtr = std::shared_ptr<YODA::BinnedEstimate<AxisT...>>;
39 using ScatterNDPtr = std::shared_ptr<YODA::ScatterND<N>>;
41 using AnalysisObjectPtr = std::shared_ptr<YODA::AnalysisObject>;
42 using CounterPtr = std::shared_ptr<YODA::Counter>;
43 using Estimate0DPtr = std::shared_ptr<YODA::Estimate0D>;
44 using Histo1DPtr = BinnedHistoPtr<double>;
45 using Histo2DPtr = BinnedHistoPtr<double,double>;
46 using Histo3DPtr = BinnedHistoPtr<double,double,double>;
47 using Profile1DPtr = BinnedProfilePtr<double>;
48 using Profile2DPtr = BinnedProfilePtr<double,double>;
49 using Profile3DPtr = BinnedProfilePtr<double,double,double>;
50 using Estimate1DPtr = BinnedEstimatePtr<double>;
51 using Estimate2DPtr = BinnedEstimatePtr<double,double>;
52 using Estimate3DPtr = BinnedEstimatePtr<double,double,double>;
53 using Scatter1DPtr = ScatterNDPtr<1>;
54 using Scatter2DPtr = ScatterNDPtr<2>;
55 using Scatter3DPtr = ScatterNDPtr<3>;
64 bool copyAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst,
const double scale=1.0) {
65 if (dst->hasAnnotation(
"Type") && src->type() != dst->type()) {
66 throw YODA::LogicError(
"Operation requries types to be the same!");
68 for (
const std::string& a : src->annotations()) {
69 dst->setAnnotation(a, src->annotation(a));
71 shared_ptr<T> dstPtr = std::static_pointer_cast<T>(dst);
72 *dstPtr = *std::static_pointer_cast<T>(src);
73 if constexpr (isFillable<T>::value) { dstPtr->scaleW(scale); }
85 virtual bool copyAO(YODA::AnalysisObjectPtr src,
86 YODA::AnalysisObjectPtr dst,
87 const double scale = 1.0)
const = 0;
89 virtual bool addAO(YODA::AnalysisObjectPtr src,
90 YODA::AnalysisObjectPtr& dst,
91 const double scale = 1.0)
const = 0;
102 bool addAO(YODA::AnalysisObjectPtr src,
103 YODA::AnalysisObjectPtr& dst,
104 const double scale = 1.0)
const {
105 if constexpr (isFillable<T>::value) {
106 std::shared_ptr<T> srcPtr = std::static_pointer_cast<T>(src);
107 srcPtr->scaleW(scale);
108 if (dst ==
nullptr) { dst = src;
return true; }
109 try { *std::static_pointer_cast<T>(dst) += *srcPtr; }
110 catch (YODA::BinningError&) {
return false; }
113 else if (dst ==
nullptr) { dst = src;
return true; }
117 bool copyAO(YODA::AnalysisObjectPtr src,
118 YODA::AnalysisObjectPtr dst,
119 const double scale = 1.0)
const {
120 return ::Rivet::copyAO<T>(src, dst, scale);
139 using Fill = pair<typename T::FillType, Weight>;
166 using YAO = YODA::Counter;
167 using Ptr = shared_ptr<FillCollector<YAO>>;
168 using YAO::operator =;
185 int fill(
const double weight=1.0,
const double fraction = 1.0) {
187 _fills.insert(_fills.end(), { YAO::FillType(), weight } );
205 template <
size_t DbnN,
typename... AxisT>
207 :
public YODA::BinnedDbn<DbnN, AxisT...> {
210 using YAO = YODA::BinnedDbn<DbnN, AxisT...>;
211 using Ptr = shared_ptr<FillCollector<YAO>>;
212 using YAO::operator =;
225 YAO::setPath(yao->path());
232 int fill(
typename YAO::FillType&& fillCoords,
233 const double weight=1.0,
const double fraction=1.0) {
235 if (YODA::containsNan(fillCoords)) {
236 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
241 typename YAO::BinningT::EdgeTypesTuple binCoords{};
242 auto extractBinCoords = [&binCoords, &fillCoords](
auto I) {
243 std::get<I>(binCoords) = std::get<I>(fillCoords);
245 MetaUtils::staticFor<
sizeof...(AxisT)>(extractBinCoords);
246 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
247 return (
int)YAO::_binning.globalIndexAt(binCoords);
251 void reset() noexcept { _fills.clear(); }
262 template <
typename AxisT>
264 :
public YODA::BinnedDbn<1, AxisT> {
267 using YAO = YODA::BinnedDbn<1, AxisT>;
268 using Ptr = shared_ptr<FillCollector<YAO>>;
269 using YAO::operator =;
275 YAO::setPath(yao->path());
282 int fill(
typename YAO::FillType&& fillCoords,
283 const double weight=1.0,
const double fraction=1.0) {
285 if (YODA::containsNan(fillCoords)) {
286 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
291 typename YAO::BinningT::EdgeTypesTuple binCoords{};
292 auto extractBinCoords = [&binCoords, &fillCoords](
auto I) {
293 std::get<I>(binCoords) = std::get<I>(fillCoords);
295 MetaUtils::staticFor<1>(extractBinCoords);
296 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
297 return (
int)YAO::_binning.globalIndexAt(binCoords);
300 int fill(
const AxisT x,
const double weight=1.0,
const double fraction=1.0) {
301 return fill(
typename YAO::FillType{x}, weight, fraction);
305 void reset() noexcept { _fills.clear(); }
316 template <
typename AxisT1,
typename AxisT2>
318 :
public YODA::BinnedDbn<2, AxisT1, AxisT2> {
321 using YAO = YODA::BinnedDbn<2, AxisT1, AxisT2>;
322 using Ptr = shared_ptr<FillCollector<YAO>>;
323 using YAO::operator =;
329 YAO::setPath(yao->path());
336 int fill(
typename YAO::FillType&& fillCoords,
337 const double weight=1.0,
const double fraction=1.0) {
339 if (YODA::containsNan(fillCoords)) {
340 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
345 typename YAO::BinningT::EdgeTypesTuple binCoords{};
346 auto extractBinCoords = [&binCoords, &fillCoords](
auto I) {
347 std::get<I>(binCoords) = std::get<I>(fillCoords);
349 MetaUtils::staticFor<1>(extractBinCoords);
350 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
351 return (
int)YAO::_binning.globalIndexAt(binCoords);
354 int fill(
const AxisT1 x,
const AxisT2 y,
const double weight=1.0,
const double fraction=1.0) {
355 return fill(
typename YAO::FillType{x,y}, weight, fraction);
359 void reset() noexcept { _fills.clear(); }
370 template <
typename AxisT1,
typename AxisT2,
typename AxisT3>
372 :
public YODA::BinnedDbn<3, AxisT1, AxisT2, AxisT3> {
375 using YAO = YODA::BinnedDbn<3, AxisT1, AxisT2, AxisT3>;
376 using Ptr = shared_ptr<FillCollector<YAO>>;
377 using YAO::operator =;
383 YAO::setPath(yao->path());
390 int fill(
typename YAO::FillType&& fillCoords,
391 const double weight=1.0,
const double fraction=1.0) {
393 if (YODA::containsNan(fillCoords)) {
394 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
399 typename YAO::BinningT::EdgeTypesTuple binCoords{};
400 auto extractBinCoords = [&binCoords, &fillCoords](
auto I) {
401 std::get<I>(binCoords) = std::get<I>(fillCoords);
403 MetaUtils::staticFor<1>(extractBinCoords);
404 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
405 return (
int)YAO::_binning.globalIndexAt(binCoords);
408 int fill(
const AxisT1 x,
const AxisT2 y,
const AxisT3 z,
const double weight=1.0,
const double fraction=1.0) {
409 return fill(
typename YAO::FillType{x,y,z}, weight, fraction);
413 void reset() noexcept { _fills.clear(); }
424 template <
typename AxisT>
426 :
public YODA::BinnedDbn<2, AxisT> {
429 using YAO = YODA::BinnedDbn<2, AxisT>;
430 using Ptr = shared_ptr<FillCollector<YAO>>;
431 using YAO::operator =;
437 YAO::setPath(yao->path());
444 int fill(
typename YAO::FillType&& fillCoords,
445 const double weight=1.0,
const double fraction=1.0) {
447 if (YODA::containsNan(fillCoords)) {
448 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
453 typename YAO::BinningT::EdgeTypesTuple binCoords{};
454 auto extractBinCoords = [&binCoords, &fillCoords](
auto I) {
455 std::get<I>(binCoords) = std::get<I>(fillCoords);
457 MetaUtils::staticFor<1>(extractBinCoords);
458 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
459 return (
int)YAO::_binning.globalIndexAt(binCoords);
462 int fill(
const AxisT x,
const double y,
const double weight=1.0,
const double fraction=1.0) {
463 return fill(
typename YAO::FillType{x,y}, weight, fraction);
467 void reset() noexcept { _fills.clear(); }
478 template <
typename AxisT1,
typename AxisT2>
480 :
public YODA::BinnedDbn<3, AxisT1, AxisT2> {
483 using YAO = YODA::BinnedDbn<3, AxisT1, AxisT2>;
484 using Ptr = shared_ptr<FillCollector<YAO>>;
485 using YAO::operator =;
491 YAO::setPath(yao->path());
498 int fill(
typename YAO::FillType&& fillCoords,
499 const double weight=1.0,
const double fraction=1.0) {
501 if (YODA::containsNan(fillCoords)) {
502 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
507 typename YAO::BinningT::EdgeTypesTuple binCoords{};
508 auto extractBinCoords = [&binCoords, &fillCoords](
auto I) {
509 std::get<I>(binCoords) = std::get<I>(fillCoords);
511 MetaUtils::staticFor<1>(extractBinCoords);
512 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
513 return (
int)YAO::_binning.globalIndexAt(binCoords);
516 int fill(
const AxisT1 x,
const AxisT2 y,
const double z,
const double weight=1.0,
const double fraction=1.0) {
517 return fill(
typename YAO::FillType{x,y,z}, weight, fraction);
521 void reset() noexcept { _fills.clear(); }
532 template <
typename AxisT1,
typename AxisT2,
typename AxisT3>
534 :
public YODA::BinnedDbn<4, AxisT1, AxisT2, AxisT3> {
537 using YAO = YODA::BinnedDbn<4, AxisT1, AxisT2, AxisT3>;
538 using Ptr = shared_ptr<FillCollector<YAO>>;
539 using YAO::operator =;
545 YAO::setPath(yao->path());
552 int fill(
typename YAO::FillType&& fillCoords,
553 const double weight=1.0,
const double fraction=1.0) {
555 if (YODA::containsNan(fillCoords)) {
556 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
561 typename YAO::BinningT::EdgeTypesTuple binCoords{};
562 auto extractBinCoords = [&binCoords, &fillCoords](
auto I) {
563 std::get<I>(binCoords) = std::get<I>(fillCoords);
565 MetaUtils::staticFor<1>(extractBinCoords);
566 _fills.insert(_fills.end(), { std::move(fillCoords), weight } );
567 return (
int)YAO::_binning.globalIndexAt(binCoords);
570 int fill(
const AxisT1 x,
const AxisT2 y,
const AxisT3 z,
const double zPlus,
const double weight=1.0,
const double fraction=1.0) {
571 return fill(
typename YAO::FillType{x,y,z,zPlus}, weight, fraction);
575 void reset() noexcept { _fills.clear(); }
592 using YAO = YODA::Estimate0D;
593 using Ptr = shared_ptr<FillCollector<YAO>>;
594 using YAO::operator =;
610 template<
typename ... AxisT>
612 :
public YODA::BinnedEstimate<AxisT...> {
615 using YAO = YODA::BinnedEstimate<AxisT...>;
616 using Ptr = shared_ptr<FillCollector<YAO>>;
617 using YAO::operator =;
637 using YAO = YODA::ScatterND<N>;
638 using Ptr = shared_ptr<FillCollector<YAO>>;
639 using YAO::operator =;
661 template<
typename... Args>
662 double distance(
const tuple<Args...>& a,
const tuple<Args...>& b) {
664 auto calculateDistance = [&](
auto I) {
665 if constexpr (std::is_floating_point<std::tuple_element_t<I,
666 std::tuple<Args...>>>::value) {
670 MetaUtils::staticFor<
sizeof...(Args)>(calculateDistance);
680 template <
typename T>
681 vector<Fills<T>> applyEmptyFillPaddingAndTranspose(
const vector<
typename FillCollector<T>::Ptr>& subevents) {
683 static Fill<T> EmptyFill = {
typename T::FillType{}, 0.0 };
685 vector<Fills<T>> matched;
686 matched.reserve(subevents.size());
688 size_t maxFillLen = 0;
689 size_t maxFillPos = 0;
690 for (
const auto& subevt : subevents) {
691 auto&& fills = subevt->fills();
692 if (fills.size() > maxFillLen) {
693 maxFillLen = fills.size();
694 maxFillPos = matched.size();
696 matched.push_back(std::move(fills));
700 const Fills<T>& maxFill = matched[maxFillPos];
701 for (
auto& fills : matched) {
703 if (fills.size() == maxFillLen)
continue;
707 while (fills.size() < maxFillLen) {
708 fills.push_back(EmptyFill);
714 for (
int i = maxFillLen - 1; i >= 0; --i) {
715 if (fills[i] == EmptyFill)
continue;
717 while (j+1 < (
int)maxFillLen &&
718 fills[j + 1] == EmptyFill &&
719 distance(fills[j].first,
721 distance(fills[j].first,
722 maxFill[j+1].first)) {
723 swap(fills[j], fills[j+1]);
729 vector<Fills<T>> result(maxFillLen, Fills<T>(matched.size()));
730 for (
size_t i = 0; i < matched.size(); ++i) {
731 for (
size_t j = 0; j < maxFillLen; ++j) {
732 result.at(j).at(i) = std::move(matched.at(i).at(j));
741 struct SubwindowType;
743 template<
typename... Args>
744 struct SubwindowType<tuple<Args...>> {
745 using type = YODA::Binning<std::decay_t<
decltype(std::declval<YODA::Axis<Args>>())>...>;
755 template<
typename YAO>
756 vector<std::tuple<typename YAO::FillType, valarray<double>,
double>>
757 applyFillWindows(shared_ptr<YAO> ao,
const Fills<YAO>& subevents,
758 const vector<valarray<double>>& weights,
const double fsmear=0.5) {
761 using SubwindowT =
typename SubwindowType<typename YAO::FillType>::type;
762 SubwindowT subwindows;
764 const size_t nFills = subevents.size();
765 vector<vector<double>> edgesLo, edgesHi;
766 edgesLo.resize(YAO::FillDimension::value);
767 edgesHi.resize(YAO::FillDimension::value);
769 auto constructWindows = [&](
auto I) {
771 using FillAxisT =
typename SubwindowT::template getAxisT<I>;
772 using isContinuous =
typename SubwindowT::template is_CAxis<I>;
774 if constexpr(!isContinuous::value) {
776 subwindows.template axis<I>() = FillAxisT({ std::get<I>(subevents[0].first) });
780 edgesHi[I].resize(nFills);
781 edgesLo[I].resize(nFills);
783 if constexpr(I < YAO::BinningT::Dimension::value) {
785 const auto& axis = ao->binning().template axis<I>();
786 size_t over = 0, under = 0;
787 const double edgeMax = ao->template max<I>();
788 const double edgeMin = ao->template min<I>();
789 const size_t binLast = axis.numBins();
791 for (
size_t i = 0; i < nFills; ++i) {
792 const double edge = get<I>(subevents[i].first);
793 size_t idx = axis.index(edge);
794 if (edge >= edgeMax) {
795 if (edge > edgeMax) ++over;
798 else if (edge < edgeMin) {
804 if (edge > axis.mid(idx)) {
805 if (idx != binLast) ++ibn;
812 const double ibw = axis.width(idx) < axis.width(ibn)? idx : ibn;
813 if ( fsmear > 0.0 ) {
814 const double delta = 0.5*fsmear*axis.width(ibw);
815 edgesHi[I][i] = edge + delta;
816 edgesLo[I][i] = edge - delta;
819 const double delta = 0.5*axis.width(ibw);
820 if (edge > edgeMax) {
821 edgesHi[I][i] =
max(edgeMax + 2*delta, edge + delta);
822 edgesLo[I][i] =
max(edgeMax, edge - delta);
824 else if (edge < edgeMin) {
825 edgesHi[I][i] =
min(edgeMin, edge + delta);
826 edgesLo[I][i] =
min(edgeMin - 2*delta, edge - delta);
829 edgesHi[I][i] = axis.max(idx);
830 edgesLo[I][i] = axis.min(idx);
834 for (
size_t i = 0; i < nFills; ++i) {
835 const double wsize = edgesHi[I][i] - edgesLo[I][i];
836 if (over == nFills && edgesLo[I][i] < edgeMax && edgesHi[I][i] > edgeMax) {
837 edgesHi[I][i] = edgeMax + wsize;
838 edgesLo[I][i] = edgeMax;
840 else if (over == 0 && edgesLo[I][i] < edgeMax && edgesHi[I][i] > edgeMax) {
841 edgesLo[I][i] = edgeMax - wsize;
842 edgesHi[I][i] = edgeMax;
844 else if (under == nFills && edgesLo[I][i] < edgeMin && edgesHi[I][i] > edgeMin) {
845 edgesLo[I][i] = edgeMin - wsize;
846 edgesHi[I][i] = edgeMin;
848 else if (under == 0 && edgesLo[I][i] < edgeMin && edgesHi[I][i] > edgeMin) {
849 edgesHi[I][i] = edgeMin + wsize;
850 edgesLo[I][i] = edgeMin;
856 for (
size_t i = 0; i < nFills; ++i) {
859 const double edge = get<I>(subevents[i].first);
860 const double delta = 0.1*fabs(edge);
861 edgesHi[I][i] = edge + delta;
862 edgesLo[I][i] = edge - delta;
866 vector<double> windowEdges;
867 std::copy(edgesLo[I].begin(), edgesLo[I].end(), std::back_inserter(windowEdges));
868 std::copy(edgesHi[I].begin(), edgesHi[I].end(), std::back_inserter(windowEdges));
869 std::sort(windowEdges.begin(), windowEdges.end());
870 windowEdges.erase( std::unique(windowEdges.begin(), windowEdges.end()), windowEdges.end() );
871 subwindows.template axis<I>() = FillAxisT(windowEdges);
875 MetaUtils::staticFor<YAO::FillDimension::value>(constructWindows);
878 vector<std::tuple<typename YAO::FillType, valarray<double>,
double>> rtn;
882 const vector<size_t> overflows = subwindows.calcOverflowBinsIndices();
883 const auto& itEnd = overflows.cend();
884 for (
size_t i = 0; i < subwindows.numBins(); ++i) {
885 if (std::find(overflows.cbegin(), itEnd, i) != itEnd)
continue;
887 const auto coords = subwindows.edgeTuple(i);
888 const double subwindowArea = subwindows.dVol(i);
889 size_t nSubfills = 0;
890 double windowFrac = 0.;
891 valarray<double> sumw(0.0, weights[0].size());
892 for (
size_t j = 0; j < nFills; ++j) {
894 double windowArea = 1.0;
895 auto checkSubwindowOverlap = [&](
auto I) {
896 using isContinuous =
typename SubwindowT::template is_CAxis<I>;
897 if constexpr (isContinuous::value) {
898 const double edge = std::get<I>(coords);
899 pass &= (edgesLo[I][j] <= edge && edge <= edgesHi[I][j]);
900 windowArea *= edgesHi[I][j] - edgesLo[I][j];
903 MetaUtils::staticFor<YAO::FillDimension::value>(checkSubwindowOverlap);
905 windowFrac = subwindowArea/windowArea;
906 sumw += subevents[j].second * weights[j];
911 const double fillFrac = (double)nSubfills/(
double)nFills;
912 rtn.emplace_back(coords, sumw/fillFrac, fillFrac*windowFrac);
943 virtual void collapseEventGroup(
const vector<std::valarray<double>>& weight,
const double nlowfrac=0.0) = 0;
949 virtual YODA::AnalysisObjectPtr
activeAO()
const = 0;
981 const vector<bool>& fillOutcomes()
const {
return _fillOutcomes; }
983 const vector<double>& fillFractions()
const {
return _fillFractions; }
987 vector<bool> _fillOutcomes;
989 vector<double> _fillFractions;
1010 template<
typename T>
1020 using MultiplexedAO::_fillOutcomes;
1021 using MultiplexedAO::_fillFractions;
1025 Multiplexer(
const vector<string>& weightNames,
const T& p) {
1026 _basePath = p.path();
1027 _baseName = p.name();
1028 for (
const string& weightname : weightNames) {
1029 _persistent.push_back(make_shared<T>(p));
1030 _final.push_back(make_shared<T>(p));
1032 typename T::Ptr obj = _persistent.back();
1033 obj->setPath(
"/RAW" + obj->path());
1034 typename T::Ptr
final = _final.back();
1035 if (weightname !=
"") {
1036 obj->setPath(obj->path() +
"[" + weightname +
"]");
1037 final->setPath(
final->path() +
"[" + weightname +
"]");
1048 #ifdef HAVE_BACKTRACE
1050 backtrace(buffer, 4);
1051 backtrace_symbols_fd(buffer, 4 , 1);
1053 assert(
false &&
"No active pointer set. Was this object booked in init()?");
1058 template <
typename U = T>
1059 auto binning() const ->
std::enable_if_t<hasBinning<T>::value, const typename U::BinningT&> {
1060 return _persistent.back()->binning();
1073 explicit operator bool()
const {
return static_cast<bool>(_active); }
1096 if (a._persistent.size() != b._persistent.size()) {
1100 for (
size_t i = 0; i < a._persistent.size(); ++i) {
1101 if (a._persistent.at(i) != b._persistent.at(i)) {
1115 if (a._persistent.size() >= b._persistent.size()) {
1118 for (
size_t i = 0; i < a._persistent.size(); ++i) {
1119 if (*(a._persistent.at(i)) >= *(b._persistent.at(i))) {
1134 _active = _persistent.at(iWeight);
1139 _active = _final.at(iWeight);
1156 _active = _evgroup.back();
1164 if constexpr( isFillable<T>::value ) {
1167 assert( _evgroup.size() == weights.size() );
1169 if (_evgroup.size() == 1) {
1173 std::fill(_fillOutcomes.begin(), _fillOutcomes.end(),
false);
1174 std::fill(_fillFractions.begin(), _fillFractions.end(), 0.0);
1178 for (
auto f : _evgroup[0]->fills()) {
1181 for (
size_t m = 0; m < _persistent.size(); ++m) {
1182 if (!m) frac = f.second;
1183 pos = _persistent[m]->fill( std::move(f.first), std::move(f.second) * weights[0][m] );
1186 _fillOutcomes[pos] =
true;
1187 _fillFractions[pos] += frac;
1204 if constexpr( isFillable<T>::value ) {
1205 if constexpr (!std::is_same<T, YODA::Counter>::value ) {
1211 for (
const Fills<T>& subEvents : applyEmptyFillPaddingAndTranspose<T>(_evgroup)) {
1213 for (
const auto& f : applyFillWindows(_persistent[0], subEvents, weights, nlowfrac)) {
1214 for (
size_t m = 0; m < _persistent.size(); ++m ) {
1215 _persistent[m]->fill(
typename T::FillType(get<0>(f)), get<1>(f)[m], get<2>(f) );
1221 for (
size_t m = 0; m < _persistent.size(); ++m) {
1222 vector<double> sumfw{0.0};
1223 for (
size_t n = 0; n < _evgroup.size(); ++n) {
1224 const auto& fills = _evgroup[n]->fills();
1227 if (fills.size() > sumfw.size()) {
1228 sumfw.resize(fills.size(), 0.0);
1231 for (
const auto& f : fills) {
1232 sumfw[fi++] += f.second * weights[n][m];
1235 for (
double fw : sumfw) {
1236 _persistent[m]->fill(std::move(fw));
1246 for (
size_t m = 0; m < _persistent.size(); ++m ) {
1247 _final.at(m)->clearAnnotations();
1248 copyAO<T>(_persistent.at(m), _final.at(m));
1250 if ( _final[m]->path().substr(0,4) ==
"/RAW" )
1251 _final[m]->setPath(_final[m]->path().substr(4));
1262 typename T::Ptr
persistent(
const size_t iWeight) {
return _persistent.at(iWeight); }
1266 const vector<typename T::Ptr>&
final()
const {
1271 typename T::Ptr
final(
const size_t iWeight) {
return _final.at(iWeight); }
1274 YODA::AnalysisObjectPtr
activeAO()
const {
return _active; }
1278 if constexpr( isFillable<T>::value ) {
1280 if constexpr (!std::is_same<T, YODA::Counter>::value ) {
1281 nPos = _persistent.back()->numBins(
true,
true);
1283 _fillOutcomes.resize(nPos);
1284 _fillFractions.resize(nPos);
1296 vector<typename T::Ptr> _persistent;
1299 vector<typename T::Ptr> _final;
1302 vector<typename FillCollector<T>::Ptr> _evgroup;
1306 typename T::Ptr _active;
1329 template <
typename T>
1334 using value_type = T;
1341 MultiplexPtr(
const vector<string>& weightNames,
const typename T::Inner& p)
1342 : _p( make_shared<T>(weightNames, p) ) { }
1345 template <typename U, typename = decltype(shared_ptr<T>(shared_ptr<U>{}))>
1346 MultiplexPtr(
const shared_ptr<U>& p) : _p(p) { }
1349 template <typename U, typename = decltype(shared_ptr<T>(shared_ptr<U>{}))>
1350 MultiplexPtr(
const MultiplexPtr<U>& p) : _p(p.
get()) { }
1354 if (_p ==
nullptr) {
1355 throw Error(
"Dereferencing null AnalysisObject pointer. Is there an unbooked histogram variable?");
1360 template <
typename U = T>
1361 auto binning() const ->
std::enable_if_t<hasBinning<typename T::Inner>::value, const typename U::Inner::BinningT&> {
1362 if (_p ==
nullptr) {
1363 throw Error(
"Dereferencing null AnalysisObject pointer. Is there an unbooked histogram variable?");
1365 return _p->binning();
1371 if (_p ==
nullptr) {
1372 throw Error(
"Dereferencing null AnalysisObject pointer. Is there an unbooked histogram variable?");
1379 const typename T::Inner&
operator * ()
const {
return **_p; }
1382 explicit operator bool()
const {
return _p && bool(*_p); }
1389 return _p == other._p;
1394 return _p != other._p;
1399 return _p < other._p;
1404 return _p > other._p;
1409 return _p <= other._p;
1414 return _p >= other._p;
1418 shared_ptr<T>
get()
const {
return _p; }
1437 using MultiplexAOPtr = MultiplexPtr<MultiplexedAO>;
1439 template<
size_t DbnN,
typename... AxisT>
1440 using BinnedDbnPtr = MultiplexPtr<Multiplexer<YODA::BinnedDbn<DbnN, AxisT...>>>;
1442 template<
typename... AxisT>
1443 using BinnedHistoPtr = BinnedDbnPtr<
sizeof...(AxisT), AxisT...>;
1445 template<
typename... AxisT>
1446 using BinnedProfilePtr = BinnedDbnPtr<
sizeof...(AxisT)+1, AxisT...>;
1448 template<
typename... AxisT>
1449 using BinnedEstimatePtr = MultiplexPtr<Multiplexer<YODA::BinnedEstimate<AxisT...>>>;
1452 using ScatterNDPtr = MultiplexPtr<Multiplexer<YODA::ScatterND<N>>>;
1454 using CounterPtr = MultiplexPtr<Multiplexer<YODA::Counter>>;
1455 using Estimate0DPtr = MultiplexPtr<Multiplexer<YODA::Estimate0D>>;
1456 using Histo1DPtr = BinnedHistoPtr<double>;
1457 using Histo2DPtr = BinnedHistoPtr<double,double>;
1458 using Histo3DPtr = BinnedHistoPtr<double,double,double>;
1459 using Profile1DPtr = BinnedProfilePtr<double>;
1460 using Profile2DPtr = BinnedProfilePtr<double,double>;
1461 using Profile3DPtr = BinnedProfilePtr<double,double,double>;
1462 using Estimate1DPtr = BinnedEstimatePtr<double>;
1463 using Estimate2DPtr = BinnedEstimatePtr<double,double>;
1464 using Estimate3DPtr = BinnedEstimatePtr<double,double,double>;
1465 using Scatter1DPtr = ScatterNDPtr<1>;
1466 using Scatter2DPtr = ScatterNDPtr<2>;
1467 using Scatter3DPtr = ScatterNDPtr<3>;
1469 using YODA::Counter;
1470 using YODA::Estimate0D;
1471 using YODA::Histo1D;
1472 using YODA::Histo2D;
1473 using YODA::Histo3D;
1474 using YODA::Profile1D;
1475 using YODA::Profile2D;
1476 using YODA::Profile3D;
1477 using YODA::Estimate1D;
1478 using YODA::Estimate2D;
1479 using YODA::Estimate3D;
1480 using YODA::Scatter1D;
1481 using YODA::Scatter2D;
1482 using YODA::Scatter3D;
1483 using YODA::Point1D;
1484 using YODA::Point2D;
1485 using YODA::Point3D;
1493 inline bool isTmpPath(
const std::string& path,
const bool tmp_only =
false) {
1494 if (tmp_only)
return path.find(
"/TMP/") != string::npos;
1495 return path.find(
"/TMP/") != string::npos || path.find(
"/_") != string::npos;
1500 map<string, YODA::AnalysisObjectPtr>
getRefData(
const string& papername);
1510 template<
typename T>
1511 struct ReferenceTraits { };
1514 struct ReferenceTraits<YODA::Counter> {
1515 using RefT = YODA::Estimate0D;
1519 struct ReferenceTraits<YODA::Estimate0D> {
1520 using RefT = YODA::Estimate0D;
1523 template<
typename... AxisT>
1524 struct ReferenceTraits<YODA::BinnedEstimate<AxisT...>> {
1525 using RefT = YODA::BinnedEstimate<AxisT...>;
1528 template<
size_t DbnN,
typename... AxisT>
1529 struct ReferenceTraits<YODA::BinnedDbn<DbnN, AxisT...>> {
1530 using RefT = YODA::ScatterND<
sizeof...(AxisT)+1>;
1534 struct ReferenceTraits<YODA::ScatterND<N>> {
1535 using RefT = YODA::ScatterND<N>;
1540 template <
typename TPtr>
1563 return a->numPoints() == b->numPoints();
1567 inline bool bookingCompatible(YODA::ScatterNDPtr<N> a, YODA::ScatterNDPtr<N> b) {
1568 return a->numPoints() == b->numPoints();
1571 inline bool beamInfoCompatible(YODA::AnalysisObjectPtr a, YODA::AnalysisObjectPtr b) {
1572 YODA::BinnedEstimatePtr<string> beamsA = std::dynamic_pointer_cast<YODA::BinnedEstimate<string>>(a);
1573 YODA::BinnedEstimatePtr<string> beamsB = std::dynamic_pointer_cast<YODA::BinnedEstimate<string>>(b);
1574 return beamsA && beamsB && (*beamsA == *beamsB) && beamsA->numBins() == 2 &&
1575 fuzzyEquals(beamsA->bin(1).val(), beamsB->bin(1).val()) &&
1576 fuzzyEquals(beamsA->bin(2).val(), beamsB->bin(2).val());
1588 AOPath(
string fullpath)
1589 : _valid(false), _path(fullpath), _raw(false), _tmp(false), _ref(false) {
1590 _valid = init(fullpath);
1594 string path()
const {
return _path; }
1597 string analysis()
const {
return _analysis; }
1600 string analysisWithOptions()
const {
return _analysis + _optionstring; }
1603 string name()
const {
return _name; }
1606 string weight()
const {
return _weight; }
1609 string weightComponent()
const {
1610 if (_weight ==
"")
return _weight;
1611 return "[" + _weight +
"]";
1615 bool isRaw()
const {
return _raw; }
1618 bool isTmp()
const {
return _tmp; }
1621 bool isRef()
const {
return _ref; }
1624 string optionString()
const {
return _optionstring; }
1627 bool hasOptions()
const {
return !_options.empty(); }
1630 void removeOption(
string opt) { _options.erase(opt); fixOptionString(); }
1633 void setOption(
string opt,
string val) { _options[opt] = val; fixOptionString(); }
1636 bool hasOption(
string opt)
const {
return _options.find(opt) != _options.end(); }
1639 string getOption(
string opt)
const {
1640 auto it = _options.find(opt);
1641 if ( it != _options.end() )
return it->second;
1646 void fixOptionString();
1649 string mkPath()
const;
1650 string setPath() {
return _path = mkPath(); }
1656 bool operator<(
const AOPath & other)
const {
1657 return _path < other._path;
1661 bool valid()
const {
return _valid; };
1662 bool operator!()
const {
return !valid(); }
1667 bool init(
string fullpath);
1668 bool chopweight(
string & fullpath);
1669 bool chopoptions(
string & anal);
1674 string _optionstring;
1680 map<string,string> _options;
This is the base class of all analysis classes in Rivet.
Definition Analysis.hh:67
FillCollector(typename YAO::Ptr yao)
Constructor.
Definition RivetYODA.hh:274
const Fills< YAO > & fills() const
Access the fill info subevent stack.
Definition RivetYODA.hh:308
void reset() noexcept
Empty the subevent stack (for start of new event group).
Definition RivetYODA.hh:305
int fill(typename YAO::FillType &&fillCoords, const double weight=1.0, const double fraction=1.0)
Definition RivetYODA.hh:282
void reset() noexcept
Empty the subevent stack (for start of new event group).
Definition RivetYODA.hh:359
FillCollector(typename YAO::Ptr yao)
Constructor.
Definition RivetYODA.hh:328
int fill(typename YAO::FillType &&fillCoords, const double weight=1.0, const double fraction=1.0)
Definition RivetYODA.hh:336
const Fills< YAO > & fills() const
Access the fill info subevent stack.
Definition RivetYODA.hh:362
void reset() noexcept
Empty the subevent stack (for start of new event group).
Definition RivetYODA.hh:467
const Fills< YAO > & fills() const
Access the fill info subevent stack.
Definition RivetYODA.hh:470
int fill(typename YAO::FillType &&fillCoords, const double weight=1.0, const double fraction=1.0)
Definition RivetYODA.hh:444
FillCollector(typename YAO::Ptr yao)
Constructor.
Definition RivetYODA.hh:436
void reset() noexcept
Empty the subevent stack (for start of new event group).
Definition RivetYODA.hh:413
FillCollector(typename YAO::Ptr yao)
Constructor.
Definition RivetYODA.hh:382
const Fills< YAO > & fills() const
Access the fill info subevent stack.
Definition RivetYODA.hh:416
int fill(typename YAO::FillType &&fillCoords, const double weight=1.0, const double fraction=1.0)
Definition RivetYODA.hh:390
const Fills< YAO > & fills() const
Access the fill info subevent stack.
Definition RivetYODA.hh:524
FillCollector(typename YAO::Ptr yao)
Constructor.
Definition RivetYODA.hh:490
void reset() noexcept
Empty the subevent stack (for start of new event group).
Definition RivetYODA.hh:521
int fill(typename YAO::FillType &&fillCoords, const double weight=1.0, const double fraction=1.0)
Definition RivetYODA.hh:498
int fill(typename YAO::FillType &&fillCoords, const double weight=1.0, const double fraction=1.0)
Definition RivetYODA.hh:552
const Fills< YAO > & fills() const
Access the fill info subevent stack.
Definition RivetYODA.hh:578
FillCollector(typename YAO::Ptr yao)
Constructor.
Definition RivetYODA.hh:544
void reset() noexcept
Empty the subevent stack (for start of new event group).
Definition RivetYODA.hh:575
void reset() noexcept
Empty the subevent stack (for start of new event group).
Definition RivetYODA.hh:251
FillCollector(typename YAO::Ptr yao)
Definition RivetYODA.hh:224
const Fills< YAO > & fills() const
Access the fill info subevent stack.
Definition RivetYODA.hh:254
int fill(typename YAO::FillType &&fillCoords, const double weight=1.0, const double fraction=1.0)
Definition RivetYODA.hh:232
FillCollector(typename YAO::Ptr yao)
Definition RivetYODA.hh:628
int fill(const double weight=1.0, const double fraction=1.0)
Definition RivetYODA.hh:185
void reset()
Empty the subevent stack (for start of new event group).
Definition RivetYODA.hh:192
const Fills< YAO > & fills() const
Access the fill info subevent stack.
Definition RivetYODA.hh:195
FillCollector(typename YAO::Ptr yao)
Definition RivetYODA.hh:179
FillCollector(typename YAO::Ptr yao)
Definition RivetYODA.hh:605
FillCollector(typename YAO::Ptr yao)
Definition RivetYODA.hh:650
FillCollectors which are used to temporarily cache unaggregated fills until collapsed by the Multiple...
Definition RivetYODA.hh:158
Definition RivetYODA.hh:1330
T & operator->()
Goes right through to the active Multiplexer<YODA> object's members.
Definition RivetYODA.hh:1353
bool operator>(const MultiplexPtr other) const
Greater-than for ptr ordering.
Definition RivetYODA.hh:1403
shared_ptr< T > get() const
Get the internal shared ptr.
Definition RivetYODA.hh:1418
T::Inner & operator*()
The active YODA object.
Definition RivetYODA.hh:1378
bool operator<(const MultiplexPtr &other) const
Less-than for ptr ordering.
Definition RivetYODA.hh:1398
MultiplexPtr(const vector< string > &weightNames, const typename T::Inner &p)
Convenience constructor, pass through to the Multiplexer constructor.
Definition RivetYODA.hh:1341
bool operator!() const
Object invalidity check.
Definition RivetYODA.hh:1385
bool operator<=(const MultiplexPtr &other) const
Less-equals for ptr ordering.
Definition RivetYODA.hh:1408
bool operator!=(const MultiplexPtr &other) const
Object invalidity check.
Definition RivetYODA.hh:1393
bool operator==(const MultiplexPtr &other) const
Object validity check.
Definition RivetYODA.hh:1388
bool operator>=(const MultiplexPtr &other) const
Greater-equals for ptr ordering.
Definition RivetYODA.hh:1413
Multiplexer base class.
Definition RivetYODA.hh:931
virtual string basePath() const =0
The histogram path, without a variation suffix.
virtual YODA::AnalysisObject * operator->()=0
Access the active analysis object for function calls.
virtual void unsetActiveWeight()=0
Unset the active-object pointer.
bool operator==(const MultiplexedAO &p)
Test for equality.
Definition RivetYODA.hh:976
virtual void pushToFinal()=0
Sync the persistent histograms to the final collection.
virtual YODA::AnalysisObjectPtr activeAO() const =0
A shared pointer to the active YODA AO.
virtual void newSubEvent()=0
Add a new layer of subevent fill staging.
virtual void collapseEventGroup(const vector< std::valarray< double > > &weight, const double nlowfrac=0.0)=0
Sync the fill proxies to the persistent histogram.
YODA::AnalysisObject Inner
The type being represented is a generic AO.
Definition RivetYODA.hh:937
bool operator!=(const MultiplexedAO &p)
Test for inequality.
Definition RivetYODA.hh:979
virtual const YODA::AnalysisObject & operator*() const =0
Access the active analysis object as a reference.
virtual void initBootstrap()=0
Set the size of the bootstrap vectors.
virtual void setActiveWeightIdx(size_t iWeight)=0
Set active object for analyze.
virtual void setActiveFinalWeightIdx(size_t iWeight)=0
Set active object for finalize.
Type-specific multiplexed YODA analysis object.
Definition RivetYODA.hh:1011
void initBootstrap()
Helper method to resize aux vectors to AO size.
Definition RivetYODA.hh:1277
friend bool operator<(const Multiplexer a, const Multiplexer &b)
Less-than operator.
Definition RivetYODA.hh:1114
string baseName() const
Get the AO name of the object, without variation suffix.
Definition RivetYODA.hh:1067
T & operator*()
Forwarding dereference operator.
Definition RivetYODA.hh:1088
void newSubEvent()
Create a new FillCollector for the next sub-event.
Definition RivetYODA.hh:1154
void pushToFinal()
Definition RivetYODA.hh:1245
void reset()
Clear the active object pointer.
Definition RivetYODA.hh:1146
T::Ptr active() const
Definition RivetYODA.hh:1046
const vector< typename T::Ptr > & final() const
Definition RivetYODA.hh:1266
bool operator!() const
Definition RivetYODA.hh:1078
string basePath() const
Get the AO path of the object, without variation suffix.
Definition RivetYODA.hh:1064
void collapseEventGroup(const vector< std::valarray< double > > &weights, const double nlowfrac=0.0)
Pushes the (possibly collapsed) fill(s) into the persistent objects.
Definition RivetYODA.hh:1161
friend bool operator!=(const Multiplexer &a, const Multiplexer &b)
Inequality operator.
Definition RivetYODA.hh:1109
friend bool operator==(const Multiplexer &a, const Multiplexer &b)
Equality operator.
Definition RivetYODA.hh:1095
void setActiveFinalWeightIdx(size_t iWeight)
Set the active-object pointer to point at a variation in the final set.
Definition RivetYODA.hh:1138
void unsetActiveWeight()
Unset the active-object pointer.
Definition RivetYODA.hh:1143
void collapseSubevents(const vector< std::valarray< double > > &weights, const double nlowfrac)
Definition RivetYODA.hh:1203
const vector< typename T::Ptr > & persistent() const
Definition RivetYODA.hh:1257
T::Ptr persistent(const size_t iWeight)
Direct access to the persistent object in weight stream iWeight.
Definition RivetYODA.hh:1262
void setActiveWeightIdx(size_t iWeight)
Set the active-object pointer to point at a variation in the persistent set.
Definition RivetYODA.hh:1133
T * operator->()
Forwarding dereference-call operator.
Definition RivetYODA.hh:1082
YODA::AnalysisObjectPtr activeAO() const
Get the currently active analysis object.
Definition RivetYODA.hh:1274
T Inner
Typedef for the YODA type being represented.
Definition RivetYODA.hh:1019
double Weight
Typedef for weights.
Definition RivetYODA.hh:135
vector< Fill< T > > Fills
A collection of several Fill objects.
Definition RivetYODA.hh:143
pair< typename T::FillType, Weight > Fill
A single fill is a (FillType, Weight) pair.
Definition RivetYODA.hh:139
map< string, YODA::AnalysisObjectPtr > getRefData(const string &papername)
bool bookingCompatible(TPtr a, TPtr b)
Definition RivetYODA.hh:1541
string getDatafilePath(const string &papername)
Get the file system path to the reference file for this paper.
Definition MC_CENT_PPB_Projections.hh:10
Cut operator!(const Cut &cptr)
Logical NOT operation on a cut.
std::enable_if< std::is_arithmetic< N1 >::value &&std::is_arithmetic< N2 >::value, typenamestd::common_type< N1, N2 >::type >::type max(N1 a, N2 b)
Get the maximum of two numbers.
Definition MathUtils.hh:111
bool copyAO(YODA::AnalysisObjectPtr src, YODA::AnalysisObjectPtr dst, const double scale=1.0)
Definition RivetYODA.hh:64
std::enable_if< std::is_arithmetic< N1 >::value &&std::is_arithmetic< N2 >::value, typenamestd::common_type< N1, N2 >::type >::type min(N1 a, N2 b)
Get the minimum of two numbers.
Definition MathUtils.hh:102
std::enable_if< std::is_arithmetic< N1 >::value &&std::is_arithmetic< N2 >::value &&(std::is_floating_point< N1 >::value||std::is_floating_point< N2 >::value), bool >::type fuzzyEquals(N1 a, N2 b, double tolerance=1e-5)
Compare two numbers for equality with a degree of fuzziness.
Definition MathUtils.hh:57
std::enable_if< std::is_arithmetic< NUM >::value, NUM >::type sqr(NUM a)
Named number-type squaring operation.
Definition MathUtils.hh:219
Generic runtime Rivet error.
Definition Exceptions.hh:12
A polymorphic base type for the AO type handles.
Definition RivetYODA.hh:79
The type-specific handle that can perform type-specific operations for objects of type T.
Definition RivetYODA.hh:100