Logo ROOT   6.30.04
Reference Guide
 All Namespaces Files Pages
RooMomentMorph.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * *
4  * This code was autogenerated by RooClassFactory *
5  *****************************************************************************/
6 
7 /** \class RooMomentMorph
8  \ingroup Roofit
9 
10 **/
11 
12 #include "Riostream.h"
13 
14 #include "RooMomentMorph.h"
15 #include "RooAbsCategory.h"
16 #include "RooRealIntegral.h"
17 #include "RooRealConstant.h"
18 #include "RooRealVar.h"
19 #include "RooFormulaVar.h"
20 #include "RooCustomizer.h"
21 #include "RooAddPdf.h"
22 #include "RooAddition.h"
23 #include "RooMoment.h"
24 #include "RooLinearVar.h"
25 #include "RooChangeTracker.h"
26 
27 #include "TMath.h"
28 #include "TH1.h"
29 
30 using namespace std;
31 
32 ClassImp(RooMomentMorph);
33 
34 ////////////////////////////////////////////////////////////////////////////////
35 /// coverity[UNINIT_CTOR]
36 
37 RooMomentMorph::RooMomentMorph() : _curNormSet(0), _mref(0), _M(0), _useHorizMorph(true)
38 {
39  _varItr = _varList.createIterator() ;
40  _pdfItr = _pdfList.createIterator() ;
41 }
42 
43 ////////////////////////////////////////////////////////////////////////////////
44 /// CTOR
45 
46 RooMomentMorph::RooMomentMorph(const char *name, const char *title,
47  RooAbsReal& _m,
48  const RooArgList& varList,
49  const RooArgList& pdfList,
50  const TVectorD& mrefpoints,
51  Setting setting) :
52  RooAbsPdf(name,title),
53  _cacheMgr(this,10,kTRUE,kTRUE),
54  m("m","m",this,_m),
55  _varList("varList","List of variables",this),
56  _pdfList("pdfList","List of pdfs",this),
57  _setting(setting),
58  _useHorizMorph(true)
59 {
60  // observables
61  TIterator* varItr = varList.createIterator() ;
62  RooAbsArg* var ;
63  for (Int_t i=0; (var = (RooAbsArg*)varItr->Next()); ++i) {
64  if (!dynamic_cast<RooAbsReal*>(var)) {
65  coutE(InputArguments) << "RooMomentMorph::ctor(" << GetName() << ") ERROR: variable " << var->GetName() << " is not of type RooAbsReal" << endl ;
66  throw string("RooPolyMorh::ctor() ERROR variable is not of type RooAbsReal") ;
67  }
68  _varList.add(*var) ;
69  }
70  delete varItr ;
71 
72  // reference p.d.f.s
73  TIterator* pdfItr = pdfList.createIterator() ;
74  RooAbsPdf* pdf ;
75  for (Int_t i=0; (pdf = dynamic_cast<RooAbsPdf*>(pdfItr->Next())); ++i) {
76  if (!pdf) {
77  coutE(InputArguments) << "RooMomentMorph::ctor(" << GetName() << ") ERROR: pdf " << pdf->GetName() << " is not of type RooAbsPdf" << endl ;
78  throw string("RooPolyMorh::ctor() ERROR pdf is not of type RooAbsPdf") ;
79  }
80  _pdfList.add(*pdf) ;
81  }
82  delete pdfItr ;
83 
84  _mref = new TVectorD(mrefpoints);
85  _varItr = _varList.createIterator() ;
86  _pdfItr = _pdfList.createIterator() ;
87 
88  // initialization
89  initialize();
90 }
91 
92 ////////////////////////////////////////////////////////////////////////////////
93 /// CTOR
94 
95 RooMomentMorph::RooMomentMorph(const char *name, const char *title,
96  RooAbsReal& _m,
97  const RooArgList& varList,
98  const RooArgList& pdfList,
99  const RooArgList& mrefList,
100  Setting setting) :
101  RooAbsPdf(name,title),
102  _cacheMgr(this,10,kTRUE,kTRUE),
103  m("m","m",this,_m),
104  _varList("varList","List of variables",this),
105  _pdfList("pdfList","List of pdfs",this),
106  _setting(setting),
107  _useHorizMorph(true)
108 {
109  // observables
110  TIterator* varItr = varList.createIterator() ;
111  RooAbsArg* var ;
112  for (Int_t i=0; (var = (RooAbsArg*)varItr->Next()); ++i) {
113  if (!dynamic_cast<RooAbsReal*>(var)) {
114  coutE(InputArguments) << "RooMomentMorph::ctor(" << GetName() << ") ERROR: variable " << var->GetName() << " is not of type RooAbsReal" << endl ;
115  throw string("RooPolyMorh::ctor() ERROR variable is not of type RooAbsReal") ;
116  }
117  _varList.add(*var) ;
118  }
119  delete varItr ;
120 
121  // reference p.d.f.s
122  TIterator* pdfItr = pdfList.createIterator() ;
123  RooAbsPdf* pdf ;
124  for (Int_t i=0; (pdf = dynamic_cast<RooAbsPdf*>(pdfItr->Next())); ++i) {
125  if (!pdf) {
126  coutE(InputArguments) << "RooMomentMorph::ctor(" << GetName() << ") ERROR: pdf " << pdf->GetName() << " is not of type RooAbsPdf" << endl ;
127  throw string("RooPolyMorh::ctor() ERROR pdf is not of type RooAbsPdf") ;
128  }
129  _pdfList.add(*pdf) ;
130  }
131  delete pdfItr ;
132 
133  // reference points in m
134  _mref = new TVectorD(mrefList.getSize());
135  TIterator* mrefItr = mrefList.createIterator() ;
136  RooAbsReal* mref ;
137  for (Int_t i=0; (mref = dynamic_cast<RooAbsReal*>(mrefItr->Next())); ++i) {
138  if (!mref) {
139  coutE(InputArguments) << "RooMomentMorph::ctor(" << GetName() << ") ERROR: mref " << mref->GetName() << " is not of type RooAbsReal" << endl ;
140  throw string("RooPolyMorh::ctor() ERROR mref is not of type RooAbsReal") ;
141  }
142  if (!dynamic_cast<RooConstVar*>(mref)) {
143  coutW(InputArguments) << "RooMomentMorph::ctor(" << GetName() << ") WARNING mref point " << i << " is not a constant, taking a snapshot of its value" << endl ;
144  }
145  (*_mref)[i] = mref->getVal() ;
146  }
147  delete mrefItr ;
148 
149  _varItr = _varList.createIterator() ;
150  _pdfItr = _pdfList.createIterator() ;
151 
152  // initialization
153  initialize();
154 }
155 
156 ////////////////////////////////////////////////////////////////////////////////
157 
158 RooMomentMorph::RooMomentMorph(const RooMomentMorph& other, const char* name) :
159  RooAbsPdf(other,name),
160  _cacheMgr(other._cacheMgr,this),
161  _curNormSet(0),
162  m("m",this,other.m),
163  _varList("varList",this,other._varList),
164  _pdfList("pdfList",this,other._pdfList),
165  _setting(other._setting),
166  _useHorizMorph(other._useHorizMorph)
167 {
168  _mref = new TVectorD(*other._mref) ;
169  _varItr = _varList.createIterator() ;
170  _pdfItr = _pdfList.createIterator() ;
171 
172  // initialization
173  initialize();
174 }
175 
176 ////////////////////////////////////////////////////////////////////////////////
177 
178 RooMomentMorph::~RooMomentMorph()
179 {
180  if (_mref) delete _mref;
181  if (_varItr) delete _varItr;
182  if (_pdfItr) delete _pdfItr;
183  if (_M) delete _M;
184 }
185 
186 ////////////////////////////////////////////////////////////////////////////////
187 
188 void RooMomentMorph::initialize()
189 {
190  Int_t nPdf = _pdfList.getSize();
191 
192  // other quantities needed
193  if (nPdf!=_mref->GetNrows()) {
194  coutE(InputArguments) << "RooMomentMorph::initialize(" << GetName() << ") ERROR: nPdf != nRefPoints" << endl ;
195  assert(0) ;
196  }
197 
198  TVectorD* dm = new TVectorD(nPdf);
199  _M = new TMatrixD(nPdf,nPdf);
200 
201  // transformation matrix for non-linear extrapolation, needed in evaluate()
202  TMatrixD M(nPdf,nPdf);
203  for (Int_t i=0; i<_mref->GetNrows(); ++i) {
204  (*dm)[i] = (*_mref)[i]-(*_mref)[0];
205  M(i,0) = 1.;
206  if (i>0) M(0,i) = 0.;
207  }
208  for (Int_t i=1; i<_mref->GetNrows(); ++i) {
209  for (Int_t j=1; j<_mref->GetNrows(); ++j) {
210  M(i,j) = TMath::Power((*dm)[i],(double)j);
211  }
212  }
213  (*_M) = M.Invert();
214 
215  delete dm ;
216 }
217 
218 ////////////////////////////////////////////////////////////////////////////////
219 
220 RooMomentMorph::CacheElem* RooMomentMorph::getCache(const RooArgSet* /*nset*/) const
221 {
222  CacheElem* cache = (CacheElem*) _cacheMgr.getObj(0,(RooArgSet*)0) ;
223  if (cache) {
224  return cache ;
225  }
226  Int_t nVar = _varList.getSize();
227  Int_t nPdf = _pdfList.getSize();
228 
229  RooAbsReal* null = 0 ;
230  vector<RooAbsReal*> meanrv(nPdf*nVar,null);
231  vector<RooAbsReal*> sigmarv(nPdf*nVar,null);
232  vector<RooAbsReal*> myrms(nVar,null);
233  vector<RooAbsReal*> mypos(nVar,null);
234  vector<RooAbsReal*> slope(nPdf*nVar,null);
235  vector<RooAbsReal*> offs(nPdf*nVar,null);
236  vector<RooAbsReal*> transVar(nPdf*nVar,null);
237  vector<RooAbsReal*> transPdf(nPdf,null);
238 
239  RooArgSet ownedComps ;
240 
241  RooArgList fracl ;
242 
243  // fraction parameters
244  RooArgList coefList("coefList");
245  RooArgList coefList2("coefList2");
246  for (Int_t i=0; i<2*nPdf; ++i) {
247  std::string fracName = Form("frac_%d",i);
248 
249  RooRealVar* frac = new RooRealVar(fracName.c_str(),fracName.c_str(),1.) ;
250 
251  fracl.add(*frac); // to be set later
252  if (i<nPdf) coefList.add(*(RooRealVar*)(fracl.at(i))) ;
253  else coefList2.add(*(RooRealVar*)(fracl.at(i))) ;
254  ownedComps.add(*(RooRealVar*)(fracl.at(i))) ;
255  }
256 
257  RooAddPdf* theSumPdf = 0;
258  std::string sumpdfName = Form("%s_sumpdf",GetName());
259 
260  if (_useHorizMorph) {
261  // mean and sigma
262  RooArgList varList(_varList) ;
263  for (Int_t i=0; i<nPdf; ++i) {
264  for (Int_t j=0; j<nVar; ++j) {
265 
266  std::string meanName = Form("%s_mean_%d_%d",GetName(),i,j);
267  std::string sigmaName = Form("%s_sigma_%d_%d",GetName(),i,j);
268 
269  RooAbsMoment* mom = nVar==1 ?
270  ((RooAbsPdf*)_pdfList.at(i))->sigma((RooRealVar&)*varList.at(j)) :
271  ((RooAbsPdf*)_pdfList.at(i))->sigma((RooRealVar&)*varList.at(j),varList) ;
272 
273  mom->setLocalNoDirtyInhibit(kTRUE) ;
274  mom->mean()->setLocalNoDirtyInhibit(kTRUE) ;
275 
276  sigmarv[ij(i,j)] = mom ;
277  meanrv[ij(i,j)] = mom->mean() ;
278 
279  ownedComps.add(*sigmarv[ij(i,j)]) ;
280  }
281  }
282 
283  // slope and offset (to be set later, depend on m)
284  for (Int_t j=0; j<nVar; ++j) {
285  RooArgList meanList("meanList");
286  RooArgList rmsList("rmsList");
287  for (Int_t i=0; i<nPdf; ++i) {
288  meanList.add(*meanrv[ij(i,j)]);
289  rmsList.add(*sigmarv[ij(i,j)]);
290  }
291  std::string myrmsName = Form("%s_rms_%d",GetName(),j);
292  std::string myposName = Form("%s_pos_%d",GetName(),j);
293  myrms[j] = new RooAddition(myrmsName.c_str(),myrmsName.c_str(),rmsList,coefList2);
294  mypos[j] = new RooAddition(myposName.c_str(),myposName.c_str(),meanList,coefList2);
295  ownedComps.add(RooArgSet(*myrms[j],*mypos[j])) ;
296  }
297 
298  // construction of unit pdfs
299  _pdfItr->Reset();
300  RooAbsPdf* pdf;
301  RooArgList transPdfList;
302 
303  for (Int_t i=0; i<nPdf; ++i) {
304  _varItr->Reset() ;
305  RooRealVar* var ;
306 
307  pdf = (RooAbsPdf*)_pdfItr->Next();
308  std::string pdfName = Form("pdf_%d",i);
309  RooCustomizer cust(*pdf,pdfName.c_str());
310 
311  for (Int_t j=0; j<nVar; ++j) {
312  // slope and offset formulas
313  std::string slopeName = Form("%s_slope_%d_%d",GetName(),i,j);
314  std::string offsetName = Form("%s_offset_%d_%d",GetName(),i,j);
315  slope[ij(i,j)] = new RooFormulaVar(slopeName.c_str(),"@0/@1",RooArgList(*sigmarv[ij(i,j)],*myrms[j]));
316  offs[ij(i,j)] = new RooFormulaVar(offsetName.c_str(),"@0-(@1*@2)",RooArgList(*meanrv[ij(i,j)],*mypos[j],*slope[ij(i,j)]));
317  ownedComps.add(RooArgSet(*slope[ij(i,j)],*offs[ij(i,j)])) ;
318  // linear transformations, so pdf can be renormalized
319  var = (RooRealVar*)(_varItr->Next());
320  std::string transVarName = Form("%s_transVar_%d_%d",GetName(),i,j);
321  //transVar[ij(i,j)] = new RooFormulaVar(transVarName.c_str(),transVarName.c_str(),"@0*@1+@2",RooArgList(*var,*slope[ij(i,j)],*offs[ij(i,j)]));
322 
323  transVar[ij(i,j)] = new RooLinearVar(transVarName.c_str(),transVarName.c_str(),*var,*slope[ij(i,j)],*offs[ij(i,j)]);
324 
325  // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
326  // This will prevent the likelihood optimizers from erroneously declaring terms constant
327  transVar[ij(i,j)]->addServer((RooAbsArg&)m.arg()) ;
328 
329  ownedComps.add(*transVar[ij(i,j)]) ;
330  cust.replaceArg(*var,*transVar[ij(i,j)]);
331  }
332  transPdf[i] = (RooAbsPdf*) cust.build() ;
333  transPdfList.add(*transPdf[i]);
334  ownedComps.add(*transPdf[i]) ;
335  }
336  // sum pdf
337  theSumPdf = new RooAddPdf(sumpdfName.c_str(),sumpdfName.c_str(),transPdfList,coefList);
338  }
339  else {
340  theSumPdf = new RooAddPdf(sumpdfName.c_str(),sumpdfName.c_str(),_pdfList,coefList);
341  }
342 
343  // *** WVE this is important *** this declares that frac effectively depends on the morphing parameters
344  // This will prevent the likelihood optimizers from erroneously declaring terms constant
345  theSumPdf->addServer((RooAbsArg&)m.arg()) ;
346  theSumPdf->addOwnedComponents(ownedComps) ;
347 
348  // change tracker for fraction parameters
349  std::string trackerName = Form("%s_frac_tracker",GetName()) ;
350  RooChangeTracker* tracker = new RooChangeTracker(trackerName.c_str(),trackerName.c_str(),m.arg(),kTRUE) ;
351 
352  // Store it in the cache
353  cache = new CacheElem(*theSumPdf,*tracker,fracl) ;
354  _cacheMgr.setObj(0,0,cache,0) ;
355 
356  cache->calculateFractions(*this, kFALSE);
357  return cache ;
358 }
359 
360 ////////////////////////////////////////////////////////////////////////////////
361 
362 RooArgList RooMomentMorph::CacheElem::containedArgs(Action)
363 {
364  return RooArgList(*_sumPdf,*_tracker) ;
365 }
366 
367 ////////////////////////////////////////////////////////////////////////////////
368 
369 RooMomentMorph::CacheElem::~CacheElem()
370 {
371  delete _sumPdf ;
372  delete _tracker ;
373 }
374 
375 ////////////////////////////////////////////////////////////////////////////////
376 /// Special version of getVal() overrides RooAbsReal::getVal() to save value of current normalization set
377 
378 Double_t RooMomentMorph::getVal(const RooArgSet* set) const
379 {
380  _curNormSet = set ? (RooArgSet*)set : (RooArgSet*)&_varList ;
381  return RooAbsPdf::getVal(set) ;
382 }
383 
384 ////////////////////////////////////////////////////////////////////////////////
385 
386 RooAbsPdf* RooMomentMorph::sumPdf(const RooArgSet* nset)
387 {
388  CacheElem* cache = getCache(nset ? nset : _curNormSet) ;
389 
390  if (cache->_tracker->hasChanged(kTRUE)) {
391  cache->calculateFractions(*this,kFALSE); // verbose turned off
392  }
393 
394  return cache->_sumPdf ;
395 }
396 
397 ////////////////////////////////////////////////////////////////////////////////
398 
399 Double_t RooMomentMorph::evaluate() const
400 {
401  CacheElem* cache = getCache(_curNormSet) ;
402 
403  if (cache->_tracker->hasChanged(kTRUE)) {
404  cache->calculateFractions(*this,kFALSE); // verbose turned off
405  }
406 
407  Double_t ret = cache->_sumPdf->getVal(_pdfList.nset());
408  return ret ;
409 }
410 
411 ////////////////////////////////////////////////////////////////////////////////
412 
413 RooRealVar* RooMomentMorph::CacheElem::frac(Int_t i )
414 {
415  return (RooRealVar*)(_frac.at(i)) ;
416 }
417 
418 ////////////////////////////////////////////////////////////////////////////////
419 
420 const RooRealVar* RooMomentMorph::CacheElem::frac(Int_t i ) const
421 {
422  return (RooRealVar*)(_frac.at(i)) ;
423 }
424 
425 ////////////////////////////////////////////////////////////////////////////////
426 
427 void RooMomentMorph::CacheElem::calculateFractions(const RooMomentMorph& self, Bool_t verbose) const
428 {
429  Int_t nPdf = self._pdfList.getSize();
430 
431  Double_t dm = self.m - (*self._mref)[0];
432 
433  // fully non-linear
434  double sumposfrac=0.;
435  for (Int_t i=0; i<nPdf; ++i) {
436  double ffrac=0.;
437  for (Int_t j=0; j<nPdf; ++j) { ffrac += (*self._M)(j,i) * (j==0?1.:TMath::Power(dm,(double)j)); }
438  if (ffrac>=0) sumposfrac+=ffrac;
439  // fractions for pdf
440  ((RooRealVar*)frac(i))->setVal(ffrac);
441  // fractions for rms and mean
442  ((RooRealVar*)frac(nPdf+i))->setVal(ffrac);
443  if (verbose) { cout << ffrac << endl; }
444  }
445 
446  // various mode settings
447  int imin = self.idxmin(self.m);
448  int imax = self.idxmax(self.m);
449  double mfrac = (self.m-(*self._mref)[imin])/((*self._mref)[imax]-(*self._mref)[imin]);
450  switch (self._setting) {
451  case NonLinear:
452  // default already set above
453  break;
454 
455  case SineLinear:
456  mfrac = TMath::Sin( TMath::PiOver2()*mfrac ); // this gives a continuous differentiable transition between grid points.
457 
458  // now fall through to Linear case
459 
460  case Linear:
461  for (Int_t i=0; i<2*nPdf; ++i)
462  ((RooRealVar*)frac(i))->setVal(0.);
463  if (imax>imin) { // m in between mmin and mmax
464  ((RooRealVar*)frac(imin))->setVal(1.-mfrac);
465  ((RooRealVar*)frac(nPdf+imin))->setVal(1.-mfrac);
466  ((RooRealVar*)frac(imax))->setVal(mfrac);
467  ((RooRealVar*)frac(nPdf+imax))->setVal(mfrac);
468  } else if (imax==imin) { // m outside mmin and mmax
469  ((RooRealVar*)frac(imin))->setVal(1.);
470  ((RooRealVar*)frac(nPdf+imin))->setVal(1.);
471  }
472  break;
473  case NonLinearLinFractions:
474  for (Int_t i=0; i<nPdf; ++i)
475  ((RooRealVar*)frac(i))->setVal(0.);
476  if (imax>imin) { // m in between mmin and mmax
477  ((RooRealVar*)frac(imin))->setVal(1.-mfrac);
478  ((RooRealVar*)frac(imax))->setVal(mfrac);
479  } else if (imax==imin) { // m outside mmin and mmax
480  ((RooRealVar*)frac(imin))->setVal(1.);
481  }
482  break;
483  case NonLinearPosFractions:
484  for (Int_t i=0; i<nPdf; ++i) {
485  if (((RooRealVar*)frac(i))->getVal()<0) ((RooRealVar*)frac(i))->setVal(0.);
486  ((RooRealVar*)frac(i))->setVal(((RooRealVar*)frac(i))->getVal()/sumposfrac);
487  }
488  break;
489  }
490 
491 }
492 
493 ////////////////////////////////////////////////////////////////////////////////
494 
495 int RooMomentMorph::idxmin(const double& mval) const
496 {
497  int imin(0);
498  Int_t nPdf = _pdfList.getSize();
499  double mmin=-DBL_MAX;
500  for (Int_t i=0; i<nPdf; ++i)
501  if ( (*_mref)[i]>mmin && (*_mref)[i]<=mval ) { mmin=(*_mref)[i]; imin=i; }
502  return imin;
503 }
504 
505 
506 ////////////////////////////////////////////////////////////////////////////////
507 
508 int RooMomentMorph::idxmax(const double& mval) const
509 {
510  int imax(0);
511  Int_t nPdf = _pdfList.getSize();
512  double mmax=DBL_MAX;
513  for (Int_t i=0; i<nPdf; ++i)
514  if ( (*_mref)[i]<mmax && (*_mref)[i]>=mval ) { mmax=(*_mref)[i]; imax=i; }
515  return imax;
516 }