ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
RooFitResult.cxx
Go to the documentation of this file.
1 /*****************************************************************************
2  * Project: RooFit *
3  * Package: RooFitCore *
4  * @(#)root/roofitcore:$Id$
5  * Authors: *
6  * WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu *
7  * DK, David Kirkby, UC Irvine, dkirkby@uci.edu *
8  * *
9  * Copyright (c) 2000-2005, Regents of the University of California *
10  * and Stanford University. All rights reserved. *
11  * *
12  * Redistribution and use in source and binary forms, *
13  * with or without modification, are permitted according to the terms *
14  * listed in LICENSE (http://roofit.sourceforge.net/license.txt) *
15  *****************************************************************************/
16 
17 //////////////////////////////////////////////////////////////////////////////
18 //
19 // RooFitResult is a container class to hold the input and output
20 // of a PDF fit to a dataset. It contains:
21 //
22 // - Values of all constant parameters
23 // - Initial and final values of floating parameters with error
24 // - Correlation matrix and global correlation coefficients
25 // - NLL and EDM at mininum
26 //
27 // No references to the fitted PDF and dataset are stored
28 //
29 
30 #include "RooFit.h"
31 #include "Riostream.h"
32 
33 #include <iomanip>
34 #include "TMinuit.h"
35 #include "TMath.h"
36 #include "TMarker.h"
37 #include "TLine.h"
38 #include "TBox.h"
39 #include "TGaxis.h"
40 #include "TMatrix.h"
41 #include "TVector.h"
42 #include "TDirectory.h"
43 #include "TClass.h"
44 #include "RooFitResult.h"
45 #include "RooArgSet.h"
46 #include "RooArgList.h"
47 #include "RooRealVar.h"
48 #include "RooPlot.h"
49 #include "RooEllipse.h"
50 #include "RooRandom.h"
51 #include "RooMsgService.h"
52 #include "TH2D.h"
53 #include "TText.h"
54 #include "TMatrixDSym.h"
55 #include "RooMultiVarGaussian.h"
56 
57 
58 
59 using namespace std;
60 
62 ;
63 
64 
65 
66 ////////////////////////////////////////////////////////////////////////////////
67 /// Constructor with name and title
68 /// coverity[UNINIT_CTOR]
69 
70 RooFitResult::RooFitResult(const char* name, const char* title) :
71  TNamed(name,title), _constPars(0), _initPars(0), _finalPars(0), _globalCorr(0), _randomPars(0), _Lt(0),
72  _CM(0), _VM(0), _GC(0)
73 {
74  if (name) appendToDir(this,kTRUE) ;
75 }
76 
77 
78 ////////////////////////////////////////////////////////////////////////////////
79 /// Copy constructor
80 
82  TNamed(other),
83  RooPrintable(other),
84  RooDirItem(other),
85  _status(other._status),
86  _covQual(other._covQual),
87  _numBadNLL(other._numBadNLL),
88  _minNLL(other._minNLL),
89  _edm(other._edm),
90  _globalCorr(0),
91  _randomPars(0),
92  _Lt(0),
93  _CM(0),
94  _VM(0),
95  _GC(0)
96 {
97  _constPars = (RooArgList*) other._constPars->snapshot() ;
98  _initPars = (RooArgList*) other._initPars->snapshot() ;
99  _finalPars = (RooArgList*) other._finalPars->snapshot() ;
100  if (other._randomPars) _randomPars = (RooArgList*) other._randomPars->snapshot() ;
101  if (other._Lt) _Lt = new TMatrix(*other._Lt);
102  if (other._VM) _VM = new TMatrixDSym(*other._VM) ;
103  if (other._CM) _CM = new TMatrixDSym(*other._CM) ;
104  if (other._GC) _GC = new TVectorD(*other._GC) ;
105 
106  if (GetName())
107  appendToDir(this, kTRUE);
108 }
109 
110 
111 
112 ////////////////////////////////////////////////////////////////////////////////
113 /// Destructor
114 
116 {
117  if (_constPars) delete _constPars ;
118  if (_initPars) delete _initPars ;
119  if (_finalPars) delete _finalPars ;
120  if (_globalCorr) delete _globalCorr;
121  if (_randomPars) delete _randomPars;
122  if (_Lt) delete _Lt;
123  if (_CM) delete _CM ;
124  if (_VM) delete _VM ;
125  if (_GC) delete _GC ;
126 
128 
129  removeFromDir(this) ;
130 }
131 
132 
133 ////////////////////////////////////////////////////////////////////////////////
134 /// Fill the list of constant parameters
135 
137 {
138  if (_constPars) delete _constPars ;
139  _constPars = (RooArgList*) list.snapshot() ;
140  TIterator* iter = _constPars->createIterator() ;
141  RooAbsArg* arg ;
142  while((arg=(RooAbsArg*)iter->Next())) {
143  RooRealVar* rrv = dynamic_cast<RooRealVar*>(arg) ;
144  if (rrv) {
145  rrv->deleteSharedProperties() ;
146  }
147  }
148  delete iter ;
149 }
150 
151 
152 
153 ////////////////////////////////////////////////////////////////////////////////
154 /// Fill the list of initial values of the floating parameters
155 
157 {
158  if (_initPars) delete _initPars ;
159  _initPars = (RooArgList*) list.snapshot() ;
160  TIterator* iter = _initPars->createIterator() ;
161  RooAbsArg* arg ;
162  while((arg=(RooAbsArg*)iter->Next())) {
163  RooRealVar* rrv = dynamic_cast<RooRealVar*>(arg) ;
164  if (rrv) {
165  rrv->deleteSharedProperties() ;
166  }
167  }
168  delete iter ;
169 }
170 
171 
172 
173 ////////////////////////////////////////////////////////////////////////////////
174 /// Fill the list of final values of the floating parameters
175 
177 {
178  if (_finalPars) delete _finalPars ;
179  _finalPars = (RooArgList*) list.snapshot() ;
180 
181  TIterator* iter = _finalPars->createIterator() ;
182  RooAbsArg* arg ;
183  while((arg=(RooAbsArg*)iter->Next())) {
184  RooRealVar* rrv = dynamic_cast<RooRealVar*>(arg) ;
185  if (rrv) {
186  rrv->deleteSharedProperties() ;
187  }
188  }
189  delete iter ;
190 }
191 
192 
193 
194 ////////////////////////////////////////////////////////////////////////////////
195 
197 {
198  if (icycle>=_statusHistory.size()) {
199  coutE(InputArguments) << "RooFitResult::statusCodeHistory(" << GetName()
200  << " ERROR request for status history slot "
201  << icycle << " exceeds history count of " << _statusHistory.size() << endl ;
202  }
203  return _statusHistory[icycle].second ;
204 }
205 
206 
207 
208 ////////////////////////////////////////////////////////////////////////////////
209 
210 const char* RooFitResult::statusLabelHistory(UInt_t icycle) const
211 {
212  if (icycle>=_statusHistory.size()) {
213  coutE(InputArguments) << "RooFitResult::statusLabelHistory(" << GetName()
214  << " ERROR request for status history slot "
215  << icycle << " exceeds history count of " << _statusHistory.size() << endl ;
216  }
217  return _statusHistory[icycle].first.c_str() ;
218 }
219 
220 
221 
222 ////////////////////////////////////////////////////////////////////////////////
223 /// Add objects to a 2D plot that represent the fit results for the
224 /// two named parameters. The input frame with the objects added is
225 /// returned, or zero in case of an error. Which objects are added
226 /// are determined by the options string which should be a concatenation
227 /// of the following (not case sensitive):
228 ///
229 /// M - a marker at the best fit result
230 /// E - an error ellipse calculated at 1-sigma using the error matrix at the minimum
231 /// 1 - the 1-sigma error bar for parameter 1
232 /// 2 - the 1-sigma error bar for parameter 2
233 /// B - the bounding box for the error ellipse
234 /// H - a line and horizontal axis for reading off the correlation coefficient
235 /// V - a line and vertical axis for reading off the correlation coefficient
236 /// A - draw axes for reading off the correlation coefficients with the H or V options
237 ///
238 /// You can change the attributes of objects in the returned RooPlot using the
239 /// various RooPlot::getAttXxx(name) member functions, e.g.
240 ///
241 /// plot->getAttLine("contour")->SetLineStyle(kDashed);
242 ///
243 /// Use plot->Print() for a list of all objects and their names (unfortunately most
244 /// of the ROOT builtin graphics objects like TLine are unnamed). Drag the left mouse
245 /// button along the labels of either axis button to interactively zoom in a plot.
246 
247 RooPlot *RooFitResult::plotOn(RooPlot *frame, const char *parName1, const char *parName2,
248  const char *options) const
249 {
250  // lookup the input parameters by name: we require that they were floated in our fit
251  const RooRealVar *par1= dynamic_cast<const RooRealVar*>(floatParsFinal().find(parName1));
252  if(0 == par1) {
253  coutE(InputArguments) << "RooFitResult::correlationPlot: parameter not floated in fit: " << parName1 << endl;
254  return 0;
255  }
256  const RooRealVar *par2= dynamic_cast<const RooRealVar*>(floatParsFinal().find(parName2));
257  if(0 == par2) {
258  coutE(InputArguments) << "RooFitResult::correlationPlot: parameter not floated in fit: " << parName2 << endl;
259  return 0;
260  }
261 
262  // options are not case sensitive
263  TString opt(options);
264  opt.ToUpper();
265 
266  // lookup the 2x2 covariance matrix elements for these variables
267  Double_t x1= par1->getVal();
268  Double_t x2= par2->getVal();
269  Double_t s1= par1->getError();
270  Double_t s2= par2->getError();
271  Double_t rho= correlation(parName1, parName2);
272 
273  // add a 1-sigma error ellipse, if requested
274  if(opt.Contains("E")) {
275  RooEllipse *contour= new RooEllipse("contour",x1,x2,s1,s2,rho);
276  contour->SetLineWidth(2) ;
277  frame->addPlotable(contour);
278  }
279 
280  // add the error bar for parameter 1, if requested
281  if(opt.Contains("1")) {
282  TLine *hline= new TLine(x1-s1,x2,x1+s1,x2);
283  hline->SetLineColor(kRed);
284  frame->addObject(hline);
285  }
286 
287  if(opt.Contains("2")) {
288  TLine *vline= new TLine(x1,x2-s2,x1,x2+s2);
289  vline->SetLineColor(kRed);
290  frame->addObject(vline);
291  }
292 
293  if(opt.Contains("B")) {
294  TBox *box= new TBox(x1-s1,x2-s2,x1+s1,x2+s2);
295  box->SetLineStyle(kDashed);
296  box->SetLineColor(kRed);
297  box->SetFillStyle(0);
298  frame->addObject(box);
299  }
300 
301  if(opt.Contains("H")) {
302  TLine *line= new TLine(x1-rho*s1,x2-s2,x1+rho*s1,x2+s2);
303  line->SetLineStyle(kDashed);
304  line->SetLineColor(kBlue);
305  line->SetLineWidth(2) ;
306  frame->addObject(line);
307  if(opt.Contains("A")) {
308  TGaxis *axis= new TGaxis(x1-s1,x2-s2,x1+s1,x2-s2,-1.,+1.,502,"-=");
309  axis->SetLineColor(kBlue);
310  frame->addObject(axis);
311  }
312  }
313 
314  if(opt.Contains("V")) {
315  TLine *line= new TLine(x1-s1,x2-rho*s2,x1+s1,x2+rho*s2);
316  line->SetLineStyle(kDashed);
317  line->SetLineColor(kBlue);
318  line->SetLineWidth(2) ;
319  frame->addObject(line);
320  if(opt.Contains("A")) {
321  TGaxis *axis= new TGaxis(x1-s1,x2-s2,x1-s1,x2+s2,-1.,+1.,502,"-=");
322  axis->SetLineColor(kBlue);
323  frame->addObject(axis);
324  }
325  }
326 
327  // add a marker at the fitted value, if requested
328  if(opt.Contains("M")) {
329  TMarker *marker= new TMarker(x1,x2,20);
330  marker->SetMarkerColor(kBlack);
331  frame->addObject(marker);
332  }
333 
334  return frame;
335 }
336 
337 
338 ////////////////////////////////////////////////////////////////////////////////
339 /// Return a list of floating parameter values that are perturbed from the final
340 /// fit values by random amounts sampled from the covariance matrix. The returned
341 /// object is overwritten with each call and belongs to the RooFitResult. Uses
342 /// the "square root method" to decompose the covariance matrix, which makes inverting
343 /// it unnecessary.
344 
346 {
347  Int_t nPar= _finalPars->getSize();
348  if(0 == _randomPars) { // first-time initialization
349  assert(0 != _finalPars);
350  // create the list of random values to fill
351  _randomPars= (RooArgList*)_finalPars->snapshot();
352  // calculate the elements of the upper-triangular matrix L that gives Lt*L = C
353  // where Lt is the transpose of L (the "square-root method")
354  TMatrix L(nPar,nPar);
355  for(Int_t iPar= 0; iPar < nPar; iPar++) {
356  // calculate the diagonal term first
357  L(iPar,iPar)= covariance(iPar,iPar);
358  for(Int_t k= 0; k < iPar; k++) {
359  Double_t tmp= L(k,iPar);
360  L(iPar,iPar)-= tmp*tmp;
361  }
362  L(iPar,iPar)= sqrt(L(iPar,iPar));
363  // then the off-diagonal terms
364  for(Int_t jPar= iPar+1; jPar < nPar; jPar++) {
365  L(iPar,jPar)= covariance(iPar,jPar);
366  for(Int_t k= 0; k < iPar; k++) {
367  L(iPar,jPar)-= L(k,iPar)*L(k,jPar);
368  }
369  L(iPar,jPar)/= L(iPar,iPar);
370  }
371  }
372  // remember Lt
374  }
375  else {
376  // reset to the final fit values
378  }
379 
380  // create a vector of unit Gaussian variables
381  TVector g(nPar);
382  for(Int_t k= 0; k < nPar; k++) g(k)= RooRandom::gaussian();
383  // multiply this vector by Lt to introduce the appropriate correlations
384  g*= (*_Lt);
385  // add the mean value offsets and store the results
387  RooRealVar *par(0);
388  Int_t index(0);
389  while((0 != (par= (RooRealVar*)iter->Next()))) {
390  par->setVal(par->getVal() + g(index++));
391  }
392  delete iter;
393 
394  return *_randomPars;
395 }
396 
397 
398 ////////////////////////////////////////////////////////////////////////////////
399 /// Return the correlation between parameters 'par1' and 'par2'
400 
401 Double_t RooFitResult::correlation(const char* parname1, const char* parname2) const
402 {
403  Int_t idx1 = _finalPars->index(parname1) ;
404  Int_t idx2 = _finalPars->index(parname2) ;
405  if (idx1<0) {
406  coutE(InputArguments) << "RooFitResult::correlation(" << GetName() << ") parameter " << parname1 << " is not a floating fit parameter" << endl ;
407  return 0 ;
408  }
409  if (idx2<0) {
410  coutE(InputArguments) << "RooFitResult::correlation(" << GetName() << ") parameter " << parname2 << " is not a floating fit parameter" << endl ;
411  return 0 ;
412  }
413  return correlation(idx1,idx2) ;
414 }
415 
416 
417 
418 ////////////////////////////////////////////////////////////////////////////////
419 /// Return the set of correlation coefficients of parameter 'par' with
420 /// all other floating parameters
421 
422 const RooArgList* RooFitResult::correlation(const char* parname) const
423 {
424  if (_globalCorr==0) {
426  }
427 
428  RooAbsArg* arg = _initPars->find(parname) ;
429  if (!arg) {
430  coutE(InputArguments) << "RooFitResult::correlation: variable " << parname << " not a floating parameter in fit" << endl ;
431  return 0 ;
432  }
433  return (RooArgList*)_corrMatrix.At(_initPars->index(arg)) ;
434 }
435 
436 
437 
438 ////////////////////////////////////////////////////////////////////////////////
439 /// Return the global correlation of the named parameter
440 
441 Double_t RooFitResult::globalCorr(const char* parname)
442 {
443  if (_globalCorr==0) {
445  }
446 
447  RooAbsArg* arg = _initPars->find(parname) ;
448  if (!arg) {
449  coutE(InputArguments) << "RooFitResult::globalCorr: variable " << parname << " not a floating parameter in fit" << endl ;
450  return 0 ;
451  }
452 
453  if (_globalCorr) {
454  return ((RooAbsReal*)_globalCorr->at(_initPars->index(arg)))->getVal() ;
455  } else {
456  return 1.0 ;
457  }
458 }
459 
460 
461 
462 ////////////////////////////////////////////////////////////////////////////////
463 /// Return the list of all global correlations
464 
466 {
467  if (_globalCorr==0) {
469  }
470 
471  return _globalCorr ;
472 }
473 
474 
475 
476 ////////////////////////////////////////////////////////////////////////////////
477 /// Return a correlation matrix element addressed with numeric indices.
478 
480 {
481  return (*_CM)(row,col) ;
482 }
483 
484 
485 ////////////////////////////////////////////////////////////////////////////////
486 /// Return the covariance matrix element addressed with numeric indices.
487 
489 {
490  return (*_VM)(row,col) ;
491 }
492 
493 
494 
495 ////////////////////////////////////////////////////////////////////////////////
496 /// Print fit result to stream 'os'. In Verbose mode, the contant parameters and
497 /// the initial and final values of the floating parameters are printed.
498 /// Standard mode only the final values of the floating parameters are printed
499 
500 void RooFitResult::printMultiline(ostream& os, Int_t /*contents*/, Bool_t verbose, TString indent) const
501 {
502 
503  os << endl
504  << indent << " RooFitResult: minimized FCN value: " << _minNLL << ", estimated distance to minimum: " << _edm << endl
505  << indent << " covariance matrix quality: " ;
506  switch(_covQual) {
507  case -1 : os << "Unknown, matrix was externally provided" ; break ;
508  case 0 : os << "Not calculated at all" ; break ;
509  case 1 : os << "Approximation only, not accurate" ; break ;
510  case 2 : os << "Full matrix, but forced positive-definite" ; break ;
511  case 3 : os << "Full, accurate covariance matrix" ; break ;
512  }
513  os << endl ;
514  os << indent << " Status : " ;
515  for (vector<pair<string,int> >::const_iterator iter = _statusHistory.begin() ; iter != _statusHistory.end() ; ++iter) {
516  os << iter->first << "=" << iter->second << " " ;
517  }
518  os << endl << endl ;;
519 
520  Int_t i ;
521  if (verbose) {
522  if (_constPars->getSize()>0) {
523  os << indent << " Constant Parameter Value " << endl
524  << indent << " -------------------- ------------" << endl ;
525 
526  for (i=0 ; i<_constPars->getSize() ; i++) {
527  os << indent << " " << setw(20) << ((RooAbsArg*)_constPars->at(i))->GetName()
528  << " " << setw(12) << Form("%12.4e",((RooRealVar*)_constPars->at(i))->getVal())
529  << endl ;
530  }
531 
532  os << endl ;
533  }
534 
535  // Has any parameter asymmetric errors?
536  Bool_t doAsymErr(kFALSE) ;
537  for (i=0 ; i<_finalPars->getSize() ; i++) {
538  if (((RooRealVar*)_finalPars->at(i))->hasAsymError()) {
539  doAsymErr=kTRUE ;
540  break ;
541  }
542  }
543 
544  if (doAsymErr) {
545  os << indent << " Floating Parameter InitialValue FinalValue (+HiError,-LoError) GblCorr." << endl
546  << indent << " -------------------- ------------ ---------------------------------- --------" << endl ;
547  } else {
548  os << indent << " Floating Parameter InitialValue FinalValue +/- Error GblCorr." << endl
549  << indent << " -------------------- ------------ -------------------------- --------" << endl ;
550  }
551 
552  for (i=0 ; i<_finalPars->getSize() ; i++) {
553  os << indent << " " << setw(20) << ((RooAbsArg*)_finalPars->at(i))->GetName() ;
554  os << indent << " " << setw(12) << Form("%12.4e",((RooRealVar*)_initPars->at(i))->getVal())
555  << indent << " " << setw(12) << Form("%12.4e",((RooRealVar*)_finalPars->at(i))->getVal()) ;
556 
557  if (((RooRealVar*)_finalPars->at(i))->hasAsymError()) {
558  os << setw(21) << Form(" (+%8.2e,-%8.2e)",((RooRealVar*)_finalPars->at(i))->getAsymErrorHi(),
559  -1*((RooRealVar*)_finalPars->at(i))->getAsymErrorLo()) ;
560  } else {
561  Double_t err = ((RooRealVar*)_finalPars->at(i))->getError() ;
562  os << (doAsymErr?" ":"") << " +/- " << setw(9) << Form("%9.2e",err) ;
563  }
564 
565  if (_globalCorr) {
566  os << " " << setw(8) << Form("%8.6f" ,((RooRealVar*)_globalCorr->at(i))->getVal()) ;
567  } else {
568  os << " <none>" ;
569  }
570 
571  os << endl ;
572  }
573 
574  } else {
575  os << indent << " Floating Parameter FinalValue +/- Error " << endl
576  << indent << " -------------------- --------------------------" << endl ;
577 
578  for (i=0 ; i<_finalPars->getSize() ; i++) {
579  Double_t err = ((RooRealVar*)_finalPars->at(i))->getError() ;
580  os << indent << " " << setw(20) << ((RooAbsArg*)_finalPars->at(i))->GetName()
581  << " " << setw(12) << Form("%12.4e",((RooRealVar*)_finalPars->at(i))->getVal())
582  << " +/- " << setw(9) << Form("%9.2e",err)
583  << endl ;
584  }
585  }
586 
587 
588  os << endl ;
589 }
590 
591 
592 ////////////////////////////////////////////////////////////////////////////////
593 /// Function called by RooMinimizer
594 
595 void RooFitResult::fillCorrMatrix(const std::vector<double>& globalCC, const TMatrixDSym& corrs, const TMatrixDSym& covs)
596 {
597  // Sanity check
598  if (globalCC.empty() || corrs.GetNoElements() < 1 || covs.GetNoElements() < 1) {
599  coutI(Minimization) << "RooFitResult::fillCorrMatrix: number of floating parameters is zero, correlation matrix not filled" << endl ;
600  return ;
601  }
602 
603  if (!_initPars) {
604  coutE(Minimization) << "RooFitResult::fillCorrMatrix: ERROR: list of initial parameters must be filled first" << endl ;
605  return ;
606  }
607 
608  // Delete eventual prevous correlation data holders
609  if (_CM) delete _CM ;
610  if (_VM) delete _VM ;
611  if (_GC) delete _GC ;
612 
613  // Build holding arrays for correlation coefficients
614  _CM = new TMatrixDSym(corrs) ;
615  _VM = new TMatrixDSym(covs) ;
616  _GC = new TVectorD(_CM->GetNcols()) ;
617  for(int i=0 ; i<_CM->GetNcols() ; i++) {
618  (*_GC)[i] = globalCC[i] ;
619  }
620  //fillLegacyCorrMatrix() ;
621 }
622 
623 
624 
625 
626 
627 ////////////////////////////////////////////////////////////////////////////////
628 /// Sanity check
629 
631 {
632  if (!_CM) return ;
633 
634  // Delete eventual prevous correlation data holders
635  if (_globalCorr) delete _globalCorr ;
637 
638  // Build holding arrays for correlation coefficients
639  _globalCorr = new RooArgList("globalCorrelations") ;
640 
641  TIterator* vIter = _initPars->createIterator() ;
642  RooAbsArg* arg ;
643  Int_t idx(0) ;
644  while((arg=(RooAbsArg*)vIter->Next())) {
645  // Create global correlation value holder
646  TString gcName("GC[") ;
647  gcName.Append(arg->GetName()) ;
648  gcName.Append("]") ;
649  TString gcTitle(arg->GetTitle()) ;
650  gcTitle.Append(" Global Correlation") ;
651  _globalCorr->addOwned(*(new RooRealVar(gcName.Data(),gcTitle.Data(),0.))) ;
652 
653  // Create array with correlation holders for this parameter
654  TString name("C[") ;
655  name.Append(arg->GetName()) ;
656  name.Append(",*]") ;
657  RooArgList* corrMatrixRow = new RooArgList(name.Data()) ;
658  _corrMatrix.Add(corrMatrixRow) ;
659  TIterator* vIter2 = _initPars->createIterator() ;
660  RooAbsArg* arg2 ;
661  while((arg2=(RooAbsArg*)vIter2->Next())) {
662 
663  TString cName("C[") ;
664  cName.Append(arg->GetName()) ;
665  cName.Append(",") ;
666  cName.Append(arg2->GetName()) ;
667  cName.Append("]") ;
668  TString cTitle("Correlation between ") ;
669  cTitle.Append(arg->GetName()) ;
670  cTitle.Append(" and ") ;
671  cTitle.Append(arg2->GetName()) ;
672  corrMatrixRow->addOwned(*(new RooRealVar(cName.Data(),cTitle.Data(),0.))) ;
673  }
674  delete vIter2 ;
675  idx++ ;
676  }
677  delete vIter ;
678 
679  TIterator *gcIter = _globalCorr->createIterator() ;
680  TIterator *parIter = _finalPars->createIterator() ;
681  RooRealVar* gcVal = 0;
682  for (unsigned int i = 0; i < (unsigned int)_CM->GetNcols() ; ++i) {
683 
684  // Find the next global correlation slot to fill, skipping fixed parameters
685  gcVal = (RooRealVar*) gcIter->Next() ;
686  gcVal->setVal((*_GC)(i)) ; // WVE FIX THIS
687 
688  // Fill a row of the correlation matrix
689  TIterator* cIter = ((RooArgList*)_corrMatrix.At(i))->createIterator() ;
690  for (unsigned int it = 0; it < (unsigned int)_CM->GetNcols() ; ++it) {
691  RooRealVar* cVal = (RooRealVar*) cIter->Next() ;
692  double value = (*_CM)(i,it) ;
693  cVal->setVal(value);
694  (*_CM)(i,it) = value;
695  }
696  delete cIter ;
697  }
698 
699  delete gcIter ;
700  delete parIter ;
701 
702 }
703 
704 
705 
706 
707 
708 ////////////////////////////////////////////////////////////////////////////////
709 /// Internal utility method to extract the correlation matrix and the
710 /// global correlation coefficients from the MINUIT memory buffer and
711 /// fill the internal arrays.
712 
714 {
715  // Sanity check
716  if (gMinuit->fNpar < 1) {
717  coutI(Minimization) << "RooFitResult::fillCorrMatrix: number of floating parameters is zero, correlation matrix not filled" << endl ;
718  return ;
719  }
720 
721  if (!_initPars) {
722  coutE(Minimization) << "RooFitResult::fillCorrMatrix: ERROR: list of initial parameters must be filled first" << endl ;
723  return ;
724  }
725 
726  // Delete eventual prevous correlation data holders
727  if (_CM) delete _CM ;
728  if (_VM) delete _VM ;
729  if (_GC) delete _GC ;
730 
731  // Build holding arrays for correlation coefficients
732  _CM = new TMatrixDSym(_initPars->getSize()) ;
733  _VM = new TMatrixDSym(_initPars->getSize()) ;
734  _GC = new TVectorD(_initPars->getSize()) ;
735 
736  // Extract correlation information for MINUIT (code taken from TMinuit::mnmatu() )
737 
738  // WVE: This code directly manipulates minuit internal workspace,
739  // if TMinuit code changes this may need updating
740  Int_t ndex, i, j, m, n, it /* nparm,id,ix */ ;
741  Int_t ndi, ndj /*, iso, isw2, isw5*/;
742  for (i = 1; i <= gMinuit->fNpar; ++i) {
743  ndi = i*(i + 1) / 2;
744  for (j = 1; j <= gMinuit->fNpar; ++j) {
745  m = TMath::Max(i,j);
746  n = TMath::Min(i,j);
747  ndex = m*(m-1) / 2 + n;
748  ndj = j*(j + 1) / 2;
749  gMinuit->fMATUvline[j-1] = gMinuit->fVhmat[ndex-1] / TMath::Sqrt(TMath::Abs(gMinuit->fVhmat[ndi-1]*gMinuit->fVhmat[ndj-1]));
750  }
751 
752  (*_GC)(i-1) = gMinuit->fGlobcc[i-1] ;
753 
754  // Fill a row of the correlation matrix
755  for (it = 1; it <= gMinuit->fNpar ; ++it) {
756  (*_CM)(i-1,it-1) = gMinuit->fMATUvline[it-1] ;
757  }
758  }
759 
760  for (int ii=0 ; ii<_finalPars->getSize() ; ii++) {
761  for (int jj=0 ; jj<_finalPars->getSize() ; jj++) {
762  (*_VM)(ii,jj) = (*_CM)(ii,jj) * ((RooRealVar*)_finalPars->at(ii))->getError() * ((RooRealVar*)_finalPars->at(jj))->getError() ;
763  }
764  }
765 }
766 
767 
768 
769 ////////////////////////////////////////////////////////////////////////////////
770 /// Return true if this fit result is identical to other within tolerance 'tol' on fitted values
771 /// and tolerance 'tolCor' on correlation coefficients
772 
773 Bool_t RooFitResult::isIdentical(const RooFitResult& other, Double_t tol, Double_t tolCorr, Bool_t /*verbose*/) const
774 {
775  Bool_t ret = kTRUE ;
776 
777  if (fabs(_minNLL-other._minNLL)>=tol) {
778  cout << "RooFitResult::isIdentical: minimized value of -log(L) is different " << _minNLL << " vs. " << other._minNLL << endl ;
779  ret = kFALSE ;
780  }
781 
782  for (Int_t i=0 ; i<_constPars->getSize() ; i++) {
783  RooAbsReal* ov = static_cast<RooAbsReal*>(other._constPars->find(_constPars->at(i)->GetName())) ;
784  if (!ov) {
785  cout << "RooFitResult::isIdentical: cannot find constant parameter " << _constPars->at(i)->GetName() << " in reference" << endl ;
786  ret = kFALSE ;
787  }
788  if (ov && fabs(static_cast<RooAbsReal*>(_constPars->at(i))->getVal()-ov->getVal())>=tol) {
789  cout << "RooFitResult::isIdentical: constant parameter " << _constPars->at(i)->GetName()
790  << " differs in value: " << static_cast<RooAbsReal*>(_constPars->at(i))->getVal() << " vs. " << ov->getVal() << endl ;
791  ret = kFALSE ;
792  }
793  }
794 
795  for (Int_t i=0 ; i<_initPars->getSize() ; i++) {
796  RooAbsReal* ov = static_cast<RooAbsReal*>(other._initPars->find(_initPars->at(i)->GetName())) ;
797  if (!ov) {
798  cout << "RooFitResult::isIdentical: cannot find initial parameter " << _initPars->at(i)->GetName() << " in reference" << endl ;
799  ret = kFALSE ;
800  }
801  if (ov && fabs(static_cast<RooAbsReal*>(_initPars->at(i))->getVal()-ov->getVal())>=tol) {
802  cout << "RooFitResult::isIdentical: initial parameter " << _initPars->at(i)->GetName()
803  << " differs in value: " << static_cast<RooAbsReal*>(_initPars->at(i))->getVal() << " vs. " << ov->getVal() << endl ;
804  ret = kFALSE ;
805  }
806  }
807 
808  for (Int_t i=0 ; i<_finalPars->getSize() ; i++) {
809  RooAbsReal* ov = static_cast<RooAbsReal*>(other._finalPars->find(_finalPars->at(i)->GetName())) ;
810  if (!ov) {
811  cout << "RooFitResult::isIdentical: cannot find final parameter " << _finalPars->at(i)->GetName() << " in reference" << endl ;
812  ret = kFALSE ;
813  }
814  if (ov && fabs(static_cast<RooAbsReal*>(_finalPars->at(i))->getVal()-ov->getVal())>=tol) {
815  cout << "RooFitResult::isIdentical: final parameter " << _finalPars->at(i)->GetName()
816  << " differs in value: " << static_cast<RooAbsReal*>(_finalPars->at(i))->getVal() << " vs. " << ov->getVal() << endl ;
817  ret = kFALSE ;
818  }
819  }
820 
821  // Only examine correlations for cases with >1 floating parameter
822  if (_finalPars->getSize()>1) {
823 
825  other.fillLegacyCorrMatrix() ;
826 
827  for (Int_t i=0 ; i<_globalCorr->getSize() ; i++) {
828  RooAbsReal* ov = static_cast<RooAbsReal*>(other._globalCorr->find(_globalCorr->at(i)->GetName())) ;
829  if (!ov) {
830  cout << "RooFitResult::isIdentical: cannot find global correlation coefficient " << _globalCorr->at(i)->GetName() << " in reference" << endl ;
831  ret = kFALSE ;
832  }
833  if (ov && fabs(static_cast<RooAbsReal*>(_globalCorr->at(i))->getVal()-ov->getVal())>=tolCorr) {
834  cout << "RooFitResult::isIdentical: global correlation coefficient " << _globalCorr->at(i)->GetName()
835  << " differs in value: " << static_cast<RooAbsReal*>(_globalCorr->at(i))->getVal() << " vs. " << ov->getVal() << endl ;
836  ret = kFALSE ;
837  }
838  }
839 
840  for (Int_t j=0 ; j<_corrMatrix.GetSize() ; j++) {
842  RooArgList* orow = (RooArgList*) other._corrMatrix.At(j) ;
843  for (Int_t i=0 ; i<row->getSize() ; i++) {
844  RooAbsReal* ov = static_cast<RooAbsReal*>(orow->find(row->at(i)->GetName())) ;
845  if (!ov) {
846  cout << "RooFitResult::isIdentical: cannot find correlation coefficient " << row->at(i)->GetName() << " in reference" << endl ;
847  ret = kFALSE ;
848  }
849  if (ov && fabs(static_cast<RooAbsReal*>(row->at(i))->getVal()-ov->getVal())>=tolCorr) {
850  cout << "RooFitResult::isIdentical: correlation coefficient " << row->at(i)->GetName()
851  << " differs in value: " << static_cast<RooAbsReal*>(row->at(i))->getVal() << " vs. " << ov->getVal() << endl ;
852  ret = kFALSE ;
853  }
854  }
855  }
856  }
857 
858  return ret ;
859 }
860 
861 
862 
863 ////////////////////////////////////////////////////////////////////////////////
864 /// Import the results of the last fit performed by gMinuit, interpreting
865 /// the fit parameters as the given varList of parameters.
866 
868 {
869  // Verify length of supplied varList
870  if (varList.getSize()>0 && varList.getSize()!=gMinuit->fNu) {
871  oocoutE((TObject*)0,InputArguments) << "RooFitResult::lastMinuitFit: ERROR: supplied variable list must be either empty " << endl
872  << " or match the number of variables of the last fit (" << gMinuit->fNu << ")" << endl ;
873  return 0 ;
874  }
875 
876  // Verify that all members of varList are of type RooRealVar
877  TIterator* iter = varList.createIterator() ;
878  RooAbsArg* arg ;
879  while((arg=(RooAbsArg*)iter->Next())) {
880  if (!dynamic_cast<RooRealVar*>(arg)) {
881  oocoutE((TObject*)0,InputArguments) << "RooFitResult::lastMinuitFit: ERROR: variable '" << arg->GetName() << "' is not of type RooRealVar" << endl ;
882  return 0 ;
883  }
884  }
885  delete iter ;
886 
887  RooFitResult* r = new RooFitResult("lastMinuitFit","Last MINUIT fit") ;
888 
889  // Extract names of fit parameters from MINUIT
890  // and construct corresponding RooRealVars
891  RooArgList constPars("constPars") ;
892  RooArgList floatPars("floatPars") ;
893 
894  Int_t i ;
895  for (i = 1; i <= gMinuit->fNu; ++i) {
896  if (gMinuit->fNvarl[i-1] < 0) continue;
897  Int_t l = gMinuit->fNiofex[i-1];
898  TString varName(gMinuit->fCpnam[i-1]) ;
899  Bool_t isConst(l==0) ;
900 
901  Double_t xlo = gMinuit->fAlim[i-1];
902  Double_t xhi = gMinuit->fBlim[i-1];
903  Double_t xerr = gMinuit->fWerr[l-1];
904  Double_t xval = gMinuit->fU[i-1] ;
905 
906  RooRealVar* var ;
907  if (varList.getSize()==0) {
908 
909  if ((xlo<xhi) && !isConst) {
910  var = new RooRealVar(varName,varName,xval,xlo,xhi) ;
911  } else {
912  var = new RooRealVar(varName,varName,xval) ;
913  }
914  var->setConstant(isConst) ;
915  } else {
916 
917  var = (RooRealVar*) varList.at(i-1)->Clone() ;
918  var->setConstant(isConst) ;
919  var->setVal(xval) ;
920  if (xlo<xhi) {
921  var->setRange(xlo,xhi) ;
922  }
923  if (varName.CompareTo(var->GetName())) {
924  oocoutI((TObject*)0,Eval) << "RooFitResult::lastMinuitFit: fit parameter '" << varName
925  << "' stored in variable '" << var->GetName() << "'" << endl ;
926  }
927 
928  }
929 
930  if (isConst) {
931  constPars.addOwned(*var) ;
932  } else {
933  var->setError(xerr) ;
934  floatPars.addOwned(*var) ;
935  }
936  }
937 
938  Int_t icode,npari,nparx ;
939  Double_t fmin,edm,errdef ;
940  gMinuit->mnstat(fmin,edm,errdef,npari,nparx,icode) ;
941 
942  r->setConstParList(constPars) ;
943  r->setInitParList(floatPars) ;
944  r->setFinalParList(floatPars) ;
945  r->setMinNLL(fmin) ;
946  r->setEDM(edm) ;
947  r->setCovQual(icode) ;
948  r->setStatus(gMinuit->fStatus) ;
949  r->fillCorrMatrix() ;
950 
951  return r ;
952 }
953 
954 
955 
956 ////////////////////////////////////////////////////////////////////////////////
957 /// Store externally provided correlation matrix in this RooFitResult ;
958 
960 {
961  // Delete any previous matrices
962  if (_VM) {
963  delete _VM ;
964  }
965  if (_CM) {
966  delete _CM ;
967  }
968 
969  // Clone input covariance matrix ;
970  _VM = (TMatrixDSym*) V.Clone() ;
971 
972  // Now construct correlation matrix from it
973  _CM = (TMatrixDSym*) _VM->Clone() ;
974  for (Int_t i=0 ; i<_CM->GetNrows() ; i++) {
975  for (Int_t j=0 ; j<_CM->GetNcols() ; j++) {
976  if (i!=j) {
977  (*_CM)(i,j) = (*_CM)(i,j) / sqrt((*_CM)(i,i)*(*_CM)(j,j)) ;
978  }
979  }
980  }
981  for (Int_t i=0 ; i<_CM->GetNrows() ; i++) {
982  (*_CM)(i,i) = 1.0 ;
983  }
984 
985  _covQual = -1 ;
986 }
987 
988 
989 
990 ////////////////////////////////////////////////////////////////////////////////
991 /// Return TH2D of correlation matrix
992 
994 {
995  Int_t n = _CM->GetNcols() ;
996 
997  TH2D* hh = new TH2D(name,name,n,0,n,n,0,n) ;
998 
999  for (Int_t i = 0 ; i<n ; i++) {
1000  for (Int_t j = 0 ; j<n; j++) {
1001  hh->Fill(i+0.5,n-j-0.5,(*_CM)(i,j)) ;
1002  }
1003  hh->GetXaxis()->SetBinLabel(i+1,_finalPars->at(i)->GetName()) ;
1004  hh->GetYaxis()->SetBinLabel(n-i,_finalPars->at(i)->GetName()) ;
1005  }
1006  hh->SetMinimum(-1) ;
1007  hh->SetMaximum(+1) ;
1008 
1009 
1010  return hh ;
1011 }
1012 
1013 
1014 
1015 
1016 ////////////////////////////////////////////////////////////////////////////////
1017 /// Return covariance matrix
1018 
1020 {
1021  return *_VM ;
1022 }
1023 
1024 
1025 
1026 
1027 ////////////////////////////////////////////////////////////////////////////////
1028 /// Return a reduced covariance matrix (Note that Vred _is_ a simple sub-matrix of V,
1029 /// row/columns are ordered to matched the convention given in input argument 'params'
1030 
1032 {
1033  const TMatrixDSym& V = covarianceMatrix() ;
1034 
1035  // Handle case where V==Vred here
1036  if (V.GetNcols()==params.getSize()) {
1037  return V ;
1038  }
1039 
1040 
1041  // Make sure that all given params were floating parameters in the represented fit
1042  RooArgList params2 ;
1043  TIterator* iter = params.createIterator() ;
1044  RooAbsArg* arg ;
1045  while((arg=(RooAbsArg*)iter->Next())) {
1046  if (_finalPars->find(arg->GetName())) {
1047  params2.add(*arg) ;
1048  } else {
1049  coutW(InputArguments) << "RooFitResult::reducedCovarianceMatrix(" << GetName() << ") WARNING input variable "
1050  << arg->GetName() << " was not a floating parameters in fit result and is ignored" << endl ;
1051  }
1052  }
1053  delete iter ;
1054 
1055  // Need to order params in vector in same order as in covariance matrix
1056  RooArgList params3 ;
1057  iter = _finalPars->createIterator() ;
1058  while((arg=(RooAbsArg*)iter->Next())) {
1059  if (params2.find(arg->GetName())) {
1060  params3.add(*arg) ;
1061  }
1062  }
1063  delete iter ;
1064 
1065  // Find (subset) of parameters that are stored in the covariance matrix
1066  vector<int> map1, map2 ;
1067  for (int i=0 ; i<_finalPars->getSize() ; i++) {
1068  if (params3.find(_finalPars->at(i)->GetName())) {
1069  map1.push_back(i) ;
1070  } else {
1071  map2.push_back(i) ;
1072  }
1073  }
1074 
1075  TMatrixDSym S11, S22 ;
1076  TMatrixD S12, S21 ;
1077  RooMultiVarGaussian::blockDecompose(V,map1,map2,S11,S12,S21,S22) ;
1078 
1079  return S11 ;
1080 }
1081 
1082 
1083 
1084 ////////////////////////////////////////////////////////////////////////////////
1085 /// Return a reduced covariance matrix, which is calculated as
1086 /// ___ -1
1087 /// Vred = V22 = V11 - V12 * V22 * V21
1088 ///
1089 /// Where V11,V12,V21,V22 represent a block decomposition of the covariance matrix into observables that
1090 /// are propagated (labeled by index '1') and that are not propagated (labeled by index '2'), and V22bar
1091 /// is the Shur complement of V22, calculated as shown above
1092 ///
1093 /// (Note that Vred is _not_ a simple sub-matrix of V)
1094 
1096 {
1097  const TMatrixDSym& V = covarianceMatrix() ;
1098 
1099  // Handle case where V==Vred here
1100  if (V.GetNcols()==params.getSize()) {
1101  return V ;
1102  }
1103 
1104  Double_t det = V.Determinant() ;
1105 
1106  if (det<=0) {
1107  coutE(Eval) << "RooFitResult::conditionalCovarianceMatrix(" << GetName() << ") ERROR: covariance matrix is not positive definite (|V|="
1108  << det << ") cannot reduce it" << endl ;
1109  throw string("RooFitResult::conditionalCovarianceMatrix() ERROR, input covariance matrix is not positive definite") ;
1110  }
1111 
1112  // Make sure that all given params were floating parameters in the represented fit
1113  RooArgList params2 ;
1114  TIterator* iter = params.createIterator() ;
1115  RooAbsArg* arg ;
1116  while((arg=(RooAbsArg*)iter->Next())) {
1117  if (_finalPars->find(arg->GetName())) {
1118  params2.add(*arg) ;
1119  } else {
1120  coutW(InputArguments) << "RooFitResult::conditionalCovarianceMatrix(" << GetName() << ") WARNING input variable "
1121  << arg->GetName() << " was not a floating parameters in fit result and is ignored" << endl ;
1122  }
1123  }
1124  delete iter ;
1125 
1126  // Need to order params in vector in same order as in covariance matrix
1127  RooArgList params3 ;
1128  iter = _finalPars->createIterator() ;
1129  while((arg=(RooAbsArg*)iter->Next())) {
1130  if (params2.find(arg->GetName())) {
1131  params3.add(*arg) ;
1132  }
1133  }
1134  delete iter ;
1135 
1136  // Find (subset) of parameters that are stored in the covariance matrix
1137  vector<int> map1, map2 ;
1138  for (int i=0 ; i<_finalPars->getSize() ; i++) {
1139  if (params3.find(_finalPars->at(i)->GetName())) {
1140  map1.push_back(i) ;
1141  } else {
1142  map2.push_back(i) ;
1143  }
1144  }
1145 
1146  // Rearrange matrix in block form with 'params' first and 'others' last
1147  // (preserving relative order)
1148  TMatrixDSym S11, S22 ;
1149  TMatrixD S12, S21 ;
1150  RooMultiVarGaussian::blockDecompose(V,map1,map2,S11,S12,S21,S22) ;
1151 
1152  // Constructed conditional matrix form -1
1153  // F(X1|X2) --> CovI --> S22bar = S11 - S12 S22 S21
1154 
1155  // Do eigenvalue decomposition
1156  TMatrixD S22Inv(TMatrixD::kInverted,S22) ;
1157  TMatrixD S22bar = S11 - S12 * (S22Inv * S21) ;
1158 
1159  // Convert explicitly to symmetric form
1160  TMatrixDSym Vred(S22bar.GetNcols()) ;
1161  for (int i=0 ; i<Vred.GetNcols() ; i++) {
1162  for (int j=i ; j<Vred.GetNcols() ; j++) {
1163  Vred(i,j) = (S22bar(i,j) + S22bar(j,i))/2 ;
1164  Vred(j,i) = Vred(i,j) ;
1165  }
1166  }
1167 
1168  return Vred ;
1169 }
1170 
1171 
1172 
1173 ////////////////////////////////////////////////////////////////////////////////
1174 /// Return correlation matrix ;
1175 
1177 {
1178  return *_CM ;
1179 }
1180 
1181 
1182 
1183 ////////////////////////////////////////////////////////////////////////////////
1184 /// Return a p.d.f that represents the fit result as a multi-variate probability densisty
1185 /// function on the floating fit parameters, including correlations
1186 
1188 {
1189  const TMatrixDSym& V = covarianceMatrix() ;
1190  Double_t det = V.Determinant() ;
1191 
1192  if (det<=0) {
1193  coutE(Eval) << "RooFitResult::createHessePdf(" << GetName() << ") ERROR: covariance matrix is not positive definite (|V|="
1194  << det << ") cannot construct p.d.f" << endl ;
1195  return 0 ;
1196  }
1197 
1198  // Make sure that all given params were floating parameters in the represented fit
1199  RooArgList params2 ;
1200  TIterator* iter = params.createIterator() ;
1201  RooAbsArg* arg ;
1202  while((arg=(RooAbsArg*)iter->Next())) {
1203  if (_finalPars->find(arg->GetName())) {
1204  params2.add(*arg) ;
1205  } else {
1206  coutW(InputArguments) << "RooFitResult::createHessePdf(" << GetName() << ") WARNING input variable "
1207  << arg->GetName() << " was not a floating parameters in fit result and is ignored" << endl ;
1208  }
1209  }
1210  delete iter ;
1211 
1212  // Need to order params in vector in same order as in covariance matrix
1213  RooArgList params3 ;
1214  iter = _finalPars->createIterator() ;
1215  while((arg=(RooAbsArg*)iter->Next())) {
1216  if (params2.find(arg->GetName())) {
1217  params3.add(*arg) ;
1218  }
1219  }
1220  delete iter ;
1221 
1222 
1223  // Handle special case of representing full covariance matrix here
1224  if (params3.getSize()==_finalPars->getSize()) {
1225 
1226  RooArgList mu ;
1227  for (Int_t i=0 ; i<_finalPars->getSize() ; i++) {
1228  RooRealVar* parclone = (RooRealVar*) _finalPars->at(i)->Clone(Form("%s_centralvalue",_finalPars->at(i)->GetName())) ;
1229  parclone->setConstant(kTRUE) ;
1230  mu.add(*parclone) ;
1231  }
1232 
1233  string name = Form("pdf_%s",GetName()) ;
1234  string title = Form("P.d.f of %s",GetTitle()) ;
1235 
1236  // Create p.d.f.
1237  RooAbsPdf* mvg = new RooMultiVarGaussian(name.c_str(),title.c_str(),params3,mu,V) ;
1238  mvg->addOwnedComponents(mu) ;
1239  return mvg ;
1240  }
1241 
1242  // -> ->
1243  // Handle case of conditional p.d.f. MVG(p1|p2) here
1244 
1245  // Find (subset) of parameters that are stored in the covariance matrix
1246  vector<int> map1, map2 ;
1247  for (int i=0 ; i<_finalPars->getSize() ; i++) {
1248  if (params3.find(_finalPars->at(i)->GetName())) {
1249  map1.push_back(i) ;
1250  } else {
1251  map2.push_back(i) ;
1252  }
1253  }
1254 
1255  // Rearrange matrix in block form with 'params' first and 'others' last
1256  // (preserving relative order)
1257  TMatrixDSym S11, S22 ;
1258  TMatrixD S12, S21 ;
1259  RooMultiVarGaussian::blockDecompose(V,map1,map2,S11,S12,S21,S22) ;
1260 
1261  // Calculate offset vectors mu1 and mu2
1262  RooArgList mu1 ;
1263  for (UInt_t i=0 ; i<map1.size() ; i++) {
1264  RooRealVar* parclone = (RooRealVar*) _finalPars->at(map1[i])->Clone(Form("%s_centralvalue",_finalPars->at(map1[i])->GetName())) ;
1265  parclone->setConstant(kTRUE) ;
1266  mu1.add(*parclone) ;
1267  }
1268 
1269  // Constructed conditional matrix form -1
1270  // F(X1|X2) --> CovI --> S22bar = S11 - S12 S22 S21
1271 
1272  // Do eigenvalue decomposition
1273  TMatrixD S22Inv(TMatrixD::kInverted,S22) ;
1274  TMatrixD S22bar = S11 - S12 * (S22Inv * S21) ;
1275 
1276  // Convert explicitly to symmetric form
1277  TMatrixDSym Vred(S22bar.GetNcols()) ;
1278  for (int i=0 ; i<Vred.GetNcols() ; i++) {
1279  for (int j=i ; j<Vred.GetNcols() ; j++) {
1280  Vred(i,j) = (S22bar(i,j) + S22bar(j,i))/2 ;
1281  Vred(j,i) = Vred(i,j) ;
1282  }
1283  }
1284  string name = Form("pdf_%s",GetName()) ;
1285  string title = Form("P.d.f of %s",GetTitle()) ;
1286 
1287  // Create p.d.f.
1288  RooAbsPdf* ret = new RooMultiVarGaussian(name.c_str(),title.c_str(),params3,mu1,Vred) ;
1289  ret->addOwnedComponents(mu1) ;
1290  return ret ;
1291 }
1292 
1293 
1294 
1295 ////////////////////////////////////////////////////////////////////////////////
1296 /// Change name of RooFitResult object
1297 
1298 void RooFitResult::SetName(const char *name)
1299 {
1300  if (_dir) _dir->GetList()->Remove(this);
1301  TNamed::SetName(name) ;
1302  if (_dir) _dir->GetList()->Add(this);
1303 }
1304 
1305 
1306 ////////////////////////////////////////////////////////////////////////////////
1307 /// Change name and title of RooFitResult object
1308 
1309 void RooFitResult::SetNameTitle(const char *name, const char* title)
1310 {
1311  if (_dir) _dir->GetList()->Remove(this);
1312  TNamed::SetNameTitle(name,title) ;
1313  if (_dir) _dir->GetList()->Add(this);
1314 }
1315 
1316 
1317 ////////////////////////////////////////////////////////////////////////////////
1318 /// Print name of fit result
1319 
1320 void RooFitResult::printName(ostream& os) const
1321 {
1322  os << GetName() ;
1323 }
1324 
1325 
1326 ////////////////////////////////////////////////////////////////////////////////
1327 /// Print title of fit result
1328 
1329 void RooFitResult::printTitle(ostream& os) const
1330 {
1331  os << GetTitle() ;
1332 }
1333 
1334 
1335 ////////////////////////////////////////////////////////////////////////////////
1336 /// Print class name of fit result
1337 
1338 void RooFitResult::printClassName(ostream& os) const
1339 {
1340  os << IsA()->GetName() ;
1341 }
1342 
1343 
1344 ////////////////////////////////////////////////////////////////////////////////
1345 /// Print arguments of fit result, i.e. the parameters of the fit
1346 
1347 void RooFitResult::printArgs(ostream& os) const
1348 {
1349  os << "[constPars=" << *_constPars << ",floatPars=" << *_finalPars << "]" ;
1350 }
1351 
1352 
1353 
1354 ////////////////////////////////////////////////////////////////////////////////
1355 /// Print the value of the fit result, i.e.g the status, minimized FCN, edm and covariance quality code
1356 
1357 void RooFitResult::printValue(ostream& os) const
1358 {
1359  os << "(status=" << _status << ",FCNmin=" << _minNLL << ",EDM=" << _edm << ",covQual=" << _covQual << ")" ;
1360 }
1361 
1362 
1363 ////////////////////////////////////////////////////////////////////////////////
1364 /// Configure default contents to be printed
1365 
1367 {
1368  return kName|kClassName|kArgs|kValue ;
1369 }
1370 
1371 
1372 ////////////////////////////////////////////////////////////////////////////////
1373 /// Configure mapping of Print() arguments to RooPrintable print styles
1374 
1376 {
1377  if (!opt || strlen(opt)==0) {
1378  return kStandard ;
1379  }
1380  return RooPrintable::defaultPrintStyle(opt) ;
1381 }
1382 
1383 
1384 ////////////////////////////////////////////////////////////////////////////////
1385 /// Stream an object of class RooFitResult.
1386 
1387 void RooFitResult::Streamer(TBuffer &R__b)
1388 {
1389  if (R__b.IsReading()) {
1390  UInt_t R__s, R__c;
1391  Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
1392  if (R__v>3) {
1393  R__b.ReadClassBuffer(RooFitResult::Class(),this,R__v,R__s,R__c);
1394  } else {
1395  // backward compatibitily streaming
1396  TNamed::Streamer(R__b);
1397  RooPrintable::Streamer(R__b);
1398  RooDirItem::Streamer(R__b);
1399  R__b >> _status;
1400  R__b >> _covQual;
1401  R__b >> _numBadNLL;
1402  R__b >> _minNLL;
1403  R__b >> _edm;
1404  R__b >> _constPars;
1405  R__b >> _initPars;
1406  R__b >> _finalPars;
1407  R__b >> _globalCorr;
1408  _corrMatrix.Streamer(R__b);
1409  R__b.CheckByteCount(R__s, R__c, RooFitResult::IsA());
1410 
1411  // Now fill new-style covariance and correlation matrix information
1412  // from legacy form
1413  _CM = new TMatrixDSym(_finalPars->getSize()) ;
1414  _VM = new TMatrixDSym(_CM->GetNcols()) ;
1415  _GC = new TVectorD(_CM->GetNcols()) ;
1416 
1417  TIterator *gcIter = _globalCorr->createIterator() ;
1418  TIterator *parIter = _finalPars->createIterator() ;
1419  RooRealVar* gcVal = 0;
1420  for (unsigned int i = 0; i < (unsigned int)_CM->GetNcols() ; ++i) {
1421 
1422  // Find the next global correlation slot to fill, skipping fixed parameters
1423  gcVal = (RooRealVar*) gcIter->Next() ;
1424  (*_GC)(i) = gcVal->getVal() ;
1425 
1426  // Fill a row of the correlation matrix
1427  TIterator* cIter = ((RooArgList*)_corrMatrix.At(i))->createIterator() ;
1428  for (unsigned int it = 0; it < (unsigned int)_CM->GetNcols() ; ++it) {
1429  RooRealVar* cVal = (RooRealVar*) cIter->Next() ;
1430  double value = cVal->getVal() ;
1431  (*_CM)(it,i) = value ;
1432  (*_CM)(i,it) = value;
1433  (*_VM)(it,i) = value*((RooRealVar*)_finalPars->at(i))->getError()*((RooRealVar*)_finalPars->at(it))->getError() ;
1434  (*_VM)(i,it) = (*_VM)(it,i) ;
1435  }
1436  delete cIter ;
1437  }
1438 
1439  delete gcIter ;
1440  delete parIter ;
1441  }
1442 
1443  } else {
1445  }
1446 }
1447 
Double_t * fBlim
Definition: TMinuit.h:77
virtual void printArgs(std::ostream &os) const
Print arguments of fit result, i.e. the parameters of the fit.
for(Int_t i=0;i< n;i++)
Definition: legend1.C:18
tuple row
Definition: mrt.py:26
virtual const char * GetTitle() const
Returns title of object.
Definition: TNamed.h:52
tuple par1
Definition: multifit.py:48
virtual void SetLineWidth(Width_t lwidth)
Definition: TAttLine.h:57
double par[1]
Definition: unuranDistr.cxx:38
virtual void mnstat(Double_t &fmin, Double_t &fedm, Double_t &errdef, Int_t &npari, Int_t &nparx, Int_t &istat)
*-*-*-*-*Returns concerning the current status of the minimization*-*-*-*-* *-* =====================...
Definition: TMinuit.cxx:7699
Double_t covariance(Int_t row, Int_t col) const
Return the covariance matrix element addressed with numeric indices.
Double_t _minNLL
Definition: RooFitResult.h:182
#define coutE(a)
Definition: RooMsgService.h:35
RooPlot * plotOn(RooPlot *frame, const RooAbsArg &par1, const RooAbsArg &par2, const char *options="ME") const
Definition: RooFitResult.h:142
Double_t * fGlobcc
Definition: TMinuit.h:81
virtual Bool_t add(const RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
virtual Int_t WriteClassBuffer(const TClass *cl, void *pointer)=0
virtual void Delete(Option_t *option="")
Remove all objects from the list AND delete all heap based objects.
Definition: TList.cxx:405
RooAbsCollection * snapshot(Bool_t deepCopy=kTRUE) const
Take a snap shot of current collection contents: An owning collection is returned containing clones o...
virtual void SetMaximum(Double_t maximum=-1111)
Definition: TH1.h:394
const RooArgList & floatParsFinal() const
Definition: RooFitResult.h:109
const RooArgList & constPars() const
Definition: RooFitResult.h:101
Bool_t IsReading() const
Definition: TBuffer.h:83
static Double_t gaussian(TRandom *generator=randomGenerator())
Return a Gaussian random variable with mean 0 and variance 1.
Definition: RooRandom.cxx:111
short Version_t
Definition: RtypesCore.h:61
static void blockDecompose(const TMatrixD &input, const std::vector< int > &map1, const std::vector< int > &map2, TMatrixDSym &S11, TMatrixD &S12, TMatrixD &S21, TMatrixDSym &S22)
Block decomposition of covI according to given maps of observables.
TLine * line
const char Option_t
Definition: RtypesCore.h:62
RooArgList L(const RooAbsArg &v1)
void addObject(TObject *obj, Option_t *drawOptions="", Bool_t invisible=kFALSE)
Add a generic object to this plot.
Definition: RooPlot.cxx:392
Double_t _edm
Definition: RooFitResult.h:183
virtual Bool_t addOwned(RooAbsArg &var, Bool_t silent=kFALSE)
Add the specified argument to list.
Definition: Rtypes.h:61
#define coutI(a)
Definition: RooMsgService.h:32
Create a Box.
Definition: TBox.h:44
#define assert(cond)
Definition: unittest.h:542
virtual void SetName(const char *name)
Change (i.e.
Definition: TNamed.cxx:128
virtual TList * GetList() const
Definition: TDirectory.h:157
Definition: Rtypes.h:60
TMatrixDSym * _VM
Definition: RooFitResult.h:195
virtual void printTitle(std::ostream &os) const
Print title of fit result.
void SetNameTitle(const char *name, const char *title)
Change name and title of RooFitResult object.
#define oocoutI(o, a)
Definition: RooMsgService.h:45
void addPlotable(RooPlotable *plotable, Option_t *drawOptions="", Bool_t invisible=kFALSE, Bool_t refreshNorm=kFALSE)
Add the specified plotable object to our plot.
Definition: RooPlot.cxx:447
RooAbsPdf * createHessePdf(const RooArgSet &params) const
Return a p.d.f that represents the fit result as a multi-variate probability densisty function on the...
void ToUpper()
Change string to upper case.
Definition: TString.cxx:1088
Buffer base class used for serializing objects.
Definition: TBuffer.h:42
Double_t * fVhmat
Definition: TMinuit.h:96
void setStatus(Int_t val)
Definition: RooFitResult.h:168
Int_t GetNoElements() const
Definition: TMatrixTBase.h:138
TMatrixDSym conditionalCovarianceMatrix(const RooArgList &params) const
Return a reduced covariance matrix, which is calculated as ___ -1 Vred = V22 = V11 - V12 * V22 * V21...
virtual void SetMinimum(Double_t minimum=-1111)
Definition: TH1.h:395
virtual Int_t CheckByteCount(UInt_t startpos, UInt_t bcnt, const TClass *clss)=0
TH2 * correlationHist(const char *name="correlation_matrix") const
Return TH2D of correlation matrix.
Basic string class.
Definition: TString.h:137
Manages Markers.
Definition: TMarker.h:40
Double_t * fMATUvline
Definition: TMinuit.h:114
RooArgList * _globalCorr
Definition: RooFitResult.h:188
Short_t Min(Short_t a, Short_t b)
Definition: TMathBase.h:170
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
virtual void SetFillStyle(Style_t fstyle)
Definition: TAttFill.h:52
Bool_t addOwnedComponents(const RooArgSet &comps)
Take ownership of the contents of 'comps'.
Definition: RooAbsArg.cxx:2282
#define coutW(a)
Definition: RooMsgService.h:34
A RooEllipse is a two-dimensional ellipse that can be used to represent an error contour.
Definition: RooEllipse.h:22
void box(Int_t pat, Double_t x1, Double_t y1, Double_t x2, Double_t y2)
Definition: fillpatterns.C:1
void setInitParList(const RooArgList &list)
Fill the list of initial values of the floating parameters.
virtual StyleOption defaultPrintStyle(Option_t *opt) const
virtual void SetNameTitle(const char *name, const char *title)
Change (i.e. set) all the TNamed parameters (name and title).
Definition: TNamed.cxx:142
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
virtual TObject * At(Int_t idx) const
Returns the object at position idx. Returns 0 if idx is out of range.
Definition: TList.cxx:311
Iterator abstract base class.
Definition: TIterator.h:32
TVectorD * _GC
Definition: RooFitResult.h:196
TObject Int_t at
void Minimization()
Definition: Minimization.C:30
void setEDM(Double_t val)
Definition: RooFitResult.h:167
TSocket * s1
Definition: hserv2.C:36
const char * Data() const
Definition: TString.h:349
double sqrt(double)
void appendToDir(TObject *obj, Bool_t forceMemoryResident=kFALSE)
Append object to directory.
Definition: RooDirItem.cxx:86
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TObject.cxx:203
static const double x2[5]
R__EXTERN TMinuit * gMinuit
Definition: TMinuit.h:279
virtual void printValue(std::ostream &os) const
Print the value of the fit result, i.e.g the status, minimized FCN, edm and covariance quality code...
RooPlotable is a 'mix-in' base class that define the standard RooFit plotting and printing methods...
Definition: RooPrintable.h:26
Int_t fNpar
Definition: TMinuit.h:48
virtual Double_t Determinant() const
The TNamed class is the base class for all named ROOT classes.
Definition: TNamed.h:33
void Class()
Definition: Class.C:29
#define oocoutE(o, a)
Definition: RooMsgService.h:48
void deleteSharedProperties()
std::map< std::string, std::string >::const_iterator iter
Definition: TAlienJob.cxx:54
TIterator * createIterator(Bool_t dir=kIterForward) const
ClassImp(RooFitResult)
virtual void SetMarkerColor(Color_t mcolor=1)
Definition: TAttMarker.h:51
TString & Append(const char *cs)
Definition: TString.h:492
Double_t * fAlim
Definition: TMinuit.h:76
Double_t getVal(const RooArgSet *set=0) const
Definition: RooAbsReal.h:64
const RooArgList & randomizePars() const
Return a list of floating parameter values that are perturbed from the final fit values by random amo...
Bool_t isIdentical(const RooFitResult &other, Double_t tol=5e-5, Double_t tolCorr=1e-4, Bool_t verbose=kTRUE) const
Return true if this fit result is identical to other within tolerance 'tol' on fitted values and tole...
void printMultiline(std::ostream &os, Int_t contents, Bool_t verbose=kFALSE, TString indent="") const
Print fit result to stream 'os'.
RooArgList * _finalPars
Definition: RooFitResult.h:186
TVectorT< Double_t > TVectorD
Definition: TVectorDfwd.h:24
Int_t * fNvarl
Definition: TMinuit.h:133
Double_t correlation(const RooAbsArg &par1, const RooAbsArg &par2) const
Definition: RooFitResult.h:116
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
RooRealVar represents a fundamental (non-derived) real valued object.
Definition: RooRealVar.h:37
void setMinNLL(Double_t val)
Definition: RooFitResult.h:166
TList _corrMatrix
List of global correlation coefficients.
Definition: RooFitResult.h:189
virtual void setVal(Double_t value)
Set value of variable to 'value'.
Definition: RooRealVar.cxx:203
virtual void SetLineColor(Color_t lcolor)
Definition: TAttLine.h:54
TMatrixDSym * _CM
triangular matrix used for generate random perturbations
Definition: RooFitResult.h:194
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
const double tol
Int_t statusCodeHistory(UInt_t icycle) const
ROOT::R::TRInterface & r
Definition: Object.C:4
Service class for 2-Dim histogram classes.
Definition: TH2.h:36
TString * fCpnam
Character to be plotted at the X,Y contour positions.
Definition: TMinuit.h:172
RooAbsArg * find(const char *name) const
Find object with given name in list.
return
Definition: TBase64.cxx:62
TPaveLabel title(3, 27.1, 15, 28.7,"ROOT Environment and Tools")
void setConstant(Bool_t value=kTRUE)
static RooFitResult * lastMinuitFit(const RooArgList &varList=RooArgList())
Import the results of the last fit performed by gMinuit, interpreting the fit parameters as the given...
Int_t fStatus
Definition: TMinuit.h:161
virtual TObject * Remove(TObject *obj)
Remove object from the list.
Definition: TList.cxx:675
Double_t * fWerr
Definition: TMinuit.h:80
TClass * IsA() const
unsigned int UInt_t
Definition: RtypesCore.h:42
TMarker * m
Definition: textangle.C:8
bool verbose
char * Form(const char *fmt,...)
virtual ~RooFitResult()
Destructor.
void setFinalParList(const RooArgList &list)
Fill the list of final values of the floating parameters.
const TMatrixDSym & covarianceMatrix() const
Return covariance matrix.
TLine * l
Definition: textangle.C:4
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
The axis painter class.
Definition: TGaxis.h:39
TAxis * GetYaxis()
Definition: TH1.h:320
RooFitResult(const char *name=0, const char *title=0)
Constructor with name and title coverity[UNINIT_CTOR].
static void indent(ostringstream &buf, int indent_level)
void removeFromDir(TObject *obj)
Remove object from directory it was added to.
Definition: RooDirItem.cxx:71
TMatrixF * _Lt
List of floating parameters with most recent random perturbation applied.
Definition: RooFitResult.h:192
void SetName(const char *name)
Change name of RooFitResult object.
TDirectory * _dir
Definition: RooDirItem.h:33
Int_t * fNiofex
Definition: TMinuit.h:134
RooArgList * _initPars
Definition: RooFitResult.h:185
A RooPlot is a plot frame and a container for graphics objects within that frame. ...
Definition: RooPlot.h:41
virtual Int_t ReadClassBuffer(const TClass *cl, void *pointer, const TClass *onfile_class=0)=0
Multivariate Gaussian p.d.f.
RooArgList * _randomPars
Correlation matrix (list of RooArgLists)
Definition: RooFitResult.h:191
void setConstParList(const RooArgList &list)
Fill the list of constant parameters.
virtual Int_t GetSize() const
Definition: TCollection.h:95
TH1F * s2
Definition: threadsh2.C:15
static const double x1[5]
RooAbsArg * at(Int_t idx) const
Definition: RooArgList.h:84
double Double_t
Definition: RtypesCore.h:55
RooAbsReal is the common abstract base class for objects that represent a real value and implements f...
Definition: RooAbsReal.h:53
TMatrixDSym reducedCovarianceMatrix(const RooArgList &params) const
Return a reduced covariance matrix (Note that Vred is a simple sub-matrix of V, row/columns are order...
RooDirItem is a utility base class for RooFit objects that are to be attached to ROOT directories...
Definition: RooDirItem.h:22
virtual TObject * Clone(const char *newname=0) const
Make a clone of an object using the Streamer facility.
Definition: RooAbsArg.h:75
TMatrixTSym< Double_t > TMatrixDSym
virtual void printName(std::ostream &os) const
Print name of fit result.
const char * statusLabelHistory(UInt_t icycle) const
Int_t GetNcols() const
Definition: TMatrixTBase.h:137
Double_t * fU
Definition: TMinuit.h:75
tuple par2
Definition: multifit.py:49
virtual StyleOption defaultPrintStyle(Option_t *opt) const
Configure mapping of Print() arguments to RooPrintable print styles.
virtual void SetBinLabel(Int_t bin, const char *label)
Set label for bin.
Definition: TAxis.cxx:793
TMatrixT< Float_t > TMatrix
Definition: TMatrix.h:26
virtual void SetLineStyle(Style_t lstyle)
Definition: TAttLine.h:56
#define name(a, b)
Definition: linkTestLib0.cpp:5
void setRange(const char *name, Double_t min, Double_t max)
Set range named 'name to [min,max].
Definition: RooRealVar.cxx:447
Mother of all ROOT objects.
Definition: TObject.h:58
else
Definition: TBase64.cxx:55
Bool_t axis
Definition: geodemo.C:37
RooAbsPdf is the abstract interface for all probability density functions The class provides hybrid a...
Definition: RooAbsPdf.h:41
RooArgList * _constPars
Definition: RooFitResult.h:184
virtual void Add(TObject *obj)
Definition: TList.h:81
Short_t Max(Short_t a, Short_t b)
Definition: TMathBase.h:202
Bool_t Contains(const char *pat, ECaseCompare cmp=kExact) const
Definition: TString.h:567
virtual TObject * Next()=0
const RooArgList * globalCorr()
Return the list of all global correlations.
void setCovQual(Int_t val)
Definition: RooFitResult.h:169
Definition: Rtypes.h:61
Int_t getSize() const
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
RooAbsArg is the common abstract base class for objects that represent a value (of arbitrary type) an...
Definition: RooAbsArg.h:66
const TMatrixDSym & correlationMatrix() const
Return correlation matrix ;.
void fillCorrMatrix()
Internal utility method to extract the correlation matrix and the global correlation coefficients fro...
virtual void printClassName(std::ostream &os) const
Print class name of fit result.
Double_t edm() const
Definition: RooFitResult.h:93
const Bool_t kTRUE
Definition: Rtypes.h:91
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:287
Int_t _numBadNLL
Definition: RooFitResult.h:181
virtual Int_t defaultPrintContents(Option_t *opt) const
Configure default contents to be printed.
Int_t fNu
Definition: TMinuit.h:137
void setCovarianceMatrix(TMatrixDSym &V)
Store externally provided correlation matrix in this RooFitResult ;.
const Int_t n
Definition: legend1.C:16
int CompareTo(const char *cs, ECaseCompare cmp=kExact) const
Compare a string to char *cs2.
Definition: TString.cxx:372
TAxis * GetXaxis()
Definition: TH1.h:319
virtual Version_t ReadVersion(UInt_t *start=0, UInt_t *bcnt=0, const TClass *cl=0)=0
void fillLegacyCorrMatrix() const
Sanity check.
std::vector< std::pair< std::string, int > > _statusHistory
Definition: RooFitResult.h:198
void setError(Double_t value)
Definition: RooRealVar.h:56
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:297
Double_t getError() const
Definition: RooRealVar.h:54
int ii
Definition: hprod.C:34