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{
456 theFolder->SetOwner();
457 // And also need to set up ownership for nested folders
458 auto subdirs = theFolder->GetListOfFolders();
459 for (auto *subdir : *subdirs) {
460 auto thisfolder = dynamic_cast<TFolder *>(subdir);
461 if (thisfolder) {
462 // no explicit deletion here, will be handled by parent
463 setOwnerRecursive(thisfolder);
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{
478 std::unique_ptr<TFolder> theFolder(inFile->Get<TFolder>(folderName.c_str()));
479 if (!theFolder) {
480 std::cerr << "Error: unable to access data from folder '" << folderName << "' from file '" << inFile->GetName()
481 << "'!" << std::endl;
482 return nullptr;
483 }
484 setOwnerRecursive(theFolder.get());
485 return theFolder;
486}
487
488///////////////////////////////////////////////////////////////////////////////
489/// Helper to load a single object from a file-resident TFolder, while
490/// avoiding memory leaks.
491/// @tparam AObjType Type of object to load.
492/// @param inFile input file to load from. Expected to be a valid pointer
493/// @param folderName Name of the TFolder to load from the file
494/// @param objName Name of the object to load
495/// @param notFoundError If set, print a detailed error if we didn't find something
496/// @return Returns a pointer to a clone of the loaded object. Ownership assigned to the caller.
497template <class AObjType>
498AObjType *loadFromFileResidentFolder(TDirectory *inFile, const std::string &folderName, const std::string &objName,
499 bool notFoundError = true)
500{
501 auto folder = readOwningFolderFromFile(inFile, folderName);
502 if (!folder) {
503 return nullptr;
504 }
505 AObjType *loadedObject = dynamic_cast<AObjType *>(folder->FindObject(objName.c_str()));
506 if (!loadedObject) {
507 if (notFoundError) {
508 std::stringstream errstr;
509 errstr << "Error: unable to retrieve object '" << objName << "' from folder '" << folderName
510 << "'. contents are:";
511 TIter next(folder->GetListOfFolders()->begin());
512 TFolder *f;
513 while ((f = (TFolder *)next())) {
514 errstr << " " << f->GetName();
515 }
516 std::cerr << errstr.str() << std::endl;
517 }
518 return nullptr;
519 }
520 // replace the loaded object by a clone, as the loaded folder will delete the original
521 return static_cast<AObjType *>(
522 loadedObject->Clone()); // can use a static_cast - confirmed validity by initial cast above.
523}
524
525///////////////////////////////////////////////////////////////////////////////
526/// retrieve a ParamSet from a certain subfolder 'name' of the file
527
528template <class T>
529void readValues(std::map<const std::string, T> &myMap, TDirectory *file, const std::string &name,
530 const std::string &key = "param_card", bool notFoundError = true)
531{
532 TH1F *h_pc = loadFromFileResidentFolder<TH1F>(file, name, key, notFoundError);
533 readValues(myMap, h_pc);
534}
535
536///////////////////////////////////////////////////////////////////////////////
537/// retrieve the param_hists file and return a map of the parameter values
538/// by providing a list of names, only the param_hists of those subfolders are
539/// read leaving the list empty is interpreted as meaning 'read everything'
540
541template <class T>
542void readValues(std::map<const std::string, std::map<const std::string, T>> &inputParameters, TDirectory *f,
543 const std::vector<std::string> &names, const std::string &key = "param_card", bool notFoundError = true)
544{
545 inputParameters.clear();
546 // if the list of names is empty, we assume that this means 'all'
547 // loop over all folders in the file
548 for (size_t i = 0; i < names.size(); i++) {
549 const std::string name(names[i]);
550 // actually read an individual param_hist
551 readValues(inputParameters[name], f, name, key, notFoundError);
552 }
553
554 // now the map is filled with all parameter values found for all samples
555}
556
557///////////////////////////////////////////////////////////////////////////////
558/// open the file and return a file pointer
559
560inline TDirectory *openFile(const std::string &filename)
561{
562 if (filename.empty()) {
563 return gDirectory;
564 } else {
565 TFile *file = TFile::Open(filename.c_str(), "READ");
566 if (!file || !file->IsOpen()) {
567 if (file)
568 delete file;
569 std::cerr << "could not open file '" << filename << "'!" << std::endl;
570 }
571 return file;
572 }
573}
574
575///////////////////////////////////////////////////////////////////////////////
576/// open the file and return a file pointer
577
578inline void closeFile(TDirectory *d)
579{
580 TFile *f = dynamic_cast<TFile *>(d);
581 if (f) {
582 f->Close();
583 delete f;
584 d = nullptr;
585 }
586}
587
588///////////////////////////////////////////////////////////////////////////////
589/// extract the operators from a single coupling
590template <class T2>
591inline void extractServers(const RooAbsArg &coupling, T2 &operators)
592{
593 int nservers = 0;
594 for (const auto server : coupling.servers()) {
595 extractServers(*server, operators);
596 nservers++;
597 }
598 if (nservers == 0) {
599 operators.add(coupling);
600 }
601}
602
603///////////////////////////////////////////////////////////////////////////////
604/// extract the operators from a list of couplings
605
607inline void extractOperators(const T1 &couplings, T2 &operators)
608{
609 // coutD(InputArguments) << "extracting operators from
610 // "<<couplings.getSize()<<" couplings" << std::endl;
611 for (auto itr : couplings) {
612 extractServers(*itr, operators);
613 }
614}
615
616///////////////////////////////////////////////////////////////////////////////
617/// extract the operators from a list of vertices
618
620inline void extractOperators(const T1 &vec, T2 &operators)
621{
622 for (const auto &v : vec) {
623 extractOperators(v, operators);
624 }
625}
626
627///////////////////////////////////////////////////////////////////////////////
628/// extract the couplings from a given set and copy them to a new one
629
630template <class T1, class T2>
631inline void extractCouplings(const T1 &inCouplings, T2 &outCouplings)
632{
633 for (auto itr : inCouplings) {
634 if (!outCouplings.find(itr->GetName())) {
635 // coutD(InputArguments) << "adding parameter "<< obj->GetName() <<
636 // std::endl;
637 outCouplings.add(*itr);
638 }
639 }
640}
641
642///////////////////////////////////////////////////////////////////////////////
643/// find and, if necessary, create a parameter from a list
644
645template <class T>
646inline RooAbsArg &get(T &operators, const char *name, double defaultval = 0)
647{
648 RooAbsArg *kappa = operators.find(name);
649 if (kappa)
650 return *kappa;
651 RooRealVar *newKappa = new RooRealVar(name, name, defaultval);
652 newKappa->setConstant(false);
653 operators.add(*newKappa);
654 return *newKappa;
655}
656
657///////////////////////////////////////////////////////////////////////////////
658/// find and, if necessary, create a parameter from a list
659
660template <class T>
661inline RooAbsArg &get(T &operators, const std::string &name, double defaultval = 0)
662{
663 return get(operators, name.c_str(), defaultval);
664}
665
666////////////////////////////////////////////////////////////////////////////////
667/// create a new coupling and add it to the set
668
669template <class T>
670inline void addCoupling(T &set, const TString &name, const TString &formula, const RooArgList &components, bool isNP)
671{
672 if (!set.find(name)) {
673 RooFormulaVar *c = new RooFormulaVar(name, formula, components);
674 c->setAttribute("NewPhysics", isNP);
675 set.add(*c);
676 }
677}
678
679////////////////////////////////////////////////////////////////////////////////
680/// set parameter values first set all values to defaultVal (if value not
681/// present in param_card then it should be 0)
682
683inline bool setParam(RooRealVar *p, double val, bool force)
684{
685 bool ok = true;
686 if (val > p->getMax()) {
687 if (force) {
688 p->setMax(val);
689 } else {
690 std::cerr << ": parameter " << p->GetName() << " out of bounds: " << val << " > " << p->getMax() << std::endl;
691 ok = false;
692 }
693 } else if (val < p->getMin()) {
694 if (force) {
695 p->setMin(val);
696 } else {
697 std::cerr << ": parameter " << p->GetName() << " out of bounds: " << val << " < " << p->getMin() << std::endl;
698 ok = false;
699 }
700 }
701 if (ok)
702 p->setVal(val);
703 return ok;
704}
705
706////////////////////////////////////////////////////////////////////////////////
707/// set parameter values first set all values to defaultVal (if value not
708/// present in param_card then it should be 0)
709
710template <class T1, class T2>
711inline bool setParams(const T2 &args, T1 val)
712{
713 for (auto itr : args) {
714 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
715 if (!param)
716 continue;
717 setParam(param, val, true);
718 }
719 return true;
720}
721
722////////////////////////////////////////////////////////////////////////////////
723/// set parameter values first set all values to defaultVal (if value not
724/// present in param_card then it should be 0)
725
726template <class T1, class T2>
727inline bool
728setParams(const std::map<const std::string, T1> &point, const T2 &args, bool force = false, T1 defaultVal = 0)
729{
730 bool ok = true;
731 for (auto itr : args) {
732 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
733 if (!param || param->isConstant())
734 continue;
735 ok = setParam(param, defaultVal, force) && ok;
736 }
737 // set all parameters to the values in the param_card histogram
738 for (auto paramit : point) {
739 // loop over all the parameters
740 const std::string param(paramit.first);
741 // retrieve them from the map
742 RooRealVar *p = dynamic_cast<RooRealVar *>(args.find(param.c_str()));
743 if (!p)
744 continue;
745 // set them to their nominal value
746 ok = setParam(p, paramit.second, force) && ok;
747 }
748 return ok;
749}
750
751////////////////////////////////////////////////////////////////////////////////
752/// set parameter values first set all values to defaultVal (if value not
753/// present in param_card then it should be 0)
754
755template <class T>
756inline bool setParams(TH1 *hist, const T &args, bool force = false)
757{
758 bool ok = true;
759
760 for (auto itr : args) {
761 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
762 if (!param)
763 continue;
764 ok = setParam(param, 0., force) && ok;
765 }
766
767 // set all parameters to the values in the param_card histogram
768 TAxis *ax = hist->GetXaxis();
769 for (int i = 1; i <= ax->GetNbins(); ++i) {
770 // loop over all the parameters
771 RooRealVar *p = dynamic_cast<RooRealVar *>(args.find(ax->GetBinLabel(i)));
772 if (!p)
773 continue;
774 // set them to their nominal value
775 ok = setParam(p, hist->GetBinContent(i), force) && ok;
776 }
777 return ok;
778}
779
780////////////////////////////////////////////////////////////////////////////////
781/// create a set of parameters
782
783template <class T>
784inline RooLagrangianMorphFunc::ParamSet getParams(const T &parameters)
785{
787 for (auto itr : parameters) {
788 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
789 if (!param)
790 continue;
791 retval[param->GetName()] = param->getVal();
792 }
793 return retval;
794}
795
796////////////////////////////////////////////////////////////////////////////////
797/// collect the histograms from the input file and convert them to RooFit
798/// objects
799
800void collectHistograms(const char *name, TDirectory *file, std::map<std::string, int> &list_hf, RooArgList &physics,
801 RooRealVar &var, const std::string &varname,
802 const RooLagrangianMorphFunc::ParamMap &inputParameters)
803{
804 bool binningOK = false;
805 for (auto sampleit : inputParameters) {
806 const std::string sample(sampleit.first);
807 TH1 *hist = loadFromFileResidentFolder<TH1>(file, sample, varname, true);
808 if (!hist)
809 return;
810
811 auto it = list_hf.find(sample);
812 if (it != list_hf.end()) {
813 RooHistFunc *hf = (RooHistFunc *)(physics.at(it->second));
814 hf->setValueDirty();
815 // commenting out To-be-resolved
816 // RooDataHist* dh = &(hf->dataHist());
817 // RooLagrangianMorphFunc::setDataHistogram(hist,&var,dh);
818 // RooArgSet vars;
819 // vars.add(var);
820 // dh->importTH1(vars,*hist,1.,false);
821 } else {
822 if (!binningOK) {
823 int n = hist->GetNbinsX();
824 std::vector<double> bins;
825 for (int i = 1; i < n + 1; ++i) {
826 bins.push_back(hist->GetBinLowEdge(i));
827 }
828 bins.push_back(hist->GetBinLowEdge(n) + hist->GetBinWidth(n));
829 var.setBinning(RooBinning(n, &(bins[0])));
830 }
831
832 // generate the mean value
833 TString histname = makeValidName(Form("dh_%s_%s", sample.c_str(), name));
834 TString funcname = makeValidName(Form("phys_%s_%s", sample.c_str(), name));
835 RooArgSet vars;
836 vars.add(var);
837
838 RooDataHist *dh = new RooDataHist(histname, histname, vars, hist);
839 // add it to the list
840 RooHistFunc *hf = new RooHistFunc(funcname, funcname, var, *dh);
841 int idx = physics.getSize();
842 list_hf[sample] = idx;
843 physics.add(*hf);
844 assert(hf = static_cast<RooHistFunc *>(physics.at(idx)));
845 }
846 // std::cout << "found histogram " << hist->GetName() << " with integral "
847 // << hist->Integral() << std::endl;
848 }
849}
850
851////////////////////////////////////////////////////////////////////////////////
852/// collect the RooAbsReal objects from the input directory
853
854void collectRooAbsReal(const char * /*name*/, TDirectory *file, std::map<std::string, int> &list_hf,
855 RooArgList &physics, const std::string &varname,
856 const RooLagrangianMorphFunc::ParamMap &inputParameters)
857{
858 for (auto sampleit : inputParameters) {
859 const std::string sample(sampleit.first);
860 RooAbsReal *obj = loadFromFileResidentFolder<RooAbsReal>(file, sample, varname, true);
861 if (!obj)
862 return;
863 auto it = list_hf.find(sample);
864 if (it == list_hf.end()) {
865 int idx = physics.getSize();
866 list_hf[sample] = idx;
867 physics.add(*obj);
868 assert(obj == physics.at(idx));
869 }
870 }
871}
872
873////////////////////////////////////////////////////////////////////////////////
874/// collect the TParameter objects from the input file and convert them to
875/// RooFit objects
876
877template <class T>
878void collectCrosssections(const char *name, TDirectory *file, std::map<std::string, int> &list_xs, RooArgList &physics,
879 const std::string &varname, const RooLagrangianMorphFunc::ParamMap &inputParameters)
880{
881 for (auto sampleit : inputParameters) {
882 const std::string sample(sampleit.first);
883 TObject *obj = loadFromFileResidentFolder<TObject>(file, sample, varname, false);
884 TParameter<T> *xsection = nullptr;
885 TParameter<T> *error = nullptr;
886 TParameter<T> *p = dynamic_cast<TParameter<T> *>(obj);
887 if (p) {
888 xsection = p;
889 }
890 TPair *pair = dynamic_cast<TPair *>(obj);
891 if (pair) {
892 xsection = dynamic_cast<TParameter<T> *>(pair->Key());
893 error = dynamic_cast<TParameter<T> *>(pair->Value());
894 }
895 if (!xsection) {
896 std::stringstream errstr;
897 errstr << "Error: unable to retrieve cross section '" << varname << "' from folder '" << sample;
898 return;
899 }
900
901 auto it = list_xs.find(sample.c_str());
902 RooRealVar *xs;
903 if (it != list_xs.end()) {
904 xs = (RooRealVar *)(physics.at(it->second));
905 xs->setVal(xsection->GetVal());
906 } else {
907 std::string objname = Form("phys_%s_%s", name, sample.c_str());
908 xs = new RooRealVar(objname.c_str(), objname.c_str(), xsection->GetVal());
909 xs->setConstant(true);
910 int idx = physics.getSize();
911 list_xs[sample] = idx;
912 physics.add(*xs);
913 assert(physics.at(idx) == xs);
914 }
915 if (error)
916 xs->setError(error->GetVal());
917 }
918}
919
920////////////////////////////////////////////////////////////////////////////////
921/// collect the TPair<TParameter,TParameter> objects from the input file and
922/// convert them to RooFit objects
923
924void collectCrosssectionsTPair(const char *name, TDirectory *file, std::map<std::string, int> &list_xs,
925 RooArgList &physics, const std::string &varname, const std::string &basefolder,
926 const RooLagrangianMorphFunc::ParamMap &inputParameters)
927{
928 TPair *pair = loadFromFileResidentFolder<TPair>(file, basefolder, varname, false);
929 if (!pair)
930 return;
931 TParameter<double> *xsec_double = dynamic_cast<TParameter<double> *>(pair->Key());
932 if (xsec_double) {
933 collectCrosssections<double>(name, file, list_xs, physics, varname, inputParameters);
934 } else {
935 TParameter<float> *xsec_float = dynamic_cast<TParameter<float> *>(pair->Key());
936 if (xsec_float) {
937 collectCrosssections<float>(name, file, list_xs, physics, varname, inputParameters);
938 } else {
939 std::cerr << "cannot morph objects of class 'TPair' if parameter is not "
940 "double or float!"
941 << std::endl;
942 }
943 }
944}
945
946///////////////////////////////////////////////////////////////////////////////
947// FORMULA CALCULATION ////////////////////////////////////////////////////////
948///////////////////////////////////////////////////////////////////////////////
949///////////////////////////////////////////////////////////////////////////////
950
951////////////////////////////////////////////////////////////////////////////////
952/// recursive function to determine polynomials
953
954void collectPolynomialsHelper(const FeynmanDiagram &diagram, MorphFuncPattern &morphfunc, std::vector<int> &term,
955 int vertexid, bool first)
956{
957 if (vertexid > 0) {
958 for (size_t i = 0; i < diagram[vertexid - 1].size(); ++i) {
959 if (!diagram[vertexid - 1][i])
960 continue;
961 std::vector<int> newterm(term);
962 newterm[i]++;
963 if (first) {
964 ::collectPolynomialsHelper(diagram, morphfunc, newterm, vertexid, false);
965 } else {
966 ::collectPolynomialsHelper(diagram, morphfunc, newterm, vertexid - 1, true);
967 }
968 }
969 } else {
970 bool found = false;
971 for (size_t i = 0; i < morphfunc.size(); ++i) {
972 bool thisfound = true;
973 for (size_t j = 0; j < morphfunc[i].size(); ++j) {
974 if (morphfunc[i][j] != term[j]) {
975 thisfound = false;
976 break;
977 }
978 }
979 if (thisfound) {
980 found = true;
981 break;
982 }
983 }
984 if (!found) {
985 morphfunc.push_back(term);
986 }
987 }
988}
989
990////////////////////////////////////////////////////////////////////////////////
991/// calculate the morphing function pattern based on a vertex map
992
993void collectPolynomials(MorphFuncPattern &morphfunc, const FeynmanDiagram &diagram)
994{
995 int nvtx(diagram.size());
996 std::vector<int> term(diagram[0].size(), 0);
997
998 ::collectPolynomialsHelper(diagram, morphfunc, term, nvtx, true);
999}
1000
1001////////////////////////////////////////////////////////////////////////////////
1002/// build a vertex map based on vertices and couplings appearing
1003
1004template <class List>
1005inline void fillFeynmanDiagram(FeynmanDiagram &diagram, const std::vector<List *> &vertices, RooArgList &couplings)
1006{
1007 const int ncouplings = couplings.getSize();
1008 // std::cout << "Number of couplings " << ncouplings << std::endl;
1009 for (auto const &vertex : vertices) {
1010 std::vector<bool> vertexCouplings(ncouplings, false);
1011 int idx = -1;
1012 RooAbsReal *coupling;
1013 for (auto citr : couplings) {
1014 coupling = dynamic_cast<RooAbsReal *>(citr);
1015 idx++;
1016 if (!coupling) {
1017 std::cerr << "encountered invalid list of couplings in vertex!" << std::endl;
1018 return;
1019 }
1020 if (vertex->find(coupling->GetName())) {
1021 vertexCouplings[idx] = true;
1022 }
1023 }
1024 diagram.push_back(vertexCouplings);
1025 }
1026}
1027
1028////////////////////////////////////////////////////////////////////////////////
1029/// fill the matrix of coefficients
1030
1031template <class MatrixT, class T1, class T2>
1032inline MatrixT buildMatrixT(const RooLagrangianMorphFunc::ParamMap &inputParameters, const FormulaList &formulas,
1033 const T1 &args, const RooLagrangianMorphFunc::FlagMap &flagValues, const T2 &flags)
1034{
1035 const size_t dim = inputParameters.size();
1036 MatrixT matrix(dim, dim);
1037 int row = 0;
1038 for (auto sampleit : inputParameters) {
1039 const std::string sample(sampleit.first);
1040 // set all vars to value stored in input file
1041 if (!setParams<double>(sampleit.second, args, true, 0)) {
1042 std::cout << "unable to set parameters for sample " << sample << "!" << std::endl;
1043 }
1044 auto flagit = flagValues.find(sample);
1045 if (flagit != flagValues.end() && !setParams<int>(flagit->second, flags, true, 1)) {
1046 std::cout << "unable to set parameters for sample " << sample << "!" << std::endl;
1047 }
1048 // loop over all the formulas
1049 int col = 0;
1050 for (auto const &formula : formulas) {
1051 if (!formula.second) {
1052 std::cerr << "Error: invalid formula encountered!" << std::endl;
1053 }
1054 matrix(row, col) = formula.second->getVal();
1055 col++;
1056 }
1057 row++;
1058 }
1059 return matrix;
1060}
1061
1062////////////////////////////////////////////////////////////////////////////////
1063/// check if the matrix is square
1064
1065inline void checkMatrix(const RooLagrangianMorphFunc::ParamMap &inputParameters, const FormulaList &formulas)
1066{
1067 if (inputParameters.size() != formulas.size()) {
1068 std::stringstream ss;
1069 ss << "matrix is not square, consistency check failed: " << inputParameters.size() << " samples, "
1070 << formulas.size() << " expressions:" << std::endl;
1071 ss << "formulas: " << std::endl;
1072 for (auto const &formula : formulas) {
1073 ss << formula.second->GetTitle() << std::endl;
1074 }
1075 ss << "samples: " << std::endl;
1076 for (auto sample : inputParameters) {
1077 ss << sample.first << std::endl;
1078 }
1079 std::cerr << ss.str() << std::endl;
1080 }
1081}
1082
1083////////////////////////////////////////////////////////////////////////////////
1084/// check if the entries in the inverted matrix are sensible
1085
1086inline void inverseSanity(const Matrix &matrix, const Matrix &inverse, double &unityDeviation, double &largestWeight)
1087{
1088 Matrix unity(inverse * matrix);
1089
1090 unityDeviation = 0.;
1091 largestWeight = 0.;
1092 const size_t dim = size(unity);
1093 for (size_t i = 0; i < dim; ++i) {
1094 for (size_t j = 0; j < dim; ++j) {
1095 if (inverse(i, j) > largestWeight) {
1096 largestWeight = (double)inverse(i, j);
1097 }
1098 if (fabs(unity(i, j) - static_cast<int>(i == j)) > unityDeviation) {
1099 unityDeviation = fabs((double)unity(i, j)) - static_cast<int>(i == j);
1100 }
1101 }
1102 }
1103}
1104
1105////////////////////////////////////////////////////////////////////////////////
1106/// check for name conflicts between the input samples and an argument set
1107template <class List>
1108inline void checkNameConflict(const RooLagrangianMorphFunc::ParamMap &inputParameters, List &args)
1109{
1110 for (auto sampleit : inputParameters) {
1111 const std::string sample(sampleit.first);
1112 RooAbsArg *arg = args.find(sample.c_str());
1113 if (arg) {
1114 std::cerr << "detected name conflict: cannot use sample '" << sample
1115 << "' - a parameter with the same name of type '" << arg->ClassName() << "' is present in set '"
1116 << args.GetName() << "'!" << std::endl;
1117 }
1118 }
1119}
1120
1121////////////////////////////////////////////////////////////////////////////////
1122/// build the formulas corresponding to the given set of input files and
1123/// the physics process
1124
1125FormulaList buildFormulas(const char *mfname, const RooLagrangianMorphFunc::ParamMap &inputParameters,
1126 const RooLagrangianMorphFunc::FlagMap &inputFlags, const MorphFuncPattern &morphfunc,
1127 const RooArgList &couplings, const RooArgList &flags,
1128 const std::vector<RooArgList *> &nonInterfering)
1129{
1130 // example vbf hww:
1131 // Operators kSM, kHww, kAww, kHdwR,kHzz, kAzz
1132 // std::vector<bool> vertexProd = {true, true, true, true, true, true };
1133 // std::vector<bool> vertexDecay = {true, true, true, true, false,false};
1134 // diagram.push_back(vertexProd);
1135 // diagram.push_back(vertexDecay);
1136
1137 const int ncouplings = couplings.getSize();
1138 std::vector<bool> couplingsZero(ncouplings, true);
1139 std::map<TString, bool> flagsZero;
1140
1141 RooArgList operators;
1142 extractOperators(couplings, operators);
1143 size_t nOps = operators.getSize();
1144
1145 for (auto sampleit : inputParameters) {
1146 const std::string sample(sampleit.first);
1147 if (!setParams(sampleit.second, operators, true)) {
1148 std::cerr << "unable to set parameters for sample '" << sample << "'!" << std::endl;
1149 }
1150
1151 if ((int)nOps != (operators.getSize())) {
1152 std::cerr << "internal error, number of operators inconsistent!" << std::endl;
1153 }
1154
1155 RooAbsReal *obj0;
1156 int idx = 0;
1157
1158 for (auto itr1 : couplings) {
1159 obj0 = dynamic_cast<RooAbsReal *>(itr1);
1160 if (obj0->getVal() != 0) {
1161 couplingsZero[idx] = false;
1162 }
1163 idx++;
1164 }
1165 }
1166
1167 for (auto itr2 : flags) {
1168 auto obj1 = dynamic_cast<RooAbsReal *>(itr2);
1169 int nZero = 0;
1170 int nNonZero = 0;
1171 for (auto sampleit : inputFlags) {
1172 const auto &flag = sampleit.second.find(obj1->GetName());
1173 if (flag != sampleit.second.end()) {
1174 if (flag->second == 0.)
1175 nZero++;
1176 else
1177 nNonZero++;
1178 }
1179 }
1180 if (nZero > 0 && nNonZero == 0)
1181 flagsZero[obj1->GetName()] = true;
1182 else
1183 flagsZero[obj1->GetName()] = false;
1184 }
1185
1186 FormulaList formulas;
1187 for (size_t i = 0; i < morphfunc.size(); ++i) {
1188 RooArgList ss;
1189 bool isZero = false;
1190 std::string reason;
1191 // check if this is a blacklisted interference term
1192 for (const auto &group : nonInterfering) {
1193 int nInterferingOperators = 0;
1194 for (size_t j = 0; j < morphfunc[i].size(); ++j) {
1195 if (morphfunc[i][j] % 2 == 0)
1196 continue; // even exponents are not interference terms
1197 if (group->find(couplings.at(j)->GetName())) { // if the coupling is part of a
1198 // "pairwise non-interfering group"
1199 nInterferingOperators++;
1200 }
1201 }
1202 if (nInterferingOperators > 1) {
1203 isZero = true;
1204 reason = "blacklisted interference term!";
1205 }
1206 }
1207 int nNP = 0;
1208 if (!isZero) {
1209 // prepare the term
1210 for (size_t j = 0; j < morphfunc[i].size(); ++j) {
1211 const int exponent = morphfunc[i][j];
1212 if (exponent == 0)
1213 continue;
1214 RooAbsReal *coupling = dynamic_cast<RooAbsReal *>(couplings.at(j));
1215 for (int k = 0; k < exponent; ++k) {
1216 ss.add(*coupling);
1217 if (coupling->getAttribute("NewPhysics")) {
1218 nNP++;
1219 }
1220 }
1221 std::string cname(coupling->GetName());
1222 if (coupling->getAttribute("LO") && exponent > 1) {
1223 isZero = true;
1224 reason = "coupling " + cname + " was listed as leading-order-only";
1225 }
1226 // mark the term as zero if any of the couplings are zero
1227 if (!isZero && couplingsZero[j]) {
1228 isZero = true;
1229 reason = "coupling " + cname + " is zero!";
1230 }
1231 }
1232 }
1233 // check and apply flags
1234 bool removedByFlag = false;
1235
1236 for (auto itr : flags) {
1237 auto obj = dynamic_cast<RooAbsReal *>(itr);
1238 if (!obj)
1239 continue;
1240 TString sval(obj->getStringAttribute("NewPhysics"));
1241 int val = atoi(sval);
1242 if (val == nNP) {
1243 if (flagsZero.find(obj->GetName()) != flagsZero.end() && flagsZero.at(obj->GetName())) {
1244 removedByFlag = true;
1245 reason = Form("flag %s is zero", obj->GetName());
1246 }
1247 ss.add(*obj);
1248 }
1249 }
1250
1251 // create and add the formula
1252 if (!isZero && !removedByFlag) {
1253 // build the name
1254 const auto name = std::string(mfname) + "_pol" + std::to_string(i);
1255 formulas[i] = std::make_unique<RooProduct>(name.c_str(), ::concatNames(ss, " * ").c_str(), ss);
1256 }
1257 }
1258 return formulas;
1259}
1260
1261////////////////////////////////////////////////////////////////////////////////
1262/// create the weight formulas required for the morphing
1263
1264FormulaList createFormulas(const char *name, const RooLagrangianMorphFunc::ParamMap &inputs,
1265 const RooLagrangianMorphFunc::FlagMap &inputFlags,
1266 const std::vector<std::vector<RooArgList *>> &diagrams, RooArgList &couplings,
1267 const RooArgList &flags, const std::vector<RooArgList *> &nonInterfering)
1268{
1269 MorphFuncPattern morphfuncpattern;
1270
1271 for (const auto &vertices : diagrams) {
1272 FeynmanDiagram d;
1273 ::fillFeynmanDiagram(d, vertices, couplings);
1274 ::collectPolynomials(morphfuncpattern, d);
1275 }
1276 FormulaList retval = buildFormulas(name, inputs, inputFlags, morphfuncpattern, couplings, flags, nonInterfering);
1277 if (retval.empty()) {
1278 std::stringstream errorMsgStream;
1279 errorMsgStream
1280 << "no formulas are non-zero, check if any if your couplings is floating and missing from your param_cards!"
1281 << std::endl;
1282 const auto errorMsg = errorMsgStream.str();
1283 throw std::runtime_error(errorMsg);
1284 }
1285 checkMatrix(inputs, retval);
1286 return retval;
1287}
1288
1289////////////////////////////////////////////////////////////////////////////////
1290/// build the sample weights required for the input templates
1291//
1292template <class T1>
1293inline void buildSampleWeights(T1 &weights, const char *fname, const RooLagrangianMorphFunc::ParamMap &inputParameters,
1294 FormulaList &formulas, const Matrix &inverse)
1295{
1296 int sampleidx = 0;
1297
1298 for (auto sampleit : inputParameters) {
1299 const std::string sample(sampleit.first);
1300 std::stringstream title;
1301 TString name_full(makeValidName(sample.c_str()));
1302 if (fname) {
1303 name_full.Append("_");
1304 name_full.Append(fname);
1305 name_full.Prepend("w_");
1306 }
1307
1308 int formulaidx = 0;
1309 // build the formula with the correct normalization
1310 RooLinearCombination *sampleformula = new RooLinearCombination(name_full.Data());
1311 for (auto const &formulait : formulas) {
1312 const RooFit::SuperFloat val(inverse(formulaidx, sampleidx));
1313 sampleformula->add(val, formulait.second.get());
1314 formulaidx++;
1315 }
1316 weights.add(*sampleformula);
1317 sampleidx++;
1318 }
1319}
1320
1321inline std::map<std::string, std::string>
1322buildSampleWeightStrings(const RooLagrangianMorphFunc::ParamMap &inputParameters, const FormulaList &formulas,
1323 const Matrix &inverse)
1324{
1325 int sampleidx = 0;
1326 std::map<std::string, std::string> weights;
1327 for (auto sampleit : inputParameters) {
1328 const std::string sample(sampleit.first);
1329 std::stringstream str;
1330 int formulaidx = 0;
1331 // build the formula with the correct normalization
1332 for (auto const &formulait : formulas) {
1333 double val(inverse(formulaidx, sampleidx));
1334 if (val != 0.) {
1335 if (formulaidx > 0 && val > 0)
1336 str << " + ";
1337 str << val << "*(" << formulait.second->GetTitle() << ")";
1338 }
1339 formulaidx++;
1340 }
1341 weights[sample] = str.str();
1342 sampleidx++;
1343 }
1344 return weights;
1345}
1346} // namespace
1347
1348///////////////////////////////////////////////////////////////////////////////
1349// CacheElem magic ////////////////////////////////////////////////////////////
1350///////////////////////////////////////////////////////////////////////////////
1351
1353public:
1354 std::unique_ptr<RooRealSumFunc> _sumFunc = nullptr;
1356
1357 FormulaList _formulas;
1359
1363
1366
1367 //////////////////////////////////////////////////////////////////////////////
1368 /// retrieve the list of contained args
1369
1371 {
1372 RooArgList args(*_sumFunc);
1373 args.add(_weights);
1374 args.add(_couplings);
1375 for (auto const &it : _formulas) {
1376 args.add(*(it.second));
1377 }
1378 return args;
1379 }
1380
1381 //////////////////////////////////////////////////////////////////////////////
1382 // default destructor
1383
1384 ~CacheElem() override {}
1385
1386 //////////////////////////////////////////////////////////////////////////////
1387 /// create the basic objects required for the morphing
1388
1389 inline void createComponents(const RooLagrangianMorphFunc::ParamMap &inputParameters,
1390 const RooLagrangianMorphFunc::FlagMap &inputFlags, const char *funcname,
1391 const std::vector<std::vector<RooListProxy *>> &diagramProxyList,
1392 const std::vector<RooArgList *> &nonInterfering, const RooArgList &flags)
1393 {
1394 RooArgList operators;
1395 std::vector<std::vector<RooArgList *>> diagrams;
1396 for (const auto &diagram : diagramProxyList) {
1397 diagrams.emplace_back();
1398 for (RooArgList *vertex : diagram) {
1399 extractCouplings(*vertex, _couplings);
1400 diagrams.back().emplace_back(vertex);
1401 }
1402 }
1403 extractOperators(_couplings, operators);
1404 _formulas =
1405 ::createFormulas(funcname, inputParameters, inputFlags, diagrams, _couplings, flags, nonInterfering);
1406 }
1407
1408 //////////////////////////////////////////////////////////////////////////////
1409 /// build and invert the morphing matrix
1410 template <class List>
1411 inline void buildMatrix(const RooLagrangianMorphFunc::ParamMap &inputParameters,
1412 const RooLagrangianMorphFunc::FlagMap &inputFlags, const List &flags)
1413 {
1414 RooArgList operators;
1415 extractOperators(_couplings, operators);
1416 Matrix matrix(buildMatrixT<Matrix>(inputParameters, _formulas, operators, inputFlags, flags));
1417 if (size(matrix) < 1) {
1418 std::cerr << "input matrix is empty, please provide suitable input samples!" << std::endl;
1419 }
1420 Matrix inverse(diagMatrix(size(matrix)));
1421
1422 double condition = (double)(invertMatrix(matrix, inverse));
1423 double unityDeviation, largestWeight;
1424 inverseSanity(matrix, inverse, unityDeviation, largestWeight);
1425 bool weightwarning(largestWeight > morphLargestWeight ? true : false);
1426 bool unitywarning(unityDeviation > morphUnityDeviation ? true : false);
1427
1428 if (false) {
1429 if (unitywarning) {
1430 oocxcoutW((TObject *)0, Eval) << "Warning: The matrix inversion seems to be unstable. This can "
1431 "be a result to input samples that are not sufficiently "
1432 "different to provide any morphing power."
1433 << std::endl;
1434 } else if (weightwarning) {
1435 oocxcoutW((TObject *)0, Eval) << "Warning: Some weights are excessively large. This can be a "
1436 "result to input samples that are not sufficiently different to "
1437 "provide any morphing power."
1438 << std::endl;
1439 }
1440 oocxcoutW((TObject *)0, Eval) << " Please consider the couplings "
1441 "encoded in your samples to cross-check:"
1442 << std::endl;
1443 for (auto sampleit : inputParameters) {
1444 const std::string sample(sampleit.first);
1445 oocxcoutW((TObject *)0, Eval) << " " << sample << ": ";
1446 // set all vars to value stored in input file
1447 setParams(sampleit.second, operators, true);
1448 bool first = true;
1449 RooAbsReal *obj;
1450
1451 for (auto itr : _couplings) {
1452 obj = dynamic_cast<RooAbsReal *>(itr);
1453 if (!first)
1454 std::cerr << ", ";
1455 oocxcoutW((TObject *)0, Eval) << obj->GetName() << "=" << obj->getVal();
1456 first = false;
1457 }
1458 oocxcoutW((TObject *)0, Eval) << std::endl;
1459 }
1460 }
1461#ifndef USE_UBLAS
1462 _matrix.ResizeTo(matrix.GetNrows(), matrix.GetNrows());
1463 _inverse.ResizeTo(matrix.GetNrows(), matrix.GetNrows());
1464#endif
1465 _matrix = matrix;
1466 _inverse = inverse;
1467 _condition = condition;
1468 }
1469
1470 ////////////////////////////////////////////////////////////////////////////////
1471 /// build the final morphing function
1472
1473 inline void buildMorphingFunction(const char *name, const RooLagrangianMorphFunc::ParamMap &inputParameters,
1474 const std::map<std::string, int> &storage, const RooArgList &physics,
1475 bool allowNegativeYields, RooRealVar *observable, RooRealVar *binWidth)
1476 {
1477 if (!binWidth) {
1478 std::cerr << "invalid bin width given!" << std::endl;
1479 return;
1480 }
1481 if (!observable) {
1482 std::cerr << "invalid observable given!" << std::endl;
1483 return;
1484 }
1485
1486 RooArgList operators;
1487 extractOperators(_couplings, operators);
1488
1489 // retrieve the weights
1490 ::buildSampleWeights(_weights, name, inputParameters, _formulas, _inverse);
1491
1492 // build the products of element and weight for each sample
1493 size_t i = 0;
1494 RooArgList sumElements;
1495 RooArgList scaleElements;
1496 for (auto sampleit : inputParameters) {
1497 // for now, we assume all the lists are nicely ordered
1498 TString prodname(makeValidName(sampleit.first.c_str()));
1499
1500 RooAbsReal *obj = static_cast<RooAbsReal *>(physics.at(storage.at(prodname.Data())));
1501
1502 if (!obj) {
1503 std::cerr << "unable to access physics object for " << prodname << std::endl;
1504 return;
1505 }
1506
1507 RooAbsReal *weight = static_cast<RooAbsReal *>(_weights.at(i));
1508
1509 if (!weight) {
1510 std::cerr << "unable to access weight object for " << prodname << std::endl;
1511 return;
1512 }
1513 prodname.Append("_");
1514 prodname.Append(name);
1515 RooArgList prodElems(*weight, *obj);
1516
1517 allowNegativeYields = true;
1518 auto prod = std::make_unique<RooProduct>(prodname, prodname, prodElems);
1519 if (!allowNegativeYields) {
1520 auto maxname = std::string(prodname) + "_max0";
1521 RooArgSet prodset(*prod);
1522
1523 auto max = std::make_unique<RooFormulaVar>(maxname.c_str(), "max(0," + prodname + ")", prodset);
1524 max->addOwnedComponents(std::move(prod));
1525 sumElements.addOwned(std::move(max));
1526 } else {
1527 sumElements.addOwned(std::move(prod));
1528 }
1529 scaleElements.add(*(binWidth));
1530 i++;
1531 }
1532
1533 // put everything together
1534 _sumFunc = make_unique<RooRealSumFunc>(Form("%s_morphfunc", name), name, sumElements, scaleElements);
1535
1536 if (!observable)
1537 std::cerr << "unable to access observable" << std::endl;
1538 _sumFunc.get()->addServer(*observable);
1539 if (!binWidth)
1540 std::cerr << "unable to access bin width" << std::endl;
1541 _sumFunc.get()->addServer(*binWidth);
1542 if (operators.getSize() < 1)
1543 std::cerr << "no operators listed" << std::endl;
1544 _sumFunc.get()->addServerList(operators);
1545 if (_weights.getSize() < 1)
1546 std::cerr << "unable to access weight objects" << std::endl;
1547 _sumFunc.get()->addOwnedComponents(_weights);
1548 _sumFunc.get()->addOwnedComponents(std::move(sumElements));
1549 _sumFunc.get()->addServerList(sumElements);
1550 _sumFunc.get()->addServerList(scaleElements);
1551
1552#ifdef USE_UBLAS
1553 std::cout.precision(std::numeric_limits<double>::digits);
1554#endif
1555 }
1556 //////////////////////////////////////////////////////////////////////////////
1557 /// create all the temporary objects required by the class
1558
1560 {
1561 std::string obsName = func->getObservable()->GetName();
1562 RooLagrangianMorphFunc::ParamSet values = getParams(func->_operators);
1563
1565
1566 cache->createComponents(func->_config.paramCards, func->_config.flagValues, func->GetName(), func->_diagrams,
1567 {func->_nonInterfering.begin(), func->_nonInterfering.end()}, func->_flags);
1568
1569 cache->buildMatrix(func->_config.paramCards, func->_config.flagValues, func->_flags);
1570 if (obsName.empty()) {
1571 std::cerr << "Matrix inversion succeeded, but no observable was "
1572 "supplied. quitting..."
1573 << std::endl;
1574 return cache;
1575 }
1576
1577 oocxcoutP((TObject *)0, ObjectHandling) << "observable: " << func->getObservable()->GetName() << std::endl;
1578 oocxcoutP((TObject *)0, ObjectHandling) << "binWidth: " << func->getBinWidth()->GetName() << std::endl;
1579
1580 setParams(func->_flags, 1);
1581 cache->buildMorphingFunction(func->GetName(), func->_config.paramCards, func->_sampleMap, func->_physics,
1582 func->_config.allowNegativeYields, func->getObservable(), func->getBinWidth());
1583 setParams(values, func->_operators, true);
1584 setParams(func->_flags, 1);
1585 return cache;
1586 }
1587
1588 //////////////////////////////////////////////////////////////////////////////
1589 /// create all the temporary objects required by the class
1590 /// function variant with precomputed inverse matrix
1591
1593 {
1594 RooLagrangianMorphFunc::ParamSet values = getParams(func->_operators);
1595
1597
1598 cache->createComponents(func->_config.paramCards, func->_config.flagValues, func->GetName(), func->_diagrams,
1599 {func->_nonInterfering.begin(), func->_nonInterfering.end()}, func->_flags);
1600
1601#ifndef USE_UBLAS
1602 cache->_inverse.ResizeTo(inverse.GetNrows(), inverse.GetNrows());
1603#endif
1604 cache->_inverse = inverse;
1605 cache->_condition = NaN;
1606
1607 setParams(func->_flags, 1);
1608 cache->buildMorphingFunction(func->GetName(), func->_config.paramCards, func->_sampleMap, func->_physics,
1609 func->_config.allowNegativeYields, func->getObservable(), func->getBinWidth());
1610 setParams(values, func->_operators, true);
1611 setParams(func->_flags, 1);
1612 return cache;
1613 }
1614};
1615
1616///////////////////////////////////////////////////////////////////////////////
1617// Class Implementation ///////////////////////////////////////////////////////
1618///////////////////////////////////////////////////////////////////////////////
1619
1620////////////////////////////////////////////////////////////////////////////////
1621/// write a matrix to a file
1622
1623void RooLagrangianMorphFunc::writeMatrixToFile(const TMatrixD &matrix, const char *fname)
1624{
1625 writeMatrixToFileT(matrix, fname);
1626}
1627
1628////////////////////////////////////////////////////////////////////////////////
1629/// write a matrix to a stream
1630
1631void RooLagrangianMorphFunc::writeMatrixToStream(const TMatrixD &matrix, std::ostream &stream)
1632{
1633 writeMatrixToStreamT(matrix, stream);
1634}
1635
1636////////////////////////////////////////////////////////////////////////////////
1637/// read a matrix from a text file
1638
1640{
1641 return readMatrixFromFileT<TMatrixD>(fname);
1642}
1643
1644////////////////////////////////////////////////////////////////////////////////
1645/// read a matrix from a stream
1646
1648{
1649 return readMatrixFromStreamT<TMatrixD>(stream);
1650}
1651
1652////////////////////////////////////////////////////////////////////////////////
1653/// setup observable, recycle existing observable if defined
1654
1656{
1657 // cxcoutP(ObjectHandling) << "setting up observable" << std::endl;
1658 RooRealVar *obs = nullptr;
1659 bool obsExists(false);
1660 if (_observables.at(0) != 0) {
1661 obs = (RooRealVar *)_observables.at(0);
1662 obsExists = true;
1663 }
1664
1665 if (mode && mode->InheritsFrom(RooHistFunc::Class())) {
1666 obs = (RooRealVar *)(dynamic_cast<RooHistFunc *>(inputExample)->getHistObsList().first());
1667 obsExists = true;
1668 _observables.add(*obs);
1669 } else if (mode && mode->InheritsFrom(RooParamHistFunc::Class())) {
1670 obs = (RooRealVar *)(dynamic_cast<RooParamHistFunc *>(inputExample)->paramList().first());
1671 obsExists = true;
1672 _observables.add(*obs);
1673 }
1674
1675 // Note: "found!" will be printed if s2 is a substring of s1, both s1 and s2
1676 // are of type std::string. s1.find(s2)
1677 // obtain the observable
1678 if (!obsExists) {
1679 if (mode && mode->InheritsFrom(TH1::Class())) {
1680 TH1 *hist = (TH1 *)(inputExample);
1681 obs = new RooRealVar(obsname, obsname, hist->GetXaxis()->GetXmin(), hist->GetXaxis()->GetXmax());
1682 obs->setBins(hist->GetNbinsX());
1683 } else {
1684 obs = new RooRealVar(obsname, obsname, 0, 1);
1685 obs->setBins(1);
1686 }
1687 _observables.add(*obs);
1688 } else {
1689 if (strcmp(obsname, obs->GetName()) != 0) {
1690 coutW(ObjectHandling) << " name of existing observable " << _observables.at(0)->GetName()
1691 << " does not match expected name " << obsname << std::endl;
1692 }
1693 }
1694
1695 TString sbw = Form("binWidth_%s", makeValidName(obs->GetName()).Data());
1696 RooRealVar *binWidth = new RooRealVar(sbw.Data(), sbw.Data(), 1.);
1697 double bw = obs->numBins() / (obs->getMax() - obs->getMin());
1698 binWidth->setVal(bw);
1699 binWidth->setConstant(true);
1700 _binWidths.add(*binWidth);
1701
1702 return obs;
1703}
1704
1705//#ifndef USE_MULTIPRECISION_LC
1706//#pragma GCC diagnostic push
1707//#pragma GCC diagnostic ignored "-Wunused-parameter"
1708//#endif
1709
1710////////////////////////////////////////////////////////////////////////////////
1711/// update sample weight (-?-)
1712
1714{
1715 //#ifdef USE_MULTIPRECISION_LC
1716 int sampleidx = 0;
1717 auto cache = this->getCache(_curNormSet);
1718 const size_t n(size(cache->_inverse));
1719 for (auto sampleit : _config.paramCards) {
1720 const std::string sample(sampleit.first);
1721 // build the formula with the correct normalization
1722 RooLinearCombination *sampleformula = dynamic_cast<RooLinearCombination *>(this->getSampleWeight(sample.c_str()));
1723 if (!sampleformula) {
1724 coutE(ObjectHandling) << Form("unable to access formula for sample '%s'!", sample.c_str()) << std::endl;
1725 return;
1726 }
1727 cxcoutP(ObjectHandling) << "updating formula for sample '" << sample << "'" << std::endl;
1728 for (size_t formulaidx = 0; formulaidx < n; ++formulaidx) {
1729 const RooFit::SuperFloat val(cache->_inverse(formulaidx, sampleidx));
1730#ifdef USE_UBLAS
1731 if (val != val) {
1732#else
1733 if (std::isnan(val)) {
1734#endif
1735 coutE(ObjectHandling) << "refusing to propagate NaN!" << std::endl;
1736 }
1737 cxcoutP(ObjectHandling) << " " << formulaidx << ":" << sampleformula->getCoefficient(formulaidx) << " -> "
1738 << val << std::endl;
1739 sampleformula->setCoefficient(formulaidx, val);
1740 assert(sampleformula->getCoefficient(formulaidx) == val);
1741 }
1742 sampleformula->setValueDirty();
1743 ++sampleidx;
1744 }
1745 //#else
1746 // ERROR("updating sample weights currently not possible without boost!");
1747 //#endif
1748}
1749//#ifndef USE_MULTIPRECISION_LC
1750//#pragma GCC diagnostic pop
1751//#endif
1752
1753////////////////////////////////////////////////////////////////////////////////
1754/// read the parameters from the input file
1755
1757{
1758 readValues<double>(_config.paramCards, f, _config.folderNames, "param_card", true);
1759 readValues<int>(_config.flagValues, f, _config.folderNames, "flags", false);
1760}
1761
1762////////////////////////////////////////////////////////////////////////////////
1763/// retrieve the physics inputs
1764
1766{
1767 std::string obsName = _config.observableName;
1768 cxcoutP(InputArguments) << "initializing physics inputs from file " << file->GetName() << " with object name(s) '"
1769 << obsName << "'" << std::endl;
1770 auto folderNames = _config.folderNames;
1771 TObject *obj = loadFromFileResidentFolder<TObject>(file, folderNames.front(), obsName, true);
1772 if (!obj) {
1773 std::cerr << "unable to locate object '" << obsName << "' in folder '" << folderNames.front() << "'!"
1774 << std::endl;
1775 return;
1776 }
1777 std::string classname = obj->ClassName();
1779
1780 RooRealVar *observable = this->setupObservable(obsName.c_str(), mode, obj);
1781 if (classname.find("TH1") != std::string::npos) {
1782 collectHistograms(this->GetName(), file, _sampleMap, _physics, *observable, obsName,
1784 } else if (classname.find("RooHistFunc") != std::string::npos ||
1785 classname.find("RooParamHistFunc") != std::string::npos ||
1786 classname.find("PiecewiseInterpolation") != std::string::npos) {
1787 collectRooAbsReal(this->GetName(), file, _sampleMap, _physics, obsName, _config.paramCards);
1788 } else if (classname.find("TParameter<double>") != std::string::npos) {
1789 collectCrosssections<double>(this->GetName(), file, _sampleMap, _physics, obsName,
1791 } else if (classname.find("TParameter<float>") != std::string::npos) {
1792 collectCrosssections<float>(this->GetName(), file, _sampleMap, _physics, obsName,
1794 } else if (classname.find("TPair") != std::string::npos) {
1795 collectCrosssectionsTPair(this->GetName(), file, _sampleMap, _physics, obsName, folderNames[0],
1797 } else {
1798 std::cerr << "cannot morph objects of class '" << mode->GetName() << "'!" << std::endl;
1799 }
1800}
1801
1802////////////////////////////////////////////////////////////////////////////////
1803/// convert the RooArgList folders into a simple vector of std::string
1804
1806{
1807 for (auto const &folder : folders) {
1808 RooStringVar *var = dynamic_cast<RooStringVar *>(folder);
1809 const std::string sample(var ? var->getVal() : folder->GetName());
1810 if (sample.empty())
1811 continue;
1812 _config.folderNames.push_back(sample);
1813 }
1814
1815 TDirectory *file = openFile(_config.fileName);
1816 TIter next(file->GetList());
1817 TObject *obj = nullptr;
1818 while ((obj = (TObject *)next())) {
1819 auto f = readOwningFolderFromFile(file, obj->GetName());
1820 if (!f)
1821 continue;
1822 std::string name(f->GetName());
1823 if (name.empty())
1824 continue;
1825 _config.folderNames.push_back(name);
1826 }
1827 closeFile(file);
1828}
1829
1830////////////////////////////////////////////////////////////////////////////////
1831/// print all the parameters and their values in the given sample to the console
1832
1833void RooLagrangianMorphFunc::printParameters(const char *samplename) const
1834{
1835 for (const auto &param : _config.paramCards.at(samplename)) {
1836 if (this->hasParameter(param.first.c_str())) {
1837 std::cout << param.first << " = " << param.second;
1838 if (this->isParameterConstant(param.first.c_str()))
1839 std::cout << " (const)";
1840 std::cout << std::endl;
1841 }
1842 }
1843}
1844
1845////////////////////////////////////////////////////////////////////////////////
1846/// print all the known samples to the console
1847
1849{
1850 // print all the known samples to the console
1851 for (auto folder : _config.folderNames) {
1852 std::cout << folder << std::endl;
1853 }
1854}
1855
1856////////////////////////////////////////////////////////////////////////////////
1857/// print the current physics values
1858
1860{
1861 for (const auto &sample : _sampleMap) {
1862 RooAbsArg *phys = _physics.at(sample.second);
1863 if (!phys)
1864 continue;
1865 phys->Print();
1866 }
1867}
1868
1869////////////////////////////////////////////////////////////////////////////////
1870/// constructor with proper arguments
1871
1872RooLagrangianMorphFunc::RooLagrangianMorphFunc(const char *name, const char *title, const Config &config)
1873 : RooAbsReal(name, title), _cacheMgr(this, 10, true, true), _physics("physics", "physics", this),
1874 _operators("operators", "set of operators", this), _observables("observables", "morphing observables", this),
1875 _binWidths("binWidths", "set of binWidth objects", this), _flags("flags", "flags", this), _config(config)
1876{
1877 this->init();
1878 this->setup(false);
1879
1881}
1882
1883////////////////////////////////////////////////////////////////////////////////
1884/// constructor with proper arguments
1885RooLagrangianMorphFunc::RooLagrangianMorphFunc(const char *name, const char *title, const char *filename,
1886 const char *observableName, const RooArgSet &couplings,
1887 const RooArgSet &folders)
1888 : RooAbsReal(name, title), _cacheMgr(this, 10, true, true), _physics("physics", "physics", this),
1889 _operators("operators", "set of operators", this), _observables("observables", "morphing observables", this),
1890 _binWidths("binWidths", "set of binWidth objects", this), _flags("flags", "flags", this)
1891{
1893 _config.observableName = observableName;
1894 _config.couplings.add(couplings);
1895 this->addFolders(folders);
1896 this->init();
1897 this->setup(false);
1898
1900}
1901
1902////////////////////////////////////////////////////////////////////////////////
1903/// setup this instance with the given set of operators and vertices
1904/// if own=true, the class will own the operators template `<class Base>`
1905
1907{
1908 _ownParameters = own;
1909 auto diagrams = _diagrams;
1910
1911 if (diagrams.size() > 0) {
1912 // TODO: The setup() function is a protected function that is only called
1913 // in the RooLagrangianMorphFunc constructor when _diagrams is still
1914 // empty. As this code branch is unreachable, can't it be removed?
1915 RooArgList operators;
1916 for (auto const &v : diagrams) {
1917 for (const RooArgList *t : v) {
1918 extractOperators(*t, operators);
1919 }
1920 }
1921
1922 if (own) {
1923 _operators.addOwned(operators);
1924 } else {
1925 _operators.add(operators);
1926 }
1927
1928 for (size_t j = 0; j < diagrams.size(); ++j) {
1929 std::vector<RooListProxy *> diagram;
1930 for (size_t i = 0; i < diagrams[j].size(); ++i) {
1931 std::stringstream name;
1932 name << "!vertex" << i;
1933 std::stringstream title;
1934 title << "set of couplings in the vertex " << i;
1935 diagram.push_back(new RooListProxy(name.str().c_str(), title.str().c_str(), this, true, false));
1936 if (own) {
1937 diagram[i]->addOwned(*diagrams[j][i]);
1938 } else {
1939 diagram[i]->add(*diagrams[j][i]);
1940 }
1941 }
1942 _diagrams.push_back(diagram);
1943 }
1944 }
1945
1946 else if (_config.couplings.size() > 0) {
1947 RooArgList operators;
1948 std::vector<RooListProxy *> vertices;
1949 extractOperators(_config.couplings, operators);
1950 vertices.push_back(new RooListProxy("!couplings", "set of couplings in the vertex", this, true, false));
1951 if (own) {
1952 _operators.addOwned(operators);
1953 vertices[0]->addOwned(_config.couplings);
1954 } else {
1955 _operators.add(operators);
1956 vertices[0]->add(_config.couplings);
1957 }
1958 _diagrams.push_back(vertices);
1959 }
1960
1961 else if (_config.prodCouplings.size() > 0 && _config.decCouplings.size() > 0) {
1962 std::vector<RooListProxy *> vertices;
1963 RooArgList operators;
1964 cxcoutP(InputArguments) << "prod/dec couplings provided" << std::endl;
1965 extractOperators(_config.prodCouplings, operators);
1966 extractOperators(_config.decCouplings, operators);
1967 vertices.push_back(
1968 new RooListProxy("!production", "set of couplings in the production vertex", this, true, false));
1969 vertices.push_back(new RooListProxy("!decay", "set of couplings in the decay vertex", this, true, false));
1970 if (own) {
1971 _operators.addOwned(operators);
1972 vertices[0]->addOwned(_config.prodCouplings);
1973 vertices[1]->addOwned(_config.decCouplings);
1974 } else {
1975 cxcoutP(InputArguments) << "adding non-own operators" << std::endl;
1976 _operators.add(operators);
1977 vertices[0]->add(_config.prodCouplings);
1978 vertices[1]->add(_config.decCouplings);
1979 }
1980 _diagrams.push_back(vertices);
1981 }
1982}
1983
1984////////////////////////////////////////////////////////////////////////////////
1985/// (-?-)
1986
1988{
1989 std::string filename = _config.fileName;
1990 TDirectory *file = openFile(filename.c_str());
1991 if (!file) {
1992 coutE(InputArguments) << "unable to open file '" << filename << "'!" << std::endl;
1993 return;
1994 }
1995 this->readParameters(file);
1996 checkNameConflict(_config.paramCards, _operators);
1997 this->collectInputs(file);
1998 closeFile(file);
1999 RooRealVar *nNP0 = new RooRealVar("nNP0", "nNP0", 1., 0, 1.);
2000 nNP0->setStringAttribute("NewPhysics", "0");
2001 nNP0->setConstant(true);
2002 _flags.add(*nNP0);
2003 RooRealVar *nNP1 = new RooRealVar("nNP1", "nNP1", 1., 0, 1.);
2004 nNP1->setStringAttribute("NewPhysics", "1");
2005 nNP1->setConstant(true);
2006 _flags.add(*nNP1);
2007 RooRealVar *nNP2 = new RooRealVar("nNP2", "nNP2", 1., 0, 1.);
2008 nNP2->setStringAttribute("NewPhysics", "2");
2009 nNP2->setConstant(true);
2010 _flags.add(*nNP2);
2011 RooRealVar *nNP3 = new RooRealVar("nNP3", "nNP3", 1., 0, 1.);
2012 nNP3->setStringAttribute("NewPhysics", "3");
2013 nNP3->setConstant(true);
2014 _flags.add(*nNP3);
2015 RooRealVar *nNP4 = new RooRealVar("nNP4", "nNP4", 1., 0, 1.);
2016 nNP4->setStringAttribute("NewPhysics", "4");
2017 nNP4->setConstant(true);
2018 _flags.add(*nNP4);
2019 // we can't use `addOwned` before, because the RooListProxy doesn't overload
2020 // `addOwned` correctly (it might in the future, then this can be changed).
2022}
2023
2024////////////////////////////////////////////////////////////////////////////////
2025/// copy constructor
2026
2028 : RooAbsReal(other, name), _cacheMgr(other._cacheMgr, this), _scale(other._scale), _sampleMap(other._sampleMap),
2029 _physics(other._physics.GetName(), this, other._physics),
2030 _operators(other._operators.GetName(), this, other._operators),
2031 _observables(other._observables.GetName(), this, other._observables),
2032 _binWidths(other._binWidths.GetName(), this, other._binWidths), _flags{other._flags.GetName(), this, other._flags},
2033 _config(other._config)
2034{
2035 for (size_t j = 0; j < other._diagrams.size(); ++j) {
2036 std::vector<RooListProxy *> diagram;
2037 for (size_t i = 0; i < other._diagrams[j].size(); ++i) {
2038 RooListProxy *list = new RooListProxy(other._diagrams[j][i]->GetName(), this, *(other._diagrams[j][i]));
2039 diagram.push_back(list);
2040 }
2041 _diagrams.push_back(diagram);
2042 }
2044}
2045
2046////////////////////////////////////////////////////////////////////////////////
2047/// set energy scale of the EFT expansion
2048
2050{
2051 _scale = val;
2052}
2053
2054////////////////////////////////////////////////////////////////////////////////
2055/// get energy scale of the EFT expansion
2056
2058{
2059 return _scale;
2060}
2061
2062////////////////////////////////////////////////////////////////////////////////
2063// default constructor
2064
2066 : _cacheMgr(this, 10, true, true), _operators("operators", "set of operators", this, true, false),
2067 _observables("observable", "morphing observable", this, true, false),
2068 _binWidths("binWidths", "set of bin width objects", this, true, false)
2069{
2070 static int counter(0);
2071 counter++;
2073}
2074
2075////////////////////////////////////////////////////////////////////////////////
2076/// default destructor
2077
2079{
2080 for (auto const &diagram : _diagrams) {
2081 for (RooListProxy *vertex : diagram) {
2082 delete vertex;
2083 }
2084 }
2086}
2087
2088////////////////////////////////////////////////////////////////////////////////
2089/// cloning method
2090
2091TObject *RooLagrangianMorphFunc::clone(const char *newname) const
2092{
2093 return new RooLagrangianMorphFunc(*this, newname);
2094}
2095
2096////////////////////////////////////////////////////////////////////////////////
2097/// calculate the number of samples needed to morph a bivertex, 2-2 physics
2098/// process
2099
2100int RooLagrangianMorphFunc::countSamples(int nprod, int ndec, int nboth)
2101{
2102 FeynmanDiagram diagram;
2103 std::vector<bool> prod;
2104 std::vector<bool> dec;
2105 for (int i = 0; i < nboth; ++i) {
2106 prod.push_back(true);
2107 dec.push_back(true);
2108 }
2109 for (int i = 0; i < nprod; ++i) {
2110 prod.push_back(true);
2111 dec.push_back(false);
2112 }
2113 for (int i = 0; i < ndec; ++i) {
2114 prod.push_back(false);
2115 dec.push_back(true);
2116 }
2117 diagram.push_back(prod);
2118 diagram.push_back(dec);
2119 MorphFuncPattern morphfuncpattern;
2120 ::collectPolynomials(morphfuncpattern, diagram);
2121 return morphfuncpattern.size();
2122}
2123
2124////////////////////////////////////////////////////////////////////////////////
2125/// calculate the number of samples needed to morph a certain physics process
2126/// usage:
2127/// countSamples ( { RooLagrangianMorphFunc::makeHCggfCouplings(),
2128/// RooLagrangianMorphFunc::makeHCHZZCouplings() } )
2129
2130int RooLagrangianMorphFunc::countSamples(std::vector<RooArgList *> &vertices)
2131{
2132 RooArgList operators, couplings;
2133 for (auto vertex : vertices) {
2134 extractOperators(*vertex, operators);
2135 extractCouplings(*vertex, couplings);
2136 }
2137 FeynmanDiagram diagram;
2138 ::fillFeynmanDiagram(diagram, vertices, couplings);
2139 MorphFuncPattern morphfuncpattern;
2140 ::collectPolynomials(morphfuncpattern, diagram);
2141 return morphfuncpattern.size();
2142}
2143
2144////////////////////////////////////////////////////////////////////////////////
2145/// create TPair containers of the type expected by the RooLagrangianMorph
2146
2148{
2149 TPair *v = new TPair(new TParameter<double>("xsection", xs), new TParameter<double>("uncertainty", unc));
2150 return v;
2151}
2152
2153////////////////////////////////////////////////////////////////////////////////
2154/// create only the weight formulas. static function for external usage.
2155
2156std::map<std::string, std::string>
2158 const std::vector<std::vector<std::string>> &vertices_str)
2159{
2160 std::stack<RooArgList> ownedVertices;
2161 std::vector<RooArgList *> vertices;
2162 RooArgList couplings;
2163 for (const auto &vtx : vertices_str) {
2164 ownedVertices.emplace();
2165 auto &vertex = ownedVertices.top();
2166 for (const auto &c : vtx) {
2167 RooRealVar *coupling = (RooRealVar *)(couplings.find(c.c_str()));
2168 if (!coupling) {
2169 coupling = new RooRealVar(c.c_str(), c.c_str(), 1., 0., 10.);
2170 couplings.add(*coupling);
2171 }
2172 vertex.add(*coupling);
2173 }
2174 vertices.push_back(&vertex);
2175 }
2176 auto retval = RooLagrangianMorphFunc::createWeightStrings(inputs, vertices, couplings);
2177 return retval;
2178}
2179
2180////////////////////////////////////////////////////////////////////////////////
2181/// create only the weight formulas. static function for external usage.
2182
2183std::map<std::string, std::string>
2185 const std::vector<RooArgList *> &vertices, RooArgList &couplings)
2186{
2187 return createWeightStrings(inputs, vertices, couplings, {}, {}, {});
2188}
2189
2190////////////////////////////////////////////////////////////////////////////////
2191/// create only the weight formulas. static function for external usage.
2192
2193std::map<std::string, std::string>
2195 const std::vector<RooArgList *> &vertices, RooArgList &couplings,
2196 const RooLagrangianMorphFunc::FlagMap &flagValues, const RooArgList &flags,
2197 const std::vector<RooArgList *> &nonInterfering)
2198{
2199 FormulaList formulas = ::createFormulas("", inputs, flagValues, {vertices}, couplings, flags, nonInterfering);
2200 RooArgSet operators;
2201 extractOperators(couplings, operators);
2202 Matrix matrix(::buildMatrixT<Matrix>(inputs, formulas, operators, flagValues, flags));
2203 if (size(matrix) < 1) {
2204 std::cerr << "input matrix is empty, please provide suitable input samples!" << std::endl;
2205 }
2206 Matrix inverse(::diagMatrix(size(matrix)));
2207 double condition __attribute__((unused)) = (double)(invertMatrix(matrix, inverse));
2208 auto retval = buildSampleWeightStrings(inputs, formulas, inverse);
2209 return retval;
2210}
2211
2212////////////////////////////////////////////////////////////////////////////////
2213/// create only the weight formulas. static function for external usage.
2214
2216 const std::vector<RooArgList *> &vertices, RooArgList &couplings,
2217 const RooLagrangianMorphFunc::FlagMap &flagValues,
2218 const RooArgList &flags,
2219 const std::vector<RooArgList *> &nonInterfering)
2220{
2221 FormulaList formulas = ::createFormulas("", inputs, flagValues, {vertices}, couplings, flags, nonInterfering);
2222 RooArgSet operators;
2223 extractOperators(couplings, operators);
2224 Matrix matrix(::buildMatrixT<Matrix>(inputs, formulas, operators, flagValues, flags));
2225 if (size(matrix) < 1) {
2226 std::cerr << "input matrix is empty, please provide suitable input samples!" << std::endl;
2227 }
2228 Matrix inverse(::diagMatrix(size(matrix)));
2229 double condition __attribute__((unused)) = (double)(invertMatrix(matrix, inverse));
2230 RooArgSet retval;
2231 ::buildSampleWeights(retval, (const char *)nullptr /* name */, inputs, formulas, inverse);
2232 return retval;
2233}
2234
2235////////////////////////////////////////////////////////////////////////////////
2236/// create only the weight formulas. static function for external usage.
2237
2239 const std::vector<RooArgList *> &vertices, RooArgList &couplings)
2240{
2241 std::vector<RooArgList *> nonInterfering;
2242 RooArgList flags;
2243 FlagMap flagValues;
2244 return RooLagrangianMorphFunc::createWeights(inputs, vertices, couplings, flagValues, flags, nonInterfering);
2245}
2246
2247////////////////////////////////////////////////////////////////////////////////
2248/// return the RooProduct that is the element of the RooRealSumPdfi
2249/// corresponding to the given sample name
2250
2252{
2253 auto mf = this->getFunc();
2254 if (!mf) {
2255 coutE(Eval) << "unable to retrieve morphing function" << std::endl;
2256 return nullptr;
2257 }
2258 RooArgSet *args = mf->getComponents();
2259 TString prodname(name);
2260 prodname.Append("_");
2261 prodname.Append(this->GetName());
2262
2263 for (auto itr : *args) {
2264 RooProduct *prod = dynamic_cast<RooProduct *>(itr);
2265 if (!prod)
2266 continue;
2267 TString sname(prod->GetName());
2268 if (sname.CompareTo(prodname) == 0) {
2269 return prod;
2270 }
2271 }
2272 return nullptr;
2273}
2274////////////////////////////////////////////////////////////////////////////////
2275/// return the vector of sample names, used to build the morph func
2276
2277std::vector<std::string> RooLagrangianMorphFunc::getSamples() const
2278{
2279 return _config.folderNames;
2280}
2281
2282////////////////////////////////////////////////////////////////////////////////
2283/// retrieve the weight (prefactor) of a sample with the given name
2284
2286{
2287 auto cache = this->getCache(_curNormSet);
2288 auto wname = std::string("w_") + name + "_" + this->GetName();
2289 return dynamic_cast<RooAbsReal *>(cache->_weights.find(wname.c_str()));
2290}
2291
2292////////////////////////////////////////////////////////////////////////////////
2293/// print the current sample weights
2294
2296{
2297 this->printSampleWeights();
2298}
2299
2300////////////////////////////////////////////////////////////////////////////////
2301/// print the current sample weights
2302
2304{
2305 auto *cache = this->getCache(_curNormSet);
2306 for (const auto &sample : _sampleMap) {
2307 auto weightName = std::string("w_") + sample.first + "_" + this->GetName();
2308 auto weight = static_cast<RooAbsReal *>(cache->_weights.find(weightName.c_str()));
2309 if (!weight)
2310 continue;
2311 }
2312}
2313
2314////////////////////////////////////////////////////////////////////////////////
2315/// randomize the parameters a bit
2316/// useful to test and debug fitting
2317
2319{
2320 RooRealVar *obj;
2321 TRandom3 r;
2322
2323 for (auto itr : _operators) {
2324 obj = dynamic_cast<RooRealVar *>(itr);
2325 double val = obj->getVal();
2326 if (obj->isConstant())
2327 continue;
2328 double variation = r.Gaus(1, z);
2329 obj->setVal(val * variation);
2330 }
2331}
2332
2333////////////////////////////////////////////////////////////////////////////////
2334/// retrive the new physics objects and update the weights in the morphing
2335/// function
2336
2338{
2339 auto cache = this->getCache(_curNormSet);
2340
2341 std::string filename = _config.fileName;
2342 TDirectory *file = openFile(filename.c_str());
2343 if (!file) {
2344 coutE(InputArguments) << "unable to open file '" << filename << "'!" << std::endl;
2345 return false;
2346 }
2347
2348 this->readParameters(file);
2349
2350 checkNameConflict(_config.paramCards, _operators);
2351 this->collectInputs(file);
2352
2353 cache->buildMatrix(_config.paramCards, _config.flagValues, _flags);
2354 this->updateSampleWeights();
2355
2356 closeFile(file);
2357 return true;
2358}
2359
2360////////////////////////////////////////////////////////////////////////////////
2361/// setup the morphing function with a predefined inverse matrix
2362/// call this function *before* any other after creating the object
2363
2365{
2366 auto cache = static_cast<RooLagrangianMorphFunc::CacheElem *>(_cacheMgr.getObj(0, (RooArgSet *)0));
2367 Matrix m = makeSuperMatrix(inverse);
2368 if (cache) {
2369 std::string filename = _config.fileName;
2370 cache->_inverse = m;
2371 TDirectory *file = openFile(filename.c_str());
2372 if (!file) {
2373 coutE(InputArguments) << "unable to open file '" << filename << "'!" << std::endl;
2374 return false;
2375 }
2376
2377 this->readParameters(file);
2378 checkNameConflict(_config.paramCards, _operators);
2379 this->collectInputs(file);
2380
2381 // then, update the weights in the morphing function
2382 this->updateSampleWeights();
2383
2384 closeFile(file);
2385 } else {
2387 if (!cache)
2388 coutE(Caching) << "unable to create cache!" << std::endl;
2389 _cacheMgr.setObj(0, 0, cache, 0);
2390 }
2391 return true;
2392}
2393
2394////////////////////////////////////////////////////////////////////////////////
2395// setup the morphing function with a predefined inverse matrix
2396// call this function *before* any other after creating the object
2397
2399{
2400 auto cache = static_cast<RooLagrangianMorphFunc::CacheElem *>(_cacheMgr.getObj(0, (RooArgSet *)0));
2401 if (cache) {
2402 return false;
2403 }
2404 cache = RooLagrangianMorphFunc::CacheElem::createCache(this, readMatrixFromFileT<Matrix>(filename));
2405 if (!cache)
2406 coutE(Caching) << "unable to create cache!" << std::endl;
2407 _cacheMgr.setObj(0, 0, cache, 0);
2408 return true;
2409}
2410
2411////////////////////////////////////////////////////////////////////////////////
2412/// write the inverse matrix to a file
2413
2415{
2416 auto cache = this->getCache(_curNormSet);
2417 if (!cache)
2418 return false;
2419 writeMatrixToFileT(cache->_inverse, filename);
2420 return true;
2421}
2422
2423////////////////////////////////////////////////////////////////////////////////
2424/// retrieve the cache object
2425
2427{
2428 auto cache = static_cast<RooLagrangianMorphFunc::CacheElem *>(_cacheMgr.getObj(0, (RooArgSet *)0));
2429 if (!cache) {
2430 cxcoutP(Caching) << "creating cache from getCache function for " << this << std::endl;
2431 cxcoutP(Caching) << "current storage has size " << _sampleMap.size() << std::endl;
2433 if (cache)
2434 _cacheMgr.setObj(0, 0, cache, 0);
2435 else
2436 coutE(Caching) << "unable to create cache!" << std::endl;
2437 }
2438 return cache;
2439}
2440
2441////////////////////////////////////////////////////////////////////////////////
2442/// return true if a cache object is present, false otherwise
2443
2445{
2446 return (bool)(_cacheMgr.getObj(0, (RooArgSet *)0));
2447}
2448
2449////////////////////////////////////////////////////////////////////////////////
2450/// set one parameter to a specific value
2451
2453{
2454 RooRealVar *param = this->getParameter(name);
2455 if (!param) {
2456 return;
2457 }
2458 if (value > param->getMax())
2459 param->setMax(value);
2460 if (value < param->getMin())
2461 param->setMin(value);
2462 param->setVal(value);
2463}
2464
2465////////////////////////////////////////////////////////////////////////////////
2466/// set one flag to a specific value
2467
2469{
2470 RooRealVar *param = this->getFlag(name);
2471 if (!param) {
2472 return;
2473 }
2474 param->setVal(value);
2475}
2476
2477////////////////////////////////////////////////////////////////////////////////
2478/// set one parameter to a specific value and range
2479
2480void RooLagrangianMorphFunc::setParameter(const char *name, double value, double min, double max)
2481{
2482 RooRealVar *param = this->getParameter(name);
2483 if (!param) {
2484 return;
2485 }
2486 param->setMin(min);
2487 param->setMax(max);
2488 param->setVal(value);
2489}
2490
2491////////////////////////////////////////////////////////////////////////////////
2492/// set one parameter to a specific value and range
2493void RooLagrangianMorphFunc::setParameter(const char *name, double value, double min, double max, double error)
2494{
2495 RooRealVar *param = this->getParameter(name);
2496 if (!param) {
2497 return;
2498 }
2499 param->setMin(min);
2500 param->setMax(max);
2501 param->setVal(value);
2502 param->setError(error);
2503}
2504
2505////////////////////////////////////////////////////////////////////////////////
2506/// return true if the parameter with the given name is set constant, false
2507/// otherwise
2508
2510{
2511 RooRealVar *param = this->getParameter(name);
2512 if (param) {
2513 return param->isConstant();
2514 }
2515 return true;
2516}
2517
2518////////////////////////////////////////////////////////////////////////////////
2519/// retrieve the RooRealVar object incorporating the parameter with the given
2520/// name
2522{
2523
2524 return dynamic_cast<RooRealVar *>(_operators.find(name));
2525}
2526
2527////////////////////////////////////////////////////////////////////////////////
2528/// retrieve the RooRealVar object incorporating the flag with the given name
2529
2531{
2532 return dynamic_cast<RooRealVar *>(_flags.find(name));
2533}
2534
2535////////////////////////////////////////////////////////////////////////////////
2536/// check if a parameter of the given name is contained in the list of known
2537/// parameters
2538
2540{
2541 return this->getParameter(name);
2542}
2543
2544////////////////////////////////////////////////////////////////////////////////
2545/// call setConstant with the boolean argument provided on the parameter with
2546/// the given name
2547
2548void RooLagrangianMorphFunc::setParameterConstant(const char *name, bool constant) const
2549{
2550 RooRealVar *param = this->getParameter(name);
2551 if (param) {
2552 return param->setConstant(constant);
2553 }
2554}
2555
2556////////////////////////////////////////////////////////////////////////////////
2557/// set one parameter to a specific value
2558
2560{
2561 RooRealVar *param = this->getParameter(name);
2562 if (param) {
2563 return param->getVal();
2564 }
2565 return 0;
2566}
2567
2568////////////////////////////////////////////////////////////////////////////////
2569/// set the morphing parameters to those supplied in the given param hist
2570
2572{
2573 setParams(paramhist, _operators, false);
2574}
2575
2576////////////////////////////////////////////////////////////////////////////////
2577/// set the morphing parameters to those supplied in the sample with the given
2578/// name
2579
2580void RooLagrangianMorphFunc::setParameters(const char *foldername)
2581{
2582 std::string filename = _config.fileName;
2583 TDirectory *file = openFile(filename.c_str());
2584 TH1 *paramhist = loadFromFileResidentFolder<TH1>(file, foldername, "param_card");
2585 setParams(paramhist, _operators, false);
2586 closeFile(file);
2587}
2588
2589/////////////////////////////////////////////////////////////////////////////////
2590/// retrieve the morphing parameters associated to the sample with the given
2591/// name
2592
2594{
2595 const std::string name(foldername);
2596 return _config.paramCards.at(name);
2597}
2598
2599////////////////////////////////////////////////////////////////////////////////
2600/// set the morphing parameters to those supplied in the list with the given
2601/// name
2602
2604{
2605 for (auto itr : *list) {
2606 RooRealVar *param = dynamic_cast<RooRealVar *>(itr);
2607 if (!param)
2608 continue;
2609 this->setParameter(param->GetName(), param->getVal());
2610 }
2611}
2612
2613////////////////////////////////////////////////////////////////////////////////
2614/// retrieve the histogram observable
2615
2617{
2618 if (_observables.getSize() < 1) {
2619 coutE(InputArguments) << "observable not available!" << std::endl;
2620 return nullptr;
2621 }
2622 return static_cast<RooRealVar *>(_observables.at(0));
2623}
2624
2625////////////////////////////////////////////////////////////////////////////////
2626/// retrieve the histogram observable
2627
2629{
2630 if (_binWidths.getSize() < 1) {
2631 coutE(InputArguments) << "bin width not available!" << std::endl;
2632 return nullptr;
2633 }
2634 return static_cast<RooRealVar *>(_binWidths.at(0));
2635}
2636
2637////////////////////////////////////////////////////////////////////////////////
2638/// retrieve a histogram output of the current morphing settings
2639
2641 /*RooFitResult *r*/)
2642{
2643 return this->createTH1(name, false /*r*/);
2644}
2645
2646////////////////////////////////////////////////////////////////////////////////
2647/// retrieve a histogram output of the current morphing settings
2648
2649TH1 *RooLagrangianMorphFunc::createTH1(const std::string &name, bool correlateErrors
2650 /* RooFitResult *r*/)
2651{
2652 auto mf = std::make_unique<RooRealSumFunc>(*(this->getFunc()));
2653 RooRealVar *observable = this->getObservable();
2654
2655 const int nbins = observable->getBins();
2656
2657 TH1 *hist = new TH1F{name.c_str(), name.c_str(), nbins, observable->getBinning().array()};
2658
2659 RooArgSet *args = mf->getComponents();
2660 for (int i = 0; i < nbins; ++i) {
2661 observable->setBin(i);
2662 double val = 0;
2663 double unc2 = 0;
2664 double unc = 0;
2665 for (auto itr : *args) {
2666 RooProduct *prod = dynamic_cast<RooProduct *>(itr);
2667 if (!prod)
2668 continue;
2669 RooAbsArg *phys = prod->components().find(Form("phys_%s", prod->GetName()));
2670 RooHistFunc *hf = dynamic_cast<RooHistFunc *>(phys);
2671 if (!hf) {
2672 continue;
2673 }
2674 const RooDataHist &dhist = hf->dataHist();
2675 dhist.get(i);
2676 RooAbsReal *formula = dynamic_cast<RooAbsReal *>(prod->components().find(Form("w_%s", prod->GetName())));
2677 double weight = formula->getVal();
2678 unc2 += dhist.weightSquared() * weight * weight;
2679 unc += sqrt(dhist.weightSquared()) * weight;
2680 val += dhist.weight() * weight;
2681 }
2682 hist->SetBinContent(i + 1, val);
2683 hist->SetBinError(i + 1, correlateErrors ? unc : sqrt(unc2));
2684 }
2685 return hist;
2686}
2687
2688////////////////////////////////////////////////////////////////////////////////
2689/// count the number of formulas that correspond to the current parameter set
2690
2692{
2693 int nFormulas = 0;
2694 auto mf = std::make_unique<RooRealSumFunc>(*(this->getFunc()));
2695 if (!mf)
2696 coutE(InputArguments) << "unable to retrieve morphing function" << std::endl;
2697 RooArgSet *args = mf->getComponents();
2698 for (auto itr : *args) {
2699 RooProduct *prod = dynamic_cast<RooProduct *>(itr);
2700 if (prod->getVal() != 0) {
2701 nFormulas++;
2702 }
2703 }
2704 return nFormulas;
2705}
2706
2707////////////////////////////////////////////////////////////////////////////////
2708/// check if there is any morphing power provided for the given parameter
2709/// morphing power is provided as soon as any two samples provide different,
2710/// non-zero values for this parameter
2711
2712bool RooLagrangianMorphFunc::isParameterUsed(const char *paramname) const
2713{
2714 std::string pname(paramname);
2715 double val = 0;
2716 bool isUsed = false;
2717 for (const auto &sample : _config.paramCards) {
2718 double thisval = sample.second.at(pname);
2719 if (thisval != val) {
2720 if (val != 0)
2721 isUsed = true;
2722 val = thisval;
2723 }
2724 }
2725 return isUsed;
2726}
2727
2728////////////////////////////////////////////////////////////////////////////////
2729/// check if there is any morphing power provided for the given coupling
2730/// morphing power is provided as soon as any two samples provide
2731/// different, non-zero values for this coupling
2732
2734{
2735 std::string cname(couplname);
2736 const RooArgList *args = this->getCouplingSet();
2737 RooAbsReal *coupling = dynamic_cast<RooAbsReal *>(args->find(couplname));
2738 if (!coupling)
2739 return false;
2741 double val = 0;
2742 bool isUsed = false;
2743 for (const auto &sample : _config.paramCards) {
2744 this->setParameters(sample.second);
2745 double thisval = coupling->getVal();
2746 if (thisval != val) {
2747 if (val != 0)
2748 isUsed = true;
2749 val = thisval;
2750 }
2751 }
2752 this->setParameters(params);
2753 return isUsed;
2754}
2755
2756////////////////////////////////////////////////////////////////////////////////
2757/// return the number of parameters in this morphing function
2758
2760{
2761 return this->getParameterSet()->getSize();
2762}
2763
2764////////////////////////////////////////////////////////////////////////////////
2765/// return the number of samples in this morphing function
2766
2768{
2769 // return the number of samples in this morphing function
2770 auto cache = getCache(_curNormSet);
2771 return cache->_formulas.size();
2772}
2773
2774////////////////////////////////////////////////////////////////////////////////
2775/// print the contributing samples and their respective weights
2776
2778{
2779 auto mf = std::make_unique<RooRealSumFunc>(*(this->getFunc()));
2780 if (!mf) {
2781 std::cerr << "Error: unable to retrieve morphing function" << std::endl;
2782 return;
2783 }
2784 RooArgSet *args = mf->getComponents();
2785 for (auto itr : *args) {
2786 RooAbsReal *formula = dynamic_cast<RooAbsReal *>(itr);
2787 if (formula) {
2788 TString name(formula->GetName());
2789 name.Remove(0, 2);
2790 name.Prepend("phys_");
2791 if (!args->find(name.Data())) {
2792 continue;
2793 }
2794 double val = formula->getVal();
2795 if (val != 0) {
2796 std::cout << formula->GetName() << ": " << val << " = " << formula->GetTitle() << std::endl;
2797 }
2798 }
2799 }
2800}
2801
2802////////////////////////////////////////////////////////////////////////////////
2803/// get the set of parameters
2804
2806{
2807 return &(_operators);
2808}
2809
2810////////////////////////////////////////////////////////////////////////////////
2811/// get the set of couplings
2812
2814{
2815 auto cache = getCache(_curNormSet);
2816 return &(cache->_couplings);
2817}
2818
2819////////////////////////////////////////////////////////////////////////////////
2820/// retrieve a set of couplings (-?-)
2821
2823{
2825 for (auto obj : *(this->getCouplingSet())) {
2826 RooAbsReal *var = dynamic_cast<RooAbsReal *>(obj);
2827 if (!var)
2828 continue;
2829 const std::string name(var->GetName());
2830 double val = var->getVal();
2831 couplings[name] = val;
2832 }
2833 return couplings;
2834}
2835
2836////////////////////////////////////////////////////////////////////////////////
2837/// retrieve the parameter set
2838
2840{
2841 return getParams(_operators);
2842}
2843
2844////////////////////////////////////////////////////////////////////////////////
2845/// retrieve a set of couplings (-?-)
2846
2848{
2849 setParams(params, _operators, false);
2850}
2851
2852////////////////////////////////////////////////////////////////////////////////
2853/// (currently similar to cloning the Pdf
2854
2855std::unique_ptr<RooWrapperPdf> RooLagrangianMorphFunc::createPdf() const
2856{
2857 auto cache = getCache(_curNormSet);
2858 auto func = std::make_unique<RooRealSumFunc>(*(cache->_sumFunc));
2859
2860 // create a wrapper on the roorealsumfunc
2861 return std::make_unique<RooWrapperPdf>(Form("pdf_%s", func->GetName()), Form("pdf of %s", func->GetTitle()), *func);
2862}
2863
2864////////////////////////////////////////////////////////////////////////////////
2865/// get the func
2866
2868{
2869 auto cache = getCache(_curNormSet);
2870 return cache->_sumFunc.get();
2871}
2872
2873////////////////////////////////////////////////////////////////////////////////
2874/// return extended mored capabilities
2875
2877{
2878 return this->createPdf()->extendMode();
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 this->createPdf()->expectedEvents(nset);
2888}
2889
2890////////////////////////////////////////////////////////////////////////////////
2891/// return the number of expected events for the current parameter set
2892
2894{
2895 RooArgSet set;
2896 set.add(*this->getObservable());
2897 return this->createPdf()->expectedEvents(set);
2898}
2899
2900////////////////////////////////////////////////////////////////////////////////
2901/// return expected number of events for extended likelihood calculation,
2902/// this is the sum of all coefficients
2903
2905{
2906 return createPdf()->expectedEvents(&nset);
2907}
2908
2909////////////////////////////////////////////////////////////////////////////////
2910/// return the expected uncertainty for the current parameter set
2911
2913{
2914 RooRealVar *observable = this->getObservable();
2915 auto cache = this->getCache(_curNormSet);
2916 double unc2 = 0;
2917 for (const auto &sample : _sampleMap) {
2918 RooAbsArg *phys = _physics.at(sample.second);
2919 auto weightName = std::string("w_") + sample.first + "_" + this->GetName();
2920 auto weight = static_cast<RooAbsReal *>(cache->_weights.find(weightName.c_str()));
2921 if (!weight) {
2922 coutE(InputArguments) << "unable to find object " + weightName << std::endl;
2923 return 0.0;
2924 }
2925 double newunc2 = 0;
2926 RooHistFunc *hf = dynamic_cast<RooHistFunc *>(phys);
2927 RooRealVar *rv = dynamic_cast<RooRealVar *>(phys);
2928 if (hf) {
2929 const RooDataHist &hist = hf->dataHist();
2930 for (Int_t j = 0; j < observable->getBins(); ++j) {
2931 hist.get(j);
2932 newunc2 += hist.weightSquared();
2933 }
2934 } else if (rv) {
2935 newunc2 = pow(rv->getError(), 2);
2936 }
2937 double w = weight->getVal();
2938 unc2 += newunc2 * w * w;
2939 // std::cout << phys->GetName() << " : " << weight->GetName() << "
2940 // thisweight: " << w << " thisxsec2: " << newunc2 << " weight " << weight
2941 // << std::endl;
2942 }
2943 return sqrt(unc2);
2944}
2945
2946////////////////////////////////////////////////////////////////////////////////
2947/// print the parameters and their current values
2948
2950{
2951 // print the parameters and their current values
2952 for (auto obj : _operators) {
2953 RooRealVar *param = static_cast<RooRealVar *>(obj);
2954 if (!param)
2955 continue;
2956 param->Print();
2957 }
2958}
2959
2960////////////////////////////////////////////////////////////////////////////////
2961/// print the flags and their current values
2962
2964{
2965 for (auto flag : _flags) {
2966 RooRealVar *param = static_cast<RooRealVar *>(flag);
2967 if (!param)
2968 continue;
2969 param->Print();
2970 }
2971}
2972
2973////////////////////////////////////////////////////////////////////////////////
2974/// print a set of couplings
2975
2977{
2979 for (auto c : couplings) {
2980 std::cout << c.first << ": " << c.second << std::endl;
2981 }
2982}
2983
2984////////////////////////////////////////////////////////////////////////////////
2985/// retrieve the list of bin boundaries
2986
2987std::list<double> *RooLagrangianMorphFunc::binBoundaries(RooAbsRealLValue &obs, double xlo, double xhi) const
2988{
2989 return this->getFunc()->binBoundaries(obs, xlo, xhi);
2990}
2991
2992////////////////////////////////////////////////////////////////////////////////
2993/// retrieve the sample Hint
2994
2995std::list<double> *RooLagrangianMorphFunc::plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const
2996{
2997 return this->getFunc()->plotSamplingHint(obs, xlo, xhi);
2998}
2999
3000////////////////////////////////////////////////////////////////////////////////
3001/// call getVal on the internal function
3002
3004{
3005 // cout << "XX RooLagrangianMorphFunc::getValV(" << this << ") set = " << set
3006 // << std::endl ;
3007 _curNormSet = set;
3008 return RooAbsReal::getValV(set);
3009}
3010
3011////////////////////////////////////////////////////////////////////////////////
3012/// call getVal on the internal function
3013
3015{
3016 // call getVal on the internal function
3017 RooRealSumFunc *pdf = this->getFunc();
3018 if (pdf)
3019 return _scale * pdf->getVal(_curNormSet);
3020 else
3021 std::cerr << "unable to acquire in-built function!" << std::endl;
3022 return 0.;
3023}
3024
3025////////////////////////////////////////////////////////////////////////////////
3026/// check if this PDF is a binned distribution in the given observable
3027
3029{
3030 return this->getFunc()->isBinnedDistribution(obs);
3031}
3032
3033////////////////////////////////////////////////////////////////////////////////
3034/// check if observable exists in the RooArgSet (-?-)
3035
3037{
3038 return this->getFunc()->checkObservables(nset);
3039}
3040
3041////////////////////////////////////////////////////////////////////////////////
3042/// Force analytical integration for the given observable
3043
3045{
3046 return this->getFunc()->forceAnalyticalInt(arg);
3047}
3048
3049////////////////////////////////////////////////////////////////////////////////
3050/// Retrieve the mat
3051
3053 const char *rangeName) const
3054{
3055 return this->getFunc()->getAnalyticalIntegralWN(allVars, numVars, normSet, rangeName);
3056}
3057
3058////////////////////////////////////////////////////////////////////////////////
3059/// Retrieve the matrix of coefficients
3060
3061double RooLagrangianMorphFunc::analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName) const
3062{
3063 return this->getFunc()->analyticalIntegralWN(code, normSet, rangeName);
3064}
3065
3066////////////////////////////////////////////////////////////////////////////////
3067/// Retrieve the matrix of coefficients
3068
3069void RooLagrangianMorphFunc::printMetaArgs(std::ostream &os) const
3070{
3071 return this->getFunc()->printMetaArgs(os);
3072}
3073
3074////////////////////////////////////////////////////////////////////////////////
3075/// Retrieve the matrix of coefficients
3076
3078{
3079 return this->getFunc()->canNodeBeCached();
3080}
3081
3082////////////////////////////////////////////////////////////////////////////////
3083/// Retrieve the matrix of coefficients
3084
3086{
3087 this->getFunc()->setCacheAndTrackHints(arg);
3088}
3089
3090////////////////////////////////////////////////////////////////////////////////
3091/// Retrieve the matrix of coefficients
3092
3094{
3095 auto cache = getCache(_curNormSet);
3096 if (!cache)
3097 coutE(Caching) << "unable to retrieve cache!" << std::endl;
3098 return makeRootMatrix(cache->_matrix);
3099}
3100
3101////////////////////////////////////////////////////////////////////////////////
3102/// Retrieve the matrix of coefficients after inversion
3103
3105{
3106 auto cache = getCache(_curNormSet);
3107 if (!cache)
3108 coutE(Caching) << "unable to retrieve cache!" << std::endl;
3109 return makeRootMatrix(cache->_inverse);
3110}
3111
3112////////////////////////////////////////////////////////////////////////////////
3113/// Retrieve the condition of the coefficient matrix. If the condition number
3114/// is very large, then the matrix is ill-conditioned and is almost singular.
3115/// The computation of the inverse is prone to large numerical errors
3116
3118{
3119 auto cache = getCache(_curNormSet);
3120 if (!cache)
3121 coutE(Caching) << "unable to retrieve cache!" << std::endl;
3122 return cache->_condition;
3123}
3124
3125////////////////////////////////////////////////////////////////////////////////
3126/// Return the RooRatio form of products and denominators of morphing functions
3127
3128std::unique_ptr<RooRatio>
3129RooLagrangianMorphFunc::makeRatio(const char *name, const char *title, RooArgList &nr, RooArgList &dr)
3130{
3131 RooArgList num, denom;
3132 for (auto it : nr) {
3133 num.add(*it);
3134 }
3135 for (auto it : dr) {
3136 denom.add(*it);
3137 }
3138 // same for denom
3139 return make_unique<RooRatio>(name, title, num, denom);
3140}
#define d(i)
Definition: RSha256.hxx:102
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
#define e(i)
Definition: RSha256.hxx:103
RooCollectionProxy< RooArgList > RooListProxy
Definition: RooAbsArg.h:50
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:99
#define coutW(a)
Definition: RooMsgService.h:36
#define oocxcoutP(o, a)
Definition: RooMsgService.h:95
#define coutE(a)
Definition: RooMsgService.h:37
#define cxcoutP(a)
Definition: RooMsgService.h:93
#define TRACE_DESTROY
Definition: RooTrace.h:24
#define TRACE_CREATE
Definition: RooTrace.h:23
#define ClassImp(name)
Definition: Rtypes.h:375
#define gDirectory
Definition: TDirectory.h:348
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void input
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char filename
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char cname
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void value
Option_t Option_t TPoint TPoint const char mode
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t TString Int_t GCValues_t GetPrimarySelectionOwner GetDisplay GetScreen GetColormap GetNativeEvent const char const char dpyName wid window const char font_name cursor keysym reg const char only_if_exist regb h Point_t winding char text const char depth char const char Int_t count const char ColorStruct_t color const char Pixmap_t Pixmap_t PictureAttributes_t attr const char char ret_data h unsigned char height h Atom_t Int_t ULong_t ULong_t unsigned char prop_list Atom_t Atom_t Atom_t Time_t type
char name[80]
Definition: TGX11.cxx:110
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:23
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2447
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:72
void setStringAttribute(const Text_t *key, const Text_t *value)
Associate string 'value' to this object under key 'key'.
Definition: RooAbsArg.cxx:310
bool isConstant() const
Check if the "Constant" attribute is set.
Definition: RooAbsArg.h:380
void Print(Option_t *options=0) const override
Print the object to the defaultPrintStream().
Definition: RooAbsArg.h:340
const RefCountList_t & servers() const
List of all servers of this object.
Definition: RooAbsArg.h:199
void setValueDirty()
Mark the element dirty. This forces a re-evaluation when a value is requested.
Definition: RooAbsArg.h:508
bool getAttribute(const Text_t *name) const
Check if a named attribute is set. By default, all attributes are unset.
Definition: RooAbsArg.cxx:301
virtual double * array() const =0
RooAbsCacheElement is the abstract base class for objects to be stored in RooAbsCache cache manager o...
virtual bool addOwned(RooAbsArg &var, bool silent=false)
Add an argument and transfer the ownership to the collection.
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * first() const
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 getMax(const char *name=0) const
Get maximum of currently defined range.
void setConstant(bool value=true)
void setBin(Int_t ibin, const char *rangeName=0) override
Set value to center of bin 'ibin' of binning 'rangeName' (or of default binning if no range is specif...
virtual double getMin(const char *name=0) const
Get minimum of currently defined range.
Int_t numBins(const char *rangeName=0) const override
virtual Int_t getBins(const char *name=0) const
Get number of bins of currently defined range.
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:64
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:94
virtual double getValV(const RooArgSet *normalisationSet=nullptr) const
Return value of object.
Definition: RooAbsReal.cxx:275
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:57
Class RooBinning is an implements RooAbsBinning in terms of an array of boundary values,...
Definition: RooBinning.h:27
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=0, const TNamed *isetRangeName=0)
Getter function without integration set.
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=0)
Setter function without integration set.
bool addOwned(RooAbsArg &var, bool silent=false) override
Overloaded RooCollection_t::addOwned() method insert object into owning set and registers object as s...
bool add(const RooAbsArg &var, bool valueServer, bool shapeServer, bool silent)
Overloaded RooCollection_t::add() method insert object into set and registers object as server to own...
The RooDataHist is a container class to hold N-dimensional binned data.
Definition: RooDataHist.h:45
double weight(std::size_t i) const
Return weight of i-th bin.
Definition: RooDataHist.h:110
double weightSquared(std::size_t i) const
Return squared weight sum of i-th bin.
Definition: RooDataHist.h:114
const RooArgSet * get() const override
Get bin centre of current bin.
Definition: RooDataHist.h:82
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
static TClass * Class()
static RooLagrangianMorphFunc::CacheElem * createCache(const RooLagrangianMorphFunc *func)
create all the temporary objects required by the class
void buildMatrix(const RooLagrangianMorphFunc::ParamMap &inputParameters, const RooLagrangianMorphFunc::FlagMap &inputFlags, const List &flags)
build and invert the morphing matrix
static RooLagrangianMorphFunc::CacheElem * createCache(const RooLagrangianMorphFunc *func, const Matrix &inverse)
create all the temporary objects required by the class function variant with precomputed inverse matr...
std::unique_ptr< RooRealSumFunc > _sumFunc
RooArgList containedArgs(Action) override
retrieve the list of contained args
void operModeHook(RooAbsArg::OperMode) override
Interface for changes of operation mode.
void createComponents(const RooLagrangianMorphFunc::ParamMap &inputParameters, const RooLagrangianMorphFunc::FlagMap &inputFlags, const char *funcname, const std::vector< std::vector< RooListProxy * > > &diagramProxyList, const std::vector< RooArgList * > &nonInterfering, const RooArgList &flags)
create the basic objects required for the morphing
void buildMorphingFunction(const char *name, const RooLagrangianMorphFunc::ParamMap &inputParameters, const std::map< std::string, int > &storage, const RooArgList &physics, bool allowNegativeYields, RooRealVar *observable, RooRealVar *binWidth)
build the final morphing function
Class RooLagrangianMorphing is a implementation of the method of Effective Lagrangian Morphing,...
bool isParameterConstant(const char *paramname) const
return true if the parameter with the given name is set constant, false otherwise
TPair * makeCrosssectionContainer(double xs, double unc)
create TPair containers of the type expected by the RooLagrangianMorph
bool isBinnedDistribution(const RooArgSet &obs) const override
check if this PDF is a binned distribution in the given observable
int nPolynomials() const
return the number of samples in this morphing function
void setParameter(const char *name, double value)
set one parameter to a specific value
ParamSet getMorphParameters() const
retrieve the parameter set
double evaluate() const override
call getVal on the internal function
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...
void setCacheAndTrackHints(RooArgSet &) override
Retrieve the matrix of coefficients.
ParamSet getCouplings() const
retrieve a set of couplings (-?-)
void printSampleWeights() const
print the current sample weights
std::map< const std::string, double > ParamSet
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.
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
std::list< double > * plotSamplingHint(RooAbsRealLValue &, double, double) const override
retrieve the sample Hint
RooRealVar * getFlag(const char *name) const
retrieve the RooRealVar object incorporating the flag with the given name
void randomizeParameters(double z)
randomize the parameters a bit useful to test and debug fitting
bool isCouplingUsed(const char *couplname)
check if there is any morphing power provided for the given coupling morphing power is provided as so...
void readParameters(TDirectory *f)
read the parameters from the input file
double getScale()
get energy scale of the EFT expansion
double getCondition() const
Retrieve the condition of the coefficient matrix.
double analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=0) const override
Retrieve the matrix of coefficients.
TMatrixD getMatrix() const
Retrieve the matrix of coefficients.
void printWeights() const
print the current sample weights
void printCouplings() const
print a set of couplings
TMatrixD readMatrixFromFile(const char *fname)
read a matrix from a text file
~RooLagrangianMorphFunc() override
default destructor
void printParameters() const
print the parameters and their current values
void printPhysics() const
print the current physics values
static std::unique_ptr< RooRatio > makeRatio(const char *name, const char *title, RooArgList &nr, RooArgList &dr)
Return the RooRatio form of products and denominators of morphing functions.
void setFlag(const char *name, double value)
set one flag to a specific value
TH1 * createTH1(const std::string &name)
retrieve a histogram output of the current morphing settings
double expectedUncertainty() const
return the expected uncertainty for the current parameter set
int nParameters() const
return the number of parameters in this morphing function
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
bool checkObservables(const RooArgSet *nset) const override
check if observable exists in the RooArgSet (-?-)
bool hasCache() const
return true if a cache object is present, false otherwise
void printFlags() const
print the flags and their current values
void setScale(double val)
set energy scale of the EFT expansion
TMatrixD getInvertedMatrix() const
Retrieve the matrix of coefficients after inversion.
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
bool forceAnalyticalInt(const RooAbsArg &arg) const override
Force analytical integration for the given observable.
void setParameterConstant(const char *paramname, bool constant) const
call setConstant with the boolean argument provided on the parameter with the given name
void printSamples() const
print all the known samples to the console
double getParameterValue(const char *name) const
set one parameter to a specific value
double getValV(const RooArgSet *set=0) const override
call getVal on the internal function
void setup(bool ownParams=true)
setup this instance with the given set of operators and vertices if own=true, the class will own the ...
void printMetaArgs(std::ostream &os) const override
Retrieve the matrix of coefficients.
std::unique_ptr< RooWrapperPdf > createPdf() const
(currently similar to cloning the Pdf
RooAbsReal * getSampleWeight(const char *name)
retrieve the weight (prefactor) of a sample with the given name
std::map< std::string, std::string > createWeightStrings(const ParamMap &inputs, const std::vector< std::vector< std::string > > &vertices)
create only the weight formulas. static function for external usage.
std::vector< std::vector< RooListProxy * > > _diagrams
double expectedEvents() const
return the number of expected events for the current parameter set
std::map< std::string, int > _sampleMap
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
RooRealSumFunc * getFunc() const
get the func
std::list< double > * binBoundaries(RooAbsRealLValue &, double, double) const override
retrieve the list of bin boundaries
void printEvaluation() const
print the contributing samples and their respective weights
bool writeCoefficients(const char *filename)
write the inverse matrix to a file
void collectInputs(TDirectory *f)
retrieve the physics inputs
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)
A histogram function that assigns scale parameters to every bin.
static TClass * Class()
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:48
std::list< double > * binBoundaries(RooAbsRealLValue &, double, double) const override
Retrieve bin boundaries if this distribution is binned in obs.
void printMetaArgs(std::ostream &os) const override
Customized printing of arguments of a RooRealSumFunc to more intuitively reflect the contents of the ...
Int_t getAnalyticalIntegralWN(RooArgSet &allVars, RooArgSet &numVars, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Variant of getAnalyticalIntegral that is also passed the normalization set that should be applied to ...
CacheMode canNodeBeCached() const override
double analyticalIntegralWN(Int_t code, const RooArgSet *normSet, const char *rangeName=nullptr) const override
Implements the actual analytical integral(s) advertised by getAnalyticalIntegral.
bool isBinnedDistribution(const RooArgSet &obs) const override
Tests if the distribution is binned. Unless overridden by derived classes, this always returns false.
bool checkObservables(const RooArgSet *nset) const override
Overloadable function in which derived classes can implement consistency checks of the variables.
void setCacheAndTrackHints(RooArgSet &) override
std::list< double > * plotSamplingHint(RooAbsRealLValue &, double, double) const override
Interface for returning an optional hint for initial sampling points when constructing a curve projec...
bool forceAnalyticalInt(const RooAbsArg &arg) const override
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
void setVal(double value) override
Set value of variable to 'value'.
Definition: RooRealVar.cxx:281
void setError(double value)
Definition: RooRealVar.h:65
void setBins(Int_t nBins, const char *name=0)
Create a uniform binning under name 'name' for this variable.
Definition: RooRealVar.cxx:434
void setMin(const char *name, double value)
Set minimum of name range to given value.
Definition: RooRealVar.cxx:487
double getError() const
Definition: RooRealVar.h:63
const RooAbsBinning & getBinning(const char *name=0, bool verbose=true, bool createOnTheFly=false) const override
Return binning definition with name.
Definition: RooRealVar.cxx:344
void setMax(const char *name, double value)
Set maximum of name range to given value.
Definition: RooRealVar.cxx:517
void setBinning(const RooAbsBinning &binning, const char *name=0)
Add given binning under name 'name' with this variable.
Definition: RooRealVar.cxx:441
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:135
const char * GetBinLabel(Int_t bin) const
Return label for bin.
Definition: TAxis.cxx:441
Double_t GetXmin() const
Definition: TAxis.h:134
Int_t GetNbins() const
Definition: TAxis.h:121
TClass instances represent classes, structs and namespaces in the ROOT type system.
Definition: TClass.h:80
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:2968
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:4019
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:481
1-D histogram with a float per channel (see TH1 documentation)}
Definition: TH1.h:574
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
static TClass * Class()
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:319
virtual Int_t GetNbinsX() const
Definition: TH1.h:295
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:9083
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:9099
virtual Double_t GetBinLowEdge(Int_t bin) const
Return bin lower edge for 1D histogram.
Definition: TH1.cxx:9029
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:5035
virtual Double_t GetBinWidth(Int_t bin) const
Return bin width for 1D histogram.
Definition: TH1.cxx:9040
Int_t GetNrows() const
Definition: TMatrixTBase.h:123
virtual TMatrixTBase< Element > & UnitMatrix()
Make a unit matrix (matrix need not be a square one).
TMatrixTBase< Element > & ResizeTo(Int_t nrows, Int_t ncols, Int_t=-1) override
Set size of the matrix to nrows x ncols New dynamic elements are created, the overlapping part of the...
Definition: TMatrixT.cxx:1211
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
const char * GetTitle() const override
Returns title of object.
Definition: TNamed.h:48
Mother of all ROOT objects.
Definition: TObject.h:37
virtual const char * GetName() const
Returns name of object.
Definition: TObject.cxx:359
virtual const char * ClassName() const
Returns name of class to which the object belongs.
Definition: TObject.cxx:130
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
TString & Append(const char *cs)
Definition: TString.h:564
TLine * line
RVec< PromoteTypes< T0, T1 > > pow(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1753
const Int_t n
Definition: legend1.C:16
#define T2
Definition: md5.inl:147
#define T1
Definition: md5.inl:146
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)
@ InputArguments
Definition: RooGlobalFunc.h:64
@ ObjectHandling
Definition: RooGlobalFunc.h:64
Definition: file.py:1
Definition: first.py:1
const char * Ref
Definition: TXMLSetup.cxx:53
std::vector< std::string > folderNames
Definition: civetweb.c:1856
TMarker m
Definition: textangle.C:8
TArc a
Definition: textangle.C:12
double * vertex
Definition: triangle.c:513