@@ -769,9 +769,12 @@ struct HistoSys {
769769 RooAbsReal const *param = nullptr ;
770770 std::vector<double > low;
771771 std::vector<double > high;
772+ // Used to validate duplicate RooFit modifiers. This is intentionally not serialized until HS3 defines the
773+ // structured histosys interpolation representation.
774+ int interpolationCode = 4 ;
772775 RooAbsPdf const *constraint = nullptr ;
773- HistoSys (const std::string &n, RooAbsReal *const p, RooHistFunc *l, RooHistFunc *h, const RooAbsPdf *c)
774- : name(n), param(p), constraint(c)
776+ HistoSys (const std::string &n, RooAbsReal *const p, RooHistFunc *l, RooHistFunc *h, int i, const RooAbsPdf *c)
777+ : name(n), param(p), interpolationCode(i), constraint(c)
775778 {
776779 low.assign (l->dataHist ().weightArray (), l->dataHist ().weightArray () + l->dataHist ().numEntries ());
777780 high.assign (h->dataHist ().weightArray (), h->dataHist ().weightArray () + h->dataHist ().numEntries ());
@@ -904,7 +907,7 @@ NormSys parseOverallModifierFormula(const std::string &s, RooFormulaVar *formula
904907 return sys;
905908}
906909
907- void collectElements (RooArgSet &elems, RooAbsArg *arg)
910+ void collectElements (RooArgList &elems, RooAbsArg *arg)
908911{
909912 if (auto prod = dynamic_cast <RooProduct *>(arg)) {
910913 for (const auto &e : prod->components ()) {
@@ -1085,7 +1088,7 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons
10851088 }
10861089 };
10871090
1088- RooArgSet elems;
1091+ RooArgList elems;
10891092 collectElements (elems, func);
10901093 collectElements (elems, sumpdf->coefList ().at (sampleidx));
10911094 processElements (elems, processElements);
@@ -1113,7 +1116,7 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons
11131116 if (!constraint && !var->isConstant ()) {
11141117 RooJSONFactoryWSTool::error (" cannot find constraint for " + std::string (var->GetName ()));
11151118 } else {
1116- sample.histosys .emplace_back (sysname, var, lo, hi, constraint);
1119+ sample.histosys .emplace_back (sysname, var, lo, hi, pip-> interpolationCodes ()[i], constraint);
11171120 }
11181121 }
11191122 }
@@ -1194,6 +1197,154 @@ Channel readChannel(RooJSONFactoryWSTool *tool, const std::string &pdfname, cons
11941197 return channel;
11951198}
11961199
1200+ bool hasSameMetadata (const RooAbsArg *lhs, const RooAbsArg *rhs)
1201+ {
1202+ if (!lhs || !rhs) {
1203+ return lhs == rhs;
1204+ }
1205+ return std::string{lhs->GetName ()} == rhs->GetName () && lhs->IsA () == rhs->IsA ();
1206+ }
1207+
1208+ [[noreturn]] void duplicateModifierError (const Channel &channel, const Sample &sample, std::string_view type,
1209+ std::string_view name, std::string_view reason)
1210+ {
1211+ std::stringstream ss;
1212+ ss << " cannot combine duplicate modifier '" << name << " ' of type '" << type << " ' in sample '" << sample.name
1213+ << " ' of channel '" << channel.name << " ': " << reason;
1214+ RooJSONFactoryWSTool::error (ss.str ().c_str ());
1215+ }
1216+
1217+ void warnDuplicateModifiersCombined (const Channel &channel, const Sample &sample, std::string_view type,
1218+ std::string_view name, std::size_t count)
1219+ {
1220+ std::stringstream ss;
1221+ ss << " combined " << count << " duplicate modifiers named '" << name << " ' of type '" << type << " ' in sample '"
1222+ << sample.name << " ' of channel '" << channel.name << " '" ;
1223+ RooJSONFactoryWSTool::warning (ss.str ());
1224+ }
1225+
1226+ // Multiplicatively combining two normsys is only faithful when the interpolation is done in log-space, so that
1227+ // f1(alpha) * f2(alpha) is again representable by a single normsys with the multiplied lo/hi factors. This holds for
1228+ // the piecewise-exponential code 1 (exact everywhere) and for the default code 4 (exact at the +-1 sigma anchors and in
1229+ // the exponential extrapolation region). The linear-space codes (e.g. 0 and 2) would turn the product into a shape that
1230+ // cannot be represented by a single normsys, so those must not be merged.
1231+ bool normSysSupportsMultiplicativeMerge (int interpolationCode)
1232+ {
1233+ return interpolationCode == 1 || interpolationCode == 4 ;
1234+ }
1235+
1236+ // Combines runs of adjacent modifiers that share the same name (the container is sorted by name beforehand) into a
1237+ // single modifier. The shared metadata (constraint, parameter and interpolation code) must be identical across the
1238+ // duplicates; the type-specific `combine` callable performs the actual merge and any additional validation.
1239+ template <class Modifiers , class CombineFn >
1240+ void mergeDuplicateModifiers (const Channel &channel, const Sample &sample, Modifiers &modifiers,
1241+ std::string_view type, CombineFn combine)
1242+ {
1243+ Modifiers mergedModifiers;
1244+ mergedModifiers.reserve (modifiers.size ());
1245+
1246+ for (std::size_t begin = 0 ; begin < modifiers.size ();) {
1247+ std::size_t end = begin + 1 ;
1248+ while (end < modifiers.size () && modifiers[end].name == modifiers[begin].name ) {
1249+ ++end;
1250+ }
1251+
1252+ auto merged = modifiers[begin];
1253+ for (std::size_t i = begin + 1 ; i < end; ++i) {
1254+ const auto &modifier = modifiers[i];
1255+ if (!hasSameMetadata (merged.constraint , modifier.constraint )) {
1256+ duplicateModifierError (channel, sample, type, merged.name , " constraint metadata differs" );
1257+ }
1258+ if (!hasSameMetadata (merged.param , modifier.param )) {
1259+ duplicateModifierError (channel, sample, type, merged.name , " parameter metadata differs" );
1260+ }
1261+ if (merged.interpolationCode != modifier.interpolationCode ) {
1262+ duplicateModifierError (channel, sample, type, merged.name , " interpolation codes differ" );
1263+ }
1264+ combine (merged, modifier);
1265+ }
1266+
1267+ if (end - begin > 1 ) {
1268+ warnDuplicateModifiersCombined (channel, sample, type, merged.name , end - begin);
1269+ }
1270+ mergedModifiers.emplace_back (std::move (merged));
1271+ begin = end;
1272+ }
1273+
1274+ modifiers = std::move (mergedModifiers);
1275+ }
1276+
1277+ void mergeDuplicateNormSys (const Channel &channel, Sample &sample)
1278+ {
1279+ mergeDuplicateModifiers (channel, sample, sample.normsys , " normsys" ,
1280+ [&](NormSys &merged, const NormSys &modifier) {
1281+ if (!normSysSupportsMultiplicativeMerge (merged.interpolationCode )) {
1282+ duplicateModifierError (
1283+ channel, sample, " normsys" , merged.name ,
1284+ " multiplicative combination is only valid for log-space interpolation codes" );
1285+ }
1286+ merged.low *= modifier.low ;
1287+ merged.high *= modifier.high ;
1288+ });
1289+ }
1290+
1291+ void mergeDuplicateHistoSys (const Channel &channel, Sample &sample)
1292+ {
1293+ const std::size_t nBins = sample.hist .size ();
1294+ mergeDuplicateModifiers (channel, sample, sample.histosys , " histosys" ,
1295+ [&](HistoSys &merged, const HistoSys &modifier) {
1296+ if (merged.interpolationCode != 4 ) {
1297+ duplicateModifierError (
1298+ channel, sample, " histosys" , merged.name ,
1299+ " non-default interpolation cannot currently be represented by the HS3 exporter" );
1300+ }
1301+ if (merged.low .size () != nBins || merged.high .size () != nBins ||
1302+ modifier.low .size () != nBins || modifier.high .size () != nBins) {
1303+ duplicateModifierError (channel, sample, " histosys" , merged.name ,
1304+ " histogram binning differs" );
1305+ }
1306+ for (std::size_t bin = 0 ; bin < nBins; ++bin) {
1307+ merged.low [bin] += modifier.low [bin] - sample.hist [bin];
1308+ merged.high [bin] += modifier.high [bin] - sample.hist [bin];
1309+ }
1310+ });
1311+ }
1312+
1313+ void ensureUniqueModifiers (const Channel &channel, const Sample &sample)
1314+ {
1315+ std::set<std::pair<std::string, std::string>> seen;
1316+ auto add = [&](std::string type, const std::string &name) {
1317+ if (!seen.emplace (type, name).second ) {
1318+ duplicateModifierError (channel, sample, type, name,
1319+ " this modifier type cannot be combined without changing its meaning" );
1320+ }
1321+ };
1322+
1323+ for (const auto &modifier : sample.normfactors )
1324+ add (" normfactor" , modifier.name );
1325+ for (const auto &modifier : sample.normsys )
1326+ add (" normsys" , modifier.name );
1327+ for (const auto &modifier : sample.histosys )
1328+ add (" histosys" , modifier.name );
1329+ for (const auto &modifier : sample.shapesys )
1330+ add (" shapesys" , modifier.name );
1331+ for (const auto &modifier : sample.otherElements )
1332+ add (" custom" , modifier.name );
1333+ for (const auto &modifier : sample.tmpElements )
1334+ add (" custom" , modifier.name );
1335+ if (sample.useBarlowBeestonLight )
1336+ add (::Literals::staterror, ::Literals::staterror);
1337+ }
1338+
1339+ void canonicalizeModifiers (Channel &channel)
1340+ {
1341+ for (auto &sample : channel.samples ) {
1342+ mergeDuplicateNormSys (channel, sample);
1343+ mergeDuplicateHistoSys (channel, sample);
1344+ ensureUniqueModifiers (channel, sample);
1345+ }
1346+ }
1347+
11971348void configureStatError (Channel &channel)
11981349{
11991350 for (auto &sample : channel.samples ) {
@@ -1424,6 +1575,8 @@ bool tryExportHistFactory(RooJSONFactoryWSTool *tool, const std::string &pdfname
14241575 }
14251576 }
14261577
1578+ canonicalizeModifiers (channel);
1579+
14271580 // stat error handling
14281581 configureStatError (channel);
14291582
0 commit comments