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