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 namespace std;
34
36
37//_____________________________________________________________________________
39 : _cacheMgr(this, 10, true, true), _curNormSet(nullptr), _mref(nullptr), _M(nullptr), _useHorizMorph(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 _useHorizMorph(true)
54{
55 // observables
57
58 // reference p.d.f.s
59 _pdfList.addTyped<RooAbsPdf>(pdfList);
60
61 // initialization
62 initialize();
63}
64
65//_____________________________________________________________________________
66RooMomentMorphFunc::RooMomentMorphFunc(const char *name, const char *title, RooAbsReal &_m, const RooArgList &varList,
67 const RooArgList &pdfList, const RooArgList &mrefList, Setting setting)
68 : RooAbsReal(name, title),
69 _cacheMgr(this, 10, true, true),
70 m("m", "m", this, _m),
71 _varList("varList", "List of variables", this),
72 _pdfList("pdfList", "List of pdfs", this),
73 _mref(new TVectorD(mrefList.size())),
74 _setting(setting),
75 _useHorizMorph(true)
76{
77 // observables
79
80 // reference p.d.f.s
81 _pdfList.addTyped<RooAbsPdf>(pdfList);
82
83 // reference points in m
84
85 Int_t i = 0;
86 for (auto *mref : mrefList) {
87 if (!dynamic_cast<RooAbsReal *>(mref)) {
88 coutE(InputArguments) << "RooMomentMorphFunc::ctor(" << GetName() << ") ERROR: mref " << mref->GetName()
89 << " is not of type RooAbsReal" << endl;
90 throw string("RooPolyMorh::ctor() ERROR mref is not of type RooAbsReal");
91 }
92 if (!dynamic_cast<RooConstVar *>(mref)) {
93 coutW(InputArguments) << "RooMomentMorphFunc::ctor(" << GetName() << ") WARNING mref point " << i
94 << " is not a constant, taking a snapshot of its value" << endl;
95 }
96 (*_mref)[i] = static_cast<RooAbsReal *>(mref)->getVal();
97 ++i;
98 }
99
100 // initialization
101 initialize();
102}
103
104//_____________________________________________________________________________
106 : RooAbsReal(other, name),
107 _cacheMgr(other._cacheMgr, this),
108 _curNormSet(nullptr),
109 m("m", this, other.m),
110 _varList("varList", this, other._varList),
111 _pdfList("pdfList", this, other._pdfList),
112 _mref(new TVectorD(*other._mref)),
113 _setting(other._setting),
114 _useHorizMorph(other._useHorizMorph)
115{
116
117 // initialization
118 initialize();
119}
120
121//_____________________________________________________________________________
123{
124 if (_mref)
125 delete _mref;
126 if (_M)
127 delete _M;
128}
129
130//_____________________________________________________________________________
132{
133
134 Int_t nPdf = _pdfList.size();
135
136 // other quantities needed
137 if (nPdf != _mref->GetNrows()) {
138 coutE(InputArguments) << "RooMomentMorphFunc::initialize(" << GetName() << ") ERROR: nPdf != nRefPoints" << endl;
139 assert(0);
140 }
141
142 TVectorD *dm = new TVectorD(nPdf);
143 _M = new TMatrixD(nPdf, nPdf);
144
145 // transformation matrix for non-linear extrapolation, needed in evaluate()
146 TMatrixD M(nPdf, nPdf);
147 for (Int_t i = 0; i < _mref->GetNrows(); ++i) {
148 (*dm)[i] = (*_mref)[i] - (*_mref)[0];
149 M(i, 0) = 1.;
150 if (i > 0)
151 M(0, i) = 0.;
152 }
153 for (Int_t i = 1; i < _mref->GetNrows(); ++i) {
154 for (Int_t j = 1; j < _mref->GetNrows(); ++j) {
155 M(i, j) = TMath::Power((*dm)[i], (double)j);
156 }
157 }
158 (*_M) = M.Invert();
159
160 delete dm;
161}
162
163//_____________________________________________________________________________
165{
166 auto cache = static_cast<CacheElem *>(_cacheMgr.getObj(nullptr, static_cast<RooArgSet const*>(nullptr)));
167 if (cache) {
168 return cache;
169 }
170 Int_t nVar = _varList.size();
171 Int_t nPdf = _pdfList.size();
172
173 RooAbsReal *null = nullptr;
174 vector<RooAbsReal *> meanrv(nPdf * nVar, null);
175 vector<RooAbsReal *> sigmarv(nPdf * nVar, null);
176 vector<RooAbsReal *> myrms(nVar, null);
177 vector<RooAbsReal *> mypos(nVar, null);
178 vector<RooAbsReal *> slope(nPdf * nVar, null);
179 vector<RooAbsReal *> offs(nPdf * nVar, null);
180 vector<RooAbsReal *> transVar(nPdf * nVar, null);
181 vector<RooAbsReal *> transPdf(nPdf, null);
182
183 RooArgSet ownedComps;
184
185 RooArgList fracl;
186
187 // fraction parameters
188 RooArgList coefList("coefList");
189 RooArgList coefList2("coefList2");
190 for (Int_t i = 0; i < 2 * nPdf; ++i) {
191 std::string fracName = Form("frac_%d", i);
192
193 RooRealVar *frac = new RooRealVar(fracName.c_str(), fracName.c_str(), 1.);
194
195 fracl.add(*frac); // to be set later
196 if (i < nPdf) {
197 coefList.add(*static_cast<RooRealVar *>(fracl.at(i)));
198 } else {
199 coefList2.add(*static_cast<RooRealVar *>(fracl.at(i)));
200 }
201 ownedComps.add(*static_cast<RooRealVar *>(fracl.at(i)));
202 }
203
204 RooRealSumFunc *theSumFunc = nullptr;
205 std::string sumfuncName = Form("%s_sumfunc", GetName());
206
207 if (_useHorizMorph) {
208 // mean and sigma
209 RooArgList varList(_varList);
210 for (Int_t i = 0; i < nPdf; ++i) {
211 for (Int_t j = 0; j < nVar; ++j) {
212
213 std::string meanName = Form("%s_mean_%d_%d", GetName(), i, j);
214 std::string sigmaName = Form("%s_sigma_%d_%d", GetName(), i, j);
215
216 RooAbsMoment *mom = nVar == 1 ? (static_cast<RooAbsPdf *>(_pdfList.at(i)))->sigma(static_cast<RooRealVar &>(*varList.at(j)))
217 : (static_cast<RooAbsPdf *>(_pdfList.at(i)))->sigma(static_cast<RooRealVar &>(*varList.at(j)), varList);
218
219 mom->setLocalNoDirtyInhibit(true);
220 mom->mean()->setLocalNoDirtyInhibit(true);
221
222 sigmarv[ij(i, j)] = mom;
223 meanrv[ij(i, j)] = mom->mean();
224
225 ownedComps.add(*sigmarv[ij(i, j)]);
226 }
227 }
228
229 // slope and offset (to be set later, depend on m)
230 for (Int_t j = 0; j < nVar; ++j) {
231 RooArgList meanList("meanList");
232 RooArgList rmsList("rmsList");
233 for (Int_t i = 0; i < nPdf; ++i) {
234 meanList.add(*meanrv[ij(i, j)]);
235 rmsList.add(*sigmarv[ij(i, j)]);
236 }
237 std::string myrmsName = Form("%s_rms_%d", GetName(), j);
238 std::string myposName = Form("%s_pos_%d", GetName(), j);
239 myrms[j] = new RooAddition(myrmsName.c_str(), myrmsName.c_str(), rmsList, coefList2);
240 mypos[j] = new RooAddition(myposName.c_str(), myposName.c_str(), meanList, coefList2);
241 ownedComps.add(RooArgSet(*myrms[j], *mypos[j]));
242 }
243
244 // construction of unit pdfs
245 RooArgList transPdfList;
246
247 for (Int_t i = 0; i < nPdf; ++i) {
248 auto& pdf = static_cast<RooAbsPdf&>(_pdfList[i]);
249 std::string pdfName = Form("pdf_%d", i);
250 RooCustomizer cust(pdf, pdfName.c_str());
251
252 for (Int_t j = 0; j < nVar; ++j) {
253 // slope and offset formulas
254 std::string slopeName = Form("%s_slope_%d_%d", GetName(), i, j);
255 std::string offsetName = Form("%s_offset_%d_%d", GetName(), i, j);
256 slope[ij(i, j)] = new RooFormulaVar(slopeName.c_str(), "@0/@1", RooArgList(*sigmarv[ij(i, j)], *myrms[j]));
257 offs[ij(i, j)] = new RooFormulaVar(offsetName.c_str(), "@0-(@1*@2)",
258 RooArgList(*meanrv[ij(i, j)], *mypos[j], *slope[ij(i, j)]));
259 ownedComps.add(RooArgSet(*slope[ij(i, j)], *offs[ij(i, j)]));
260 // linear transformations, so pdf can be renormalized
261 auto& var = static_cast<RooRealVar&>(*_varList[j]);
262 std::string transVarName = Form("%s_transVar_%d_%d", GetName(), i, j);
263 // transVar[ij(i,j)] = new
264 // RooFormulaVar(transVarName.c_str(),transVarName.c_str(),"@0*@1+@2",RooArgList(*var,*slope[ij(i,j)],*offs[ij(i,j)]));
265
266 transVar[ij(i, j)] =
267 new RooLinearVar(transVarName.c_str(), transVarName.c_str(), var, *slope[ij(i, j)], *offs[ij(i, j)]);
268
269 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
270 // This will prevent the likelihood optimizers from erroneously declaring terms constant
271 transVar[ij(i, j)]->addServer((RooAbsArg &)m.arg());
272
273 ownedComps.add(*transVar[ij(i, j)]);
274 cust.replaceArg(var, *transVar[ij(i, j)]);
275 }
276 transPdf[i] = static_cast<RooAbsPdf *>(cust.build());
277 transPdfList.add(*transPdf[i]);
278 ownedComps.add(*transPdf[i]);
279 }
280 // sum pdf
281 theSumFunc = new RooRealSumFunc(sumfuncName.c_str(), sumfuncName.c_str(), transPdfList, coefList);
282 } else {
283 theSumFunc = new RooRealSumFunc(sumfuncName.c_str(), sumfuncName.c_str(), _pdfList, coefList);
284 }
285
286 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
287 // This will prevent the likelihood optimizers from erroneously declaring terms constant
288 theSumFunc->addServer((RooAbsArg &)m.arg());
289 theSumFunc->addOwnedComponents(ownedComps);
290
291 // change tracker for fraction parameters
292 std::string trackerName = Form("%s_frac_tracker", GetName());
293 RooChangeTracker *tracker = new RooChangeTracker(trackerName.c_str(), trackerName.c_str(), m.arg(), true);
294
295 // Store it in the cache
296 cache = new CacheElem(*theSumFunc, *tracker, fracl);
297 _cacheMgr.setObj(nullptr, nullptr, cache, nullptr);
298
299 return cache;
300}
301
302//_____________________________________________________________________________
304{
305 return RooArgList(*_sumFunc, *_tracker);
306}
307
308//_____________________________________________________________________________
310{
311 delete _sumFunc;
312 delete _tracker;
313}
314
315//_____________________________________________________________________________
316double RooMomentMorphFunc::getVal(const RooArgSet *set) const
317{
318 // Special version of getVal() overrides RooAbsReal::getVal() to save value of current normalization set
319 _curNormSet = set ? const_cast<RooArgSet *>(set) : const_cast<RooArgSet *>(static_cast<RooArgSet const*>(&_varList));
320 return RooAbsReal::getVal(set);
321}
322
323//_____________________________________________________________________________
325{
326 CacheElem *cache = getCache(nset ? nset : _curNormSet);
327
328 if (cache->_tracker->hasChanged(true)) {
329 cache->calculateFractions(*this, false); // verbose turned off
330 }
331
332 return cache->_sumFunc;
333}
334
335//_____________________________________________________________________________
337{
338 CacheElem *cache = getCache(nset ? nset : _curNormSet);
339
340 if (cache->_tracker->hasChanged(true)) {
341 cache->calculateFractions(*this, false); // verbose turned off
342 }
343
344 return cache->_sumFunc;
345}
346
347//_____________________________________________________________________________
349{
351
352 if (cache->_tracker->hasChanged(true)) {
353 cache->calculateFractions(*this, false); // verbose turned off
354 }
355
356 double ret = cache->_sumFunc->getVal(_pdfList.nset());
357 return ret;
358}
359
360//_____________________________________________________________________________
362{
363 return static_cast<RooRealVar *>(_frac.at(i));
364}
365
366//_____________________________________________________________________________
368{
369 return static_cast<RooRealVar *>(_frac.at(i));
370}
371
372//_____________________________________________________________________________
374{
375 Int_t nPdf = self._pdfList.size();
376
377 double dm = self.m - (*self._mref)[0];
378
379 // fully non-linear
380 double sumposfrac = 0.;
381 for (Int_t i = 0; i < nPdf; ++i) {
382 double ffrac = 0.;
383 for (Int_t j = 0; j < nPdf; ++j) {
384 ffrac += (*self._M)(j, i) * (j == 0 ? 1. : TMath::Power(dm, (double)j));
385 }
386 if (ffrac >= 0)
387 sumposfrac += ffrac;
388 // fractions for pdf
389 const_cast<RooRealVar *>(frac(i))->setVal(ffrac);
390 // fractions for rms and mean
391 const_cast<RooRealVar *>(frac(nPdf + i))->setVal(ffrac);
392 if (verbose) {
393 cout << ffrac << endl;
394 }
395 }
396
397 // various mode settings
398 int imin = self.idxmin(self.m);
399 int imax = self.idxmax(self.m);
400 double mfrac = (self.m - (*self._mref)[imin]) / ((*self._mref)[imax] - (*self._mref)[imin]);
401 switch (self._setting) {
402 case NonLinear:
403 // default already set above
404 break;
405
406 case SineLinear:
407 mfrac =
408 TMath::Sin(TMath::PiOver2() * mfrac); // this gives a continuous differentiable transition between grid points.
409
410 // now fall through to Linear case
411
412 case Linear:
413 for (Int_t i = 0; i < 2 * nPdf; ++i) const_cast<RooRealVar *>(frac(i))->setVal(0.);
414 if (imax > imin) { // m in between mmin and mmax
415 const_cast<RooRealVar *>(frac(imin))->setVal(1. - mfrac);
416 const_cast<RooRealVar *>(frac(nPdf + imin))->setVal(1. - mfrac);
417 const_cast<RooRealVar *>(frac(imax))->setVal(mfrac);
418 const_cast<RooRealVar *>(frac(nPdf + imax))->setVal(mfrac);
419 } else if (imax == imin) { // m outside mmin and mmax
420 const_cast<RooRealVar *>(frac(imin))->setVal(1.);
421 const_cast<RooRealVar *>(frac(nPdf + imin))->setVal(1.);
422 }
423 break;
425 for (Int_t i = 0; i < nPdf; ++i) const_cast<RooRealVar *>(frac(i))->setVal(0.);
426 if (imax > imin) { // m in between mmin and mmax
427 const_cast<RooRealVar *>(frac(imin))->setVal(1. - mfrac);
428 const_cast<RooRealVar *>(frac(imax))->setVal(mfrac);
429 } else if (imax == imin) { // m outside mmin and mmax
430 const_cast<RooRealVar *>(frac(imin))->setVal(1.);
431 }
432 break;
434 for (Int_t i = 0; i < nPdf; ++i) {
435 if (frac(i)->getVal() < 0)
436 const_cast<RooRealVar *>(frac(i))->setVal(0.);
437 const_cast<RooRealVar *>(frac(i))->setVal(frac(i)->getVal() / sumposfrac);
438 }
439 break;
440 }
441}
442
443//_____________________________________________________________________________
444int RooMomentMorphFunc::idxmin(const double &mval) const
445{
446 int imin(0);
447 Int_t nPdf = _pdfList.size();
448 double mmin = -DBL_MAX;
449 for (Int_t i = 0; i < nPdf; ++i) {
450 if ((*_mref)[i] > mmin && (*_mref)[i] <= mval) {
451 mmin = (*_mref)[i];
452 imin = i;
453 }
454 }
455 return imin;
456}
457
458//_____________________________________________________________________________
459int RooMomentMorphFunc::idxmax(const double &mval) const
460{
461 int imax(0);
462 Int_t nPdf = _pdfList.size();
463 double mmax = DBL_MAX;
464 for (Int_t i = 0; i < nPdf; ++i) {
465 if ((*_mref)[i] < mmax && (*_mref)[i] >= mval) {
466 mmax = (*_mref)[i];
467 imax = i;
468 }
469 }
470 return imax;
471}
472
473//_____________________________________________________________________________
474std::list<double> *RooMomentMorphFunc::plotSamplingHint(RooAbsRealLValue &obs, double xlo, double xhi) const
475{
476 return sumFunc(nullptr)->plotSamplingHint(obs, xlo, xhi);
477}
478
479//_____________________________________________________________________________
480std::list<double> *RooMomentMorphFunc::binBoundaries(RooAbsRealLValue &obs, double xlo, double xhi) const
481{
482 return sumFunc(nullptr)->binBoundaries(obs, xlo, xhi);
483}
484
485//_____________________________________________________________________________
487{
488 return sumFunc(nullptr)->isBinnedDistribution(obs);
489}
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:377
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:77
void setLocalNoDirtyInhibit(bool flag) const
Definition RooAbsArg.h:700
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
friend class RooRealSumFunc
Definition RooAbsReal.h:410
virtual bool isBinnedDistribution(const RooArgSet &) const
Tests if the distribution is binned. Unless overridden by derived classes, this always returns false.
Definition RooAbsReal.h:353
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:368
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:55
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
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...
virtual double getVal(const RooArgSet *set=nullptr) const
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
const T & arg() const
Return reference to object held in proxy.
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:75
constexpr Double_t PiOver2()
Definition TMath.h:51
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:721
Double_t Sin(Double_t)
Returns the sine of an angle of x radians.
Definition TMath.h:588
TMarker m
Definition textangle.C:8