85#define NaN std::numeric_limits<double>::quiet_NaN()
94template <
typename Test,
template <
typename...>
class Ref>
98template <
template <
typename...>
class Ref,
typename... Args>
109template <
class MatrixT>
110inline size_t size(
const MatrixT &matrix);
120template <
class MatrixT>
123 if (!stream.good()) {
126 for (
size_t i = 0; i <
size(matrix); ++i) {
127 for (
size_t j = 0; j <
size(matrix); ++j) {
129 stream << std::setprecision(RooFit::SuperFloatPrecision::digits10) << matrix(i, j) <<
"\t";
131 stream << matrix(i, j) <<
"\t";
141template <
class MatrixT>
144 std::ofstream of(fname);
146 cerr <<
"unable to read file '" << fname <<
"'!" << std::endl;
155#pragma GCC diagnostic push
156#pragma GCC diagnostic ignored "-Wshadow"
157#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
158#include <boost/numeric/ublas/io.hpp>
159#include <boost/numeric/ublas/lu.hpp>
160#include <boost/numeric/ublas/matrix.hpp>
161#include <boost/numeric/ublas/matrix_expression.hpp>
162#include <boost/numeric/ublas/symmetric.hpp>
163#include <boost/numeric/ublas/triangular.hpp>
164#include <boost/operators.hpp>
166#pragma GCC diagnostic pop
168typedef boost::numeric::ublas::matrix<RooFit::SuperFloat>
Matrix;
175 for (
size_t i = 0; i < mat.size1(); ++i) {
176 for (
size_t j = 0; j < mat.size2(); ++j) {
177 std::cout << std::setprecision(RooFit::SuperFloatPrecision::digits10) << mat(i, j) <<
" ,\t";
179 std::cout << std::endl;
187inline size_t size<Matrix>(
const Matrix &matrix)
189 return matrix.size1();
197 return boost::numeric::ublas::identity_matrix<RooFit::SuperFloat>(
n);
207 for (
size_t i = 0; i <
n; ++i) {
208 for (
size_t j = 0; j <
n; ++j) {
209 mat(i, j) =
double(in(i, j));
222 for (
size_t i = 0; i <
n; ++i) {
223 for (
size_t j = 0; j <
n; ++j) {
224 mat(i, j) =
double(in(i, j));
236 return prod(
m, otherM);
244 boost::numeric::ublas::permutation_matrix<size_t> pm(
size(matrix));
248 int res = lu_factorize(lu, pm);
250 std::stringstream ss;
255 lu_substitute(lu, pm, inverse);
256 }
catch (boost::numeric::ublas::internal_logic &error) {
310 bool status = lu.
Invert(inverse);
313 std::cerr <<
" matrix is not invertible!" << std::endl;
316 const size_t n =
size(inverse);
318 for (
size_t i = 0; i <
n; ++i)
319 for (
size_t j = 0; j <
n; ++j)
337typedef std::vector<std::vector<bool>> FeynmanDiagram;
338typedef std::vector<std::vector<int>> MorphFuncPattern;
339typedef std::map<int, std::unique_ptr<RooAbsReal>> FormulaList;
347 retval.ReplaceAll(
"/",
"_");
348 retval.ReplaceAll(
"^",
"");
349 retval.ReplaceAll(
"*",
"X");
350 retval.ReplaceAll(
"[",
"");
351 retval.ReplaceAll(
"]",
"");
359std::string concatNames(
const List &
c,
const char *
sep)
361 std::stringstream ss;
366 ss << itr->GetName();
376template <
class A,
class B>
377inline void assignElement(
A &
a,
const B &
b)
379 a =
static_cast<A>(
b);
384template <
class MatrixT>
385inline MatrixT readMatrixFromStreamT(std::istream &stream)
387 std::vector<std::vector<RooFit::SuperFloat>> matrix;
388 std::vector<RooFit::SuperFloat>
line;
389 while (!stream.eof()) {
390 if (stream.peek() ==
'\n') {
398 while (stream.peek() ==
' ' || stream.peek() ==
'\t') {
401 if (stream.peek() ==
'\n') {
402 matrix.push_back(
line);
406 MatrixT retval(matrix.size(), matrix.size());
407 for (
size_t i = 0; i < matrix.size(); ++i) {
408 if (matrix[i].
size() != matrix.size()) {
409 std::cerr <<
"matrix read from stream doesn't seem to be square!" << std::endl;
411 for (
size_t j = 0; j < matrix[i].size(); ++j) {
412 assignElement(retval(i, j), matrix[i][j]);
421template <
class MatrixT>
422inline MatrixT readMatrixFromFileT(
const char *fname)
424 std::ifstream in(fname);
426 std::cerr <<
"unable to read file '" << fname <<
"'!" << std::endl;
428 MatrixT mat = readMatrixFromStreamT<MatrixT>(in);
437void readValues(std::map<const std::string, T> &myMap,
TH1 *h_pc)
441 for (
int ibx = 1; ibx <= h_pc->
GetNbinsX(); ++ibx) {
446 if (!s_coup.empty()) {
447 myMap[s_coup] =
T(coup_val);
456void setOwnerRecursive(
TFolder *theFolder)
461 for (
auto *subdir : *subdirs) {
462 auto thisfolder =
dynamic_cast<TFolder *
>(subdir);
465 setOwnerRecursive(thisfolder);
478std::unique_ptr<TFolder> readOwningFolderFromFile(
TDirectory *inFile,
const std::string &folderName)
480 std::unique_ptr<TFolder> theFolder(inFile->
Get<
TFolder>(folderName.c_str()));
482 std::cerr <<
"Error: unable to access data from folder '" << folderName <<
"' from file '" << inFile->
GetName()
483 <<
"'!" << std::endl;
486 setOwnerRecursive(theFolder.get());
499template <
class AObjType>
500std::unique_ptr<AObjType> loadFromFileResidentFolder(
TDirectory *inFile,
const std::string &folderName,
501 const std::string &objName,
bool notFoundError =
true)
503 auto folder = readOwningFolderFromFile(inFile, folderName);
507 AObjType *loadedObject =
dynamic_cast<AObjType *
>(folder->FindObject(objName.c_str()));
510 std::stringstream errstr;
511 errstr <<
"Error: unable to retrieve object '" << objName <<
"' from folder '" << folderName
512 <<
"'. contents are:";
513 TIter next(folder->GetListOfFolders()->begin());
516 errstr <<
" " <<
f->GetName();
518 std::cerr << errstr.str() << std::endl;
524 return std::unique_ptr<AObjType>{
static_cast<AObjType *
>(loadedObject->Clone())};
531void readValues(std::map<const std::string, T> &myMap,
TDirectory *
file,
const std::string &
name,
532 const std::string &key =
"param_card",
bool notFoundError =
true)
534 auto h_pc = loadFromFileResidentFolder<TH1F>(
file,
name, key, notFoundError);
535 readValues(myMap, h_pc.get());
544void readValues(std::map<
const std::string, std::map<const std::string, T>> &inputParameters,
TDirectory *
f,
545 const std::vector<std::string> &names,
const std::string &key =
"param_card",
bool notFoundError =
true)
547 inputParameters.clear();
550 for (
size_t i = 0; i < names.size(); i++) {
551 const std::string
name(names[i]);
553 readValues(inputParameters[
name],
f,
name, key, notFoundError);
571 std::cerr <<
"could not open file '" <<
filename <<
"'!" << std::endl;
593inline void extractServers(
const RooAbsArg &coupling,
T2 &operators)
596 for (
const auto server : coupling.
servers()) {
597 extractServers(*server, operators);
601 operators.add(coupling);
609inline void extractOperators(
const T1 &couplings,
T2 &operators)
613 for (
auto itr : couplings) {
614 extractServers(*itr, operators);
622inline void extractOperators(
const T1 &
vec,
T2 &operators)
624 for (
const auto &
v :
vec) {
625 extractOperators(
v, operators);
632template <
class T1,
class T2>
633inline void extractCouplings(
const T1 &inCouplings,
T2 &outCouplings)
635 for (
auto itr : inCouplings) {
636 if (!outCouplings.find(itr->GetName())) {
639 outCouplings.add(*itr);
648inline bool setParam(
RooRealVar *
p,
double val,
bool force)
651 if (val >
p->getMax()) {
655 std::cerr <<
": parameter " <<
p->
GetName() <<
" out of bounds: " << val <<
" > " <<
p->getMax() << std::endl;
658 }
else if (val < p->getMin()) {
662 std::cerr <<
": parameter " <<
p->
GetName() <<
" out of bounds: " << val <<
" < " <<
p->getMin() << std::endl;
675template <
class T1,
class T2>
676inline bool setParams(
const T2 &args,
T1 val)
678 for (
auto itr : args) {
682 setParam(param, val,
true);
691template <
class T1,
class T2>
693setParams(
const std::map<const std::string, T1> &point,
const T2 &args,
bool force =
false,
T1 defaultVal = 0)
696 for (
auto itr : args) {
700 ok = setParam(param, defaultVal, force) && ok;
703 for (
auto paramit : point) {
705 const std::string param(paramit.first);
711 ok = setParam(
p, paramit.second, force) && ok;
721inline bool setParams(
TH1 *hist,
const T &args,
bool force =
false)
725 for (
auto itr : args) {
729 ok = setParam(param, 0., force) && ok;
734 for (
int i = 1; i <= ax->
GetNbins(); ++i) {
752 for (
auto itr : parameters) {
769 bool binningOK =
false;
770 for (
auto sampleit : inputParameters) {
771 const std::string sample(sampleit.first);
772 auto hist = loadFromFileResidentFolder<TH1>(
file, sample, varname,
true);
780 auto it = list_hf.find(sample);
781 if (it != list_hf.end()) {
793 std::vector<double> bins;
794 for (
int i = 1; i <
n + 1; ++i) {
802 TString histname = makeValidName(
Form(
"dh_%s_%s", sample.c_str(),
name));
803 TString funcname = makeValidName(
Form(
"phys_%s_%s", sample.c_str(),
name));
807 auto dh = std::make_unique<RooDataHist>(histname, histname, vars, hist.get());
809 auto hf = std::make_unique<RooHistFunc>(funcname, funcname, var, std::move(dh));
811 list_hf[sample] = idx;
822void collectRooAbsReal(
const char * ,
TDirectory *
file, std::map<std::string, int> &list_hf,
823 RooArgList &physics,
const std::string &varname,
826 for (
auto sampleit : inputParameters) {
827 const std::string sample(sampleit.first);
828 auto obj = loadFromFileResidentFolder<RooAbsReal>(
file, sample, varname,
true);
831 auto it = list_hf.find(sample);
832 if (it == list_hf.end()) {
834 list_hf[sample] = idx;
848 for (
auto sampleit : inputParameters) {
849 const std::string sample(sampleit.first);
850 auto obj = loadFromFileResidentFolder<TObject>(
file, sample, varname,
false);
857 TPair *pair =
dynamic_cast<TPair *
>(obj.get());
863 std::stringstream errstr;
864 errstr <<
"Error: unable to retrieve cross section '" << varname <<
"' from folder '" << sample;
868 auto it = list_xs.find(sample.c_str());
870 if (it != list_xs.end()) {
874 std::string objname =
Form(
"phys_%s_%s",
name, sample.c_str());
875 auto xsOwner = std::make_unique<RooRealVar>(objname.c_str(), objname.c_str(), xsection->
GetVal());
879 list_xs[sample] = idx;
880 physics.
addOwned(std::move(xsOwner));
881 assert(physics.
at(idx) == xs);
893void collectCrosssectionsTPair(
const char *
name,
TDirectory *
file, std::map<std::string, int> &list_xs,
894 RooArgList &physics,
const std::string &varname,
const std::string &basefolder,
897 auto pair = loadFromFileResidentFolder<TPair>(
file, basefolder, varname,
false);
901 collectCrosssections<double>(
name,
file, list_xs, physics, varname, inputParameters);
903 collectCrosssections<float>(
name,
file, list_xs, physics, varname, inputParameters);
905 std::cerr <<
"cannot morph objects of class 'TPair' if parameter is not "
919void collectPolynomialsHelper(
const FeynmanDiagram &diagram, MorphFuncPattern &morphfunc, std::vector<int> &term,
920 int vertexid,
bool first)
923 for (
size_t i = 0; i < diagram[vertexid - 1].size(); ++i) {
924 if (!diagram[vertexid - 1][i])
926 std::vector<int> newterm(term);
929 ::collectPolynomialsHelper(diagram, morphfunc, newterm, vertexid,
false);
931 ::collectPolynomialsHelper(diagram, morphfunc, newterm, vertexid - 1,
true);
936 for (
size_t i = 0; i < morphfunc.size(); ++i) {
937 bool thisfound =
true;
938 for (
size_t j = 0; j < morphfunc[i].size(); ++j) {
939 if (morphfunc[i][j] != term[j]) {
950 morphfunc.push_back(term);
958void collectPolynomials(MorphFuncPattern &morphfunc,
const FeynmanDiagram &diagram)
960 int nvtx(diagram.size());
961 std::vector<int> term(diagram[0].
size(), 0);
963 ::collectPolynomialsHelper(diagram, morphfunc, term, nvtx,
true);
970inline void fillFeynmanDiagram(FeynmanDiagram &diagram,
const std::vector<List *> &vertices,
RooArgList &couplings)
972 const int ncouplings = couplings.
getSize();
974 for (
auto const &
vertex : vertices) {
975 std::vector<bool> vertexCouplings(ncouplings,
false);
978 for (
auto citr : couplings) {
982 std::cerr <<
"encountered invalid list of couplings in vertex!" << std::endl;
986 vertexCouplings[idx] =
true;
989 diagram.push_back(vertexCouplings);
996template <
class MatrixT,
class T1,
class T2>
1000 const size_t dim = inputParameters.size();
1001 MatrixT matrix(dim, dim);
1003 for (
auto sampleit : inputParameters) {
1004 const std::string sample(sampleit.first);
1006 if (!setParams<double>(sampleit.second, args,
true, 0)) {
1007 std::cout <<
"unable to set parameters for sample " << sample <<
"!" << std::endl;
1009 auto flagit = flagValues.find(sample);
1010 if (flagit != flagValues.end() && !setParams<int>(flagit->second, flags,
true, 1)) {
1011 std::cout <<
"unable to set parameters for sample " << sample <<
"!" << std::endl;
1015 for (
auto const &formula : formulas) {
1016 if (!formula.second) {
1017 std::cerr <<
"Error: invalid formula encountered!" << std::endl;
1019 matrix(row, col) = formula.second->getVal();
1032 if (inputParameters.size() != formulas.size()) {
1033 std::stringstream ss;
1034 ss <<
"matrix is not square, consistency check failed: " << inputParameters.size() <<
" samples, "
1035 << formulas.size() <<
" expressions:" << std::endl;
1036 ss <<
"formulas: " << std::endl;
1037 for (
auto const &formula : formulas) {
1038 ss << formula.second->GetTitle() << std::endl;
1040 ss <<
"samples: " << std::endl;
1041 for (
auto sample : inputParameters) {
1042 ss << sample.first << std::endl;
1044 std::cerr << ss.str() << std::endl;
1051inline void inverseSanity(
const Matrix &matrix,
const Matrix &inverse,
double &unityDeviation,
double &largestWeight)
1053 Matrix unity(inverse * matrix);
1055 unityDeviation = 0.;
1057 const size_t dim =
size(unity);
1058 for (
size_t i = 0; i < dim; ++i) {
1059 for (
size_t j = 0; j < dim; ++j) {
1060 if (inverse(i, j) > largestWeight) {
1061 largestWeight = (
double)inverse(i, j);
1063 if (
std::abs(unity(i, j) -
static_cast<int>(i == j)) > unityDeviation) {
1064 unityDeviation =
std::abs((
double)unity(i, j)) -
static_cast<int>(i == j);
1072template <
class List>
1075 for (
auto sampleit : inputParameters) {
1076 const std::string sample(sampleit.first);
1077 RooAbsArg *arg = args.find(sample.c_str());
1079 std::cerr <<
"detected name conflict: cannot use sample '" << sample
1080 <<
"' - a parameter with the same name of type '" << arg->
ClassName() <<
"' is present in set '"
1081 << args.GetName() <<
"'!" << std::endl;
1093 const std::vector<std::vector<std::string>> &nonInterfering)
1102 const int ncouplings = couplings.
getSize();
1103 std::vector<bool> couplingsZero(ncouplings,
true);
1104 std::map<TString, bool> flagsZero;
1107 extractOperators(couplings, operators);
1108 size_t nOps = operators.
getSize();
1110 for (
auto sampleit : inputParameters) {
1111 const std::string sample(sampleit.first);
1112 if (!setParams(sampleit.second, operators,
true)) {
1113 std::cerr <<
"unable to set parameters for sample '" << sample <<
"'!" << std::endl;
1116 if ((
int)nOps != (operators.
getSize())) {
1117 std::cerr <<
"internal error, number of operators inconsistent!" << std::endl;
1123 for (
auto itr1 : couplings) {
1125 if (obj0->
getVal() != 0) {
1126 couplingsZero[idx] =
false;
1132 for (
auto itr2 : flags) {
1133 auto obj1 =
dynamic_cast<RooAbsReal *
>(itr2);
1136 for (
auto sampleit : inputFlags) {
1137 const auto &flag = sampleit.second.find(obj1->GetName());
1138 if (flag != sampleit.second.end()) {
1139 if (flag->second == 0.)
1145 if (nZero > 0 && nNonZero == 0)
1146 flagsZero[obj1->GetName()] =
true;
1148 flagsZero[obj1->GetName()] =
false;
1151 FormulaList formulas;
1152 for (
size_t i = 0; i < morphfunc.size(); ++i) {
1154 bool isZero =
false;
1157 for (
const auto &
group : nonInterfering) {
1158 int nInterferingOperators = 0;
1159 for (
size_t j = 0; j < morphfunc[i].size(); ++j) {
1160 if (morphfunc[i][j] % 2 == 0)
1163 if (std::find(
group.begin(),
group.end(), couplings.at(j)->GetName()) !=
group.end()) {
1164 nInterferingOperators++;
1167 if (nInterferingOperators > 1) {
1169 reason =
"blacklisted interference term!";
1175 for (
size_t j = 0; j < morphfunc[i].size(); ++j) {
1176 const int exponent = morphfunc[i][j];
1180 for (
int k = 0; k < exponent; ++k) {
1189 reason =
"coupling " +
cname +
" was listed as leading-order-only";
1192 if (!isZero && couplingsZero[j]) {
1194 reason =
"coupling " +
cname +
" is zero!";
1199 bool removedByFlag =
false;
1201 for (
auto itr : flags) {
1205 TString sval(obj->getStringAttribute(
"NewPhysics"));
1206 int val = atoi(sval);
1208 if (flagsZero.find(obj->GetName()) != flagsZero.end() && flagsZero.at(obj->GetName())) {
1209 removedByFlag =
true;
1210 reason =
Form(
"flag %s is zero", obj->GetName());
1217 if (!isZero && !removedByFlag) {
1219 const auto name = std::string(mfname) +
"_pol" + std::to_string(i);
1220 formulas[i] = std::make_unique<RooProduct>(
name.c_str(), ::concatNames(ss,
" * ").c_str(), ss);
1231 const std::vector<std::vector<RooArgList *>> &diagrams,
RooArgList &couplings,
1232 const RooArgList &flags,
const std::vector<std::vector<std::string>> &nonInterfering)
1234 MorphFuncPattern morphfuncpattern;
1236 for (
const auto &vertices : diagrams) {
1238 ::fillFeynmanDiagram(
d, vertices, couplings);
1239 ::collectPolynomials(morphfuncpattern,
d);
1241 FormulaList retval = buildFormulas(
name, inputs, inputFlags, morphfuncpattern, couplings, flags, nonInterfering);
1242 if (retval.empty()) {
1243 std::stringstream errorMsgStream;
1245 <<
"no formulas are non-zero, check if any if your couplings is floating and missing from your param_cards!"
1247 const auto errorMsg = errorMsgStream.str();
1248 throw std::runtime_error(errorMsg);
1250 checkMatrix(inputs, retval);
1259 FormulaList &formulas,
const Matrix &inverse)
1263 for (
auto sampleit : inputParameters) {
1264 const std::string sample(sampleit.first);
1265 std::stringstream title;
1266 TString name_full(makeValidName(sample.c_str()));
1268 name_full.Append(
"_");
1269 name_full.Append(fname);
1270 name_full.Prepend(
"w_");
1275 auto sampleformula = std::make_unique<RooLinearCombination>(name_full.Data());
1276 for (
auto const &formulait : formulas) {
1278 sampleformula->add(val, formulait.second.get());
1281 weights.addOwned(std::move(sampleformula));
1286inline std::map<std::string, std::string>
1291 std::map<std::string, std::string> weights;
1292 for (
auto sampleit : inputParameters) {
1293 const std::string sample(sampleit.first);
1294 std::stringstream str;
1297 for (
auto const &formulait : formulas) {
1298 double val(inverse(formulaidx, sampleidx));
1300 if (formulaidx > 0 && val > 0)
1302 str << val <<
"*(" << formulait.second->GetTitle() <<
")";
1306 weights[sample] = str.str();
1341 args.
add(*(it.second));
1356 const std::vector<std::vector<RooListProxy *>> &diagramProxyList,
1357 const std::vector<std::vector<std::string>> &nonInterfering,
const RooArgList &flags)
1360 std::vector<std::vector<RooArgList *>> diagrams;
1361 for (
const auto &diagram : diagramProxyList) {
1362 diagrams.emplace_back();
1365 diagrams.back().emplace_back(
vertex);
1369 _formulas = ::createFormulas(funcname, inputParameters, inputFlags, diagrams,
_couplings, flags, nonInterfering);
1374 template <
class List>
1380 Matrix matrix(buildMatrixT<Matrix>(inputParameters,
_formulas, operators, inputFlags, flags));
1381 if (
size(matrix) < 1) {
1382 std::cerr <<
"input matrix is empty, please provide suitable input samples!" << std::endl;
1387 double unityDeviation, largestWeight;
1388 inverseSanity(matrix, inverse, unityDeviation, largestWeight);
1394 oocxcoutW((
TObject *)0,
Eval) <<
"Warning: The matrix inversion seems to be unstable. This can "
1395 "be a result to input samples that are not sufficiently "
1396 "different to provide any morphing power."
1398 }
else if (weightwarning) {
1400 "result to input samples that are not sufficiently different to "
1401 "provide any morphing power."
1405 "encoded in your samples to cross-check:"
1407 for (
auto sampleit : inputParameters) {
1408 const std::string sample(sampleit.first);
1411 setParams(sampleit.second, operators,
true);
1438 const std::map<std::string, int> &storage,
const RooArgList &physics,
1442 std::cerr <<
"invalid bin width given!" << std::endl;
1446 std::cerr <<
"invalid observable given!" << std::endl;
1460 for (
auto sampleit : inputParameters) {
1462 TString prodname(makeValidName(sampleit.first.c_str()));
1467 std::cerr <<
"unable to access physics object for " << prodname << std::endl;
1474 std::cerr <<
"unable to access weight object for " << prodname << std::endl;
1481 allowNegativeYields =
true;
1482 auto prod = std::make_unique<RooProduct>(prodname, prodname, prodElems);
1483 if (!allowNegativeYields) {
1484 auto maxname = std::string(prodname) +
"_max0";
1487 auto max = std::make_unique<RooFormulaVar>(maxname.c_str(),
"max(0," + prodname +
")", prodset);
1488 max->addOwnedComponents(std::move(prod));
1489 sumElements.
addOwned(std::move(max));
1491 sumElements.
addOwned(std::move(prod));
1493 scaleElements.
add(*(binWidth));
1498 _sumFunc = make_unique<RooRealSumFunc>(
Form(
"%s_morphfunc",
name),
name, sumElements, scaleElements);
1501 std::cerr <<
"unable to access observable" << std::endl;
1502 _sumFunc.get()->addServer(*observable);
1504 std::cerr <<
"unable to access bin width" << std::endl;
1505 _sumFunc.get()->addServer(*binWidth);
1507 std::cerr <<
"no operators listed" << std::endl;
1508 _sumFunc.get()->addServerList(operators);
1510 std::cerr <<
"unable to access weight objects" << std::endl;
1511 _sumFunc.get()->addOwnedComponents(std::move(sumElements));
1512 _sumFunc.get()->addServerList(sumElements);
1513 _sumFunc.get()->addServerList(scaleElements);
1516 std::cout.precision(std::numeric_limits<double>::digits);
1533 if (obsName.empty()) {
1534 std::cerr <<
"Matrix inversion succeeded, but no observable was "
1535 "supplied. quitting..."
1543 setParams(func->
_flags, 1);
1547 setParams(func->
_flags, 1);
1570 setParams(func->
_flags, 1);
1574 setParams(func->
_flags, 1);
1604 return readMatrixFromFileT<TMatrixD>(fname);
1612 return readMatrixFromStreamT<TMatrixD>(stream);
1622 bool obsExists(
false);
1643 TH1 *hist = (
TH1 *)(inputExample);
1646 obs = obsOwner.get();
1650 auto obsOwner = std::make_unique<RooRealVar>(obsname, obsname, 0, 1);
1651 obs = obsOwner.get();
1657 if (strcmp(obsname, obs->
GetName()) != 0) {
1659 <<
" does not match expected name " << obsname << std::endl;
1664 auto binWidth = std::make_unique<RooRealVar>(sbw.
Data(), sbw.
Data(), 1.);
1666 binWidth->setVal(bw);
1667 binWidth->setConstant(
true);
1686 const size_t n(
size(cache->_inverse));
1688 const std::string sample(sampleit.first);
1691 if (!sampleformula) {
1696 for (
size_t formulaidx = 0; formulaidx <
n; ++formulaidx) {
1701 if (std::isnan(val)) {
1706 << val << std::endl;
1735 std::string obsName;
1748 << obsName <<
"'" << std::endl;
1750 auto obj = loadFromFileResidentFolder<TObject>(
file, folderNames.front(), obsName,
true);
1752 std::cerr <<
"unable to locate object '" << obsName <<
"' in folder '" << folderNames.front() <<
"'!"
1756 std::string classname = obj->ClassName();
1760 if (classname.find(
"TH1") != std::string::npos) {
1763 }
else if (classname.find(
"RooHistFunc") != std::string::npos ||
1764 classname.find(
"RooParamHistFunc") != std::string::npos ||
1765 classname.find(
"PiecewiseInterpolation") != std::string::npos) {
1767 }
else if (classname.find(
"TParameter<double>") != std::string::npos) {
1769 }
else if (classname.find(
"TParameter<float>") != std::string::npos) {
1771 }
else if (classname.find(
"TPair") != std::string::npos) {
1775 std::cerr <<
"cannot morph objects of class '" <<
mode->GetName() <<
"'!" << std::endl;
1786 std::cout << param.first <<
" = " << param.second;
1788 std::cout <<
" (const)";
1789 std::cout << std::endl;
1801 std::cout << folder << std::endl;
1823 _operators(
"operators",
"set of operators", this),
_observables(
"observables",
"morphing observables", this),
1841 std::vector<RooListProxy *> vertices;
1843 vertices.push_back(
new RooListProxy(
"!couplings",
"set of couplings in the vertex",
this,
true,
false));
1855 std::vector<RooListProxy *> vertices;
1861 new RooListProxy(
"!production",
"set of couplings in the production vertex",
this,
true,
false));
1862 vertices.push_back(
new RooListProxy(
"!decay",
"set of couplings in the decay vertex",
this,
true,
false));
1883 std::stringstream
name;
1884 name <<
"noInteference";
1885 for (
auto c : nonInterfering) {
1889 for (
auto c : nonInterfering) {
1900 for (
size_t i = 0; i < nonInterfering.size(); ++i) {
1949 :
RooAbsReal(other,
name), _cacheMgr(other._cacheMgr, this), _scale(other._scale), _sampleMap(other._sampleMap),
1950 _physics(other._physics.GetName(), this, other._physics),
1951 _operators(other._operators.GetName(), this, other._operators),
1952 _observables(other._observables.GetName(), this, other._observables),
1953 _binWidths(other._binWidths.GetName(), this, other._binWidths), _flags{other._flags.GetName(), this, other._flags},
1954 _config(other._config)
1956 for (
size_t j = 0; j < other.
_diagrams.size(); ++j) {
1957 std::vector<RooListProxy *> diagram;
1958 for (
size_t i = 0; i < other.
_diagrams[j].size(); ++i) {
1960 diagram.push_back(list);
1987 : _cacheMgr(this, 10, true, true), _operators(
"operators",
"set of operators", this, true, false),
1988 _observables(
"observable",
"morphing observable", this, true, false),
1989 _binWidths(
"binWidths",
"set of bin width objects", this, true, false)
2013 FeynmanDiagram diagram;
2014 std::vector<bool> prod;
2015 std::vector<bool> dec;
2016 for (
int i = 0; i < nboth; ++i) {
2017 prod.push_back(
true);
2018 dec.push_back(
true);
2020 for (
int i = 0; i < nprod; ++i) {
2021 prod.push_back(
true);
2022 dec.push_back(
false);
2024 for (
int i = 0; i < ndec; ++i) {
2025 prod.push_back(
false);
2026 dec.push_back(
true);
2028 diagram.push_back(prod);
2029 diagram.push_back(dec);
2030 MorphFuncPattern morphfuncpattern;
2031 ::collectPolynomials(morphfuncpattern, diagram);
2032 return morphfuncpattern.size();
2041 for (
auto vertex : vertices) {
2042 extractOperators(*
vertex, operators);
2043 extractCouplings(*
vertex, couplings);
2045 FeynmanDiagram diagram;
2046 ::fillFeynmanDiagram(diagram, vertices, couplings);
2047 MorphFuncPattern morphfuncpattern;
2048 ::collectPolynomials(morphfuncpattern, diagram);
2049 return morphfuncpattern.size();
2055std::map<std::string, std::string>
2057 const std::vector<std::vector<std::string>> &vertices_str)
2059 std::stack<RooArgList> ownedVertices;
2060 std::vector<RooArgList *> vertices;
2062 for (
const auto &vtx : vertices_str) {
2063 ownedVertices.emplace();
2064 auto &
vertex = ownedVertices.top();
2065 for (
const auto &
c : vtx) {
2068 auto couplingOwner = std::make_unique<RooRealVar>(
c.c_str(),
c.c_str(), 1., 0., 10.);
2069 coupling = couplingOwner.get();
2070 couplings.
addOwned(std::move(couplingOwner));
2074 vertices.push_back(&
vertex);
2083std::map<std::string, std::string>
2085 const std::vector<RooArgList *> &vertices,
RooArgList &couplings)
2093std::map<std::string, std::string>
2095 const std::vector<RooArgList *> &vertices,
RooArgList &couplings,
2097 const std::vector<std::vector<std::string>> &nonInterfering)
2099 FormulaList formulas = ::createFormulas(
"", inputs, flagValues, {vertices}, couplings, flags, nonInterfering);
2101 extractOperators(couplings, operators);
2102 Matrix matrix(::buildMatrixT<Matrix>(inputs, formulas, operators, flagValues, flags));
2103 if (
size(matrix) < 1) {
2104 std::cerr <<
"input matrix is empty, please provide suitable input samples!" << std::endl;
2108 auto retval = buildSampleWeightStrings(inputs, formulas, inverse);
2116 const std::vector<RooArgList *> &vertices,
RooArgList &couplings,
2119 const std::vector<std::vector<std::string>> &nonInterfering)
2121 FormulaList formulas = ::createFormulas(
"", inputs, flagValues, {vertices}, couplings, flags, nonInterfering);
2123 extractOperators(couplings, operators);
2124 Matrix matrix(::buildMatrixT<Matrix>(inputs, formulas, operators, flagValues, flags));
2125 if (
size(matrix) < 1) {
2126 std::cerr <<
"input matrix is empty, please provide suitable input samples!" << std::endl;
2131 ::buildSampleWeights(retval, (
const char *)
nullptr , inputs, formulas, inverse);
2139 const std::vector<RooArgList *> &vertices,
RooArgList &couplings)
2154 coutE(
Eval) <<
"unable to retrieve morphing function" << std::endl;
2162 for (
auto itr : *args) {
2187 auto wname = std::string(
"w_") +
name +
"_" + this->
GetName();
2188 return dynamic_cast<RooAbsReal *
>(cache->_weights.find(wname.c_str()));
2206 auto weightName = std::string(
"w_") + sample.first +
"_" + this->
GetName();
2207 auto weight =
static_cast<RooAbsReal *
>(cache->_weights.find(weightName.c_str()));
2224 double val = obj->
getVal();
2227 double variation =
r.Gaus(1, z);
2228 obj->
setVal(val * variation);
2269 cache->_inverse =
m;
2287 coutE(
Caching) <<
"unable to create cache!" << std::endl;
2305 coutE(
Caching) <<
"unable to create cache!" << std::endl;
2329 cxcoutP(
Caching) <<
"creating cache from getCache function for " <<
this << std::endl;
2335 coutE(
Caching) <<
"unable to create cache!" << std::endl;
2359 if (value < param->getMin())
2483 auto paramhist = loadFromFileResidentFolder<TH1>(
file, foldername,
"param_card");
2484 setParams(paramhist.get(),
_operators,
false);
2494 const std::string
name(foldername);
2504 for (
auto itr : *list) {
2549 auto mf = std::make_unique<RooRealSumFunc>(*(this->
getFunc()));
2552 const int nbins = observable->
getBins();
2557 for (
int i = 0; i < nbins; ++i) {
2562 for (
auto itr : *args) {
2574 double weight = formula->
getVal();
2577 val += dhist.
weight() * weight;
2582 return hist.release();
2591 auto mf = std::make_unique<RooRealSumFunc>(*(this->
getFunc()));
2595 for (
auto itr : *args) {
2597 if (prod->
getVal() != 0) {
2611 std::string pname(paramname);
2613 bool isUsed =
false;
2615 double thisval = sample.second.at(pname);
2616 if (thisval != val) {
2632 std::string
cname(couplname);
2639 bool isUsed =
false;
2642 double thisval = coupling->
getVal();
2643 if (thisval != val) {
2668 return cache->_formulas.size();
2676 auto mf = std::make_unique<RooRealSumFunc>(*(this->
getFunc()));
2678 std::cerr <<
"Error: unable to retrieve morphing function" << std::endl;
2682 for (
auto itr : *args) {
2687 name.Prepend(
"phys_");
2691 double val = formula->
getVal();
2693 std::cout << formula->
GetName() <<
": " << val <<
" = " << formula->
GetTitle() << std::endl;
2713 return &(cache->_couplings);
2727 double val = var->
getVal();
2728 couplings[
name] = val;
2755 auto func = std::make_unique<RooRealSumFunc>(*(cache->_sumFunc));
2758 return std::make_unique<RooWrapperPdf>(
Form(
"pdf_%s", func->GetName()),
Form(
"pdf of %s", func->GetTitle()), *func);
2767 return cache->_sumFunc.get();
2784 return this->
createPdf()->expectedEvents(nset);
2794 return this->
createPdf()->expectedEvents(set);
2803 return createPdf()->expectedEvents(&nset);
2816 auto weightName = std::string(
"w_") + sample.first +
"_" + this->
GetName();
2817 auto weight =
static_cast<RooAbsReal *
>(cache->_weights.find(weightName.c_str()));
2834 double w = weight->getVal();
2835 unc2 += newunc2 *
w *
w;
2862 for (
auto flag :
_flags) {
2876 for (
auto c : couplings) {
2877 std::cout <<
c.first <<
": " <<
c.second << std::endl;
2911 std::cerr <<
"unable to acquire in-built function!" << std::endl;
2943 const char *rangeName)
const
2987 coutE(
Caching) <<
"unable to retrieve cache!" << std::endl;
2998 coutE(
Caching) <<
"unable to retrieve cache!" << std::endl;
3011 coutE(
Caching) <<
"unable to retrieve cache!" << std::endl;
3012 return cache->_condition;
3018std::unique_ptr<RooRatio>
3022 for (
auto it : nr) {
3025 for (
auto it : dr) {
3029 return make_unique<RooRatio>(
name, title, num, denom);
3037std::vector<std::string> asStringV(std::string
const &arg)
3039 std::vector<std::string> out;
3041 for (std::string &tok :
ROOT::Split(arg,
",{}",
true)) {
3042 if (tok[0] ==
'\'') {
3043 out.emplace_back(tok.substr(1, tok.size() - 2));
3045 throw std::runtime_error(
"Strings in factory expressions need to be in single quotes!");
3055 create(
RooFactoryWSTool &,
const char *typeName,
const char *instName, std::vector<std::string> args)
override;
3058std::string LMIFace::create(
RooFactoryWSTool &ft,
const char * ,
const char *instanceName,
3059 std::vector<std::string> args)
3062 const std::array<std::string, 4> funcArgs{{
"fileName",
"observableName",
"couplings",
"folders"}};
3063 std::map<string, string> mappedInputs;
3065 for (
unsigned int i = 1; i < args.size(); i++) {
3066 if (args[i].find(
"$fileName(") != 0 && args[i].find(
"$observableName(") != 0 &&
3067 args[i].find(
"$couplings(") != 0 && args[i].find(
"$folders(") != 0 && args[i].find(
"$NewPhysics(") != 0) {
3068 throw std::string(
Form(
"%s::create() ERROR: unknown token %s encountered", instanceName, args[i].c_str()));
3072 for (
unsigned int i = 0; i < args.size(); i++) {
3073 if (args[i].find(
"$NewPhysics(") == 0) {
3075 for (
const auto &subarg : subargs) {
3076 std::vector<std::string> parts =
ROOT::Split(subarg,
"=");
3077 if (parts.size() == 2) {
3078 ft.
ws().
arg(parts[0].c_str())->
setAttribute(
"NewPhysics", atoi(parts[1].c_str()));
3080 throw std::string(
Form(
"%s::create() ERROR: unknown token %s encountered, check input provided for %s",
3081 instanceName, subarg.c_str(), args[i].c_str()));
3085 if (subargs.size() == 1) {
3087 for (
auto const ¶m : funcArgs) {
3088 if (args[i].find(param) != string::npos)
3089 mappedInputs[param] = subargs[0];
3093 Form(
"Incorrect number of arguments in %s, have %d, expect 1", args[i].c_str(), (
Int_t)subargs.size()));
3098 config.
fileName = asStringV(mappedInputs[
"fileName"])[0];
3099 config.
observableName = asStringV(mappedInputs[
"observableName"])[0];
3100 config.
folderNames = asStringV(mappedInputs[
"folders"]);
3105 return instanceName;
RooCollectionProxy< RooArgList > RooListProxy
size_t size< TMatrixD >(const TMatrixD &mat)
void writeMatrixToStreamT(const MatrixT &matrix, std::ostream &stream)
write a matrix to a stream
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
static constexpr double morphUnityDeviation
Matrix makeSuperMatrix(const TMatrixD &in)
convert a TMatrixD into a Matrix
static constexpr double morphLargestWeight
void writeMatrixToFileT(const MatrixT &matrix, const char *fname)
write a matrix to a text file
double invertMatrix(const Matrix &matrix, Matrix &inverse)
TMatrixD makeRootMatrix(const Matrix &in)
convert a matrix into a TMatrixD
Matrix diagMatrix(size_t n)
create a new diagonal matrix of size n
void printMatrix(const TMatrixD &mat)
write a matrix
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void input
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char cname
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Option_t Option_t TPoint TPoint const char mode
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
TMatrixT< Double_t > TMatrixD
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
std::string & operator+=(std::string &left, const TString &right)
TTime operator*(const TTime &t1, const TTime &t2)
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
void Print(Option_t *options=nullptr) const override
Print the object to the defaultPrintStream().
void setStringAttribute(const Text_t *key, const Text_t *value)
Associate string 'value' to this object under key 'key'.
bool isConstant() const
Check if the "Constant" attribute is set.
bool addOwnedComponents(const RooAbsCollection &comps)
Take ownership of the contents of 'comps'.
const RefCountList_t & servers() const
List of all servers of this object.
void setValueDirty()
Mark the element dirty. This forces a re-evaluation when a value is requested.
bool getAttribute(const Text_t *name) const
Check if a named attribute is set. By default, all attributes are unset.
void setAttribute(const Text_t *name, bool value=true)
Set (default) or clear a named boolean attribute of this object.
virtual double * array() const =0
RooAbsCacheElement is the abstract base class for objects to be stored in RooAbsCache cache manager o...
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * first() const
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
Int_t numBins(const char *rangeName=nullptr) const override
virtual Int_t getBins(const char *name=nullptr) const
Get number of bins of currently defined range.
void setConstant(bool value=true)
virtual double getMax(const char *name=nullptr) const
Get maximum of currently defined range.
void setBin(Int_t ibin, const char *rangeName=nullptr) override
Set value to center of bin 'ibin' of binning 'rangeName' (or of default binning if no range is specif...
virtual double getMin(const char *name=nullptr) const
Get minimum of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
RooArgList is a container object that can hold multiple RooAbsArg objects.
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Class RooBinning is an implements RooAbsBinning in terms of an array of boundary values,...
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=nullptr)
Setter function without integration set.
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=nullptr, const TNamed *isetRangeName=nullptr)
Getter function without integration set.
bool addOwned(RooAbsArg &var, bool silent=false) override
Overloaded RooCollection_t::addOwned() method insert object into owning set and registers object as s...
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
The RooDataHist is a container class to hold N-dimensional binned data.
double weight(std::size_t i) const
Return weight of i-th bin.
double weightSquared(std::size_t i) const
Return squared weight sum of i-th bin.
const RooArgSet * get() const override
Get bin centre of current bin.
RooHistFunc implements a real-valued function sampled from a multidimensional histogram.
RooArgSet const & getHistObsList() const
RooDataHist & dataHist()
Return RooDataHist that is represented.
static RooLagrangianMorphFunc::CacheElem * createCache(const RooLagrangianMorphFunc *func)
create all the temporary objects required by the class
void buildMatrix(const RooLagrangianMorphFunc::ParamMap &inputParameters, const RooLagrangianMorphFunc::FlagMap &inputFlags, const List &flags)
build and invert the morphing matrix
static RooLagrangianMorphFunc::CacheElem * createCache(const RooLagrangianMorphFunc *func, const Matrix &inverse)
create all the temporary objects required by the class function variant with precomputed inverse matr...
std::unique_ptr< RooRealSumFunc > _sumFunc
RooArgList containedArgs(Action) override
retrieve the list of contained args
void operModeHook(RooAbsArg::OperMode) override
Interface for changes of operation mode.
void createComponents(const RooLagrangianMorphFunc::ParamMap &inputParameters, const RooLagrangianMorphFunc::FlagMap &inputFlags, const char *funcname, const std::vector< std::vector< RooListProxy * > > &diagramProxyList, const std::vector< std::vector< std::string > > &nonInterfering, const RooArgList &flags)
create the basic objects required for the morphing
void buildMorphingFunction(const char *name, const RooLagrangianMorphFunc::ParamMap &inputParameters, const std::map< std::string, int > &storage, const RooArgList &physics, bool allowNegativeYields, RooRealVar *observable, RooRealVar *binWidth)
build the final morphing function
Class RooLagrangianMorphing is a implementation of the method of Effective Lagrangian Morphing,...
bool isParameterConstant(const char *paramname) const
return true if the parameter with the given name is set constant, false otherwise
bool isBinnedDistribution(const RooArgSet &obs) const override
check if this PDF is a binned distribution in the given observable
int nPolynomials() const
return the number of samples in this morphing function
void setParameter(const char *name, double value)
set one parameter to a specific value
RooArgSet createWeights(const ParamMap &inputs, const std::vector< RooArgList * > &vertices, RooArgList &couplings, const FlagMap &inputFlags, const RooArgList &flags, const std::vector< std::vector< std::string > > &nonInterfering)
create only the weight formulas. static function for external usage.
ParamSet getMorphParameters() const
retrieve the parameter set
double evaluate() const override
call getVal on the internal function
void disableInterference(const std::vector< const char * > &nonInterfering)
disable interference between terms
RooProduct * getSumElement(const char *name) const
return the RooProduct that is the element of the RooRealSumPdfi corresponding to the given sample nam...
RooRealVar * getBinWidth() const
retrieve the histogram observable
void writeMatrixToFile(const TMatrixD &matrix, const char *fname)
write a matrix to a file
RooRealVar * getParameter(const char *name) const
retrieve the RooRealVar object incorporating the parameter with the given name
bool useCoefficients(const TMatrixD &inverse)
setup the morphing function with a predefined inverse matrix call this function before any other afte...
const RooArgSet * getParameterSet() const
get the set of parameters
TMatrixD readMatrixFromStream(std::istream &stream)
read a matrix from a stream
std::vector< std::vector< std::string > > _nonInterfering
std::vector< std::string > getSamples() const
return the vector of sample names, used to build the morph func
int countSamples(std::vector< RooArgList * > &vertices)
calculate the number of samples needed to morph a certain physics process
void setCacheAndTrackHints(RooArgSet &) override
Retrieve the matrix of coefficients.
ParamSet getCouplings() const
retrieve a set of couplings (-?-)
void printSampleWeights() const
print the current sample weights
std::map< const std::string, double > ParamSet
void writeMatrixToStream(const TMatrixD &matrix, std::ostream &stream)
write a matrix to a stream
std::map< const std::string, ParamSet > ParamMap
bool updateCoefficients()
retrive the new physics objects and update the weights in the morphing function
double _scale
The cache manager.
RooRealVar * getObservable() const
retrieve the histogram observable
int countContributingFormulas() const
count the number of formulas that correspond to the current parameter set
std::list< double > * plotSamplingHint(RooAbsRealLValue &, double, double) const override
retrieve the sample Hint
RooRealVar * getFlag(const char *name) const
retrieve the RooRealVar object incorporating the flag with the given name
void randomizeParameters(double z)
randomize the parameters a bit useful to test and debug fitting
bool isCouplingUsed(const char *couplname)
check if there is any morphing power provided for the given coupling morphing power is provided as so...
void readParameters(TDirectory *f)
read the parameters from the input file
double getScale()
get energy scale of the EFT expansion
double getCondition() const
Retrieve the condition of the coefficient matrix.
TMatrixD getMatrix() const
Retrieve the matrix of coefficients.
void printWeights() const
print the current sample weights
void printCouplings() const
print a set of couplings
TMatrixD readMatrixFromFile(const char *fname)
read a matrix from a text file
~RooLagrangianMorphFunc() override
default destructor
void printParameters() const
print the parameters and their current values
void printPhysics() const
print the current physics values
static std::unique_ptr< RooRatio > makeRatio(const char *name, const char *title, RooArgList &nr, RooArgList &dr)
Return the RooRatio form of products and denominators of morphing functions.
void setFlag(const char *name, double value)
set one flag to a specific value
TH1 * createTH1(const std::string &name)
retrieve a histogram output of the current morphing settings
double expectedUncertainty() const
return the expected uncertainty for the current parameter set
int nParameters() const
return the number of parameters in this morphing function
bool hasParameter(const char *paramname) const
check if a parameter of the given name is contained in the list of known parameters
bool checkObservables(const RooArgSet *nset) const override
check if observable exists in the RooArgSet (-?-)
bool hasCache() const
return true if a cache object is present, false otherwise
void printFlags() const
print the flags and their current values
void setScale(double val)
set energy scale of the EFT expansion
TMatrixD getInvertedMatrix() const
Retrieve the matrix of coefficients after inversion.
double analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Retrieve the matrix of coefficients.
RooAbsArg::CacheMode canNodeBeCached() const override
Retrieve the matrix of coefficients.
void updateSampleWeights()
update sample weight (-?-)
void setParameters(const char *foldername)
set the morphing parameters to those supplied in the sample with the given name
RooObjCacheManager _cacheMgr
bool isParameterUsed(const char *paramname) const
check if there is any morphing power provided for the given parameter morphing power is provided as s...
RooListProxy _observables
RooAbsPdf::ExtendMode extendMode() const
return extended mored capabilities
std::map< const std::string, FlagSet > FlagMap
bool forceAnalyticalInt(const RooAbsArg &arg) const override
Force analytical integration for the given observable.
void setParameterConstant(const char *paramname, bool constant) const
call setConstant with the boolean argument provided on the parameter with the given name
void disableInterferences(const std::vector< std::vector< const char * > > &nonInterfering)
disable interference between terms
void printSamples() const
print all the known samples to the console
double getParameterValue(const char *name) const
set one parameter to a specific value
void setup(bool ownParams=true)
setup this instance with the given set of operators and vertices if own=true, the class will own the ...
void printMetaArgs(std::ostream &os) const override
Retrieve the matrix of coefficients.
std::unique_ptr< RooWrapperPdf > createPdf() const
(currently similar to cloning the Pdf
RooAbsReal * getSampleWeight(const char *name)
retrieve the weight (prefactor) of a sample with the given name
std::map< std::string, std::string > createWeightStrings(const ParamMap &inputs, const std::vector< std::vector< std::string > > &vertices)
create only the weight formulas. static function for external usage.
Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &numVars, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Retrieve the mat.
std::vector< std::vector< RooListProxy * > > _diagrams
double expectedEvents() const
return the number of expected events for the current parameter set
std::map< std::string, int > _sampleMap
RooLagrangianMorphFunc::CacheElem * getCache() const
retrieve the cache object
RooRealVar * setupObservable(const char *obsname, TClass *mode, TObject *inputExample)
setup observable, recycle existing observable if defined
const RooArgList * getCouplingSet() const
get the set of couplings
RooRealSumFunc * getFunc() const
get the func
std::list< double > * binBoundaries(RooAbsRealLValue &, double, double) const override
retrieve the list of bin boundaries
void printEvaluation() const
print the contributing samples and their respective weights
bool writeCoefficients(const char *filename)
write the inverse matrix to a file
void collectInputs(TDirectory *f)
retrieve the physics inputs
void init()
initialise inputs required for the morphing function
RooLinearCombination is a class that helps perform linear combination of floating point numbers and p...
void setCoefficient(size_t idx, RooFit::SuperFloat c)
RooFit::SuperFloat getCoefficient(size_t idx)
A histogram function that assigns scale parameters to every bin.
const RooArgList & paramList() const
A RooProduct represents the product of a given set of RooAbsReal objects.
std::list< double > * binBoundaries(RooAbsRealLValue &, double, double) const override
Retrieve bin boundaries if this distribution is binned in obs.
void printMetaArgs(std::ostream &os) const override
Customized printing of arguments of a RooRealSumFunc to more intuitively reflect the contents of the ...
Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &numVars, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Variant of getAnalyticalIntegral that is also passed the normalization set that should be applied to ...
CacheMode canNodeBeCached() const override
double analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
bool isBinnedDistribution(const RooArgSet &obs) const override
Tests if the distribution is binned. Unless overridden by derived classes, this always returns false.
bool checkObservables(const RooArgSet *nset) const override
Overloadable function in which derived classes can implement consistency checks of the variables.
void setCacheAndTrackHints(RooArgSet &) override
std::list< double > * plotSamplingHint(RooAbsRealLValue &, double, double) const override
Interface for returning an optional hint for initial sampling points when constructing a curve projec...
bool forceAnalyticalInt(const RooAbsArg &arg) const override
RooRealVar represents a variable that can be changed from the outside.
void setVal(double value) override
Set value of variable to 'value'.
void setError(double value)
void setMin(const char *name, double value)
Set minimum of name range to given value.
const RooAbsBinning & getBinning(const char *name=nullptr, bool verbose=true, bool createOnTheFly=false) const override
Return binning definition with name.
void setBinning(const RooAbsBinning &binning, const char *name=nullptr)
Add given binning under name 'name' with this variable.
void setBins(Int_t nBins, const char *name=nullptr)
Create a uniform binning under name 'name' for this variable.
void setMax(const char *name, double value)
Set maximum of name range to given value.
bool import(const RooAbsArg &arg, const RooCmdArg &arg1=RooCmdArg(), const RooCmdArg &arg2=RooCmdArg(), const RooCmdArg &arg3=RooCmdArg(), const RooCmdArg &arg4=RooCmdArg(), const RooCmdArg &arg5=RooCmdArg(), const RooCmdArg &arg6=RooCmdArg(), const RooCmdArg &arg7=RooCmdArg(), const RooCmdArg &arg8=RooCmdArg(), const RooCmdArg &arg9=RooCmdArg())
Import a RooAbsArg object, e.g.
RooAbsArg * arg(RooStringView name) const
Return RooAbsArg with given name. A null pointer is returned if none is found.
Class to manage histogram axis.
const char * GetBinLabel(Int_t bin) const
Return label for bin.
TClass instances represent classes, structs and namespaces in the ROOT type system.
static TClass * GetClass(const char *name, Bool_t load=kTRUE, Bool_t silent=kFALSE)
Static method returning pointer to TClass of the specified class name.
Double_t GetCondition() const
Bool_t Invert(TMatrixD &inv)
For a matrix A(m,m), its inverse A_inv is defined as A * A_inv = A_inv * A = unit (m x m) Ainv is ret...
Describe directory structure in memory.
virtual TObject * Get(const char *namecycle)
Return pointer to object identified by namecycle.
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
<div class="legacybox"><h2>Legacy Code</h2> TFolder is a legacy interface: there will be no bug fixes...
TCollection * GetListOfFolders() const
virtual void SetOwner(Bool_t owner=kTRUE)
Set ownership.
TH1 is the base class of all histogram classes in ROOT.
virtual Int_t GetNbinsX() const
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1) override
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
const char * GetName() const override
Returns name of object.
const char * GetTitle() const override
Returns title of object.
Mother of all ROOT objects.
virtual const char * GetName() const
Returns name of object.
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Class used by TMap to store (key,value) pairs.
Named parameter, streamable and storable.
const AParamType & GetVal() const
Random number generator class based on M.
int CompareTo(const char *cs, ECaseCompare cmp=kExact) const
Compare a string to char *cs2.
const char * Data() const
TString & Append(const char *cs)
RooCmdArg Silence(bool flag=true)
RVec< PromoteType< T > > abs(const RVec< T > &v)
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
void(off) SmallVectorTemplateBase< T
VecExpr< UnaryOp< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
std::vector< std::string > Split(std::string_view str, std::string_view delims, bool skipEmpty=false)
Splits a string at each character in delims.
void init()
Inspect hardware capabilities, and load the optimal library for RooFit computations.
std::string observableName
std::vector< std::string > folderNames
std::vector< std::vector< const char * > > nonInterfering