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