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