Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooMomentMorphFunc.cxx
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 *
4 * Copyright (c) 2023, CERN
5 *
6 * Redistribution and use in source and binary forms,
7 * with or without modification, are permitted according to the terms
8 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
9 */
10
11/** \class RooMomentMorphFunc
12 \ingroup Roofit
13
14**/
15
16#include "Riostream.h"
17
18#include "RooMomentMorphFunc.h"
19#include "RooAbsCategory.h"
20#include "RooRealConstant.h"
21#include "RooRealVar.h"
22#include "RooFormulaVar.h"
23#include "RooCustomizer.h"
24#include "RooRealSumFunc.h"
25#include "RooAddition.h"
26#include "RooMoment.h"
27#include "RooLinearVar.h"
28#include "RooChangeTracker.h"
29
30#include "TMath.h"
31#include "TH1.h"
32
33using std::cout, std::endl, std::string, std::vector;
34
36
37//_____________________________________________________________________________
39 : _cacheMgr(this, 10, true, true)
40{
41}
42
43//_____________________________________________________________________________
44RooMomentMorphFunc::RooMomentMorphFunc(const char *name, const char *title, RooAbsReal &_m, const RooArgList &varList,
45 const RooArgList &pdfList, const TVectorD &mrefpoints, Setting setting)
46 : RooAbsReal(name, title),
47 _cacheMgr(this, 10, true, true),
48 m("m", "m", this, _m),
49 _varList("varList", "List of variables", this),
50 _pdfList("pdfList", "List of pdfs", this),
51 _mref(new TVectorD(mrefpoints)),
52 _setting(setting)
53{
54 // observables
56
57 // reference p.d.f.s
58 _pdfList.addTyped<RooAbsPdf>(pdfList);
59
60 // initialization
61 initialize();
62}
63
64//_____________________________________________________________________________
65RooMomentMorphFunc::RooMomentMorphFunc(const char *name, const char *title, RooAbsReal &_m, const RooArgList &varList,
66 const RooArgList &pdfList, const RooArgList &mrefList, Setting setting)
67 : RooAbsReal(name, title),
68 _cacheMgr(this, 10, true, true),
69 m("m", "m", this, _m),
70 _varList("varList", "List of variables", this),
71 _pdfList("pdfList", "List of pdfs", this),
72 _mref(new TVectorD(mrefList.size())),
73 _setting(setting)
74{
75 // observables
77
78 // reference p.d.f.s
79 _pdfList.addTyped<RooAbsPdf>(pdfList);
80
81 // reference points in m
82
83 Int_t i = 0;
84 for (auto *mref : mrefList) {
85 if (!dynamic_cast<RooAbsReal *>(mref)) {
86 coutE(InputArguments) << "RooMomentMorphFunc::ctor(" << GetName() << ") ERROR: mref " << mref->GetName()
87 << " is not of type RooAbsReal" << endl;
88 throw string("RooPolyMorh::ctor() ERROR mref is not of type RooAbsReal");
89 }
90 if (!dynamic_cast<RooConstVar *>(mref)) {
91 coutW(InputArguments) << "RooMomentMorphFunc::ctor(" << GetName() << ") WARNING mref point " << i
92 << " is not a constant, taking a snapshot of its value" << endl;
93 }
94 (*_mref)[i] = static_cast<RooAbsReal *>(mref)->getVal();
95 ++i;
96 }
97
98 // initialization
99 initialize();
100}
101
102//_____________________________________________________________________________
104 : RooAbsReal(other, name),
105 _cacheMgr(other._cacheMgr, this),
106 m("m", this, other.m),
107 _varList("varList", this, other._varList),
108 _pdfList("pdfList", this, other._pdfList),
109 _mref(new TVectorD(*other._mref)),
110 _setting(other._setting),
111 _useHorizMorph(other._useHorizMorph)
112{
113
114 // initialization
115 initialize();
116}
117
118//_____________________________________________________________________________
120{
121 if (_mref)
122 delete _mref;
123 if (_M)
124 delete _M;
125}
126
127//_____________________________________________________________________________
129{
130
131 Int_t nPdf = _pdfList.size();
132
133 // other quantities needed
134 if (nPdf != _mref->GetNrows()) {
135 coutE(InputArguments) << "RooMomentMorphFunc::initialize(" << GetName() << ") ERROR: nPdf != nRefPoints" << endl;
136 assert(0);
137 }
138
139 TVectorD *dm = new TVectorD(nPdf);
140 _M = new TMatrixD(nPdf, nPdf);
141
142 // transformation matrix for non-linear extrapolation, needed in evaluate()
143 TMatrixD M(nPdf, nPdf);
144 for (Int_t i = 0; i < _mref->GetNrows(); ++i) {
145 (*dm)[i] = (*_mref)[i] - (*_mref)[0];
146 M(i, 0) = 1.;
147 if (i > 0)
148 M(0, i) = 0.;
149 }
150 for (Int_t i = 1; i < _mref->GetNrows(); ++i) {
151 for (Int_t j = 1; j < _mref->GetNrows(); ++j) {
152 M(i, j) = std::pow((*dm)[i], (double)j);
153 }
154 }
155 (*_M) = M.Invert();
156
157 delete dm;
158}
159
160//_____________________________________________________________________________
162{
163 auto cache = static_cast<CacheElem *>(_cacheMgr.getObj(nullptr, static_cast<RooArgSet const*>(nullptr)));
164 if (cache) {
165 return cache;
166 }
167 Int_t nVar = _varList.size();
168 Int_t nPdf = _pdfList.size();
169
170 RooAbsReal *null = nullptr;
171 vector<RooAbsReal *> meanrv(nPdf * nVar, null);
172 vector<RooAbsReal *> sigmarv(nPdf * nVar, null);
173 vector<RooAbsReal *> myrms(nVar, null);
174 vector<RooAbsReal *> mypos(nVar, null);
175 vector<RooAbsReal *> slope(nPdf * nVar, null);
176 vector<RooAbsReal *> offs(nPdf * nVar, null);
177 vector<RooAbsReal *> transVar(nPdf * nVar, null);
178 vector<RooAbsReal *> transPdf(nPdf, null);
179
180 RooArgSet ownedComps;
181
182 RooArgList fracl;
183
184 // fraction parameters
185 RooArgList coefList("coefList");
186 RooArgList coefList2("coefList2");
187 for (Int_t i = 0; i < 2 * nPdf; ++i) {
188 std::string fracName = Form("frac_%d", i);
189
190 RooRealVar *frac = new RooRealVar(fracName.c_str(), fracName.c_str(), 1.);
191
192 fracl.add(*frac); // to be set later
193 if (i < nPdf) {
194 coefList.add(*static_cast<RooRealVar *>(fracl.at(i)));
195 } else {
196 coefList2.add(*static_cast<RooRealVar *>(fracl.at(i)));
197 }
198 ownedComps.add(*static_cast<RooRealVar *>(fracl.at(i)));
199 }
200
201 RooRealSumFunc *theSumFunc = nullptr;
202 std::string sumfuncName = Form("%s_sumfunc", GetName());
203
204 if (_useHorizMorph) {
205 // mean and sigma
206 RooArgList varList(_varList);
207 for (Int_t i = 0; i < nPdf; ++i) {
208 for (Int_t j = 0; j < nVar; ++j) {
209
210 std::string meanName = Form("%s_mean_%d_%d", GetName(), i, j);
211 std::string sigmaName = Form("%s_sigma_%d_%d", GetName(), i, j);
212
213 RooAbsMoment *mom = nVar == 1 ? (static_cast<RooAbsPdf *>(_pdfList.at(i)))->sigma(static_cast<RooRealVar &>(*varList.at(j)))
214 : (static_cast<RooAbsPdf *>(_pdfList.at(i)))->sigma(static_cast<RooRealVar &>(*varList.at(j)), varList);
215
216 mom->setLocalNoDirtyInhibit(true);
217 mom->mean()->setLocalNoDirtyInhibit(true);
218
219 sigmarv[ij(i, j)] = mom;
220 meanrv[ij(i, j)] = mom->mean();
221
222 ownedComps.add(*sigmarv[ij(i, j)]);
223 }
224 }
225
226 // slope and offset (to be set later, depend on m)
227 for (Int_t j = 0; j < nVar; ++j) {
228 RooArgList meanList("meanList");
229 RooArgList rmsList("rmsList");
230 for (Int_t i = 0; i < nPdf; ++i) {
231 meanList.add(*meanrv[ij(i, j)]);
232 rmsList.add(*sigmarv[ij(i, j)]);
233 }
234 std::string myrmsName = Form("%s_rms_%d", GetName(), j);
235 std::string myposName = Form("%s_pos_%d", GetName(), j);
236 myrms[j] = new RooAddition(myrmsName.c_str(), myrmsName.c_str(), rmsList, coefList2);
237 mypos[j] = new RooAddition(myposName.c_str(), myposName.c_str(), meanList, coefList2);
238 ownedComps.add(RooArgSet(*myrms[j], *mypos[j]));
239 }
240
241 // construction of unit pdfs
242 RooArgList transPdfList;
243
244 for (Int_t i = 0; i < nPdf; ++i) {
245 auto& pdf = static_cast<RooAbsPdf&>(_pdfList[i]);
246 std::string pdfName = Form("pdf_%d", i);
247 RooCustomizer cust(pdf, pdfName.c_str());
248
249 for (Int_t j = 0; j < nVar; ++j) {
250 // slope and offset formulas
251 std::string slopeName = Form("%s_slope_%d_%d", GetName(), i, j);
252 std::string offsetName = Form("%s_offset_%d_%d", GetName(), i, j);
253 slope[ij(i, j)] = new RooFormulaVar(slopeName.c_str(), "@0/@1", RooArgList(*sigmarv[ij(i, j)], *myrms[j]));
254 offs[ij(i, j)] = new RooFormulaVar(offsetName.c_str(), "@0-(@1*@2)",
255 RooArgList(*meanrv[ij(i, j)], *mypos[j], *slope[ij(i, j)]));
256 ownedComps.add(RooArgSet(*slope[ij(i, j)], *offs[ij(i, j)]));
257 // linear transformations, so pdf can be renormalized
258 auto& var = static_cast<RooRealVar&>(*_varList[j]);
259 std::string transVarName = Form("%s_transVar_%d_%d", GetName(), i, j);
260 // transVar[ij(i,j)] = new
261 // RooFormulaVar(transVarName.c_str(),transVarName.c_str(),"@0*@1+@2",RooArgList(*var,*slope[ij(i,j)],*offs[ij(i,j)]));
262
263 transVar[ij(i, j)] =
264 new RooLinearVar(transVarName.c_str(), transVarName.c_str(), var, *slope[ij(i, j)], *offs[ij(i, j)]);
265
266 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
267 // This will prevent the likelihood optimizers from erroneously declaring terms constant
268 transVar[ij(i, j)]->addServer((RooAbsArg &)m.arg());
269
270 ownedComps.add(*transVar[ij(i, j)]);
271 cust.replaceArg(var, *transVar[ij(i, j)]);
272 }
273 transPdf[i] = static_cast<RooAbsPdf *>(cust.build());
274 transPdfList.add(*transPdf[i]);
275 ownedComps.add(*transPdf[i]);
276 }
277 // sum pdf
278 theSumFunc = new RooRealSumFunc(sumfuncName.c_str(), sumfuncName.c_str(), transPdfList, coefList);
279 } else {
280 theSumFunc = new RooRealSumFunc(sumfuncName.c_str(), sumfuncName.c_str(), _pdfList, coefList);
281 }
282
283 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
284 // This will prevent the likelihood optimizers from erroneously declaring terms constant
285 theSumFunc->addServer((RooAbsArg &)m.arg());
286 theSumFunc->addOwnedComponents(ownedComps);
287
288 // change tracker for fraction parameters
289 std::string trackerName = Form("%s_frac_tracker", GetName());
290 RooChangeTracker *tracker = new RooChangeTracker(trackerName.c_str(), trackerName.c_str(), m.arg(), true);
291
292 // Store it in the cache
293 cache = new CacheElem(*theSumFunc, *tracker, fracl);
294 _cacheMgr.setObj(nullptr, nullptr, cache, nullptr);
295
296 return cache;
297}
298
299//_____________________________________________________________________________
301{
302 return RooArgList(*_sumFunc, *_tracker);
303}
304
305//_____________________________________________________________________________
307{
308 delete _sumFunc;
309 delete _tracker;
310}
311
312//_____________________________________________________________________________
314{
315 // Special version of getValV() overrides RooAbsReal::getVal() to save value of current normalization set
316 _curNormSet = set ? const_cast<RooArgSet *>(set) : const_cast<RooArgSet *>(static_cast<RooArgSet const*>(&_varList));
317 return RooAbsReal::getValV(set);
318}
319
320//_____________________________________________________________________________
322{
323 CacheElem *cache = getCache(nset ? nset : _curNormSet);
324
325 if (cache->_tracker->hasChanged(true)) {
326 cache->calculateFractions(*this, false); // verbose turned off
327 }
328
329 return cache->_sumFunc;
330}
331
332//_____________________________________________________________________________
334{
335 CacheElem *cache = getCache(nset ? nset : _curNormSet);
336
337 if (cache->_tracker->hasChanged(true)) {
338 cache->calculateFractions(*this, false); // verbose turned off
339 }
340
341 return cache->_sumFunc;
342}
343
344//_____________________________________________________________________________
346{
348
349 if (cache->_tracker->hasChanged(true)) {
350 cache->calculateFractions(*this, false); // verbose turned off
351 }
352
353 double ret = cache->_sumFunc->getVal(_pdfList.nset());
354 return ret;
355}
356
357//_____________________________________________________________________________
359{
360 return static_cast<RooRealVar *>(_frac.at(i));
361}
362
363//_____________________________________________________________________________
365{
366 return static_cast<RooRealVar *>(_frac.at(i));
367}
368
369//_____________________________________________________________________________
371{
372 Int_t nPdf = self._pdfList.size();
373
374 double dm = self.m - (*self._mref)[0];
375
376 // fully non-linear
377 double sumposfrac = 0.;
378 for (Int_t i = 0; i < nPdf; ++i) {
379 double ffrac = 0.;
380 for (Int_t j = 0; j < nPdf; ++j) {
381 ffrac += (*self._M)(j, i) * (j == 0 ? 1. : std::pow(dm, (double)j));
382 }
383 if (ffrac >= 0)
384 sumposfrac += ffrac;
385 // fractions for pdf
386 const_cast<RooRealVar *>(frac(i))->setVal(ffrac);
387 // fractions for rms and mean
388 const_cast<RooRealVar *>(frac(nPdf + i))->setVal(ffrac);
389 if (verbose) {
390 cout << ffrac << endl;
391 }
392 }
393
394 // various mode settings
395 int imin = self.idxmin(self.m);
396 int imax = self.idxmax(self.m);
397 double mfrac = (self.m - (*self._mref)[imin]) / ((*self._mref)[imax] - (*self._mref)[imin]);
398 switch (self._setting) {
399 case NonLinear:
400 // default already set above
401 break;
402
403 case SineLinear:
404 mfrac =
405 std::sin(TMath::PiOver2() * mfrac); // this gives a continuous differentiable transition between grid points.
406
407 // now fall through to Linear case
408
409 case Linear:
410 for (Int_t i = 0; i < 2 * nPdf; ++i) const_cast<RooRealVar *>(frac(i))->setVal(0.);
411 if (imax > imin) { // m in between mmin and mmax
412 const_cast<RooRealVar *>(frac(imin))->setVal(1. - mfrac);
413 const_cast<RooRealVar *>(frac(nPdf + imin))->setVal(1. - mfrac);
414 const_cast<RooRealVar *>(frac(imax))->setVal(mfrac);
415 const_cast<RooRealVar *>(frac(nPdf + imax))->setVal(mfrac);
416 } else if (imax == imin) { // m outside mmin and mmax
417 const_cast<RooRealVar *>(frac(imin))->setVal(1.);
418 const_cast<RooRealVar *>(frac(nPdf + imin))->setVal(1.);
419 }
420 break;
422 for (Int_t i = 0; i < nPdf; ++i) const_cast<RooRealVar *>(frac(i))->setVal(0.);
423 if (imax > imin) { // m in between mmin and mmax
424 const_cast<RooRealVar *>(frac(imin))->setVal(1. - mfrac);
425 const_cast<RooRealVar *>(frac(imax))->setVal(mfrac);
426 } else if (imax == imin) { // m outside mmin and mmax
427 const_cast<RooRealVar *>(frac(imin))->setVal(1.);
428 }
429 break;
431 for (Int_t i = 0; i < nPdf; ++i) {
432 if (frac(i)->getVal() < 0)
433 const_cast<RooRealVar *>(frac(i))->setVal(0.);
434 const_cast<RooRealVar *>(frac(i))->setVal(frac(i)->getVal() / sumposfrac);
435 }
436 break;
437 }
438}
439
440//_____________________________________________________________________________
441int RooMomentMorphFunc::idxmin(const double &mval) const
442{
443 int imin(0);
444 Int_t nPdf = _pdfList.size();
445 double mmin = -DBL_MAX;
446 for (Int_t i = 0; i < nPdf; ++i) {
447 if ((*_mref)[i] > mmin && (*_mref)[i] <= mval) {
448 mmin = (*_mref)[i];
449 imin = i;
450 }
451 }
452 return imin;
453}
454
455//_____________________________________________________________________________
456int RooMomentMorphFunc::idxmax(const double &mval) const
457{
458 int imax(0);
459 Int_t nPdf = _pdfList.size();
460 double mmax = DBL_MAX;
461 for (Int_t i = 0; i < nPdf; ++i) {
462 if ((*_mref)[i] < mmax && (*_mref)[i] >= mval) {
463 mmax = (*_mref)[i];
464 imax = i;
465 }
466 }
467 return imax;
468}
469
470//_____________________________________________________________________________
471std::list<double> *RooMomentMorphFunc::plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const
472{
473 return sumFunc(nullptr)->plotSamplingHint(obs, xlo, xhi);
474}
475
476//_____________________________________________________________________________
477std::list<double> *RooMomentMorphFunc::binBoundaries(RooAbsRealLValue &obs, double xlo, double xhi) const
478{
479 return sumFunc(nullptr)->binBoundaries(obs, xlo, xhi);
480}
481
482//_____________________________________________________________________________
484{
485 return sumFunc(nullptr)->isBinnedDistribution(obs);
486}
size_t size(const MatrixT &matrix)
retrieve the size of a square matrix
#define coutW(a)
#define coutE(a)
#define ClassImp(name)
Definition Rtypes.h:382
char name[80]
Definition TGX11.cxx:110
TMatrixT< Double_t > TMatrixD
Definition TMatrixDfwd.h:23
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
TVectorT< Double_t > TVectorD
Definition TVectorDfwd.h:23
Common abstract base class for objects that represent a value and a "shape" in RooFit.
Definition RooAbsArg.h:79
void setLocalNoDirtyInhibit(bool flag) const
Definition RooAbsArg.h:674
R__DEPRECATED(6, 36, "Use getObservables().") RooFit R__DEPRECATED(6, 36, "Use getObservables().") RooFit R__DEPRECATED(6, 36, "Use getObservables().") RooFit const RooAbsArg &testArg const
Definition RooAbsArg.h:145
bool addOwnedComponents(const RooAbsCollection &comps)
Take ownership of the contents of 'comps'.
void addServer(RooAbsArg &server, bool valueProp=true, bool shapeProp=false, std::size_t refCount=1)
Register another RooAbsArg as a server to us, ie, declare that we depend on it.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
bool addTyped(const RooAbsCollection &list, bool silent=false)
Adds elements of a given RooAbsCollection to the container if they match the specified type.
RooAbsReal * mean()
Abstract interface for all probability density functions.
Definition RooAbsPdf.h:40
const RooArgSet * nset() const
Definition RooAbsProxy.h:52
Abstract base class for objects that represent a real value that may appear on the left hand side of ...
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
virtual std::list< double > * binBoundaries(RooAbsRealLValue &obs, double xlo, double xhi) const
Retrieve bin boundaries if this distribution is binned in obs.
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
virtual double getValV(const RooArgSet *normalisationSet=nullptr) const
Return value of object.
friend class RooRealSumFunc
Definition RooAbsReal.h:403
virtual bool isBinnedDistribution(const RooArgSet &) const
Tests if the distribution is binned. Unless overridden by derived classes, this always returns false.
Definition RooAbsReal.h:346
virtual std::list< double > * plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const
Interface for returning an optional hint for initial sampling points when constructing a curve projec...
RooAbsMoment * sigma(RooRealVar &obs)
Definition RooAbsReal.h:361
Calculates the sum of a set of RooAbsReal terms, or when constructed with two sets,...
Definition RooAddition.h:27
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:24
Int_t setObj(const RooArgSet *nset, T *obj, const TNamed *isetRangeName=nullptr)
Setter function without integration set.
T * getObj(const RooArgSet *nset, Int_t *sterileIndex=nullptr, const TNamed *isetRangeName=nullptr)
Getter function without integration set.
Meta object that tracks value changes in a given set of RooAbsArgs by registering itself as value cli...
bool hasChanged(bool clearState)
Returns true if state has changed since last call with clearState=true.
Represents a constant real-valued object.
Definition RooConstVar.h:23
RooCustomizer is a factory class to produce clones of a prototype composite PDF object with the same ...
void replaceArg(const RooAbsArg &orig, const RooAbsArg &subst)
Replace any occurrence of arg 'orig' with arg 'subst'.
RooAbsArg * build(const char *masterCatState, bool verbose=false)
Build a clone of the prototype executing all registered 'replace' rules and 'split' rules for the mas...
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
RooLinearVar is the most general form of a derived real-valued object that can be used by RooRealInte...
void calculateFractions(const RooMomentMorphFunc &self, bool verbose=true) const
RooArgList containedArgs(Action) override
double getValV(const RooArgSet *set=nullptr) const override
Return value of object.
RooArgSet * _curNormSet
The cache manager.
bool isBinnedDistribution(const RooArgSet &obs) const override
Tests if the distribution is binned. Unless overridden by derived classes, this always returns false.
RooObjCacheManager _cacheMgr
CacheElem * getCache(const RooArgSet *nset) const
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooAbsReal * sumFunc(const RooArgSet *nset)
int idxmin(const double &m) const
std::list< double > * binBoundaries(RooAbsRealLValue &, double, double) const override
Retrieve bin boundaries if this distribution is binned in obs.
std::list< double > * plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const override
Interface for returning an optional hint for initial sampling points when constructing a curve projec...
Int_t ij(const Int_t &i, const Int_t &j) const
int idxmax(const double &m) const
Variable that can be changed from the outside.
Definition RooRealVar.h:37
TMatrixT< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant.
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
Int_t GetNrows() const
Definition TVectorT.h:73
constexpr Double_t PiOver2()
Definition TMath.h:51
TMarker m
Definition textangle.C:8