Logo ROOT  
Reference Guide
RooLagrangianMorphFunc.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * Package: RooLagrangianMorphing *
4 * @(#)root/roofit:$Id$
5 * Authors: *
6 * Lydia Brenner (lbrenner@cern.ch), Carsten Burgard (cburgard@cern.ch) *
7 * Katharina Ecker (kecker@cern.ch), Adam Kaluza (akaluza@cern.ch) *
8 * *
9 * Copyright (c) 2000-2005, Regents of the University of California *
10 * and Stanford University. All rights reserved. *
11 * *
12 * Redistribution and use in source and binary forms, *
13 * with or without modification, are permitted according to the terms *
14 * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15 *****************************************************************************/
16
17/** \class RooLagrangianMorphFunc
18 \ingroup Roofit
19Class RooLagrangianMorphing is a implementation of the method of Effective
20Lagrangian Morphing, described in ATL-PHYS-PUB-2015-047.
21Effective Lagrangian Morphing is a method to construct a continuous signal
22model in the coupling parameter space. Basic assumption is that shape and
23cross section of a physical distribution is proportional to it's
24squared matrix element.
25The signal model is constructed by a weighted sum over N input distributions.
26The calculation of the weights is based on Matrix Elements evaluated for the
27different input scenarios.
28The number of input files depends on the number of couplings in production
29and decay vertices, and also whether the decay and production vertices
30describe the same process or not.
31**/
32
33// uncomment to force UBLAS multiprecision matrices
34// #define USE_UBLAS 1
35// #undef USE_UBLAS
36
37#include "Riostream.h"
38
39#include "RooAbsCollection.h"
40#include "RooArgList.h"
41#include "RooArgProxy.h"
42#include "RooArgSet.h"
43#include "RooBinning.h"
44#include "RooDataHist.h"
45#include "RooFormulaVar.h"
46#include "RooHistFunc.h"
49#include "RooParamHistFunc.h"
50#include "RooProduct.h"
51#include "RooRealVar.h"
52#include "RooStringVar.h"
53#include "RooWorkspace.h"
54#include "TFile.h"
55#include "TFolder.h"
56#include "TH1.h"
57#include "TMap.h"
58#include "TParameter.h"
59#include "TRandom3.h"
60// stl includes
61#include <algorithm>
62#include <cmath>
63#include <cstddef>
64#include <iostream>
65#include <limits>
66#include <map>
67#include <memory>
68#include <sstream>
69#include <stdexcept>
70#include <type_traits>
71#include <typeinfo>
72
73using namespace std;
75
76//#define _DEBUG_
77
78///////////////////////////////////////////////////////////////////////////////
79// PREPROCESSOR MAGIC /////////////////////////////////////////////////////////
80///////////////////////////////////////////////////////////////////////////////
81
82// various preprocessor helpers
83#define NaN std::numeric_limits<double>::quiet_NaN()
84
85constexpr static double morphLargestWeight = 10e7;
86constexpr static double morphUnityDeviation = 10e-6;
87
88///////////////////////////////////////////////////////////////////////////////
89// TEMPLATE MAGIC /////////////////////////////////////////////////////////////
90///////////////////////////////////////////////////////////////////////////////
91
92template <typename Test, template <typename...> class Ref>
93struct is_specialization : std::false_type {};
94
95template <template <typename...> class Ref, typename... Args>
96struct is_specialization<Ref<Args...>, Ref> : std::true_type {};
97
98///////////////////////////////////////////////////////////////////////////////
99// LINEAR ALGEBRA HELPERS /////////////////////////////////////////////////////
100///////////////////////////////////////////////////////////////////////////////
101
102////////////////////////////////////////////////////////////////////////////////
103/// retrieve the size of a square matrix
104
105template <class MatrixT> inline size_t size(const MatrixT &matrix);
106template <> inline size_t size<TMatrixD>(const TMatrixD &mat) {
107 return mat.GetNrows();
108}
109
110////////////////////////////////////////////////////////////////////////////////
111/// write a matrix to a stream
112
113template <class MatrixT>
114void writeMatrixToStreamT(const MatrixT &matrix, std::ostream &stream) {
115 if(!stream.good()) {
116 return;
117 }
118 for (size_t i = 0; i < size(matrix); ++i) {
119 for (size_t j = 0; j < size(matrix); ++j) {
120#ifdef USE_UBLAS
121 stream << std::setprecision(SuperFloatPrecision::digits10) << matrix(i, j)
122 << "\t";
123#else
124 stream << matrix(i, j) << "\t";
125#endif
126 }
127 stream << std::endl;
128 }
129}
130
131////////////////////////////////////////////////////////////////////////////////
132/// write a matrix to a text file
133
134template <class MatrixT>
135inline void writeMatrixToFileT(const MatrixT &matrix, const char *fname) {
136 std::ofstream of(fname);
137 if (!of.good()) {
138 cerr << "unable to read file '" << fname << "'!" << std::endl;
139 }
140 writeMatrixToStreamT(matrix, of);
141 of.close();
142}
143
144#ifdef USE_UBLAS
145
146// boost includes
147#pragma GCC diagnostic push
148#pragma GCC diagnostic ignored "-Wshadow"
149#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
150#include <boost/numeric/ublas/io.hpp>
151#include <boost/numeric/ublas/lu.hpp>
152#include <boost/numeric/ublas/matrix.hpp>
153#include <boost/numeric/ublas/matrix_expression.hpp>
154#include <boost/numeric/ublas/symmetric.hpp> //inc diag
155#include <boost/numeric/ublas/triangular.hpp>
156#include <boost/operators.hpp>
157
158#pragma GCC diagnostic pop
159
160typedef boost::numeric::ublas::matrix<SuperFloat> Matrix;
161
162////////////////////////////////////////////////////////////////////////////////
163/// write a matrix
164
165inline void printMatrix(const Matrix &mat) {
166 for (size_t i = 0; i < mat.size1(); ++i) {
167 for (size_t j = 0; j < mat.size2(); ++j) {
168 // std::cout << std::setprecision(SuperFloatPrecision::digits10) <<
169 // mat(i,j) << " ,\t";
170 }
171 std::cout << std::endl;
172 }
173}
174
175////////////////////////////////////////////////////////////////////////////////
176/// retrieve the size of a square matrix
177
178template <> inline size_t size<Matrix>(const Matrix &matrix) {
179 return matrix.size1();
180}
181
182////////////////////////////////////////////////////////////////////////////////
183/// create a new diagonal matrix of size n
184
185inline Matrix diagMatrix(size_t n) {
186 return boost::numeric::ublas::identity_matrix<SuperFloat>(n);
187}
188
189////////////////////////////////////////////////////////////////////////////////
190/// convert a matrix into a TMatrixD
191
192inline TMatrixD makeRootMatrix(const Matrix &in) {
193 size_t n = size(in);
194 TMatrixD mat(n, n);
195 for (size_t i = 0; i < n; ++i) {
196 for (size_t j = 0; j < n; ++j) {
197 mat(i, j) = double(in(i, j));
198 }
199 }
200 return mat;
201}
202
203////////////////////////////////////////////////////////////////////////////////
204/// convert a TMatrixD into a matrix
205
206inline Matrix makeSuperMatrix(const TMatrixD &in) {
207 size_t n = in.GetNrows();
208 Matrix mat(n, n);
209 for (size_t i = 0; i < n; ++i) {
210 for (size_t j = 0; j < n; ++j) {
211 mat(i, j) = double(in(i, j));
212 }
213 }
214 return mat;
215}
216
217inline Matrix operator+=(const Matrix &rhs) { return add(rhs); }
218inline Matrix operator*(const Matrix &m, const Matrix &otherM) {
219 return prod(m, otherM);
220}
221
222////////////////////////////////////////////////////////////////////////////////
223/// calculate the inverse of a matrix, returning the condition
224
225inline SuperFloat invertMatrix(const Matrix &matrix, Matrix &inverse) {
226 boost::numeric::ublas::permutation_matrix<size_t> pm(size(matrix));
227 SuperFloat mnorm = norm_inf(matrix);
228 Matrix lu(matrix);
229 try {
230 int res = lu_factorize(lu, pm);
231 if (res != 0) {
232 std::stringstream ss;
233 ::writeMatrixToStreamT(matrix, ss);
234 cxcoutP(Eval) << ss.str << std::endl;
235 }
236 // back-substitute to get the inverse
237 lu_substitute(lu, pm, inverse);
238 } catch (boost::numeric::ublas::internal_logic &error) {
239 // coutE(Eval) << "boost::numberic::ublas error: matrix is not invertible!"
240 // << std::endl;
241 }
242 SuperFloat inorm = norm_inf(inverse);
243 SuperFloat condition = mnorm * inorm;
244 return condition;
245}
246
247#else
248
249#include "TDecompLU.h"
251
252////////////////////////////////////////////////////////////////////////////////
253/// convert a matrix into a TMatrixD
254
255inline TMatrixD makeRootMatrix(const Matrix &in) { return TMatrixD(in); }
256
257////////////////////////////////////////////////////////////////////////////////
258/// convert a TMatrixD into a Matrix
259
260inline Matrix makeSuperMatrix(const TMatrixD &in) { return in; }
261
262////////////////////////////////////////////////////////////////////////////////
263/// create a new diagonal matrix of size n
264
265inline Matrix diagMatrix(size_t n) {
266 TMatrixD mat(n, n);
267 mat.UnitMatrix();
268 return mat;
269}
270
271////////////////////////////////////////////////////////////////////////////////
272/// write a matrix
273
274inline void printMatrix(const TMatrixD &mat) {
275 writeMatrixToStreamT(mat, std::cout);
276}
277
278////////////////////////////////////////////////////////////////////////////////
279// calculate the inverse of a matrix, returning the condition
280
281inline double invertMatrix(const Matrix &matrix, Matrix &inverse) {
282 TDecompLU lu(matrix);
283 bool status = lu.Invert(inverse);
284 // check if the matrix is invertible
285 if (!status) {
286 std::cerr << " matrix is not invertible!" << std::endl;
287 }
288 double condition = lu.GetCondition();
289 const size_t n = size(inverse);
290 // sanitize numeric problems
291 for (size_t i = 0; i < n; ++i)
292 for (size_t j = 0; j < n; ++j)
293 if (fabs(inverse(i, j)) < 1e-9)
294 inverse(i, j) = 0;
295 return condition;
296}
297#endif
298
299/////////////////////////////////////////////////////////////////////////////////
300// LOCAL FUNCTIONS AND DEFINITIONS
301// //////////////////////////////////////////////
302/////////////////////////////////////////////////////////////////////////////////
303/// anonymous namespace to prohibit use of these functions outside the class
304/// itself
305namespace {
306///////////////////////////////////////////////////////////////////////////////
307// HELPERS ////////////////////////////////////////////////////////////////////
308///////////////////////////////////////////////////////////////////////////////
309
310typedef std::vector<std::vector<bool>> FeynmanDiagram;
311typedef std::vector<std::vector<int>> MorphFuncPattern;
312typedef std::map<int, std::unique_ptr<RooAbsReal>> FormulaList;
313
314///////////////////////////////////////////////////////////////////////////////
315/// check if a std::string begins with the given character set
316
317// inline bool begins_with(const std::string &input, const std::string &match) {
318// return input.size() >= match.size() &&
319// equal(match.begin(), match.end(), input.begin());
320// }
321
322///////////////////////////////////////////////////////////////////////////////
323/// (-?-)
324
325inline TString makeValidName(const char *input) {
326 TString retval(input);
327 retval.ReplaceAll("/", "_");
328 retval.ReplaceAll("^", "");
329 retval.ReplaceAll("*", "X");
330 retval.ReplaceAll("[", "");
331 retval.ReplaceAll("]", "");
332 return retval;
333}
334
335//////////////////////////////////////////////////////////////////////////////
336/// concatenate the names of objects in a collection to a single string
337
338template <class List> std::string concatNames(const List &c, const char *sep) {
339 std::stringstream ss;
340 bool first = true;
341 for (auto itr : c) {
342 if (!first)
343 ss << sep;
344 ss << itr->GetName();
345 first = false;
346 }
347 return ss.str();
348}
349
350///////////////////////////////////////////////////////////////////////////////
351/// this is a workaround for the missing implicit conversion from
352/// SuperFloat<>double
353
354template <class A, class B> inline void assignElement(A &a, const B &b) {
355 a = static_cast<A>(b);
356}
357///////////////////////////////////////////////////////////////////////////////
358// read a matrix from a stream
359
360template <class MatrixT>
361inline MatrixT readMatrixFromStreamT(std::istream &stream) {
362 std::vector<std::vector<SuperFloat>> matrix;
363 std::vector<SuperFloat> line;
364 while (!stream.eof()) {
365 if (stream.peek() == '\n') {
366 stream.get();
367 stream.peek();
368 continue;
369 }
370 SuperFloat val;
371 stream >> val;
372 line.push_back(val);
373 while (stream.peek() == ' ' || stream.peek() == '\t') {
374 stream.get();
375 }
376 if (stream.peek() == '\n') {
377 matrix.push_back(line);
378 line.clear();
379 }
380 }
381 MatrixT retval(matrix.size(), matrix.size());
382 for (size_t i = 0; i < matrix.size(); ++i) {
383 if (matrix[i].size() != matrix.size()) {
384 std::cerr << "matrix read from stream doesn't seem to be square!"
385 << std::endl;
386 }
387 for (size_t j = 0; j < matrix[i].size(); ++j) {
388 assignElement(retval(i, j), matrix[i][j]);
389 }
390 }
391 return retval;
392}
393
394///////////////////////////////////////////////////////////////////////////////
395/// read a matrix from a text file
396
397template <class MatrixT> inline MatrixT readMatrixFromFileT(const char *fname) {
398 std::ifstream in(fname);
399 if (!in.good()) {
400 std::cerr << "unable to read file '" << fname << "'!" << std::endl;
401 }
402 MatrixT mat = readMatrixFromStreamT<MatrixT>(in);
403 in.close();
404 return mat;
405}
406
407///////////////////////////////////////////////////////////////////////////////
408/// convert a TH1* param hist into the corresponding ParamSet object
409
410template <class T>
411void readValues(std::map<const std::string, T> & myMap, TH1 *h_pc) {
412 if (h_pc) {
413 // loop over all bins of the param_card histogram
414 for (int ibx = 1; ibx <= h_pc->GetNbinsX(); ++ibx) {
415 // read the value of one parameter
416 const std::string s_coup(h_pc->GetXaxis()->GetBinLabel(ibx));
417 double coup_val = h_pc->GetBinContent(ibx);
418 // add it to the map
419 if (!s_coup.empty()) {
420 myMap[s_coup] = T(coup_val);
421 }
422 }
423 }
424}
425
426///////////////////////////////////////////////////////////////////////////////
427/// retrieve a param_hist from a certain subfolder 'name' of the file
428
429inline TH1F *getParamHist(TDirectory *file, const std::string &name,
430 const std::string &objkey = "param_card",
431 bool notFoundError = true) {
432 auto f_tmp = file->Get<TFolder>(name.c_str());
433 if (!f_tmp) {
434 std::cerr << "unable to retrieve folder '" << name << "' from file '"
435 << file->GetName() << "'!" << std::endl;
436 return nullptr;
437}
438 // retrieve the histogram param_card which should live directly in the folder
439 TH1F *h_pc = dynamic_cast<TH1F *>(f_tmp->FindObject(objkey.c_str()));
440 if (h_pc) {
441 return h_pc;
442 }
443 if (notFoundError) {
444 std::cerr << "unable to retrieve " << objkey << " histogram from folder '"
445 << name << "'" << std::endl;
446 }
447 return nullptr;
448}
449
450///////////////////////////////////////////////////////////////////////////////
451/// retrieve a ParamSet from a certain subfolder 'name' of the file
452
453template <class T>
454void readValues(std::map<const std::string, T> & myMap,
455 TDirectory *file, const std::string &name,
456 const std::string &key = "param_card", bool notFoundError = true) {
457 TH1F *h_pc = getParamHist(file, name, key, notFoundError);
458 readValues(myMap, h_pc);
459}
460
461///////////////////////////////////////////////////////////////////////////////
462/// retrieve the param_hists file and return a map of the parameter values
463/// by providing a list of names, only the param_hists of those subfolders are
464/// read leaving the list empty is interpreted as meaning 'read everything'
465
466template <class T>
467void readValues(std::map<const std::string, std::map<const std::string, T>> & inputParameters,
468 TDirectory *f, const std::vector<std::string> &names,
469 const std::string &key = "param_card", bool notFoundError = true) {
470 inputParameters.clear();
471 // if the list of names is empty, we assume that this means 'all'
472 // loop over all folders in the file
473 for (size_t i = 0; i < names.size(); i++) {
474 const std::string name(names[i]);
475 // actually read an individual param_hist
476 // std::cout << "reading " << key << " '" << name << "'!" << std::endl;
477 readValues(inputParameters[name], f, name, key, notFoundError);
478 }
479
480 // now the map is filled with all parameter values found for all samples
481}
482
483///////////////////////////////////////////////////////////////////////////////
484/// open the file and return a file pointer
485
486inline TDirectory *openFile(const std::string &filename) {
487 if (filename.empty()) {
488 return gDirectory;
489 } else {
490 TFile *file = TFile::Open(filename.c_str(), "READ");
491 if (!file || !file->IsOpen()) {
492 if (file)
493 delete file;
494 std::cerr << "could not open file '" << filename << "'!" << std::endl;
495 }
496 return file;
497 }
498}
499
500///////////////////////////////////////////////////////////////////////////////
501/// open the file and return a file pointer
502
503inline void closeFile(TDirectory * d) {
504 TFile *f = dynamic_cast<TFile *>(d);
505 if (f) {
506 f->Close();
507 delete f;
508 d = nullptr;
509 }
510}
511
512///////////////////////////////////////////////////////////////////////////////
513/// extract the operators from a single coupling
514template <class T2>
515inline void extractServers(const RooAbsArg &coupling, T2 &operators) {
516 int nservers = 0;
517 for (const auto server : coupling.servers()) {
518 extractServers(*server, operators);
519 nservers++;
520 }
521 if (nservers == 0) {
522 operators.add(coupling);
523 }
524}
525
526///////////////////////////////////////////////////////////////////////////////
527/// extract the operators from a list of couplings
528
529template <class T1, class T2,
530 typename std::enable_if<!is_specialization<T1, std::vector>::value,
531 T1>::type * = nullptr>
532inline void extractOperators(const T1 &couplings, T2 &operators) {
533 // coutD(InputArguments) << "extracting operators from
534 // "<<couplings.getSize()<<" couplings" << std::endl;
535 for (auto itr : couplings) {
536 extractServers(*itr, operators);
537 }
538}
539
540///////////////////////////////////////////////////////////////////////////////
541/// extract the operators from a list of vertices
542
543template <class T1, class T2,
544 typename std::enable_if<is_specialization<T1, std::vector>::value,
545 T1>::type * = nullptr>
546inline void extractOperators(const T1 &vec, T2 &operators) {
547 for (const auto &v : vec) {
548 extractOperators(v, operators);
549 }
550}
551
552///////////////////////////////////////////////////////////////////////////////
553/// extract the couplings from a given set and copy them to a new one
554
555template <class T1, class T2>
556inline void extractCouplings(const T1 &inCouplings, T2 &outCouplings) {
557 for (auto itr : inCouplings) {
558 if (!outCouplings.find(itr->GetName())) {
559 // coutD(InputArguments) << "adding parameter "<< obj->GetName() <<
560 // std::endl;
561 outCouplings.add(*itr);
562 }
563 }
564}
565
566///////////////////////////////////////////////////////////////////////////////
567/// find and, if necessary, create a parameter from a list
568
569template <class T>
570inline RooAbsArg &get(T &operators, const char *name, double defaultval = 0) {
571 RooAbsArg *kappa = operators.find(name);
572 if (kappa)
573 return *kappa;
574 RooRealVar *newKappa = new RooRealVar(name, name, defaultval);
575 newKappa->setConstant(false);
576 operators.add(*newKappa);
577 return *newKappa;
578}
579
580///////////////////////////////////////////////////////////////////////////////
581/// find and, if necessary, create a parameter from a list
582
583template <class T>
584inline RooAbsArg &get(T &operators, const std::string &name,
585 double defaultval = 0) {
586 return get(operators, name.c_str(), defaultval);
587}
588
589////////////////////////////////////////////////////////////////////////////////
590/// create a new coupling and add it to the set
591
592template <class T>
593inline void addCoupling(T &set, const TString &name, const TString &formula,
594 const RooArgList &components, bool isNP) {
595 if (!set.find(name)) {
596 RooFormulaVar *c = new RooFormulaVar(name, formula, components);
597 c->setAttribute("NP", isNP);
598 set.add(*c);
599 }
600}
601
602////////////////////////////////////////////////////////////////////////////////
603/// set parameter values first set all values to defaultVal (if value not
604/// present in param_card then it should be 0)
605
606inline bool setParam(RooRealVar *p, double val, bool force) {
607 // coutD(InputArguments) << "setparam for "<<p->GetName()<<" to "<<val <<
608 // std::endl;
609 bool ok = true;
610 if (val > p->getMax()) {
611 if (force) {
612 p->setMax(val);
613 } else {
614 std::cerr << ": parameter " << p->GetName() << " out of bounds: " << val
615 << " > " << p->getMax() << std::endl;
616 ok = false;
617 }
618 } else if (val < p->getMin()) {
619 if (force) {
620 p->setMin(val);
621 } else {
622 std::cerr << ": parameter " << p->GetName() << " out of bounds: " << val
623 << " < " << p->getMin() << std::endl;
624 ok = false;
625 }
626 }
627 if (ok)
628 p->setVal(val);
629 return ok;
630}
631
632////////////////////////////////////////////////////////////////////////////////
633/// set parameter values first set all values to defaultVal (if value not
634/// present in param_card then it should be 0)
635
636template <class T1, class T2> inline bool setParams(const T2 &args, T1 val) {
637 for (auto itr : args) {
638 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
639 if (!param)
640 continue;
641 setParam(param, val, true);
642 }
643 return true;
644}
645
646////////////////////////////////////////////////////////////////////////////////
647/// set parameter values first set all values to defaultVal (if value not
648/// present in param_card then it should be 0)
649
650template <class T1, class T2>
651inline bool setParams(const std::map<const std::string, T1> &point,
652 const T2 &args, bool force = false, T1 defaultVal = 0) {
653 bool ok = true;
654 for (auto itr : args) {
655 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
656 if (!param || param->isConstant())
657 continue;
658 ok = setParam(param, defaultVal, force) && ok;
659 }
660 // set all parameters to the values in the param_card histogram
661 for (auto paramit : point) {
662 // loop over all the parameters
663 const std::string param(paramit.first);
664 // retrieve them from the map
665 RooRealVar *p = dynamic_cast<RooRealVar *>(args.find(param.c_str()));
666 if (!p)
667 continue;
668 // set them to their nominal value
669 ok = setParam(p, paramit.second, force) && ok;
670 }
671 return ok;
672}
673
674////////////////////////////////////////////////////////////////////////////////
675/// set parameter values first set all values to defaultVal (if value not
676/// present in param_card then it should be 0)
677
678template <class T>
679inline bool setParams(TH1 *hist, const T &args, bool force = false) {
680 bool ok = true;
681
682 for (auto itr : args) {
683 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
684 if (!param)
685 continue;
686 ok = setParam(param, 0., force) && ok;
687 }
688
689 // set all parameters to the values in the param_card histogram
690 TAxis *ax = hist->GetXaxis();
691 for (int i = 1; i <= ax->GetNbins(); ++i) {
692 // loop over all the parameters
693 RooRealVar *p = dynamic_cast<RooRealVar *>(args.find(ax->GetBinLabel(i)));
694 if (!p)
695 continue;
696 // set them to their nominal value
697 ok = setParam(p, hist->GetBinContent(i), force) && ok;
698 }
699 return ok;
700}
701
702////////////////////////////////////////////////////////////////////////////////
703/// create a set of parameters
704
705template <class T>
706inline RooLagrangianMorphFunc::ParamSet getParams(const T &parameters) {
708 for (auto itr : parameters) {
709 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
710 if (!param)
711 continue;
712 retval[param->GetName()] = param->getVal();
713 }
714 return retval;
715}
716
717////////////////////////////////////////////////////////////////////////////////
718/// collect the histograms from the input file and convert them to RooFit
719/// objects
720
721void collectHistograms(
722 const char *name, TDirectory *file, std::map<std::string, int> &list_hf,
723 RooArgList &physics, RooRealVar &var, const std::string &varname,
724 const RooLagrangianMorphFunc::ParamMap &inputParameters) {
725 // oocoutD((TObject*)0,InputArguments) << "building list of histogram
726 // functions" << std::endl;
727 bool binningOK = false;
728 for (auto sampleit : inputParameters) {
729 const std::string sample(sampleit.first);
730 auto folder = file->Get<TFolder>(sample.c_str());
731 if (!folder) {
732 std::cerr << "Error: unable to access data from folder '" << sample
733 << "'!" << std::endl;
734 continue;
735 }
736 TH1 *hist = dynamic_cast<TH1 *>(folder->FindObject(varname.c_str()));
737 if (!hist) {
738 std::stringstream errstr;
739 errstr << "Error: unable to retrieve histogram '" << varname
740 << "' from folder '" << sample << "'. contents are:";
741 TIter next(folder->GetListOfFolders()->begin());
742 TFolder *f;
743 while ((f = (TFolder *)next())) {
744 errstr << " " << f->GetName();
745 }
746 std::cerr << errstr.str() << std::endl;
747 return;
748 }
749
750 auto it = list_hf.find(sample);
751 if (it != list_hf.end()) {
752 RooHistFunc *hf = (RooHistFunc *)(physics.at(it->second));
753 hf->setValueDirty();
754 // commenting out To-be-resolved
755 // RooDataHist* dh = &(hf->dataHist());
756 // RooLagrangianMorphFunc::setDataHistogram(hist,&var,dh);
757 // RooArgSet vars;
758 // vars.add(var);
759 // dh->importTH1(vars,*hist,1.,kFALSE);
760 } else {
761 if (!binningOK) {
762 int n = hist->GetNbinsX();
763 std::vector<double> bins;
764 for (int i = 1; i < n + 1; ++i) {
765 bins.push_back(hist->GetBinLowEdge(i));
766 }
767 bins.push_back(hist->GetBinLowEdge(n) + hist->GetBinWidth(n));
768 var.setBinning(RooBinning(n, &(bins[0])));
769 }
770
771 // generate the mean value
772 TString histname = makeValidName(Form("dh_%s_%s", sample.c_str(), name));
773 TString funcname =
774 makeValidName(Form("phys_%s_%s", sample.c_str(), name));
775 RooArgSet vars;
776 vars.add(var);
777
778 RooDataHist *dh = new RooDataHist(histname.View(), histname.View(), vars, hist);
779 // add it to the list
780 RooHistFunc *hf = new RooHistFunc(funcname, funcname, var, *dh);
781 int idx = physics.getSize();
782 list_hf[sample] = idx;
783 physics.add(*hf);
784 assert(hf = static_cast<RooHistFunc *>(physics.at(idx)));
785 }
786 // std::cout << "found histogram " << hist->GetName() << " with integral "
787 // << hist->Integral() << std::endl;
788 }
789}
790
791////////////////////////////////////////////////////////////////////////////////
792/// collect the RooAbsReal objects from the input directory
793
794void collectRooAbsReal(
795 const char * /*name*/, TDirectory *file,
796 std::map<std::string, int> &list_hf, RooArgList &physics,
797 const std::string &varname,
798 const RooLagrangianMorphFunc::ParamMap &inputParameters) {
799 // cxcoutP(ObjectHandling) << "building list of RooAbsReal objects" <<
800 // std::endl;
801 for (auto sampleit : inputParameters) {
802 const std::string sample(sampleit.first);
803 auto folder = file->Get<TFolder>(sample.c_str());
804 if (!folder) {
805 std::cerr << "Error: unable to access data from folder '" << sample
806 << "'!" << std::endl;
807 continue;
808 }
809
810 RooAbsReal *obj =
811 dynamic_cast<RooAbsReal *>(folder->FindObject(varname.c_str()));
812 if (!obj) {
813 std::stringstream errstr;
814 std::cerr << "Error: unable to retrieve RooAbsArg '" << varname
815 << "' from folder '" << sample
816 << "'. contents are:" << std::endl;
817 TIter next(folder->GetListOfFolders()->begin());
818 TFolder *f;
819 while ((f = (TFolder *)next())) {
820 errstr << " " << f->GetName();
821 }
822 std::cerr << errstr.str() << std::endl;
823 return;
824 }
825 auto it = list_hf.find(sample);
826 if (it == list_hf.end()) {
827 int idx = physics.getSize();
828 list_hf[sample] = idx;
829 physics.add(*obj);
830 assert(obj == physics.at(idx));
831 }
832 }
833}
834
835////////////////////////////////////////////////////////////////////////////////
836/// collect the TParameter objects from the input file and convert them to
837/// RooFit objects
838
839template <class T>
840void collectCrosssections(
841 const char *name, TDirectory *file, std::map<std::string, int> &list_xs,
842 RooArgList &physics, const std::string &varname,
843 const RooLagrangianMorphFunc::ParamMap &inputParameters) {
844 // cxcoutP(ObjectHandling) << "building list of histogram functions" <<
845 // std::endl;
846 for (auto sampleit : inputParameters) {
847 const std::string sample(sampleit.first);
848 auto folder = file->Get<TFolder>(sample.c_str());
849 if (!folder) {
850 std::cerr << "unable to access data from folder '" << sample << "'!"
851 << std::endl;
852 return;
853 }
854 TObject *obj = folder->FindObject(varname.c_str());
855 TParameter<T> *xsection = nullptr;
856 TParameter<T> *error = nullptr;
857 TParameter<T> *p = dynamic_cast<TParameter<T> *>(obj);
858 if (p) {
859 xsection = p;
860 }
861 TPair *pair = dynamic_cast<TPair *>(obj);
862 if (pair) {
863 xsection = dynamic_cast<TParameter<T> *>(pair->Key());
864 error = dynamic_cast<TParameter<T> *>(pair->Value());
865 }
866 if (!xsection) {
867 std::stringstream errstr;
868 errstr << "Error: unable to retrieve cross section '" << varname
869 << "' from folder '" << sample << "'. contents are:";
870 TIter next(folder->GetListOfFolders()->begin());
871 TFolder *f;
872 while ((f = (TFolder *)next())) {
873 errstr << " " << f->GetName();
874 }
875 std::cerr << errstr.str() << std::endl;
876 return;
877 }
878
879 auto it = list_xs.find(sample.c_str());
880 RooRealVar *xs;
881 if (it != list_xs.end()) {
882 xs = (RooRealVar *)(physics.at(it->second));
883 xs->setVal(xsection->GetVal());
884 } else {
885 std::string objname = Form("phys_%s_%s", name, sample.c_str());
886 xs = new RooRealVar(objname.c_str(), objname.c_str(), xsection->GetVal());
887 xs->setConstant(true);
888 int idx = physics.getSize();
889 list_xs[sample] = idx;
890 physics.add(*xs);
891 assert(physics.at(idx) == xs);
892 }
893 if (error)
894 xs->setError(error->GetVal());
895 }
896}
897
898////////////////////////////////////////////////////////////////////////////////
899/// collect the TPair<TParameter,TParameter> objects from the input file and
900/// convert them to RooFit objects
901
902void collectCrosssectionsTPair(
903 const char *name, TDirectory *file, std::map<std::string, int> &list_xs,
904 RooArgList &physics, const std::string &varname,
905 const std::string &basefolder,
906 const RooLagrangianMorphFunc::ParamMap &inputParameters) {
907 auto folder = file->Get<TFolder>(basefolder.c_str());
908 TPair *pair = dynamic_cast<TPair *>(folder->FindObject(varname.c_str()));
909 TParameter<double> *xsec_double =
910 dynamic_cast<TParameter<double> *>(pair->Key());
911 if (xsec_double) {
912 collectCrosssections<double>(name, file, list_xs, physics, varname,
913 inputParameters);
914 } else {
915 TParameter<float> *xsec_float =
916 dynamic_cast<TParameter<float> *>(pair->Key());
917 if (xsec_float) {
918 collectCrosssections<float>(name, file, list_xs, physics, varname,
919 inputParameters);
920 } else {
921 std::cerr << "cannot morph objects of class 'TPair' if parameter is not "
922 "double or float!"
923 << std::endl;
924 }
925 }
926}
927
928///////////////////////////////////////////////////////////////////////////////
929// FORMULA CALCULATION ////////////////////////////////////////////////////////
930///////////////////////////////////////////////////////////////////////////////
931///////////////////////////////////////////////////////////////////////////////
932
933////////////////////////////////////////////////////////////////////////////////
934/// recursive function to determine polynomials
935
936void collectPolynomialsHelper(const FeynmanDiagram &diagram,
937 MorphFuncPattern &morphfunc,
938 std::vector<int> &term, int vertexid,
939 bool first) {
940 if (vertexid > 0) {
941 for (size_t i = 0; i < diagram[vertexid - 1].size(); ++i) {
942 if (!diagram[vertexid - 1][i])
943 continue;
944 std::vector<int> newterm(term);
945 newterm[i]++;
946 if (first) {
947 ::collectPolynomialsHelper(diagram, morphfunc, newterm, vertexid,
948 false);
949 } else {
950 ::collectPolynomialsHelper(diagram, morphfunc, newterm, vertexid - 1,
951 true);
952 }
953 }
954 } else {
955 bool found = false;
956 for (size_t i = 0; i < morphfunc.size(); ++i) {
957 bool thisfound = true;
958 for (size_t j = 0; j < morphfunc[i].size(); ++j) {
959 if (morphfunc[i][j] != term[j]) {
960 thisfound = false;
961 break;
962 }
963 }
964 if (thisfound) {
965 found = true;
966 break;
967 }
968 }
969 if (!found) {
970 morphfunc.push_back(term);
971 }
972 }
973}
974
975////////////////////////////////////////////////////////////////////////////////
976/// calculate the morphing function pattern based on a vertex map
977
978void collectPolynomials(MorphFuncPattern &morphfunc,
979 const FeynmanDiagram &diagram) {
980 int nvtx(diagram.size());
981 std::vector<int> term(diagram[0].size(), 0);
982
983 ::collectPolynomialsHelper(diagram, morphfunc, term, nvtx, true);
984}
985
986////////////////////////////////////////////////////////////////////////////////
987/// build a vertex map based on vertices and couplings appearing
988
989template <class List>
990inline void fillFeynmanDiagram(FeynmanDiagram &diagram,
991 const std::vector<List *> &vertices,
992 RooArgList &couplings) {
993 const int ncouplings = couplings.getSize();
994 // std::cout << "Number of couplings " << ncouplings << std::endl;
995 for (auto const& vertex : vertices) {
996 std::vector<bool> vertexCouplings(ncouplings, false);
997 int idx = -1;
998 RooAbsReal *coupling;
999 for (auto citr : couplings) {
1000 coupling = dynamic_cast<RooAbsReal *>(citr);
1001 idx++;
1002 if (!coupling) {
1003 std::cerr << "encountered invalid list of couplings in vertex!"
1004 << std::endl;
1005 return;
1006 }
1007 if (vertex->find(coupling->GetName())) {
1008 vertexCouplings[idx] = true;
1009 }
1010 }
1011 diagram.push_back(vertexCouplings);
1012 }
1013}
1014
1015////////////////////////////////////////////////////////////////////////////////
1016/// fill the matrix of coefficients
1017
1018template <class MatrixT, class T1, class T2>
1019inline MatrixT
1020buildMatrixT(const RooLagrangianMorphFunc::ParamMap &inputParameters,
1021 const FormulaList &formulas, const T1 &args,
1022 const RooLagrangianMorphFunc::FlagMap &flagValues,
1023 const T2 &flags) {
1024 const size_t dim = inputParameters.size();
1025 MatrixT matrix(dim, dim);
1026 int row = 0;
1027 for (auto sampleit : inputParameters) {
1028 const std::string sample(sampleit.first);
1029 // set all vars to value stored in input file
1030 if (!setParams<double>(sampleit.second, args, true, 0)) {
1031 // std::cout << "unable to set parameters for sample "<<sample<<"!" <<
1032 // std::endl;
1033 }
1034 auto flagit = flagValues.find(sample);
1035 if (flagit != flagValues.end() &&
1036 !setParams<int>(flagit->second, flags, true, 1)) {
1037 // std::cout << "unable to set parameters for sample "<<sample<<"!" <<
1038 // std::endl;
1039 }
1040 // loop over all the formulas
1041 int col = 0;
1042 for (auto const &formula : formulas) {
1043 if (!formula.second) {
1044 std::cerr << "Error: invalid formula encountered!" << std::endl;
1045 }
1046 matrix(row, col) = formula.second->getVal();
1047 // std::cout << formula.second->getVal() << " = " <<
1048 // formula.second->GetTitle() << " for " << sample << std::endl;
1049 col++;
1050 }
1051 row++;
1052 }
1053 return matrix;
1054}
1055
1056////////////////////////////////////////////////////////////////////////////////
1057/// check if the matrix is square
1058
1059inline void checkMatrix(const RooLagrangianMorphFunc::ParamMap &inputParameters,
1060 const FormulaList &formulas) {
1061 if (inputParameters.size() != formulas.size()) {
1062 std::stringstream ss;
1063 ss << "matrix is not square, consistency check failed: "
1064 << inputParameters.size() << " samples, " << formulas.size()
1065 << " expressions:" << std::endl;
1066 ss << "formulas: " << std::endl;
1067 for (auto const &formula : formulas) {
1068 ss << formula.second->GetTitle() << std::endl;
1069 }
1070 ss << "samples: " << std::endl;
1071 for (auto sample : inputParameters) {
1072 ss << sample.first << std::endl;
1073 }
1074 std::cerr << ss.str() << std::endl;
1075 }
1076}
1077
1078////////////////////////////////////////////////////////////////////////////////
1079/// check if the entries in the inverted matrix are sensible
1080
1081inline void inverseSanity(const Matrix &matrix, const Matrix &inverse,
1082 double &unityDeviation, double &largestWeight) {
1083 // cxcoutP(Eval) << "multiplying for sanity check" << std::endl;
1084 Matrix unity(inverse * matrix);
1085
1086 unityDeviation = 0.;
1087 largestWeight = 0.;
1088 const size_t dim = size(unity);
1089 for (size_t i = 0; i < dim; ++i) {
1090 for (size_t j = 0; j < dim; ++j) {
1091 if (inverse(i, j) > largestWeight) {
1092 largestWeight = (double)inverse(i, j);
1093 }
1094 if (fabs(unity(i, j) - static_cast<int>(i == j)) > unityDeviation) {
1095 unityDeviation = fabs((double)unity(i, j)) - static_cast<int>(i == j);
1096 }
1097 }
1098 }
1099 // std::cout << "found deviation of " << unityDeviation << ", largest weight
1100 // is " << largestWeight << "." << std::endl;
1101}
1102
1103////////////////////////////////////////////////////////////////////////////////
1104/// check for name conflicts between the input samples and an argument set
1105template <class List>
1106inline void
1107checkNameConflict(const RooLagrangianMorphFunc::ParamMap &inputParameters,
1108 List &args) {
1109 for (auto sampleit : inputParameters) {
1110 const std::string sample(sampleit.first);
1111 RooAbsArg *arg = args.find(sample.c_str());
1112 if (arg) {
1113 std::cerr << "detected name conflict: cannot use sample '" << sample
1114 << "' - a parameter with the same name of type '"
1115 << arg->ClassName() << "' is present in set '" << args.GetName()
1116 << "'!" << std::endl;
1117 }
1118 }
1119}
1120
1121////////////////////////////////////////////////////////////////////////////////
1122/// build the formulas corresponding to the given set of input files and
1123/// the physics process
1124
1125template <class List>
1126inline FormulaList
1127buildFormulas(const char *mfname,
1128 const RooLagrangianMorphFunc::ParamMap &inputParameters,
1129 const RooLagrangianMorphFunc::FlagMap &inputFlags,
1130 const MorphFuncPattern &morphfunc, const RooArgList &couplings,
1131 const List &flags, const std::vector<List *> &nonInterfering) {
1132 // example vbf hww:
1133 // Operators kSM, kHww, kAww, kHdwR,kHzz, kAzz
1134 // std::vector<bool> vertexProd = {true, true, true, true, true, true };
1135 // std::vector<bool> vertexDecay = {true, true, true, true, false,false};
1136 // diagram.push_back(vertexProd);
1137 // diagram.push_back(vertexDecay);
1138
1139 const int ncouplings = couplings.getSize();
1140 std::vector<bool> couplingsZero(ncouplings, true);
1141 std::map<TString, bool> flagsZero;
1142
1143 // std::cout << "number of couplings : " << ncouplings << std::endl;
1144 RooArgList operators;
1145 extractOperators(couplings, operators);
1146 size_t nOps = operators.getSize();
1147 // std::cout << "number of operators : " << nOps << std::endl;
1148
1149 for (auto sampleit : inputParameters) {
1150 // std::cout << "inside inputParameters loop" << std::endl;
1151 const std::string sample(sampleit.first);
1152 if (!setParams(sampleit.second, operators, true)) {
1153 std::cerr << "unable to set parameters for sample '" << sample << "'!"
1154 << std::endl;
1155 }
1156
1157 if ((int)nOps != (operators.getSize())) {
1158 std::cerr << "internal error, number of operators inconsistent!"
1159 << std::endl;
1160 }
1161
1162 RooAbsReal *obj0;
1163 int idx = 0;
1164
1165 for (auto itr1 : couplings) {
1166 obj0 = dynamic_cast<RooAbsReal *>(itr1);
1167 if (obj0->getVal() != 0) {
1168 // std::cout << obj0->GetName() << " is non-zero for sample " << sample
1169 // << " (idx=" << idx << ")!" << std::endl;
1170 couplingsZero[idx] = false;
1171 } else {
1172 // std::cout << obj0->GetName() << " is zero for sample " << sample << "
1173 // (idx=" << idx << ")!" << std::endl;
1174 }
1175 idx++;
1176 }
1177 }
1178
1179 for (auto itr2 : flags) {
1180 auto obj1 = dynamic_cast<RooAbsReal *>(itr2);
1181 int nZero = 0;
1182 int nNonZero = 0;
1183 for (auto sampleit : inputFlags) {
1184 const auto &flag = sampleit.second.find(obj1->GetName());
1185 if (flag != sampleit.second.end()) {
1186 if (flag->second == 0.)
1187 nZero++;
1188 else
1189 nNonZero++;
1190 // std::cout << "flag found " << obj1->GetName() << ", value = " <<
1191 // flag->second << std::endl;
1192 } else {
1193 // std::cout << "flag not found " << obj1->GetName() << std::endl;
1194 }
1195 }
1196 if (nZero > 0 && nNonZero == 0)
1197 flagsZero[obj1->GetName()] = true;
1198 else
1199 flagsZero[obj1->GetName()] = false;
1200 }
1201
1202 // {
1203 // int idx = 0;
1204 // RooAbsReal* obj2;
1205 // for(auto itr = couplings.begin(); itr != couplings.end(); ++itr){
1206 // obj2 = dynamic_cast<RooAbsReal*>(*itr);
1207 // if(couplingsZero[idx]){
1208 // //coutD(ObjectHandling) << obj2->GetName() << " is zero (idx=" <<
1209 // idx << ")" << std::endl;
1210 // } else {
1211 // //coutD(ObjectHandling) << obj2->GetName() << " is non-zero (idx="
1212 // << idx << ")" << std::endl;
1213 // }
1214 // idx++;
1215 // }
1216 // }
1217
1218 FormulaList formulas;
1219 for (size_t i = 0; i < morphfunc.size(); ++i) {
1220 RooArgList ss;
1221 bool isZero = false;
1222 std::string reason;
1223 // check if this is a blacklisted interference term
1224 for (const auto &group : nonInterfering) {
1225 int nInterferingOperators = 0;
1226 for (size_t j = 0; j < morphfunc[i].size(); ++j) {
1227 if (morphfunc[i][j] % 2 == 0)
1228 continue; // even exponents are not interference terms
1229 if (group->find(couplings.at(j)
1230 ->GetName())) { // if the coupling is part of a
1231 // "pairwise non-interfering group"
1232 nInterferingOperators++;
1233 }
1234 }
1235 if (nInterferingOperators > 1) {
1236 isZero = true;
1237 reason = "blacklisted interference term!";
1238 }
1239 }
1240 int nNP = 0;
1241 if (!isZero) {
1242 // prepare the term
1243 for (size_t j = 0; j < morphfunc[i].size(); ++j) {
1244 const int exponent = morphfunc[i][j];
1245 if (exponent == 0)
1246 continue;
1247 RooAbsReal *coupling = dynamic_cast<RooAbsReal *>(couplings.at(j));
1248 for (int k = 0; k < exponent; ++k) {
1249 ss.add(*coupling);
1250 if (coupling->getAttribute("NP")) {
1251 nNP++;
1252 }
1253 }
1254 std::string cname(coupling->GetName());
1255 if (coupling->getAttribute("LO") && exponent > 1) {
1256 isZero = true;
1257 reason = "coupling " + cname + " was listed as leading-order-only";
1258 }
1259 // mark the term as zero if any of the couplings are zero
1260 if (!isZero && couplingsZero[j]) {
1261 isZero = true;
1262 reason = "coupling " + cname + " is zero!";
1263 }
1264 }
1265 }
1266 // check and apply flags
1267 bool removedByFlag = false;
1268
1269 for (auto itr : flags) {
1270 auto obj = dynamic_cast<RooAbsReal *>(itr);
1271 if (!obj)
1272 continue;
1273 TString sval(obj->getStringAttribute("NP"));
1274 int val = atoi(sval);
1275 if (val == nNP) {
1276 if (flagsZero.find(obj->GetName()) != flagsZero.end() &&
1277 flagsZero.at(obj->GetName())) {
1278 removedByFlag = true;
1279 reason = Form("flag %s is zero", obj->GetName());
1280 }
1281 ss.add(*obj);
1282 }
1283 }
1284
1285 // create and add the formula
1286 if (!isZero && !removedByFlag) {
1287 // build the name
1288 const auto name = std::string(mfname) + "_pol" + std::to_string(i);
1289 // formulas[i] = new RooProduct(name.Data(),::concatNames(ss," *
1290 // ").c_str(),ss);
1291 formulas[i] = std::make_unique<RooProduct>(
1292 name.c_str(), ::concatNames(ss, " * ").c_str(), ss);
1293 std::cout << "creating formula " << name << ": "
1294 << formulas[i]->GetTitle() << std::endl;
1295 } else {
1296 // print a message and continue without doing anything
1297 // std::cout << "killing formula " << ::concatNames(ss," * ") << " because
1298 // " << reason << std::endl;
1299 }
1300 }
1301 return formulas;
1302}
1303
1304////////////////////////////////////////////////////////////////////////////////
1305/// create the weight formulas required for the morphing
1306
1307template <class T>
1308inline FormulaList
1309createFormulas(const char *name, const RooLagrangianMorphFunc::ParamMap &inputs,
1310 const RooLagrangianMorphFunc::FlagMap &inputFlags,
1311 const std::vector<std::vector<T *>> &diagrams,
1312 RooArgList &couplings, const T &flags,
1313 const std::vector<T *> &nonInterfering) {
1314 MorphFuncPattern morphfuncpattern;
1315 // std::cout << "building vertex map" << std::endl;
1316
1317 for (const auto &vertices : diagrams) {
1318 FeynmanDiagram d;
1319 // std::cout << "building vertex map" << std::endl;
1320 ::fillFeynmanDiagram(d, vertices, couplings);
1321 // std::cout << "collecting polynomials for diagram of size " << d.size() <<
1322 // std::endl;
1323 ::collectPolynomials(morphfuncpattern, d);
1324 }
1325 // coutD(ObjectHandling) << "building formulas" << std::endl;
1326 FormulaList retval = buildFormulas(name, inputs, inputFlags, morphfuncpattern,
1327 couplings, flags, nonInterfering);
1328 if (retval.size() == 0) {
1329 // std::cout << "no formulas are non-zero, check if any if your couplings is
1330 // floating and missing from your param_cards!" << std::endl;
1331 }
1332 // coutD(ObjectHandling) << "checking matrix consistency" << std::endl;
1333 checkMatrix(inputs, retval);
1334 return retval;
1335}
1336
1337////////////////////////////////////////////////////////////////////////////////
1338/// (-?-)
1339//
1340template <class T1>
1341inline void
1342buildSampleWeights(T1 &weights, const char *fname,
1343 const RooLagrangianMorphFunc::ParamMap &inputParameters,
1344 FormulaList &formulas, const Matrix &inverse) {
1345 int sampleidx = 0;
1346
1347 for (auto sampleit : inputParameters) {
1348 const std::string sample(sampleit.first);
1349 std::stringstream title;
1350 TString name_full(makeValidName(sample.c_str()));
1351 if (fname) {
1352 name_full.Append("_");
1353 name_full.Append(fname);
1354 name_full.Prepend("w_");
1355 }
1356
1357 int formulaidx = 0;
1358 // build the formula with the correct normalization
1359 RooLinearCombination *sampleformula = new RooLinearCombination(name_full.Data());
1360 for (auto const &formulait : formulas) {
1361 const SuperFloat val(inverse(formulaidx, sampleidx));
1362 sampleformula->add(val, formulait.second.get());
1363 formulaidx++;
1364 }
1365 weights.add(*sampleformula);
1366 sampleidx++;
1367 }
1368 // coutD(ObjectHandling) << "done building sample weights" << std::endl;
1369}
1370
1371inline std::map<std::string, std::string> buildSampleWeightStrings(
1372 const RooLagrangianMorphFunc::ParamMap &inputParameters,
1373 const FormulaList &formulas, const Matrix &inverse) {
1374 int sampleidx = 0;
1375 std::map<std::string, std::string> weights;
1376 for (auto sampleit : inputParameters) {
1377 const std::string sample(sampleit.first);
1378 std::stringstream str;
1379 int formulaidx = 0;
1380 // build the formula with the correct normalization
1381 for (auto const &formulait : formulas) {
1382 double val(inverse(formulaidx, sampleidx));
1383 if (val != 0.) {
1384 if (formulaidx > 0 && val > 0)
1385 str << " + ";
1386 str << val << "*(" << formulait.second->GetTitle() << ")";
1387 }
1388 formulaidx++;
1389 }
1390 weights[sample] = str.str();
1391 // std::cout << str.str() << std::endl;
1392 sampleidx++;
1393 }
1394 return weights;
1395}
1396} // namespace
1397
1398///////////////////////////////////////////////////////////////////////////////
1399// CacheElem magic ////////////////////////////////////////////////////////////
1400///////////////////////////////////////////////////////////////////////////////
1401
1403public:
1404 std::unique_ptr<RooRealSumFunc> _sumFunc = nullptr;
1406
1407 FormulaList _formulas;
1409
1413
1415 virtual void operModeHook(RooAbsArg::OperMode) override{};
1416
1417 //////////////////////////////////////////////////////////////////////////////
1418 /// retrieve the list of contained args
1419
1420 virtual RooArgList containedArgs(Action) override {
1421 RooArgList args(*_sumFunc);
1422 args.add(_weights);
1423 args.add(_couplings);
1424 for (auto const &it : _formulas) {
1425 args.add(*(it.second));
1426 }
1427 return args;
1428 }
1429
1430 //////////////////////////////////////////////////////////////////////////////
1431 // default destructor
1432
1433 virtual ~CacheElem() {}
1434
1435 //////////////////////////////////////////////////////////////////////////////
1436 /// create the basic objects required for the morphing
1437
1438 inline void
1440 const RooLagrangianMorphFunc::FlagMap &inputFlags,
1441 const char *funcname,
1442 const std::vector<std::vector<RooListProxy *>> &diagrams,
1443 const std::vector<RooListProxy *> &nonInterfering,
1444 const RooListProxy &flags) {
1445 RooArgList operators;
1446 // coutD(ObjectHandling) << "collecting couplings" << std::endl;
1447 for (const auto &diagram : diagrams) {
1448 for (const auto &vertex : diagram) {
1449 extractCouplings(*vertex, this->_couplings);
1450 }
1451 }
1452 extractOperators(this->_couplings, operators);
1453 this->_formulas =
1454 ::createFormulas(funcname, inputParameters, inputFlags, diagrams,
1455 this->_couplings, flags, nonInterfering);
1456 }
1457
1458 //////////////////////////////////////////////////////////////////////////////
1459 /// build and invert the morphing matrix
1460 template <class List>
1461 inline void
1463 const RooLagrangianMorphFunc::FlagMap &inputFlags,
1464 const List &flags) {
1465 RooArgList operators;
1466 extractOperators(this->_couplings, operators);
1467 Matrix matrix(buildMatrixT<Matrix>(inputParameters, this->_formulas,
1468 operators, inputFlags, flags));
1469 if (size(matrix) < 1) {
1470 std::cerr
1471 << "input matrix is empty, please provide suitable input samples!"
1472 << std::endl;
1473 }
1474 Matrix inverse(diagMatrix(size(matrix)));
1475
1476 double condition = (double)(invertMatrix(matrix, inverse));
1477 // coutD(Eval) << "Condition of the matrix :" << condition << std::endl;
1478
1479 double unityDeviation, largestWeight;
1480 inverseSanity(matrix, inverse, unityDeviation, largestWeight);
1481 bool weightwarning(largestWeight > morphLargestWeight ? true : false);
1482 bool unitywarning(unityDeviation > morphUnityDeviation ? true : false);
1483
1484 if (false) {
1485 if (unitywarning) {
1486 oocxcoutW((TObject *)0, Eval)
1487 << "Warning: The matrix inversion seems to be unstable. This can "
1488 "be a result to input samples that are not sufficiently "
1489 "different to provide any morphing power."
1490 << std::endl;
1491 } else if (weightwarning) {
1492 oocxcoutW((TObject *)0, Eval)
1493 << "Warning: Some weights are excessively large. This can be a "
1494 "result to input samples that are not sufficiently different to "
1495 "provide any morphing power."
1496 << std::endl;
1497 }
1498 oocxcoutW((TObject *)0, Eval) << " Please consider the couplings "
1499 "encoded in your samples to cross-check:"
1500 << std::endl;
1501 for (auto sampleit : inputParameters) {
1502 const std::string sample(sampleit.first);
1503 oocxcoutW((TObject *)0, Eval) << " " << sample << ": ";
1504 // set all vars to value stored in input file
1505 setParams(sampleit.second, operators, true);
1506 bool first = true;
1507 RooAbsReal *obj;
1508
1509 for (auto itr : this->_couplings) {
1510 obj = dynamic_cast<RooAbsReal *>(itr);
1511 if (!first)
1512 std::cerr << ", ";
1513 oocxcoutW((TObject *)0, Eval)
1514 << obj->GetName() << "=" << obj->getVal();
1515 first = false;
1516 }
1517 oocxcoutW((TObject *)0, Eval) << std::endl;
1518 }
1519 }
1520#ifndef USE_UBLAS
1521 this->_matrix.ResizeTo(matrix.GetNrows(), matrix.GetNrows());
1522 this->_inverse.ResizeTo(matrix.GetNrows(), matrix.GetNrows());
1523#endif
1524 this->_matrix = matrix;
1525 this->_inverse = inverse;
1526 this->_condition = condition;
1527 }
1528
1529 ////////////////////////////////////////////////////////////////////////////////
1530 /// build the final morphing function
1531
1533 const char *name, const RooLagrangianMorphFunc::ParamMap &inputParameters,
1534 const std::map<std::string, int> &storage, const RooArgList &physics,
1535 bool allowNegativeYields, RooRealVar *observable, RooRealVar *binWidth) {
1536 if (!binWidth) {
1537 std::cerr << "invalid bin width given!" << std::endl;
1538 return;
1539 }
1540 if (!observable) {
1541 std::cerr << "invalid observable given!" << std::endl;
1542 return;
1543 }
1544
1545 RooArgList operators;
1546 extractOperators(this->_couplings, operators);
1547
1548 // retrieve the weights
1549 ::buildSampleWeights(this->_weights, name, inputParameters, this->_formulas,
1550 this->_inverse);
1551
1552 // coutD(ObjectHandling) << "creating RooProducts" << std::endl;
1553 // build the products of element and weight for each sample
1554 size_t i = 0;
1555 RooArgList sumElements;
1556 RooArgList scaleElements;
1557 for (auto sampleit : inputParameters) {
1558 // for now, we assume all the lists are nicely ordered
1559 TString prodname(makeValidName(sampleit.first.c_str()));
1560
1561 RooAbsReal *obj =
1562 static_cast<RooAbsReal *>(physics.at(storage.at(prodname.Data())));
1563
1564 // obj->Print();
1565 if (!obj) {
1566 std::cerr << "unable to access physics object for " << prodname
1567 << std::endl;
1568 return;
1569 }
1570
1571 RooAbsReal *weight = static_cast<RooAbsReal *>(this->_weights.at(i));
1572
1573 // weight->Print();
1574 if (!weight) {
1575 std::cerr << "unable to access weight object for " << prodname
1576 << std::endl;
1577 return;
1578 }
1579 prodname.Append("_");
1580 prodname.Append(name);
1581 RooArgList prodElems(*weight, *obj);
1582
1583 allowNegativeYields = true;
1584 RooProduct *prod = new RooProduct(prodname, prodname, prodElems);
1585 if (!allowNegativeYields) {
1586 auto maxname = std::string(prodname) + "_max0";
1587 RooArgSet prodset(*prod);
1588
1589 RooFormulaVar *max =
1590 new RooFormulaVar(maxname.c_str(), "max(0," + prodname + ")", prodset);
1591 sumElements.add(*max);
1592 } else {
1593 sumElements.add(*prod);
1594 }
1595 scaleElements.add(*(binWidth));
1596 i++;
1597 }
1598
1599 // put everything together
1600 this->_sumFunc = make_unique<RooRealSumFunc>(
1601 Form("%s_morphfunc", name), name, sumElements, scaleElements);
1602
1603 // coutD(ObjectHandling) << "ownership handling" << std::endl;
1604 if (!observable)
1605 std::cerr << "unable to access observable" << std::endl;
1606 this->_sumFunc.get()->addServer(*observable);
1607 if (!binWidth)
1608 std::cerr << "unable to access bin width" << std::endl;
1609 this->_sumFunc.get()->addServer(*binWidth);
1610 if (operators.getSize() < 1)
1611 std::cerr << "no operators listed" << std::endl;
1612 this->_sumFunc.get()->addServerList(operators);
1613 if (this->_weights.getSize() < 1)
1614 std::cerr << "unable to access weight objects" << std::endl;
1615 this->_sumFunc.get()->addOwnedComponents(this->_weights);
1616 this->_sumFunc.get()->addOwnedComponents(sumElements);
1617 this->_sumFunc.get()->addServerList(sumElements);
1618 this->_sumFunc.get()->addServerList(scaleElements);
1619
1620#ifdef USE_UBLAS
1621 std::cout.precision(std::numeric_limits<double>::digits);
1622#endif
1623
1624 // oocxcoutP((TObject*)0,ObjectHandling) << morphfunc->Print() << std::endl;
1625
1626 // fill the this
1627 // this->_sumFunc = morphfunc;
1628 }
1629 //////////////////////////////////////////////////////////////////////////////
1630 /// create all the temporary objects required by the class
1631
1634 std::string obsName = func->getObservable()->GetName();
1635 // std::cout << "creating cache for RooLagrangianMorphFunc " << func <<
1636 // std::endl;
1637 RooLagrangianMorphFunc::ParamSet values = getParams(func->_operators);
1638
1641 cache->createComponents(
1642 func->_config.getParamCards(), func->_config.getFlagValues(),
1643 func->GetName(), func->_diagrams, func->_nonInterfering, func->_flags);
1644
1645 // cxcoutP(Eval) << "performing matrix operations" << std::endl;
1646 cache->buildMatrix(func->_config.getParamCards(),
1647 func->_config.getFlagValues(), func->_flags);
1648 if (obsName.size() == 0) {
1649 std::cerr << "Matrix inversion succeeded, but no observable was "
1650 "supplied. quitting..."
1651 << std::endl;
1652 return cache;
1653 }
1654
1655 // cxcoutP(ObjectHandling) << "building morphing function" << std::endl;
1657 << "observable: " << func->getObservable()->GetName() << std::endl;
1659 << "binWidth: " << func->getBinWidth()->GetName() << std::endl;
1660
1661 setParams(func->_flags, 1);
1662 cache->buildMorphingFunction(func->GetName(), func->_config.getParamCards(),
1663 func->_sampleMap, func->_physics,
1665 func->getObservable(), func->getBinWidth());
1666 setParams(values, func->_operators, true);
1667 setParams(func->_flags, 1);
1668 return cache;
1669 }
1670
1671 //////////////////////////////////////////////////////////////////////////////
1672 /// create all the temporary objects required by the class
1673 /// function variant with precomputed inverse matrix
1674
1676 createCache(const RooLagrangianMorphFunc *func, const Matrix &inverse) {
1677 // cxcoutP(Caching) << "creating cache for RooLagrangianMorphFunc = " <<
1678 // func << " with matrix" << std::endl;
1679 RooLagrangianMorphFunc::ParamSet values = getParams(func->_operators);
1680
1683 cache->createComponents(
1684 func->_config.getParamCards(), func->_config.getFlagValues(),
1685 func->GetName(), func->_diagrams, func->_nonInterfering, func->_flags);
1686
1687#ifndef USE_UBLAS
1688 cache->_inverse.ResizeTo(inverse.GetNrows(), inverse.GetNrows());
1689#endif
1690 cache->_inverse = inverse;
1691 cache->_condition = NaN;
1692
1693 // cxcoutP(ObjectHandling) << "building morphing function" << std::endl;
1694 setParams(func->_flags, 1);
1695 cache->buildMorphingFunction(func->GetName(), func->_config.getParamCards(),
1696 func->_sampleMap, func->_physics,
1698 func->getObservable(), func->getBinWidth());
1699 setParams(values, func->_operators, true);
1700 setParams(func->_flags, 1);
1701 return cache;
1702 }
1703};
1704
1705///////////////////////////////////////////////////////////////////////////////
1706// Class Implementation ///////////////////////////////////////////////////////
1707///////////////////////////////////////////////////////////////////////////////
1708
1709////////////////////////////////////////////////////////////////////////////////
1710/// insert an object into a workspace (wrapper for RooWorkspace::import)
1711
1713 const RooAbsReal *object) {
1714 if (ws && object)
1715 ws->import(*object, RooFit::RecycleConflictNodes());
1716}
1717
1718////////////////////////////////////////////////////////////////////////////////
1719/// insert an object into a workspace (wrapper for RooWorkspace::import)
1720
1722 RooAbsData *object) {
1723 if (ws && object)
1724 ws->import(*object, RooFit::RecycleConflictNodes());
1725}
1726
1727////////////////////////////////////////////////////////////////////////////////
1728/// insert this object into a workspace
1729
1732}
1733
1734////////////////////////////////////////////////////////////////////////////////
1735/// write a matrix to a file
1736
1738 const char *fname) {
1739 writeMatrixToFileT(matrix, fname);
1740}
1741
1742////////////////////////////////////////////////////////////////////////////////
1743/// write a matrix to a stream
1744
1746 std::ostream &stream) {
1747 writeMatrixToStreamT(matrix, stream);
1748}
1749
1750////////////////////////////////////////////////////////////////////////////////
1751/// read a matrix from a text file
1752
1754 return readMatrixFromFileT<TMatrixD>(fname);
1755}
1756
1757////////////////////////////////////////////////////////////////////////////////
1758/// read a matrix from a stream
1759
1761 return readMatrixFromStreamT<TMatrixD>(stream);
1762}
1763
1764////////////////////////////////////////////////////////////////////////////////
1765/// setup observable, recycle existing observable if defined
1766
1768 TClass *mode,
1769 TObject *inputExample) {
1770 // cxcoutP(ObjectHandling) << "setting up observable" << std::endl;
1771 RooRealVar *obs = nullptr;
1772 Bool_t obsExists(false);
1773 if (this->_observables.at(0) != 0) {
1774 obs = (RooRealVar *)this->_observables.at(0);
1775 obsExists = true;
1776 }
1777
1778 if (mode && mode->InheritsFrom(RooHistFunc::Class())) {
1779 obs = (RooRealVar *)(dynamic_cast<RooHistFunc *>(inputExample)
1780 ->getHistObsList()
1781 .first());
1782 obsExists = true;
1783 this->_observables.add(*obs);
1784 // ClassName()
1785 } else if (mode && mode->InheritsFrom(RooParamHistFunc::Class())) {
1786 obs = (RooRealVar *)(dynamic_cast<RooParamHistFunc *>(inputExample)
1787 ->paramList()
1788 .first());
1789 obsExists = true;
1790 this->_observables.add(*obs);
1791 }
1792
1793 // Note: "found!" will be printed if s2 is a substring of s1, both s1 and s2
1794 // are of type std::string. s1.find(s2)
1795 // obtain the observable
1796 if (!obsExists) {
1797 if (mode && mode->InheritsFrom(TH1::Class())) {
1798 // cxcoutP(ObjectHandling) << "getObservable: creating new multi-bin
1799 // observable object " << obsname << std::endl;
1800 TH1 *hist = (TH1 *)(inputExample);
1801 obs = new RooRealVar(obsname, obsname, hist->GetXaxis()->GetXmin(),
1802 hist->GetXaxis()->GetXmax());
1803 obs->setBins(hist->GetNbinsX());
1804 } else {
1805 // cxcoutP(ObjectHandling) << "getObservable: creating new single-bin
1806 // observable object " << obsname << std::endl;
1807 obs = new RooRealVar(obsname, obsname, 0, 1);
1808 obs->setBins(1);
1809 }
1810 this->_observables.add(*obs);
1811 } else {
1812 // cxcoutP(ObjectHandling) << "getobservable: recycling existing observable
1813 // object " << this->_observables.at(0) << std::endl;
1814 if (strcmp(obsname, obs->GetName()) != 0) {
1815 // coutW(ObjectHandling) << " name of existing observable " <<
1816 // this->_observables.at(0)->GetName() << " does not match expected name
1817 // " << obsname << std::endl ;
1818 }
1819 }
1820
1821 TString sbw = Form("binWidth_%s", makeValidName(obs->GetName()).Data());
1822 RooRealVar *binWidth = new RooRealVar(sbw.Data(), sbw.Data(), 1.);
1823 double bw = obs->numBins() / (obs->getMax() - obs->getMin());
1824 binWidth->setVal(bw);
1825 binWidth->setConstant(true);
1826 this->_binWidths.add(*binWidth);
1827
1828 return obs;
1829}
1830
1831//#ifndef USE_MULTIPRECISION_LC
1832//#pragma GCC diagnostic push
1833//#pragma GCC diagnostic ignored "-Wunused-parameter"
1834//#endif
1835
1836////////////////////////////////////////////////////////////////////////////////
1837/// update sample weight (-?-)
1838
1840 //#ifdef USE_MULTIPRECISION_LC
1841 int sampleidx = 0;
1842 auto cache = this->getCache(_curNormSet);
1843 const size_t n(size(cache->_inverse));
1844 for (auto sampleit : this->_config.getParamCards()) {
1845 const std::string sample(sampleit.first);
1846 // build the formula with the correct normalization
1847 RooLinearCombination *sampleformula = dynamic_cast<RooLinearCombination *>(
1848 this->getSampleWeight(sample.c_str()));
1849 if (!sampleformula) {
1850 // coutE(ObjectHandling) << Form("unable to access formula for sample
1851 // '%s'!",sample.c_str()) << std::endl;
1852 return;
1853 }
1855 << "updating formula for sample '" << sample << "'" << std::endl;
1856 for (size_t formulaidx = 0; formulaidx < n; ++formulaidx) {
1857 const SuperFloat val(cache->_inverse(formulaidx, sampleidx));
1858#ifdef USE_UBLAS
1859 if (val != val) {
1860#else
1861 if (std::isnan(val)) {
1862#endif
1863 coutE(ObjectHandling) << "refusing to propagate NaN!" << std::endl;
1864 }
1865 cxcoutP(ObjectHandling) << " " << formulaidx << ":"
1866 << sampleformula->getCoefficient(formulaidx)
1867 << " -> " << val << std::endl;
1868 sampleformula->setCoefficient(formulaidx, val);
1869 assert(sampleformula->getCoefficient(formulaidx) == val);
1870 }
1871 sampleformula->setValueDirty();
1872 ++sampleidx;
1873 }
1874 //#else
1875 // ERROR("updating sample weights currently not possible without boost!");
1876 //#endif
1877}
1878//#ifndef USE_MULTIPRECISION_LC
1879//#pragma GCC diagnostic pop
1880//#endif
1881
1882////////////////////////////////////////////////////////////////////////////////
1883/// read the parameters from the input file
1884
1886 readValues<double>(_paramCards, f, _folderNames, "param_card", true);
1887 readValues<int>(_flagValues, f, _folderNames, "flags", false);
1888}
1889
1890////////////////////////////////////////////////////////////////////////////////
1891/// retrieve the physics inputs
1892
1894 std::string obsName = this->_config.getObservableName();
1895 cxcoutP(InputArguments) << "initializing physics inputs from file "
1896 << file->GetName() << " with object name(s) '"
1897 << obsName << "'" << std::endl;
1898 auto foldernames = this->_config.getFolderNames();
1899 auto base = file->Get<TFolder>(foldernames[0].c_str());
1900 TObject *obj = base->FindObject(obsName.c_str());
1901 if (!obj) {
1902 std::cerr << "unable to locate object '" << obsName << "' in folder '"
1903 << base << "'!" << std::endl;
1904 return;
1905 }
1906 std::string classname = obj->ClassName();
1907 TClass *mode = TClass::GetClass(obj->ClassName());
1908
1909 RooRealVar *observable = this->setupObservable(obsName.c_str(), mode, obj);
1910 if (classname.find("TH1") != std::string::npos) {
1911 collectHistograms(this->GetName(), file, this->_sampleMap, this->_physics,
1912 *observable, obsName, this->_config.getParamCards());
1913 } else if (classname.find("RooHistFunc") != std::string::npos ||
1914 classname.find("RooParamHistFunc") != std::string::npos ||
1915 classname.find("PiecewiseInterpolation") != std::string::npos) {
1916 collectRooAbsReal(this->GetName(), file, this->_sampleMap, this->_physics,
1917 obsName, this->_config.getParamCards());
1918 } else if (classname.find("TParameter<double>") != std::string::npos) {
1919 collectCrosssections<double>(this->GetName(), file, this->_sampleMap,
1920 this->_physics, obsName,
1921 this->_config.getParamCards());
1922 } else if (classname.find("TParameter<float>") != std::string::npos) {
1923 collectCrosssections<float>(this->GetName(), file, this->_sampleMap,
1924 this->_physics, obsName,
1925 this->_config.getParamCards());
1926 } else if (classname.find("TPair") != std::string::npos) {
1927 collectCrosssectionsTPair(this->GetName(), file, this->_sampleMap,
1928 this->_physics, obsName, foldernames[0],
1929 this->_config.getParamCards());
1930 } else {
1931 std::cerr << "cannot morph objects of class '" << mode->GetName() << "'!"
1932 << std::endl;
1933 }
1934}
1935
1936////////////////////////////////////////////////////////////////////////////////
1937/// convert the RooArgList folders into a simple vector of std::string
1938
1940 for (auto const& folder : folders) {
1941 RooStringVar *var = dynamic_cast<RooStringVar *>(folder);
1942 const std::string sample(var ? var->getVal() : folder->GetName());
1943 if (sample.empty()) continue;
1944 this->_folderNames.push_back(sample);
1945 }
1946
1947 TDirectory *file = openFile(this->getFileName());
1948 TIter next(file->GetList());
1949 TObject *obj = nullptr;
1950 while ((obj = (TObject *)next())) {
1951 auto f = file->Get<TFolder>(obj->GetName());
1952 if (!f)
1953 continue;
1954 std::string name(f->GetName());
1955 if (name.empty())
1956 continue;
1957 this->_folderNames.push_back(name);
1958 }
1959 closeFile(file);
1960}
1961
1962////////////////////////////////////////////////////////////////////////////////
1963/// parameterised constructor
1965 extractCouplings(couplings, this->_couplings);
1966}
1967
1968////////////////////////////////////////////////////////////////////////////////
1969/// parameterised constructor
1971 const RooAbsCollection &decCouplings) {
1972 extractCouplings(prodCouplings, this->_prodCouplings);
1973 extractCouplings(decCouplings, this->_decCouplings);
1974}
1975
1976////////////////////////////////////////////////////////////////////////////////
1977/// config setter for couplings
1979 const RooAbsCollection &couplings) {
1980 extractCouplings(couplings, this->_couplings);
1981}
1982
1983////////////////////////////////////////////////////////////////////////////////
1984/// config setter for production and decay couplings
1986 const RooAbsCollection &prodCouplings,
1987 const RooAbsCollection &decCouplings) {
1988 extractCouplings(prodCouplings, this->_prodCouplings);
1989 extractCouplings(decCouplings, this->_decCouplings);
1990}
1991
1992////////////////////////////////////////////////////////////////////////////////
1993/// config setter for input file name
1995 this->_fileName = filename;
1996}
1997
1998////////////////////////////////////////////////////////////////////////////////
1999/// config setter for input folders
2001 for (auto itr : folderlist)
2002 this->_folderlist.add(*itr);
2003}
2004
2005////////////////////////////////////////////////////////////////////////////////
2006/// config setter for name of the the observable to be morphed
2008 this->_obsName = obsname;
2009}
2010
2011////////////////////////////////////////////////////////////////////////////////
2012/// config setter for name of the the observable to be morphed
2014 bool isallowNegativeYields) {
2015 this->_allowNegativeYields = isallowNegativeYields;
2016}
2017
2018////////////////////////////////////////////////////////////////////////////////
2019/// config setter for diagrams
2020template <class T>
2022 const std::vector<std::vector<T>> &diagrams) {
2023 for (size_t j = 0; j < diagrams.size(); ++j) {
2024 std::vector<RooArgList *> vertices;
2025 for (size_t i = 0; i < diagrams[j].size(); i++) {
2026 vertices.push_back(new RooArgList());
2027 vertices[i]->add(diagrams[j][i]);
2028 }
2029 this->_configDiagrams.push_back(vertices);
2030 }
2031}
2032
2033////////////////////////////////////////////////////////////////////////////////
2034/// config setter for vertices
2035template <class T>
2037 const std::vector<T> &vertices) {
2038 std::vector<std::vector<T>> diagrams;
2039 diagrams.push_back(vertices);
2040 this->setDiagrams(diagrams);
2041}
2042
2043////////////////////////////////////////////////////////////////////////////////
2044/// append the parameter map with a parmeter set
2045
2046//void RooLagrangianMorphFunc::Config::append(
2047// RooLagrangianMorphFunc::ParamMap &map, const char *str,
2048// RooLagrangianMorphFunc::ParamSet &set) {
2049// map[str] = set;
2050//}
2051
2052////////////////////////////////////////////////////////////////////////////////
2053/// set values to parameter set (-?-)
2054
2056 RooLagrangianMorphFunc::ParamSet &set, const char *str, double val) {
2057 set[str] = val;
2058}
2059
2060template void RooLagrangianMorphFunc::Config::setVertices<RooArgSet>(
2061 const std::vector<RooArgSet> &);
2062template void RooLagrangianMorphFunc::Config::setDiagrams<RooArgSet>(
2063 const std::vector<std::vector<RooArgSet>> &);
2064template void RooLagrangianMorphFunc::Config::setVertices<RooArgList>(
2065 const std::vector<RooArgList> &);
2066template void RooLagrangianMorphFunc::Config::setDiagrams<RooArgList>(
2067 const std::vector<std::vector<RooArgList>> &);
2068
2069/*
2070////////////////////////////////////////////////////////////////////////////////
2071/// config setter for vertices
2072void RooLagrangianMorphFunc::Config::disableInterference(const std::vector<const
2073char*>& nonInterfering){
2074 // disable interference between the listed operators
2075 std::stringstream name;
2076 name << "noInteference";
2077 for(auto c:nonInterfering){
2078 name << c;
2079 }
2080 RooArgList* p = new RooArgList(name.str().c_str());
2081 this->_nonInterfering.push_back(new RooArgList());
2082 for(auto c:nonInterfering){
2083 p->add(*(new RooStringVar(c,c,c)));
2084 }
2085}
2086
2087////////////////////////////////////////////////////////////////////////////////
2088/// config setter for vertices
2089void RooLagrangianMorphFunc::Config::disableInterferences(const
2090std::vector<std::vector<const char*> >& nonInterfering){
2091 // disable interferences between the listed groups of operators
2092 for(size_t i=0;i<nonInterfering.size();++i){
2093 this->disableInterference(nonInterfering[i]);
2094 }
2095}
2096*/
2097
2098
2099////////////////////////////////////////////////////////////////////////////////
2100/// print all the parameters and their values in the given sample to the console
2101
2102void RooLagrangianMorphFunc::printParameters(const char *samplename) const {
2103 for (const auto &param : this->_config.getParamCards().at(samplename)) {
2104 if (this->hasParameter(param.first.c_str())) {
2105 std::cout << param.first << " = " << param.second;
2106 if (this->isParameterConstant(param.first.c_str()))
2107 std::cout << " (const)";
2108 std::cout << std::endl;
2109 }
2110 }
2111}
2112
2113////////////////////////////////////////////////////////////////////////////////
2114/// print all the known samples to the console
2115
2117 // print all the known samples to the console
2118 for (auto folder : this->_config.getFolderNames()) {
2119 std::cout << folder << std::endl;
2120 }
2121}
2122
2123////////////////////////////////////////////////////////////////////////////////
2124/// print the current physics values
2125
2127 for (const auto &sample : this->_sampleMap) {
2128 RooAbsArg *phys = this->_physics.at(sample.second);
2129 if (!phys)
2130 continue;
2131 phys->Print();
2132 }
2133}
2134
2135////////////////////////////////////////////////////////////////////////////////
2136/// protected constructor with proper arguments
2137
2139 const char *title,
2140 const Config &config)
2141 : RooAbsReal(name, title), _cacheMgr(this, 10, kTRUE, kTRUE),
2142 _operators("operators", "set of operators", this, kTRUE, kFALSE),
2143 _observables("observables", "morphing observables", this, kTRUE, kFALSE),
2144 _binWidths("binWidths", "set of binWidth objects", this, kTRUE, kFALSE),
2145 _config(config), _curNormSet(0) {
2147 this->init();
2148 this->setup(false);
2149
2151}
2152
2153/*////////////////////////////////////////////////////////////////////////////////
2154/// protected constructor with proper arguments
2155
2156RooLagrangianMorphFunc::RooLagrangianMorphFunc(const char *name, const char
2157*title, const Config& config, bool allowNegativeYields) :
2158 RooLagrangianMorphFunc(name,title,config,"",folders,allowNegativeYields) {
2159// this->disableInterferences(this->_nonInterfering);
2160 this->setup(false);
2161 TRACE_CREATE
2162}
2163*/
2164////////////////////////////////////////////////////////////////////////////////
2165/// setup this instance with the given set of operators and vertices
2166/// if own=true, the class will own the operatorsemplate <class Base>
2167
2169 cxcoutP(Eval) << "setup(" << own << ") called" << std::endl;
2170 this->_ownParameters = own;
2171 auto diagrams = this->_config.getDiagrams();
2172
2173 if (diagrams.size() > 0) {
2174 RooArgList operators;
2175 for (auto const& v : diagrams) {
2176 for (const RooArgList *t : v) {
2177 extractOperators(*t, operators);
2178 }
2179 }
2180
2181 if (own) {
2182 this->_operators.addOwned(operators);
2183 } else {
2184 this->_operators.add(operators);
2185 }
2186
2187 for (size_t j = 0; j < diagrams.size(); ++j) {
2188 std::vector<RooListProxy *> diagram;
2189 for (size_t i = 0; i < diagrams[j].size(); ++i) {
2190 std::stringstream name;
2191 name << "!vertex" << i;
2192 std::stringstream title;
2193 title << "set of couplings in the vertex " << i;
2194 diagram.push_back(new RooListProxy(
2195 name.str().c_str(), title.str().c_str(), this, kTRUE, kFALSE));
2196 if (own) {
2197 diagram[i]->addOwned(*diagrams[j][i]);
2198 } else {
2199 diagram[i]->add(*diagrams[j][i]);
2200 }
2201 }
2202 // ownership of contents of diagrams[i][j]
2203 this->_diagrams.push_back(diagram);
2204 }
2205 }
2206
2207 else if (this->_config.getCouplings().size() > 0) {
2208 RooArgList operators;
2209 std::vector<RooListProxy *> vertices;
2210 cxcoutP(InputArguments) << "couplings provided" << std::endl;
2211 extractOperators(this->_config.getCouplings(), operators);
2212 vertices.push_back(new RooListProxy(
2213 "!couplings", "set of couplings in the vertex", this, kTRUE, kFALSE));
2214 if (own) {
2215 cxcoutP(InputArguments) << "adding own operators" << std::endl;
2216 this->_operators.addOwned(operators);
2217 vertices[0]->addOwned(this->_config.getCouplings());
2218 } else {
2219 cxcoutP(InputArguments) << "adding non-own operators" << std::endl;
2220 this->_operators.add(operators);
2221 vertices[0]->add(this->_config.getCouplings());
2222 }
2223 this->_diagrams.push_back(vertices);
2224 }
2225
2226 else if (this->_config.getProdCouplings().size() > 0 &&
2227 this->_config.getDecCouplings().size() > 0) {
2228 std::vector<RooListProxy *> vertices;
2229 RooArgList operators;
2230 cxcoutP(InputArguments) << "prod/dec couplings provided" << std::endl;
2231 extractOperators(this->_config.getProdCouplings(), operators);
2232 extractOperators(this->_config.getDecCouplings(), operators);
2233 vertices.push_back(new RooListProxy(
2234 "!production", "set of couplings in the production vertex", this, kTRUE,
2235 kFALSE));
2236 vertices.push_back(new RooListProxy(
2237 "!decay", "set of couplings in the decay vertex", this, kTRUE, kFALSE));
2238 if (own) {
2239 cxcoutP(InputArguments) << "adding own operators" << std::endl;
2240 this->_operators.addOwned(operators);
2241 vertices[0]->addOwned(this->_config.getProdCouplings());
2242 vertices[1]->addOwned(this->_config.getDecCouplings());
2243 } else {
2244 cxcoutP(InputArguments) << "adding non-own operators" << std::endl;
2245 this->_operators.add(operators);
2246 vertices[0]->add(this->_config.getProdCouplings());
2247 vertices[1]->add(this->_config.getDecCouplings());
2248 }
2249 this->_diagrams.push_back(vertices);
2250 }
2251}
2252
2253////////////////////////////////////////////////////////////////////////////////
2254/// (-?-)
2255
2257 std::string filename = this->_config.getFileName();
2258 TDirectory *file = openFile(filename.c_str());
2259 if (!file) {
2260 coutE(InputArguments) << "unable to open file '" << filename << "'!"
2261 << std::endl;
2262 return;
2263 }
2264 this->readParameters(file);
2265 checkNameConflict(this->_config.getParamCards(), this->_operators);
2266 this->collectInputs(file);
2267 closeFile(file);
2268 this->addServerList(this->_physics);
2269 cxcoutP(InputArguments) << "adding flags" << std::endl;
2270 RooRealVar *nNP0 = new RooRealVar("nNP0", "nNP0", 1., 0, 1.);
2271 nNP0->setStringAttribute("NP", "0");
2272 nNP0->setConstant(true);
2273 this->_flags.add(*nNP0);
2274 RooRealVar *nNP1 = new RooRealVar("nNP1", "nNP1", 1., 0, 1.);
2275 nNP1->setStringAttribute("NP", "1");
2276 nNP1->setConstant(true);
2277 this->_flags.add(*nNP1);
2278 RooRealVar *nNP2 = new RooRealVar("nNP2", "nNP2", 1., 0, 1.);
2279 nNP2->setStringAttribute("NP", "2");
2280 nNP2->setConstant(true);
2281 this->_flags.add(*nNP2);
2282 RooRealVar *nNP3 = new RooRealVar("nNP3", "nNP3", 1., 0, 1.);
2283 nNP3->setStringAttribute("NP", "3");
2284 nNP3->setConstant(true);
2285 this->_flags.add(*nNP3);
2286 RooRealVar *nNP4 = new RooRealVar("nNP4", "nNP4", 1., 0, 1.);
2287 nNP4->setStringAttribute("NP", "4");
2288 nNP4->setConstant(true);
2289 this->_flags.add(*nNP4);
2290}
2291
2292////////////////////////////////////////////////////////////////////////////////
2293/// copy constructor
2294
2296 const RooLagrangianMorphFunc &other, const char *name)
2297 : RooAbsReal(other, name), _cacheMgr(other._cacheMgr, this),
2298 _scale(other._scale), _sampleMap(other._sampleMap),
2299 _physics(other._physics.GetName(), this, other._physics),
2300 _operators(other._operators.GetName(), this, other._operators),
2301 _observables(other._observables.GetName(), this, other._observables),
2302 _binWidths(other._binWidths.GetName(), this, other._binWidths),
2303 _flags(other._flags.GetName(), this, other._flags),
2304 _config(other._config), _curNormSet(0) {
2305 std::cout << "COPY CONSTRUCTOR" << std::endl;
2306 for (size_t j = 0; j < other._diagrams.size(); ++j) {
2307 std::vector<RooListProxy *> diagram;
2308 for (size_t i = 0; i < other._diagrams[j].size(); ++i) {
2309 RooListProxy *list = new RooListProxy(other._diagrams[j][i]->GetName(),
2310 this, *(other._diagrams[j][i]));
2311 diagram.push_back(list);
2312 }
2313 this->_diagrams.push_back(diagram);
2314 }
2316}
2317
2318////////////////////////////////////////////////////////////////////////////////
2319/// set energy scale of the EFT expansion
2320
2321void RooLagrangianMorphFunc::setScale(double val) { this->_scale = val; }
2322
2323////////////////////////////////////////////////////////////////////////////////
2324/// get energy scale of the EFT expansion
2325
2327
2328////////////////////////////////////////////////////////////////////////////////
2329// default constructor
2330
2332 : _operators("operators", "set of operators", this, kTRUE, kFALSE),
2333 _observables("observable", "morphing observable", this, kTRUE, kFALSE),
2334 _binWidths("binWidths", "set of bin width objects", this, kTRUE, kFALSE) {
2335 static int counter(0);
2336 counter++;
2338}
2339
2340////////////////////////////////////////////////////////////////////////////////
2341/// default destructor
2342
2344
2345////////////////////////////////////////////////////////////////////////////////
2346/// cloning method
2347
2348TObject *RooLagrangianMorphFunc::clone(const char *newname) const {
2349 return new RooLagrangianMorphFunc(*this, newname);
2350}
2351
2352////////////////////////////////////////////////////////////////////////////////
2353/// print the authors information
2354
2356 std::cout << "\033[1mRooLagrangianMorphFunc\033[0m: a RooFit class for "
2357 "morphing physics distributions between configurations. authors:"
2358 << std::endl;
2359 std::cout << " "
2360 << "Lydia Brenner (lbrenner@cern.ch)" << std::endl;
2361 std::cout << " "
2362 << "Carsten Burgard (cburgard@cern.ch)" << std::endl;
2363 std::cout << " "
2364 << "Katharina Ecker (kecker@cern.ch)" << std::endl;
2365 std::cout << " "
2366 << "Adam Kaluza (akaluza@cern.ch)" << std::endl;
2367 std::cout << "please feel free to contact with questions and suggestions."
2368 << std::endl;
2369}
2370
2371////////////////////////////////////////////////////////////////////////////////
2372/// calculate the number of samples needed to morph a bivertex, 2-2 physics
2373/// process
2374
2375int RooLagrangianMorphFunc::countSamples(int nprod, int ndec, int nboth) {
2376 FeynmanDiagram diagram;
2377 std::vector<bool> prod;
2378 std::vector<bool> dec;
2379 for (int i = 0; i < nboth; ++i) {
2380 prod.push_back(true);
2381 dec.push_back(true);
2382 }
2383 for (int i = 0; i < nprod; ++i) {
2384 prod.push_back(true);
2385 dec.push_back(false);
2386 }
2387 for (int i = 0; i < ndec; ++i) {
2388 prod.push_back(false);
2389 dec.push_back(true);
2390 }
2391 diagram.push_back(prod);
2392 diagram.push_back(dec);
2393 MorphFuncPattern morphfuncpattern;
2394 ::collectPolynomials(morphfuncpattern, diagram);
2395 return morphfuncpattern.size();
2396}
2397
2398////////////////////////////////////////////////////////////////////////////////
2399/// calculate the number of samples needed to morph a certain physics process
2400/// usage:
2401/// countSamples ( { RooLagrangianMorphFunc::makeHCggfCouplings(),
2402/// RooLagrangianMorphFunc::makeHCHZZCouplings() } )
2403
2404int RooLagrangianMorphFunc::countSamples(std::vector<RooArgList *> &vertices) {
2405 RooArgList operators, couplings;
2406 for (auto vertex : vertices) {
2407 extractOperators(*vertex, operators);
2408 extractCouplings(*vertex, couplings);
2409 }
2410 FeynmanDiagram diagram;
2411 ::fillFeynmanDiagram(diagram, vertices, couplings);
2412 MorphFuncPattern morphfuncpattern;
2413 ::collectPolynomials(morphfuncpattern, diagram);
2414 return morphfuncpattern.size();
2415}
2416
2417////////////////////////////////////////////////////////////////////////////////
2418/// create TPair containers of the type expected by the RooLagrangianMorph
2419
2421 double unc) {
2422 TPair *v = new TPair(new TParameter<double>("xsection", xs),
2423 new TParameter<double>("uncertainty", unc));
2424 return v;
2425}
2426
2427////////////////////////////////////////////////////////////////////////////////
2428/// create only the weight formulas. static function for external usage.
2429
2430std::map<std::string, std::string> RooLagrangianMorphFunc::createWeightStrings(
2432 const std::vector<std::vector<std::string>> &vertices_str) {
2433 std::stack<RooArgList> ownedVertices;
2434 std::vector<RooArgList *> vertices;
2435 RooArgList couplings;
2436 for (const auto &vtx : vertices_str) {
2437 ownedVertices.emplace();
2438 auto& vertex = ownedVertices.top();
2439 for (const auto &c : vtx) {
2440 RooRealVar *coupling = (RooRealVar *)(couplings.find(c.c_str()));
2441 if (!coupling) {
2442 coupling = new RooRealVar(c.c_str(), c.c_str(), 1., 0., 10.);
2443 couplings.add(*coupling);
2444 }
2445 vertex.add(*coupling);
2446 }
2447 vertices.push_back(&vertex);
2448 }
2449 auto retval =
2450 RooLagrangianMorphFunc::createWeightStrings(inputs, vertices, couplings);
2451 return retval;
2452}
2453
2454////////////////////////////////////////////////////////////////////////////////
2455/// create only the weight formulas. static function for external usage.
2456
2457std::map<std::string, std::string> RooLagrangianMorphFunc::createWeightStrings(
2459 const std::vector<RooArgList *> &vertices, RooArgList &couplings) {
2460 return createWeightStrings(inputs, vertices, couplings, {}, {}, {});
2461}
2462
2463////////////////////////////////////////////////////////////////////////////////
2464/// create only the weight formulas. static function for external usage.
2465
2466std::map<std::string, std::string> RooLagrangianMorphFunc::createWeightStrings(
2468 const std::vector<RooArgList *> &vertices, RooArgList &couplings,
2469 const RooLagrangianMorphFunc::FlagMap &flagValues, const RooArgList &flags,
2470 const std::vector<RooArgList *> &nonInterfering) {
2471 FormulaList formulas = ::createFormulas("", inputs, flagValues, {vertices},
2472 couplings, flags, nonInterfering);
2473 RooArgSet operators;
2474 extractOperators(couplings, operators);
2475 Matrix matrix(
2476 ::buildMatrixT<Matrix>(inputs, formulas, operators, flagValues, flags));
2477 if (size(matrix) < 1) {
2478 std::cerr << "input matrix is empty, please provide suitable input samples!"
2479 << std::endl;
2480 }
2481 Matrix inverse(::diagMatrix(size(matrix)));
2482 double condition __attribute__((unused)) =
2483 (double)(invertMatrix(matrix, inverse));
2484 auto retval = buildSampleWeightStrings(inputs, formulas, inverse);
2485 return retval;
2486}
2487
2488////////////////////////////////////////////////////////////////////////////////
2489/// create only the weight formulas. static function for external usage.
2490
2493 const std::vector<RooArgList *> &vertices, RooArgList &couplings,
2494 const RooLagrangianMorphFunc::FlagMap &flagValues, const RooArgList &flags,
2495 const std::vector<RooArgList *> &nonInterfering) {
2496 FormulaList formulas = ::createFormulas("", inputs, flagValues, {vertices},
2497 couplings, flags, nonInterfering);
2498 RooArgSet operators;
2499 extractOperators(couplings, operators);
2500 Matrix matrix(
2501 ::buildMatrixT<Matrix>(inputs, formulas, operators, flagValues, flags));
2502 if (size(matrix) < 1) {
2503 std::cerr << "input matrix is empty, please provide suitable input samples!"
2504 << std::endl;
2505 }
2506 Matrix inverse(::diagMatrix(size(matrix)));
2507 double condition __attribute__((unused)) =
2508 (double)(invertMatrix(matrix, inverse));
2509 RooArgSet retval;
2510 ::buildSampleWeights(retval, (const char *)nullptr /* name */, inputs, formulas,
2511 inverse);
2512 return retval;
2513}
2514
2515////////////////////////////////////////////////////////////////////////////////
2516/// create only the weight formulas. static function for external usage.
2517
2520 const std::vector<RooArgList *> &vertices, RooArgList &couplings) {
2521 std::vector<RooArgList *> nonInterfering;
2522 RooArgList flags;
2523 FlagMap flagValues;
2525 inputs, vertices, couplings, flagValues, flags, nonInterfering);
2526}
2527
2528////////////////////////////////////////////////////////////////////////////////
2529/// return the RooProduct that is the element of the RooRealSumPdfi
2530/// corresponding to the given sample name
2531
2533 auto mf = this->getFunc();
2534 if (!mf) {
2535 coutE(Eval) << "unable to retrieve morphing function" << std::endl;
2536 return nullptr;
2537 }
2538 RooArgSet *args = mf->getComponents();
2539 TString prodname(name);
2540 prodname.Append("_");
2541 prodname.Append(this->GetName());
2542
2543 for (auto itr : *args) {
2544 RooProduct *prod = dynamic_cast<RooProduct *>(itr);
2545 if (!prod)
2546 continue;
2547 TString sname(prod->GetName());
2548 if (sname.CompareTo(prodname) == 0) {
2549 return prod;
2550 }
2551 }
2552 return nullptr;
2553}
2554////////////////////////////////////////////////////////////////////////////////
2555/// return the vector of sample names, used to build the morph func
2556
2557std::vector<std::string> RooLagrangianMorphFunc::getSamples() const {
2558 return this->_config.getFolderNames();
2559}
2560
2561////////////////////////////////////////////////////////////////////////////////
2562/// retrieve the weight (prefactor) of a sample with the given name
2563
2565 auto cache = this->getCache(_curNormSet);
2566 auto wname = std::string("w_") + name + "_" + this->GetName();
2567 return dynamic_cast<RooAbsReal *>(cache->_weights.find(wname.c_str()));
2568}
2569
2570////////////////////////////////////////////////////////////////////////////////
2571/// print the current sample weights
2572
2574 this->printSampleWeights();
2575}
2576
2577////////////////////////////////////////////////////////////////////////////////
2578/// print the current sample weights
2579
2581 auto *cache = this->getCache(this->_curNormSet);
2582 for (const auto &sample : this->_sampleMap) {
2583 auto weightName = std::string("w_") + sample.first + "_" + this->GetName();
2584 auto weight = static_cast<RooAbsReal *>(cache->_weights.find(weightName.c_str()));
2585 if (!weight)
2586 continue;
2587 std::cout << weight->GetName() << " = " << weight->GetTitle() << " = "
2588 << weight->getVal() << std::endl;
2589 }
2590}
2591
2592////////////////////////////////////////////////////////////////////////////////
2593/// randomize the parameters a bit
2594/// useful to test and debug fitting
2595
2597 RooRealVar *obj;
2598 TRandom3 r;
2599
2600 for (auto itr : _operators) {
2601 obj = dynamic_cast<RooRealVar *>(itr);
2602 double val = obj->getVal();
2603 if (obj->isConstant())
2604 continue;
2605 double variation = r.Gaus(1, z);
2606 obj->setVal(val * variation);
2607 }
2608}
2609
2610////////////////////////////////////////////////////////////////////////////////
2611/// retrive the new physics objects and update the weights in the morphing
2612/// function
2613
2615 auto cache = this->getCache(_curNormSet);
2616
2617 std::string filename = this->_config.getFileName();
2618 TDirectory *file = openFile(filename.c_str());
2619 if (!file) {
2620 coutE(InputArguments) << "unable to open file '" << filename << "'!"
2621 << std::endl;
2622 return false;
2623 }
2624 cxcoutP(InputArguments) << "reading parameter sets." << std::endl;
2625
2626 this->readParameters(file);
2627
2628 checkNameConflict(this->_config.getParamCards(), this->_operators);
2629 this->collectInputs(file);
2630
2631 cache->buildMatrix(this->_config.getParamCards(),
2632 this->_config.getFlagValues(), this->_flags);
2633 this->updateSampleWeights();
2634
2635 closeFile(file);
2636 return true;
2637}
2638
2639////////////////////////////////////////////////////////////////////////////////
2640/// setup the morphing function with a predefined inverse matrix
2641/// call this function *before* any other after creating the object
2642
2644 auto cache = static_cast<RooLagrangianMorphFunc::CacheElem *>(_cacheMgr.getObj(0, (RooArgSet *)0));
2645 Matrix m = makeSuperMatrix(inverse);
2646 if (cache) {
2647 std::string filename = this->_config.getFileName();
2648 cache->_inverse = m;
2649 TDirectory *file = openFile(filename.c_str());
2650 if (!file) {
2652 << "unable to open file '" << filename << "'!" << std::endl;
2653 return false;
2654 }
2655 cxcoutP(InputArguments) << "reading parameter sets." << std::endl;
2656
2657 this->readParameters(file);
2658 checkNameConflict(this->_config.getParamCards(), this->_operators);
2659 this->collectInputs(file);
2660
2661 // then, update the weights in the morphing function
2662 this->updateSampleWeights();
2663
2664 closeFile(file);
2665 } else {
2667 if (!cache)
2668 coutE(Caching) << "unable to create cache!" << std::endl;
2669 this->_cacheMgr.setObj(0, 0, cache, 0);
2670 }
2671 return true;
2672}
2673
2674////////////////////////////////////////////////////////////////////////////////
2675// setup the morphing function with a predefined inverse matrix
2676// call this function *before* any other after creating the object
2677
2679 auto cache = static_cast<RooLagrangianMorphFunc::CacheElem *>(_cacheMgr.getObj(0, (RooArgSet *)0));
2680 if (cache) {
2681 return false;
2682 }
2684 this, readMatrixFromFileT<Matrix>(filename));
2685 if (!cache)
2686 coutE(Caching) << "unable to create cache!" << std::endl;
2687 this->_cacheMgr.setObj(0, 0, cache, 0);
2688 return true;
2689}
2690
2691////////////////////////////////////////////////////////////////////////////////
2692/// write the inverse matrix to a file
2693
2695 auto cache = this->getCache(_curNormSet);
2696 if (!cache)
2697 return false;
2698 writeMatrixToFileT(cache->_inverse, filename);
2699 return true;
2700}
2701
2702////////////////////////////////////////////////////////////////////////////////
2703/// retrieve the cache object
2704
2707 auto cache = static_cast<RooLagrangianMorphFunc::CacheElem *>(_cacheMgr.getObj(0, (RooArgSet *)0));
2708 if (!cache) {
2709 cxcoutP(Caching) << "creating cache from getCache function for " << this
2710 << std::endl;
2711 cxcoutP(Caching) << "current storage has size " << this->_sampleMap.size()
2712 << std::endl;
2714 if (cache)
2715 this->_cacheMgr.setObj(0, 0, cache, 0);
2716 else
2717 coutE(Caching) << "unable to create cache!" << std::endl;
2718 }
2719 return cache;
2720}
2721
2722////////////////////////////////////////////////////////////////////////////////
2723/// return true if a cache object is present, false otherwise
2724
2726 return (bool)(_cacheMgr.getObj(0, (RooArgSet *)0));
2727}
2728
2729////////////////////////////////////////////////////////////////////////////////
2730/// set one parameter to a specific value
2731
2732void RooLagrangianMorphFunc::setParameter(const char *name, double value) {
2733 RooRealVar *param = this->getParameter(name);
2734 if (!param) {
2735 return;
2736 }
2737 if (value > param->getMax())
2738 param->setMax(value);
2739 if (value < param->getMin())
2740 param->setMin(value);
2741 param->setVal(value);
2742}
2743
2744////////////////////////////////////////////////////////////////////////////////
2745/// set one flag to a specific value
2746
2747void RooLagrangianMorphFunc::setFlag(const char *name, double value) {
2748 RooRealVar *param = this->getFlag(name);
2749 if (!param) {
2750 return;
2751 }
2752 param->setVal(value);
2753}
2754
2755////////////////////////////////////////////////////////////////////////////////
2756/// set one parameter to a specific value and range
2757
2758void RooLagrangianMorphFunc::setParameter(const char *name, double value,
2759 double min, double max) {
2760 RooRealVar *param = this->getParameter(name);
2761 if (!param) {
2762 return;
2763 }
2764 param->setMin(min);
2765 param->setMax(max);
2766 param->setVal(value);
2767}
2768
2769////////////////////////////////////////////////////////////////////////////////
2770/// set one parameter to a specific value and range
2771void RooLagrangianMorphFunc::setParameter(const char *name, double value,
2772 double min, double max,
2773 double error) {
2774 RooRealVar *param = this->getParameter(name);
2775 if (!param) {
2776 return;
2777 }
2778 param->setMin(min);
2779 param->setMax(max);
2780 param->setVal(value);
2781 param->setError(error);
2782}
2783
2784////////////////////////////////////////////////////////////////////////////////
2785/// return true if the parameter with the given name is set constant, false
2786/// otherwise
2787
2789 RooRealVar *param = this->getParameter(name);
2790 if (param) {
2791 return param->isConstant();
2792 }
2793 return true;
2794}
2795
2796////////////////////////////////////////////////////////////////////////////////
2797/// retrieve the RooRealVar object incorporating the parameter with the given
2798/// name
2800
2801 return dynamic_cast<RooRealVar *>(this->_operators.find(name));
2802}
2803
2804////////////////////////////////////////////////////////////////////////////////
2805/// retrieve the RooRealVar object incorporating the flag with the given name
2806
2808 return dynamic_cast<RooRealVar *>(this->_flags.find(name));
2809}
2810
2811////////////////////////////////////////////////////////////////////////////////
2812/// check if a parameter of the given name is contained in the list of known
2813/// parameters
2814
2816 return this->getParameter(name);
2817}
2818
2819////////////////////////////////////////////////////////////////////////////////
2820/// call setConstant with the boolean argument provided on the parameter with
2821/// the given name
2822
2824 bool constant) const {
2825 RooRealVar *param = this->getParameter(name);
2826 if (param) {
2827 return param->setConstant(constant);
2828 }
2829}
2830
2831////////////////////////////////////////////////////////////////////////////////
2832/// set one parameter to a specific value
2833
2835 RooRealVar *param = this->getParameter(name);
2836 if (param) {
2837 return param->getVal();
2838 }
2839 return 0;
2840}
2841
2842////////////////////////////////////////////////////////////////////////////////
2843/// set the morphing parameters to those supplied in the given param hist
2844
2846 setParams(paramhist, this->_operators, false);
2847}
2848
2849////////////////////////////////////////////////////////////////////////////////
2850/// set the morphing parameters to those supplied in the sample with the given
2851/// name
2852
2853void RooLagrangianMorphFunc::setParameters(const char *foldername) {
2854 std::string filename = this->_config.getFileName();
2855 TDirectory *file = openFile(filename.c_str());
2856 TH1 *paramhist = getParamHist(file, foldername);
2857 setParams(paramhist, this->_operators, false);
2858 closeFile(file);
2859}
2860
2861/////////////////////////////////////////////////////////////////////////////////
2862/// retrieve the morphing parameters associated to the sample with the given
2863/// name
2864
2866RooLagrangianMorphFunc::getMorphParameters(const char *foldername) const {
2867 const std::string name(foldername);
2868 return this->_config.getParamCards().at(name);
2869}
2870
2871////////////////////////////////////////////////////////////////////////////////
2872/// set the morphing parameters to those supplied in the list with the given
2873/// name
2874
2876 for (auto itr : *list) {
2877 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
2878 if (!param)
2879 continue;
2880 this->setParameter(param->GetName(), param->getVal());
2881 }
2882}
2883
2884////////////////////////////////////////////////////////////////////////////////
2885/// retrieve the histogram observable
2886
2888 if (this->_observables.getSize() < 1) {
2889 coutE(InputArguments) << "observable not available!" << std::endl;
2890 return nullptr;
2891 }
2892 return static_cast<RooRealVar *>(this->_observables.at(0));
2893}
2894
2895////////////////////////////////////////////////////////////////////////////////
2896/// retrieve the histogram observable
2897
2899 if (this->_binWidths.getSize() < 1) {
2900 coutE(InputArguments) << "bin width not available!" << std::endl;
2901 return nullptr;
2902 }
2903 return static_cast<RooRealVar *>(this->_binWidths.at(0));
2904}
2905
2906////////////////////////////////////////////////////////////////////////////////
2907/// retrieve a histogram output of the current morphing settings
2908
2910 /*RooFitResult *r*/) {
2911 return this->createTH1(name, false /*r*/);
2912}
2913
2914////////////////////////////////////////////////////////////////////////////////
2915/// retrieve a histogram output of the current morphing settings
2916
2918 bool correlateErrors
2919 /* RooFitResult *r*/) {
2920 auto mf = std::make_unique<RooRealSumFunc>(*(this->getFunc()));
2921 RooRealVar *observable = this->getObservable();
2922
2923 const int nbins = observable->getBins();
2924
2925 TH1 *hist = new TH1F(name.c_str(), name.c_str(), nbins,
2926 observable->getBinning().array());
2927
2928 //bool ownResult = !(bool)(r);
2929 RooArgSet *args = mf->getComponents();
2930 for (int i = 0; i < nbins; ++i) {
2931 observable->setBin(i);
2932 double val = 0;
2933 double unc2 = 0;
2934 double unc = 0;
2935 for (auto itr : *args) {
2936 RooProduct *prod = dynamic_cast<RooProduct *>(itr);
2937 if (!prod)
2938 continue;
2939 RooAbsArg *phys =
2940 prod->components().find(Form("phys_%s", prod->GetName()));
2941 RooHistFunc *hf = dynamic_cast<RooHistFunc *>(phys);
2942 if (!hf) {
2943 continue;
2944 }
2945 const RooDataHist &dhist = hf->dataHist();
2946 dhist.get(i);
2947 RooAbsReal *formula = dynamic_cast<RooAbsReal *>(
2948 prod->components().find(Form("w_%s", prod->GetName())));
2949 double weight = formula->getVal();
2950 unc2 += dhist.weightSquared() * weight * weight;
2951 unc += sqrt(dhist.weightSquared()) * weight;
2952 val += dhist.weight() * weight;
2953 }
2954 hist->SetBinContent(i + 1, val);
2955 hist->SetBinError(i + 1, correlateErrors ? unc : sqrt(unc2));
2956 }
2957 /*if (ownResult)
2958 delete r;*/
2959 return hist;
2960}
2961
2962////////////////////////////////////////////////////////////////////////////////
2963/// count the number of formulas that correspond to the current parameter set
2964
2966 int nFormulas = 0;
2967 auto mf = std::make_unique<RooRealSumFunc>(*(this->getFunc()));
2968 if (!mf)
2969 coutE(InputArguments) << "unable to retrieve morphing function"
2970 << std::endl;
2971 RooArgSet *args = mf->getComponents();
2972 for (auto itr : *args) {
2973 RooProduct *prod = dynamic_cast<RooProduct *>(itr);
2974 if (prod->getVal() != 0) {
2975 nFormulas++;
2976 }
2977 }
2978 return nFormulas;
2979}
2980
2981////////////////////////////////////////////////////////////////////////////////
2982/// check if there is any morphing power provided for the given parameter
2983/// morphing power is provided as soon as any two samples provide different,
2984/// non-zero values for this parameter
2985
2986bool RooLagrangianMorphFunc::isParameterUsed(const char *paramname) const {
2987 std::string pname(paramname);
2988 double val = 0;
2989 bool isUsed = false;
2990 for (const auto &sample : this->_config.getParamCards()) {
2991 double thisval = sample.second.at(pname);
2992 if (thisval != val) {
2993 if (val != 0)
2994 isUsed = true;
2995 val = thisval;
2996 }
2997 }
2998 return isUsed;
2999}
3000
3001////////////////////////////////////////////////////////////////////////////////
3002/// check if there is any morphing power provided for the given coupling
3003/// morphing power is provided as soon as any two samples provide
3004/// different, non-zero values for this coupling
3005
3006bool RooLagrangianMorphFunc::isCouplingUsed(const char *couplname) {
3007 std::string cname(couplname);
3008 const RooArgList *args = this->getCouplingSet();
3009 RooAbsReal *coupling = dynamic_cast<RooAbsReal *>(args->find(couplname));
3010 if (!coupling)
3011 return false;
3013 double val = 0;
3014 bool isUsed = false;
3015 for (const auto &sample : this->_config.getParamCards()) {
3016 this->setParameters(sample.second);
3017 double thisval = coupling->getVal();
3018 if (thisval != val) {
3019 if (val != 0)
3020 isUsed = true;
3021 val = thisval;
3022 }
3023 }
3024 this->setParameters(params);
3025 return isUsed;
3026}
3027
3028////////////////////////////////////////////////////////////////////////////////
3029/// return the number of parameters in this morphing function
3030
3032 return this->getParameterSet()->getSize();
3033}
3034
3035////////////////////////////////////////////////////////////////////////////////
3036/// return the number of samples in this morphing function
3037
3039 // return the number of samples in this morphing function
3040 auto cache = getCache(_curNormSet);
3041 return cache->_formulas.size();
3042}
3043
3044////////////////////////////////////////////////////////////////////////////////
3045/// print the contributing samples and their respective weights
3046
3048 auto mf = std::make_unique<RooRealSumFunc>(*(this->getFunc()));
3049 if (!mf) {
3050 std::cerr << "Error: unable to retrieve morphing function" << std::endl;
3051 return;
3052 }
3053 RooArgSet *args = mf->getComponents();
3054 for (auto itr : *args) {
3055 RooAbsReal *formula = dynamic_cast<RooAbsReal *>(itr);
3056 if (formula) {
3057 TString name(formula->GetName());
3058 name.Remove(0, 2);
3059 name.Prepend("phys_");
3060 if (!args->find(name.Data())) {
3061 continue;
3062 }
3063 double val = formula->getVal();
3064 if (val != 0) {
3065 std::cout << formula->GetName() << ": " << val << " = "
3066 << formula->GetTitle() << std::endl;
3067 }
3068 }
3069 }
3070}
3071
3072////////////////////////////////////////////////////////////////////////////////
3073/// get the set of parameters
3074
3076 return &(this->_operators);
3077}
3078
3079////////////////////////////////////////////////////////////////////////////////
3080/// get the set of couplings
3081
3083 auto cache = getCache(_curNormSet);
3084 return &(cache->_couplings);
3085}
3086
3087////////////////////////////////////////////////////////////////////////////////
3088/// retrieve a set of couplings (-?-)
3089
3092 for (auto obj : *(this->getCouplingSet())) {
3093 RooAbsReal *var = dynamic_cast<RooAbsReal *>(obj);
3094 if (!var)
3095 continue;
3096 const std::string name(var->GetName());
3097 double val = var->getVal();
3098 couplings[name] = val;
3099 }
3100 return couplings;
3101}
3102
3103////////////////////////////////////////////////////////////////////////////////
3104/// retrieve the parameter set
3105
3107 return getParams(this->_operators);
3108}
3109
3110////////////////////////////////////////////////////////////////////////////////
3111/// retrieve a set of couplings (-?-)
3112
3114 setParams(params, this->_operators, false);
3115}
3116
3117////////////////////////////////////////////////////////////////////////////////
3118/// (currently similar to cloning the Pdf
3119
3120std::unique_ptr<RooWrapperPdf> RooLagrangianMorphFunc::createPdf() const {
3121 auto cache = getCache(_curNormSet);
3122 auto func = std::make_unique<RooRealSumFunc>(*(cache->_sumFunc));
3123
3124 // create a wrapper on the roorealsumfunc
3125 return std::make_unique<RooWrapperPdf>(Form("pdf_%s", func->GetName()),
3126 Form("pdf of %s", func->GetTitle()),
3127 *func);
3128}
3129/*
3130////////////////////////////////////////////////////////////////////////////////
3131/// get a standalone clone of the pdf that does not depend on this object
3132
3133RooWrapperPdf* RooLagrangianMorphFunc::clonePdf() const
3134{
3135 auto cache = getCache(_curNormSet);
3136 RooRealSumFunc* func = cache->_sumFunc;
3137 // create a wrapper on the roorealsumfunc
3138 RooWrapperPdf* x = new RooWrapperPdf(Form("pdf_%s",func->GetName()),Form("pdf
3139of %s",func->GetTitle()), *func); return x;
3140}
3141*/
3142////////////////////////////////////////////////////////////////////////////////
3143/// get the func
3144
3146 auto cache = getCache(_curNormSet);
3147 return cache->_sumFunc.get();
3148}
3149
3150////////////////////////////////////////////////////////////////////////////////
3151/// get a standalone clone of the func that does not depend on this object
3152/* do we need a cloneFunc method ?
3153 i get segfault when i do,
3154 std::unique_ptr<RooRealSumFunc> morphfunc_ptjj1 = morphfunc_ptjj.cloneFunc();
3155 morphfunc_ptjj1.get()->Print();
3156
3157std::unique_ptr<RooRealSumFunc> RooLagrangianMorphFunc::cloneFunc() const
3158{
3159 auto cache = getCache(_curNormSet);
3160 auto orig = cache->_sumFunc.get();
3161 return
3162std::make_unique<RooRealSumFunc>(orig->GetName(),orig->GetTitle(),orig->funcList(),orig->coefList());
3163}
3164*/
3165////////////////////////////////////////////////////////////////////////////////
3166/// return extended mored capabilities
3167
3169 return this->createPdf()->extendMode();
3170}
3171
3172////////////////////////////////////////////////////////////////////////////////
3173/// return expected number of events for extended likelihood calculation,
3174/// this is the sum of all coefficients
3175
3177 return this->createPdf()->expectedEvents(nset);
3178}
3179
3180////////////////////////////////////////////////////////////////////////////////
3181/// return the number of expected events for the current parameter set
3182
3184 RooArgSet set;
3185 set.add(*this->getObservable());
3186 return this->createPdf()->expectedEvents(set);
3187}
3188
3189////////////////////////////////////////////////////////////////////////////////
3190/// return expected number of events for extended likelihood calculation,
3191/// this is the sum of all coefficients
3192
3194 return createPdf()->expectedEvents(&nset);
3195}
3196////////////////////////////////////////////////////////////////////////////////
3197/// return the expected uncertainty for the current parameter set
3198
3200 RooRealVar *observable = this->getObservable();
3201 auto cache = this->getCache(_curNormSet);
3202 double unc2 = 0;
3203 for (const auto &sample : this->_sampleMap) {
3204 RooAbsArg *phys = this->_physics.at(sample.second);
3205 auto weightName = std::string("w_") + sample.first + "_" + this->GetName();
3206 auto weight = static_cast<RooAbsReal *>(cache->_weights.find(weightName.c_str()));
3207 if (!weight) {
3209 << "unable to find object " + weightName << std::endl;
3210 return 0.0;
3211 }
3212 double newunc2 = 0;
3213 RooHistFunc *hf = dynamic_cast<RooHistFunc *>(phys);
3214 RooRealVar *rv = dynamic_cast<RooRealVar *>(phys);
3215 if (hf) {
3216 const RooDataHist &hist = hf->dataHist();
3217 for (Int_t j = 0; j < observable->getBins(); ++j) {
3218 hist.get(j);
3219 newunc2 += hist.weightSquared();
3220 }
3221 } else if (rv) {
3222 newunc2 = pow(rv->getError(), 2);
3223 }
3224 double w = weight->getVal();
3225 unc2 += newunc2 * w * w;
3226 // std::cout << phys->GetName() << " : " << weight->GetName() << "
3227 // thisweight: " << w << " thisxsec2: " << newunc2 << " weight " << weight
3228 // << std::endl;
3229 }
3230 return sqrt(unc2);
3231}
3232
3233////////////////////////////////////////////////////////////////////////////////
3234/// print the parameters and their current values
3235
3237 // print the parameters and their current values
3238 for (auto obj : this->_operators) {
3239 RooRealVar *param = static_cast<RooRealVar *>(obj);
3240 if (!param)
3241 continue;
3242 param->Print();
3243 }
3244}
3245
3246////////////////////////////////////////////////////////////////////////////////
3247/// print the flags and their current values
3248
3250 for (auto flag : this->_flags) {
3251 RooRealVar *param = static_cast<RooRealVar *>(flag);
3252 if (!param)
3253 continue;
3254 param->Print();
3255 }
3256}
3257
3258////////////////////////////////////////////////////////////////////////////////
3259/// print a set of couplings
3260
3263 for (auto c : couplings) {
3264 std::cout << c.first << ": " << c.second << std::endl;
3265 }
3266}
3267
3268////////////////////////////////////////////////////////////////////////////////
3269/// retrieve the list of bin boundaries
3270
3271std::list<Double_t> *
3273 Double_t xhi) const {
3274 return this->getFunc()->binBoundaries(obs, xlo, xhi);
3275}
3276
3277////////////////////////////////////////////////////////////////////////////////
3278/// retrieve the sample Hint
3279
3280std::list<Double_t> *
3282 Double_t xhi) const {
3283 return this->getFunc()->plotSamplingHint(obs, xlo, xhi);
3284}
3285
3286////////////////////////////////////////////////////////////////////////////////
3287/// call getVal on the internal function
3288
3290 // cout << "XX RooLagrangianMorphFunc::getValV(" << this << ") set = " << set
3291 // << std::endl ;
3292 this->_curNormSet = set;
3293 return RooAbsReal::getValV(set);
3294}
3295
3296////////////////////////////////////////////////////////////////////////////////
3297/// call getVal on the internal function
3298
3300 // call getVal on the internal function
3301 RooRealSumFunc *pdf = this->getFunc();
3302 if (pdf)
3303 return this->_scale * pdf->getVal(_curNormSet);
3304 else
3305 std::cerr << "unable to acquire in-built function!" << std::endl;
3306 return 0.;
3307}
3308
3309////////////////////////////////////////////////////////////////////////////////
3310/// check if this PDF is a binned distribution in the given observable
3311
3312Bool_t
3314 return this->getFunc()->isBinnedDistribution(obs);
3315}
3316
3317////////////////////////////////////////////////////////////////////////////////
3318/// check if observable exists in the RooArgSet (-?-)
3319
3321 return this->getFunc()->checkObservables(nset);
3322}
3323
3324////////////////////////////////////////////////////////////////////////////////
3325/// Force analytical integration for the given observable
3326
3328 return this->getFunc()->forceAnalyticalInt(arg);
3329}
3330
3331////////////////////////////////////////////////////////////////////////////////
3332/// Retrieve the mat
3333
3335 RooArgSet &allVars, RooArgSet &numVars, const RooArgSet *normSet,
3336 const char *rangeName) const {
3337 return this->getFunc()->getAnalyticalIntegralWN(allVars, numVars, normSet,
3338 rangeName);
3339}
3340
3341////////////////////////////////////////////////////////////////////////////////
3342/// Retrieve the matrix of coefficients
3343
3345 Int_t code, const RooArgSet *normSet, const char *rangeName) const {
3346 return this->getFunc()->analyticalIntegralWN(code, normSet, rangeName);
3347}
3348
3349////////////////////////////////////////////////////////////////////////////////
3350/// Retrieve the matrix of coefficients
3351
3352void RooLagrangianMorphFunc::printMetaArgs(std::ostream &os) const {
3353 return this->getFunc()->printMetaArgs(os);
3354}
3355
3356////////////////////////////////////////////////////////////////////////////////
3357/// Retrieve the matrix of coefficients
3358
3360 return this->getFunc()->canNodeBeCached();
3361}
3362
3363////////////////////////////////////////////////////////////////////////////////
3364/// Retrieve the matrix of coefficients
3365
3367 this->getFunc()->setCacheAndTrackHints(arg);
3368}
3369
3370////////////////////////////////////////////////////////////////////////////////
3371/// Retrieve the matrix of coefficients
3372
3374 auto cache = getCache(_curNormSet);
3375 if (!cache)
3376 coutE(Caching) << "unable to retrieve cache!" << std::endl;
3377 return makeRootMatrix(cache->_matrix);
3378}
3379
3380////////////////////////////////////////////////////////////////////////////////
3381/// Retrieve the matrix of coefficients after inversion
3382
3384 auto cache = getCache(_curNormSet);
3385 if (!cache)
3386 coutE(Caching) << "unable to retrieve cache!" << std::endl;
3387 return makeRootMatrix(cache->_inverse);
3388}
3389
3390////////////////////////////////////////////////////////////////////////////////
3391/// Retrieve the condition of the coefficient matrix. If the condition number
3392/// is very large, then the matrix is ill-conditioned and is almost singular.
3393/// The computation of the inverse is prone to large numerical errors
3394
3396 auto cache = getCache(_curNormSet);
3397 if (!cache)
3398 coutE(Caching) << "unable to retrieve cache!" << std::endl;
3399 return cache->_condition;
3400}
3401
3402////////////////////////////////////////////////////////////////////////////////
3403/// Return the RooRatio form of products and denominators of morphing functions
3404
3405std::unique_ptr<RooRatio> RooLagrangianMorphFunc::makeRatio(const char *name,
3406 const char *title,
3407 RooArgList &nr,
3408 RooArgList &dr) {
3409 RooArgList num, denom;
3410 for (auto it : nr) {
3411 num.add(*it);
3412 }
3413 for (auto it : dr) {
3414 denom.add(*it);
3415 }
3416 // same for denom
3417 return make_unique<RooRatio>(name, title, num, denom);
3418}
int Int_t
Definition: CPyCppyy.h:43
void Class()
Definition: Class.C:29
double
Definition: Converters.cxx:939
double SuperFloat
Definition: Floats.h:22
ROOT::R::TRInterface & r
Definition: Object.C:4
#define d(i)
Definition: RSha256.hxx:102
#define b(i)
Definition: RSha256.hxx:100
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
size_t size< TMatrixD >(const TMatrixD &mat)
void writeMatrixToStreamT(const MatrixT &matrix, std::ostream &stream)
write a matrix to a stream
TMatrixD Matrix
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
Matrix makeSuperMatrix(const TMatrixD &in)
convert a TMatrixD into a Matrix
void writeMatrixToFileT(const MatrixT &matrix, const char *fname)
write a matrix to a text file
double invertMatrix(const Matrix &matrix, Matrix &inverse)
#define NaN
constexpr static double morphUnityDeviation
TMatrixD makeRootMatrix(const Matrix &in)
convert a matrix into a TMatrixD
constexpr static double morphLargestWeight
Matrix diagMatrix(size_t n)
create a new diagonal matrix of size n
void printMatrix(const TMatrixD &mat)
write a matrix
#define oocxcoutW(o, a)
Definition: RooMsgService.h:95
#define oocxcoutP(o, a)
Definition: RooMsgService.h:91
#define coutE(a)
Definition: RooMsgService.h:33
#define cxcoutP(a)
Definition: RooMsgService.h:89
#define TRACE_DESTROY
Definition: RooTrace.h:24
#define TRACE_CREATE
Definition: RooTrace.h:23
const Bool_t kFALSE
Definition: RtypesCore.h:101
bool Bool_t
Definition: RtypesCore.h:63
double Double_t
Definition: RtypesCore.h:59
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:364
#define gDirectory
Definition: TDirectory.h:327
int type
Definition: TGX11.cxx:121
int isnan(double)
double pow(double, double)
double sqrt(double)
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:22
std::string & operator+=(std::string &left, const TString &right)
Definition: TString.h:478
char * Form(const char *fmt,...)
TTime operator*(const TTime &t1, const TTime &t2)
Definition: TTime.h:85
RooAbsArg is the common abstract base class for objects that represent a value and a "shape" in RooFi...
Definition: RooAbsArg.h:72
void setStringAttribute(const Text_t *key, const Text_t *value)
Associate string 'value' to this object under key 'key'.
Definition: RooAbsArg.cxx:323
const RefCountList_t & servers() const
List of all servers of this object.
Definition: RooAbsArg.h:199
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:323
friend class RooListProxy
Definition: RooAbsArg.h:621
void setValueDirty()
Mark the element dirty. This forces a re-evaluation when a value is requested.
Definition: RooAbsArg.h:491
Bool_t getAttribute(const Text_t *name) const
Check if a named attribute is set. By default, all attributes are unset.
Definition: RooAbsArg.cxx:314
void addServerList(RooAbsCollection &serverList, Bool_t valueProp=kTRUE, Bool_t shapeProp=kFALSE)
Register a list of RooAbsArg as servers to us by calling addServer() for each arg in the list.
Definition: RooAbsArg.cxx:424
Bool_t isConstant() const
Check if the "Constant" attribute is set.
Definition: RooAbsArg.h:363
virtual Double_t * array() const =0
RooAbsCacheElement is the abstract base class for objects to be stored in RooAbsCache cache manager o...
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects.
RooAbsArg * first() const
Int_t getSize() const
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsData is the common abstract base class for binned and unbinned datasets.
Definition: RooAbsData.h:67
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
virtual Double_t getMax(const char *name=0) const
Get maximum of currently defined range.
virtual void setBin(Int_t ibin, const char *rangeName=0)
Set value to center of bin 'ibin' of binning 'rangeName' (or of default binning if no range is specif...
virtual Int_t numBins(const char *rangeName=0) const
virtual Int_t getBins(const char *name=0) const
Get number of bins of currently defined range.
void setConstant(Bool_t value=kTRUE)
virtual Double_t getMin(const char *name=0) const
Get miniminum of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:61
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
virtual Double_t getValV(const RooArgSet *normalisationSet=nullptr) const
Return value of object.
Definition: RooAbsReal.cxx:274
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition: RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition: RooArgSet.h:35
Class RooBinning is an implements RooAbsBinning in terms of an array of boundary values,...
Definition: RooBinning.h:28
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=0, const TNamed *isetRangeName=0)
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=0)
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:39
double weight(std::size_t i) const
Return weight of i-th bin.
Definition: RooDataHist.h:106
const RooArgSet * get() const override
Get bin centre of current bin.
Definition: RooDataHist.h:78
double weightSquared(std::size_t i) const
Return squared weight sum of i-th bin.
Definition: RooDataHist.h:110
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
Definition: RooFormulaVar.h:30
RooHistFunc implements a real-valued function sampled from a multidimensional histogram.
Definition: RooHistFunc.h:30
RooDataHist & dataHist()
Return RooDataHist that is represented.
Definition: RooHistFunc.h:40
RooArgSet const & getHistObsList() const
Definition: RooHistFunc.h:85
virtual RooArgList containedArgs(Action) override
retrieve the list of contained args
void buildMatrix(const RooLagrangianMorphFunc::ParamMap &inputParameters, const RooLagrangianMorphFunc::FlagMap &inputFlags, const List &flags)
build and invert the morphing matrix
std::unique_ptr< RooRealSumFunc > _sumFunc
void createComponents(const RooLagrangianMorphFunc::ParamMap &inputParameters, const RooLagrangianMorphFunc::FlagMap &inputFlags, const char *funcname, const std::vector< std::vector< RooListProxy * >> &diagrams, const std::vector< RooListProxy * > &nonInterfering, const RooListProxy &flags)
create the basic objects required for the morphing
virtual void operModeHook(RooAbsArg::OperMode) override
Interface for changes of operation mode.
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...
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
static RooLagrangianMorphFunc::CacheElem * createCache(const RooLagrangianMorphFunc *func)
create all the temporary objects required by the class
RooArgList const & getCouplings() const
void setCouplings(const RooAbsCollection &couplings)
config setter for couplings
FlagMap const & getFlagValues() const
std::vector< std::vector< RooArgList * > > const & getDiagrams() const
void append(ParamSet &set, const char *str, double val)
append the parameter map with a parmeter set
ParamMap const & getParamCards() const
void setFileName(const char *filename)
config setter for input file name
std::vector< std::string > const & getFolderNames() const
std::string const & getFileName() const
void readParameters(TDirectory *f)
read the parameters from the input file
RooArgList const & getDecCouplings() const
void setFolders(const RooArgList &folderlist)
config setter for input folders
void allowNegativeYields(Bool_t allowNegativeYields)
config setter for name of the the observable to be morphed
RooArgList const & getFolders() const
void setVertices(const std::vector< T > &vertices)
config setter for vertices
std::string const & getObservableName() const
void setDiagrams(const std::vector< std::vector< T >> &diagrams)
config setter for diagrams
void addFolders(const RooArgList &folders)
convert the RooArgList folders into a simple vector of std::string
void setObservableName(const char *obsname)
config setter for name of the the observable to be morphed
RooArgList const & getProdCouplings() const
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
TPair * makeCrosssectionContainer(double xs, double unc)
create TPair containers of the type expected by the RooLagrangianMorph
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
ParamSet getMorphParameters() const
retrieve the parameter set
RooProduct * getSumElement(const char *name) const
return the RooProduct that is the element of the RooRealSumPdfi corresponding to the given sample nam...
void insert(RooWorkspace *ws)
insert this object into a workspace
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::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 usage: countSamples ( { Roo...
virtual Double_t analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=0) const override
Retrieve the matrix of coefficients.
virtual void setCacheAndTrackHints(RooArgSet &) override
Retrieve the matrix of coefficients.
virtual Double_t getValV(const RooArgSet *set=0) const override
call getVal on the internal function
ParamSet getCouplings() const
retrieve a set of couplings (-?-)
void printSampleWeights() const
print the current sample weights
std::map< const std::string, double > ParamSet
RooArgSet createWeights(const ParamMap &inputs, const std::vector< RooArgList * > &vertices, RooArgList &couplings, const FlagMap &inputFlags, const RooArgList &flags, const std::vector< RooArgList * > &nonInterfering)
create only the weight formulas. static function for external usage.
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.
virtual Double_t evaluate() const override
call getVal on the internal function
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
RooRealVar * getObservable() const
retrieve the histogram observable
int countContributingFormulas() const
count the number of formulas that correspond to the current parameter set
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
virtual ~RooLagrangianMorphFunc()
default destructor
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)
double getScale()
get energy scale of the EFT expansion
std::vector< RooListProxy * > _nonInterfering
virtual std::list< Double_t > * plotSamplingHint(RooAbsRealLValue &, Double_t, Double_t) const override
retrieve the sample Hint
virtual Bool_t checkObservables(const RooArgSet *nset) const override
check if observable exists in the RooArgSet (-?-)
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
void printParameters() const
print the parameters and their current values
void printAuthors() const
print the authors information
void printPhysics() const
print the current physics values
virtual std::list< Double_t > * binBoundaries(RooAbsRealLValue &, Double_t, Double_t) const override
retrieve the list of bin boundaries
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
RooLagrangianMorphFunc::CacheElem * getCache(const RooArgSet *nset) const
retrieve the cache object
bool hasParameter(const char *paramname) const
check if a parameter of the given name is contained in the list of known parameters
virtual Bool_t isBinnedDistribution(const RooArgSet &obs) const override
check if this PDF is a binned distribution in the given observable
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
Double_t expectedEvents() const
return the number of expected events for the current parameter set
TMatrixD getInvertedMatrix() const
Retrieve the matrix of coefficients after inversion.
virtual 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
bool isParameterUsed(const char *paramname) const
check if there is any morphing power provided for the given parameter morphing power is provided as s...
RooAbsPdf::ExtendMode extendMode() const
get a standalone clone of the func that does not depend on this object
std::map< const std::string, FlagSet > FlagMap
void setParameterConstant(const char *paramname, bool constant) const
call setConstant with the boolean argument provided on the parameter with the given name
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 ...
virtual 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::vector< std::vector< RooListProxy * > > _diagrams
void importToWorkspace(RooWorkspace *ws, const RooAbsReal *object)
insert an object into a workspace (wrapper for RooWorkspace::import)
std::map< std::string, int > _sampleMap
virtual TObject * clone(const char *newname) const override
cloning method
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
virtual Bool_t forceAnalyticalInt(const RooAbsArg &arg) const override
Force analytical integration for the given observable.
RooRealSumFunc * getFunc() const
get the func
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
virtual Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &numVars, const RooArgSet *normSet, const char *rangeName=0) const override
Retrieve the mat.
RooLinearCombination is a class that helps perform linear combination of floating point numbers and p...
void add(SuperFloat c, RooAbsReal *t)
void setCoefficient(size_t idx, SuperFloat c)
SuperFloat getCoefficient(size_t idx)
RooListProxy is the concrete proxy for RooArgList objects.
Definition: RooListProxy.h:24
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Reimplementation of standard RooArgList::add()
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.
Definition: RooProduct.h:29
RooArgList components()
Definition: RooProduct.h:47
virtual Bool_t forceAnalyticalInt(const RooAbsArg &arg) const
Bool_t isBinnedDistribution(const RooArgSet &obs) const
Tests if the distribution is binned. Unless overridden by derived classes, this always returns false.
virtual std::list< Double_t > * binBoundaries(RooAbsRealLValue &, Double_t, Double_t) const
Retrieve bin boundaries if this distribution is binned in obs.
Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &numVars, const RooArgSet *normSet, const char *rangeName=0) const
Variant of getAnalyticalIntegral that is also passed the normalization set that should be applied to ...
virtual std::list< Double_t > * plotSamplingHint(RooAbsRealLValue &, Double_t, Double_t) const
Interface for returning an optional hint for initial sampling points when constructing a curve projec...
virtual CacheMode canNodeBeCached() const
virtual Bool_t checkObservables(const RooArgSet *nset) const
Overloadable function in which derived classes can implement consistency checks of the variables.
Double_t analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=0) const
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
virtual void setCacheAndTrackHints(RooArgSet &)
void printMetaArgs(std::ostream &os) const
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:39
void setMin(const char *name, Double_t value)
Set minimum of name range to given value.
Definition: RooRealVar.cxx:461
void setBins(Int_t nBins, const char *name=0)
Create a uniform binning under name 'name' for this variable.
Definition: RooRealVar.cxx:408
void setError(Double_t value)
Definition: RooRealVar.h:62
void setMax(const char *name, Double_t value)
Set maximum of name range to given value.
Definition: RooRealVar.cxx:491
Double_t getError() const
Definition: RooRealVar.h:60
const RooAbsBinning & getBinning(const char *name=0, Bool_t verbose=kTRUE, Bool_t createOnTheFly=kFALSE) const
Return binning definition with name.
Definition: RooRealVar.cxx:318
void setBinning(const RooAbsBinning &binning, const char *name=0)
Add given binning under name 'name' with this variable.
Definition: RooRealVar.cxx:415
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:257
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE) override
Overloaded RooArgSet::add() method inserts 'var' into set and registers 'var' as server to owner with...
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE) override
Overloaded RooArgSet::addOwned() method insert object into owning set and registers object as server ...
RooStringVar is a RooAbsArg implementing string values.
Definition: RooStringVar.h:23
const char * getVal() const
Definition: RooStringVar.h:34
The RooWorkspace is a persistable container for RooFit projects.
Definition: RooWorkspace.h:43
Class to manage histogram axis.
Definition: TAxis.h:30
Double_t GetXmax() const
Definition: TAxis.h:134
const char * GetBinLabel(Int_t bin) const
Return label for bin.
Definition: TAxis.cxx:440
Double_t GetXmin() const
Definition: TAxis.h:133
Int_t GetNbins() const
Definition: TAxis.h:121
TClass instances represent classes, structs and namespaces in the ROOT type system.
Definition: TClass.h:80
Bool_t InheritsFrom(const char *cl) const
Return kTRUE if this class inherits from a class with name "classname".
Definition: TClass.cxx:4867
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.
Definition: TClass.cxx:2955
Double_t GetCondition() const
Definition: TDecompBase.h:70
LU Decomposition class.
Definition: TDecompLU.h:24
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...
Definition: TDecompLU.cxx:526
Describe directory structure in memory.
Definition: TDirectory.h:45
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:54
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.
Definition: TFile.cxx:4011
A TFolder object is a collection of objects and folders.
Definition: TFolder.h:30
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:575
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:320
virtual Int_t GetNbinsX() const
Definition: TH1.h:296
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...
Definition: TH1.cxx:9040
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...
Definition: TH1.cxx:9056
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition: TH1.cxx:8986
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4993
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition: TH1.cxx:8997
Int_t GetNrows() const
Definition: TMatrixTBase.h:123
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
virtual TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1)
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1211
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
Mother of all ROOT objects.
Definition: TObject.h:37
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:359
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:130
virtual TObject * FindObject(const char *name) const
Must be redefined in derived classes.
Definition: TObject.cxx:323
Class used by TMap to store (key,value) pairs.
Definition: TMap.h:102
TObject * Key() const
Definition: TMap.h:120
TObject * Value() const
Definition: TMap.h:121
Named parameter, streamable and storable.
Definition: TParameter.h:35
const AParamType & GetVal() const
Definition: TParameter.h:67
Random number generator class based on M.
Definition: TRandom3.h:27
Basic string class.
Definition: TString.h:136
int CompareTo(const char *cs, ECaseCompare cmp=kExact) const
Compare a string to char *cs2.
Definition: TString.cxx:442
std::string_view View() const
Definition: TString.h:444
const char * Data() const
Definition: TString.h:369
TString & Append(const char *cs)
Definition: TString.h:564
TLine * line
RooCmdArg RecycleConflictNodes(Bool_t flag=kTRUE)
const Int_t n
Definition: legend1.C:16
#define T2
Definition: md5.inl:146
#define T1
Definition: md5.inl:145
static double B[]
static double A[]
double T(double x)
Definition: ChebyshevPol.h:34
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
static const std::string name("name")
@ InputArguments
Definition: RooGlobalFunc.h:61
@ ObjectHandling
Definition: RooGlobalFunc.h:61
Definition: file.py:1
Definition: first.py:1
const char * Ref
Definition: TXMLSetup.cxx:53
Definition: civetweb.c:2228
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12
REAL * vertex
Definition: triangle.c:513
void ws()
Definition: ws.C:66