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