Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
RooMomentMorphFuncND.cxx
Go to the documentation of this file.
1/*
2 * Project: RooFit
3 *
4 * Copyright (c) 2023, CERN
5 *
6 * Redistribution and use in source and binary forms,
7 * with or without modification, are permitted according to the terms
8 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)
9 */
10
11/** \class RooMomentMorphFuncND
12 \ingroup Roofit
13
14**/
15
17
18#include <RooAbsCategory.h>
19#include <RooAbsMoment.h>
20#include <RooAddPdf.h>
21#include <RooAddition.h>
22#include <RooChangeTracker.h>
23#include <RooConstVar.h>
24#include <RooCustomizer.h>
25#include <RooFormulaVar.h>
26#include <RooLinearVar.h>
27#include <RooMoment.h>
28#include <RooNumIntConfig.h>
29#include <RooRealSumFunc.h>
30#include <RooRealVar.h>
31
33
34#include <Riostream.h>
35
36#include <TMap.h>
37#include <TMath.h>
38#include <TVector.h>
39
40#include <algorithm>
41#include <map>
42
43using namespace std;
44
47
48//_____________________________________________________________________________
49RooMomentMorphFuncND::RooMomentMorphFuncND() : _cacheMgr(this, 10, true, true), _setting(RooMomentMorphFuncND::Linear), _useHorizMorph(true)
50{
52}
53
54//_____________________________________________________________________________
55RooMomentMorphFuncND::RooMomentMorphFuncND(const char *name, const char *title, const RooArgList &parList, const RooArgList &obsList,
56 const Grid2 &referenceGrid, const Setting &setting)
58 _cacheMgr(this, 10, true, true),
59 _parList("parList", "List of morph parameters", this),
60 _obsList("obsList", "List of observables", this),
61 _referenceGrid(referenceGrid),
62 _pdfList("pdfList", "List of pdfs", this),
63 _setting(setting),
64 _useHorizMorph(true)
65{
66 // morph parameters
68
69 // observables
71
73
74 // general initialization
75 initialize();
76
78}
79
80//_____________________________________________________________________________
81RooMomentMorphFuncND::RooMomentMorphFuncND(const char *name, const char *title, RooAbsReal &_m, const RooArgList &varList,
82 const RooArgList &pdfList, const TVectorD &mrefpoints, Setting setting)
84 _cacheMgr(this, 10, true, true),
85 _parList("parList", "List of morph parameters", this),
86 _obsList("obsList", "List of observables", this),
87 _pdfList("pdfList", "List of pdfs", this),
88 _setting(setting),
89 _useHorizMorph(true)
90{
91 // make reference grid
92 RooBinning grid(mrefpoints.GetNrows() - 1, mrefpoints.GetMatrixArray());
94
95 for (int i = 0; i < mrefpoints.GetNrows(); ++i) {
96 for (int j = 0; j < grid.numBoundaries(); ++j) {
97 if (mrefpoints[i] == grid.array()[j]) {
98 _referenceGrid.addPdf(*static_cast<Base_t *>(pdfList.at(i)), j);
99 break;
100 }
101 }
102 }
103
105
106 // morph parameters
107 RooArgList parList;
108 parList.add(_m);
109 _parList.addTyped<RooAbsReal>(parList);
110
111 // observables
112 _obsList.addTyped<RooAbsReal>(varList);
113
114 // general initialization
115 initialize();
116
118}
119
120//_____________________________________________________________________________
121RooMomentMorphFuncND::RooMomentMorphFuncND(const char *name, const char *title, RooAbsReal &_m, const RooArgList &varList,
122 const RooArgList &pdfList, const RooArgList &mrefList, Setting setting)
124 _cacheMgr(this, 10, true, true),
125 _parList("parList", "List of morph parameters", this),
126 _obsList("obsList", "List of observables", this),
127 _pdfList("pdfList", "List of pdfs", this),
128 _setting(setting),
129 _useHorizMorph(true)
130{
131 // make reference grid
132 TVectorD mrefpoints(mrefList.size());
133 Int_t i = 0;
134 for (auto *mref : mrefList) {
135 if (!dynamic_cast<RooAbsReal *>(mref)) {
136 coutE(InputArguments) << "RooMomentMorphFuncND::ctor(" << GetName() << ") ERROR: mref " << mref->GetName()
137 << " is not of type RooAbsReal" << endl;
138 throw string("RooMomentMorphFuncND::ctor() ERROR mref is not of type RooAbsReal");
139 }
140 if (!dynamic_cast<RooConstVar *>(mref)) {
141 coutW(InputArguments) << "RooMomentMorphFuncND::ctor(" << GetName() << ") WARNING mref point " << i
142 << " is not a constant, taking a snapshot of its value" << endl;
143 }
144 mrefpoints[i] = static_cast<RooAbsReal *>(mref)->getVal();
145 i++;
146 }
147
148 RooBinning grid(mrefpoints.GetNrows() - 1, mrefpoints.GetMatrixArray());
150 for (i = 0; i < mrefpoints.GetNrows(); ++i) {
151 for (int j = 0; j < grid.numBoundaries(); ++j) {
152 if (mrefpoints[i] == grid.array()[j]) {
153 _referenceGrid.addPdf(static_cast<Base_t &>(pdfList[i]), j);
154 break;
155 }
156 }
157 }
158
160
161 // morph parameters
162 RooArgList parList;
163 parList.add(_m);
164 _parList.addTyped<RooAbsReal>(parList);
165
166 // observables
167 _obsList.addTyped<RooAbsReal>(varList);
168
169 // general initialization
170 initialize();
171
173}
174
175//_____________________________________________________________________________
178 _cacheMgr(other._cacheMgr, this),
179 _parList("parList", this, other._parList),
180 _obsList("obsList", this, other._obsList),
181 _referenceGrid(other._referenceGrid),
182 _pdfList("pdfList", this, other._pdfList),
183 _setting(other._setting),
184 _useHorizMorph(other._useHorizMorph)
185{
186 // general initialization
187 initialize();
188
190}
191
192//_____________________________________________________________________________
194{
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.size();
207 int nDim = _referenceGrid._grid.size();
208 int nPdf = _referenceGrid._pdfList.size();
209 int nRef = _referenceGrid._nref.size();
210 int depth = TMath::Power(2, nPar);
211
212 if (nPar != nDim) {
213 coutE(InputArguments) << "RooMomentMorphFuncND::initialize(" << GetName() << ") ERROR: nPar != nDim"
214 << ": " << nPar << " !=" << nDim << endl;
215 assert(0);
216 }
217
218 if (nPdf != nRef) {
219 coutE(InputArguments) << "RooMomentMorphFuncND::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 xtmp.reserve(_referenceGrid._nnuis[idim]);
244 for (int ix = 0; ix < _referenceGrid._nnuis[idim]; ix++) {
245 xtmp.push_back(ix);
246 }
247 powers.push_back(xtmp);
248 }
249
250 vector<vector<int>> output;
252 int nCombs = output.size();
253
254 for (int k = 0; k < nPdf; ++k) {
255 int nperm = 0;
256 for (int i = 0; i < nCombs; i++) {
257 double tmpDm = 1.0;
258 for (int ix = 0; ix < nPar; ix++) {
259 double delta = dm[k][ix];
260 tmpDm *= TMath::Power(delta, static_cast<double>(output[i][ix]));
261 }
262 M(k, nperm) = tmpDm;
263 nperm++;
264 }
265 }
266
267 // M.Print();
268 (*_M) = M.Invert();
269 }
270
271 // Resize transformation vectors
272 _squareVec.resize(TMath::Power(2, nPar));
273 _squareIdx.resize(TMath::Power(2, nPar));
274}
275
276//_____________________________________________________________________________
278 : _pdfList(other._pdfList), _pdfMap(other._pdfMap), _nref(other._nref)
279{
280 for (unsigned int i = 0; i < other._grid.size(); i++) {
281 _grid.push_back(other._grid[i]->clone());
282 }
283}
284
285//_____________________________________________________________________________
287{
288 for (RooAbsBinning *binning : _grid) {
289 delete binning;
290 }
291}
292
293//_____________________________________________________________________________
295{
296 vector<int> thisBoundaries;
297 vector<double> thisBoundaryCoordinates;
298 thisBoundaries.push_back(bin_x);
299 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
300 _pdfList.add(pdf);
301 _pdfMap[thisBoundaries] = _pdfList.size() - 1;
302 _nref.push_back(thisBoundaryCoordinates);
303}
304
305//_____________________________________________________________________________
307{
308 vector<int> thisBoundaries;
309 vector<double> thisBoundaryCoordinates;
310 thisBoundaries.push_back(bin_x);
311 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
312 thisBoundaries.push_back(bin_y);
313 thisBoundaryCoordinates.push_back(_grid[1]->array()[bin_y]);
314 _pdfList.add(pdf);
315 _pdfMap[thisBoundaries] = _pdfList.size() - 1;
316 _nref.push_back(thisBoundaryCoordinates);
317}
318
319//_____________________________________________________________________________
320void RooMomentMorphFuncND::Grid2::addPdf(const RooMomentMorphFuncND::Base_t &pdf, int bin_x, int bin_y, int bin_z)
321{
322 vector<int> thisBoundaries;
323 vector<double> thisBoundaryCoordinates;
324 thisBoundaries.push_back(bin_x);
325 thisBoundaryCoordinates.push_back(_grid[0]->array()[bin_x]);
326 thisBoundaries.push_back(bin_y);
327 thisBoundaryCoordinates.push_back(_grid[1]->array()[bin_y]);
328 thisBoundaries.push_back(bin_z);
329 thisBoundaryCoordinates.push_back(_grid[2]->array()[bin_z]);
330 _pdfList.add(pdf);
331 _pdfMap[thisBoundaries] = _pdfList.size() - 1;
332 _nref.push_back(thisBoundaryCoordinates);
333}
334
335//_____________________________________________________________________________
337{
338 vector<double> thisBoundaryCoordinates;
339 int nBins = bins.size();
340 thisBoundaryCoordinates.reserve(nBins);
341 for (int i = 0; i < nBins; i++) {
342 thisBoundaryCoordinates.push_back(_grid[i]->array()[bins[i]]);
343 }
344 _pdfList.add(pdf);
345 _pdfMap[bins] = _pdfList.size() - 1;
346 _nref.push_back(thisBoundaryCoordinates);
347}
348
349//_____________________________________________________________________________
351{
352 auto cache = static_cast<CacheElem *>(_cacheMgr.getObj(nullptr, static_cast<RooArgSet const*>(nullptr)));
353 if (cache) {
354 return cache;
355 }
356
357 int nObs = _obsList.size();
358 int nPdf = _referenceGrid._pdfList.size();
359
360 RooAbsReal *null = nullptr;
361 vector<RooAbsReal *> meanrv(nPdf * nObs, null);
362 vector<RooAbsReal *> sigmarv(nPdf * nObs, null);
363 vector<RooAbsReal *> myrms(nObs, null);
364 vector<RooAbsReal *> mypos(nObs, null);
365 vector<RooAbsReal *> slope(nPdf * nObs, null);
366 vector<RooAbsReal *> offsets(nPdf * nObs, null);
367 vector<RooAbsReal *> transVar(nPdf * nObs, null);
368 vector<RooAbsReal *> transPdf(nPdf, null);
369
370 RooArgSet ownedComps;
371 RooArgList fracl;
372
373 // fraction parameters
374 RooArgList coefList("coefList"); // fractions multiplied with input pdfs
375 RooArgList coefList2("coefList2"); // fractions multiplied with mean position of observable contribution
376 RooArgList coefList3("coefList3"); // fractions multiplied with rms position of observable contribution
377
378 for (int i = 0; i < 3 * nPdf; ++i) {
379 string fracName = Form("frac_%d", i);
380 double initval = _isPdfMode ? 1.0 : 0.0;
381 RooRealVar *frac = new RooRealVar(fracName.c_str(), fracName.c_str(), initval); // to be set later
382
383 fracl.add(*frac);
384 if (i < nPdf) {
385 coefList.add(*static_cast<RooRealVar *>(fracl.at(i)));
386 } else if (i < 2 * nPdf) {
387 coefList2.add(*static_cast<RooRealVar *>(fracl.at(i)));
388 } else {
389 coefList3.add(*static_cast<RooRealVar *>(fracl.at(i)));
390 }
391 ownedComps.add(*static_cast<RooRealVar *>(fracl.at(i)));
392 }
393
394 std::unique_ptr<RooAbsReal> theSum;
395 string sumName = Form("%s_sum", GetName());
396
397 RooArgList transPdfList;
398 if (_useHorizMorph) {
399 // mean and sigma
400 RooArgList obsList(_obsList);
401 for (int i = 0; i < nPdf; ++i) {
402 for (int j = 0; j < nObs; ++j) {
403 RooAbsMoment *mom = nObs == 1 ? (static_cast<Base_t *>(_pdfList.at(i)))->sigma(static_cast<RooRealVar &>(*obsList.at(j)))
404 : (static_cast<Base_t *>(_pdfList.at(i)))->sigma(static_cast<RooRealVar &>(*obsList.at(j)), obsList);
405
406 mom->setLocalNoDirtyInhibit(true);
407 mom->mean()->setLocalNoDirtyInhibit(true);
408
409 sigmarv[sij(i, j)] = mom;
410 meanrv[sij(i, j)] = mom->mean();
411
412 ownedComps.add(*sigmarv[sij(i, j)]);
413 }
414 }
415
416 // slope and offset (to be set later, depend on nuisance parameters)
417 for (int j = 0; j < nObs; ++j) {
418 RooArgList meanList("meanList");
419 RooArgList rmsList("rmsList");
420 for (int i = 0; i < nPdf; ++i) {
421 meanList.add(*meanrv[sij(i, j)]);
422 rmsList.add(*sigmarv[sij(i, j)]);
423 }
424 string myrmsName = Form("%s_rms_%d", GetName(), j);
425 string myposName = Form("%s_pos_%d", GetName(), j);
426 mypos[j] = new RooAddition(myposName.c_str(), myposName.c_str(), meanList, coefList2);
427 myrms[j] = new RooAddition(myrmsName.c_str(), myrmsName.c_str(), rmsList, coefList3);
428 ownedComps.add(RooArgSet(*myrms[j], *mypos[j]));
429 }
430
431 // construction of unit pdfs
432
433 Int_t i = 0;
434 for (auto const *pdf : static_range_cast<Base_t *>(_pdfList)) {
435
436 string pdfName = Form("pdf_%d", i);
437 RooCustomizer cust(*pdf, pdfName.c_str());
438
439 Int_t j = 0;
440 for (auto *var : static_range_cast<RooRealVar *>(obsList)) {
441 // slope and offset formulas
442 string slopeName = Form("%s_slope_%d_%d", GetName(), i, j);
443 string offsetName = Form("%s_offset_%d_%d", GetName(), i, j);
444
445 slope[sij(i, j)] =
446 new RooFormulaVar(slopeName.c_str(), "@0/@1", RooArgList(*sigmarv[sij(i, j)], *myrms[j]));
447 offsets[sij(i, j)] = new RooFormulaVar(offsetName.c_str(), "@0-(@1*@2)",
448 RooArgList(*meanrv[sij(i, j)], *mypos[j], *slope[sij(i, j)]));
449 ownedComps.add(RooArgSet(*slope[sij(i, j)], *offsets[sij(i, j)]));
450
451 // linear transformations, so pdf can be renormalized easily
452 string transVarName = Form("%s_transVar_%d_%d", GetName(), i, j);
453 transVar[sij(i, j)] = new RooLinearVar(transVarName.c_str(), transVarName.c_str(), *var, *slope[sij(i, j)],
454 *offsets[sij(i, j)]);
455
456 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
457 // This will prevent the likelihood optimizers from erroneously declaring terms constant
458 transVar[sij(i, j)]->addServerList((RooAbsCollection &)_parList);
459
460 ownedComps.add(*transVar[sij(i, j)]);
461 cust.replaceArg(*var, *transVar[sij(i, j)]);
462 ++j;
463 }
464 transPdf[i] = static_cast<Base_t *>(cust.build());
465 transPdfList.add(*transPdf[i]);
466 ownedComps.add(*transPdf[i]);
467 ++i;
468 }
469 }
470
471 // sum pdf
472 RooArgList const &pdfList = _useHorizMorph ? transPdfList : static_cast<RooArgList const &>(_pdfList);
473 if (_isPdfMode) {
474 theSum = std::make_unique<RooAddPdf>(sumName.c_str(), sumName.c_str(), pdfList, coefList);
475 } else {
476 theSum = std::make_unique<RooRealSumFunc>(sumName.c_str(), sumName.c_str(), pdfList, coefList);
477 }
478
479 // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
480 // This will prevent the likelihood optimizers from erroneously declaring terms constant
481 theSum->addServerList((RooAbsCollection &)_parList);
482 theSum->addOwnedComponents(ownedComps);
483
484 // change tracker for fraction parameters
485 std::string trackerName = std::string(GetName()) + "_frac_tracker";
486
487 // Store it in the cache
488 cache = new CacheElem(std::move(theSum),
489 std::make_unique<RooChangeTracker>(trackerName.c_str(), trackerName.c_str(), _parList, true),
490 fracl);
491 _cacheMgr.setObj(nullptr, nullptr, cache, nullptr);
492
493 return cache;
494}
495
497 std::unique_ptr<RooChangeTracker> &&tracker, const RooArgList &flist)
498 : _sum(std::move(sumFunc)), _tracker(std::move(tracker))
499{
500 _frac.add(flist);
501}
502
503//_____________________________________________________________________________
505{
506 return RooArgList(*_sum, *_tracker);
507}
508
509//_____________________________________________________________________________
511
512//_____________________________________________________________________________
514{
515 // Special version of getValV() overrides Base_t::getValV() to save value of current normalization set
516 _curNormSet = set ? const_cast<RooArgSet *>(set) : const_cast<RooArgSet *>(static_cast<RooArgSet const*>(&_obsList));
517 return Base_t::getValV(set);
518}
519
520//_____________________________________________________________________________
522{
523 CacheElem *cache = getCache(nset ? nset : _curNormSet);
524
525 if (cache->_tracker->hasChanged(true)) {
526 cache->calculateFractions(*this, false); // verbose turned off
527 }
528 return cache->_sum.get();
529}
530
531//_____________________________________________________________________________
533{
535
536 if (cache->_tracker->hasChanged(true)) {
537 cache->calculateFractions(*this, false); // verbose turned off
538 }
539
540 double ret = cache->_sum->getVal(_obsList.nset());
541
542 return ret;
543}
544
545//_____________________________________________________________________________
547{
548 return static_cast<RooRealVar *>(_frac.at(i));
549}
550
551//_____________________________________________________________________________
553{
554 return static_cast<RooRealVar *>(_frac.at(i));
555}
556
557//_____________________________________________________________________________
559{
560 int nPdf = self._pdfList.size();
561 int nPar = self._parList.size();
562
563 double fracLinear(1.);
564 double fracNonLinear(1.);
565
567 // Calculate the delta vector
568 vector<double> dm2;
569 for (int idim = 0; idim < nPar; idim++) {
570 double delta = (static_cast<RooRealVar *>(self._parList.at(idim)))->getVal() - self._referenceGrid._nref[0][idim];
571 dm2.push_back(delta);
572 }
573
574 vector<vector<int>> powers;
575 for (int idim = 0; idim < nPar; idim++) {
576 vector<int> xtmp;
577 xtmp.reserve(self._referenceGrid._nnuis[idim]);
578 for (int ix = 0; ix < self._referenceGrid._nnuis[idim]; ix++) {
579 xtmp.push_back(ix);
580 }
581 powers.push_back(xtmp);
582 }
583
584 vector<vector<int>> output;
586 int nCombs = output.size();
587
588 vector<double> deltavec(nPdf, 1.0);
589
590 int nperm = 0;
591 for (int i = 0; i < nCombs; i++) {
592 double tmpDm = 1.0;
593 for (int ix = 0; ix < nPar; ix++) {
594 double delta = dm2[ix];
595 tmpDm *= TMath::Power(delta, static_cast<double>(output[i][ix]));
596 }
597 deltavec[nperm] = tmpDm;
598 nperm++;
599 }
600
601 double sumposfrac = 0.0;
602 for (int i = 0; i < nPdf; ++i) {
603 double ffrac = 0.0;
604
605 for (int j = 0; j < nPdf; ++j) {
606 ffrac += (*self._M)(j, i) * deltavec[j] * fracNonLinear;
607 }
608
609 if (ffrac >= 0) {
610 sumposfrac += ffrac;
611 }
612
613 // fractions for pdf
614 if (self._setting != NonLinearLinFractions) {
615 const_cast<RooRealVar *>(frac(i))->setVal(ffrac);
616 }
617
618 // fractions for rms and mean
619 const_cast<RooRealVar *>(frac(nPdf + i))->setVal(ffrac); // need to add up
620 const_cast<RooRealVar *>(frac(2 * nPdf + i))->setVal(ffrac); // need to add up
621
622 if (verbose) {
623 cout << "NonLinear fraction " << ffrac << endl;
624 frac(i)->Print();
625 frac(nPdf + i)->Print();
626 frac(2 * nPdf + i)->Print();
627 }
628 }
629
630 if (self._setting == NonLinearPosFractions) {
631 for (int i = 0; i < nPdf; ++i) {
632 if (frac(i)->getVal() < 0)
633 const_cast<RooRealVar *>(frac(i))->setVal(0.);
634 const_cast<RooRealVar *>(frac(i))->setVal(frac(i)->getVal() / sumposfrac);
635 }
636 }
637 }
638
639 if (self._setting == Linear || self._setting == NonLinearLinFractions) {
640 // zero all fractions
641 // for (int i = 0; i < 3*nPdf; ++i) {
642 for (int i = 0; i < nPdf; ++i) {
643 double initval = 0;
644 const_cast<RooRealVar *>(frac(i))->setVal(initval);
645 const_cast<RooRealVar *>(frac(nPdf + i))->setVal(initval);
646 const_cast<RooRealVar *>(frac(2 * nPdf + i))->setVal(initval);
647 }
648
649 std::vector<double> mtmp;
650
651 // loop over parList
652 for (auto *m : static_range_cast<RooRealVar *>(self._parList)) {
653 mtmp.push_back(m->getVal());
654 }
655
656 self.findShape(mtmp); // this sets _squareVec and _squareIdx quantities
657
658 int depth = TMath::Power(2, nPar);
659 vector<double> deltavec(depth, 1.0);
660
661 int nperm = 0;
662
663 vector<int> xtmp;
664 xtmp.reserve(nPar);
665 for (int ix = 0; ix < nPar; ix++) {
666 xtmp.push_back(ix);
667 }
668
669 for (int iperm = 1; iperm <= nPar; ++iperm) {
670 do {
671 double dtmp = mtmp[xtmp[0]] - self._squareVec[0][xtmp[0]];
672 for (int itmp = 1; itmp < iperm; ++itmp) {
673 dtmp *= mtmp[xtmp[itmp]] - self._squareVec[0][xtmp[itmp]];
674 }
675 deltavec[nperm + 1] = dtmp;
676 nperm++;
677 } while (RooFit::Detail::nextCombination(xtmp.begin(), xtmp.begin() + iperm, xtmp.end()));
678 }
679
680 double origFrac1(0.);
681 double origFrac2(0.);
682 for (int i = 0; i < depth; ++i) {
683 double ffrac = 0.;
684 for (int j = 0; j < depth; ++j) {
685 ffrac += (*self._MSqr)(j, i) * deltavec[j] * fracLinear;
686 }
687
688 // set fractions for pdf
689 origFrac1 = frac(self._squareIdx[i])->getVal(); // already set in case of smoothlinear
690 const_cast<RooRealVar *>(frac(self._squareIdx[i]))->setVal(origFrac1 + ffrac); // need to add up
691
692 // set fractions for rms and mean
693 if (self._setting != NonLinearLinFractions) {
694 origFrac2 =
695 frac(nPdf + self._squareIdx[i])->getVal(); // already set in case of smoothlinear
696 const_cast<RooRealVar *>(frac(nPdf + self._squareIdx[i]))->setVal(origFrac2 + ffrac); // need to add up
697 const_cast<RooRealVar *>(frac(2 * nPdf + self._squareIdx[i]))->setVal(origFrac2 + ffrac); // need to add up
698 }
699
700 if (verbose) {
701 cout << "Linear fraction " << ffrac << endl;
702 frac(self._squareIdx[i])->Print();
703 frac(nPdf + self._squareIdx[i])->Print();
704 frac(2 * nPdf + self._squareIdx[i])->Print();
705 }
706 }
707 }
708}
709
710//_____________________________________________________________________________
711void RooMomentMorphFuncND::findShape(const vector<double> &x) const
712{
713 int nPar = _parList.size();
714 int nRef = _referenceGrid._nref.size();
715
716 // Find hypercube enclosing the location to morph to
717 // bool isEnclosed = true;
718 // for (int i = 0; i < nPar; i++) {
719 // if (x[i] < _referenceGrid._grid[i]->lowBound())
720 // isEnclosed = false;
721 // if (x[i] > _referenceGrid._grid[i]->highBound())
722 // isEnclosed = false;
723 // }
724
725 // cout << "isEnclosed = " << isEnclosed << endl;
726
727 int depth = TMath::Power(2, nPar);
728
729 vector<vector<double>> boundaries(nPar);
730 for (int idim = 0; idim < nPar; idim++) {
731 int bin = _referenceGrid._grid[idim]->binNumber(x[idim]);
732 double lo = _referenceGrid._grid[idim]->binLow(bin);
733 double hi = _referenceGrid._grid[idim]->binHigh(bin);
734 boundaries[idim].push_back(lo);
735 boundaries[idim].push_back(hi);
736 }
737
738 vector<vector<double>> output;
741
742 for (int isq = 0; isq < depth; isq++) {
743 for (int iref = 0; iref < nRef; iref++) {
744 if (_squareVec[isq] == _referenceGrid._nref[iref]) {
745 _squareIdx[isq] = iref;
746 break;
747 }
748 }
749 }
750
751 // cout << endl;
752
753 // for (int isq = 0; isq < _squareVec.size(); isq++) {
754 // cout << _squareIdx[isq];
755 // cout << " (";
756 // for (int isqq = 0; isqq < _squareVec[isq].size(); isqq++) {
757 // cout << _squareVec[isq][isqq] << ((isqq<_squareVec[isq].size()-1)?",":"");
758 // }
759 // cout << ") ";
760 // }
761
762 // construct transformation matrix for linear extrapolation
763 TMatrixD M(depth, depth);
764
765 vector<int> xtmp;
766 xtmp.reserve(nPar);
767 for (int ix = 0; ix < nPar; ix++) {
768 xtmp.push_back(ix);
769 }
770
771 for (int k = 0; k < depth; ++k) {
772 M(k, 0) = 1.0;
773
774 int nperm = 0;
775 vector<double> squareBase = _squareVec[0];
776
777 for (int iperm = 1; iperm <= nPar; ++iperm) {
778 do {
779 double dtmp = _squareVec[k][xtmp[0]] - squareBase[xtmp[0]];
780 for (int itmp = 1; itmp < iperm; ++itmp) {
781 dtmp *= _squareVec[k][xtmp[itmp]] - squareBase[xtmp[itmp]];
782 }
783 M(k, nperm + 1) = dtmp;
784 nperm++;
785 } while (RooFit::Detail::nextCombination(xtmp.begin(), xtmp.begin() + iperm, xtmp.end()));
786 }
787 }
788
789 // M.Print();
790 (*_MSqr) = M.Invert();
791}
792
793//_____________________________________________________________________________
795{
796 if (allVars.size() == 1) {
797 RooAbsReal *temp = const_cast<RooMomentMorphFuncND *>(this);
798 temp->specialIntegratorConfig(true)->method1D().setLabel("RooBinIntegrator");
799 int nbins = (static_cast<RooRealVar *>(allVars.first()))->numBins();
800 temp->specialIntegratorConfig(true)->getConfigSection("RooBinIntegrator").setRealValue("numBins", nbins);
801 return true;
802 } else {
803 cout << "Currently BinIntegrator only knows how to deal with 1-d " << endl;
804 return false;
805 }
806 return false;
807}
#define coutW(a)
#define coutE(a)
#define TRACE_DESTROY
Definition RooTrace.h:24
#define TRACE_CREATE
Definition RooTrace.h:23
#define ClassImp(name)
Definition Rtypes.h:377
char name[80]
Definition TGX11.cxx:110
#define hi
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
void setLocalNoDirtyInhibit(bool flag) const
Definition RooAbsArg.h:700
Abstract base class for RooRealVar binning definitions.
Abstract container object that can hold multiple RooAbsArg objects.
virtual bool add(const RooAbsArg &var, bool silent=false)
Add the specified argument to list.
Storage_t::size_type size() const
RooAbsArg * first() const
bool addTyped(const RooAbsCollection &list, bool silent=false)
Adds elements of a given RooAbsCollection to the container if they match the specified type.
bool setRealValue(const char *name, double newVal=0.0, bool verbose=false)
Set value of a RooAbsRealLValue stored in set with given name to newVal No error messages are printed...
RooAbsReal * mean()
const RooArgSet * nset() const
Definition RooAbsProxy.h:52
Abstract base class for objects that represent a real value and implements functionality common to al...
Definition RooAbsReal.h:59
double getVal(const RooArgSet *normalisationSet=nullptr) const
Evaluate object.
Definition RooAbsReal.h:103
virtual double getValV(const RooArgSet *normalisationSet=nullptr) const
Return value of object.
RooNumIntConfig * specialIntegratorConfig() const
Returns the specialized integrator configuration for this RooAbsReal.
RooAbsMoment * sigma(RooRealVar &obs)
Definition RooAbsReal.h:368
Calculates the sum of a set of RooAbsReal terms, or when constructed with two sets,...
Definition RooAddition.h:27
RooArgList is a container object that can hold multiple RooAbsArg objects.
Definition RooArgList.h:22
RooAbsArg * at(Int_t idx) const
Return object at given index, or nullptr if index is out of range.
Definition RooArgList.h:110
RooArgSet is a container object that can hold multiple RooAbsArg objects.
Definition RooArgSet.h:55
Implements a RooAbsBinning in terms of an array of boundary values, posing no constraints on the choi...
Definition RooBinning.h:27
double * array() const override
Return array of boundary values.
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.
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...
Represents a constant real-valued object.
Definition RooConstVar.h:23
RooCustomizer is a factory class to produce clones of a prototype composite PDF object with the same ...
void replaceArg(const RooAbsArg &orig, const RooAbsArg &subst)
Replace any occurrence of arg 'orig' with arg 'subst'.
RooAbsArg * build(const char *masterCatState, bool verbose=false)
Build a clone of the prototype executing all registered 'replace' rules and 'split' rules for the mas...
A RooFormulaVar is a generic implementation of a real-valued object, which takes a RooArgList of serv...
RooLinearVar is the most general form of a derived real-valued object that can be used by RooRealInte...
void calculateFractions(const RooMomentMorphFuncND &self, bool verbose=true) const
std::unique_ptr< RooChangeTracker > _tracker
std::unique_ptr< RooAbsReal > _sum
RooArgList containedArgs(Action) override
CacheElem(std::unique_ptr< RooAbsReal > &&sumFunc, std::unique_ptr< RooChangeTracker > &&tracker, const RooArgList &flist)
void addBinning(const RooAbsBinning &binning)
std::vector< RooAbsBinning * > _grid
std::vector< std::vector< double > > _nref
void addPdf(const RooAbsReal &func, int bin_x)
RooObjCacheManager _cacheMgr
! Transient cache manager
std::unique_ptr< TMatrixD > _MSqr
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.
RooAbsReal * sumFunc(const RooArgSet *nset)
CacheElem * getCache(const RooArgSet *nset) const
std::unique_ptr< TMatrixD > _M
RooArgSet * _curNormSet
! Transient cache manager
bool setBinIntegrator(RooArgSet &allVars)
double getValV(const RooArgSet *set=nullptr) const override
Return value of object.
std::vector< std::vector< double > > _squareVec
int sij(const int &i, const int &j) const
std::vector< int > _squareIdx
const RooArgSet & getConfigSection(const char *name) const
Retrieve configuration information specific to integrator with given name.
RooCategory & method1D()
Variable that can be changed from the outside.
Definition RooRealVar.h:37
TMatrixT< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant.
const char * GetName() const override
Returns name of object.
Definition TNamed.h:47
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
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:721
TMarker m
Definition textangle.C:8
static void output()