Logo ROOT  
Reference Guide
RooMomentMorphND.cxx
Go to the documentation of this file.
1/*****************************************************************************
2 * Project: RooFit *
3 * *
4 * This code was autogenerated by RooClassFactory *
5 *****************************************************************************/
6
7#include "Riostream.h"
8
9#include "RooMomentMorphND.h"
10#include "RooAbsCategory.h"
11#include "RooRealIntegral.h"
12#include "RooRealConstant.h"
13#include "RooRealVar.h"
14#include "RooFormulaVar.h"
15#include "RooCustomizer.h"
16#include "RooAddPdf.h"
17#include "RooAddition.h"
18#include "RooAbsMoment.h"
19#include "RooMoment.h"
20#include "RooLinearVar.h"
21#include "RooChangeTracker.h"
22#include "RooNumIntConfig.h"
23#include "RooHistPdf.h"
24
25#include "TMath.h"
26#include "TVector.h"
27#include "TMap.h"
28
29#include <map>
30#include <algorithm>
31
32using namespace std;
33
35
36//_____________________________________________________________________________
38 : _cacheMgr(this, 10, true, true), _curNormSet(0), _M(0), _MSqr(0), _setting(RooMomentMorphND::Linear), _useHorizMorph(true)
39{
41}
42
43//_____________________________________________________________________________
44RooMomentMorphND::RooMomentMorphND(const char *name, const char *title, const RooArgList &parList,
45 const RooArgList &obsList, const Grid &referenceGrid, const Setting &setting)
46 : RooAbsPdf(name, title), _cacheMgr(this, 10, true, true), _parList("parList", "List of morph parameters", this),
47 _obsList("obsList", "List of observables", this), _referenceGrid(referenceGrid),
48 _pdfList("pdfList", "List of pdfs", this), _setting(setting), _useHorizMorph(true)
49{
50 // morph parameters
51 initializeParameters(parList);
52
53 // observables
54 initializeObservables(obsList);
55
57
58 // general initialization
59 initialize();
60
62}
63
64//_____________________________________________________________________________
65RooMomentMorphND::RooMomentMorphND(const char *name, const char *title, RooAbsReal &_m, const RooArgList &varList,
66 const RooArgList &pdfList, const TVectorD &mrefpoints, Setting setting)
67 : RooAbsPdf(name, title), _cacheMgr(this, 10, true, true), _parList("parList", "List of morph parameters", this),
68 _obsList("obsList", "List of observables", this), _pdfList("pdfList", "List of pdfs", this), _setting(setting),
69 _useHorizMorph(true)
70{
71 // make reference grid
72 RooBinning grid(mrefpoints.GetNrows() - 1, mrefpoints.GetMatrixArray());
74
75 for (int i = 0; i < mrefpoints.GetNrows(); ++i) {
76 for (int j = 0; j < grid.numBoundaries(); ++j) {
77 if (mrefpoints[i] == grid.array()[j]) {
78 _referenceGrid.addPdf(*(RooAbsPdf *)pdfList.at(i), j);
79 break;
80 }
81 }
82 }
83
85
86 // morph parameters
87 RooArgList parList;
88 parList.add(_m);
89 initializeParameters(parList);
90
91 // observables
92 initializeObservables(varList);
93
94 // general initialization
95 initialize();
96
98}
99
100//_____________________________________________________________________________
101RooMomentMorphND::RooMomentMorphND(const char *name, const char *title, RooAbsReal &_m, const RooArgList &varList,
102 const RooArgList &pdfList, const RooArgList &mrefList, Setting setting)
103 : RooAbsPdf(name, title), _cacheMgr(this, 10, true, true), _parList("parList", "List of morph parameters", this),
104 _obsList("obsList", "List of observables", this), _pdfList("pdfList", "List of pdfs", this), _setting(setting),
105 _useHorizMorph(true)
106{
107 // make reference grid
108 TVectorD mrefpoints(mrefList.getSize());
109 Int_t i = 0;
110 for (auto *mref : mrefList) {
111 if (!dynamic_cast<RooAbsReal*>(mref)) {
112 coutE(InputArguments) << "RooMomentMorphND::ctor(" << GetName() << ") ERROR: mref " << mref->GetName()
113 << " is not of type RooAbsReal" << endl;
114 throw string("RooMomentMorphND::ctor() ERROR mref is not of type RooAbsReal");
115 }
116 if (!dynamic_cast<RooConstVar *>(mref)) {
117 coutW(InputArguments) << "RooMomentMorphND::ctor(" << GetName() << ") WARNING mref point " << i
118 << " is not a constant, taking a snapshot of its value" << endl;
119 }
120 mrefpoints[i] = static_cast<RooAbsReal *>(mref)->getVal();
121 i++;
122 }
123
124 RooBinning grid(mrefpoints.GetNrows() - 1, mrefpoints.GetMatrixArray());
126
127 for (i = 0; i < mrefpoints.GetNrows(); ++i) {
128 for (int j = 0; j < grid.numBoundaries(); ++j) {
129 if (mrefpoints[i] == grid.array()[j]) {
130 _referenceGrid.addPdf(static_cast<RooAbsPdf &>(pdfList[i]), j);
131 break;
132 }
133 }
134 }
135
137
138 // morph parameters
139 RooArgList parList;
140 parList.add(_m);
141 initializeParameters(parList);
142
143 // observables
144 initializeObservables(varList);
145
146 // general initialization
147 initialize();
148
150}
151
152//_____________________________________________________________________________
154 : RooAbsPdf(other, name), _cacheMgr(other._cacheMgr, this), _curNormSet(0),
155 _parList("parList", this, other._parList), _obsList("obsList", this, other._obsList),
156 _referenceGrid(other._referenceGrid), _pdfList("pdfList", this, other._pdfList), _M(0), _MSqr(0),
157 _setting(other._setting), _useHorizMorph(other._useHorizMorph)
158{
159 // general initialization
160 initialize();
161
163}
164
165//_____________________________________________________________________________
167{
168 if (_M)
169 delete _M;
170 if (_MSqr)
171 delete _MSqr;
172
174}
175
176//_____________________________________________________________________________
178{
179 for (auto *par : parList) {
180 if (!dynamic_cast<RooAbsReal *>(par)) {
181 coutE(InputArguments) << "RooMomentMorphND::ctor(" << GetName() << ") ERROR: parameter " << par->GetName()
182 << " is not of type RooAbsReal" << endl;
183 throw string("RooMomentMorphND::initializeParameters() ERROR parameter is not of type RooAbsReal");
184 }
185 _parList.add(*par);
186 }
187}
188
189//_____________________________________________________________________________
191{
192 for (auto *var : obsList){
193 if (!dynamic_cast<RooAbsReal *>(var)) {
194 coutE(InputArguments) << "RooMomentMorphND::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
195 << " is not of type RooAbsReal" << endl;
196 throw string("RooMomentMorphND::initializeObservables() ERROR variable is not of type RooAbsReal");
197 }
198 _obsList.add(*var);
199 }
200}
201
202//_____________________________________________________________________________
203// from http://stackoverflow.com/a/5279601
204template <typename T>
205struct Digits {
206 typename vector<T>::const_iterator begin;
207 typename vector<T>::const_iterator end;
208 typename vector<T>::const_iterator me;
209};
210
211template <typename T>
212inline void cartesian_product(vector<vector<T>> &out, vector<vector<T>> &in)
213{
214 vector<Digits<T>> vd;
215
216 for (typename vector<vector<T>>::const_iterator it = in.begin(); it != in.end(); ++it) {
217 Digits<T> d = {(*it).begin(), (*it).end(), (*it).begin()};
218 vd.push_back(d);
219 }
220
221 while (1) {
222 vector<T> result;
223 for (typename vector<Digits<T>>::const_iterator it = vd.begin(); it != vd.end(); ++it) {
224 result.push_back(*(it->me));
225 }
226 out.push_back(result);
227
228 for (typename vector<Digits<T>>::iterator it = vd.begin();;) {
229 ++(it->me);
230 if (it->me == it->end) {
231 if (it + 1 == vd.end()) {
232 return;
233 } else {
234 it->me = it->begin;
235 ++it;
236 }
237 } else {
238 break;
239 }
240 }
241 }
242}
243
244//_____________________________________________________________________________
246{
247 for (vector<RooAbsBinning *>::iterator itr = _referenceGrid._grid.begin(); itr != _referenceGrid._grid.end();
248 ++itr) {
249 _referenceGrid._nnuis.push_back((*itr)->numBins() + 1);
250 }
251
252 int nPar = _parList.getSize();
253 int nDim = _referenceGrid._grid.size();
254 int nPdf = _referenceGrid._pdfList.getSize();
255 int nRef = _referenceGrid._nref.size();
256 int depth = TMath::Power(2, nPar);
257
258 if (nPar != nDim) {
259 coutE(InputArguments) << "RooMomentMorphND::initialize(" << GetName() << ") ERROR: nPar != nDim"
260 << ": " << nPar << " !=" << nDim << endl;
261 assert(0);
262 }
263
264 if (nPdf != nRef) {
265 coutE(InputArguments) << "RooMomentMorphND::initialize(" << GetName() << ") ERROR: nPdf != nRef"
266 << ": " << nPdf << " !=" << nRef << endl;
267 assert(0);
268 }
269
270 // Transformation matrix for NonLinear settings
271 _M = new TMatrixD(nPdf, nPdf);
272 _MSqr = new TMatrixD(depth, depth);
274 TMatrixD M(nPdf, nPdf);
275
276 vector<vector<double>> dm(nPdf);
277 for (int k = 0; k < nPdf; ++k) {
278 vector<double> dm2;
279 for (int idim = 0; idim < nPar; idim++) {
280 double delta = _referenceGrid._nref[k][idim] - _referenceGrid._nref[0][idim];
281 dm2.push_back(delta);
282 }
283 dm[k] = dm2;
284 }
285
286 vector<vector<int>> powers;
287 for (int idim = 0; idim < nPar; idim++) {
288 vector<int> xtmp;
289 for (int ix = 0; ix < _referenceGrid._nnuis[idim]; ix++) {
290 xtmp.push_back(ix);
291 }
292 powers.push_back(xtmp);
293 }
294
295 vector<vector<int>> output;
296 cartesian_product(output, powers);
297 int nCombs = output.size();
298
299 for (int k = 0; k < nPdf; ++k) {
300 int nperm = 0;
301 for (int i = 0; i < nCombs; i++) {
302 double tmpDm = 1.0;
303 for (int ix = 0; ix < nPar; ix++) {
304 double delta = dm[k][ix];
305 tmpDm *= TMath::Power(delta, static_cast<double>(output[i][ix]));
306 }
307 M(k, nperm) = tmpDm;
308 nperm++;
309 }
310 }
311
312 // M.Print();
313 (*_M) = M.Invert();
314 }
315
316 // Resize transformation vectors
317 _squareVec.resize(TMath::Power(2, nPar));
318 _squareIdx.resize(TMath::Power(2, nPar));
319}
320
321//_____________________________________________________________________________
323 : _grid(other._grid), _pdfList(other._pdfList), _pdfMap(other._pdfMap), _nref(other._nref)
324{
325}
326
327//_____________________________________________________________________________
329{
330}
331
332//_____________________________________________________________________________
333void RooMomentMorphND::Grid::addPdf(const RooAbsPdf &pdf, int bin_x)
334{
335 vector<int> thisBoundaries;
336 vector<double> thisBoundaryCoordinates;
337 thisBoundaries.push_back(bin_x);
338 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
339 _pdfList.add(pdf);
340 _pdfMap[thisBoundaries] = _pdfList.getSize() - 1;
341 _nref.push_back(thisBoundaryCoordinates);
342}
343
344//_____________________________________________________________________________
345void RooMomentMorphND::Grid::addPdf(const RooAbsPdf &pdf, int bin_x, int bin_y)
346{
347 vector<int> thisBoundaries;
348 vector<double> thisBoundaryCoordinates;
349 thisBoundaries.push_back(bin_x);
350 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
351 thisBoundaries.push_back(bin_y);
352 thisBoundaryCoordinates.push_back(_grid[1]->array()[bin_y]);
353 _pdfList.add(pdf);
354 _pdfMap[thisBoundaries] = _pdfList.getSize() - 1;
355 _nref.push_back(thisBoundaryCoordinates);
356}
357
358//_____________________________________________________________________________
359void RooMomentMorphND::Grid::addPdf(const RooAbsPdf &pdf, int bin_x, int bin_y, int bin_z)
360{
361 vector<int> thisBoundaries;
362 vector<double> thisBoundaryCoordinates;
363 thisBoundaries.push_back(bin_x);
364 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
365 thisBoundaries.push_back(bin_y);
366 thisBoundaryCoordinates.push_back(_grid[1]->array()[bin_y]);
367 thisBoundaries.push_back(bin_z);
368 thisBoundaryCoordinates.push_back(_grid[2]->array()[bin_z]);
369 _pdfList.add(pdf);
370 _pdfMap[thisBoundaries] = _pdfList.getSize() - 1;
371 _nref.push_back(thisBoundaryCoordinates);
372}
373
374//_____________________________________________________________________________
375void RooMomentMorphND::Grid::addPdf(const RooAbsPdf &pdf, vector<int> bins)
376{
377 vector<double> thisBoundaryCoordinates;
378 int nBins = bins.size();
379 for (int i = 0; i < nBins; i++) {
380 thisBoundaryCoordinates.push_back(_grid[i]->array()[bins[i]]);
381 }
382 _pdfList.add(pdf);
383 _pdfMap[bins] = _pdfList.getSize() - 1;
384 _nref.push_back(thisBoundaryCoordinates);
385}
386
387//_____________________________________________________________________________
389{
390 CacheElem *cache = static_cast<CacheElem *>(_cacheMgr.getObj(0, (RooArgSet *)0));
391 if (cache) {
392 return cache;
393 }
394
395 int nObs = _obsList.getSize();
396 int nPdf = _referenceGrid._pdfList.getSize();
397
398 RooAbsReal *null = 0;
399 vector<RooAbsReal *> meanrv(nPdf * nObs, null);
400 vector<RooAbsReal *> sigmarv(nPdf * nObs, null);
401 vector<RooAbsReal *> myrms(nObs, null);
402 vector<RooAbsReal *> mypos(nObs, null);
403 vector<RooAbsReal *> slope(nPdf * nObs, null);
404 vector<RooAbsReal *> offsetrv(nPdf * nObs, null);
405 vector<RooAbsReal *> transVar(nPdf * nObs, null);
406 vector<RooAbsReal *> transPdf(nPdf, null);
407
408 RooArgSet ownedComps;
409 RooArgList fracl;
410
411 // fraction parameters
412 RooArgList coefList("coefList"); // fractions multiplied with input pdfs
413 RooArgList coefList2("coefList2"); // fractions multiplied with mean position of observable contribution
414 RooArgList coefList3("coefList3"); // fractions multiplied with rms position of observable contribution
415
416 for (int i = 0; i < 3 * nPdf; ++i) {
417 string fracName = Form("frac_%d", i);
418 RooRealVar *frac = new RooRealVar(fracName.c_str(), fracName.c_str(), /*value=*/1.); // to be set later
419
420 fracl.add(*frac);
421 if (i < nPdf)
422 coefList.add(*(RooRealVar *)(fracl.at(i)));
423 else if (i < 2 * nPdf)
424 coefList2.add(*(RooRealVar *)(fracl.at(i)));
425 else
426 coefList3.add(*(RooRealVar *)(fracl.at(i)));
427 ownedComps.add(*(RooRealVar *)(fracl.at(i)));
428 }
429
430 RooAddPdf *theSumPdf = 0;
431 string sumpdfName = Form("%s_sumpdf", GetName());
432
433 if (_useHorizMorph) {
434 // mean and sigma
435 RooArgList obsList(_obsList);
436 for (int i = 0; i < nPdf; ++i) {
437 for (int j = 0; j < nObs; ++j) {
438 RooAbsMoment *mom = nObs == 1 ? ((RooAbsPdf *)_pdfList.at(i))->sigma((RooRealVar &)*obsList.at(j))
439 : ((RooAbsPdf *)_pdfList.at(i))->sigma((RooRealVar &)*obsList.at(j), obsList);
440
441 mom->setLocalNoDirtyInhibit(true);
442 mom->mean()->setLocalNoDirtyInhibit(true);
443
444 sigmarv[sij(i, j)] = mom;
445 meanrv[sij(i, j)] = mom->mean();
446
447 ownedComps.add(*sigmarv[sij(i, j)]);
448 }
449 }
450
451 // slope and offset (to be set later, depend on nuisance parameters)
452 for (int j = 0; j < nObs; ++j) {
453 RooArgList meanList("meanList");
454 RooArgList rmsList("rmsList");
455 for (int i = 0; i < nPdf; ++i) {
456 meanList.add(*meanrv[sij(i, j)]);
457 rmsList.add(*sigmarv[sij(i, j)]);
458 }
459 string myrmsName = Form("%s_rms_%d", GetName(), j);
460 string myposName = Form("%s_pos_%d", GetName(), j);
461 mypos[j] = new RooAddition(myposName.c_str(), myposName.c_str(), meanList, coefList2);
462 myrms[j] = new RooAddition(myrmsName.c_str(), myrmsName.c_str(), rmsList, coefList3);
463 ownedComps.add(RooArgSet(*myrms[j], *mypos[j]));
464 }
465
466 // construction of unit pdfs
467 RooArgList transPdfList;
468
469 Int_t i = 0;
470 for (auto const *pdf : static_range_cast<RooRealVar *>(_pdfList)) {
471
472 string pdfName = Form("pdf_%d", i);
473 RooCustomizer cust(*pdf, pdfName.c_str());
474
475 Int_t j = 0;
476 for (auto *var : static_range_cast<RooRealVar *>(obsList)) {
477 // slope and offset formulas
478 string slopeName = Form("%s_slope_%d_%d", GetName(), i, j);
479 string offsetName = Form("%s_offset_%d_%d", GetName(), i, j);
480
481 slope[sij(i, j)] =
482 new RooFormulaVar(slopeName.c_str(), "@0/@1", RooArgList(*sigmarv[sij(i, j)], *myrms[j]));
483 offsetrv[sij(i, j)] = new RooFormulaVar(offsetName.c_str(), "@0-(@1*@2)",
484 RooArgList(*meanrv[sij(i, j)], *mypos[j], *slope[sij(i, j)]));
485 ownedComps.add(RooArgSet(*slope[sij(i, j)], *offsetrv[sij(i, j)]));
486
487 // linear transformations, so pdf can be renormalized easily
488 string transVarName = Form("%s_transVar_%d_%d", GetName(), i, j);
489 transVar[sij(i, j)] = new RooLinearVar(transVarName.c_str(), transVarName.c_str(), *var, *slope[sij(i, j)],
490 *offsetrv[sij(i, j)]);
491
492 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
493 // This will prevent the likelihood optimizers from erroneously declaring terms constant
494 transVar[sij(i, j)]->addServerList((RooAbsCollection &)_parList);
495
496 ownedComps.add(*transVar[sij(i, j)]);
497 cust.replaceArg(*var, *transVar[sij(i, j)]);
498 ++j;
499 }
500 transPdf[i] = static_cast<RooAbsPdf *>(cust.build());
501 transPdfList.add(*transPdf[i]);
502 ownedComps.add(*transPdf[i]);
503 ++i;
504 }
505
506 // sum pdf
507 theSumPdf = new RooAddPdf(sumpdfName.c_str(), sumpdfName.c_str(), transPdfList, coefList);
508 } else {
509 theSumPdf = new RooAddPdf(sumpdfName.c_str(), sumpdfName.c_str(), _pdfList, coefList);
510 }
511
512 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
513 // This will prevent the likelihood optimizers from erroneously declaring terms constant
515 theSumPdf->addOwnedComponents(ownedComps);
516
517 // change tracker for fraction parameters
518 string trackerName = Form("%s_frac_tracker", GetName());
519 RooChangeTracker *tracker = new RooChangeTracker(trackerName.c_str(), trackerName.c_str(), _parList, true);
520
521 // Store it in the cache
522 cache = new CacheElem(*theSumPdf, *tracker, fracl);
523 _cacheMgr.setObj(0, 0, cache, 0);
524
525 cache->calculateFractions(*this, false);
526 return cache;
527}
528
529//_____________________________________________________________________________
531{
532 return RooArgList(*_sumPdf, *_tracker);
533}
534
535//_____________________________________________________________________________
537{
538 delete _sumPdf;
539 delete _tracker;
540}
541
542//_____________________________________________________________________________
543double RooMomentMorphND::getVal(const RooArgSet *set) const
544{
545 // Special version of getVal() overrides RooAbsReal::getVal() to save value of current normalization set
546 _curNormSet = set ? (RooArgSet *)set : (RooArgSet *)&_obsList;
547 return RooAbsPdf::getVal(set);
548}
549
550//_____________________________________________________________________________
552{
553 CacheElem *cache = getCache(nset ? nset : _curNormSet);
554
555 if (cache->_tracker->hasChanged(true)) {
556 cache->calculateFractions(*this, false); // verbose turned off
557 }
558 return cache->_sumPdf;
559}
560
561//_____________________________________________________________________________
563{
565
566 if (cache->_tracker->hasChanged(true)) {
567 cache->calculateFractions(*this, false); // verbose turned off
568 }
569
570 double ret = cache->_sumPdf->getVal(_obsList.nset());
571
572 return ret;
573}
574
575//_____________________________________________________________________________
577{
578 return (RooRealVar *)(_frac.at(i));
579}
580
581//_____________________________________________________________________________
583{
584 return (RooRealVar *)(_frac.at(i));
585}
586
587//_____________________________________________________________________________
588// from http://stackoverflow.com/a/5097100/8747
589template <typename Iterator>
590inline bool next_combination(const Iterator first, Iterator k, const Iterator last)
591{
592 if ((first == last) || (first == k) || (last == k)) {
593 return false;
594 }
595 Iterator itr1 = first;
596 Iterator itr2 = last;
597 ++itr1;
598 if (last == itr1) {
599 return false;
600 }
601 itr1 = last;
602 --itr1;
603 itr1 = k;
604 --itr2;
605 while (first != itr1) {
606 if (*--itr1 < *itr2) {
607 Iterator j = k;
608 while (!(*itr1 < *j)) ++j;
609 iter_swap(itr1, j);
610 ++itr1;
611 ++j;
612 itr2 = k;
613 rotate(itr1, j, last);
614 while (last != j) {
615 ++j;
616 ++itr2;
617 }
618 rotate(k, itr2, last);
619 return true;
620 }
621 }
622 rotate(first, k, last);
623 return false;
624}
625
626//_____________________________________________________________________________
628{
629 int nPdf = self._pdfList.getSize();
630 int nPar = self._parList.getSize();
631
632 double fracLinear(1.);
633 double fracNonLinear(1.);
634
636 // Calculate the delta vector
637 vector<double> dm2;
638 for (int idim = 0; idim < nPar; idim++) {
639 double delta = ((RooRealVar *)self._parList.at(idim))->getVal() - self._referenceGrid._nref[0][idim];
640 dm2.push_back(delta);
641 }
642
643 vector<vector<int>> powers;
644 for (int idim = 0; idim < nPar; idim++) {
645 vector<int> xtmp;
646 for (int ix = 0; ix < self._referenceGrid._nnuis[idim]; ix++) {
647 xtmp.push_back(ix);
648 }
649 powers.push_back(xtmp);
650 }
651
652 vector<vector<int>> output;
653 cartesian_product(output, powers);
654 int nCombs = output.size();
655
656 vector<double> deltavec(nPdf, 1.0);
657
658 int nperm = 0;
659 for (int i = 0; i < nCombs; i++) {
660 double tmpDm = 1.0;
661 for (int ix = 0; ix < nPar; ix++) {
662 double delta = dm2[ix];
663 tmpDm *= TMath::Power(delta, static_cast<double>(output[i][ix]));
664 }
665 deltavec[nperm] = tmpDm;
666 nperm++;
667 }
668
669 double sumposfrac = 0.0;
670 for (int i = 0; i < nPdf; ++i) {
671 double ffrac = 0.0;
672
673 for (int j = 0; j < nPdf; ++j) {
674 ffrac += (*self._M)(j, i) * deltavec[j] * fracNonLinear;
675 }
676
677 if (ffrac >= 0) {
678 sumposfrac += ffrac;
679 }
680
681 // fractions for pdf
682 if (self._setting != NonLinearLinFractions) {
683 ((RooRealVar *)frac(i))->setVal(ffrac);
684 }
685
686 // fractions for rms and mean
687 ((RooRealVar *)frac(nPdf + i))->setVal(ffrac); // need to add up
688 ((RooRealVar *)frac(2 * nPdf + i))->setVal(ffrac); // need to add up
689
690 if (verbose) {
691 cout << "NonLinear fraction " << ffrac << endl;
692 frac(i)->Print();
693 frac(nPdf + i)->Print();
694 frac(2 * nPdf + i)->Print();
695 }
696 }
697
698 if (self._setting == NonLinearPosFractions) {
699 for (int i = 0; i < nPdf; ++i) {
700 if (((RooRealVar *)frac(i))->getVal() < 0)
701 ((RooRealVar *)frac(i))->setVal(0.);
702 ((RooRealVar *)frac(i))->setVal(((RooRealVar *)frac(i))->getVal() / sumposfrac);
703 }
704 }
705 }
706
707 if (self._setting == Linear || self._setting == NonLinearLinFractions) {
708 // zero all fractions
709 // for (int i = 0; i < 3*nPdf; ++i) {
710 for (int i = 0; i < nPdf; ++i) {
711 double initval = 0;
712 ((RooRealVar *)frac(i))->setVal(initval);
713 ((RooRealVar *)frac(nPdf + i))->setVal(initval);
714 ((RooRealVar *)frac(2 * nPdf + i))->setVal(initval);
715 }
716
717 std::vector<double> mtmp;
718
719 // loop over parList
720 for (auto * m : static_range_cast<RooRealVar*>(self._parList)) {
721 mtmp.push_back(m->getVal());
722 }
723
724 self.findShape(mtmp); // this sets _squareVec and _squareIdx quantities
725
726 int depth = TMath::Power(2, nPar);
727 vector<double> deltavec(depth, 1.0);
728
729 int nperm = 0;
730
731 vector<int> xtmp;
732 for (int ix = 0; ix < nPar; ix++) {
733 xtmp.push_back(ix);
734 }
735
736 for (int iperm = 1; iperm <= nPar; ++iperm) {
737 do {
738 double dtmp = mtmp[xtmp[0]] - self._squareVec[0][xtmp[0]];
739 for (int itmp = 1; itmp < iperm; ++itmp) {
740 dtmp *= mtmp[xtmp[itmp]] - self._squareVec[0][xtmp[itmp]];
741 }
742 deltavec[nperm + 1] = dtmp;
743 nperm++;
744 } while (next_combination(xtmp.begin(), xtmp.begin() + iperm, xtmp.end()));
745 }
746
747 double origFrac1(0.), origFrac2(0.);
748 for (int i = 0; i < depth; ++i) {
749 double ffrac = 0.;
750 for (int j = 0; j < depth; ++j) {
751 ffrac += (*self._MSqr)(j, i) * deltavec[j] * fracLinear;
752 }
753
754 // set fractions for pdf
755 origFrac1 = ((RooRealVar *)frac(self._squareIdx[i]))->getVal(); // already set in case of smoothlinear
756 ((RooRealVar *)frac(self._squareIdx[i]))->setVal(origFrac1 + ffrac); // need to add up
757
758 // set fractions for rms and mean
759 if (self._setting != NonLinearLinFractions) {
760 origFrac2 =
761 ((RooRealVar *)frac(nPdf + self._squareIdx[i]))->getVal(); // already set in case of smoothlinear
762 ((RooRealVar *)frac(nPdf + self._squareIdx[i]))->setVal(origFrac2 + ffrac); // need to add up
763 ((RooRealVar *)frac(2 * nPdf + self._squareIdx[i]))->setVal(origFrac2 + ffrac); // need to add up
764 }
765
766 if (verbose) {
767 cout << "Linear fraction " << ffrac << endl;
768 frac(self._squareIdx[i])->Print();
769 frac(nPdf + self._squareIdx[i])->Print();
770 frac(2 * nPdf + self._squareIdx[i])->Print();
771 }
772 }
773 }
774}
775
776//_____________________________________________________________________________
777void RooMomentMorphND::findShape(const vector<double> &x) const
778{
779 int nPar = _parList.getSize();
780 int nRef = _referenceGrid._nref.size();
781
782 // Find hypercube enclosing the location to morph to
783 // bool isEnclosed = true;
784 // for (int i = 0; i < nPar; i++) {
785 // if (x[i] < _referenceGrid._grid[i]->lowBound())
786 // isEnclosed = false;
787 // if (x[i] > _referenceGrid._grid[i]->highBound())
788 // isEnclosed = false;
789 // }
790
791 // cout << "isEnclosed = " << isEnclosed << endl;
792
793 int depth = TMath::Power(2, nPar);
794
795 vector<vector<double>> boundaries(nPar);
796 for (int idim = 0; idim < nPar; idim++) {
797 int bin = _referenceGrid._grid[idim]->binNumber(x[idim]);
798 double lo = _referenceGrid._grid[idim]->binLow(bin);
799 double hi = _referenceGrid._grid[idim]->binHigh(bin);
800 boundaries[idim].push_back(lo);
801 boundaries[idim].push_back(hi);
802 }
803
804 vector<vector<double>> output;
805 cartesian_product(output, boundaries);
807
808 for (int isq = 0; isq < depth; isq++) {
809 for (int iref = 0; iref < nRef; iref++) {
810 if (_squareVec[isq] == _referenceGrid._nref[iref]) {
811 _squareIdx[isq] = iref;
812 break;
813 }
814 }
815 }
816
817 // cout << endl;
818
819 // for (int isq = 0; isq < _squareVec.size(); isq++) {
820 // cout << _squareIdx[isq];
821 // cout << " (";
822 // for (int isqq = 0; isqq < _squareVec[isq].size(); isqq++) {
823 // cout << _squareVec[isq][isqq] << ((isqq<_squareVec[isq].size()-1)?",":"");
824 // }
825 // cout << ") ";
826 // }
827
828 // construct transformation matrix for linear extrapolation
829 TMatrixD M(depth, depth);
830
831 vector<int> xtmp;
832 for (int ix = 0; ix < nPar; ix++) {
833 xtmp.push_back(ix);
834 }
835
836 for (int k = 0; k < depth; ++k) {
837 M(k, 0) = 1.0;
838
839 int nperm = 0;
840 vector<double> squareBase = _squareVec[0];
841
842 for (int iperm = 1; iperm <= nPar; ++iperm) {
843 do {
844 double dtmp = _squareVec[k][xtmp[0]] - squareBase[xtmp[0]];
845 for (int itmp = 1; itmp < iperm; ++itmp) {
846 dtmp *= _squareVec[k][xtmp[itmp]] - squareBase[xtmp[itmp]];
847 }
848 M(k, nperm + 1) = dtmp;
849 nperm++;
850 } while (next_combination(xtmp.begin(), xtmp.begin() + iperm, xtmp.end()));
851 }
852 }
853
854 // M.Print();
855 (*_MSqr) = M.Invert();
856}
857
858//_____________________________________________________________________________
860{
861 if (allVars.getSize() == 1) {
862 RooAbsReal *temp = const_cast<RooMomentMorphND *>(this);
863 temp->specialIntegratorConfig(true)->method1D().setLabel("RooBinIntegrator");
864 int nbins = ((RooRealVar *)allVars.first())->numBins();
865 temp->specialIntegratorConfig(true)->getConfigSection("RooBinIntegrator").setRealValue("numBins", nbins);
866 return true;
867 } else {
868 cout << "Currently BinIntegrator only knows how to deal with 1-d " << endl;
869 return false;
870 }
871 return false;
872}
#define d(i)
Definition: RSha256.hxx:102
bool next_combination(const Iterator first, Iterator k, const Iterator last)
void cartesian_product(vector< vector< T > > &out, vector< vector< T > > &in)
#define coutW(a)
Definition: RooMsgService.h:36
#define coutE(a)
Definition: RooMsgService.h:37
#define TRACE_DESTROY
Definition: RooTrace.h:24
#define TRACE_CREATE
Definition: RooTrace.h:23
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t result
char name[80]
Definition: TGX11.cxx:110
#define hi
Definition: THbookFile.cxx:128
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:23
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2456
void setLocalNoDirtyInhibit(bool flag) const
Definition: RooAbsArg.h:698
bool addOwnedComponents(const RooAbsCollection &comps)
Take ownership of the contents of 'comps'.
Definition: RooAbsArg.cxx:2185
void addServerList(RooAbsCollection &serverList, bool valueProp=true, bool shapeProp=false)
Register a list of RooAbsArg as servers to us by calling addServer() for each arg in the list.
Definition: RooAbsArg.cxx:387
RooAbsCollection is an abstract container object that can hold multiple RooAbsArg objects.
Int_t getSize() const
Return the number of elements in the collection.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
RooAbsArg * first() const
bool setRealValue(const char *name, double newVal=0.0, bool verbose=false)
Set value of a RooAbsRealLValye stored in set with given name to newVal No error messages are printed...
RooAbsMoment represents the first, second, or third order derivative of any RooAbsReal as calculated ...
Definition: RooAbsMoment.h:27
RooAbsReal * mean()
Definition: RooAbsMoment.h:37
const RooArgSet * nset() const
Definition: RooAbsProxy.h:48
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:62
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:91
friend class RooAddPdf
Definition: RooAbsReal.h:551
RooNumIntConfig * specialIntegratorConfig() const
Returns the specialized integrator configuration for this RooAbsReal.
RooAbsMoment * sigma(RooRealVar &obs)
Definition: RooAbsReal.h:360
RooAddPdf is an efficient implementation of a sum of PDFs of the form.
Definition: RooAddPdf.h:34
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:56
Class RooBinning is an implements RooAbsBinning in terms of an array of boundary values,...
Definition: RooBinning.h:27
double * array() const override
Return array of boundary values.
Definition: RooBinning.cxx:219
Int_t numBoundaries() const override
Return the number boundaries.
Definition: RooBinning.h:38
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 setLabel(const char *label, bool printError=true) override
Set value by specifying the name of the desired state.
RooChangeTracker is a meta object that tracks value changes in a given set of RooAbsArgs by registeri...
bool hasChanged(bool clearState)
Returns true if state has changed since last call with clearState=true.
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:26
RooCustomizer is a factory class to produce clones of a prototype composite PDF object with the same ...
Definition: RooCustomizer.h:35
void replaceArg(const RooAbsArg &orig, const RooAbsArg &subst)
Replace any occurence 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...
Definition: RooFormulaVar.h:30
RooLinearVar is the most general form of a derived real-valued object that can be used by RooRealInte...
Definition: RooLinearVar.h:30
RooChangeTracker * _tracker
void calculateFractions(const RooMomentMorphND &self, bool verbose=true) const
RooArgList containedArgs(Action) override
void addBinning(const RooAbsBinning &binning)
std::vector< std::vector< double > > _nref
std::vector< RooAbsBinning * > _grid
void addPdf(const RooAbsPdf &pdf, int bin_x)
std::vector< int > _nnuis
std::vector< int > _squareIdx
CacheElem * getCache(const RooArgSet *nset) const
bool setBinIntegrator(RooArgSet &allVars)
RooObjCacheManager _cacheMgr
! Transient cache manager
RooArgSet * _curNormSet
friend class CacheElem
std::vector< std::vector< double > > _squareVec
void initializeObservables(const RooArgList &obsList)
double evaluate() const override
Evaluate this PDF / function / constant. Needs to be overridden by all derived classes.
int sij(const int &i, const int &j) const
virtual double getVal(const RooArgSet *set=nullptr) const
void findShape(const std::vector< double > &x) const
RooSetProxy _obsList
RooAbsPdf * sumPdf(const RooArgSet *nset)
void initializeParameters(const RooArgList &parList)
~RooMomentMorphND() override
RooListProxy _parList
RooListProxy _pdfList
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
RooCategory & method1D()
RooRealVar represents a variable that can be changed from the outside.
Definition: RooRealVar.h:40
TMatrixT< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant.
Definition: TMatrixT.cxx:1397
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
Int_t GetNrows() const
Definition: TVectorT.h:75
Element * GetMatrixArray()
Definition: TVectorT.h:78
Double_t x[n]
Definition: legend1.C:17
@ InputArguments
Definition: RooGlobalFunc.h:62
std::shared_ptr< std::function< double(double)> > Linear
Definition: NeuralNet.cxx:24
null_t< F > null()
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:719
Definition: first.py:1
vector< T >::const_iterator begin
vector< T >::const_iterator end
vector< T >::const_iterator me
TMarker m
Definition: textangle.C:8
static void output()