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