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