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