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