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};
95
96template <template <typename...> class Ref, typename... Args>
97struct is_specialization<Ref<Args...>, Ref> : std::true_type {
98};
99
100///////////////////////////////////////////////////////////////////////////////
101// LINEAR ALGEBRA HELPERS /////////////////////////////////////////////////////
102///////////////////////////////////////////////////////////////////////////////
103
104////////////////////////////////////////////////////////////////////////////////
105/// retrieve the size of a square matrix
106
107template <class MatrixT>
108inline size_t size(const MatrixT &matrix);
109template <>
110inline size_t size<TMatrixD>(const TMatrixD &mat)
111{
112 return mat.GetNrows();
113}
114
115////////////////////////////////////////////////////////////////////////////////
116/// write a matrix to a stream
117
118template <class MatrixT>
119void writeMatrixToStreamT(const MatrixT &matrix, std::ostream &stream)
120{
121 if (!stream.good()) {
122 return;
123 }
124 for (size_t i = 0; i < size(matrix); ++i) {
125 for (size_t j = 0; j < size(matrix); ++j) {
126#ifdef USE_UBLAS
127 stream << std::setprecision(RooFit::SuperFloatPrecision::digits10) << matrix(i, j) << "\t";
128#else
129 stream << matrix(i, j) << "\t";
130#endif
131 }
132 stream << std::endl;
133 }
134}
135
136////////////////////////////////////////////////////////////////////////////////
137/// write a matrix to a text file
138
139template <class MatrixT>
140inline void writeMatrixToFileT(const MatrixT &matrix, const char *fname)
141{
142 std::ofstream of(fname);
143 if (!of.good()) {
144 cerr << "unable to read file '" << fname << "'!" << std::endl;
145 }
146 writeMatrixToStreamT(matrix, of);
147 of.close();
148}
149
150#ifdef USE_UBLAS
151
152// boost includes
153#pragma GCC diagnostic push
154#pragma GCC diagnostic ignored "-Wshadow"
155#pragma GCC diagnostic ignored "-Wunused-local-typedefs"
156#include <boost/numeric/ublas/io.hpp>
157#include <boost/numeric/ublas/lu.hpp>
158#include <boost/numeric/ublas/matrix.hpp>
159#include <boost/numeric/ublas/matrix_expression.hpp>
160#include <boost/numeric/ublas/symmetric.hpp> //inc diag
161#include <boost/numeric/ublas/triangular.hpp>
162#include <boost/operators.hpp>
163
164#pragma GCC diagnostic pop
165
166typedef boost::numeric::ublas::matrix<RooFit::SuperFloat> Matrix;
167
168////////////////////////////////////////////////////////////////////////////////
169/// write a matrix
170
171inline void printMatrix(const Matrix &mat)
172{
173 for (size_t i = 0; i < mat.size1(); ++i) {
174 for (size_t j = 0; j < mat.size2(); ++j) {
175 std::cout << std::setprecision(RooFit::SuperFloatPrecision::digits10) << mat(i, j) << " ,\t";
176 }
177 std::cout << std::endl;
178 }
179}
180
181////////////////////////////////////////////////////////////////////////////////
182/// retrieve the size of a square matrix
183
184template <>
185inline size_t size<Matrix>(const Matrix &matrix)
186{
187 return matrix.size1();
188}
189
190////////////////////////////////////////////////////////////////////////////////
191/// create a new diagonal matrix of size n
192
193inline Matrix diagMatrix(size_t n)
194{
195 return boost::numeric::ublas::identity_matrix<RooFit::SuperFloat>(n);
196}
197
198////////////////////////////////////////////////////////////////////////////////
199/// convert a matrix into a TMatrixD
200
201inline TMatrixD makeRootMatrix(const Matrix &in)
202{
203 size_t n = size(in);
204 TMatrixD mat(n, n);
205 for (size_t i = 0; i < n; ++i) {
206 for (size_t j = 0; j < n; ++j) {
207 mat(i, j) = double(in(i, j));
208 }
209 }
210 return mat;
211}
212
213////////////////////////////////////////////////////////////////////////////////
214/// convert a TMatrixD into a matrix
215
216inline Matrix makeSuperMatrix(const TMatrixD &in)
217{
218 size_t n = in.GetNrows();
219 Matrix mat(n, n);
220 for (size_t i = 0; i < n; ++i) {
221 for (size_t j = 0; j < n; ++j) {
222 mat(i, j) = double(in(i, j));
223 }
224 }
225 return mat;
226}
227
228inline Matrix operator+=(const Matrix &rhs)
229{
230 return add(rhs);
231}
232inline Matrix operator*(const Matrix &m, const Matrix &otherM)
233{
234 return prod(m, otherM);
235}
236
237////////////////////////////////////////////////////////////////////////////////
238/// calculate the inverse of a matrix, returning the condition
239
240inline RooFit::SuperFloat invertMatrix(const Matrix &matrix, Matrix &inverse)
241{
242 boost::numeric::ublas::permutation_matrix<size_t> pm(size(matrix));
243 RooFit::SuperFloat mnorm = norm_inf(matrix);
244 Matrix lu(matrix);
245 try {
246 int res = lu_factorize(lu, pm);
247 if (res != 0) {
248 std::stringstream ss;
249 ::writeMatrixToStreamT(matrix, ss);
250 cxcoutP(Eval) << ss.str << std::endl;
251 }
252 // back-substitute to get the inverse
253 lu_substitute(lu, pm, inverse);
254 } catch (boost::numeric::ublas::internal_logic &error) {
255 // coutE(Eval) << "boost::numberic::ublas error: matrix is not invertible!"
256 // << std::endl;
257 }
258 RooFit::SuperFloat inorm = norm_inf(inverse);
259 RooFit::SuperFloat condition = mnorm * inorm;
260 return condition;
261}
262
263#else
264
265#include "TDecompLU.h"
267
268////////////////////////////////////////////////////////////////////////////////
269/// convert a matrix into a TMatrixD
270
272{
273 return TMatrixD(in);
274}
275
276////////////////////////////////////////////////////////////////////////////////
277/// convert a TMatrixD into a Matrix
278
280{
281 return in;
282}
283
284////////////////////////////////////////////////////////////////////////////////
285/// create a new diagonal matrix of size n
286
287inline Matrix diagMatrix(size_t n)
288{
289 TMatrixD mat(n, n);
290 mat.UnitMatrix();
291 return mat;
292}
293
294////////////////////////////////////////////////////////////////////////////////
295/// write a matrix
296
297inline void printMatrix(const TMatrixD &mat)
298{
299 writeMatrixToStreamT(mat, std::cout);
300}
301
302////////////////////////////////////////////////////////////////////////////////
303// calculate the inverse of a matrix, returning the condition
304
305inline double invertMatrix(const Matrix &matrix, Matrix &inverse)
306{
307 TDecompLU lu(matrix);
308 bool status = lu.Invert(inverse);
309 // check if the matrix is invertible
310 if (!status) {
311 std::cerr << " matrix is not invertible!" << std::endl;
312 }
313 double condition = lu.GetCondition();
314 const size_t n = size(inverse);
315 // sanitize numeric problems
316 for (size_t i = 0; i < n; ++i)
317 for (size_t j = 0; j < n; ++j)
318 if (fabs(inverse(i, j)) < 1e-9)
319 inverse(i, j) = 0;
320 return condition;
321}
322#endif
323
324/////////////////////////////////////////////////////////////////////////////////
325// LOCAL FUNCTIONS AND DEFINITIONS
326// //////////////////////////////////////////////
327/////////////////////////////////////////////////////////////////////////////////
328/// anonymous namespace to prohibit use of these functions outside the class
329/// itself
330namespace {
331///////////////////////////////////////////////////////////////////////////////
332// HELPERS ////////////////////////////////////////////////////////////////////
333///////////////////////////////////////////////////////////////////////////////
334
335typedef std::vector<std::vector<bool>> FeynmanDiagram;
336typedef std::vector<std::vector<int>> MorphFuncPattern;
337typedef std::map<int, std::unique_ptr<RooAbsReal>> FormulaList;
338
339///////////////////////////////////////////////////////////////////////////////
340/// (-?-)
341
342inline TString makeValidName(const char *input)
343{
344 TString retval(input);
345 retval.ReplaceAll("/", "_");
346 retval.ReplaceAll("^", "");
347 retval.ReplaceAll("*", "X");
348 retval.ReplaceAll("[", "");
349 retval.ReplaceAll("]", "");
350 return retval;
351}
352
353//////////////////////////////////////////////////////////////////////////////
354/// concatenate the names of objects in a collection to a single string
355
356template <class List>
357std::string concatNames(const List &c, const char *sep)
358{
359 std::stringstream ss;
360 bool first = true;
361 for (auto itr : c) {
362 if (!first)
363 ss << sep;
364 ss << itr->GetName();
365 first = false;
366 }
367 return ss.str();
368}
369
370///////////////////////////////////////////////////////////////////////////////
371/// this is a workaround for the missing implicit conversion from
372/// SuperFloat<>double
373
374template <class A, class B>
375inline void assignElement(A &a, const B &b)
376{
377 a = static_cast<A>(b);
378}
379///////////////////////////////////////////////////////////////////////////////
380// read a matrix from a stream
381
382template <class MatrixT>
383inline MatrixT readMatrixFromStreamT(std::istream &stream)
384{
385 std::vector<std::vector<RooFit::SuperFloat>> matrix;
386 std::vector<RooFit::SuperFloat> line;
387 while (!stream.eof()) {
388 if (stream.peek() == '\n') {
389 stream.get();
390 stream.peek();
391 continue;
392 }
394 stream >> val;
395 line.push_back(val);
396 while (stream.peek() == ' ' || stream.peek() == '\t') {
397 stream.get();
398 }
399 if (stream.peek() == '\n') {
400 matrix.push_back(line);
401 line.clear();
402 }
403 }
404 MatrixT retval(matrix.size(), matrix.size());
405 for (size_t i = 0; i < matrix.size(); ++i) {
406 if (matrix[i].size() != matrix.size()) {
407 std::cerr << "matrix read from stream doesn't seem to be square!" << std::endl;
408 }
409 for (size_t j = 0; j < matrix[i].size(); ++j) {
410 assignElement(retval(i, j), matrix[i][j]);
411 }
412 }
413 return retval;
414}
415
416///////////////////////////////////////////////////////////////////////////////
417/// read a matrix from a text file
418
419template <class MatrixT>
420inline MatrixT readMatrixFromFileT(const char *fname)
421{
422 std::ifstream in(fname);
423 if (!in.good()) {
424 std::cerr << "unable to read file '" << fname << "'!" << std::endl;
425 }
426 MatrixT mat = readMatrixFromStreamT<MatrixT>(in);
427 in.close();
428 return mat;
429}
430
431///////////////////////////////////////////////////////////////////////////////
432/// convert a TH1* param hist into the corresponding ParamSet object
433
434template <class T>
435void readValues(std::map<const std::string, T> &myMap, TH1 *h_pc)
436{
437 if (h_pc) {
438 // loop over all bins of the param_card histogram
439 for (int ibx = 1; ibx <= h_pc->GetNbinsX(); ++ibx) {
440 // read the value of one parameter
441 const std::string s_coup(h_pc->GetXaxis()->GetBinLabel(ibx));
442 double coup_val = h_pc->GetBinContent(ibx);
443 // add it to the map
444 if (!s_coup.empty()) {
445 myMap[s_coup] = T(coup_val);
446 }
447 }
448 }
449}
450
451///////////////////////////////////////////////////////////////////////////////
452/// Set up folder ownership over its children, and treat likewise any subfolders.
453/// @param theFolder: folder to update. Assumed to be a valid pointer
454void setOwnerRecursive(TFolder* theFolder){
455 theFolder->SetOwner();
456 // And also need to set up ownership for nested folders
457 auto subdirs = theFolder->GetListOfFolders();
458 for (auto* subdir : *subdirs){
459 auto thisfolder = dynamic_cast<TFolder*>(subdir);
460 if (thisfolder){
461 // no explicit deletion here, will be handled by parent
462 setOwnerRecursive(thisfolder);
463 }
464 }
465}
466
467
468///////////////////////////////////////////////////////////////////////////////
469/// Load a TFolder from a file while ensuring it owns its content.
470/// This avoids memory leaks. Note that when fetching objects
471/// from this folder, you need to clone them to prevent deletion.
472/// Also recursively updates nested subfolders accordingly
473/// @param inFile: Input file to read - assumed to be a valid pointer
474/// @param folderName: Name of the folder to read from the file
475/// @return a unique_ptr to the folder. Nullptr if not found.
476std::unique_ptr<TFolder> readOwningFolderFromFile(TDirectory* inFile, const std::string & folderName){
477 std::unique_ptr<TFolder> theFolder(inFile->Get<TFolder>(folderName.c_str()));
478 if (!theFolder) {
479 std::cerr << "Error: unable to access data from folder '" << folderName << "' from file '"<<inFile->GetName()<<"'!" << std::endl;
480 return nullptr;
481 }
482 setOwnerRecursive(theFolder.get());
483 return theFolder;
484}
485
486///////////////////////////////////////////////////////////////////////////////
487/// Helper to load a single object from a file-resident TFolder, while
488/// avoiding memory leaks.
489/// @tparam AObjType Type of object to load.
490/// @param inFile input file to load from. Expected to be a valid pointer
491/// @param folderName Name of the TFolder to load from the file
492/// @param objName Name of the object to load
493/// @param notFoundError If set, print a detailed error if we didn't find something
494/// @return Returns a pointer to a clone of the loaded object. Ownership assigned to the caller.
495template <class AObjType> AObjType* loadFromFileResidentFolder(TDirectory* inFile,
496 const std::string & folderName,
497 const std::string & objName,
498 bool notFoundError = true)
499{
500 auto folder = readOwningFolderFromFile(inFile, folderName);
501 if (!folder) {
502 return nullptr;
503 }
504 AObjType *loadedObject = dynamic_cast<AObjType *>(folder->FindObject(objName.c_str()));
505 if (!loadedObject) {
506 if (notFoundError){
507 std::stringstream errstr;
508 errstr << "Error: unable to retrieve object '" << objName << "' from folder '" << folderName
509 << "'. contents are:";
510 TIter next(folder->GetListOfFolders()->begin());
511 TFolder *f;
512 while ((f = (TFolder *)next())) {
513 errstr << " " << f->GetName();
514 }
515 std::cerr << errstr.str() << std::endl;
516 }
517 return nullptr;
518 }
519 // replace the loaded object by a clone, as the loaded folder will delete the original
520 return static_cast<AObjType *>(loadedObject->Clone()); // can use a static_cast - confirmed validity by initial cast above.
521}
522
523
524///////////////////////////////////////////////////////////////////////////////
525/// retrieve a ParamSet from a certain subfolder 'name' of the file
526
527template <class T>
528void readValues(std::map<const std::string, T> &myMap, TDirectory *file, const std::string &name,
529 const std::string &key = "param_card", bool notFoundError = true)
530{
531 TH1F *h_pc = loadFromFileResidentFolder<TH1F>(file, name, key, notFoundError);
532 readValues(myMap, h_pc);
533}
534
535///////////////////////////////////////////////////////////////////////////////
536/// retrieve the param_hists file and return a map of the parameter values
537/// by providing a list of names, only the param_hists of those subfolders are
538/// read leaving the list empty is interpreted as meaning 'read everything'
539
540template <class T>
541void readValues(std::map<const std::string, std::map<const std::string, T>> &inputParameters, TDirectory *f,
542 const std::vector<std::string> &names, const std::string &key = "param_card", bool notFoundError = true)
543{
544 inputParameters.clear();
545 // if the list of names is empty, we assume that this means 'all'
546 // loop over all folders in the file
547 for (size_t i = 0; i < names.size(); i++) {
548 const std::string name(names[i]);
549 // actually read an individual param_hist
550 readValues(inputParameters[name], f, name, key, notFoundError);
551 }
552
553 // now the map is filled with all parameter values found for all samples
554}
555
556///////////////////////////////////////////////////////////////////////////////
557/// open the file and return a file pointer
558
559inline TDirectory *openFile(const std::string &filename)
560{
561 if (filename.empty()) {
562 return gDirectory;
563 } else {
564 TFile *file = TFile::Open(filename.c_str(), "READ");
565 if (!file || !file->IsOpen()) {
566 if (file)
567 delete file;
568 std::cerr << "could not open file '" << filename << "'!" << std::endl;
569 }
570 return file;
571 }
572}
573
574///////////////////////////////////////////////////////////////////////////////
575/// open the file and return a file pointer
576
577inline void closeFile(TDirectory *d)
578{
579 TFile *f = dynamic_cast<TFile *>(d);
580 if (f) {
581 f->Close();
582 delete f;
583 d = nullptr;
584 }
585}
586
587///////////////////////////////////////////////////////////////////////////////
588/// extract the operators from a single coupling
589template <class T2>
590inline void extractServers(const RooAbsArg &coupling, T2 &operators)
591{
592 int nservers = 0;
593 for (const auto server : coupling.servers()) {
594 extractServers(*server, operators);
595 nservers++;
596 }
597 if (nservers == 0) {
598 operators.add(coupling);
599 }
600}
601
602///////////////////////////////////////////////////////////////////////////////
603/// extract the operators from a list of couplings
604
605template <class T1, class T2, typename std::enable_if<!is_specialization<T1, std::vector>::value, T1>::type * = nullptr>
606inline void extractOperators(const T1 &couplings, T2 &operators)
607{
608 // coutD(InputArguments) << "extracting operators from
609 // "<<couplings.getSize()<<" couplings" << std::endl;
610 for (auto itr : couplings) {
611 extractServers(*itr, operators);
612 }
613}
614
615///////////////////////////////////////////////////////////////////////////////
616/// extract the operators from a list of vertices
617
618template <class T1, class T2, typename std::enable_if<is_specialization<T1, std::vector>::value, T1>::type * = nullptr>
619inline void extractOperators(const T1 &vec, T2 &operators)
620{
621 for (const auto &v : vec) {
622 extractOperators(v, operators);
623 }
624}
625
626///////////////////////////////////////////////////////////////////////////////
627/// extract the couplings from a given set and copy them to a new one
628
629template <class T1, class T2>
630inline void extractCouplings(const T1 &inCouplings, T2 &outCouplings)
631{
632 for (auto itr : inCouplings) {
633 if (!outCouplings.find(itr->GetName())) {
634 // coutD(InputArguments) << "adding parameter "<< obj->GetName() <<
635 // std::endl;
636 outCouplings.add(*itr);
637 }
638 }
639}
640
641///////////////////////////////////////////////////////////////////////////////
642/// find and, if necessary, create a parameter from a list
643
644template <class T>
645inline RooAbsArg &get(T &operators, const char *name, double defaultval = 0)
646{
647 RooAbsArg *kappa = operators.find(name);
648 if (kappa)
649 return *kappa;
650 RooRealVar *newKappa = new RooRealVar(name, name, defaultval);
651 newKappa->setConstant(false);
652 operators.add(*newKappa);
653 return *newKappa;
654}
655
656///////////////////////////////////////////////////////////////////////////////
657/// find and, if necessary, create a parameter from a list
658
659template <class T>
660inline RooAbsArg &get(T &operators, const std::string &name, double defaultval = 0)
661{
662 return get(operators, name.c_str(), defaultval);
663}
664
665////////////////////////////////////////////////////////////////////////////////
666/// create a new coupling and add it to the set
667
668template <class T>
669inline void addCoupling(T &set, const TString &name, const TString &formula, const RooArgList &components, bool isNP)
670{
671 if (!set.find(name)) {
672 RooFormulaVar *c = new RooFormulaVar(name, formula, components);
673 c->setAttribute("NewPhysics", isNP);
674 set.add(*c);
675 }
676}
677
678////////////////////////////////////////////////////////////////////////////////
679/// set parameter values first set all values to defaultVal (if value not
680/// present in param_card then it should be 0)
681
682inline bool setParam(RooRealVar *p, double val, bool force)
683{
684 bool ok = true;
685 if (val > p->getMax()) {
686 if (force) {
687 p->setMax(val);
688 } else {
689 std::cerr << ": parameter " << p->GetName() << " out of bounds: " << val << " > " << p->getMax() << std::endl;
690 ok = false;
691 }
692 } else if (val < p->getMin()) {
693 if (force) {
694 p->setMin(val);
695 } else {
696 std::cerr << ": parameter " << p->GetName() << " out of bounds: " << val << " < " << p->getMin() << std::endl;
697 ok = false;
698 }
699 }
700 if (ok)
701 p->setVal(val);
702 return ok;
703}
704
705////////////////////////////////////////////////////////////////////////////////
706/// set parameter values first set all values to defaultVal (if value not
707/// present in param_card then it should be 0)
708
709template <class T1, class T2>
710inline bool setParams(const T2 &args, T1 val)
711{
712 for (auto itr : args) {
713 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
714 if (!param)
715 continue;
716 setParam(param, val, true);
717 }
718 return true;
719}
720
721////////////////////////////////////////////////////////////////////////////////
722/// set parameter values first set all values to defaultVal (if value not
723/// present in param_card then it should be 0)
724
725template <class T1, class T2>
726inline bool
727setParams(const std::map<const std::string, T1> &point, const T2 &args, bool force = false, T1 defaultVal = 0)
728{
729 bool ok = true;
730 for (auto itr : args) {
731 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
732 if (!param || param->isConstant())
733 continue;
734 ok = setParam(param, defaultVal, force) && ok;
735 }
736 // set all parameters to the values in the param_card histogram
737 for (auto paramit : point) {
738 // loop over all the parameters
739 const std::string param(paramit.first);
740 // retrieve them from the map
741 RooRealVar *p = dynamic_cast<RooRealVar *>(args.find(param.c_str()));
742 if (!p)
743 continue;
744 // set them to their nominal value
745 ok = setParam(p, paramit.second, force) && ok;
746 }
747 return ok;
748}
749
750////////////////////////////////////////////////////////////////////////////////
751/// set parameter values first set all values to defaultVal (if value not
752/// present in param_card then it should be 0)
753
754template <class T>
755inline bool setParams(TH1 *hist, const T &args, bool force = false)
756{
757 bool ok = true;
758
759 for (auto itr : args) {
760 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
761 if (!param)
762 continue;
763 ok = setParam(param, 0., force) && ok;
764 }
765
766 // set all parameters to the values in the param_card histogram
767 TAxis *ax = hist->GetXaxis();
768 for (int i = 1; i <= ax->GetNbins(); ++i) {
769 // loop over all the parameters
770 RooRealVar *p = dynamic_cast<RooRealVar *>(args.find(ax->GetBinLabel(i)));
771 if (!p)
772 continue;
773 // set them to their nominal value
774 ok = setParam(p, hist->GetBinContent(i), force) && ok;
775 }
776 return ok;
777}
778
779////////////////////////////////////////////////////////////////////////////////
780/// create a set of parameters
781
782template <class T>
783inline RooLagrangianMorphFunc::ParamSet getParams(const T &parameters)
784{
786 for (auto itr : parameters) {
787 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
788 if (!param)
789 continue;
790 retval[param->GetName()] = param->getVal();
791 }
792 return retval;
793}
794
795////////////////////////////////////////////////////////////////////////////////
796/// collect the histograms from the input file and convert them to RooFit
797/// objects
798
799void collectHistograms(const char *name, TDirectory *file, std::map<std::string, int> &list_hf, RooArgList &physics,
800 RooRealVar &var, const std::string &varname,
801 const RooLagrangianMorphFunc::ParamMap &inputParameters)
802{
803 bool binningOK = false;
804 for (auto sampleit : inputParameters) {
805 const std::string sample(sampleit.first);
806 TH1* hist = loadFromFileResidentFolder<TH1>(file, sample, varname, true);
807 if (!hist) return;
808
809 auto it = list_hf.find(sample);
810 if (it != list_hf.end()) {
811 RooHistFunc *hf = (RooHistFunc *)(physics.at(it->second));
812 hf->setValueDirty();
813 // commenting out To-be-resolved
814 // RooDataHist* dh = &(hf->dataHist());
815 // RooLagrangianMorphFunc::setDataHistogram(hist,&var,dh);
816 // RooArgSet vars;
817 // vars.add(var);
818 // dh->importTH1(vars,*hist,1.,kFALSE);
819 } else {
820 if (!binningOK) {
821 int n = hist->GetNbinsX();
822 std::vector<double> bins;
823 for (int i = 1; i < n + 1; ++i) {
824 bins.push_back(hist->GetBinLowEdge(i));
825 }
826 bins.push_back(hist->GetBinLowEdge(n) + hist->GetBinWidth(n));
827 var.setBinning(RooBinning(n, &(bins[0])));
828 }
829
830 // generate the mean value
831 TString histname = makeValidName(Form("dh_%s_%s", sample.c_str(), name));
832 TString funcname = makeValidName(Form("phys_%s_%s", sample.c_str(), name));
833 RooArgSet vars;
834 vars.add(var);
835
836 RooDataHist *dh = new RooDataHist(histname.View(), histname.View(), vars, hist);
837 // add it to the list
838 RooHistFunc *hf = new RooHistFunc(funcname, funcname, var, *dh);
839 int idx = physics.getSize();
840 list_hf[sample] = idx;
841 physics.add(*hf);
842 assert(hf = static_cast<RooHistFunc *>(physics.at(idx)));
843 }
844 // std::cout << "found histogram " << hist->GetName() << " with integral "
845 // << hist->Integral() << std::endl;
846 }
847}
848
849////////////////////////////////////////////////////////////////////////////////
850/// collect the RooAbsReal objects from the input directory
851
852void collectRooAbsReal(const char * /*name*/, TDirectory *file, std::map<std::string, int> &list_hf,
853 RooArgList &physics, const std::string &varname,
854 const RooLagrangianMorphFunc::ParamMap &inputParameters)
855{
856 for (auto sampleit : inputParameters) {
857 const std::string sample(sampleit.first);
858 RooAbsReal* obj = loadFromFileResidentFolder<RooAbsReal>(file,sample,varname,true);
859 if (!obj) return;
860 auto it = list_hf.find(sample);
861 if (it == list_hf.end()) {
862 int idx = physics.getSize();
863 list_hf[sample] = idx;
864 physics.add(*obj);
865 assert(obj == physics.at(idx));
866 }
867 }
868}
869
870////////////////////////////////////////////////////////////////////////////////
871/// collect the TParameter objects from the input file and convert them to
872/// RooFit objects
873
874template <class T>
875void collectCrosssections(const char *name, TDirectory *file, std::map<std::string, int> &list_xs, RooArgList &physics,
876 const std::string &varname, const RooLagrangianMorphFunc::ParamMap &inputParameters)
877{
878 for (auto sampleit : inputParameters) {
879 const std::string sample(sampleit.first);
880 TObject* obj = loadFromFileResidentFolder<TObject>(file, sample, varname,false);
881 TParameter<T> *xsection = nullptr;
882 TParameter<T> *error = nullptr;
883 TParameter<T> *p = dynamic_cast<TParameter<T> *>(obj);
884 if (p) {
885 xsection = p;
886 }
887 TPair *pair = dynamic_cast<TPair *>(obj);
888 if (pair) {
889 xsection = dynamic_cast<TParameter<T> *>(pair->Key());
890 error = dynamic_cast<TParameter<T> *>(pair->Value());
891 }
892 if (!xsection) {
893 std::stringstream errstr;
894 errstr << "Error: unable to retrieve cross section '" << varname << "' from folder '" << sample;
895 return;
896 }
897
898 auto it = list_xs.find(sample.c_str());
899 RooRealVar *xs;
900 if (it != list_xs.end()) {
901 xs = (RooRealVar *)(physics.at(it->second));
902 xs->setVal(xsection->GetVal());
903 } else {
904 std::string objname = Form("phys_%s_%s", name, sample.c_str());
905 xs = new RooRealVar(objname.c_str(), objname.c_str(), xsection->GetVal());
906 xs->setConstant(true);
907 int idx = physics.getSize();
908 list_xs[sample] = idx;
909 physics.add(*xs);
910 assert(physics.at(idx) == xs);
911 }
912 if (error)
913 xs->setError(error->GetVal());
914 }
915}
916
917////////////////////////////////////////////////////////////////////////////////
918/// collect the TPair<TParameter,TParameter> objects from the input file and
919/// convert them to RooFit objects
920
921void collectCrosssectionsTPair(const char *name, TDirectory *file, std::map<std::string, int> &list_xs,
922 RooArgList &physics, const std::string &varname, const std::string &basefolder,
923 const RooLagrangianMorphFunc::ParamMap &inputParameters)
924{
925 TPair* pair = loadFromFileResidentFolder<TPair>(file,basefolder,varname,false);
926 if (!pair) return;
927 TParameter<double> *xsec_double = dynamic_cast<TParameter<double> *>(pair->Key());
928 if (xsec_double) {
929 collectCrosssections<double>(name, file, list_xs, physics, varname, inputParameters);
930 } else {
931 TParameter<float> *xsec_float = dynamic_cast<TParameter<float> *>(pair->Key());
932 if (xsec_float) {
933 collectCrosssections<float>(name, file, list_xs, physics, varname, inputParameters);
934 } else {
935 std::cerr << "cannot morph objects of class 'TPair' if parameter is not "
936 "double or float!"
937 << std::endl;
938 }
939 }
940}
941
942///////////////////////////////////////////////////////////////////////////////
943// FORMULA CALCULATION ////////////////////////////////////////////////////////
944///////////////////////////////////////////////////////////////////////////////
945///////////////////////////////////////////////////////////////////////////////
946
947////////////////////////////////////////////////////////////////////////////////
948/// recursive function to determine polynomials
949
950void collectPolynomialsHelper(const FeynmanDiagram &diagram, MorphFuncPattern &morphfunc, std::vector<int> &term,
951 int vertexid, bool first)
952{
953 if (vertexid > 0) {
954 for (size_t i = 0; i < diagram[vertexid - 1].size(); ++i) {
955 if (!diagram[vertexid - 1][i])
956 continue;
957 std::vector<int> newterm(term);
958 newterm[i]++;
959 if (first) {
960 ::collectPolynomialsHelper(diagram, morphfunc, newterm, vertexid, false);
961 } else {
962 ::collectPolynomialsHelper(diagram, morphfunc, newterm, vertexid - 1, true);
963 }
964 }
965 } else {
966 bool found = false;
967 for (size_t i = 0; i < morphfunc.size(); ++i) {
968 bool thisfound = true;
969 for (size_t j = 0; j < morphfunc[i].size(); ++j) {
970 if (morphfunc[i][j] != term[j]) {
971 thisfound = false;
972 break;
973 }
974 }
975 if (thisfound) {
976 found = true;
977 break;
978 }
979 }
980 if (!found) {
981 morphfunc.push_back(term);
982 }
983 }
984}
985
986////////////////////////////////////////////////////////////////////////////////
987/// calculate the morphing function pattern based on a vertex map
988
989void collectPolynomials(MorphFuncPattern &morphfunc, const FeynmanDiagram &diagram)
990{
991 int nvtx(diagram.size());
992 std::vector<int> term(diagram[0].size(), 0);
993
994 ::collectPolynomialsHelper(diagram, morphfunc, term, nvtx, true);
995}
996
997////////////////////////////////////////////////////////////////////////////////
998/// build a vertex map based on vertices and couplings appearing
999
1000template <class List>
1001inline void fillFeynmanDiagram(FeynmanDiagram &diagram, const std::vector<List *> &vertices, RooArgList &couplings)
1002{
1003 const int ncouplings = couplings.getSize();
1004 // std::cout << "Number of couplings " << ncouplings << std::endl;
1005 for (auto const &vertex : vertices) {
1006 std::vector<bool> vertexCouplings(ncouplings, false);
1007 int idx = -1;
1008 RooAbsReal *coupling;
1009 for (auto citr : couplings) {
1010 coupling = dynamic_cast<RooAbsReal *>(citr);
1011 idx++;
1012 if (!coupling) {
1013 std::cerr << "encountered invalid list of couplings in vertex!" << std::endl;
1014 return;
1015 }
1016 if (vertex->find(coupling->GetName())) {
1017 vertexCouplings[idx] = true;
1018 }
1019 }
1020 diagram.push_back(vertexCouplings);
1021 }
1022}
1023
1024////////////////////////////////////////////////////////////////////////////////
1025/// fill the matrix of coefficients
1026
1027template <class MatrixT, class T1, class T2>
1028inline MatrixT buildMatrixT(const RooLagrangianMorphFunc::ParamMap &inputParameters, const FormulaList &formulas,
1029 const T1 &args, const RooLagrangianMorphFunc::FlagMap &flagValues, const T2 &flags)
1030{
1031 const size_t dim = inputParameters.size();
1032 MatrixT matrix(dim, dim);
1033 int row = 0;
1034 for (auto sampleit : inputParameters) {
1035 const std::string sample(sampleit.first);
1036 // set all vars to value stored in input file
1037 if (!setParams<double>(sampleit.second, args, true, 0)) {
1038 std::cout << "unable to set parameters for sample " << sample << "!" << std::endl;
1039 }
1040 auto flagit = flagValues.find(sample);
1041 if (flagit != flagValues.end() && !setParams<int>(flagit->second, flags, true, 1)) {
1042 std::cout << "unable to set parameters for sample " << sample << "!" << std::endl;
1043 }
1044 // loop over all the formulas
1045 int col = 0;
1046 for (auto const &formula : formulas) {
1047 if (!formula.second) {
1048 std::cerr << "Error: invalid formula encountered!" << std::endl;
1049 }
1050 matrix(row, col) = formula.second->getVal();
1051 col++;
1052 }
1053 row++;
1054 }
1055 return matrix;
1056}
1057
1058////////////////////////////////////////////////////////////////////////////////
1059/// check if the matrix is square
1060
1061inline void checkMatrix(const RooLagrangianMorphFunc::ParamMap &inputParameters, const FormulaList &formulas)
1062{
1063 if (inputParameters.size() != formulas.size()) {
1064 std::stringstream ss;
1065 ss << "matrix is not square, consistency check failed: " << inputParameters.size() << " samples, "
1066 << formulas.size() << " expressions:" << std::endl;
1067 ss << "formulas: " << std::endl;
1068 for (auto const &formula : formulas) {
1069 ss << formula.second->GetTitle() << std::endl;
1070 }
1071 ss << "samples: " << std::endl;
1072 for (auto sample : inputParameters) {
1073 ss << sample.first << std::endl;
1074 }
1075 std::cerr << ss.str() << std::endl;
1076 }
1077}
1078
1079////////////////////////////////////////////////////////////////////////////////
1080/// check if the entries in the inverted matrix are sensible
1081
1082inline void inverseSanity(const Matrix &matrix, const Matrix &inverse, double &unityDeviation, double &largestWeight)
1083{
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}
1100
1101////////////////////////////////////////////////////////////////////////////////
1102/// check for name conflicts between the input samples and an argument set
1103template <class List>
1104inline void checkNameConflict(const RooLagrangianMorphFunc::ParamMap &inputParameters, List &args)
1105{
1106 for (auto sampleit : inputParameters) {
1107 const std::string sample(sampleit.first);
1108 RooAbsArg *arg = args.find(sample.c_str());
1109 if (arg) {
1110 std::cerr << "detected name conflict: cannot use sample '" << sample
1111 << "' - a parameter with the same name of type '" << arg->ClassName() << "' is present in set '"
1112 << args.GetName() << "'!" << std::endl;
1113 }
1114 }
1115}
1116
1117////////////////////////////////////////////////////////////////////////////////
1118/// build the formulas corresponding to the given set of input files and
1119/// the physics process
1120
1121template <class List>
1122inline FormulaList
1123buildFormulas(const char *mfname, const RooLagrangianMorphFunc::ParamMap &inputParameters,
1124 const RooLagrangianMorphFunc::FlagMap &inputFlags, const MorphFuncPattern &morphfunc,
1125 const RooArgList &couplings, const List &flags, const std::vector<List *> &nonInterfering)
1126{
1127 // example vbf hww:
1128 // Operators kSM, kHww, kAww, kHdwR,kHzz, kAzz
1129 // std::vector<bool> vertexProd = {true, true, true, true, true, true };
1130 // std::vector<bool> vertexDecay = {true, true, true, true, false,false};
1131 // diagram.push_back(vertexProd);
1132 // diagram.push_back(vertexDecay);
1133
1134 const int ncouplings = couplings.getSize();
1135 std::vector<bool> couplingsZero(ncouplings, true);
1136 std::map<TString, bool> flagsZero;
1137
1138 RooArgList operators;
1139 extractOperators(couplings, operators);
1140 size_t nOps = operators.getSize();
1141
1142 for (auto sampleit : inputParameters) {
1143 const std::string sample(sampleit.first);
1144 if (!setParams(sampleit.second, operators, true)) {
1145 std::cerr << "unable to set parameters for sample '" << sample << "'!" << std::endl;
1146 }
1147
1148 if ((int)nOps != (operators.getSize())) {
1149 std::cerr << "internal error, number of operators inconsistent!" << std::endl;
1150 }
1151
1152 RooAbsReal *obj0;
1153 int idx = 0;
1154
1155 for (auto itr1 : couplings) {
1156 obj0 = dynamic_cast<RooAbsReal *>(itr1);
1157 if (obj0->getVal() != 0) {
1158 couplingsZero[idx] = false;
1159 }
1160 idx++;
1161 }
1162 }
1163
1164 for (auto itr2 : flags) {
1165 auto obj1 = dynamic_cast<RooAbsReal *>(itr2);
1166 int nZero = 0;
1167 int nNonZero = 0;
1168 for (auto sampleit : inputFlags) {
1169 const auto &flag = sampleit.second.find(obj1->GetName());
1170 if (flag != sampleit.second.end()) {
1171 if (flag->second == 0.)
1172 nZero++;
1173 else
1174 nNonZero++;
1175 }
1176 }
1177 if (nZero > 0 && nNonZero == 0)
1178 flagsZero[obj1->GetName()] = true;
1179 else
1180 flagsZero[obj1->GetName()] = false;
1181 }
1182
1183 FormulaList formulas;
1184 for (size_t i = 0; i < morphfunc.size(); ++i) {
1185 RooArgList ss;
1186 bool isZero = false;
1187 std::string reason;
1188 // check if this is a blacklisted interference term
1189 for (const auto &group : nonInterfering) {
1190 int nInterferingOperators = 0;
1191 for (size_t j = 0; j < morphfunc[i].size(); ++j) {
1192 if (morphfunc[i][j] % 2 == 0)
1193 continue; // even exponents are not interference terms
1194 if (group->find(couplings.at(j)->GetName())) { // if the coupling is part of a
1195 // "pairwise non-interfering group"
1196 nInterferingOperators++;
1197 }
1198 }
1199 if (nInterferingOperators > 1) {
1200 isZero = true;
1201 reason = "blacklisted interference term!";
1202 }
1203 }
1204 int nNP = 0;
1205 if (!isZero) {
1206 // prepare the term
1207 for (size_t j = 0; j < morphfunc[i].size(); ++j) {
1208 const int exponent = morphfunc[i][j];
1209 if (exponent == 0)
1210 continue;
1211 RooAbsReal *coupling = dynamic_cast<RooAbsReal *>(couplings.at(j));
1212 for (int k = 0; k < exponent; ++k) {
1213 ss.add(*coupling);
1214 if (coupling->getAttribute("NewPhysics")) {
1215 nNP++;
1216 }
1217 }
1218 std::string cname(coupling->GetName());
1219 if (coupling->getAttribute("LO") && exponent > 1) {
1220 isZero = true;
1221 reason = "coupling " + cname + " was listed as leading-order-only";
1222 }
1223 // mark the term as zero if any of the couplings are zero
1224 if (!isZero && couplingsZero[j]) {
1225 isZero = true;
1226 reason = "coupling " + cname + " is zero!";
1227 }
1228 }
1229 }
1230 // check and apply flags
1231 bool removedByFlag = false;
1232
1233 for (auto itr : flags) {
1234 auto obj = dynamic_cast<RooAbsReal *>(itr);
1235 if (!obj)
1236 continue;
1237 TString sval(obj->getStringAttribute("NewPhysics"));
1238 int val = atoi(sval);
1239 if (val == nNP) {
1240 if (flagsZero.find(obj->GetName()) != flagsZero.end() && flagsZero.at(obj->GetName())) {
1241 removedByFlag = true;
1242 reason = Form("flag %s is zero", obj->GetName());
1243 }
1244 ss.add(*obj);
1245 }
1246 }
1247
1248 // create and add the formula
1249 if (!isZero && !removedByFlag) {
1250 // build the name
1251 const auto name = std::string(mfname) + "_pol" + std::to_string(i);
1252 formulas[i] = std::make_unique<RooProduct>(name.c_str(), ::concatNames(ss, " * ").c_str(), ss);
1253 }
1254 }
1255 return formulas;
1256}
1257
1258////////////////////////////////////////////////////////////////////////////////
1259/// create the weight formulas required for the morphing
1260
1261template <class T>
1262inline FormulaList
1263createFormulas(const char *name, const RooLagrangianMorphFunc::ParamMap &inputs,
1264 const RooLagrangianMorphFunc::FlagMap &inputFlags, const std::vector<std::vector<T *>> &diagrams,
1265 RooArgList &couplings, const T &flags, const std::vector<T *> &nonInterfering)
1266{
1267 MorphFuncPattern morphfuncpattern;
1268
1269 for (const auto &vertices : diagrams) {
1270 FeynmanDiagram d;
1271 ::fillFeynmanDiagram(d, vertices, couplings);
1272 ::collectPolynomials(morphfuncpattern, d);
1273 }
1274 FormulaList retval = buildFormulas(name, inputs, inputFlags, morphfuncpattern, couplings, flags, nonInterfering);
1275 if (retval.size() == 0) {
1276 std::stringstream errorMsgStream;
1277 errorMsgStream
1278 << "no formulas are non-zero, check if any if your couplings is floating and missing from your param_cards!"
1279 << std::endl;
1280 const auto errorMsg = errorMsgStream.str();
1281 throw std::runtime_error(errorMsg);
1282 }
1283 checkMatrix(inputs, retval);
1284 return retval;
1285}
1286
1287////////////////////////////////////////////////////////////////////////////////
1288/// build the sample weights required for the input templates
1289//
1290template <class T1>
1291inline void buildSampleWeights(T1 &weights, const char *fname, const RooLagrangianMorphFunc::ParamMap &inputParameters,
1292 FormulaList &formulas, const Matrix &inverse)
1293{
1294 int sampleidx = 0;
1295
1296 for (auto sampleit : inputParameters) {
1297 const std::string sample(sampleit.first);
1298 std::stringstream title;
1299 TString name_full(makeValidName(sample.c_str()));
1300 if (fname) {
1301 name_full.Append("_");
1302 name_full.Append(fname);
1303 name_full.Prepend("w_");
1304 }
1305
1306 int formulaidx = 0;
1307 // build the formula with the correct normalization
1308 RooLinearCombination *sampleformula = new RooLinearCombination(name_full.Data());
1309 for (auto const &formulait : formulas) {
1310 const RooFit::SuperFloat val(inverse(formulaidx, sampleidx));
1311 sampleformula->add(val, formulait.second.get());
1312 formulaidx++;
1313 }
1314 weights.add(*sampleformula);
1315 sampleidx++;
1316 }
1317}
1318
1319inline std::map<std::string, std::string>
1320buildSampleWeightStrings(const RooLagrangianMorphFunc::ParamMap &inputParameters, const FormulaList &formulas,
1321 const Matrix &inverse)
1322{
1323 int sampleidx = 0;
1324 std::map<std::string, std::string> weights;
1325 for (auto sampleit : inputParameters) {
1326 const std::string sample(sampleit.first);
1327 std::stringstream str;
1328 int formulaidx = 0;
1329 // build the formula with the correct normalization
1330 for (auto const &formulait : formulas) {
1331 double val(inverse(formulaidx, sampleidx));
1332 if (val != 0.) {
1333 if (formulaidx > 0 && val > 0)
1334 str << " + ";
1335 str << val << "*(" << formulait.second->GetTitle() << ")";
1336 }
1337 formulaidx++;
1338 }
1339 weights[sample] = str.str();
1340 sampleidx++;
1341 }
1342 return weights;
1343}
1344} // namespace
1345
1346///////////////////////////////////////////////////////////////////////////////
1347// CacheElem magic ////////////////////////////////////////////////////////////
1348///////////////////////////////////////////////////////////////////////////////
1349
1351public:
1352 std::unique_ptr<RooRealSumFunc> _sumFunc = nullptr;
1354
1355 FormulaList _formulas;
1357
1361
1363 virtual void operModeHook(RooAbsArg::OperMode) override{};
1364
1365 //////////////////////////////////////////////////////////////////////////////
1366 /// retrieve the list of contained args
1367
1369 {
1370 RooArgList args(*_sumFunc);
1371 args.add(_weights);
1372 args.add(_couplings);
1373 for (auto const &it : _formulas) {
1374 args.add(*(it.second));
1375 }
1376 return args;
1377 }
1378
1379 //////////////////////////////////////////////////////////////////////////////
1380 // default destructor
1381
1382 virtual ~CacheElem() {}
1383
1384 //////////////////////////////////////////////////////////////////////////////
1385 /// create the basic objects required for the morphing
1386
1387 inline void createComponents(const RooLagrangianMorphFunc::ParamMap &inputParameters,
1388 const RooLagrangianMorphFunc::FlagMap &inputFlags, const char *funcname,
1389 const std::vector<std::vector<RooListProxy *>> &diagrams,
1390 const std::vector<RooListProxy *> &nonInterfering, const RooListProxy &flags)
1391 {
1392 RooArgList operators;
1393 // coutD(ObjectHandling) << "collecting couplings" << std::endl;
1394 for (const auto &diagram : diagrams) {
1395 for (const auto &vertex : diagram) {
1396 extractCouplings(*vertex, this->_couplings);
1397 }
1398 }
1399 extractOperators(this->_couplings, operators);
1400 this->_formulas =
1401 ::createFormulas(funcname, inputParameters, inputFlags, diagrams, this->_couplings, flags, nonInterfering);
1402 }
1403
1404 //////////////////////////////////////////////////////////////////////////////
1405 /// build and invert the morphing matrix
1406 template <class List>
1407 inline void buildMatrix(const RooLagrangianMorphFunc::ParamMap &inputParameters,
1408 const RooLagrangianMorphFunc::FlagMap &inputFlags, const List &flags)
1409 {
1410 RooArgList operators;
1411 extractOperators(this->_couplings, operators);
1412 Matrix matrix(buildMatrixT<Matrix>(inputParameters, this->_formulas, operators, inputFlags, flags));
1413 if (size(matrix) < 1) {
1414 std::cerr << "input matrix is empty, please provide suitable input samples!" << std::endl;
1415 }
1416 Matrix inverse(diagMatrix(size(matrix)));
1417
1418 double condition = (double)(invertMatrix(matrix, inverse));
1419 double unityDeviation, largestWeight;
1420 inverseSanity(matrix, inverse, unityDeviation, largestWeight);
1421 bool weightwarning(largestWeight > morphLargestWeight ? true : false);
1422 bool unitywarning(unityDeviation > morphUnityDeviation ? true : false);
1423
1424 if (false) {
1425 if (unitywarning) {
1426 oocxcoutW((TObject *)0, Eval) << "Warning: The matrix inversion seems to be unstable. This can "
1427 "be a result to input samples that are not sufficiently "
1428 "different to provide any morphing power."
1429 << std::endl;
1430 } else if (weightwarning) {
1431 oocxcoutW((TObject *)0, Eval) << "Warning: Some weights are excessively large. This can be a "
1432 "result to input samples that are not sufficiently different to "
1433 "provide any morphing power."
1434 << std::endl;
1435 }
1436 oocxcoutW((TObject *)0, Eval) << " Please consider the couplings "
1437 "encoded in your samples to cross-check:"
1438 << std::endl;
1439 for (auto sampleit : inputParameters) {
1440 const std::string sample(sampleit.first);
1441 oocxcoutW((TObject *)0, Eval) << " " << sample << ": ";
1442 // set all vars to value stored in input file
1443 setParams(sampleit.second, operators, true);
1444 bool first = true;
1445 RooAbsReal *obj;
1446
1447 for (auto itr : this->_couplings) {
1448 obj = dynamic_cast<RooAbsReal *>(itr);
1449 if (!first)
1450 std::cerr << ", ";
1451 oocxcoutW((TObject *)0, Eval) << obj->GetName() << "=" << obj->getVal();
1452 first = false;
1453 }
1454 oocxcoutW((TObject *)0, Eval) << std::endl;
1455 }
1456 }
1457#ifndef USE_UBLAS
1458 this->_matrix.ResizeTo(matrix.GetNrows(), matrix.GetNrows());
1459 this->_inverse.ResizeTo(matrix.GetNrows(), matrix.GetNrows());
1460#endif
1461 this->_matrix = matrix;
1462 this->_inverse = inverse;
1463 this->_condition = condition;
1464 }
1465
1466 ////////////////////////////////////////////////////////////////////////////////
1467 /// build the final morphing function
1468
1469 inline void buildMorphingFunction(const char *name, const RooLagrangianMorphFunc::ParamMap &inputParameters,
1470 const std::map<std::string, int> &storage, const RooArgList &physics,
1471 bool allowNegativeYields, RooRealVar *observable, RooRealVar *binWidth)
1472 {
1473 if (!binWidth) {
1474 std::cerr << "invalid bin width given!" << std::endl;
1475 return;
1476 }
1477 if (!observable) {
1478 std::cerr << "invalid observable given!" << std::endl;
1479 return;
1480 }
1481
1482 RooArgList operators;
1483 extractOperators(this->_couplings, operators);
1484
1485 // retrieve the weights
1486 ::buildSampleWeights(this->_weights, name, inputParameters, this->_formulas, this->_inverse);
1487
1488 // build the products of element and weight for each sample
1489 size_t i = 0;
1490 RooArgList sumElements;
1491 RooArgList scaleElements;
1492 for (auto sampleit : inputParameters) {
1493 // for now, we assume all the lists are nicely ordered
1494 TString prodname(makeValidName(sampleit.first.c_str()));
1495
1496 RooAbsReal *obj = static_cast<RooAbsReal *>(physics.at(storage.at(prodname.Data())));
1497
1498 if (!obj) {
1499 std::cerr << "unable to access physics object for " << prodname << std::endl;
1500 return;
1501 }
1502
1503 RooAbsReal *weight = static_cast<RooAbsReal *>(this->_weights.at(i));
1504
1505 if (!weight) {
1506 std::cerr << "unable to access weight object for " << prodname << std::endl;
1507 return;
1508 }
1509 prodname.Append("_");
1510 prodname.Append(name);
1511 RooArgList prodElems(*weight, *obj);
1512
1513 allowNegativeYields = true;
1514 RooProduct *prod = new RooProduct(prodname, prodname, prodElems);
1515 if (!allowNegativeYields) {
1516 auto maxname = std::string(prodname) + "_max0";
1517 RooArgSet prodset(*prod);
1518
1519 RooFormulaVar *max = new RooFormulaVar(maxname.c_str(), "max(0," + prodname + ")", prodset);
1520 sumElements.add(*max);
1521 } else {
1522 sumElements.add(*prod);
1523 }
1524 scaleElements.add(*(binWidth));
1525 i++;
1526 }
1527
1528 // put everything together
1529 this->_sumFunc = make_unique<RooRealSumFunc>(Form("%s_morphfunc", name), name, sumElements, scaleElements);
1530
1531 if (!observable)
1532 std::cerr << "unable to access observable" << std::endl;
1533 this->_sumFunc.get()->addServer(*observable);
1534 if (!binWidth)
1535 std::cerr << "unable to access bin width" << std::endl;
1536 this->_sumFunc.get()->addServer(*binWidth);
1537 if (operators.getSize() < 1)
1538 std::cerr << "no operators listed" << std::endl;
1539 this->_sumFunc.get()->addServerList(operators);
1540 if (this->_weights.getSize() < 1)
1541 std::cerr << "unable to access weight objects" << std::endl;
1542 this->_sumFunc.get()->addOwnedComponents(this->_weights);
1543 this->_sumFunc.get()->addOwnedComponents(sumElements);
1544 this->_sumFunc.get()->addServerList(sumElements);
1545 this->_sumFunc.get()->addServerList(scaleElements);
1546
1547#ifdef USE_UBLAS
1548 std::cout.precision(std::numeric_limits<double>::digits);
1549#endif
1550 }
1551 //////////////////////////////////////////////////////////////////////////////
1552 /// create all the temporary objects required by the class
1553
1555 {
1556 std::string obsName = func->getObservable()->GetName();
1557 RooLagrangianMorphFunc::ParamSet values = getParams(func->_operators);
1558
1560 cache->createComponents(func->_config.paramCards, func->_config.flagValues, func->GetName(), func->_diagrams,
1561 func->_nonInterfering, func->_flags);
1562
1563 cache->buildMatrix(func->_config.paramCards, func->_config.flagValues, func->_flags);
1564 if (obsName.size() == 0) {
1565 std::cerr << "Matrix inversion succeeded, but no observable was "
1566 "supplied. quitting..."
1567 << std::endl;
1568 return cache;
1569 }
1570
1571 oocxcoutP((TObject *)0, ObjectHandling) << "observable: " << func->getObservable()->GetName() << std::endl;
1572 oocxcoutP((TObject *)0, ObjectHandling) << "binWidth: " << func->getBinWidth()->GetName() << std::endl;
1573
1574 setParams(func->_flags, 1);
1575 cache->buildMorphingFunction(func->GetName(), func->_config.paramCards, func->_sampleMap, func->_physics,
1576 func->_config.allowNegativeYields, func->getObservable(), func->getBinWidth());
1577 setParams(values, func->_operators, true);
1578 setParams(func->_flags, 1);
1579 return cache;
1580 }
1581
1582 //////////////////////////////////////////////////////////////////////////////
1583 /// create all the temporary objects required by the class
1584 /// function variant with precomputed inverse matrix
1585
1587 {
1588 RooLagrangianMorphFunc::ParamSet values = getParams(func->_operators);
1589
1591 cache->createComponents(func->_config.paramCards, func->_config.flagValues, func->GetName(), func->_diagrams,
1592 func->_nonInterfering, func->_flags);
1593
1594#ifndef USE_UBLAS
1595 cache->_inverse.ResizeTo(inverse.GetNrows(), inverse.GetNrows());
1596#endif
1597 cache->_inverse = inverse;
1598 cache->_condition = NaN;
1599
1600 setParams(func->_flags, 1);
1601 cache->buildMorphingFunction(func->GetName(), func->_config.paramCards, func->_sampleMap, func->_physics,
1602 func->_config.allowNegativeYields, func->getObservable(), func->getBinWidth());
1603 setParams(values, func->_operators, true);
1604 setParams(func->_flags, 1);
1605 return cache;
1606 }
1607};
1608
1609///////////////////////////////////////////////////////////////////////////////
1610// Class Implementation ///////////////////////////////////////////////////////
1611///////////////////////////////////////////////////////////////////////////////
1612
1613////////////////////////////////////////////////////////////////////////////////
1614/// write a matrix to a file
1615
1616void RooLagrangianMorphFunc::writeMatrixToFile(const TMatrixD &matrix, const char *fname)
1617{
1618 writeMatrixToFileT(matrix, fname);
1619}
1620
1621////////////////////////////////////////////////////////////////////////////////
1622/// write a matrix to a stream
1623
1624void RooLagrangianMorphFunc::writeMatrixToStream(const TMatrixD &matrix, std::ostream &stream)
1625{
1626 writeMatrixToStreamT(matrix, stream);
1627}
1628
1629////////////////////////////////////////////////////////////////////////////////
1630/// read a matrix from a text file
1631
1633{
1634 return readMatrixFromFileT<TMatrixD>(fname);
1635}
1636
1637////////////////////////////////////////////////////////////////////////////////
1638/// read a matrix from a stream
1639
1641{
1642 return readMatrixFromStreamT<TMatrixD>(stream);
1643}
1644
1645////////////////////////////////////////////////////////////////////////////////
1646/// setup observable, recycle existing observable if defined
1647
1648RooRealVar *RooLagrangianMorphFunc::setupObservable(const char *obsname, TClass *mode, TObject *inputExample)
1649{
1650 // cxcoutP(ObjectHandling) << "setting up observable" << std::endl;
1651 RooRealVar *obs = nullptr;
1652 Bool_t obsExists(false);
1653 if (this->_observables.at(0) != 0) {
1654 obs = (RooRealVar *)this->_observables.at(0);
1655 obsExists = true;
1656 }
1657
1658 if (mode && mode->InheritsFrom(RooHistFunc::Class())) {
1659 obs = (RooRealVar *)(dynamic_cast<RooHistFunc *>(inputExample)->getHistObsList().first());
1660 obsExists = true;
1661 this->_observables.add(*obs);
1662 } else if (mode && mode->InheritsFrom(RooParamHistFunc::Class())) {
1663 obs = (RooRealVar *)(dynamic_cast<RooParamHistFunc *>(inputExample)->paramList().first());
1664 obsExists = true;
1665 this->_observables.add(*obs);
1666 }
1667
1668 // Note: "found!" will be printed if s2 is a substring of s1, both s1 and s2
1669 // are of type std::string. s1.find(s2)
1670 // obtain the observable
1671 if (!obsExists) {
1672 if (mode && mode->InheritsFrom(TH1::Class())) {
1673 TH1 *hist = (TH1 *)(inputExample);
1674 obs = new RooRealVar(obsname, obsname, hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax());
1675 obs->setBins(hist->GetNbinsX());
1676 } else {
1677 obs = new RooRealVar(obsname, obsname, 0, 1);
1678 obs->setBins(1);
1679 }
1680 this->_observables.add(*obs);
1681 } else {
1682 if (strcmp(obsname, obs->GetName()) != 0) {
1683 coutW(ObjectHandling) << " name of existing observable " << this->_observables.at(0)->GetName()
1684 << " does not match expected name " << obsname << std::endl ;
1685 }
1686 }
1687
1688 TString sbw = Form("binWidth_%s", makeValidName(obs->GetName()).Data());
1689 RooRealVar *binWidth = new RooRealVar(sbw.Data(), sbw.Data(), 1.);
1690 double bw = obs->numBins() / (obs->getMax() - obs->getMin());
1691 binWidth->setVal(bw);
1692 binWidth->setConstant(true);
1693 this->_binWidths.add(*binWidth);
1694
1695 return obs;
1696}
1697
1698//#ifndef USE_MULTIPRECISION_LC
1699//#pragma GCC diagnostic push
1700//#pragma GCC diagnostic ignored "-Wunused-parameter"
1701//#endif
1702
1703////////////////////////////////////////////////////////////////////////////////
1704/// update sample weight (-?-)
1705
1707{
1708 //#ifdef USE_MULTIPRECISION_LC
1709 int sampleidx = 0;
1710 auto cache = this->getCache(_curNormSet);
1711 const size_t n(size(cache->_inverse));
1712 for (auto sampleit : this->_config.paramCards) {
1713 const std::string sample(sampleit.first);
1714 // build the formula with the correct normalization
1715 RooLinearCombination *sampleformula = dynamic_cast<RooLinearCombination *>(this->getSampleWeight(sample.c_str()));
1716 if (!sampleformula) {
1717 coutE(ObjectHandling) << Form("unable to access formula for sample '%s'!",sample.c_str()) << std::endl;
1718 return;
1719 }
1720 cxcoutP(ObjectHandling) << "updating formula for sample '" << sample << "'" << std::endl;
1721 for (size_t formulaidx = 0; formulaidx < n; ++formulaidx) {
1722 const RooFit::SuperFloat val(cache->_inverse(formulaidx, sampleidx));
1723#ifdef USE_UBLAS
1724 if (val != val) {
1725#else
1726 if (std::isnan(val)) {
1727#endif
1728 coutE(ObjectHandling) << "refusing to propagate NaN!" << std::endl;
1729 }
1730 cxcoutP(ObjectHandling) << " " << formulaidx << ":" << sampleformula->getCoefficient(formulaidx) << " -> "
1731 << val << std::endl;
1732 sampleformula->setCoefficient(formulaidx, val);
1733 assert(sampleformula->getCoefficient(formulaidx) == val);
1734 }
1735 sampleformula->setValueDirty();
1736 ++sampleidx;
1737 }
1738 //#else
1739 // ERROR("updating sample weights currently not possible without boost!");
1740 //#endif
1741}
1742//#ifndef USE_MULTIPRECISION_LC
1743//#pragma GCC diagnostic pop
1744//#endif
1745
1746////////////////////////////////////////////////////////////////////////////////
1747/// read the parameters from the input file
1748
1750{
1751 readValues<double>(this->_config.paramCards, f, this->_config.folderNames, "param_card", true);
1752 readValues<int>(this->_config.flagValues, f, this->_config.folderNames, "flags", false);
1753}
1754
1755////////////////////////////////////////////////////////////////////////////////
1756/// retrieve the physics inputs
1757
1759{
1760 std::string obsName = this->_config.observableName;
1761 cxcoutP(InputArguments) << "initializing physics inputs from file " << file->GetName() << " with object name(s) '"
1762 << obsName << "'" << std::endl;
1763 auto folderNames = this->_config.folderNames;
1764 TObject* obj = loadFromFileResidentFolder<TObject>(file, folderNames.front(), obsName, true);
1765 if (!obj) {
1766 std::cerr << "unable to locate object '" << obsName << "' in folder '" << folderNames.front() << "'!" << std::endl;
1767 return;
1768 }
1769 std::string classname = obj->ClassName();
1770 TClass *mode = TClass::GetClass(obj->ClassName());
1771
1772 RooRealVar *observable = this->setupObservable(obsName.c_str(), mode, obj);
1773 if (classname.find("TH1") != std::string::npos) {
1774 collectHistograms(this->GetName(), file, this->_sampleMap, this->_physics, *observable, obsName,
1775 this->_config.paramCards);
1776 } else if (classname.find("RooHistFunc") != std::string::npos ||
1777 classname.find("RooParamHistFunc") != std::string::npos ||
1778 classname.find("PiecewiseInterpolation") != std::string::npos) {
1779 collectRooAbsReal(this->GetName(), file, this->_sampleMap, this->_physics, obsName, this->_config.paramCards);
1780 } else if (classname.find("TParameter<double>") != std::string::npos) {
1781 collectCrosssections<double>(this->GetName(), file, this->_sampleMap, this->_physics, obsName,
1782 this->_config.paramCards);
1783 } else if (classname.find("TParameter<float>") != std::string::npos) {
1784 collectCrosssections<float>(this->GetName(), file, this->_sampleMap, this->_physics, obsName,
1785 this->_config.paramCards);
1786 } else if (classname.find("TPair") != std::string::npos) {
1787 collectCrosssectionsTPair(this->GetName(), file, this->_sampleMap, this->_physics, obsName, folderNames[0],
1788 this->_config.paramCards);
1789 } else {
1790 std::cerr << "cannot morph objects of class '" << mode->GetName() << "'!" << std::endl;
1791 }
1792}
1793
1794////////////////////////////////////////////////////////////////////////////////
1795/// convert the RooArgList folders into a simple vector of std::string
1796
1798{
1799 for (auto const &folder : folders) {
1800 RooStringVar *var = dynamic_cast<RooStringVar *>(folder);
1801 const std::string sample(var ? var->getVal() : folder->GetName());
1802 if (sample.empty())
1803 continue;
1804 this->_config.folderNames.push_back(sample);
1805 }
1806
1807 TDirectory *file = openFile(this->_config.fileName);
1808 TIter next(file->GetList());
1809 TObject *obj = nullptr;
1810 while ((obj = (TObject *)next())) {
1811 auto f = readOwningFolderFromFile(file,obj->GetName());
1812 if (!f)
1813 continue;
1814 std::string name(f->GetName());
1815 if (name.empty())
1816 continue;
1817 this->_config.folderNames.push_back(name);
1818 }
1819 closeFile(file);
1820}
1821
1822////////////////////////////////////////////////////////////////////////////////
1823/// print all the parameters and their values in the given sample to the console
1824
1825void RooLagrangianMorphFunc::printParameters(const char *samplename) const
1826{
1827 for (const auto &param : this->_config.paramCards.at(samplename)) {
1828 if (this->hasParameter(param.first.c_str())) {
1829 std::cout << param.first << " = " << param.second;
1830 if (this->isParameterConstant(param.first.c_str()))
1831 std::cout << " (const)";
1832 std::cout << std::endl;
1833 }
1834 }
1835}
1836
1837////////////////////////////////////////////////////////////////////////////////
1838/// print all the known samples to the console
1839
1841{
1842 // print all the known samples to the console
1843 for (auto folder : this->_config.folderNames) {
1844 std::cout << folder << std::endl;
1845 }
1846}
1847
1848////////////////////////////////////////////////////////////////////////////////
1849/// print the current physics values
1850
1852{
1853 for (const auto &sample : this->_sampleMap) {
1854 RooAbsArg *phys = this->_physics.at(sample.second);
1855 if (!phys)
1856 continue;
1857 phys->Print();
1858 }
1859}
1860
1861////////////////////////////////////////////////////////////////////////////////
1862/// constructor with proper arguments
1863
1864RooLagrangianMorphFunc::RooLagrangianMorphFunc(const char *name, const char *title, const Config &config)
1865 : RooAbsReal(name, title), _cacheMgr(this, 10, kTRUE, kTRUE),
1866 _operators("operators", "set of operators", this, kTRUE, kFALSE),
1867 _observables("observables", "morphing observables", this, kTRUE, kFALSE),
1868 _binWidths("binWidths", "set of binWidth objects", this, kTRUE, kFALSE), _config(config), _curNormSet(0)
1869{
1870 this->init();
1871 this->setup(false);
1872
1874}
1875
1876////////////////////////////////////////////////////////////////////////////////
1877/// constructor with proper arguments
1878RooLagrangianMorphFunc::RooLagrangianMorphFunc(const char *name, const char *title, const char *filename,
1879 const char *observableName, const RooArgSet &couplings,
1880 const RooArgSet &folders)
1881 : RooAbsReal(name, title), _cacheMgr(this, 10, kTRUE, kTRUE),
1882 _operators("operators", "set of operators", this, kTRUE, kFALSE),
1883 _observables("observables", "morphing observables", this, kTRUE, kFALSE),
1884 _binWidths("binWidths", "set of binWidth objects", this, kTRUE, kFALSE), _curNormSet(0)
1885{
1886 this->_config.fileName = filename;
1887 this->_config.observableName = observableName;
1888 this->_config.couplings.add(couplings);
1889 this->addFolders(folders);
1890 this->init();
1891 this->setup(false);
1892
1894}
1895
1896////////////////////////////////////////////////////////////////////////////////
1897/// setup this instance with the given set of operators and vertices
1898/// if own=true, the class will own the operatorsemplate <class Base>
1899
1901{
1902 this->_ownParameters = own;
1903 auto diagrams = this->_diagrams;
1904
1905 if (diagrams.size() > 0) {
1906 RooArgList operators;
1907 for (auto const &v : diagrams) {
1908 for (const RooArgList *t : v) {
1909 extractOperators(*t, operators);
1910 }
1911 }
1912
1913 if (own) {
1914 this->_operators.addOwned(operators);
1915 } else {
1916 this->_operators.add(operators);
1917 }
1918
1919 for (size_t j = 0; j < diagrams.size(); ++j) {
1920 std::vector<RooListProxy *> diagram;
1921 for (size_t i = 0; i < diagrams[j].size(); ++i) {
1922 std::stringstream name;
1923 name << "!vertex" << i;
1924 std::stringstream title;
1925 title << "set of couplings in the vertex " << i;
1926 diagram.push_back(new RooListProxy(name.str().c_str(), title.str().c_str(), this, kTRUE, kFALSE));
1927 if (own) {
1928 diagram[i]->addOwned(*diagrams[j][i]);
1929 } else {
1930 diagram[i]->add(*diagrams[j][i]);
1931 }
1932 }
1933 this->_diagrams.push_back(diagram);
1934 }
1935 }
1936
1937 else if (this->_config.couplings.size() > 0) {
1938 RooArgList operators;
1939 std::vector<RooListProxy *> vertices;
1940 extractOperators(this->_config.couplings, operators);
1941 vertices.push_back(new RooListProxy("!couplings", "set of couplings in the vertex", this, kTRUE, kFALSE));
1942 if (own) {
1943 this->_operators.addOwned(operators);
1944 vertices[0]->addOwned(this->_config.couplings);
1945 } else {
1946 this->_operators.add(operators);
1947 vertices[0]->add(this->_config.couplings);
1948 }
1949 this->_diagrams.push_back(vertices);
1950 }
1951
1952 else if (this->_config.prodCouplings.size() > 0 && this->_config.decCouplings.size() > 0) {
1953 std::vector<RooListProxy *> vertices;
1954 RooArgList operators;
1955 cxcoutP(InputArguments) << "prod/dec couplings provided" << std::endl;
1956 extractOperators(this->_config.prodCouplings, operators);
1957 extractOperators(this->_config.decCouplings, operators);
1958 vertices.push_back(
1959 new RooListProxy("!production", "set of couplings in the production vertex", this, kTRUE, kFALSE));
1960 vertices.push_back(new RooListProxy("!decay", "set of couplings in the decay vertex", this, kTRUE, kFALSE));
1961 if (own) {
1962 this->_operators.addOwned(operators);
1963 vertices[0]->addOwned(this->_config.prodCouplings);
1964 vertices[1]->addOwned(this->_config.decCouplings);
1965 } else {
1966 cxcoutP(InputArguments) << "adding non-own operators" << std::endl;
1967 this->_operators.add(operators);
1968 vertices[0]->add(this->_config.prodCouplings);
1969 vertices[1]->add(this->_config.decCouplings);
1970 }
1971 this->_diagrams.push_back(vertices);
1972 }
1973}
1974
1975////////////////////////////////////////////////////////////////////////////////
1976/// (-?-)
1977
1979{
1980 std::string filename = this->_config.fileName;
1981 TDirectory *file = openFile(filename.c_str());
1982 if (!file) {
1983 coutE(InputArguments) << "unable to open file '" << filename << "'!" << std::endl;
1984 return;
1985 }
1986 this->readParameters(file);
1987 checkNameConflict(this->_config.paramCards, this->_operators);
1988 this->collectInputs(file);
1989 closeFile(file);
1990 this->addServerList(this->_physics);
1991 RooRealVar *nNP0 = new RooRealVar("nNP0", "nNP0", 1., 0, 1.);
1992 nNP0->setStringAttribute("NewPhysics", "0");
1993 nNP0->setConstant(true);
1994 this->_flags.add(*nNP0);
1995 RooRealVar *nNP1 = new RooRealVar("nNP1", "nNP1", 1., 0, 1.);
1996 nNP1->setStringAttribute("NewPhysics", "1");
1997 nNP1->setConstant(true);
1998 this->_flags.add(*nNP1);
1999 RooRealVar *nNP2 = new RooRealVar("nNP2", "nNP2", 1., 0, 1.);
2000 nNP2->setStringAttribute("NewPhysics", "2");
2001 nNP2->setConstant(true);
2002 this->_flags.add(*nNP2);
2003 RooRealVar *nNP3 = new RooRealVar("nNP3", "nNP3", 1., 0, 1.);
2004 nNP3->setStringAttribute("NewPhysics", "3");
2005 nNP3->setConstant(true);
2006 this->_flags.add(*nNP3);
2007 RooRealVar *nNP4 = new RooRealVar("nNP4", "nNP4", 1., 0, 1.);
2008 nNP4->setStringAttribute("NewPhysics", "4");
2009 nNP4->setConstant(true);
2010 this->_flags.add(*nNP4);
2011}
2012
2013////////////////////////////////////////////////////////////////////////////////
2014/// copy constructor
2015
2017 : RooAbsReal(other, name), _cacheMgr(other._cacheMgr, this), _scale(other._scale), _sampleMap(other._sampleMap),
2018 _physics(other._physics.GetName(), this, other._physics),
2019 _operators(other._operators.GetName(), this, other._operators),
2020 _observables(other._observables.GetName(), this, other._observables),
2021 _binWidths(other._binWidths.GetName(), this, other._binWidths), _flags(other._flags.GetName(), this, other._flags),
2022 _config(other._config), _curNormSet(0)
2023{
2024 for (size_t j = 0; j < other._diagrams.size(); ++j) {
2025 std::vector<RooListProxy *> diagram;
2026 for (size_t i = 0; i < other._diagrams[j].size(); ++i) {
2027 RooListProxy *list = new RooListProxy(other._diagrams[j][i]->GetName(), this, *(other._diagrams[j][i]));
2028 diagram.push_back(list);
2029 }
2030 this->_diagrams.push_back(diagram);
2031 }
2033}
2034
2035////////////////////////////////////////////////////////////////////////////////
2036/// set energy scale of the EFT expansion
2037
2039{
2040 this->_scale = val;
2041}
2042
2043////////////////////////////////////////////////////////////////////////////////
2044/// get energy scale of the EFT expansion
2045
2047{
2048 return this->_scale;
2049}
2050
2051////////////////////////////////////////////////////////////////////////////////
2052// default constructor
2053
2055 : _cacheMgr(this, 10, true, true), _operators("operators", "set of operators", this, kTRUE, kFALSE),
2056 _observables("observable", "morphing observable", this, kTRUE, kFALSE),
2057 _binWidths("binWidths", "set of bin width objects", this, kTRUE, kFALSE)
2058{
2059 static int counter(0);
2060 counter++;
2062}
2063
2064////////////////////////////////////////////////////////////////////////////////
2065/// default destructor
2066
2068
2069////////////////////////////////////////////////////////////////////////////////
2070/// cloning method
2071
2072TObject *RooLagrangianMorphFunc::clone(const char *newname) const
2073{
2074 return new RooLagrangianMorphFunc(*this, newname);
2075}
2076
2077////////////////////////////////////////////////////////////////////////////////
2078/// calculate the number of samples needed to morph a bivertex, 2-2 physics
2079/// process
2080
2081int RooLagrangianMorphFunc::countSamples(int nprod, int ndec, int nboth)
2082{
2083 FeynmanDiagram diagram;
2084 std::vector<bool> prod;
2085 std::vector<bool> dec;
2086 for (int i = 0; i < nboth; ++i) {
2087 prod.push_back(true);
2088 dec.push_back(true);
2089 }
2090 for (int i = 0; i < nprod; ++i) {
2091 prod.push_back(true);
2092 dec.push_back(false);
2093 }
2094 for (int i = 0; i < ndec; ++i) {
2095 prod.push_back(false);
2096 dec.push_back(true);
2097 }
2098 diagram.push_back(prod);
2099 diagram.push_back(dec);
2100 MorphFuncPattern morphfuncpattern;
2101 ::collectPolynomials(morphfuncpattern, diagram);
2102 return morphfuncpattern.size();
2103}
2104
2105////////////////////////////////////////////////////////////////////////////////
2106/// calculate the number of samples needed to morph a certain physics process
2107/// usage:
2108/// countSamples ( { RooLagrangianMorphFunc::makeHCggfCouplings(),
2109/// RooLagrangianMorphFunc::makeHCHZZCouplings() } )
2110
2111int RooLagrangianMorphFunc::countSamples(std::vector<RooArgList *> &vertices)
2112{
2113 RooArgList operators, couplings;
2114 for (auto vertex : vertices) {
2115 extractOperators(*vertex, operators);
2116 extractCouplings(*vertex, couplings);
2117 }
2118 FeynmanDiagram diagram;
2119 ::fillFeynmanDiagram(diagram, vertices, couplings);
2120 MorphFuncPattern morphfuncpattern;
2121 ::collectPolynomials(morphfuncpattern, diagram);
2122 return morphfuncpattern.size();
2123}
2124
2125////////////////////////////////////////////////////////////////////////////////
2126/// create TPair containers of the type expected by the RooLagrangianMorph
2127
2129{
2130 TPair *v = new TPair(new TParameter<double>("xsection", xs), new TParameter<double>("uncertainty", unc));
2131 return v;
2132}
2133
2134////////////////////////////////////////////////////////////////////////////////
2135/// create only the weight formulas. static function for external usage.
2136
2137std::map<std::string, std::string>
2139 const std::vector<std::vector<std::string>> &vertices_str)
2140{
2141 std::stack<RooArgList> ownedVertices;
2142 std::vector<RooArgList *> vertices;
2143 RooArgList couplings;
2144 for (const auto &vtx : vertices_str) {
2145 ownedVertices.emplace();
2146 auto &vertex = ownedVertices.top();
2147 for (const auto &c : vtx) {
2148 RooRealVar *coupling = (RooRealVar *)(couplings.find(c.c_str()));
2149 if (!coupling) {
2150 coupling = new RooRealVar(c.c_str(), c.c_str(), 1., 0., 10.);
2151 couplings.add(*coupling);
2152 }
2153 vertex.add(*coupling);
2154 }
2155 vertices.push_back(&vertex);
2156 }
2157 auto retval = RooLagrangianMorphFunc::createWeightStrings(inputs, vertices, couplings);
2158 return retval;
2159}
2160
2161////////////////////////////////////////////////////////////////////////////////
2162/// create only the weight formulas. static function for external usage.
2163
2164std::map<std::string, std::string>
2166 const std::vector<RooArgList *> &vertices, RooArgList &couplings)
2167{
2168 return createWeightStrings(inputs, vertices, couplings, {}, {}, {});
2169}
2170
2171////////////////////////////////////////////////////////////////////////////////
2172/// create only the weight formulas. static function for external usage.
2173
2174std::map<std::string, std::string>
2176 const std::vector<RooArgList *> &vertices, RooArgList &couplings,
2177 const RooLagrangianMorphFunc::FlagMap &flagValues, const RooArgList &flags,
2178 const std::vector<RooArgList *> &nonInterfering)
2179{
2180 FormulaList formulas = ::createFormulas("", inputs, flagValues, {vertices}, couplings, flags, nonInterfering);
2181 RooArgSet operators;
2182 extractOperators(couplings, operators);
2183 Matrix matrix(::buildMatrixT<Matrix>(inputs, formulas, operators, flagValues, flags));
2184 if (size(matrix) < 1) {
2185 std::cerr << "input matrix is empty, please provide suitable input samples!" << std::endl;
2186 }
2187 Matrix inverse(::diagMatrix(size(matrix)));
2188 double condition __attribute__((unused)) = (double)(invertMatrix(matrix, inverse));
2189 auto retval = buildSampleWeightStrings(inputs, formulas, inverse);
2190 return retval;
2191}
2192
2193////////////////////////////////////////////////////////////////////////////////
2194/// create only the weight formulas. static function for external usage.
2195
2197 const std::vector<RooArgList *> &vertices, RooArgList &couplings,
2198 const RooLagrangianMorphFunc::FlagMap &flagValues,
2199 const RooArgList &flags,
2200 const std::vector<RooArgList *> &nonInterfering)
2201{
2202 FormulaList formulas = ::createFormulas("", inputs, flagValues, {vertices}, couplings, flags, nonInterfering);
2203 RooArgSet operators;
2204 extractOperators(couplings, operators);
2205 Matrix matrix(::buildMatrixT<Matrix>(inputs, formulas, operators, flagValues, flags));
2206 if (size(matrix) < 1) {
2207 std::cerr << "input matrix is empty, please provide suitable input samples!" << std::endl;
2208 }
2209 Matrix inverse(::diagMatrix(size(matrix)));
2210 double condition __attribute__((unused)) = (double)(invertMatrix(matrix, inverse));
2211 RooArgSet retval;
2212 ::buildSampleWeights(retval, (const char *)nullptr /* name */, inputs, formulas, inverse);
2213 return retval;
2214}
2215
2216////////////////////////////////////////////////////////////////////////////////
2217/// create only the weight formulas. static function for external usage.
2218
2220 const std::vector<RooArgList *> &vertices, RooArgList &couplings)
2221{
2222 std::vector<RooArgList *> nonInterfering;
2223 RooArgList flags;
2224 FlagMap flagValues;
2225 return RooLagrangianMorphFunc::createWeights(inputs, vertices, couplings, flagValues, flags, nonInterfering);
2226}
2227
2228////////////////////////////////////////////////////////////////////////////////
2229/// return the RooProduct that is the element of the RooRealSumPdfi
2230/// corresponding to the given sample name
2231
2233{
2234 auto mf = this->getFunc();
2235 if (!mf) {
2236 coutE(Eval) << "unable to retrieve morphing function" << std::endl;
2237 return nullptr;
2238 }
2239 RooArgSet *args = mf->getComponents();
2240 TString prodname(name);
2241 prodname.Append("_");
2242 prodname.Append(this->GetName());
2243
2244 for (auto itr : *args) {
2245 RooProduct *prod = dynamic_cast<RooProduct *>(itr);
2246 if (!prod)
2247 continue;
2248 TString sname(prod->GetName());
2249 if (sname.CompareTo(prodname) == 0) {
2250 return prod;
2251 }
2252 }
2253 return nullptr;
2254}
2255////////////////////////////////////////////////////////////////////////////////
2256/// return the vector of sample names, used to build the morph func
2257
2258std::vector<std::string> RooLagrangianMorphFunc::getSamples() const
2259{
2260 return this->_config.folderNames;
2261}
2262
2263////////////////////////////////////////////////////////////////////////////////
2264/// retrieve the weight (prefactor) of a sample with the given name
2265
2267{
2268 auto cache = this->getCache(_curNormSet);
2269 auto wname = std::string("w_") + name + "_" + this->GetName();
2270 return dynamic_cast<RooAbsReal *>(cache->_weights.find(wname.c_str()));
2271}
2272
2273////////////////////////////////////////////////////////////////////////////////
2274/// print the current sample weights
2275
2277{
2278 this->printSampleWeights();
2279}
2280
2281////////////////////////////////////////////////////////////////////////////////
2282/// print the current sample weights
2283
2285{
2286 auto *cache = this->getCache(this->_curNormSet);
2287 for (const auto &sample : this->_sampleMap) {
2288 auto weightName = std::string("w_") + sample.first + "_" + this->GetName();
2289 auto weight = static_cast<RooAbsReal *>(cache->_weights.find(weightName.c_str()));
2290 if (!weight)
2291 continue;
2292 }
2293}
2294
2295////////////////////////////////////////////////////////////////////////////////
2296/// randomize the parameters a bit
2297/// useful to test and debug fitting
2298
2300{
2301 RooRealVar *obj;
2302 TRandom3 r;
2303
2304 for (auto itr : _operators) {
2305 obj = dynamic_cast<RooRealVar *>(itr);
2306 double val = obj->getVal();
2307 if (obj->isConstant())
2308 continue;
2309 double variation = r.Gaus(1, z);
2310 obj->setVal(val * variation);
2311 }
2312}
2313
2314////////////////////////////////////////////////////////////////////////////////
2315/// retrive the new physics objects and update the weights in the morphing
2316/// function
2317
2319{
2320 auto cache = this->getCache(_curNormSet);
2321
2322 std::string filename = this->_config.fileName;
2323 TDirectory *file = openFile(filename.c_str());
2324 if (!file) {
2325 coutE(InputArguments) << "unable to open file '" << filename << "'!" << std::endl;
2326 return false;
2327 }
2328
2329 this->readParameters(file);
2330
2331 checkNameConflict(this->_config.paramCards, this->_operators);
2332 this->collectInputs(file);
2333
2334 cache->buildMatrix(this->_config.paramCards, this->_config.flagValues, this->_flags);
2335 this->updateSampleWeights();
2336
2337 closeFile(file);
2338 return true;
2339}
2340
2341////////////////////////////////////////////////////////////////////////////////
2342/// setup the morphing function with a predefined inverse matrix
2343/// call this function *before* any other after creating the object
2344
2346{
2347 auto cache = static_cast<RooLagrangianMorphFunc::CacheElem *>(_cacheMgr.getObj(0, (RooArgSet *)0));
2348 Matrix m = makeSuperMatrix(inverse);
2349 if (cache) {
2350 std::string filename = this->_config.fileName;
2351 cache->_inverse = m;
2352 TDirectory *file = openFile(filename.c_str());
2353 if (!file) {
2354 coutE(InputArguments) << "unable to open file '" << filename << "'!" << std::endl;
2355 return false;
2356 }
2357
2358 this->readParameters(file);
2359 checkNameConflict(this->_config.paramCards, this->_operators);
2360 this->collectInputs(file);
2361
2362 // then, update the weights in the morphing function
2363 this->updateSampleWeights();
2364
2365 closeFile(file);
2366 } else {
2368 if (!cache)
2369 coutE(Caching) << "unable to create cache!" << std::endl;
2370 this->_cacheMgr.setObj(0, 0, cache, 0);
2371 }
2372 return true;
2373}
2374
2375////////////////////////////////////////////////////////////////////////////////
2376// setup the morphing function with a predefined inverse matrix
2377// call this function *before* any other after creating the object
2378
2380{
2381 auto cache = static_cast<RooLagrangianMorphFunc::CacheElem *>(_cacheMgr.getObj(0, (RooArgSet *)0));
2382 if (cache) {
2383 return false;
2384 }
2385 cache = RooLagrangianMorphFunc::CacheElem::createCache(this, readMatrixFromFileT<Matrix>(filename));
2386 if (!cache)
2387 coutE(Caching) << "unable to create cache!" << std::endl;
2388 this->_cacheMgr.setObj(0, 0, cache, 0);
2389 return true;
2390}
2391
2392////////////////////////////////////////////////////////////////////////////////
2393/// write the inverse matrix to a file
2394
2396{
2397 auto cache = this->getCache(_curNormSet);
2398 if (!cache)
2399 return false;
2400 writeMatrixToFileT(cache->_inverse, filename);
2401 return true;
2402}
2403
2404////////////////////////////////////////////////////////////////////////////////
2405/// retrieve the cache object
2406
2408{
2409 auto cache = static_cast<RooLagrangianMorphFunc::CacheElem *>(_cacheMgr.getObj(0, (RooArgSet *)0));
2410 if (!cache) {
2411 cxcoutP(Caching) << "creating cache from getCache function for " << this << std::endl;
2412 cxcoutP(Caching) << "current storage has size " << this->_sampleMap.size() << std::endl;
2414 if (cache)
2415 this->_cacheMgr.setObj(0, 0, cache, 0);
2416 else
2417 coutE(Caching) << "unable to create cache!" << std::endl;
2418 }
2419 return cache;
2420}
2421
2422////////////////////////////////////////////////////////////////////////////////
2423/// return true if a cache object is present, false otherwise
2424
2426{
2427 return (bool)(_cacheMgr.getObj(0, (RooArgSet *)0));
2428}
2429
2430////////////////////////////////////////////////////////////////////////////////
2431/// set one parameter to a specific value
2432
2433void RooLagrangianMorphFunc::setParameter(const char *name, double value)
2434{
2435 RooRealVar *param = this->getParameter(name);
2436 if (!param) {
2437 return;
2438 }
2439 if (value > param->getMax())
2440 param->setMax(value);
2441 if (value < param->getMin())
2442 param->setMin(value);
2443 param->setVal(value);
2444}
2445
2446////////////////////////////////////////////////////////////////////////////////
2447/// set one flag to a specific value
2448
2449void RooLagrangianMorphFunc::setFlag(const char *name, double value)
2450{
2451 RooRealVar *param = this->getFlag(name);
2452 if (!param) {
2453 return;
2454 }
2455 param->setVal(value);
2456}
2457
2458////////////////////////////////////////////////////////////////////////////////
2459/// set one parameter to a specific value and range
2460
2461void RooLagrangianMorphFunc::setParameter(const char *name, double value, double min, double max)
2462{
2463 RooRealVar *param = this->getParameter(name);
2464 if (!param) {
2465 return;
2466 }
2467 param->setMin(min);
2468 param->setMax(max);
2469 param->setVal(value);
2470}
2471
2472////////////////////////////////////////////////////////////////////////////////
2473/// set one parameter to a specific value and range
2474void RooLagrangianMorphFunc::setParameter(const char *name, double value, double min, double max, double error)
2475{
2476 RooRealVar *param = this->getParameter(name);
2477 if (!param) {
2478 return;
2479 }
2480 param->setMin(min);
2481 param->setMax(max);
2482 param->setVal(value);
2483 param->setError(error);
2484}
2485
2486////////////////////////////////////////////////////////////////////////////////
2487/// return true if the parameter with the given name is set constant, false
2488/// otherwise
2489
2491{
2492 RooRealVar *param = this->getParameter(name);
2493 if (param) {
2494 return param->isConstant();
2495 }
2496 return true;
2497}
2498
2499////////////////////////////////////////////////////////////////////////////////
2500/// retrieve the RooRealVar object incorporating the parameter with the given
2501/// name
2503{
2504
2505 return dynamic_cast<RooRealVar *>(this->_operators.find(name));
2506}
2507
2508////////////////////////////////////////////////////////////////////////////////
2509/// retrieve the RooRealVar object incorporating the flag with the given name
2510
2512{
2513 return dynamic_cast<RooRealVar *>(this->_flags.find(name));
2514}
2515
2516////////////////////////////////////////////////////////////////////////////////
2517/// check if a parameter of the given name is contained in the list of known
2518/// parameters
2519
2521{
2522 return this->getParameter(name);
2523}
2524
2525////////////////////////////////////////////////////////////////////////////////
2526/// call setConstant with the boolean argument provided on the parameter with
2527/// the given name
2528
2529void RooLagrangianMorphFunc::setParameterConstant(const char *name, bool constant) const
2530{
2531 RooRealVar *param = this->getParameter(name);
2532 if (param) {
2533 return param->setConstant(constant);
2534 }
2535}
2536
2537////////////////////////////////////////////////////////////////////////////////
2538/// set one parameter to a specific value
2539
2541{
2542 RooRealVar *param = this->getParameter(name);
2543 if (param) {
2544 return param->getVal();
2545 }
2546 return 0;
2547}
2548
2549////////////////////////////////////////////////////////////////////////////////
2550/// set the morphing parameters to those supplied in the given param hist
2551
2553{
2554 setParams(paramhist, this->_operators, false);
2555}
2556
2557////////////////////////////////////////////////////////////////////////////////
2558/// set the morphing parameters to those supplied in the sample with the given
2559/// name
2560
2561void RooLagrangianMorphFunc::setParameters(const char *foldername)
2562{
2563 std::string filename = this->_config.fileName;
2564 TDirectory *file = openFile(filename.c_str());
2565 TH1 *paramhist = loadFromFileResidentFolder<TH1>(file, foldername,"param_card");
2566 setParams(paramhist, this->_operators, false);
2567 closeFile(file);
2568}
2569
2570/////////////////////////////////////////////////////////////////////////////////
2571/// retrieve the morphing parameters associated to the sample with the given
2572/// name
2573
2575{
2576 const std::string name(foldername);
2577 return this->_config.paramCards.at(name);
2578}
2579
2580////////////////////////////////////////////////////////////////////////////////
2581/// set the morphing parameters to those supplied in the list with the given
2582/// name
2583
2585{
2586 for (auto itr : *list) {
2587 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
2588 if (!param)
2589 continue;
2590 this->setParameter(param->GetName(), param->getVal());
2591 }
2592}
2593
2594////////////////////////////////////////////////////////////////////////////////
2595/// retrieve the histogram observable
2596
2598{
2599 if (this->_observables.getSize() < 1) {
2600 coutE(InputArguments) << "observable not available!" << std::endl;
2601 return nullptr;
2602 }
2603 return static_cast<RooRealVar *>(this->_observables.at(0));
2604}
2605
2606////////////////////////////////////////////////////////////////////////////////
2607/// retrieve the histogram observable
2608
2610{
2611 if (this->_binWidths.getSize() < 1) {
2612 coutE(InputArguments) << "bin width not available!" << std::endl;
2613 return nullptr;
2614 }
2615 return static_cast<RooRealVar *>(this->_binWidths.at(0));
2616}
2617
2618////////////////////////////////////////////////////////////////////////////////
2619/// retrieve a histogram output of the current morphing settings
2620
2622 /*RooFitResult *r*/)
2623{
2624 return this->createTH1(name, false /*r*/);
2625}
2626
2627////////////////////////////////////////////////////////////////////////////////
2628/// retrieve a histogram output of the current morphing settings
2629
2630TH1 *RooLagrangianMorphFunc::createTH1(const std::string &name, bool correlateErrors
2631 /* RooFitResult *r*/)
2632{
2633 auto mf = std::make_unique<RooRealSumFunc>(*(this->getFunc()));
2634 RooRealVar *observable = this->getObservable();
2635
2636 const int nbins = observable->getBins();
2637
2638 TH1 *hist = new TH1F(name.c_str(), name.c_str(), nbins, observable->getBinning().array());
2639
2640 RooArgSet *args = mf->getComponents();
2641 for (int i = 0; i < nbins; ++i) {
2642 observable->setBin(i);
2643 double val = 0;
2644 double unc2 = 0;
2645 double unc = 0;
2646 for (auto itr : *args) {
2647 RooProduct *prod = dynamic_cast<RooProduct *>(itr);
2648 if (!prod)
2649 continue;
2650 RooAbsArg *phys = prod->components().find(Form("phys_%s", prod->GetName()));
2651 RooHistFunc *hf = dynamic_cast<RooHistFunc *>(phys);
2652 if (!hf) {
2653 continue;
2654 }
2655 const RooDataHist &dhist = hf->dataHist();
2656 dhist.get(i);
2657 RooAbsReal *formula = dynamic_cast<RooAbsReal *>(prod->components().find(Form("w_%s", prod->GetName())));
2658 double weight = formula->getVal();
2659 unc2 += dhist.weightSquared() * weight * weight;
2660 unc += sqrt(dhist.weightSquared()) * weight;
2661 val += dhist.weight() * weight;
2662 }
2663 hist->SetBinContent(i + 1, val);
2664 hist->SetBinError(i + 1, correlateErrors ? unc : sqrt(unc2));
2665 }
2666 return hist;
2667}
2668
2669////////////////////////////////////////////////////////////////////////////////
2670/// count the number of formulas that correspond to the current parameter set
2671
2673{
2674 int nFormulas = 0;
2675 auto mf = std::make_unique<RooRealSumFunc>(*(this->getFunc()));
2676 if (!mf)
2677 coutE(InputArguments) << "unable to retrieve morphing function" << std::endl;
2678 RooArgSet *args = mf->getComponents();
2679 for (auto itr : *args) {
2680 RooProduct *prod = dynamic_cast<RooProduct *>(itr);
2681 if (prod->getVal() != 0) {
2682 nFormulas++;
2683 }
2684 }
2685 return nFormulas;
2686}
2687
2688////////////////////////////////////////////////////////////////////////////////
2689/// check if there is any morphing power provided for the given parameter
2690/// morphing power is provided as soon as any two samples provide different,
2691/// non-zero values for this parameter
2692
2693bool RooLagrangianMorphFunc::isParameterUsed(const char *paramname) const
2694{
2695 std::string pname(paramname);
2696 double val = 0;
2697 bool isUsed = false;
2698 for (const auto &sample : this->_config.paramCards) {
2699 double thisval = sample.second.at(pname);
2700 if (thisval != val) {
2701 if (val != 0)
2702 isUsed = true;
2703 val = thisval;
2704 }
2705 }
2706 return isUsed;
2707}
2708
2709////////////////////////////////////////////////////////////////////////////////
2710/// check if there is any morphing power provided for the given coupling
2711/// morphing power is provided as soon as any two samples provide
2712/// different, non-zero values for this coupling
2713
2715{
2716 std::string cname(couplname);
2717 const RooArgList *args = this->getCouplingSet();
2718 RooAbsReal *coupling = dynamic_cast<RooAbsReal *>(args->find(couplname));
2719 if (!coupling)
2720 return false;
2722 double val = 0;
2723 bool isUsed = false;
2724 for (const auto &sample : this->_config.paramCards) {
2725 this->setParameters(sample.second);
2726 double thisval = coupling->getVal();
2727 if (thisval != val) {
2728 if (val != 0)
2729 isUsed = true;
2730 val = thisval;
2731 }
2732 }
2733 this->setParameters(params);
2734 return isUsed;
2735}
2736
2737////////////////////////////////////////////////////////////////////////////////
2738/// return the number of parameters in this morphing function
2739
2741{
2742 return this->getParameterSet()->getSize();
2743}
2744
2745////////////////////////////////////////////////////////////////////////////////
2746/// return the number of samples in this morphing function
2747
2749{
2750 // return the number of samples in this morphing function
2751 auto cache = getCache(_curNormSet);
2752 return cache->_formulas.size();
2753}
2754
2755////////////////////////////////////////////////////////////////////////////////
2756/// print the contributing samples and their respective weights
2757
2759{
2760 auto mf = std::make_unique<RooRealSumFunc>(*(this->getFunc()));
2761 if (!mf) {
2762 std::cerr << "Error: unable to retrieve morphing function" << std::endl;
2763 return;
2764 }
2765 RooArgSet *args = mf->getComponents();
2766 for (auto itr : *args) {
2767 RooAbsReal *formula = dynamic_cast<RooAbsReal *>(itr);
2768 if (formula) {
2769 TString name(formula->GetName());
2770 name.Remove(0, 2);
2771 name.Prepend("phys_");
2772 if (!args->find(name.Data())) {
2773 continue;
2774 }
2775 double val = formula->getVal();
2776 if (val != 0) {
2777 std::cout << formula->GetName() << ": " << val << " = " << formula->GetTitle() << std::endl;
2778 }
2779 }
2780 }
2781}
2782
2783////////////////////////////////////////////////////////////////////////////////
2784/// get the set of parameters
2785
2787{
2788 return &(this->_operators);
2789}
2790
2791////////////////////////////////////////////////////////////////////////////////
2792/// get the set of couplings
2793
2795{
2796 auto cache = getCache(_curNormSet);
2797 return &(cache->_couplings);
2798}
2799
2800////////////////////////////////////////////////////////////////////////////////
2801/// retrieve a set of couplings (-?-)
2802
2804{
2806 for (auto obj : *(this->getCouplingSet())) {
2807 RooAbsReal *var = dynamic_cast<RooAbsReal *>(obj);
2808 if (!var)
2809 continue;
2810 const std::string name(var->GetName());
2811 double val = var->getVal();
2812 couplings[name] = val;
2813 }
2814 return couplings;
2815}
2816
2817////////////////////////////////////////////////////////////////////////////////
2818/// retrieve the parameter set
2819
2821{
2822 return getParams(this->_operators);
2823}
2824
2825////////////////////////////////////////////////////////////////////////////////
2826/// retrieve a set of couplings (-?-)
2827
2829{
2830 setParams(params, this->_operators, false);
2831}
2832
2833////////////////////////////////////////////////////////////////////////////////
2834/// (currently similar to cloning the Pdf
2835
2836std::unique_ptr<RooWrapperPdf> RooLagrangianMorphFunc::createPdf() const
2837{
2838 auto cache = getCache(_curNormSet);
2839 auto func = std::make_unique<RooRealSumFunc>(*(cache->_sumFunc));
2840
2841 // create a wrapper on the roorealsumfunc
2842 return std::make_unique<RooWrapperPdf>(Form("pdf_%s", func->GetName()), Form("pdf of %s", func->GetTitle()), *func);
2843}
2844
2845////////////////////////////////////////////////////////////////////////////////
2846/// get the func
2847
2849{
2850 auto cache = getCache(_curNormSet);
2851 return cache->_sumFunc.get();
2852}
2853
2854////////////////////////////////////////////////////////////////////////////////
2855/// return extended mored capabilities
2856
2858{
2859 return this->createPdf()->extendMode();
2860}
2861
2862////////////////////////////////////////////////////////////////////////////////
2863/// return expected number of events for extended likelihood calculation,
2864/// this is the sum of all coefficients
2865
2867{
2868 return this->createPdf()->expectedEvents(nset);
2869}
2870
2871////////////////////////////////////////////////////////////////////////////////
2872/// return the number of expected events for the current parameter set
2873
2875{
2876 RooArgSet set;
2877 set.add(*this->getObservable());
2878 return this->createPdf()->expectedEvents(set);
2879}
2880
2881////////////////////////////////////////////////////////////////////////////////
2882/// return expected number of events for extended likelihood calculation,
2883/// this is the sum of all coefficients
2884
2886{
2887 return createPdf()->expectedEvents(&nset);
2888}
2889
2890////////////////////////////////////////////////////////////////////////////////
2891/// return the expected uncertainty for the current parameter set
2892
2894{
2895 RooRealVar *observable = this->getObservable();
2896 auto cache = this->getCache(_curNormSet);
2897 double unc2 = 0;
2898 for (const auto &sample : this->_sampleMap) {
2899 RooAbsArg *phys = this->_physics.at(sample.second);
2900 auto weightName = std::string("w_") + sample.first + "_" + this->GetName();
2901 auto weight = static_cast<RooAbsReal *>(cache->_weights.find(weightName.c_str()));
2902 if (!weight) {
2903 coutE(InputArguments) << "unable to find object " + weightName << std::endl;
2904 return 0.0;
2905 }
2906 double newunc2 = 0;
2907 RooHistFunc *hf = dynamic_cast<RooHistFunc *>(phys);
2908 RooRealVar *rv = dynamic_cast<RooRealVar *>(phys);
2909 if (hf) {
2910 const RooDataHist &hist = hf->dataHist();
2911 for (Int_t j = 0; j < observable->getBins(); ++j) {
2912 hist.get(j);
2913 newunc2 += hist.weightSquared();
2914 }
2915 } else if (rv) {
2916 newunc2 = pow(rv->getError(), 2);
2917 }
2918 double w = weight->getVal();
2919 unc2 += newunc2 * w * w;
2920 // std::cout << phys->GetName() << " : " << weight->GetName() << "
2921 // thisweight: " << w << " thisxsec2: " << newunc2 << " weight " << weight
2922 // << std::endl;
2923 }
2924 return sqrt(unc2);
2925}
2926
2927////////////////////////////////////////////////////////////////////////////////
2928/// print the parameters and their current values
2929
2931{
2932 // print the parameters and their current values
2933 for (auto obj : this->_operators) {
2934 RooRealVar *param = static_cast<RooRealVar *>(obj);
2935 if (!param)
2936 continue;
2937 param->Print();
2938 }
2939}
2940
2941////////////////////////////////////////////////////////////////////////////////
2942/// print the flags and their current values
2943
2945{
2946 for (auto flag : this->_flags) {
2947 RooRealVar *param = static_cast<RooRealVar *>(flag);
2948 if (!param)
2949 continue;
2950 param->Print();
2951 }
2952}
2953
2954////////////////////////////////////////////////////////////////////////////////
2955/// print a set of couplings
2956
2958{
2960 for (auto c : couplings) {
2961 std::cout << c.first << ": " << c.second << std::endl;
2962 }
2963}
2964
2965////////////////////////////////////////////////////////////////////////////////
2966/// retrieve the list of bin boundaries
2967
2969{
2970 return this->getFunc()->binBoundaries(obs, xlo, xhi);
2971}
2972
2973////////////////////////////////////////////////////////////////////////////////
2974/// retrieve the sample Hint
2975
2977{
2978 return this->getFunc()->plotSamplingHint(obs, xlo, xhi);
2979}
2980
2981////////////////////////////////////////////////////////////////////////////////
2982/// call getVal on the internal function
2983
2985{
2986 // cout << "XX RooLagrangianMorphFunc::getValV(" << this << ") set = " << set
2987 // << std::endl ;
2988 this->_curNormSet = set;
2989 return RooAbsReal::getValV(set);
2990}
2991
2992////////////////////////////////////////////////////////////////////////////////
2993/// call getVal on the internal function
2994
2996{
2997 // call getVal on the internal function
2998 RooRealSumFunc *pdf = this->getFunc();
2999 if (pdf)
3000 return this->_scale * pdf->getVal(_curNormSet);
3001 else
3002 std::cerr << "unable to acquire in-built function!" << std::endl;
3003 return 0.;
3004}
3005
3006////////////////////////////////////////////////////////////////////////////////
3007/// check if this PDF is a binned distribution in the given observable
3008
3010{
3011 return this->getFunc()->isBinnedDistribution(obs);
3012}
3013
3014////////////////////////////////////////////////////////////////////////////////
3015/// check if observable exists in the RooArgSet (-?-)
3016
3018{
3019 return this->getFunc()->checkObservables(nset);
3020}
3021
3022////////////////////////////////////////////////////////////////////////////////
3023/// Force analytical integration for the given observable
3024
3026{
3027 return this->getFunc()->forceAnalyticalInt(arg);
3028}
3029
3030////////////////////////////////////////////////////////////////////////////////
3031/// Retrieve the mat
3032
3034 const char *rangeName) const
3035{
3036 return this->getFunc()->getAnalyticalIntegralWN(allVars, numVars, normSet, rangeName);
3037}
3038
3039////////////////////////////////////////////////////////////////////////////////
3040/// Retrieve the matrix of coefficients
3041
3042Double_t RooLagrangianMorphFunc::analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName) const
3043{
3044 return this->getFunc()->analyticalIntegralWN(code, normSet, rangeName);
3045}
3046
3047////////////////////////////////////////////////////////////////////////////////
3048/// Retrieve the matrix of coefficients
3049
3050void RooLagrangianMorphFunc::printMetaArgs(std::ostream &os) const
3051{
3052 return this->getFunc()->printMetaArgs(os);
3053}
3054
3055////////////////////////////////////////////////////////////////////////////////
3056/// Retrieve the matrix of coefficients
3057
3059{
3060 return this->getFunc()->canNodeBeCached();
3061}
3062
3063////////////////////////////////////////////////////////////////////////////////
3064/// Retrieve the matrix of coefficients
3065
3067{
3068 this->getFunc()->setCacheAndTrackHints(arg);
3069}
3070
3071////////////////////////////////////////////////////////////////////////////////
3072/// Retrieve the matrix of coefficients
3073
3075{
3076 auto cache = getCache(_curNormSet);
3077 if (!cache)
3078 coutE(Caching) << "unable to retrieve cache!" << std::endl;
3079 return makeRootMatrix(cache->_matrix);
3080}
3081
3082////////////////////////////////////////////////////////////////////////////////
3083/// Retrieve the matrix of coefficients after inversion
3084
3086{
3087 auto cache = getCache(_curNormSet);
3088 if (!cache)
3089 coutE(Caching) << "unable to retrieve cache!" << std::endl;
3090 return makeRootMatrix(cache->_inverse);
3091}
3092
3093////////////////////////////////////////////////////////////////////////////////
3094/// Retrieve the condition of the coefficient matrix. If the condition number
3095/// is very large, then the matrix is ill-conditioned and is almost singular.
3096/// The computation of the inverse is prone to large numerical errors
3097
3099{
3100 auto cache = getCache(_curNormSet);
3101 if (!cache)
3102 coutE(Caching) << "unable to retrieve cache!" << std::endl;
3103 return cache->_condition;
3104}
3105
3106////////////////////////////////////////////////////////////////////////////////
3107/// Return the RooRatio form of products and denominators of morphing functions
3108
3109std::unique_ptr<RooRatio>
3110RooLagrangianMorphFunc::makeRatio(const char *name, const char *title, RooArgList &nr, RooArgList &dr)
3111{
3112 RooArgList num, denom;
3113 for (auto it : nr) {
3114 num.add(*it);
3115 }
3116 for (auto it : dr) {
3117 denom.add(*it);
3118 }
3119 // same for denom
3120 return make_unique<RooRatio>(name, title, num, denom);
3121}
void Class()
Definition: Class.C:29
double
Definition: Converters.cxx:939
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
static constexpr double morphUnityDeviation
Matrix makeSuperMatrix(const TMatrixD &in)
convert a TMatrixD into a Matrix
static constexpr double morphLargestWeight
void writeMatrixToFileT(const MatrixT &matrix, const char *fname)
write a matrix to a text file
double invertMatrix(const Matrix &matrix, Matrix &inverse)
#define NaN
TMatrixD makeRootMatrix(const Matrix &in)
convert a matrix into a TMatrixD
Matrix diagMatrix(size_t n)
create a new diagonal matrix of size n
void printMatrix(const TMatrixD &mat)
write a matrix
#define oocxcoutW(o, a)
Definition: RooMsgService.h:95
#define coutW(a)
Definition: RooMsgService.h:32
#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
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:364
#define gDirectory
Definition: TDirectory.h:348
char name[80]
Definition: TGX11.cxx:110
int type
Definition: TGX11.cxx:121
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:22
char * Form(const char *fmt,...)
std::string & operator+=(std::string &left, const TString &right)
Definition: TString.h:478
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:75
void setStringAttribute(const Text_t *key, const Text_t *value)
Associate string 'value' to this object under key 'key'.
Definition: RooAbsArg.cxx:316
virtual void Print(Option_t *options=0) const
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:343
friend class RooListProxy
Definition: RooAbsArg.h:666
const RefCountList_t & servers() const
List of all servers of this object.
Definition: RooAbsArg.h:202
void setValueDirty()
Mark the element dirty. This forces a re-evaluation when a value is requested.
Definition: RooAbsArg.h:511
Bool_t getAttribute(const Text_t *name) const
Check if a named attribute is set. By default, all attributes are unset.
Definition: RooAbsArg.cxx:307
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:417
Bool_t isConstant() const
Check if the "Constant" attribute is set.
Definition: RooAbsArg.h:383
virtual Double_t * array() const =0
RooAbsCacheElement is the abstract base class for objects to be stored in RooAbsCache cache manager o...
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 * first() const
RooAbsArg * find(const char *name) const
Find object with given name in list.
RooAbsRealLValue is the common abstract base class for objects that represent a real value that may a...
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:63
Double_t getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:93
virtual Double_t getValV(const RooArgSet *normalisationSet=nullptr) const
Return value of object.
Definition: RooAbsReal.cxx:276
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:45
double weight(std::size_t i) const
Return weight of i-th bin.
Definition: RooDataHist.h:112
double weightSquared(std::size_t i) const
Return squared weight sum of i-th bin.
Definition: RooDataHist.h:116
const RooArgSet * get() const override
Get bin centre of current bin.
Definition: RooDataHist.h:84
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
RooArgSet const & getHistObsList() const
Definition: RooHistFunc.h:85
RooDataHist & dataHist()
Return RooDataHist that is represented.
Definition: RooHistFunc.h:40
virtual RooArgList containedArgs(Action) override
retrieve the list of contained args
static RooLagrangianMorphFunc::CacheElem * createCache(const RooLagrangianMorphFunc *func)
create all the temporary objects required by the class
void buildMatrix(const RooLagrangianMorphFunc::ParamMap &inputParameters, const RooLagrangianMorphFunc::FlagMap &inputFlags, const List &flags)
build and invert the morphing matrix
static RooLagrangianMorphFunc::CacheElem * createCache(const RooLagrangianMorphFunc *func, const Matrix &inverse)
create all the temporary objects required by the class function variant with precomputed inverse matr...
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
std::unique_ptr< RooRealSumFunc > _sumFunc
virtual void operModeHook(RooAbsArg::OperMode) override
Interface for changes of operation mode.
void buildMorphingFunction(const char *name, const RooLagrangianMorphFunc::ParamMap &inputParameters, const std::map< std::string, int > &storage, const RooArgList &physics, bool allowNegativeYields, RooRealVar *observable, RooRealVar *binWidth)
build the final morphing function
Class RooLagrangianMorphing is a implementation of the method of Effective Lagrangian Morphing,...
bool isParameterConstant(const char *paramname) const
return true if the parameter with the given name is set constant, false otherwise
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...
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.
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
void addFolders(const RooArgList &folders)
The cache manager.
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)
read the parameters from the input file
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 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
return extended mored capabilities
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::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.
std::vector< std::vector< RooListProxy * > > _diagrams
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(RooFit::SuperFloat c, RooAbsReal *t)
void setCoefficient(size_t idx, RooFit::SuperFloat c)
RooFit::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:486
void setBins(Int_t nBins, const char *name=0)
Create a uniform binning under name 'name' for this variable.
Definition: RooRealVar.cxx:433
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:516
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:343
void setBinning(const RooAbsBinning &binning, const char *name=0)
Add given binning under name 'name' with this variable.
Definition: RooRealVar.cxx:440
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:282
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
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
virtual TObject * Get(const char *namecycle)
Return pointer to object identified by namecycle.
Definition: TDirectory.cxx:814
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:4015
A TFolder object is a collection of objects and folders.
Definition: TFolder.h:30
TCollection * GetListOfFolders() const
Definition: TFolder.h:55
virtual void SetOwner(Bool_t owner=kTRUE)
Set ownership.
Definition: TFolder.cxx:489
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:9084
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:9100
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition: TH1.cxx:9030
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:5028
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition: TH1.cxx:9041
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 * GetTitle() const
Returns title of object.
Definition: TNamed.h:48
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
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
Class used by TMap to store (key,value) pairs.
Definition: TMap.h:102
TObject * Value() const
Definition: TMap.h:121
TObject * Key() const
Definition: TMap.h:120
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
const char * Data() const
Definition: TString.h:369
std::string_view View() const
Definition: TString.h:444
TString & Append(const char *cs)
Definition: TString.h:564
TLine * line
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< Sqrt< T >, VecExpr< A, T, D >, T >, T, D > sqrt(const VecExpr< A, T, D > &rhs)
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
double SuperFloat
Definition: Floats.h:29
@ 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
std::vector< std::string > folderNames
Definition: civetweb.c:1856
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12
REAL * vertex
Definition: triangle.c:513