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