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