Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooMomentMorph.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 RooMomentMorph
12 \ingroup Roofit
13
14**/
15
16#include "Riostream.h"
17
18#include "RooMomentMorph.h"
19#include "RooAbsCategory.h"
20#include "RooRealConstant.h"
21#include "RooRealVar.h"
22#include "RooFormulaVar.h"
23#include "RooCustomizer.h"
24#include "RooAddPdf.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
34
35////////////////////////////////////////////////////////////////////////////////
36/// coverity[UNINIT_CTOR]
37
39 : _cacheMgr(this,10,true,true)
40{
41}
42
43////////////////////////////////////////////////////////////////////////////////
44/// CTOR
45
46RooMomentMorph::RooMomentMorph(const char *name, const char *title, RooAbsReal &_m, const RooArgList &varList,
47 const RooArgList &pdfList, const TVectorD &mrefpoints, Setting setting)
48 : RooAbsPdf(name, title),
49 _cacheMgr(this, 10, true, true),
50 m("m", "m", this, _m),
51 _varList("varList", "List of variables", this),
52 _pdfList("pdfList", "List of pdfs", this),
53 _mref(new TVectorD(mrefpoints)),
54 _setting(setting),
55 _useHorizMorph(true)
56{
57 // observables
59
60 // reference p.d.f.s
61 _pdfList.addTyped<RooAbsPdf>(pdfList);
62
63 // initialization
64 initialize();
65}
66
67////////////////////////////////////////////////////////////////////////////////
68/// CTOR
69
70RooMomentMorph::RooMomentMorph(const char *name, const char *title, RooAbsReal &_m, const RooArgList &varList,
71 const RooArgList &pdfList, const RooArgList &mrefList, Setting setting)
72 : RooAbsPdf(name, title),
73 _cacheMgr(this, 10, true, true),
74 m("m", "m", this, _m),
75 _varList("varList", "List of variables", this),
76 _pdfList("pdfList", "List of pdfs", this),
77 _mref(new TVectorD(mrefList.size())),
78 _setting(setting),
79 _useHorizMorph(true)
80{
81 // observables
83
84 // reference p.d.f.s
85 _pdfList.addTyped<RooAbsPdf>(pdfList);
86
87 // reference points in m
88
89 Int_t i = 0;
90 for (auto *mref : mrefList) {
91 if (!dynamic_cast<RooAbsReal*>(mref)) {
92 coutE(InputArguments) << "RooMomentMorph::ctor(" << GetName() << ") ERROR: mref " << mref->GetName() << " is not of type RooAbsReal" << std::endl ;
93 throw std::string("RooPolyMorh::ctor() ERROR mref is not of type RooAbsReal") ;
94 }
95 if (!dynamic_cast<RooConstVar*>(mref)) {
96 coutW(InputArguments) << "RooMomentMorph::ctor(" << GetName() << ") WARNING mref point " << i << " is not a constant, taking a snapshot of its value" << std::endl ;
97 }
98 (*_mref)[i] = static_cast<RooAbsReal*>(mref)->getVal() ;
99 i++;
100 }
101
102 // initialization
103 initialize();
104}
105
106////////////////////////////////////////////////////////////////////////////////
107
109 : RooAbsPdf(other, name),
110 _cacheMgr(other._cacheMgr, this),
111 m("m", this, other.m),
112 _varList("varList", this, other._varList),
113 _pdfList("pdfList", this, other._pdfList),
114 _mref(new TVectorD(*other._mref)),
115 _setting(other._setting),
116 _useHorizMorph(other._useHorizMorph)
117{
118
119 // initialization
120 initialize();
121}
122
123////////////////////////////////////////////////////////////////////////////////
124
126{
127 if (_mref) delete _mref;
128 if (_M) delete _M;
129}
130
131////////////////////////////////////////////////////////////////////////////////
132
134{
135 Int_t nPdf = _pdfList.size();
136
137 // other quantities needed
138 if (nPdf!=_mref->GetNrows()) {
139 coutE(InputArguments) << "RooMomentMorph::initialize(" << GetName() << ") ERROR: nPdf != nRefPoints" << std::endl ;
140 assert(0) ;
141 }
142
143 TVectorD dm{nPdf};
144 _M = new TMatrixD(nPdf,nPdf);
145
146 // transformation matrix for non-linear extrapolation, needed in evaluate()
147 TMatrixD M(nPdf,nPdf);
148 for (Int_t i=0; i<_mref->GetNrows(); ++i) {
149 dm[i] = (*_mref)[i]-(*_mref)[0];
150 M(i,0) = 1.;
151 if (i>0) 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) = std::pow(dm[i],(double)j);
156 }
157 }
158 (*_M) = M.Invert();
159}
160
161////////////////////////////////////////////////////////////////////////////////
162
163RooMomentMorph::CacheElem::CacheElem(std::unique_ptr<RooAbsPdf> && sumPdf,
164 std::unique_ptr<RooChangeTracker> && tracker,
165 const RooArgList& flist)
166 : _sumPdf(std::move(sumPdf)), _tracker(std::move(tracker)) {
167 _frac.add(flist);
168}
169
170////////////////////////////////////////////////////////////////////////////////
171
173{
174 if (auto* cache = static_cast<CacheElem*>(_cacheMgr.getObj(nullptr,static_cast<RooArgSet*>(nullptr)))) {
175 return cache ;
176 }
177 Int_t nVar = _varList.size();
178 Int_t nPdf = _pdfList.size();
179
180 RooAbsReal* null = nullptr;
181 std::vector<RooAbsReal*> meanrv(nPdf*nVar,null);
182 std::vector<RooAbsReal*> sigmarv(nPdf*nVar,null);
183 std::vector<RooAbsReal*> myrms(nVar,null);
184 std::vector<RooAbsReal*> mypos(nVar,null);
185 std::vector<RooAbsReal*> slope(nPdf*nVar,null);
186 std::vector<RooAbsReal*> offs(nPdf*nVar,null);
187 std::vector<RooAbsReal*> transVar(nPdf*nVar,null);
188 std::vector<RooAbsReal*> transPdf(nPdf,null);
189
190 RooArgSet ownedComps ;
191
192 RooArgList fracl ;
193
194 // fraction parameters
195 RooArgList coefList("coefList");
196 RooArgList coefList2("coefList2");
197 for (Int_t i=0; i<2*nPdf; ++i) {
198 std::string fracName = "frac_" + std::to_string(i);
199
200 RooRealVar* frac = new RooRealVar(fracName.c_str(),fracName.c_str(),1.) ;
201
202 fracl.add(*frac); // to be set later
203 if (i<nPdf) coefList.add(*frac);
204 else coefList2.add(*frac);
205 ownedComps.add(*frac);
206 }
207
208 std::unique_ptr<RooAddPdf> theSumPdf;
209 std::string sumpdfName = Form("%s_sumpdf",GetName());
210
211 if (_useHorizMorph) {
212 // mean and sigma
213 RooArgList varList(_varList) ;
214 for (Int_t i=0; i<nPdf; ++i) {
215 for (Int_t j=0; j<nVar; ++j) {
216
217 std::string meanName = Form("%s_mean_%d_%d",GetName(),i,j);
218 std::string sigmaName = Form("%s_sigma_%d_%d",GetName(),i,j);
219
220 RooAbsMoment* mom = nVar==1 ?
221 (static_cast<RooAbsPdf*>(_pdfList.at(i)))->sigma(static_cast<RooRealVar&>(*varList.at(j))) :
222 (static_cast<RooAbsPdf*>(_pdfList.at(i)))->sigma(static_cast<RooRealVar&>(*varList.at(j)),varList) ;
223
224 mom->setLocalNoDirtyInhibit(true) ;
225 mom->mean()->setLocalNoDirtyInhibit(true) ;
226
227 sigmarv[ij(i,j)] = mom ;
228 meanrv[ij(i,j)] = mom->mean() ;
229
230 ownedComps.add(*sigmarv[ij(i,j)]) ;
231 }
232 }
233
234 // slope and offset (to be set later, depend on m)
235 for (Int_t j=0; j<nVar; ++j) {
236 RooArgList meanList("meanList");
237 RooArgList rmsList("rmsList");
238 for (Int_t i=0; i<nPdf; ++i) {
239 meanList.add(*meanrv[ij(i,j)]);
240 rmsList.add(*sigmarv[ij(i,j)]);
241 }
242 std::string myrmsName = Form("%s_rms_%d",GetName(),j);
243 std::string myposName = Form("%s_pos_%d",GetName(),j);
244 myrms[j] = new RooAddition(myrmsName.c_str(),myrmsName.c_str(),rmsList,coefList2);
245 mypos[j] = new RooAddition(myposName.c_str(),myposName.c_str(),meanList,coefList2);
246 ownedComps.add(RooArgSet(*myrms[j],*mypos[j])) ;
247 }
248
249 // construction of unit pdfs
250 RooArgList transPdfList;
251
252 for (Int_t i=0; i<nPdf; ++i) {
253
254 auto& pdf = static_cast<RooAbsPdf&>(_pdfList[i]);
255 std::string pdfName = "pdf_" + std::to_string(i);
256 RooCustomizer cust(pdf,pdfName.c_str());
257
258 for (Int_t j=0; j<nVar; ++j) {
259 // slope and offset formulas
260 std::string slopeName = Form("%s_slope_%d_%d",GetName(),i,j);
261 std::string offsetName = Form("%s_offset_%d_%d",GetName(),i,j);
262 slope[ij(i,j)] = new RooFormulaVar(slopeName.c_str(),"@0/@1", {*sigmarv[ij(i,j)],*myrms[j]});
263 offs[ij(i,j)] = new RooFormulaVar(offsetName.c_str(),"@0-(@1*@2)", {*meanrv[ij(i,j)],*mypos[j],*slope[ij(i,j)]});
264 ownedComps.add(RooArgSet(*slope[ij(i,j)],*offs[ij(i,j)])) ;
265 // linear transformations, so pdf can be renormalized
266 auto* var = static_cast<RooRealVar*>(_varList[j]);
267 std::string transVarName = Form("%s_transVar_%d_%d",GetName(),i,j);
268 //transVar[ij(i,j)] = new RooFormulaVar(transVarName.c_str(),transVarName.c_str(),"@0*@1+@2",RooArgList(*var,*slope[ij(i,j)],*offs[ij(i,j)]));
269
270 transVar[ij(i,j)] = new RooLinearVar(transVarName.c_str(),transVarName.c_str(),*var,*slope[ij(i,j)],*offs[ij(i,j)]);
271
272 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
273 // This will prevent the likelihood optimizers from erroneously declaring terms constant
274 transVar[ij(i,j)]->addServer(const_cast<RooAbsReal&>(m.arg()));
275
276 ownedComps.add(*transVar[ij(i,j)]) ;
277 cust.replaceArg(*var,*transVar[ij(i,j)]);
278 }
279 transPdf[i] = static_cast<RooAbsPdf*>(cust.build()) ;
280 transPdfList.add(*transPdf[i]);
281 ownedComps.add(*transPdf[i]) ;
282 }
283 // sum pdf
284 theSumPdf = std::make_unique<RooAddPdf>(sumpdfName.c_str(),sumpdfName.c_str(),transPdfList,coefList);
285 }
286 else {
287 theSumPdf = std::make_unique<RooAddPdf>(sumpdfName.c_str(),sumpdfName.c_str(),_pdfList,coefList);
288 }
289 theSumPdf->fixCoefNormalization(*nset);
290
291 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
292 // This will prevent the likelihood optimizers from erroneously declaring terms constant
293 theSumPdf->addServer(const_cast<RooAbsReal&>(m.arg()));
294 theSumPdf->addOwnedComponents(ownedComps) ;
295
296 // change tracker for fraction parameters
297 std::string trackerName = Form("%s_frac_tracker",GetName()) ;
298 auto tracker = std::make_unique<RooChangeTracker>(trackerName.c_str(),trackerName.c_str(),m.arg(),true) ;
299
300 // Store it in the cache
301 auto cache = new CacheElem(std::move(theSumPdf),std::move(tracker),fracl) ;
302 _cacheMgr.setObj(nullptr,nullptr,cache,nullptr);
303
304 cache->calculateFractions(*this, false);
305 return cache ;
306}
307
308////////////////////////////////////////////////////////////////////////////////
309
311{
312 return RooArgList(*_sumPdf,*_tracker) ;
313}
314
315////////////////////////////////////////////////////////////////////////////////
316
318
319////////////////////////////////////////////////////////////////////////////////
320/// Special version of getVal() overrides RooAbsReal::getVal() to save value of current normalization set
321
322double RooMomentMorph::getValV(const RooArgSet* set) const
323{
324 _curNormSet = set ? const_cast<RooArgSet*>(set) : const_cast<RooArgSet*>(static_cast<RooArgSet const*>(&_varList));
325 return RooAbsPdf::getValV(set) ;
326}
327
328////////////////////////////////////////////////////////////////////////////////
329
331{
332 CacheElem* cache = getCache(nset ? nset : _curNormSet) ;
333
334 if (cache->_tracker->hasChanged(true)) {
335 cache->calculateFractions(*this,false); // verbose turned off
336 }
337
338 return cache->_sumPdf.get();
339}
340
341////////////////////////////////////////////////////////////////////////////////
342
344{
345 CacheElem* cache = getCache(_curNormSet) ;
346
347 if (cache->_tracker->hasChanged(true)) {
348 cache->calculateFractions(*this,false); // verbose turned off
349 }
350
351 return cache->_sumPdf->getVal(_pdfList.nset());
352}
353
354////////////////////////////////////////////////////////////////////////////////
355
357{
358 return static_cast<RooRealVar*>(_frac.at(i)) ;
359}
360
361////////////////////////////////////////////////////////////////////////////////
362
364{
365 return static_cast<RooRealVar*>(_frac.at(i)) ;
366}
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) { ffrac += (*self._M)(j,i) * (j==0?1.:std::pow(dm,(double)j)); }
381 if (ffrac>=0) sumposfrac+=ffrac;
382 // fractions for pdf
383 const_cast<RooRealVar*>(frac(i))->setVal(ffrac);
384 // fractions for rms and mean
385 const_cast<RooRealVar*>(frac(nPdf+i))->setVal(ffrac);
386 if (verbose) { std::cout << ffrac << std::endl; }
387 }
388
389 // various mode settings
390 int imin = self.idxmin(self.m);
391 int imax = self.idxmax(self.m);
392 double mfrac = (self.m-(*self._mref)[imin])/((*self._mref)[imax]-(*self._mref)[imin]);
393 switch (self._setting) {
394 case NonLinear:
395 // default already set above
396 break;
397
398 case SineLinear:
399 mfrac = std::sin( TMath::PiOver2()*mfrac ); // this gives a continuous differentiable transition between grid points.
400
401 // now fall through to Linear case
402
403 case Linear:
404 for (Int_t i=0; i<2*nPdf; ++i)
405 const_cast<RooRealVar*>(frac(i))->setVal(0.);
406 if (imax>imin) { // m in between mmin and mmax
407 const_cast<RooRealVar*>(frac(imin))->setVal(1.-mfrac);
408 const_cast<RooRealVar*>(frac(nPdf+imin))->setVal(1.-mfrac);
409 const_cast<RooRealVar*>(frac(imax))->setVal(mfrac);
410 const_cast<RooRealVar*>(frac(nPdf+imax))->setVal(mfrac);
411 } else if (imax==imin) { // m outside mmin and mmax
412 const_cast<RooRealVar*>(frac(imin))->setVal(1.);
413 const_cast<RooRealVar*>(frac(nPdf+imin))->setVal(1.);
414 }
415 break;
417 for (Int_t i=0; i<nPdf; ++i)
418 const_cast<RooRealVar*>(frac(i))->setVal(0.);
419 if (imax>imin) { // m in between mmin and mmax
420 const_cast<RooRealVar*>(frac(imin))->setVal(1.-mfrac);
421 const_cast<RooRealVar*>(frac(imax))->setVal(mfrac);
422 } else if (imax==imin) { // m outside mmin and mmax
423 const_cast<RooRealVar*>(frac(imin))->setVal(1.);
424 }
425 break;
427 for (Int_t i = 0; i < nPdf; ++i) {
428 if (frac(i)->getVal() < 0)
429 const_cast<RooRealVar *>(frac(i))->setVal(0.);
430 const_cast<RooRealVar *>(frac(i))->setVal(frac(i)->getVal() / sumposfrac);
431 }
432 break;
433 }
434
435}
436
437////////////////////////////////////////////////////////////////////////////////
438
439int RooMomentMorph::idxmin(const double& mval) const
440{
441 int imin(0);
442 Int_t nPdf = _pdfList.size();
443 double mmin=-DBL_MAX;
444 for (Int_t i=0; i<nPdf; ++i)
445 if ( (*_mref)[i]>mmin && (*_mref)[i]<=mval ) { mmin=(*_mref)[i]; imin=i; }
446 return imin;
447}
448
449
450////////////////////////////////////////////////////////////////////////////////
451
452int RooMomentMorph::idxmax(const double& mval) const
453{
454 int imax(0);
455 Int_t nPdf = _pdfList.size();
456 double mmax=DBL_MAX;
457 for (Int_t i=0; i<nPdf; ++i)
458 if ( (*_mref)[i]<mmax && (*_mref)[i]>=mval ) { mmax=(*_mref)[i]; imax=i; }
459 return imax;
460}
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
void setLocalNoDirtyInhibit(bool flag) const
Definition RooAbsArg.h:700
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
double getValV(const RooArgSet *set=nullptr) const override
Return current value, normalized by integrating over the observables in nset.
const RooArgSet * nset() const
Definition RooAbsProxy.h:52
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
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.
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...
RooArgList containedArgs(Action) override
CacheElem(std::unique_ptr< RooAbsPdf > &&sumPdf, std::unique_ptr< RooChangeTracker > &&tracker, const RooArgList &flist)
RooRealVar * frac(Int_t i)
void calculateFractions(const RooMomentMorph &self, bool verbose=true) const
std::unique_ptr< RooAbsPdf > _sumPdf
std::unique_ptr< RooChangeTracker > _tracker
RooObjCacheManager _cacheMgr
~RooMomentMorph() override
RooSetProxy _varList
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
RooMomentMorph()
coverity[UNINIT_CTOR]
RooRealProxy m
RooAbsPdf * sumPdf(const RooArgSet *nset)
Int_t ij(const Int_t &i, const Int_t &j) const
RooArgSet * _curNormSet
The cache manager.
RooListProxy _pdfList
int idxmin(const double &m) const
int idxmax(const double &m) const
double getValV(const RooArgSet *set=nullptr) const override
Special version of getVal() overrides RooAbsReal::getVal() to save value of current normalization set...
CacheElem * getCache(const RooArgSet *nset) const
TVectorD * _mref
Variable that can be changed from the outside.
Definition RooRealVar.h:37
void setVal(double value) override
Set value of variable to 'value'.
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:73
constexpr Double_t PiOver2()
Definition TMath.h:51
TMarker m
Definition textangle.C:8