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