RooFitResult.cxx

Go to the documentation of this file.
00001 /*****************************************************************************
00002  * Project: RooFit                                                           *
00003  * Package: RooFitCore                                                       *
00004  * @(#)root/roofitcore:$Id: RooFitResult.cxx 36274 2010-10-11 08:05:03Z wouter $
00005  * Authors:                                                                  *
00006  *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
00007  *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
00008  *                                                                           *
00009  * Copyright (c) 2000-2005, Regents of the University of California          *
00010  *                          and Stanford University. All rights reserved.    *
00011  *                                                                           *
00012  * Redistribution and use in source and binary forms,                        *
00013  * with or without modification, are permitted according to the terms        *
00014  * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
00015  *****************************************************************************/
00016 
00017 //////////////////////////////////////////////////////////////////////////////
00018 //
00019 // RooFitResult is a container class to hold the input and output
00020 // of a PDF fit to a dataset. It contains:
00021 //
00022 //   - Values of all constant parameters
00023 //   - Initial and final values of floating parameters with error
00024 //   - Correlation matrix and global correlation coefficients
00025 //   - NLL and EDM at mininum
00026 //
00027 // No references to the fitted PDF and dataset are stored
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   // Constructor with name and title
00070   // coverity[UNINIT_CTOR]
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   // Copy constructor
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   // Destructor
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   // Fill the list of constant parameters
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   // Fill the list of initial values of the floating parameters 
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   // Fill the list of final values of the floating parameters 
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   // Add objects to a 2D plot that represent the fit results for the
00192   // two named parameters.  The input frame with the objects added is
00193   // returned, or zero in case of an error.  Which objects are added
00194   // are determined by the options string which should be a concatenation
00195   // of the following (not case sensitive):
00196   //
00197   //   M - a marker at the best fit result
00198   //   E - an error ellipse calculated at 1-sigma using the error matrix at the minimum
00199   //   1 - the 1-sigma error bar for parameter 1
00200   //   2 - the 1-sigma error bar for parameter 2
00201   //   B - the bounding box for the error ellipse
00202   //   H - a line and horizontal axis for reading off the correlation coefficient
00203   //   V - a line and vertical axis for reading off the correlation coefficient
00204   //   A - draw axes for reading off the correlation coefficients with the H or V options
00205   //
00206   // You can change the attributes of objects in the returned RooPlot using the
00207   // various RooPlot::getAttXxx(name) member functions, e.g.
00208   //
00209   //   plot->getAttLine("contour")->SetLineStyle(kDashed);
00210   //
00211   // Use plot->Print() for a list of all objects and their names (unfortunately most
00212   // of the ROOT builtin graphics objects like TLine are unnamed). Drag the left mouse
00213   // button along the labels of either axis button to interactively zoom in a plot.
00214 
00215   // lookup the input parameters by name: we require that they were floated in our fit
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   // options are not case sensitive
00228   TString opt(options);
00229   opt.ToUpper();
00230 
00231   // lookup the 2x2 covariance matrix elements for these variables
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   // add a 1-sigma error ellipse, if requested
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   // add the error bar for parameter 1, if requested
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   // add a marker at the fitted value, if requested
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   // Return a list of floating parameter values that are perturbed from the final
00307   // fit values by random amounts sampled from the covariance matrix. The returned
00308   // object is overwritten with each call and belongs to the RooFitResult. Uses
00309   // the "square root method" to decompose the covariance matrix, which makes inverting
00310   // it unnecessary.
00311   
00312   Int_t nPar= _finalPars->getSize();
00313   if(0 == _randomPars) { // first-time initialization
00314     assert(0 != _finalPars);
00315     // create the list of random values to fill
00316     _randomPars= (RooArgList*)_finalPars->snapshot();
00317     // calculate the elements of the upper-triangular matrix L that gives Lt*L = C
00318     // where Lt is the transpose of L (the "square-root method")
00319     TMatrix L(nPar,nPar);
00320     for(Int_t iPar= 0; iPar < nPar; iPar++) {
00321       // calculate the diagonal term first
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       // then the off-diagonal terms
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     // remember Lt
00338     _Lt= new TMatrix(TMatrix::kTransposed,L);
00339   }
00340   else {
00341     // reset to the final fit values
00342     *_randomPars= *_finalPars;
00343   }
00344 
00345   // create a vector of unit Gaussian variables
00346   TVector g(nPar);
00347   for(Int_t k= 0; k < nPar; k++) g(k)= RooRandom::gaussian();
00348   // multiply this vector by Lt to introduce the appropriate correlations
00349   g*= (*_Lt);
00350   // add the mean value offsets and store the results
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   // Return the correlation between parameters 'par1' and 'par2'
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   // Return the set of correlation coefficients of parameter 'par' with
00386   // all other floating parameters
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   // Return the global correlation of the named parameter
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   // Return the list of all global correlations
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   // Return a correlation matrix element addressed with numeric indices.
00444   return (*_CM)(row,col) ;
00445 }
00446 
00447 
00448 //_____________________________________________________________________________
00449 Double_t RooFitResult::covariance(Int_t row, Int_t col) const 
00450 {
00451   // Return the covariance matrix element addressed with numeric indices.
00452   return (*_VM)(row,col) ;
00453 }
00454 
00455 
00456 
00457 //_____________________________________________________________________________
00458 void RooFitResult::printMultiline(ostream& os, Int_t /*contents*/, Bool_t verbose, TString indent) const
00459 {
00460   // Print fit result to stream 'os'. In Verbose mode, the contant parameters and
00461   // the initial and final values of the floating parameters are printed. 
00462   // Standard mode only the final values of the floating parameters are printed
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     // Has any parameter asymmetric errors?
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   // Function called by RooMinimizer
00554 
00555   // Sanity check
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   // Delete eventual prevous correlation data holders
00567   if (_CM) delete _CM ;
00568   if (_VM) delete _VM ;
00569   if (_GC) delete _GC ;
00570 
00571   // Build holding arrays for correlation coefficients
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   // Sanity check
00589   if (!_CM) return ;
00590 
00591   // Delete eventual prevous correlation data holders
00592   if (_globalCorr) delete _globalCorr ;
00593   _corrMatrix.Delete();
00594 
00595   // Build holding arrays for correlation coefficients
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     // Create global correlation value holder
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     // Create array with correlation holders for this parameter
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     // Find the next global correlation slot to fill, skipping fixed parameters
00642     gcVal = (RooRealVar*) gcIter->Next() ;
00643     gcVal->setVal((*_GC)(i)) ; // WVE FIX THIS 
00644 
00645     // Fill a row of the correlation matrix
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   // Internal utility method to extract the correlation matrix and the
00669   // global correlation coefficients from the MINUIT memory buffer and
00670   // fill the internal arrays.
00671 
00672   // Sanity check
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   // Delete eventual prevous correlation data holders
00684   if (_CM) delete _CM ;
00685   if (_VM) delete _VM ;
00686   if (_GC) delete _GC ;
00687 
00688   // Build holding arrays for correlation coefficients
00689   _CM = new TMatrixDSym(_initPars->getSize()) ;
00690   _VM = new TMatrixDSym(_initPars->getSize()) ;
00691   _GC = new TVectorD(_initPars->getSize()) ;
00692 
00693   // Extract correlation information for MINUIT (code taken from TMinuit::mnmatu() )
00694 
00695   // WVE: This code directly manipulates minuit internal workspace, 
00696   //      if TMinuit code changes this may need updating
00697   Int_t ndex, i, j, m, n, ncoef, nparm, /*id,*/ it, ix ;
00698   Int_t ndi, ndj /*, iso, isw2, isw5*/;
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     // Fill a row of the correlation matrix
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 /*verbose*/) const 
00734 {
00735   // Return true if this fit result is identical to other within tolerance 'tol' on fitted values
00736   // and tolerance 'tolCor' on correlation coefficients
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   // Only examine correlations for cases with >1 floating paramater
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   // Import the results of the last fit performed by gMinuit, interpreting
00830   // the fit parameters as the given varList of parameters.
00831 
00832   // Verify length of supplied varList
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   // Verify that all members of varList are of type RooRealVar
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   // Extract names of fit parameters from MINUIT 
00853   // and construct corresponding RooRealVars
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   // Store externally provided correlation matrix in his RooFitResult ;
00923 
00924   // Delete any previous matrices
00925   if (_VM) {
00926     delete _VM ;
00927   }
00928   if (_CM) {
00929     delete _CM ;
00930   }
00931   
00932   // Clone input covariance matrix ;
00933   _VM = (TMatrixDSym*) V.Clone() ;
00934 
00935   // Now construct correlation matrix from it
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   // Return TH2D of correlation matrix 
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   // Return covariance matrix 
00982   return *_VM ;
00983 }
00984 
00985 
00986 
00987 
00988 //_____________________________________________________________________________
00989 TMatrixDSym RooFitResult::reducedCovarianceMatrix(const RooArgList& params) const 
00990 {
00991   // Return a reduced covariance matrix, which is calculated as
00992   //        ___                   -1
00993   // Vred = V22  = V11 - V12 * V22   * V21
00994   //
00995   // Where V11,V12,V21,V22 represent a block decomposition of the covariance matrix into observables that
00996   // are propagated (labeled by index '1') and that are not propagated (labeled by index '2'), and V22bar
00997   // is the Shur complement of V22, calculated as shown above  
00998   //
00999   // (Note that Vred is _not_ a simple sub-matrix of V)
01000 
01001   const TMatrixDSym& V = covarianceMatrix() ;
01002 
01003   // Handle case where V==Vred here
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   // Make sure that all given params were floating parameters in the represented fit
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   // Need to order params in vector in same order as in covariance matrix
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   // Find (subset) of parameters that are stored in the covariance matrix
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   // Rearrange matrix in block form with 'params' first and 'others' last
01051   // (preserving relative order) 
01052   TMatrixDSym S11, S22 ;
01053   TMatrixD S12, S21 ;
01054   RooMultiVarGaussian::blockDecompose(V,map1,map2,S11,S12,S21,S22) ;
01055 
01056   // Calculate offset vectors mu1 and mu2
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   // Constructed conditional matrix form         -1
01063   // F(X1|X2) --> CovI --> S22bar = S11 - S12 S22  S21
01064   
01065   // Do eigenvalue decomposition
01066   TMatrixD S22Inv(TMatrixD::kInverted,S22) ;
01067   TMatrixD S22bar =  S11 - S12 * (S22Inv * S21) ;
01068 
01069   // Convert explicitly to symmetric form
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   // Return correlation matrix ;
01087   return *_CM ;
01088 }
01089 
01090 
01091 
01092 //_____________________________________________________________________________
01093 RooAbsPdf* RooFitResult::createHessePdf(const RooArgSet& params) const
01094 {
01095   // Return a p.d.f that represents the fit result as a multi-variate probability densisty
01096   // function on the floating fit parameters, including correlations
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   // Make sure that all given params were floating parameters in the represented fit
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   // Need to order params in vector in same order as in covariance matrix
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   // Handle special case of representing full covariance matrix here
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     // Create p.d.f.
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   // Handle case of conditional p.d.f. MVG(p1|p2) here
01153 
01154   // Find (subset) of parameters that are stored in the covariance matrix
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   // Rearrange matrix in block form with 'params' first and 'others' last
01165   // (preserving relative order) 
01166   TMatrixDSym S11, S22 ;
01167   TMatrixD S12, S21 ;
01168   RooMultiVarGaussian::blockDecompose(V,map1,map2,S11,S12,S21,S22) ;
01169 
01170   // Calculate offset vectors mu1 and mu2
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   // Constructed conditional matrix form         -1
01179   // F(X1|X2) --> CovI --> S22bar = S11 - S12 S22  S21
01180   
01181   // Do eigenvalue decomposition
01182   TMatrixD S22Inv(TMatrixD::kInverted,S22) ;
01183   TMatrixD S22bar =  S11 - S12 * (S22Inv * S21) ;
01184 
01185   // Convert explicitly to symmetric form
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   // Create p.d.f.
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   // Change name of RooFitResult object
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   // Change name and title of RooFitResult object
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   // Print name of fit result
01230 
01231   os << GetName() ;
01232 }
01233 
01234 
01235 //_____________________________________________________________________________
01236 void RooFitResult::printTitle(ostream& os) const 
01237 {
01238   // Print title of fit result
01239 
01240   os << GetTitle() ;
01241 }
01242 
01243 
01244 //_____________________________________________________________________________
01245 void RooFitResult::printClassName(ostream& os) const 
01246 {
01247   // Print class name of fit result
01248 
01249   os << IsA()->GetName() ;
01250 }
01251 
01252 
01253 //_____________________________________________________________________________
01254 void RooFitResult::printArgs(ostream& os) const 
01255 {
01256   // Print arguments of fit result, i.e. the parameters of the fit
01257 
01258   os << "[constPars=" << *_constPars << ",floatPars=" << *_finalPars << "]" ;
01259 }
01260 
01261 
01262 
01263 //_____________________________________________________________________________
01264 void RooFitResult::printValue(ostream& os) const 
01265 {
01266   // Print the value of the fit result, i.e.g the status, minimized FCN, edm and covariance quality code
01267 
01268   os << "(status=" << _status << ",FCNmin=" << _minNLL << ",EDM=" << _edm << ",covQual=" << _covQual << ")" ;
01269 }
01270 
01271 
01272 //_____________________________________________________________________________
01273 Int_t RooFitResult::defaultPrintContents(Option_t* /*opt*/) const 
01274 {
01275   // Configure default contents to be printed
01276 
01277   return kName|kClassName|kArgs|kValue ;
01278 }
01279 
01280 
01281 //_____________________________________________________________________________
01282 RooPrintable::StyleOption RooFitResult::defaultPrintStyle(Option_t* opt) const 
01283 {
01284   // Configure mapping of Print() arguments to RooPrintable print styles
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    // Stream an object of class RooFitResult.
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       // backward compatibitily streaming 
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       // Now fill new-style covariance and correlation matrix information
01320       // from legacy form
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         // Find the next global correlation slot to fill, skipping fixed parameters
01331         gcVal = (RooRealVar*) gcIter->Next() ;
01332         (*_GC)(i) = gcVal->getVal() ;
01333         
01334         // Fill a row of the correlation matrix
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 

Generated on Tue Jul 5 15:06:40 2011 for ROOT_528-00b_version by  doxygen 1.5.1