Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
RooMomentMorphFuncND.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * *
4  * This code was autogenerated by RooClassFactory *
5  *****************************************************************************/
6 
7 #include "Riostream.h"
8 
9 #include "RooMomentMorphFuncND.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 "RooRealSumFunc.h"
17 #include "RooAddition.h"
18 #include "RooAbsMoment.h"
19 #include "RooMoment.h"
20 #include "RooLinearVar.h"
21 #include "RooChangeTracker.h"
22 #include "RooNumIntConfig.h"
23 #include "RooHistPdf.h"
24 
25 #include "TMath.h"
26 #include "TVector.h"
27 #include "TMap.h"
28 
29 #include <map>
30 #include <algorithm>
31 
32 using namespace std;
33 
34 // ClassImp for RooMomentMorphFuncND
35 ClassImp(RooMomentMorphFuncND)
36 
37  // ClassImp for RooMomentMorphFuncND
38  ClassImp(RooMomentMorphFuncND::Grid2)
39 
40  //_____________________________________________________________________________
41  RooMomentMorphFuncND::RooMomentMorphFuncND()
42  : _curNormSet(0), _M(0), _MSqr(0), _setting(RooMomentMorphFuncND::Linear), _useHorizMorph(true)
43 {
44  _parItr = _parList.createIterator();
45  _obsItr = _obsList.createIterator();
46 
47  TRACE_CREATE
48 }
49 
50 //_____________________________________________________________________________
51 RooMomentMorphFuncND::RooMomentMorphFuncND(const char *name, const char *title, const RooArgList &parList,
52  const RooArgList &obsList, const Grid2 &referenceGrid,
53  const Setting &setting)
54  : RooAbsReal(name, title), _cacheMgr(this, 10, kTRUE, kTRUE), _parList("parList", "List of morph parameters", this),
55  _obsList("obsList", "List of observables", this), _referenceGrid(referenceGrid),
56  _pdfList("pdfList", "List of pdfs", this), _setting(setting), _useHorizMorph(true)
57 {
58  // morph parameters
59  initializeParameters(parList);
60 
61  // observables
62  initializeObservables(obsList);
63 
64  _pdfList.add(_referenceGrid._pdfList);
65 
66  // general initialization
67  initialize();
68 
69  TRACE_CREATE
70 }
71 
72 //_____________________________________________________________________________
73 RooMomentMorphFuncND::RooMomentMorphFuncND(const char *name, const char *title, RooAbsReal &_m,
74  const RooArgList &varList, const RooArgList &pdfList,
75  const TVectorD &mrefpoints, Setting setting)
76  : RooAbsReal(name, title), _cacheMgr(this, 10, kTRUE, kTRUE), _parList("parList", "List of morph parameters", this),
77  _obsList("obsList", "List of observables", this), _pdfList("pdfList", "List of pdfs", this), _setting(setting),
78  _useHorizMorph(true)
79 {
80  // make reference grid
81  RooBinning grid(mrefpoints.GetNrows() - 1, mrefpoints.GetMatrixArray());
82  _referenceGrid.addBinning(grid);
83 
84  for (int i = 0; i < mrefpoints.GetNrows(); ++i) {
85  for (int j = 0; j < grid.numBoundaries(); ++j) {
86  if (mrefpoints[i] == grid.array()[j]) {
87  _referenceGrid.addPdf(*(RooAbsReal *)pdfList.at(i), j);
88  break;
89  }
90  }
91  }
92 
93  _pdfList.add(_referenceGrid._pdfList);
94 
95  // morph parameters
96  RooArgList parList;
97  parList.add(_m);
98  initializeParameters(parList);
99 
100  // observables
101  initializeObservables(varList);
102 
103  // general initialization
104  initialize();
105 
106  TRACE_CREATE
107 }
108 
109 //_____________________________________________________________________________
110 RooMomentMorphFuncND::RooMomentMorphFuncND(const char *name, const char *title, RooAbsReal &_m,
111  const RooArgList &varList, const RooArgList &pdfList,
112  const RooArgList &mrefList, Setting setting)
113  : RooAbsReal(name, title), _cacheMgr(this, 10, kTRUE, kTRUE), _parList("parList", "List of morph parameters", this),
114  _obsList("obsList", "List of observables", this), _pdfList("pdfList", "List of pdfs", this), _setting(setting),
115  _useHorizMorph(true)
116 {
117  // make reference grid
118  TVectorD mrefpoints(mrefList.getSize());
119  TIterator *mrefItr = mrefList.createIterator();
120  RooAbsReal *mref;
121  for (int i = 0; (mref = dynamic_cast<RooAbsReal *>(mrefItr->Next())); ++i) {
122  if (!mref) {
123  coutE(InputArguments) << "RooMomentMorphFuncND::ctor(" << GetName() << ") ERROR: mref " << mref->GetName()
124  << " is not of type RooAbsReal" << endl;
125  throw string("RooMomentMorphFuncND::ctor() ERROR mref is not of type RooAbsReal");
126  }
127  if (!dynamic_cast<RooConstVar *>(mref)) {
128  coutW(InputArguments) << "RooMomentMorphFuncND::ctor(" << GetName() << ") WARNING mref point " << i
129  << " is not a constant, taking a snapshot of its value" << endl;
130  }
131  mrefpoints[i] = mref->getVal();
132  }
133  delete mrefItr;
134 
135  RooBinning grid(mrefpoints.GetNrows() - 1, mrefpoints.GetMatrixArray());
136  _referenceGrid.addBinning(grid);
137 
138  for (int i = 0; i < mrefpoints.GetNrows(); ++i) {
139  for (int j = 0; j < grid.numBoundaries(); ++j) {
140  if (mrefpoints[i] == grid.array()[j]) {
141  _referenceGrid.addPdf(*(RooAbsReal *)pdfList.at(i), j);
142  break;
143  }
144  }
145  }
146 
147  _pdfList.add(_referenceGrid._pdfList);
148 
149  // morph parameters
150  RooArgList parList;
151  parList.add(_m);
152  initializeParameters(parList);
153 
154  // observables
155  initializeObservables(varList);
156 
157  // general initialization
158  initialize();
159 
160  TRACE_CREATE
161 }
162 
163 //_____________________________________________________________________________
164 RooMomentMorphFuncND::RooMomentMorphFuncND(const RooMomentMorphFuncND &other, const char *name)
165  : RooAbsReal(other, name), _cacheMgr(other._cacheMgr, this), _curNormSet(0),
166  _parList("parList", this, other._parList), _obsList("obsList", this, other._obsList),
167  _referenceGrid(other._referenceGrid), _pdfList("pdfList", this, other._pdfList), _M(0), _MSqr(0),
168  _setting(other._setting), _useHorizMorph(other._useHorizMorph)
169 {
170  _parItr = _parList.createIterator();
171  _obsItr = _obsList.createIterator();
172 
173  // general initialization
174  initialize();
175 
176  TRACE_CREATE
177 }
178 
179 //_____________________________________________________________________________
180 RooMomentMorphFuncND::~RooMomentMorphFuncND()
181 {
182  if (_parItr)
183  delete _parItr;
184  if (_obsItr)
185  delete _obsItr;
186  if (_M)
187  delete _M;
188  if (_MSqr)
189  delete _MSqr;
190 
191  TRACE_DESTROY
192 }
193 
194 //_____________________________________________________________________________
195 void RooMomentMorphFuncND::initializeParameters(const RooArgList &parList)
196 {
197  TIterator *parItr = parList.createIterator();
198  RooAbsArg *par;
199  for (int i = 0; (par = (RooAbsArg *)parItr->Next()); ++i) {
200  if (!dynamic_cast<RooAbsReal *>(par)) {
201  coutE(InputArguments) << "RooMomentMorphFuncND::ctor(" << GetName() << ") ERROR: parameter " << par->GetName()
202  << " is not of type RooAbsReal" << endl;
203  throw string("RooMomentMorphFuncND::initializeParameters() ERROR parameter is not of type RooAbsReal");
204  }
205  _parList.add(*par);
206  }
207  delete parItr;
208 
209  _parItr = _parList.createIterator();
210 }
211 
212 //_____________________________________________________________________________
213 void RooMomentMorphFuncND::initializeObservables(const RooArgList &obsList)
214 {
215  TIterator *obsItr = obsList.createIterator();
216  RooAbsArg *var;
217  for (int i = 0; (var = (RooAbsArg *)obsItr->Next()); ++i) {
218  if (!dynamic_cast<RooAbsReal *>(var)) {
219  coutE(InputArguments) << "RooMomentMorphFuncND::ctor(" << GetName() << ") ERROR: variable " << var->GetName()
220  << " is not of type RooAbsReal" << endl;
221  throw string("RooMomentMorphFuncND::initializeObservables() ERROR variable is not of type RooAbsReal");
222  }
223  _obsList.add(*var);
224  }
225  delete obsItr;
226 
227  _obsItr = _obsList.createIterator();
228 }
229 
230 //_____________________________________________________________________________
231 // from http://stackoverflow.com/a/5279601
232 template <typename T>
233 void RooMomentMorphFuncND::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 //_____________________________________________________________________________
266 void RooMomentMorphFuncND::initialize()
267 {
268  // TIterator* pdfItr = _referenceGrid._pdfList.createIterator() ;
269  // RooAbsReal* pdf ;
270  // for (int i=0; (pdf = dynamic_cast<RooAbsReal*>(pdfItr->Next())); ++i) {
271  // if (!pdf) {
272  // coutE(InputArguments) << "RooMomentMorphFunc::ctor(" << GetName() << ") ERROR: pdf " << pdf->GetName() << " is
273  // not of type RooAbsReal" << endl ;
274  // throw string("RooPolyMorh::ctor() ERROR pdf is not of type RooAbsReal") ;
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) << "RooMomentMorphFuncND::initialize(" << GetName() << ") ERROR: nPar != nDim"
293  << ": " << nPar << " !=" << nDim << endl;
294  assert(0);
295  }
296 
297  if (nPdf != nRef) {
298  coutE(InputArguments) << "RooMomentMorphFuncND::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);
306  if (_setting == NonLinear || _setting == NonLinearPosFractions || _setting == NonLinearLinFractions) {
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 //_____________________________________________________________________________
355 RooMomentMorphFuncND::Grid2::Grid2(const RooMomentMorphFuncND::Grid2 &other)
356  : _grid(other._grid), _pdfList(other._pdfList), _pdfMap(other._pdfMap), _nref(other._nref)
357 {
358 }
359 
360 //_____________________________________________________________________________
361 RooMomentMorphFuncND::Grid2::~Grid2()
362 {
363 }
364 
365 //_____________________________________________________________________________
366 void RooMomentMorphFuncND::Grid2::addPdf(const RooAbsReal &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 //_____________________________________________________________________________
378 void RooMomentMorphFuncND::Grid2::addPdf(const RooAbsReal &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 //_____________________________________________________________________________
392 void RooMomentMorphFuncND::Grid2::addPdf(const RooAbsReal &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 //_____________________________________________________________________________
408 void RooMomentMorphFuncND::Grid2::addPdf(const RooAbsReal &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 //_____________________________________________________________________________
421 RooMomentMorphFuncND::CacheElem *RooMomentMorphFuncND::getCache(const RooArgSet * /*nset*/) const
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 
431  TIterator *pdfItr = _pdfList.createIterator();
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 *> offsets(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  RooRealSumFunc *theSumFunc = 0;
467  string sumfuncName = Form("%s_sumfunc", 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
475  ? ((RooAbsReal *)_pdfList.at(i))->sigma((RooRealVar &)*obsList.at(j))
476  : ((RooAbsReal *)_pdfList.at(i))->sigma((RooRealVar &)*obsList.at(j), obsList);
477 
478  mom->setLocalNoDirtyInhibit(kTRUE);
479  mom->mean()->setLocalNoDirtyInhibit(kTRUE);
480 
481  sigmarv[sij(i, j)] = mom;
482  meanrv[sij(i, j)] = mom->mean();
483 
484  ownedComps.add(*sigmarv[sij(i, j)]);
485  }
486  }
487 
488  // slope and offsets (to be set later, depend on nuisance parameters)
489  for (int j = 0; j < nObs; ++j) {
490  RooArgList meanList("meanList");
491  RooArgList rmsList("rmsList");
492  for (int i = 0; i < nPdf; ++i) {
493  meanList.add(*meanrv[sij(i, j)]);
494  rmsList.add(*sigmarv[sij(i, j)]);
495  }
496  string myrmsName = Form("%s_rms_%d", GetName(), j);
497  string myposName = Form("%s_pos_%d", GetName(), j);
498  mypos[j] = new RooAddition(myposName.c_str(), myposName.c_str(), meanList, coefList2);
499  myrms[j] = new RooAddition(myrmsName.c_str(), myrmsName.c_str(), rmsList, coefList3);
500  ownedComps.add(RooArgSet(*myrms[j], *mypos[j]));
501  }
502 
503  // construction of unit pdfs
504  pdfItr->Reset();
505  RooAbsReal *pdf;
506  RooArgList transPdfList;
507 
508  for (int i = 0; i < nPdf; ++i) {
509  _obsItr->Reset();
510  RooRealVar *var;
511 
512  pdf = (RooAbsReal *)pdfItr->Next();
513  string pdfName = Form("pdf_%d", i);
514  RooCustomizer cust(*pdf, pdfName.c_str());
515 
516  for (int j = 0; j < nObs; ++j) {
517  // slope and offsets formulas
518  string slopeName = Form("%s_slope_%d_%d", GetName(), i, j);
519  string offsetName = Form("%s_offset_%d_%d", GetName(), i, j);
520 
521  slope[sij(i, j)] =
522  new RooFormulaVar(slopeName.c_str(), "@0/@1", RooArgList(*sigmarv[sij(i, j)], *myrms[j]));
523  offsets[sij(i, j)] = new RooFormulaVar(offsetName.c_str(), "@0-(@1*@2)",
524  RooArgList(*meanrv[sij(i, j)], *mypos[j], *slope[sij(i, j)]));
525  ownedComps.add(RooArgSet(*slope[sij(i, j)], *offsets[sij(i, j)]));
526 
527  // linear transformations, so pdf can be renormalized easily
528  var = (RooRealVar *)(_obsItr->Next());
529  string transVarName = Form("%s_transVar_%d_%d", GetName(), i, j);
530  transVar[sij(i, j)] = new RooLinearVar(transVarName.c_str(), transVarName.c_str(), *var, *slope[sij(i, j)],
531  *offsets[sij(i, j)]);
532 
533  // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
534  // This will prevent the likelihood optimizers from erroneously declaring terms constant
535  transVar[sij(i, j)]->addServerList((RooAbsCollection &)_parList);
536 
537  ownedComps.add(*transVar[sij(i, j)]);
538  cust.replaceArg(*var, *transVar[sij(i, j)]);
539  }
540  transPdf[i] = (RooAbsReal *)cust.build();
541  transPdfList.add(*transPdf[i]);
542  ownedComps.add(*transPdf[i]);
543  }
544 
545  // sum pdf
546  theSumFunc = new RooRealSumFunc(sumfuncName.c_str(), sumfuncName.c_str(), transPdfList, coefList);
547  } else {
548  theSumFunc = new RooRealSumFunc(sumfuncName.c_str(), sumfuncName.c_str(), _pdfList, coefList);
549  }
550 
551  // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
552  // This will prevent the likelihood optimizers from erroneously declaring terms constant
553  theSumFunc->addServerList((RooAbsCollection &)_parList);
554  theSumFunc->addOwnedComponents(ownedComps);
555 
556  // change tracker for fraction parameters
557  string trackerName = Form("%s_frac_tracker", GetName());
558  RooChangeTracker *tracker = new RooChangeTracker(trackerName.c_str(), trackerName.c_str(), _parList, kTRUE);
559 
560  // Store it in the cache
561  cache = new CacheElem(*theSumFunc, *tracker, fracl);
562  _cacheMgr.setObj(0, 0, cache, 0);
563 
564  return cache;
565 }
566 
567 //_____________________________________________________________________________
568 RooArgList RooMomentMorphFuncND::CacheElem::containedArgs(Action)
569 {
570  return RooArgList(*_sumFunc, *_tracker);
571 }
572 
573 //_____________________________________________________________________________
574 RooMomentMorphFuncND::CacheElem::~CacheElem()
575 {
576  delete _sumFunc;
577  delete _tracker;
578 }
579 
580 //_____________________________________________________________________________
581 Double_t RooMomentMorphFuncND::getVal(const RooArgSet *set) const
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 RooAbsReal::getVal(set);
586 }
587 
588 //_____________________________________________________________________________
589 RooAbsReal *RooMomentMorphFuncND::sumFunc(const RooArgSet *nset)
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->_sumFunc;
597 }
598 
599 //_____________________________________________________________________________
600 Double_t RooMomentMorphFuncND::evaluate() const
601 {
602  CacheElem *cache = getCache(_curNormSet);
603 
604  if (cache->_tracker->hasChanged(kTRUE)) {
605  cache->calculateFractions(*this, kFALSE); // verbose turned off
606  }
607 
608  Double_t ret = cache->_sumFunc->getVal(_obsList.nset());
609 
610  return ret;
611 }
612 
613 //_____________________________________________________________________________
614 RooRealVar *RooMomentMorphFuncND::CacheElem::frac(int i)
615 {
616  return (RooRealVar *)(_frac.at(i));
617 }
618 
619 //_____________________________________________________________________________
620 const RooRealVar *RooMomentMorphFuncND::CacheElem::frac(int i) const
621 {
622  return (RooRealVar *)(_frac.at(i));
623 }
624 
625 //_____________________________________________________________________________
626 // from http://stackoverflow.com/a/5097100/8747
627 template <typename Iterator>
628 bool RooMomentMorphFuncND::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 //_____________________________________________________________________________
665 void RooMomentMorphFuncND::CacheElem::calculateFractions(const RooMomentMorphFuncND &self, Bool_t verbose) const
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 
673  if (self._setting == NonLinear || self._setting == NonLinearLinFractions || self._setting == NonLinearPosFractions) {
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 //_____________________________________________________________________________
818 void RooMomentMorphFuncND::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  (void)isEnclosed; // Silence warning for isEnclosed being unused.
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);
847  _squareVec = output;
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 //_____________________________________________________________________________
900 Bool_t RooMomentMorphFuncND::setBinIntegrator(RooArgSet &allVars)
901 {
902  if (allVars.getSize() == 1) {
903  RooAbsReal *temp = const_cast<RooMomentMorphFuncND *>(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 }