00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030 #include "RooFit.h"
00031 #include "Riostream.h"
00032
00033 #include <iomanip>
00034 #include "TMinuit.h"
00035 #include "TMath.h"
00036 #include "TMarker.h"
00037 #include "TLine.h"
00038 #include "TBox.h"
00039 #include "TGaxis.h"
00040 #include "TMatrix.h"
00041 #include "TVector.h"
00042 #include "TDirectory.h"
00043 #include "TClass.h"
00044 #include "RooFitResult.h"
00045 #include "RooArgSet.h"
00046 #include "RooArgList.h"
00047 #include "RooRealVar.h"
00048 #include "RooPlot.h"
00049 #include "RooEllipse.h"
00050 #include "RooRandom.h"
00051 #include "RooMsgService.h"
00052 #include "TH2D.h"
00053 #include "TText.h"
00054 #include "TMatrixDSym.h"
00055 #include "RooMultiVarGaussian.h"
00056
00057
00058
00059 ClassImp(RooFitResult)
00060 ;
00061
00062
00063
00064
00065 RooFitResult::RooFitResult(const char* name, const char* title) :
00066 TNamed(name,title), _constPars(0), _initPars(0), _finalPars(0), _globalCorr(0), _randomPars(0), _Lt(0),
00067 _CM(0), _VM(0), _GC(0)
00068 {
00069
00070
00071 if (name) appendToDir(this,kTRUE) ;
00072 }
00073
00074
00075
00076 RooFitResult::RooFitResult(const RooFitResult& other) :
00077 TNamed(other),
00078 RooPrintable(other),
00079 RooDirItem(other),
00080 _status(other._status),
00081 _covQual(other._covQual),
00082 _numBadNLL(other._numBadNLL),
00083 _minNLL(other._minNLL),
00084 _edm(other._edm),
00085 _globalCorr(0),
00086 _randomPars(0),
00087 _Lt(0),
00088 _CM(0),
00089 _VM(0),
00090 _GC(0)
00091 {
00092
00093
00094 _constPars = (RooArgList*) other._constPars->snapshot() ;
00095 _initPars = (RooArgList*) other._initPars->snapshot() ;
00096 _finalPars = (RooArgList*) other._finalPars->snapshot() ;
00097 if (other._randomPars) _randomPars = (RooArgList*) other._randomPars->snapshot() ;
00098 if (other._Lt) _Lt = new TMatrix(*other._Lt);
00099 if (other._VM) _VM = new TMatrixDSym(*other._VM) ;
00100 if (other._CM) _CM = new TMatrixDSym(*other._CM) ;
00101 if (other._GC) _GC = new TVectorD(*other._GC) ;
00102 }
00103
00104
00105
00106
00107 RooFitResult::~RooFitResult()
00108 {
00109
00110
00111 if (_constPars) delete _constPars ;
00112 if (_initPars) delete _initPars ;
00113 if (_finalPars) delete _finalPars ;
00114 if (_globalCorr) delete _globalCorr;
00115 if (_randomPars) delete _randomPars;
00116 if (_Lt) delete _Lt;
00117 if (_CM) delete _CM ;
00118 if (_VM) delete _VM ;
00119 if (_GC) delete _GC ;
00120
00121 _corrMatrix.Delete();
00122
00123 removeFromDir(this) ;
00124 }
00125
00126
00127
00128 void RooFitResult::setConstParList(const RooArgList& list)
00129 {
00130
00131
00132 if (_constPars) delete _constPars ;
00133 _constPars = (RooArgList*) list.snapshot() ;
00134 TIterator* iter = _constPars->createIterator() ;
00135 RooAbsArg* arg ;
00136 while((arg=(RooAbsArg*)iter->Next())) {
00137 RooRealVar* rrv = dynamic_cast<RooRealVar*>(arg) ;
00138 if (rrv) {
00139 rrv->deleteSharedProperties() ;
00140 }
00141 }
00142 delete iter ;
00143 }
00144
00145
00146
00147
00148 void RooFitResult::setInitParList(const RooArgList& list)
00149 {
00150
00151
00152 if (_initPars) delete _initPars ;
00153 _initPars = (RooArgList*) list.snapshot() ;
00154 TIterator* iter = _initPars->createIterator() ;
00155 RooAbsArg* arg ;
00156 while((arg=(RooAbsArg*)iter->Next())) {
00157 RooRealVar* rrv = dynamic_cast<RooRealVar*>(arg) ;
00158 if (rrv) {
00159 rrv->deleteSharedProperties() ;
00160 }
00161 }
00162 delete iter ;
00163 }
00164
00165
00166
00167
00168 void RooFitResult::setFinalParList(const RooArgList& list)
00169 {
00170
00171
00172 if (_finalPars) delete _finalPars ;
00173 _finalPars = (RooArgList*) list.snapshot() ;
00174
00175 TIterator* iter = _finalPars->createIterator() ;
00176 RooAbsArg* arg ;
00177 while((arg=(RooAbsArg*)iter->Next())) {
00178 RooRealVar* rrv = dynamic_cast<RooRealVar*>(arg) ;
00179 if (rrv) {
00180 rrv->deleteSharedProperties() ;
00181 }
00182 }
00183 delete iter ;
00184 }
00185
00186
00187
00188 RooPlot *RooFitResult::plotOn(RooPlot *frame, const char *parName1, const char *parName2,
00189 const char *options) const
00190 {
00191
00192
00193
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205
00206
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216 const RooRealVar *par1= dynamic_cast<const RooRealVar*>(floatParsFinal().find(parName1));
00217 if(0 == par1) {
00218 coutE(InputArguments) << "RooFitResult::correlationPlot: parameter not floated in fit: " << parName1 << endl;
00219 return 0;
00220 }
00221 const RooRealVar *par2= dynamic_cast<const RooRealVar*>(floatParsFinal().find(parName2));
00222 if(0 == par2) {
00223 coutE(InputArguments) << "RooFitResult::correlationPlot: parameter not floated in fit: " << parName2 << endl;
00224 return 0;
00225 }
00226
00227
00228 TString opt(options);
00229 opt.ToUpper();
00230
00231
00232 Double_t x1= par1->getVal();
00233 Double_t x2= par2->getVal();
00234 Double_t s1= par1->getError();
00235 Double_t s2= par2->getError();
00236 Double_t rho= correlation(parName1, parName2);
00237
00238
00239 if(opt.Contains("E")) {
00240 RooEllipse *contour= new RooEllipse("contour",x1,x2,s1,s2,rho);
00241 contour->SetLineWidth(2) ;
00242 frame->addPlotable(contour);
00243 }
00244
00245
00246 if(opt.Contains("1")) {
00247 TLine *hline= new TLine(x1-s1,x2,x1+s1,x2);
00248 hline->SetLineColor(kRed);
00249 frame->addObject(hline);
00250 }
00251
00252 if(opt.Contains("2")) {
00253 TLine *vline= new TLine(x1,x2-s2,x1,x2+s2);
00254 vline->SetLineColor(kRed);
00255 frame->addObject(vline);
00256 }
00257
00258 if(opt.Contains("B")) {
00259 TBox *box= new TBox(x1-s1,x2-s2,x1+s1,x2+s2);
00260 box->SetLineStyle(kDashed);
00261 box->SetLineColor(kRed);
00262 box->SetFillStyle(0);
00263 frame->addObject(box);
00264 }
00265
00266 if(opt.Contains("H")) {
00267 TLine *line= new TLine(x1-rho*s1,x2-s2,x1+rho*s1,x2+s2);
00268 line->SetLineStyle(kDashed);
00269 line->SetLineColor(kBlue);
00270 line->SetLineWidth(2) ;
00271 frame->addObject(line);
00272 if(opt.Contains("A")) {
00273 TGaxis *axis= new TGaxis(x1-s1,x2-s2,x1+s1,x2-s2,-1.,+1.,502,"-=");
00274 axis->SetLineColor(kBlue);
00275 frame->addObject(axis);
00276 }
00277 }
00278
00279 if(opt.Contains("V")) {
00280 TLine *line= new TLine(x1-s1,x2-rho*s2,x1+s1,x2+rho*s2);
00281 line->SetLineStyle(kDashed);
00282 line->SetLineColor(kBlue);
00283 line->SetLineWidth(2) ;
00284 frame->addObject(line);
00285 if(opt.Contains("A")) {
00286 TGaxis *axis= new TGaxis(x1-s1,x2-s2,x1-s1,x2+s2,-1.,+1.,502,"-=");
00287 axis->SetLineColor(kBlue);
00288 frame->addObject(axis);
00289 }
00290 }
00291
00292
00293 if(opt.Contains("M")) {
00294 TMarker *marker= new TMarker(x1,x2,20);
00295 marker->SetMarkerColor(kBlack);
00296 frame->addObject(marker);
00297 }
00298
00299 return frame;
00300 }
00301
00302
00303
00304 const RooArgList& RooFitResult::randomizePars() const
00305 {
00306
00307
00308
00309
00310
00311
00312 Int_t nPar= _finalPars->getSize();
00313 if(0 == _randomPars) {
00314 assert(0 != _finalPars);
00315
00316 _randomPars= (RooArgList*)_finalPars->snapshot();
00317
00318
00319 TMatrix L(nPar,nPar);
00320 for(Int_t iPar= 0; iPar < nPar; iPar++) {
00321
00322 L(iPar,iPar)= covariance(iPar,iPar);
00323 for(Int_t k= 0; k < iPar; k++) {
00324 Double_t tmp= L(k,iPar);
00325 L(iPar,iPar)-= tmp*tmp;
00326 }
00327 L(iPar,iPar)= sqrt(L(iPar,iPar));
00328
00329 for(Int_t jPar= iPar+1; jPar < nPar; jPar++) {
00330 L(iPar,jPar)= covariance(iPar,jPar);
00331 for(Int_t k= 0; k < iPar; k++) {
00332 L(iPar,jPar)-= L(k,iPar)*L(k,jPar);
00333 }
00334 L(iPar,jPar)/= L(iPar,iPar);
00335 }
00336 }
00337
00338 _Lt= new TMatrix(TMatrix::kTransposed,L);
00339 }
00340 else {
00341
00342 *_randomPars= *_finalPars;
00343 }
00344
00345
00346 TVector g(nPar);
00347 for(Int_t k= 0; k < nPar; k++) g(k)= RooRandom::gaussian();
00348
00349 g*= (*_Lt);
00350
00351 TIterator *iter= _randomPars->createIterator();
00352 RooRealVar *par(0);
00353 Int_t index(0);
00354 while((0 != (par= (RooRealVar*)iter->Next()))) {
00355 par->setVal(par->getVal() + g(index++));
00356 }
00357 delete iter;
00358
00359 return *_randomPars;
00360 }
00361
00362
00363
00364 Double_t RooFitResult::correlation(const char* parname1, const char* parname2) const
00365 {
00366
00367 Int_t idx1 = _finalPars->index(parname1) ;
00368 Int_t idx2 = _finalPars->index(parname2) ;
00369 if (idx1<0) {
00370 coutE(InputArguments) << "RooFitResult::correlation(" << GetName() << ") parameter " << parname1 << " is not a floating fit parameter" << endl ;
00371 return 0 ;
00372 }
00373 if (idx2<0) {
00374 coutE(InputArguments) << "RooFitResult::correlation(" << GetName() << ") parameter " << parname2 << " is not a floating fit parameter" << endl ;
00375 return 0 ;
00376 }
00377 return correlation(idx1,idx2) ;
00378 }
00379
00380
00381
00382
00383 const RooArgList* RooFitResult::correlation(const char* parname) const
00384 {
00385
00386
00387
00388 if (_globalCorr==0) {
00389 fillLegacyCorrMatrix() ;
00390 }
00391
00392 RooAbsArg* arg = _initPars->find(parname) ;
00393 if (!arg) {
00394 coutE(InputArguments) << "RooFitResult::correlation: variable " << parname << " not a floating parameter in fit" << endl ;
00395 return 0 ;
00396 }
00397 return (RooArgList*)_corrMatrix.At(_initPars->index(arg)) ;
00398 }
00399
00400
00401
00402
00403 Double_t RooFitResult::globalCorr(const char* parname)
00404 {
00405
00406
00407 if (_globalCorr==0) {
00408 fillLegacyCorrMatrix() ;
00409 }
00410
00411 RooAbsArg* arg = _initPars->find(parname) ;
00412 if (!arg) {
00413 coutE(InputArguments) << "RooFitResult::globalCorr: variable " << parname << " not a floating parameter in fit" << endl ;
00414 return 0 ;
00415 }
00416
00417 if (_globalCorr) {
00418 return ((RooAbsReal*)_globalCorr->at(_initPars->index(arg)))->getVal() ;
00419 } else {
00420 return 1.0 ;
00421 }
00422 }
00423
00424
00425
00426
00427 const RooArgList* RooFitResult::globalCorr()
00428 {
00429
00430
00431 if (_globalCorr==0) {
00432 fillLegacyCorrMatrix() ;
00433 }
00434
00435 return _globalCorr ;
00436 }
00437
00438
00439
00440
00441 Double_t RooFitResult::correlation(Int_t row, Int_t col) const
00442 {
00443
00444 return (*_CM)(row,col) ;
00445 }
00446
00447
00448
00449 Double_t RooFitResult::covariance(Int_t row, Int_t col) const
00450 {
00451
00452 return (*_VM)(row,col) ;
00453 }
00454
00455
00456
00457
00458 void RooFitResult::printMultiline(ostream& os, Int_t , Bool_t verbose, TString indent) const
00459 {
00460
00461
00462
00463
00464
00465 os << endl
00466 << indent << " RooFitResult: minimized FCN value: " << _minNLL << ", estimated distance to minimum: " << _edm << endl
00467 << indent << " covariance matrix quality: " ;
00468 switch(_covQual) {
00469 case -1 : os << "Unknown, matrix was externally provided" ; break ;
00470 case 0 : os << "Not calculated at all" ; break ;
00471 case 1 : os << "Approximation only, not accurate" ; break ;
00472 case 2 : os << "Full matrix, but forced positive-definite" ; break ;
00473 case 3 : os << "Full, accurate covariance matrix" ; break ;
00474 }
00475 os << endl
00476 << endl ;
00477
00478 Int_t i ;
00479 if (verbose) {
00480 if (_constPars->getSize()>0) {
00481 os << indent << " Constant Parameter Value " << endl
00482 << indent << " -------------------- ------------" << endl ;
00483
00484 for (i=0 ; i<_constPars->getSize() ; i++) {
00485 os << indent << " " << setw(20) << ((RooAbsArg*)_constPars->at(i))->GetName()
00486 << " " << setw(12) << Form("%12.4e",((RooRealVar*)_constPars->at(i))->getVal())
00487 << endl ;
00488 }
00489
00490 os << endl ;
00491 }
00492
00493
00494 Bool_t doAsymErr(kFALSE) ;
00495 for (i=0 ; i<_finalPars->getSize() ; i++) {
00496 if (((RooRealVar*)_finalPars->at(i))->hasAsymError()) {
00497 doAsymErr=kTRUE ;
00498 break ;
00499 }
00500 }
00501
00502 if (doAsymErr) {
00503 os << indent << " Floating Parameter InitialValue FinalValue (+HiError,-LoError) GblCorr." << endl
00504 << indent << " -------------------- ------------ ---------------------------------- --------" << endl ;
00505 } else {
00506 os << indent << " Floating Parameter InitialValue FinalValue +/- Error GblCorr." << endl
00507 << indent << " -------------------- ------------ -------------------------- --------" << endl ;
00508 }
00509
00510 for (i=0 ; i<_finalPars->getSize() ; i++) {
00511 os << indent << " " << setw(20) << ((RooAbsArg*)_finalPars->at(i))->GetName() ;
00512 os << indent << " " << setw(12) << Form("%12.4e",((RooRealVar*)_initPars->at(i))->getVal())
00513 << indent << " " << setw(12) << Form("%12.4e",((RooRealVar*)_finalPars->at(i))->getVal()) ;
00514
00515 if (((RooRealVar*)_finalPars->at(i))->hasAsymError()) {
00516 os << setw(21) << Form(" (+%8.2e,-%8.2e)",((RooRealVar*)_finalPars->at(i))->getAsymErrorHi(),
00517 -1*((RooRealVar*)_finalPars->at(i))->getAsymErrorLo()) ;
00518 } else {
00519 Double_t err = ((RooRealVar*)_finalPars->at(i))->getError() ;
00520 os << (doAsymErr?" ":"") << " +/- " << setw(9) << Form("%9.2e",err) ;
00521 }
00522
00523 if (_globalCorr) {
00524 os << " " << setw(8) << Form("%8.6f" ,((RooRealVar*)_globalCorr->at(i))->getVal()) ;
00525 } else {
00526 os << " <none>" ;
00527 }
00528
00529 os << endl ;
00530 }
00531
00532 } else {
00533 os << indent << " Floating Parameter FinalValue +/- Error " << endl
00534 << indent << " -------------------- --------------------------" << endl ;
00535
00536 for (i=0 ; i<_finalPars->getSize() ; i++) {
00537 Double_t err = ((RooRealVar*)_finalPars->at(i))->getError() ;
00538 os << indent << " " << setw(20) << ((RooAbsArg*)_finalPars->at(i))->GetName()
00539 << " " << setw(12) << Form("%12.4e",((RooRealVar*)_finalPars->at(i))->getVal())
00540 << " +/- " << setw(9) << Form("%9.2e",err)
00541 << endl ;
00542 }
00543 }
00544
00545
00546 os << endl ;
00547 }
00548
00549
00550
00551 void RooFitResult::fillCorrMatrix(const std::vector<double>& globalCC, const TMatrixDSym& corrs, const TMatrixDSym& covs)
00552 {
00553
00554
00555
00556 if (globalCC.empty() || corrs.GetNoElements() < 1 || covs.GetNoElements() < 1) {
00557 coutI(Minimization) << "RooFitResult::fillCorrMatrix: number of floating parameters is zero, correlation matrix not filled" << endl ;
00558 return ;
00559 }
00560
00561 if (!_initPars) {
00562 coutE(Minimization) << "RooFitResult::fillCorrMatrix: ERROR: list of initial parameters must be filled first" << endl ;
00563 return ;
00564 }
00565
00566
00567 if (_CM) delete _CM ;
00568 if (_VM) delete _VM ;
00569 if (_GC) delete _GC ;
00570
00571
00572 _CM = new TMatrixDSym(corrs) ;
00573 _VM = new TMatrixDSym(covs) ;
00574 _GC = new TVectorD(_CM->GetNcols()) ;
00575 for(int i=0 ; i<_CM->GetNcols() ; i++) {
00576 (*_GC)[i] = globalCC[i] ;
00577 }
00578 fillLegacyCorrMatrix() ;
00579 }
00580
00581
00582
00583
00584
00585
00586 void RooFitResult::fillLegacyCorrMatrix() const
00587 {
00588
00589 if (!_CM) return ;
00590
00591
00592 if (_globalCorr) delete _globalCorr ;
00593 _corrMatrix.Delete();
00594
00595
00596 _globalCorr = new RooArgList("globalCorrelations") ;
00597
00598 TIterator* vIter = _initPars->createIterator() ;
00599 RooAbsArg* arg ;
00600 Int_t idx(0) ;
00601 while((arg=(RooAbsArg*)vIter->Next())) {
00602
00603 TString gcName("GC[") ;
00604 gcName.Append(arg->GetName()) ;
00605 gcName.Append("]") ;
00606 TString gcTitle(arg->GetTitle()) ;
00607 gcTitle.Append(" Global Correlation") ;
00608 _globalCorr->addOwned(*(new RooRealVar(gcName.Data(),gcTitle.Data(),0.))) ;
00609
00610
00611 TString name("C[") ;
00612 name.Append(arg->GetName()) ;
00613 name.Append(",*]") ;
00614 RooArgList* corrMatrixRow = new RooArgList(name.Data()) ;
00615 _corrMatrix.Add(corrMatrixRow) ;
00616 TIterator* vIter2 = _initPars->createIterator() ;
00617 RooAbsArg* arg2 ;
00618 while((arg2=(RooAbsArg*)vIter2->Next())) {
00619
00620 TString cName("C[") ;
00621 cName.Append(arg->GetName()) ;
00622 cName.Append(",") ;
00623 cName.Append(arg2->GetName()) ;
00624 cName.Append("]") ;
00625 TString cTitle("Correlation between ") ;
00626 cTitle.Append(arg->GetName()) ;
00627 cTitle.Append(" and ") ;
00628 cTitle.Append(arg2->GetName()) ;
00629 corrMatrixRow->addOwned(*(new RooRealVar(cName.Data(),cTitle.Data(),0.))) ;
00630 }
00631 delete vIter2 ;
00632 idx++ ;
00633 }
00634 delete vIter ;
00635
00636 TIterator *gcIter = _globalCorr->createIterator() ;
00637 TIterator *parIter = _finalPars->createIterator() ;
00638 RooRealVar* gcVal = 0;
00639 for (unsigned int i = 0; i < (unsigned int)_CM->GetNcols() ; ++i) {
00640
00641
00642 gcVal = (RooRealVar*) gcIter->Next() ;
00643 gcVal->setVal((*_GC)(i)) ;
00644
00645
00646 TIterator* cIter = ((RooArgList*)_corrMatrix.At(i))->createIterator() ;
00647 for (unsigned int it = 0; it < (unsigned int)_CM->GetNcols() ; ++it) {
00648 RooRealVar* cVal = (RooRealVar*) cIter->Next() ;
00649 double value = (*_CM)(i,it) ;
00650 cVal->setVal(value);
00651 (*_CM)(i,it) = value;
00652 }
00653 delete cIter ;
00654 }
00655
00656 delete gcIter ;
00657 delete parIter ;
00658
00659 }
00660
00661
00662
00663
00664
00665
00666 void RooFitResult::fillCorrMatrix()
00667 {
00668
00669
00670
00671
00672
00673 if (gMinuit->fNpar < 1) {
00674 coutI(Minimization) << "RooFitResult::fillCorrMatrix: number of floating parameters is zero, correlation matrix not filled" << endl ;
00675 return ;
00676 }
00677
00678 if (!_initPars) {
00679 coutE(Minimization) << "RooFitResult::fillCorrMatrix: ERROR: list of initial parameters must be filled first" << endl ;
00680 return ;
00681 }
00682
00683
00684 if (_CM) delete _CM ;
00685 if (_VM) delete _VM ;
00686 if (_GC) delete _GC ;
00687
00688
00689 _CM = new TMatrixDSym(_initPars->getSize()) ;
00690 _VM = new TMatrixDSym(_initPars->getSize()) ;
00691 _GC = new TVectorD(_initPars->getSize()) ;
00692
00693
00694
00695
00696
00697 Int_t ndex, i, j, m, n, ncoef, nparm, it, ix ;
00698 Int_t ndi, ndj ;
00699 ncoef = (gMinuit->fNpagwd - 19) / 6;
00700 nparm = TMath::Min(gMinuit->fNpar,ncoef);
00701 Double_t tmp[1000] ;
00702 for (i = 1; i <= gMinuit->fNpar; ++i) {
00703 ix = gMinuit->fNexofi[i-1];
00704 ndi = i*(i + 1) / 2;
00705 for (j = 1; j <= gMinuit->fNpar; ++j) {
00706 m = TMath::Max(i,j);
00707 n = TMath::Min(i,j);
00708 ndex = m*(m-1) / 2 + n;
00709 ndj = j*(j + 1) / 2;
00710 gMinuit->fMATUvline[j-1] = gMinuit->fVhmat[ndex-1] / TMath::Sqrt(TMath::Abs(gMinuit->fVhmat[ndi-1]*gMinuit->fVhmat[ndj-1]));
00711 tmp[j-1] = gMinuit->fVhmat[ndex-1] ;
00712 }
00713 nparm = TMath::Min(gMinuit->fNpar,ncoef);
00714
00715 (*_GC)(i-1) = gMinuit->fGlobcc[i-1] ;
00716
00717
00718 for (it = 1; it <= gMinuit->fNpar ; ++it) {
00719 (*_CM)(i-1,it-1) = gMinuit->fMATUvline[it-1] ;
00720 }
00721 }
00722
00723 for (int ii=0 ; ii<_finalPars->getSize() ; ii++) {
00724 for (int jj=0 ; jj<_finalPars->getSize() ; jj++) {
00725 (*_VM)(ii,jj) = (*_CM)(ii,jj) * ((RooRealVar*)_finalPars->at(ii))->getError() * ((RooRealVar*)_finalPars->at(jj))->getError() ;
00726 }
00727 }
00728 }
00729
00730
00731
00732
00733 Bool_t RooFitResult::isIdentical(const RooFitResult& other, Double_t tol, Double_t tolCorr, Bool_t ) const
00734 {
00735
00736
00737
00738 Bool_t ret = kTRUE ;
00739
00740 if (fabs(_minNLL-other._minNLL)>=tol) {
00741 cout << "RooFitResult::isIdentical: minimized value of -log(L) is different " << _minNLL << " vs. " << other._minNLL << endl ;
00742 ret = kFALSE ;
00743 }
00744
00745 for (Int_t i=0 ; i<_constPars->getSize() ; i++) {
00746 RooAbsReal* ov = static_cast<RooAbsReal*>(other._constPars->find(_constPars->at(i)->GetName())) ;
00747 if (!ov) {
00748 cout << "RooFitResult::isIdentical: cannot find constant parameter " << _constPars->at(i)->GetName() << " in reference" << endl ;
00749 ret = kFALSE ;
00750 }
00751 if (ov && fabs(static_cast<RooAbsReal*>(_constPars->at(i))->getVal()-ov->getVal())>=tol) {
00752 cout << "RooFitResult::isIdentical: constant parameter " << _constPars->at(i)->GetName()
00753 << " differs in value: " << static_cast<RooAbsReal*>(_constPars->at(i))->getVal() << " vs. " << ov->getVal() << endl ;
00754 ret = kFALSE ;
00755 }
00756 }
00757
00758 for (Int_t i=0 ; i<_initPars->getSize() ; i++) {
00759 RooAbsReal* ov = static_cast<RooAbsReal*>(other._initPars->find(_initPars->at(i)->GetName())) ;
00760 if (!ov) {
00761 cout << "RooFitResult::isIdentical: cannot find initial parameter " << _initPars->at(i)->GetName() << " in reference" << endl ;
00762 ret = kFALSE ;
00763 }
00764 if (ov && fabs(static_cast<RooAbsReal*>(_initPars->at(i))->getVal()-ov->getVal())>=tol) {
00765 cout << "RooFitResult::isIdentical: initial parameter " << _initPars->at(i)->GetName()
00766 << " differs in value: " << static_cast<RooAbsReal*>(_initPars->at(i))->getVal() << " vs. " << ov->getVal() << endl ;
00767 ret = kFALSE ;
00768 }
00769 }
00770
00771 for (Int_t i=0 ; i<_finalPars->getSize() ; i++) {
00772 RooAbsReal* ov = static_cast<RooAbsReal*>(other._finalPars->find(_finalPars->at(i)->GetName())) ;
00773 if (!ov) {
00774 cout << "RooFitResult::isIdentical: cannot find final parameter " << _finalPars->at(i)->GetName() << " in reference" << endl ;
00775 ret = kFALSE ;
00776 }
00777 if (ov && fabs(static_cast<RooAbsReal*>(_finalPars->at(i))->getVal()-ov->getVal())>=tol) {
00778 cout << "RooFitResult::isIdentical: final parameter " << _finalPars->at(i)->GetName()
00779 << " differs in value: " << static_cast<RooAbsReal*>(_finalPars->at(i))->getVal() << " vs. " << ov->getVal() << endl ;
00780 ret = kFALSE ;
00781 }
00782 }
00783
00784
00785 if (_finalPars->getSize()>1) {
00786
00787 fillLegacyCorrMatrix() ;
00788 other.fillLegacyCorrMatrix() ;
00789
00790 for (Int_t i=0 ; i<_globalCorr->getSize() ; i++) {
00791 RooAbsReal* ov = static_cast<RooAbsReal*>(other._globalCorr->find(_globalCorr->at(i)->GetName())) ;
00792 if (!ov) {
00793 cout << "RooFitResult::isIdentical: cannot find global correlation coefficient " << _globalCorr->at(i)->GetName() << " in reference" << endl ;
00794 ret = kFALSE ;
00795 }
00796 if (ov && fabs(static_cast<RooAbsReal*>(_globalCorr->at(i))->getVal()-ov->getVal())>=tolCorr) {
00797 cout << "RooFitResult::isIdentical: global correlation coefficient " << _globalCorr->at(i)->GetName()
00798 << " differs in value: " << static_cast<RooAbsReal*>(_globalCorr->at(i))->getVal() << " vs. " << ov->getVal() << endl ;
00799 ret = kFALSE ;
00800 }
00801 }
00802
00803 for (Int_t j=0 ; j<_corrMatrix.GetSize() ; j++) {
00804 RooArgList* row = (RooArgList*) _corrMatrix.At(j) ;
00805 RooArgList* orow = (RooArgList*) other._corrMatrix.At(j) ;
00806 for (Int_t i=0 ; i<row->getSize() ; i++) {
00807 RooAbsReal* ov = static_cast<RooAbsReal*>(orow->find(row->at(i)->GetName())) ;
00808 if (!ov) {
00809 cout << "RooFitResult::isIdentical: cannot find correlation coefficient " << row->at(i)->GetName() << " in reference" << endl ;
00810 ret = kFALSE ;
00811 }
00812 if (ov && fabs(static_cast<RooAbsReal*>(row->at(i))->getVal()-ov->getVal())>=tolCorr) {
00813 cout << "RooFitResult::isIdentical: correlation coefficient " << row->at(i)->GetName()
00814 << " differs in value: " << static_cast<RooAbsReal*>(row->at(i))->getVal() << " vs. " << ov->getVal() << endl ;
00815 ret = kFALSE ;
00816 }
00817 }
00818 }
00819 }
00820
00821 return ret ;
00822 }
00823
00824
00825
00826
00827 RooFitResult* RooFitResult::lastMinuitFit(const RooArgList& varList)
00828 {
00829
00830
00831
00832
00833 if (varList.getSize()>0 && varList.getSize()!=gMinuit->fNu) {
00834 oocoutE((TObject*)0,InputArguments) << "RooFitResult::lastMinuitFit: ERROR: supplied variable list must be either empty " << endl
00835 << " or match the number of variables of the last fit (" << gMinuit->fNu << ")" << endl ;
00836 return 0 ;
00837 }
00838
00839
00840 TIterator* iter = varList.createIterator() ;
00841 RooAbsArg* arg ;
00842 while((arg=(RooAbsArg*)iter->Next())) {
00843 if (!dynamic_cast<RooRealVar*>(arg)) {
00844 oocoutE((TObject*)0,InputArguments) << "RooFitResult::lastMinuitFit: ERROR: variable '" << arg->GetName() << "' is not of type RooRealVar" << endl ;
00845 return 0 ;
00846 }
00847 }
00848 delete iter ;
00849
00850 RooFitResult* r = new RooFitResult("lastMinuitFit","Last MINUIT fit") ;
00851
00852
00853
00854 RooArgList constPars("constPars") ;
00855 RooArgList floatPars("floatPars") ;
00856
00857 Int_t i ;
00858 for (i = 1; i <= gMinuit->fNu; ++i) {
00859 if (gMinuit->fNvarl[i-1] < 0) continue;
00860 Int_t l = gMinuit->fNiofex[i-1];
00861 TString varName(gMinuit->fCpnam[i-1]) ;
00862 Bool_t isConst(l==0) ;
00863
00864 Double_t xlo = gMinuit->fAlim[i-1];
00865 Double_t xhi = gMinuit->fBlim[i-1];
00866 Double_t xerr = gMinuit->fWerr[l-1];
00867 Double_t xval = gMinuit->fU[i-1] ;
00868
00869 RooRealVar* var ;
00870 if (varList.getSize()==0) {
00871
00872 if ((xlo<xhi) && !isConst) {
00873 var = new RooRealVar(varName,varName,xval,xlo,xhi) ;
00874 } else {
00875 var = new RooRealVar(varName,varName,xval) ;
00876 }
00877 var->setConstant(isConst) ;
00878 } else {
00879
00880 var = (RooRealVar*) varList.at(i-1)->Clone() ;
00881 var->setConstant(isConst) ;
00882 var->setVal(xval) ;
00883 if (xlo<xhi) {
00884 var->setRange(xlo,xhi) ;
00885 }
00886 if (varName.CompareTo(var->GetName())) {
00887 oocoutI((TObject*)0,Eval) << "RooFitResult::lastMinuitFit: fit parameter '" << varName
00888 << "' stored in variable '" << var->GetName() << "'" << endl ;
00889 }
00890
00891 }
00892
00893 if (isConst) {
00894 constPars.addOwned(*var) ;
00895 } else {
00896 var->setError(xerr) ;
00897 floatPars.addOwned(*var) ;
00898 }
00899 }
00900
00901 Int_t icode,npari,nparx ;
00902 Double_t fmin,edm,errdef ;
00903 gMinuit->mnstat(fmin,edm,errdef,npari,nparx,icode) ;
00904
00905 r->setConstParList(constPars) ;
00906 r->setInitParList(floatPars) ;
00907 r->setFinalParList(floatPars) ;
00908 r->setMinNLL(fmin) ;
00909 r->setEDM(edm) ;
00910 r->setCovQual(icode) ;
00911 r->setStatus(gMinuit->fStatus) ;
00912 r->fillCorrMatrix() ;
00913
00914 return r ;
00915 }
00916
00917
00918
00919
00920 void RooFitResult::setCovarianceMatrix(TMatrixDSym& V)
00921 {
00922
00923
00924
00925 if (_VM) {
00926 delete _VM ;
00927 }
00928 if (_CM) {
00929 delete _CM ;
00930 }
00931
00932
00933 _VM = (TMatrixDSym*) V.Clone() ;
00934
00935
00936 _CM = (TMatrixDSym*) _VM->Clone() ;
00937 for (Int_t i=0 ; i<_CM->GetNrows() ; i++) {
00938 for (Int_t j=0 ; j<_CM->GetNcols() ; j++) {
00939 if (i!=j) {
00940 (*_CM)(i,j) = (*_CM)(i,j) / sqrt((*_CM)(i,i)*(*_CM)(j,j)) ;
00941 }
00942 }
00943 }
00944 for (Int_t i=0 ; i<_CM->GetNrows() ; i++) {
00945 (*_CM)(i,i) = 1.0 ;
00946 }
00947
00948 _covQual = -1 ;
00949 }
00950
00951
00952
00953
00954 TH2* RooFitResult::correlationHist(const char* name) const
00955 {
00956
00957 Int_t n = _CM->GetNcols() ;
00958
00959 TH2D* hh = new TH2D(name,name,n,0,n,n,0,n) ;
00960
00961 for (Int_t i = 0 ; i<n ; i++) {
00962 for (Int_t j = 0 ; j<n; j++) {
00963 hh->Fill(i+0.5,n-j-0.5,(*_CM)(i,j)) ;
00964 }
00965 hh->GetXaxis()->SetBinLabel(i+1,_finalPars->at(i)->GetName()) ;
00966 hh->GetYaxis()->SetBinLabel(n-i,_finalPars->at(i)->GetName()) ;
00967 }
00968 hh->SetMinimum(-1) ;
00969 hh->SetMaximum(+1) ;
00970
00971
00972 return hh ;
00973 }
00974
00975
00976
00977
00978
00979 const TMatrixDSym& RooFitResult::covarianceMatrix() const
00980 {
00981
00982 return *_VM ;
00983 }
00984
00985
00986
00987
00988
00989 TMatrixDSym RooFitResult::reducedCovarianceMatrix(const RooArgList& params) const
00990 {
00991
00992
00993
00994
00995
00996
00997
00998
00999
01000
01001 const TMatrixDSym& V = covarianceMatrix() ;
01002
01003
01004 if (V.GetNcols()==params.getSize()) {
01005 return V ;
01006 }
01007
01008 Double_t det = V.Determinant() ;
01009
01010 if (det<=0) {
01011 coutE(Eval) << "RooFitResult::reducedCovarianceMatrix(" << GetName() << ") ERROR: covariance matrix is not positive definite (|V|="
01012 << det << ") cannot reduce it" << endl ;
01013 throw string("RooFitResult::reducedCovarianceMatrix() ERROR, input covariance matrix is not positive definite") ;
01014 }
01015
01016
01017 RooArgList params2 ;
01018 TIterator* iter = params.createIterator() ;
01019 RooAbsArg* arg ;
01020 while((arg=(RooAbsArg*)iter->Next())) {
01021 if (_finalPars->find(arg->GetName())) {
01022 params2.add(*arg) ;
01023 } else {
01024 coutW(InputArguments) << "RooFitResult::reducedCovarianceMatrix(" << GetName() << ") WARNING input variable "
01025 << arg->GetName() << " was not a floating parameters in fit result and is ignored" << endl ;
01026 }
01027 }
01028 delete iter ;
01029
01030
01031 RooArgList params3 ;
01032 iter = _finalPars->createIterator() ;
01033 while((arg=(RooAbsArg*)iter->Next())) {
01034 if (params2.find(arg->GetName())) {
01035 params3.add(*arg) ;
01036 }
01037 }
01038 delete iter ;
01039
01040
01041 vector<int> map1, map2 ;
01042 for (int i=0 ; i<_finalPars->getSize() ; i++) {
01043 if (params3.find(_finalPars->at(i)->GetName())) {
01044 map1.push_back(i) ;
01045 } else {
01046 map2.push_back(i) ;
01047 }
01048 }
01049
01050
01051
01052 TMatrixDSym S11, S22 ;
01053 TMatrixD S12, S21 ;
01054 RooMultiVarGaussian::blockDecompose(V,map1,map2,S11,S12,S21,S22) ;
01055
01056
01057 TVectorD mu1(map1.size()) ;
01058 for (UInt_t i=0 ; i<map1.size() ; i++) {
01059 mu1(i) = ((RooAbsReal*)_finalPars->at(map1[i]))->getVal() ;
01060 }
01061
01062
01063
01064
01065
01066 TMatrixD S22Inv(TMatrixD::kInverted,S22) ;
01067 TMatrixD S22bar = S11 - S12 * (S22Inv * S21) ;
01068
01069
01070 TMatrixDSym Vred(S22bar.GetNcols()) ;
01071 for (int i=0 ; i<Vred.GetNcols() ; i++) {
01072 for (int j=i ; j<Vred.GetNcols() ; j++) {
01073 Vred(i,j) = (S22bar(i,j) + S22bar(j,i))/2 ;
01074 Vred(j,i) = Vred(i,j) ;
01075 }
01076 }
01077
01078 return Vred ;
01079 }
01080
01081
01082
01083
01084 const TMatrixDSym& RooFitResult::correlationMatrix() const
01085 {
01086
01087 return *_CM ;
01088 }
01089
01090
01091
01092
01093 RooAbsPdf* RooFitResult::createHessePdf(const RooArgSet& params) const
01094 {
01095
01096
01097
01098 const TMatrixDSym& V = covarianceMatrix() ;
01099 Double_t det = V.Determinant() ;
01100
01101 if (det<=0) {
01102 coutE(Eval) << "RooFitResult::createHessePdf(" << GetName() << ") ERROR: covariance matrix is not positive definite (|V|="
01103 << det << ") cannot construct p.d.f" << endl ;
01104 return 0 ;
01105 }
01106
01107
01108 RooArgList params2 ;
01109 TIterator* iter = params.createIterator() ;
01110 RooAbsArg* arg ;
01111 while((arg=(RooAbsArg*)iter->Next())) {
01112 if (_finalPars->find(arg->GetName())) {
01113 params2.add(*arg) ;
01114 } else {
01115 coutW(InputArguments) << "RooFitResult::createHessePdf(" << GetName() << ") WARNING input variable "
01116 << arg->GetName() << " was not a floating parameters in fit result and is ignored" << endl ;
01117 }
01118 }
01119 delete iter ;
01120
01121
01122 RooArgList params3 ;
01123 iter = _finalPars->createIterator() ;
01124 while((arg=(RooAbsArg*)iter->Next())) {
01125 if (params2.find(arg->GetName())) {
01126 params3.add(*arg) ;
01127 }
01128 }
01129 delete iter ;
01130
01131
01132
01133 if (params3.getSize()==_finalPars->getSize()) {
01134
01135 RooArgList mu ;
01136 for (Int_t i=0 ; i<_finalPars->getSize() ; i++) {
01137 RooRealVar* parclone = (RooRealVar*) _finalPars->at(i)->Clone(Form("%s_centralvalue",_finalPars->at(i)->GetName())) ;
01138 parclone->setConstant(kTRUE) ;
01139 mu.add(*parclone) ;
01140 }
01141
01142 string name = Form("pdf_%s",GetName()) ;
01143 string title = Form("P.d.f of %s",GetTitle()) ;
01144
01145
01146 RooAbsPdf* mvg = new RooMultiVarGaussian(name.c_str(),title.c_str(),params3,mu,V) ;
01147 mvg->addOwnedComponents(mu) ;
01148 return mvg ;
01149 }
01150
01151
01152
01153
01154
01155 vector<int> map1, map2 ;
01156 for (int i=0 ; i<_finalPars->getSize() ; i++) {
01157 if (params3.find(_finalPars->at(i)->GetName())) {
01158 map1.push_back(i) ;
01159 } else {
01160 map2.push_back(i) ;
01161 }
01162 }
01163
01164
01165
01166 TMatrixDSym S11, S22 ;
01167 TMatrixD S12, S21 ;
01168 RooMultiVarGaussian::blockDecompose(V,map1,map2,S11,S12,S21,S22) ;
01169
01170
01171 RooArgList mu1 ;
01172 for (UInt_t i=0 ; i<map1.size() ; i++) {
01173 RooRealVar* parclone = (RooRealVar*) _finalPars->at(map1[i])->Clone(Form("%s_centralvalue",_finalPars->at(i)->GetName())) ;
01174 parclone->setConstant(kTRUE) ;
01175 mu1.add(*parclone) ;
01176 }
01177
01178
01179
01180
01181
01182 TMatrixD S22Inv(TMatrixD::kInverted,S22) ;
01183 TMatrixD S22bar = S11 - S12 * (S22Inv * S21) ;
01184
01185
01186 TMatrixDSym Vred(S22bar.GetNcols()) ;
01187 for (int i=0 ; i<Vred.GetNcols() ; i++) {
01188 for (int j=i ; j<Vred.GetNcols() ; j++) {
01189 Vred(i,j) = (S22bar(i,j) + S22bar(j,i))/2 ;
01190 Vred(j,i) = Vred(i,j) ;
01191 }
01192 }
01193 string name = Form("pdf_%s",GetName()) ;
01194 string title = Form("P.d.f of %s",GetTitle()) ;
01195
01196
01197 RooAbsPdf* ret = new RooMultiVarGaussian(name.c_str(),title.c_str(),params3,mu1,Vred) ;
01198 ret->addOwnedComponents(mu1) ;
01199 return ret ;
01200 }
01201
01202
01203
01204
01205 void RooFitResult::SetName(const char *name)
01206 {
01207
01208
01209 if (_dir) _dir->GetList()->Remove(this);
01210 TNamed::SetName(name) ;
01211 if (_dir) _dir->GetList()->Add(this);
01212 }
01213
01214
01215
01216 void RooFitResult::SetNameTitle(const char *name, const char* title)
01217 {
01218
01219
01220 if (_dir) _dir->GetList()->Remove(this);
01221 TNamed::SetNameTitle(name,title) ;
01222 if (_dir) _dir->GetList()->Add(this);
01223 }
01224
01225
01226
01227 void RooFitResult::printName(ostream& os) const
01228 {
01229
01230
01231 os << GetName() ;
01232 }
01233
01234
01235
01236 void RooFitResult::printTitle(ostream& os) const
01237 {
01238
01239
01240 os << GetTitle() ;
01241 }
01242
01243
01244
01245 void RooFitResult::printClassName(ostream& os) const
01246 {
01247
01248
01249 os << IsA()->GetName() ;
01250 }
01251
01252
01253
01254 void RooFitResult::printArgs(ostream& os) const
01255 {
01256
01257
01258 os << "[constPars=" << *_constPars << ",floatPars=" << *_finalPars << "]" ;
01259 }
01260
01261
01262
01263
01264 void RooFitResult::printValue(ostream& os) const
01265 {
01266
01267
01268 os << "(status=" << _status << ",FCNmin=" << _minNLL << ",EDM=" << _edm << ",covQual=" << _covQual << ")" ;
01269 }
01270
01271
01272
01273 Int_t RooFitResult::defaultPrintContents(Option_t* ) const
01274 {
01275
01276
01277 return kName|kClassName|kArgs|kValue ;
01278 }
01279
01280
01281
01282 RooPrintable::StyleOption RooFitResult::defaultPrintStyle(Option_t* opt) const
01283 {
01284
01285 if (!opt || strlen(opt)==0) {
01286 return kStandard ;
01287 }
01288 return RooPrintable::defaultPrintStyle(opt) ;
01289 }
01290
01291
01292
01293 void RooFitResult::Streamer(TBuffer &R__b)
01294 {
01295
01296
01297 if (R__b.IsReading()) {
01298 UInt_t R__s, R__c;
01299 Version_t R__v = R__b.ReadVersion(&R__s, &R__c);
01300 if (R__v>3) {
01301 R__b.ReadClassBuffer(RooFitResult::Class(),this,R__v,R__s,R__c);
01302 } else {
01303
01304 TNamed::Streamer(R__b);
01305 RooPrintable::Streamer(R__b);
01306 RooDirItem::Streamer(R__b);
01307 R__b >> _status;
01308 R__b >> _covQual;
01309 R__b >> _numBadNLL;
01310 R__b >> _minNLL;
01311 R__b >> _edm;
01312 R__b >> _constPars;
01313 R__b >> _initPars;
01314 R__b >> _finalPars;
01315 R__b >> _globalCorr;
01316 _corrMatrix.Streamer(R__b);
01317 R__b.CheckByteCount(R__s, R__c, RooFitResult::IsA());
01318
01319
01320
01321 _CM = new TMatrixDSym(_finalPars->getSize()) ;
01322 _VM = new TMatrixDSym(_CM->GetNcols()) ;
01323 _GC = new TVectorD(_CM->GetNcols()) ;
01324
01325 TIterator *gcIter = _globalCorr->createIterator() ;
01326 TIterator *parIter = _finalPars->createIterator() ;
01327 RooRealVar* gcVal = 0;
01328 for (unsigned int i = 0; i < (unsigned int)_CM->GetNcols() ; ++i) {
01329
01330
01331 gcVal = (RooRealVar*) gcIter->Next() ;
01332 (*_GC)(i) = gcVal->getVal() ;
01333
01334
01335 TIterator* cIter = ((RooArgList*)_corrMatrix.At(i))->createIterator() ;
01336 for (unsigned int it = 0; it < (unsigned int)_CM->GetNcols() ; ++it) {
01337 RooRealVar* cVal = (RooRealVar*) cIter->Next() ;
01338 double value = cVal->getVal() ;
01339 (*_CM)(it,i) = value ;
01340 (*_CM)(i,it) = value;
01341 (*_VM)(it,i) = value*((RooRealVar*)_finalPars->at(i))->getError()*((RooRealVar*)_finalPars->at(it))->getError() ;
01342 (*_VM)(i,it) = (*_VM)(it,i) ;
01343 }
01344 delete cIter ;
01345 }
01346
01347 delete gcIter ;
01348 delete parIter ;
01349 }
01350
01351 } else {
01352 R__b.WriteClassBuffer(RooFitResult::Class(),this);
01353 }
01354 }
01355