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
26
27#include "TMath.h"
28#include "TVector.h"
29#include "TMap.h"
30
31#include <map>
32#include <algorithm>
33
34using namespace std;
35
37
38//_____________________________________________________________________________
40 : _cacheMgr(this, 10, true, true), _setting(RooMomentMorphND::Linear), _useHorizMorph(true)
41{
43}
44
45//_____________________________________________________________________________
46RooMomentMorphND::RooMomentMorphND(const char *name, const char *title, const RooArgList &parList,
47 const RooArgList &obsList, const Grid2 &referenceGrid, const Setting &setting)
48 : RooMomentMorphND::Base_t(name, title), _cacheMgr(this, 10, true, true),
49 _parList("parList", "List of morph parameters", this), _obsList("obsList", "List of observables", this),
50 _referenceGrid(referenceGrid), _pdfList("pdfList", "List of pdfs", this), _setting(setting), _useHorizMorph(true)
51{
52 // morph parameters
53 initializeParameters(parList);
54
55 // observables
56 initializeObservables(obsList);
57
59
60 // general initialization
61 initialize();
62
64}
65
66//_____________________________________________________________________________
67RooMomentMorphND::RooMomentMorphND(const char *name, const char *title, RooAbsReal &_m, const RooArgList &varList,
68 const RooArgList &pdfList, const TVectorD &mrefpoints, Setting setting)
69 : RooMomentMorphND::Base_t(name, title), _cacheMgr(this, 10, true, true),
70 _parList("parList", "List of morph parameters", this), _obsList("obsList", "List of observables", this),
71 _pdfList("pdfList", "List of pdfs", this), _setting(setting), _useHorizMorph(true)
72{
73 // make reference grid
74 RooBinning grid(mrefpoints.GetNrows() - 1, mrefpoints.GetMatrixArray());
76
77 for (int i = 0; i < mrefpoints.GetNrows(); ++i) {
78 for (int j = 0; j < grid.numBoundaries(); ++j) {
79 if (mrefpoints[i] == grid.array()[j]) {
80 _referenceGrid.addPdf(*(Base_t *)pdfList.at(i), j);
81 break;
82 }
83 }
84 }
85
87
88 // morph parameters
89 RooArgList parList;
90 parList.add(_m);
91 initializeParameters(parList);
92
93 // observables
94 initializeObservables(varList);
95
96 // general initialization
97 initialize();
98
100}
101
102//_____________________________________________________________________________
103RooMomentMorphND::RooMomentMorphND(const char *name, const char *title, RooAbsReal &_m, const RooArgList &varList,
104 const RooArgList &pdfList, const RooArgList &mrefList, Setting setting)
105 : RooMomentMorphND::Base_t(name, title), _cacheMgr(this, 10, true, true),
106 _parList("parList", "List of morph parameters", this), _obsList("obsList", "List of observables", this),
107 _pdfList("pdfList", "List of pdfs", this), _setting(setting), _useHorizMorph(true)
108{
109 // make reference grid
110 TVectorD mrefpoints(mrefList.getSize());
111 Int_t i = 0;
112 for (auto *mref : mrefList) {
113 if (!dynamic_cast<RooAbsReal *>(mref)) {
114 coutE(InputArguments) << "RooMomentMorphND::ctor(" << GetName() << ") ERROR: mref " << mref->GetName()
115 << " is not of type RooAbsReal" << endl;
116 throw string("RooMomentMorphND::ctor() ERROR mref is not of type RooAbsReal");
117 }
118 if (!dynamic_cast<RooConstVar *>(mref)) {
119 coutW(InputArguments) << "RooMomentMorphND::ctor(" << GetName() << ") WARNING mref point " << i
120 << " is not a constant, taking a snapshot of its value" << endl;
121 }
122 mrefpoints[i] = static_cast<RooAbsReal *>(mref)->getVal();
123 i++;
124 }
125
126 RooBinning grid(mrefpoints.GetNrows() - 1, mrefpoints.GetMatrixArray());
128
129 for (i = 0; i < mrefpoints.GetNrows(); ++i) {
130 for (int j = 0; j < grid.numBoundaries(); ++j) {
131 if (mrefpoints[i] == grid.array()[j]) {
132 _referenceGrid.addPdf(static_cast<Base_t &>(pdfList[i]), j);
133 break;
134 }
135 }
136 }
137
139
140 // morph parameters
141 RooArgList parList;
142 parList.add(_m);
143 initializeParameters(parList);
144
145 // observables
146 initializeObservables(varList);
147
148 // general initialization
149 initialize();
150
152}
153
154//_____________________________________________________________________________
156 : RooMomentMorphND::Base_t(other, name), _cacheMgr(other._cacheMgr, this), _parList("parList", this, other._parList),
157 _obsList("obsList", this, other._obsList), _referenceGrid(other._referenceGrid),
158 _pdfList("pdfList", this, other._pdfList), _setting(other._setting), _useHorizMorph(other._useHorizMorph)
159{
160 // general initialization
161 initialize();
162
164}
165
166//_____________________________________________________________________________
168{
170}
171
172//_____________________________________________________________________________
174{
175 for (auto *par : parList) {
176 if (!dynamic_cast<RooAbsReal *>(par)) {
177 coutE(InputArguments) << "RooMomentMorphND::ctor(" << GetName() << ") ERROR: parameter " << par->GetName()
178 << " is not of type RooAbsReal" << endl;
179 throw string("RooMomentMorphND::initializeParameters() ERROR parameter is not of type RooAbsReal");
180 }
181 _parList.add(*par);
182 }
183}
184
185//_____________________________________________________________________________
187{
188 for (auto *var : obsList) {
189 if (!dynamic_cast<RooAbsReal *>(var)) {
190 coutE(InputArguments) << "RooMomentMorphND::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
191 << " is not of type RooAbsReal" << endl;
192 throw string("RooMomentMorphND::initializeObservables() ERROR variable is not of type RooAbsReal");
193 }
194 _obsList.add(*var);
195 }
196}
197
198//_____________________________________________________________________________
200{
201 for (vector<RooAbsBinning *>::iterator itr = _referenceGrid._grid.begin(); itr != _referenceGrid._grid.end();
202 ++itr) {
203 _referenceGrid._nnuis.push_back((*itr)->numBins() + 1);
204 }
205
206 int nPar = _parList.getSize();
207 int nDim = _referenceGrid._grid.size();
208 int nPdf = _referenceGrid._pdfList.getSize();
209 int nRef = _referenceGrid._nref.size();
210 int depth = TMath::Power(2, nPar);
211
212 if (nPar != nDim) {
213 coutE(InputArguments) << "RooMomentMorphND::initialize(" << GetName() << ") ERROR: nPar != nDim"
214 << ": " << nPar << " !=" << nDim << endl;
215 assert(0);
216 }
217
218 if (nPdf != nRef) {
219 coutE(InputArguments) << "RooMomentMorphND::initialize(" << GetName() << ") ERROR: nPdf != nRef"
220 << ": " << nPdf << " !=" << nRef << endl;
221 assert(0);
222 }
223
224 // Transformation matrix for NonLinear settings
225 _M = std::make_unique<TMatrixD>(nPdf, nPdf);
226 _MSqr = std::make_unique<TMatrixD>(depth, depth);
228 TMatrixD M(nPdf, nPdf);
229
230 vector<vector<double>> dm(nPdf);
231 for (int k = 0; k < nPdf; ++k) {
232 vector<double> dm2;
233 for (int idim = 0; idim < nPar; idim++) {
234 double delta = _referenceGrid._nref[k][idim] - _referenceGrid._nref[0][idim];
235 dm2.push_back(delta);
236 }
237 dm[k] = dm2;
238 }
239
240 vector<vector<int>> powers;
241 for (int idim = 0; idim < nPar; idim++) {
242 vector<int> xtmp;
243 for (int ix = 0; ix < _referenceGrid._nnuis[idim]; ix++) {
244 xtmp.push_back(ix);
245 }
246 powers.push_back(xtmp);
247 }
248
249 vector<vector<int>> output;
251 int nCombs = output.size();
252
253 for (int k = 0; k < nPdf; ++k) {
254 int nperm = 0;
255 for (int i = 0; i < nCombs; i++) {
256 double tmpDm = 1.0;
257 for (int ix = 0; ix < nPar; ix++) {
258 double delta = dm[k][ix];
259 tmpDm *= TMath::Power(delta, static_cast<double>(output[i][ix]));
260 }
261 M(k, nperm) = tmpDm;
262 nperm++;
263 }
264 }
265
266 // M.Print();
267 (*_M) = M.Invert();
268 }
269
270 // Resize transformation vectors
271 _squareVec.resize(TMath::Power(2, nPar));
272 _squareIdx.resize(TMath::Power(2, nPar));
273}
274
275//_____________________________________________________________________________
277 : _pdfList(other._pdfList), _pdfMap(other._pdfMap), _nref(other._nref)
278{
279 for (unsigned int i = 0; i < other._grid.size(); i++) {
280 _grid.push_back(other._grid[i]->clone());
281 }
282}
283
284//_____________________________________________________________________________
286{
287 for (RooAbsBinning *binning : _grid) {
288 delete binning;
289 }
290}
291
292//_____________________________________________________________________________
294{
295 vector<int> thisBoundaries;
296 vector<double> thisBoundaryCoordinates;
297 thisBoundaries.push_back(bin_x);
298 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
299 _pdfList.add(pdf);
300 _pdfMap[thisBoundaries] = _pdfList.getSize() - 1;
301 _nref.push_back(thisBoundaryCoordinates);
302}
303
304//_____________________________________________________________________________
305void RooMomentMorphND::Grid2::addPdf(const RooMomentMorphND::Base_t &pdf, int bin_x, int bin_y)
306{
307 vector<int> thisBoundaries;
308 vector<double> thisBoundaryCoordinates;
309 thisBoundaries.push_back(bin_x);
310 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
311 thisBoundaries.push_back(bin_y);
312 thisBoundaryCoordinates.push_back(_grid[1]->array()[bin_y]);
313 _pdfList.add(pdf);
314 _pdfMap[thisBoundaries] = _pdfList.getSize() - 1;
315 _nref.push_back(thisBoundaryCoordinates);
316}
317
318//_____________________________________________________________________________
319void RooMomentMorphND::Grid2::addPdf(const RooMomentMorphND::Base_t &pdf, int bin_x, int bin_y, int bin_z)
320{
321 vector<int> thisBoundaries;
322 vector<double> thisBoundaryCoordinates;
323 thisBoundaries.push_back(bin_x);
324 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
325 thisBoundaries.push_back(bin_y);
326 thisBoundaryCoordinates.push_back(_grid[1]->array()[bin_y]);
327 thisBoundaries.push_back(bin_z);
328 thisBoundaryCoordinates.push_back(_grid[2]->array()[bin_z]);
329 _pdfList.add(pdf);
330 _pdfMap[thisBoundaries] = _pdfList.getSize() - 1;
331 _nref.push_back(thisBoundaryCoordinates);
332}
333
334//_____________________________________________________________________________
336{
337 vector<double> thisBoundaryCoordinates;
338 int nBins = bins.size();
339 for (int i = 0; i < nBins; i++) {
340 thisBoundaryCoordinates.push_back(_grid[i]->array()[bins[i]]);
341 }
342 _pdfList.add(pdf);
343 _pdfMap[bins] = _pdfList.getSize() - 1;
344 _nref.push_back(thisBoundaryCoordinates);
345}
346
347//_____________________________________________________________________________
349{
350 CacheElem *cache = static_cast<CacheElem *>(_cacheMgr.getObj(0, (RooArgSet *)0));
351 if (cache) {
352 return cache;
353 }
354
355 int nObs = _obsList.getSize();
356 int nPdf = _referenceGrid._pdfList.getSize();
357
358 RooAbsReal *null = 0;
359 vector<RooAbsReal *> meanrv(nPdf * nObs, null);
360 vector<RooAbsReal *> sigmarv(nPdf * nObs, null);
361 vector<RooAbsReal *> myrms(nObs, null);
362 vector<RooAbsReal *> mypos(nObs, null);
363 vector<RooAbsReal *> slope(nPdf * nObs, null);
364 vector<RooAbsReal *> offsets(nPdf * nObs, null);
365 vector<RooAbsReal *> transVar(nPdf * nObs, null);
366 vector<RooAbsReal *> transPdf(nPdf, null);
367
368 RooArgSet ownedComps;
369 RooArgList fracl;
370
371 // fraction parameters
372 RooArgList coefList("coefList"); // fractions multiplied with input pdfs
373 RooArgList coefList2("coefList2"); // fractions multiplied with mean position of observable contribution
374 RooArgList coefList3("coefList3"); // fractions multiplied with rms position of observable contribution
375
376 for (int i = 0; i < 3 * nPdf; ++i) {
377 string fracName = Form("frac_%d", i);
378 RooRealVar *frac = new RooRealVar(fracName.c_str(), fracName.c_str(), /*value=*/1.); // to be set later
379
380 fracl.add(*frac);
381 if (i < nPdf)
382 coefList.add(*(RooRealVar *)(fracl.at(i)));
383 else if (i < 2 * nPdf)
384 coefList2.add(*(RooRealVar *)(fracl.at(i)));
385 else
386 coefList3.add(*(RooRealVar *)(fracl.at(i)));
387 ownedComps.add(*(RooRealVar *)(fracl.at(i)));
388 }
389
390 Sum_t *theSum = nullptr;
391 string sumName = Form("%s_sum", GetName());
392
393 if (_useHorizMorph) {
394 // mean and sigma
395 RooArgList obsList(_obsList);
396 for (int i = 0; i < nPdf; ++i) {
397 for (int j = 0; j < nObs; ++j) {
398 RooAbsMoment *mom = nObs == 1 ? ((Base_t *)_pdfList.at(i))->sigma((RooRealVar &)*obsList.at(j))
399 : ((Base_t *)_pdfList.at(i))->sigma((RooRealVar &)*obsList.at(j), obsList);
400
401 mom->setLocalNoDirtyInhibit(true);
402 mom->mean()->setLocalNoDirtyInhibit(true);
403
404 sigmarv[sij(i, j)] = mom;
405 meanrv[sij(i, j)] = mom->mean();
406
407 ownedComps.add(*sigmarv[sij(i, j)]);
408 }
409 }
410
411 // slope and offset (to be set later, depend on nuisance parameters)
412 for (int j = 0; j < nObs; ++j) {
413 RooArgList meanList("meanList");
414 RooArgList rmsList("rmsList");
415 for (int i = 0; i < nPdf; ++i) {
416 meanList.add(*meanrv[sij(i, j)]);
417 rmsList.add(*sigmarv[sij(i, j)]);
418 }
419 string myrmsName = Form("%s_rms_%d", GetName(), j);
420 string myposName = Form("%s_pos_%d", GetName(), j);
421 mypos[j] = new RooAddition(myposName.c_str(), myposName.c_str(), meanList, coefList2);
422 myrms[j] = new RooAddition(myrmsName.c_str(), myrmsName.c_str(), rmsList, coefList3);
423 ownedComps.add(RooArgSet(*myrms[j], *mypos[j]));
424 }
425
426 // construction of unit pdfs
427 RooArgList transPdfList;
428
429 Int_t i = 0;
430 for (auto const *pdf : static_range_cast<Base_t *>(_pdfList)) {
431
432 string pdfName = Form("pdf_%d", i);
433 RooCustomizer cust(*pdf, pdfName.c_str());
434
435 Int_t j = 0;
436 for (auto *var : static_range_cast<RooRealVar *>(obsList)) {
437 // slope and offset formulas
438 string slopeName = Form("%s_slope_%d_%d", GetName(), i, j);
439 string offsetName = Form("%s_offset_%d_%d", GetName(), i, j);
440
441 slope[sij(i, j)] =
442 new RooFormulaVar(slopeName.c_str(), "@0/@1", RooArgList(*sigmarv[sij(i, j)], *myrms[j]));
443 offsets[sij(i, j)] = new RooFormulaVar(offsetName.c_str(), "@0-(@1*@2)",
444 RooArgList(*meanrv[sij(i, j)], *mypos[j], *slope[sij(i, j)]));
445 ownedComps.add(RooArgSet(*slope[sij(i, j)], *offsets[sij(i, j)]));
446
447 // linear transformations, so pdf can be renormalized easily
448 string transVarName = Form("%s_transVar_%d_%d", GetName(), i, j);
449 transVar[sij(i, j)] = new RooLinearVar(transVarName.c_str(), transVarName.c_str(), *var, *slope[sij(i, j)],
450 *offsets[sij(i, j)]);
451
452 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
453 // This will prevent the likelihood optimizers from erroneously declaring terms constant
454 transVar[sij(i, j)]->addServerList((RooAbsCollection &)_parList);
455
456 ownedComps.add(*transVar[sij(i, j)]);
457 cust.replaceArg(*var, *transVar[sij(i, j)]);
458 ++j;
459 }
460 transPdf[i] = static_cast<Base_t *>(cust.build());
461 transPdfList.add(*transPdf[i]);
462 ownedComps.add(*transPdf[i]);
463 ++i;
464 }
465
466 // sum pdf
467 theSum = new Sum_t(sumName.c_str(), sumName.c_str(), transPdfList, coefList);
468 } else {
469 theSum = new Sum_t(sumName.c_str(), sumName.c_str(), _pdfList, coefList);
470 }
471
472 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
473 // This will prevent the likelihood optimizers from erroneously declaring terms constant
475 theSum->addOwnedComponents(ownedComps);
476
477 // change tracker for fraction parameters
478 string trackerName = Form("%s_frac_tracker", GetName());
479 RooChangeTracker *tracker = new RooChangeTracker(trackerName.c_str(), trackerName.c_str(), _parList, true);
480
481 // Store it in the cache
482 cache = new CacheElem(*theSum, *tracker, fracl);
483 _cacheMgr.setObj(0, 0, cache, 0);
484
485 cache->calculateFractions(*this, false);
486 return cache;
487}
488
489//_____________________________________________________________________________
491{
492 return RooArgList(*_sum, *_tracker);
493}
494
495//_____________________________________________________________________________
497{
498 delete _sum;
499 delete _tracker;
500}
501
502//_____________________________________________________________________________
503double RooMomentMorphND::getVal(const RooArgSet *set) const
504{
505 // Special version of getVal() overrides Base_t::getVal() to save value of current normalization set
506 _curNormSet = set ? (RooArgSet *)set : (RooArgSet *)&_obsList;
507 return Base_t::getVal(set);
508}
509
510//_____________________________________________________________________________
512{
513 CacheElem *cache = getCache(nset ? nset : _curNormSet);
514
515 if (cache->_tracker->hasChanged(true)) {
516 cache->calculateFractions(*this, false); // verbose turned off
517 }
518 return cache->_sum;
519}
520
521//_____________________________________________________________________________
523{
525
526 if (cache->_tracker->hasChanged(true)) {
527 cache->calculateFractions(*this, false); // verbose turned off
528 }
529
530 double ret = cache->_sum->getVal(_obsList.nset());
531
532 return ret;
533}
534
535//_____________________________________________________________________________
537{
538 return (RooRealVar *)(_frac.at(i));
539}
540
541//_____________________________________________________________________________
543{
544 return (RooRealVar *)(_frac.at(i));
545}
546
547//_____________________________________________________________________________
549{
550 int nPdf = self._pdfList.getSize();
551 int nPar = self._parList.getSize();
552
553 double fracLinear(1.);
554 double fracNonLinear(1.);
555
557 // Calculate the delta vector
558 vector<double> dm2;
559 for (int idim = 0; idim < nPar; idim++) {
560 double delta = ((RooRealVar *)self._parList.at(idim))->getVal() - self._referenceGrid._nref[0][idim];
561 dm2.push_back(delta);
562 }
563
564 vector<vector<int>> powers;
565 for (int idim = 0; idim < nPar; idim++) {
566 vector<int> xtmp;
567 for (int ix = 0; ix < self._referenceGrid._nnuis[idim]; ix++) {
568 xtmp.push_back(ix);
569 }
570 powers.push_back(xtmp);
571 }
572
573 vector<vector<int>> output;
575 int nCombs = output.size();
576
577 vector<double> deltavec(nPdf, 1.0);
578
579 int nperm = 0;
580 for (int i = 0; i < nCombs; i++) {
581 double tmpDm = 1.0;
582 for (int ix = 0; ix < nPar; ix++) {
583 double delta = dm2[ix];
584 tmpDm *= TMath::Power(delta, static_cast<double>(output[i][ix]));
585 }
586 deltavec[nperm] = tmpDm;
587 nperm++;
588 }
589
590 double sumposfrac = 0.0;
591 for (int i = 0; i < nPdf; ++i) {
592 double ffrac = 0.0;
593
594 for (int j = 0; j < nPdf; ++j) {
595 ffrac += (*self._M)(j, i) * deltavec[j] * fracNonLinear;
596 }
597
598 if (ffrac >= 0) {
599 sumposfrac += ffrac;
600 }
601
602 // fractions for pdf
603 if (self._setting != NonLinearLinFractions) {
604 ((RooRealVar *)frac(i))->setVal(ffrac);
605 }
606
607 // fractions for rms and mean
608 ((RooRealVar *)frac(nPdf + i))->setVal(ffrac); // need to add up
609 ((RooRealVar *)frac(2 * nPdf + i))->setVal(ffrac); // need to add up
610
611 if (verbose) {
612 cout << "NonLinear fraction " << ffrac << endl;
613 frac(i)->Print();
614 frac(nPdf + i)->Print();
615 frac(2 * nPdf + i)->Print();
616 }
617 }
618
619 if (self._setting == NonLinearPosFractions) {
620 for (int i = 0; i < nPdf; ++i) {
621 if (((RooRealVar *)frac(i))->getVal() < 0)
622 ((RooRealVar *)frac(i))->setVal(0.);
623 ((RooRealVar *)frac(i))->setVal(((RooRealVar *)frac(i))->getVal() / sumposfrac);
624 }
625 }
626 }
627
628 if (self._setting == Linear || self._setting == NonLinearLinFractions) {
629 // zero all fractions
630 // for (int i = 0; i < 3*nPdf; ++i) {
631 for (int i = 0; i < nPdf; ++i) {
632 double initval = 0;
633 ((RooRealVar *)frac(i))->setVal(initval);
634 ((RooRealVar *)frac(nPdf + i))->setVal(initval);
635 ((RooRealVar *)frac(2 * nPdf + i))->setVal(initval);
636 }
637
638 std::vector<double> mtmp;
639
640 // loop over parList
641 for (auto *m : static_range_cast<RooRealVar *>(self._parList)) {
642 mtmp.push_back(m->getVal());
643 }
644
645 self.findShape(mtmp); // this sets _squareVec and _squareIdx quantities
646
647 int depth = TMath::Power(2, nPar);
648 vector<double> deltavec(depth, 1.0);
649
650 int nperm = 0;
651
652 vector<int> xtmp;
653 for (int ix = 0; ix < nPar; ix++) {
654 xtmp.push_back(ix);
655 }
656
657 for (int iperm = 1; iperm <= nPar; ++iperm) {
658 do {
659 double dtmp = mtmp[xtmp[0]] - self._squareVec[0][xtmp[0]];
660 for (int itmp = 1; itmp < iperm; ++itmp) {
661 dtmp *= mtmp[xtmp[itmp]] - self._squareVec[0][xtmp[itmp]];
662 }
663 deltavec[nperm + 1] = dtmp;
664 nperm++;
665 } while (RooFit::Detail::nextCombination(xtmp.begin(), xtmp.begin() + iperm, xtmp.end()));
666 }
667
668 double origFrac1(0.), origFrac2(0.);
669 for (int i = 0; i < depth; ++i) {
670 double ffrac = 0.;
671 for (int j = 0; j < depth; ++j) {
672 ffrac += (*self._MSqr)(j, i) * deltavec[j] * fracLinear;
673 }
674
675 // set fractions for pdf
676 origFrac1 = ((RooRealVar *)frac(self._squareIdx[i]))->getVal(); // already set in case of smoothlinear
677 ((RooRealVar *)frac(self._squareIdx[i]))->setVal(origFrac1 + ffrac); // need to add up
678
679 // set fractions for rms and mean
680 if (self._setting != NonLinearLinFractions) {
681 origFrac2 =
682 ((RooRealVar *)frac(nPdf + self._squareIdx[i]))->getVal(); // already set in case of smoothlinear
683 ((RooRealVar *)frac(nPdf + self._squareIdx[i]))->setVal(origFrac2 + ffrac); // need to add up
684 ((RooRealVar *)frac(2 * nPdf + self._squareIdx[i]))->setVal(origFrac2 + ffrac); // need to add up
685 }
686
687 if (verbose) {
688 cout << "Linear fraction " << ffrac << endl;
689 frac(self._squareIdx[i])->Print();
690 frac(nPdf + self._squareIdx[i])->Print();
691 frac(2 * nPdf + self._squareIdx[i])->Print();
692 }
693 }
694 }
695}
696
697//_____________________________________________________________________________
698void RooMomentMorphND::findShape(const vector<double> &x) const
699{
700 int nPar = _parList.getSize();
701 int nRef = _referenceGrid._nref.size();
702
703 // Find hypercube enclosing the location to morph to
704 // bool isEnclosed = true;
705 // for (int i = 0; i < nPar; i++) {
706 // if (x[i] < _referenceGrid._grid[i]->lowBound())
707 // isEnclosed = false;
708 // if (x[i] > _referenceGrid._grid[i]->highBound())
709 // isEnclosed = false;
710 // }
711
712 // cout << "isEnclosed = " << isEnclosed << endl;
713
714 int depth = TMath::Power(2, nPar);
715
716 vector<vector<double>> boundaries(nPar);
717 for (int idim = 0; idim < nPar; idim++) {
718 int bin = _referenceGrid._grid[idim]->binNumber(x[idim]);
719 double lo = _referenceGrid._grid[idim]->binLow(bin);
720 double hi = _referenceGrid._grid[idim]->binHigh(bin);
721 boundaries[idim].push_back(lo);
722 boundaries[idim].push_back(hi);
723 }
724
725 vector<vector<double>> output;
728
729 for (int isq = 0; isq < depth; isq++) {
730 for (int iref = 0; iref < nRef; iref++) {
731 if (_squareVec[isq] == _referenceGrid._nref[iref]) {
732 _squareIdx[isq] = iref;
733 break;
734 }
735 }
736 }
737
738 // cout << endl;
739
740 // for (int isq = 0; isq < _squareVec.size(); isq++) {
741 // cout << _squareIdx[isq];
742 // cout << " (";
743 // for (int isqq = 0; isqq < _squareVec[isq].size(); isqq++) {
744 // cout << _squareVec[isq][isqq] << ((isqq<_squareVec[isq].size()-1)?",":"");
745 // }
746 // cout << ") ";
747 // }
748
749 // construct transformation matrix for linear extrapolation
750 TMatrixD M(depth, depth);
751
752 vector<int> xtmp;
753 for (int ix = 0; ix < nPar; ix++) {
754 xtmp.push_back(ix);
755 }
756
757 for (int k = 0; k < depth; ++k) {
758 M(k, 0) = 1.0;
759
760 int nperm = 0;
761 vector<double> squareBase = _squareVec[0];
762
763 for (int iperm = 1; iperm <= nPar; ++iperm) {
764 do {
765 double dtmp = _squareVec[k][xtmp[0]] - squareBase[xtmp[0]];
766 for (int itmp = 1; itmp < iperm; ++itmp) {
767 dtmp *= _squareVec[k][xtmp[itmp]] - squareBase[xtmp[itmp]];
768 }
769 M(k, nperm + 1) = dtmp;
770 nperm++;
771 } while (RooFit::Detail::nextCombination(xtmp.begin(), xtmp.begin() + iperm, xtmp.end()));
772 }
773 }
774
775 // M.Print();
776 (*_MSqr) = M.Invert();
777}
778
779//_____________________________________________________________________________
781{
782 if (allVars.getSize() == 1) {
783 RooAbsReal *temp = const_cast<RooMomentMorphND *>(this);
784 temp->specialIntegratorConfig(true)->method1D().setLabel("RooBinIntegrator");
785 int nbins = ((RooRealVar *)allVars.first())->numBins();
786 temp->specialIntegratorConfig(true)->getConfigSection("RooBinIntegrator").setRealValue("numBins", nbins);
787 return true;
788 } else {
789 cout << "Currently BinIntegrator only knows how to deal with 1-d " << endl;
790 return false;
791 }
792 return false;
793}
#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
char name[80]
Definition: TGX11.cxx:110
#define hi
Definition: THbookFile.cxx:128
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition: TString.cxx:2468
void setLocalNoDirtyInhibit(bool flag) const
Definition: RooAbsArg.h:703
bool addOwnedComponents(const RooAbsCollection &comps)
Take ownership of the contents of 'comps'.
Definition: RooAbsArg.cxx:2211
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
RooAbsBinning is the abstract base class for RooRealVar binning definitions.
Definition: RooAbsBinning.h:25
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:60
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition: RooAbsReal.h:104
RooNumIntConfig * specialIntegratorConfig() const
Returns the specialized integrator configuration for this RooAbsReal.
RooAbsMoment * sigma(RooRealVar &obs)
Definition: RooAbsReal.h:370
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< RooAbsBinning * > _grid
std::vector< int > _nnuis
void addPdf(const RooAbsPdf &pdf, int bin_x)
std::vector< std::vector< double > > _nref
std::vector< int > _squareIdx
CacheElem * getCache(const RooArgSet *nset) const
bool setBinIntegrator(RooArgSet &allVars)
RooObjCacheManager _cacheMgr
! Transient cache manager
RooArgSet * _curNormSet
! Transient cache manager
friend class CacheElem
std::unique_ptr< TMatrixD > _M
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
std::unique_ptr< TMatrixD > _MSqr
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
void cartesianProduct(std::vector< std::vector< T > > &out, std::vector< std::vector< T > > &in)
Definition: Algorithms.h:22
bool nextCombination(const Iterator first, Iterator k, const Iterator last)
Definition: Algorithms.h:64
@ InputArguments
Definition: RooGlobalFunc.h:63
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
TMarker m
Definition: textangle.C:8
static void output()