RooMinimizer.cxx

Go to the documentation of this file.
00001 /*****************************************************************************
00002  * Project: RooFit                                                           *
00003  * Package: RooFitCore                                                       *
00004  * @(#)root/roofitcore:$Id: RooMinimizer.cxx 34286 2010-07-01 20:38:57Z rdm $
00005  * Authors:                                                                  *
00006  *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
00007  *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
00008  *   AL, Alfio Lazzaro,   INFN Milan,        alfio.lazzaro@mi.infn.it        *
00009  *                                                                           *
00010  * Redistribution and use in source and binary forms,                        *
00011  * with or without modification, are permitted according to the terms        *
00012  * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
00013  *****************************************************************************/
00014 
00015 //////////////////////////////////////////////////////////////////////////////
00016 //
00017 // BEGIN_HTML
00018 // RooMinimizer is a wrapper class around ROOT::Fit:Fitter that
00019 // provides a seamless interface between the minimizer functionality
00020 // and the native RooFit interface.
00021 // <p>
00022 // By default the Minimizer is MINUIT.
00023 // <p>
00024 // RooMinimizer can minimize any RooAbsReal function with respect to
00025 // its parameters. Usual choices for minimization are RooNLLVar
00026 // and RooChi2Var
00027 // <p>
00028 // RooMinimizer has methods corresponding to MINUIT functions like
00029 // hesse(), migrad(), minos() etc. In each of these function calls
00030 // the state of the MINUIT engine is synchronized with the state
00031 // of the RooFit variables: any change in variables, change
00032 // in the constant status etc is forwarded to MINUIT prior to
00033 // execution of the MINUIT call. Afterwards the RooFit objects
00034 // are resynchronized with the output state of MINUIT: changes
00035 // parameter values, errors are propagated.
00036 // <p>
00037 // Various methods are available to control verbosity, profiling,
00038 // automatic PDF optimization.
00039 // END_HTML
00040 //
00041 
00042 #ifndef __ROOFIT_NOROOMINIMIZER
00043 
00044 #include "RooFit.h"
00045 #include "Riostream.h"
00046 
00047 #include "TClass.h"
00048 
00049 #include <fstream>
00050 #include <iomanip>
00051 
00052 #include "TH1.h"
00053 #include "TH2.h"
00054 #include "TMarker.h"
00055 #include "TGraph.h"
00056 #include "TStopwatch.h"
00057 #include "TDirectory.h"
00058 #include "TMatrixDSym.h"
00059 
00060 #include "RooArgSet.h"
00061 #include "RooArgList.h"
00062 #include "RooAbsReal.h"
00063 #include "RooAbsRealLValue.h"
00064 #include "RooRealVar.h"
00065 #include "RooAbsPdf.h"
00066 #include "RooSentinel.h"
00067 #include "RooMsgService.h"
00068 #include "RooPlot.h"
00069 
00070 
00071 #include "RooMinimizer.h"
00072 #include "RooFitResult.h"
00073 
00074 #include "Math/Minimizer.h"
00075 
00076 #if (__GNUC__==3&&__GNUC_MINOR__==2&&__GNUC_PATCHLEVEL__==3)
00077 char* operator+( streampos&, char* );
00078 #endif
00079 
00080 ClassImp(RooMinimizer)
00081 ;
00082 
00083 ROOT::Fit::Fitter *RooMinimizer::_theFitter = 0 ;
00084 
00085 
00086 
00087 //_____________________________________________________________________________
00088 void RooMinimizer::cleanup()
00089 {
00090   // Cleanup method called by atexit handler installed by RooSentinel
00091   // to delete all global heap objects when the program is terminated
00092   if (_theFitter) {
00093     delete _theFitter ;
00094     _theFitter =0 ;
00095   }
00096 }
00097 
00098 
00099 
00100 //_____________________________________________________________________________
00101 RooMinimizer::RooMinimizer(RooAbsReal& function)
00102 {
00103   // Construct MINUIT interface to given function. Function can be anything,
00104   // but is typically a -log(likelihood) implemented by RooNLLVar or a chi^2
00105   // (implemented by RooChi2Var). Other frequent use cases are a RooAddition
00106   // of a RooNLLVar plus a penalty or constraint term. This class propagates
00107   // all RooFit information (floating parameters, their values and errors)
00108   // to MINUIT before each MINUIT call and propagates all MINUIT information
00109   // back to the RooFit object at the end of each call (updated parameter
00110   // values, their (asymmetric errors) etc. The default MINUIT error level
00111   // for HESSE and MINOS error analysis is taken from the defaultErrorLevel()
00112   // value of the input function.
00113 
00114   RooSentinel::activate() ;
00115 
00116   // Store function reference
00117   _extV = 0 ;
00118   _func = &function ;
00119   _optConst = kFALSE ;
00120   _verbose = kFALSE ;
00121   _profile = kFALSE ;
00122   _printLevel = 1 ;
00123   _minimizerType = "Minuit"; // default minimizer
00124 
00125   if (_theFitter) delete _theFitter ;
00126   _theFitter = new ROOT::Fit::Fitter;
00127   _fcn = new RooMinimizerFcn(_func,this,_verbose);
00128   _theFitter->Config().SetMinimizer(_minimizerType.c_str());
00129   setEps(1.0); // default tolerance
00130   // default max number of calls
00131   _theFitter->Config().MinimizerOptions().SetMaxIterations(500*_fcn->NDim());
00132   _theFitter->Config().MinimizerOptions().SetMaxFunctionCalls(500*_fcn->NDim());
00133 
00134   // Shut up for now
00135   setPrintLevel(-1) ;
00136 
00137   // Use +0.5 for 1-sigma errors
00138   setErrorLevel(_func->defaultErrorLevel()) ;
00139 
00140   // Declare our parameters to MINUIT
00141   _fcn->Synchronize(_theFitter->Config().ParamsSettings(),
00142                     _optConst,_verbose) ;
00143 
00144   // Now set default verbosity
00145   if (RooMsgService::instance().silentMode()) {
00146     setPrintLevel(-1) ;
00147   } else {
00148     setPrintLevel(1) ;
00149   }
00150 }
00151 
00152 
00153 
00154 //_____________________________________________________________________________
00155 RooMinimizer::~RooMinimizer()
00156 {
00157   // Destructor
00158 
00159   if (_extV) {
00160     delete _extV ;
00161   }
00162 
00163   if (_fcn) {
00164     delete _fcn;
00165   }
00166 
00167 }
00168 
00169 
00170 
00171 //_____________________________________________________________________________
00172 void RooMinimizer::setStrategy(Int_t istrat)
00173 {
00174   // Change MINUIT strategy to istrat. Accepted codes
00175   // are 0,1,2 and represent MINUIT strategies for dealing
00176   // most efficiently with fast FCNs (0), expensive FCNs (2)
00177   // and 'intermediate' FCNs (1)
00178 
00179   _theFitter->Config().MinimizerOptions().SetStrategy(istrat);
00180 
00181 }
00182 
00183 
00184 
00185 //_____________________________________________________________________________
00186 void RooMinimizer::setErrorLevel(Double_t level)
00187 {
00188   // Set the level for MINUIT error analysis to the given
00189   // value. This function overrides the default value
00190   // that is taken in the RooMinimizer constructor from
00191   // the defaultErrorLevel() method of the input function
00192 
00193   _theFitter->Config().MinimizerOptions().SetErrorDef(level);
00194 
00195 }
00196 
00197 
00198 
00199 //_____________________________________________________________________________
00200 void RooMinimizer::setEps(Double_t eps)
00201 {
00202   // Change MINUIT epsilon
00203 
00204   _theFitter->Config().MinimizerOptions().SetTolerance(eps);
00205 
00206 }
00207 
00208 
00209 
00210 
00211 //_____________________________________________________________________________
00212 void RooMinimizer::setMinimizerType(const char* type)
00213 {
00214   // Choose the minimzer algorithm.
00215   _minimizerType = type;
00216 }
00217 
00218 
00219 
00220 
00221 //_____________________________________________________________________________
00222 RooFitResult* RooMinimizer::fit(const char* options)
00223 {
00224   // Parse traditional RooAbsPdf::fitTo driver options
00225   //
00226   //  m - Run Migrad only
00227   //  h - Run Hesse to estimate errors
00228   //  v - Verbose mode
00229   //  l - Log parameters after each Minuit steps to file
00230   //  t - Activate profile timer
00231   //  r - Save fit result
00232   //  0 - Run Migrad with strategy 0
00233 
00234   TString opts(options) ;
00235   opts.ToLower() ;
00236 
00237   // Initial configuration
00238   if (opts.Contains("v")) setVerbose(1) ;
00239   if (opts.Contains("t")) setProfile(1) ;
00240   if (opts.Contains("l")) setLogFile(Form("%s.log",_func->GetName())) ;
00241   if (opts.Contains("c")) optimizeConst(1) ;
00242 
00243   // Fitting steps
00244   if (opts.Contains("0")) setStrategy(0) ;
00245   migrad() ;
00246   if (opts.Contains("0")) setStrategy(1) ;
00247   if (opts.Contains("h")||!opts.Contains("m")) hesse() ;
00248   if (!opts.Contains("m")) minos() ;
00249 
00250   return (opts.Contains("r")) ? save() : 0 ;
00251 }
00252 
00253 
00254 
00255 
00256 //_____________________________________________________________________________
00257 Int_t RooMinimizer::minimize(const char* type, const char* alg)
00258 {
00259   _fcn->Synchronize(_theFitter->Config().ParamsSettings(),
00260                     _optConst,_verbose) ;
00261 
00262   _theFitter->Config().SetMinimizer(type,alg);
00263 
00264   profileStart() ;
00265   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
00266   RooAbsReal::clearEvalErrorLog() ;
00267 
00268   bool ret = _theFitter->FitFCN(*_fcn);
00269   _status = ((ret) ? _theFitter->Result().Status() : -1);
00270 
00271   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
00272   profileStop() ;
00273   _fcn->BackProp(_theFitter->Result());
00274 
00275   return _status ;
00276 }
00277 
00278 
00279 
00280 //_____________________________________________________________________________
00281 Int_t RooMinimizer::migrad()
00282 {
00283   // Execute MIGRAD. Changes in parameter values
00284   // and calculated errors are automatically
00285   // propagated back the RooRealVars representing
00286   // the floating parameters in the MINUIT operation
00287 
00288   _fcn->Synchronize(_theFitter->Config().ParamsSettings(),
00289                     _optConst,_verbose) ;
00290   profileStart() ;
00291   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
00292   RooAbsReal::clearEvalErrorLog() ;
00293 
00294   _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migrad");
00295   bool ret = _theFitter->FitFCN(*_fcn);
00296   _status = ((ret) ? _theFitter->Result().Status() : -1);
00297 
00298   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
00299   profileStop() ;
00300   _fcn->BackProp(_theFitter->Result());
00301 
00302   return _status ;
00303 }
00304 
00305 
00306 
00307 //_____________________________________________________________________________
00308 Int_t RooMinimizer::hesse()
00309 {
00310   // Execute HESSE. Changes in parameter values
00311   // and calculated errors are automatically
00312   // propagated back the RooRealVars representing
00313   // the floating parameters in the MINUIT operation
00314 
00315   if (_theFitter->GetMinimizer()==0) {
00316     coutW(Minimization) << "RooMinimizer::hesse: Error, run Migrad before Hesse!"
00317                         << endl ;
00318     _status = -1;
00319   }
00320   else {
00321 
00322     _fcn->Synchronize(_theFitter->Config().ParamsSettings(),
00323                     _optConst,_verbose) ;
00324     profileStart() ;
00325     RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
00326     RooAbsReal::clearEvalErrorLog() ;
00327 
00328     _theFitter->Config().SetMinimizer(_minimizerType.c_str());
00329     bool ret = _theFitter->CalculateHessErrors();
00330     _status = ((ret) ? _theFitter->Result().Status() : -1);
00331 
00332     RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
00333     profileStop() ;
00334     _fcn->BackProp(_theFitter->Result());
00335 
00336   }
00337 
00338   return _status ;
00339 
00340 }
00341 
00342 //_____________________________________________________________________________
00343 Int_t RooMinimizer::minos()
00344 {
00345   // Execute MINOS. Changes in parameter values
00346   // and calculated errors are automatically
00347   // propagated back the RooRealVars representing
00348   // the floating parameters in the MINUIT operation
00349 
00350   if (_theFitter->GetMinimizer()==0) {
00351     coutW(Minimization) << "RooMinimizer::minos: Error, run Migrad before Minos!"
00352                         << endl ;
00353     _status = -1;
00354   }
00355   else {
00356 
00357     _fcn->Synchronize(_theFitter->Config().ParamsSettings(),
00358                       _optConst,_verbose) ;
00359     profileStart() ;
00360     RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
00361     RooAbsReal::clearEvalErrorLog() ;
00362 
00363     _theFitter->Config().SetMinimizer(_minimizerType.c_str());
00364     bool ret = _theFitter->CalculateMinosErrors();
00365     _status = ((ret) ? _theFitter->Result().Status() : -1);
00366 
00367     RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
00368     profileStop() ;
00369     _fcn->BackProp(_theFitter->Result());
00370   }
00371 
00372   return _status ;
00373 
00374 }
00375 
00376 
00377 //_____________________________________________________________________________
00378 Int_t RooMinimizer::minos(const RooArgSet& minosParamList)
00379 {
00380   // Execute MINOS for given list of parameters. Changes in parameter values
00381   // and calculated errors are automatically
00382   // propagated back the RooRealVars representing
00383   // the floating parameters in the MINUIT operation
00384 
00385   if (_theFitter->GetMinimizer()==0) {
00386     coutW(Minimization) << "RooMinimizer::minos: Error, run Migrad before Minos!"
00387                         << endl ;
00388     _status = -1;
00389   }
00390   else if (minosParamList.getSize()>0) {
00391 
00392     _fcn->Synchronize(_theFitter->Config().ParamsSettings(),
00393                       _optConst,_verbose) ;
00394     profileStart() ;
00395     RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
00396     RooAbsReal::clearEvalErrorLog() ;
00397 
00398     // get list of parameters for Minos
00399     TIterator* aIter = minosParamList.createIterator() ;
00400     RooAbsArg* arg ;
00401     std::vector<unsigned int> paramInd;
00402     while((arg=(RooAbsArg*)aIter->Next())) {
00403       RooAbsArg* par = _fcn->GetFloatParamList()->find(arg->GetName());
00404       if (par && !par->isConstant()) {
00405         Int_t index = _fcn->GetFloatParamList()->index(par);
00406         paramInd.push_back(index);
00407       }
00408     }
00409     delete aIter ;
00410 
00411     if (paramInd.size()) {
00412       // set the parameter indeces
00413       _theFitter->Config().SetMinosErrors(paramInd);
00414 
00415       _theFitter->Config().SetMinimizer(_minimizerType.c_str());
00416       bool ret = _theFitter->CalculateMinosErrors();
00417       _status = ((ret) ? _theFitter->Result().Status() : -1);
00418 
00419     }
00420 
00421     RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
00422     profileStop() ;
00423     _fcn->BackProp(_theFitter->Result());
00424 
00425   }
00426 
00427   return _status ;
00428 }
00429 
00430 
00431 
00432 //_____________________________________________________________________________
00433 Int_t RooMinimizer::seek()
00434 {
00435   // Execute SEEK. Changes in parameter values
00436   // and calculated errors are automatically
00437   // propagated back the RooRealVars representing
00438   // the floating parameters in the MINUIT operation
00439 
00440   _fcn->Synchronize(_theFitter->Config().ParamsSettings(),
00441                     _optConst,_verbose) ;
00442   profileStart() ;
00443   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
00444   RooAbsReal::clearEvalErrorLog() ;
00445 
00446   _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"seek");
00447   bool ret = _theFitter->FitFCN(*_fcn);
00448   _status = ((ret) ? _theFitter->Result().Status() : -1);
00449 
00450   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
00451   profileStop() ;
00452   _fcn->BackProp(_theFitter->Result());
00453 
00454   return _status ;
00455 }
00456 
00457 
00458 
00459 //_____________________________________________________________________________
00460 Int_t RooMinimizer::simplex()
00461 {
00462   // Execute SIMPLEX. Changes in parameter values
00463   // and calculated errors are automatically
00464   // propagated back the RooRealVars representing
00465   // the floating parameters in the MINUIT operation
00466 
00467   _fcn->Synchronize(_theFitter->Config().ParamsSettings(),
00468                     _optConst,_verbose) ;
00469   profileStart() ;
00470   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
00471   RooAbsReal::clearEvalErrorLog() ;
00472 
00473   _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"simplex");
00474   bool ret = _theFitter->FitFCN(*_fcn);
00475   _status = ((ret) ? _theFitter->Result().Status() : -1);
00476 
00477   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
00478   profileStop() ;
00479   _fcn->BackProp(_theFitter->Result());
00480 
00481   return _status ;
00482 }
00483 
00484 
00485 
00486 //_____________________________________________________________________________
00487 Int_t RooMinimizer::improve()
00488 {
00489   // Execute IMPROVE. Changes in parameter values
00490   // and calculated errors are automatically
00491   // propagated back the RooRealVars representing
00492   // the floating parameters in the MINUIT operation
00493 
00494   _fcn->Synchronize(_theFitter->Config().ParamsSettings(),
00495                     _optConst,_verbose) ;
00496   profileStart() ;
00497   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
00498   RooAbsReal::clearEvalErrorLog() ;
00499 
00500   _theFitter->Config().SetMinimizer(_minimizerType.c_str(),"migradimproved");
00501   bool ret = _theFitter->FitFCN(*_fcn);
00502   _status = ((ret) ? _theFitter->Result().Status() : -1);
00503 
00504   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
00505   profileStop() ;
00506   _fcn->BackProp(_theFitter->Result());
00507 
00508   return _status ;
00509 }
00510 
00511 
00512 
00513 //_____________________________________________________________________________
00514 Int_t RooMinimizer::setPrintLevel(Int_t newLevel)
00515 {
00516   // Change the MINUIT internal printing level
00517 
00518   Int_t ret = _printLevel ;
00519   _theFitter->Config().MinimizerOptions().SetPrintLevel(newLevel+1);
00520   _printLevel = newLevel+1 ;
00521   return ret ;
00522 }
00523 
00524 //_____________________________________________________________________________
00525 void RooMinimizer::optimizeConst(Bool_t flag)
00526 {
00527   // If flag is true, perform constant term optimization on
00528   // function being minimized.
00529 
00530   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::CollectErrors) ;
00531 
00532   if (_optConst && !flag){
00533     if (_printLevel>-1) coutI(Minimization) << "RooMinimizer::optimizeConst: deactivating const optimization" << endl ;
00534     _func->constOptimizeTestStatistic(RooAbsArg::DeActivate) ;
00535     _optConst = flag ;
00536   } else if (!_optConst && flag) {
00537     if (_printLevel>-1) coutI(Minimization) << "RooMinimizer::optimizeConst: activating const optimization" << endl ;
00538     _func->constOptimizeTestStatistic(RooAbsArg::Activate) ;
00539     _optConst = flag ;
00540   } else if (_optConst && flag) {
00541     if (_printLevel>-1) coutI(Minimization) << "RooMinimizer::optimizeConst: const optimization already active" << endl ;
00542   } else {
00543     if (_printLevel>-1) coutI(Minimization) << "RooMinimizer::optimizeConst: const optimization wasn't active" << endl ;
00544   }
00545 
00546   RooAbsReal::setEvalErrorLoggingMode(RooAbsReal::PrintErrors) ;
00547 
00548 }
00549 
00550 
00551 
00552 //_____________________________________________________________________________
00553 RooFitResult* RooMinimizer::save(const char* userName, const char* userTitle)
00554 {
00555   // Save and return a RooFitResult snaphot of current minimizer status.
00556   // This snapshot contains the values of all constant parameters,
00557   // the value of all floating parameters at RooMinimizer construction and
00558   // after the last MINUIT operation, the MINUIT status, variance quality,
00559   // EDM setting, number of calls with evaluation problems, the minimized
00560   // function value and the full correlation matrix
00561 
00562   if (_theFitter->GetMinimizer()==0) {
00563     coutW(Minimization) << "RooMinimizer::save: Error, run minimization before!"
00564                         << endl ;
00565     return 0;
00566   }
00567 
00568   TString name,title ;
00569   name = userName ? userName : Form("%s", _func->GetName()) ;
00570   title = userTitle ? userTitle : Form("%s", _func->GetTitle()) ;
00571   RooFitResult* fitRes = new RooFitResult(name,title) ;
00572 
00573   // Move eventual fixed paramaters in floatList to constList
00574   Int_t i ;
00575   RooArgList saveConstList(*(_fcn->GetConstParamList())) ;
00576   RooArgList saveFloatInitList(*(_fcn->GetInitFloatParamList())) ;
00577   RooArgList saveFloatFinalList(*(_fcn->GetFloatParamList())) ;
00578   for (i=0 ; i<_fcn->GetFloatParamList()->getSize() ; i++) {
00579     RooAbsArg* par = _fcn->GetFloatParamList()->at(i) ;
00580     if (par->isConstant()) {
00581       saveFloatInitList.remove(*saveFloatInitList.find(par->GetName()),kTRUE) ;
00582       saveFloatFinalList.remove(*par) ;
00583       saveConstList.add(*par) ;
00584     }
00585   }
00586   saveConstList.sort() ;
00587 
00588   fitRes->setConstParList(saveConstList) ;
00589   fitRes->setInitParList(saveFloatInitList) ;
00590 
00591   fitRes->setStatus(_status) ;
00592   fitRes->setCovQual(_theFitter->GetMinimizer()->CovMatrixStatus()) ;
00593   fitRes->setMinNLL(_theFitter->Result().MinFcnValue()) ;
00594   fitRes->setNumInvalidNLL(_fcn->GetNumInvalidNLL()) ;
00595   fitRes->setEDM(_theFitter->Result().Edm()) ;
00596   fitRes->setFinalParList(saveFloatFinalList) ;
00597   if (!_extV) {
00598     std::vector<double> globalCC;
00599     TMatrixDSym corrs(_theFitter->Result().Parameters().size()) ;
00600     TMatrixDSym covs(_theFitter->Result().Parameters().size()) ;
00601     for (UInt_t ic=0; ic<_theFitter->Result().Parameters().size(); ic++) {
00602       globalCC.push_back(_theFitter->Result().GlobalCC(ic));
00603       for (UInt_t ii=0; ii<_theFitter->Result().Parameters().size(); ii++) {
00604         corrs(ic,ii) = _theFitter->Result().Correlation(ic,ii);
00605         covs(ic,ii) = _theFitter->Result().CovMatrix(ic,ii);
00606       }
00607     }
00608     fitRes->fillCorrMatrix(globalCC,corrs,covs) ;
00609   } else {
00610     fitRes->setCovarianceMatrix(*_extV) ;
00611   }
00612 
00613   return fitRes ;
00614 
00615 }
00616 
00617 //_____________________________________________________________________________
00618 RooPlot* RooMinimizer::contour(RooRealVar& var1, RooRealVar& var2,
00619                                Double_t n1, Double_t n2, Double_t n3,
00620                                Double_t n4, Double_t n5, Double_t n6)
00621 {
00622   // Create and draw a TH2 with the error contours in parameters var1 and v2 at up to 6 'sigma' settings
00623   // where 'sigma' is calculated as n*n*errorLevel
00624 
00625   RooArgList* params = _fcn->GetFloatParamList() ;
00626   RooArgList* paramSave = (RooArgList*) params->snapshot() ;
00627 
00628   // Verify that both variables are floating parameters of PDF
00629   Int_t index1= _fcn->GetFloatParamList()->index(&var1);
00630   if(index1 < 0) {
00631     coutE(Minimization) << "RooMinimizer::contour(" << GetName()
00632                         << ") ERROR: " << var1.GetName()
00633                         << " is not a floating parameter of "
00634                         << _func->GetName() << endl ;
00635     return 0;
00636   }
00637 
00638   Int_t index2= _fcn->GetFloatParamList()->index(&var2);
00639   if(index2 < 0) {
00640     coutE(Minimization) << "RooMinimizer::contour(" << GetName()
00641                         << ") ERROR: " << var2.GetName()
00642                         << " is not a floating parameter of PDF "
00643                         << _func->GetName() << endl ;
00644     return 0;
00645   }
00646 
00647   // create and draw a frame
00648   RooPlot* frame = new RooPlot(var1,var2) ;
00649 
00650   // draw a point at the current parameter values
00651   TMarker *point= new TMarker(var1.getVal(), var2.getVal(), 8);
00652   frame->addObject(point) ;
00653 
00654   // remember our original value of ERRDEF
00655   Double_t errdef= _theFitter->Config().MinimizerOptions().ErrorDef();
00656 
00657   Double_t n[6] ;
00658   n[0] = n1 ; n[1] = n2 ; n[2] = n3 ; n[3] = n4 ; n[4] = n5 ; n[5] = n6 ;
00659   unsigned int npoints(50);
00660 
00661   for (Int_t ic = 0 ; ic<6 ; ic++) {
00662     if(n[ic] > 0) {
00663       // set the value corresponding to an n1-sigma contour
00664       _theFitter->Config().MinimizerOptions().SetErrorDef(n[ic]*n[ic]*errdef);
00665 
00666       // calculate and draw the contour
00667       Double_t *xcoor = new Double_t[npoints+1];
00668       Double_t *ycoor = new Double_t[npoints+1];
00669       bool ret = _theFitter->GetMinimizer()->Contour(index1,index2,npoints,xcoor,ycoor);
00670 
00671       if (!ret) {
00672         coutE(Minimization) << "RooMinimizer::contour("
00673                             << GetName()
00674                             << ") ERROR: MINUIT did not return a contour graph for n="
00675                             << n[ic] << endl ;
00676       } else {
00677         xcoor[npoints] = xcoor[0];
00678         ycoor[npoints] = ycoor[0];
00679         TGraph* graph = new TGraph(npoints+1,xcoor,ycoor);
00680 
00681         graph->SetName(Form("contour_%s_n%f",_func->GetName(),n[ic])) ;
00682         graph->SetLineStyle(ic+1) ;
00683         graph->SetLineWidth(2) ;
00684         graph->SetLineColor(kBlue) ;
00685         frame->addObject(graph,"L") ;
00686       }
00687 
00688       delete [] xcoor;
00689       delete [] ycoor;
00690     }
00691   }
00692 
00693 
00694   // restore the original ERRDEF
00695   _theFitter->Config().MinimizerOptions().SetErrorDef(errdef);
00696 
00697   // restore parameter values
00698   *params = *paramSave ;
00699   delete paramSave ;
00700 
00701   return frame ;
00702 
00703 }
00704 
00705 
00706 //_____________________________________________________________________________
00707 void RooMinimizer::profileStart()
00708 {
00709   // Start profiling timer
00710   if (_profile) {
00711     _timer.Start() ;
00712     _cumulTimer.Start(kFALSE) ;
00713   }
00714 }
00715 
00716 
00717 //_____________________________________________________________________________
00718 void RooMinimizer::profileStop()
00719 {
00720   // Stop profiling timer and report results of last session
00721   if (_profile) {
00722     _timer.Stop() ;
00723     _cumulTimer.Stop() ;
00724     coutI(Minimization) << "Command timer: " ; _timer.Print() ;
00725     coutI(Minimization) << "Session timer: " ; _cumulTimer.Print() ;
00726   }
00727 }
00728 
00729 
00730 
00731 
00732 
00733 //_____________________________________________________________________________
00734 void RooMinimizer::applyCovarianceMatrix(TMatrixDSym& V)
00735 {
00736   // Apply results of given external covariance matrix. i.e. propagate its errors
00737   // to all RRV parameter representations and give this matrix instead of the
00738   // HESSE matrix at the next save() call
00739 
00740   _extV = (TMatrixDSym*) V.Clone() ;
00741   _fcn->ApplyCovarianceMatrix(*_extV);
00742 
00743 }
00744 
00745 
00746 
00747 RooFitResult* RooMinimizer::lastMinuitFit(const RooArgList& varList)
00748 {
00749   // Import the results of the last fit performed, interpreting
00750   // the fit parameters as the given varList of parameters.
00751 
00752   if (_theFitter==0 || _theFitter->GetMinimizer()==0) {
00753     oocoutE((TObject*)0,InputArguments) << "RooMinimizer::save: Error, run minimization before!"
00754                                         << endl ;
00755     return 0;
00756   }
00757 
00758   // Verify length of supplied varList
00759   if (varList.getSize()>0 && varList.getSize()!=Int_t(_theFitter->Result().NTotalParameters())) {
00760     oocoutE((TObject*)0,InputArguments)
00761       << "RooMinimizer::lastMinuitFit: ERROR: supplied variable list must be either empty " << endl
00762       << "                             or match the number of variables of the last fit ("
00763       << _theFitter->Result().NTotalParameters() << ")" << endl ;
00764     return 0 ;
00765   }
00766 
00767 
00768   // Verify that all members of varList are of type RooRealVar
00769   TIterator* iter = varList.createIterator() ;
00770   RooAbsArg* arg  ;
00771   while((arg=(RooAbsArg*)iter->Next())) {
00772     if (!dynamic_cast<RooRealVar*>(arg)) {
00773       oocoutE((TObject*)0,InputArguments) << "RooMinimizer::lastMinuitFit: ERROR: variable '"
00774                                           << arg->GetName() << "' is not of type RooRealVar" << endl ;
00775       return 0 ;
00776     }
00777   }
00778   delete iter ;
00779 
00780   RooFitResult* res = new RooFitResult("lastMinuitFit","Last MINUIT fit") ;
00781 
00782   // Extract names of fit parameters
00783   // and construct corresponding RooRealVars
00784   RooArgList constPars("constPars") ;
00785   RooArgList floatPars("floatPars") ;
00786 
00787   UInt_t i ;
00788   for (i = 0; i < _theFitter->Result().NTotalParameters(); ++i) {
00789 
00790     TString varName(_theFitter->Result().GetParameterName(i));
00791     Bool_t isConst(_theFitter->Result().IsParameterFixed(i)) ;
00792 
00793     Double_t xlo = _theFitter->Config().ParSettings(i).LowerLimit();
00794     Double_t xhi = _theFitter->Config().ParSettings(i).UpperLimit();
00795     Double_t xerr = _theFitter->Result().Error(i);
00796     Double_t xval = _theFitter->Result().Value(i);
00797 
00798     RooRealVar* var ;
00799     if (varList.getSize()==0) {
00800 
00801       if ((xlo<xhi) && !isConst) {
00802         var = new RooRealVar(varName,varName,xval,xlo,xhi) ;
00803       } else {
00804         var = new RooRealVar(varName,varName,xval) ;
00805       }
00806       var->setConstant(isConst) ;
00807     } else {
00808 
00809       var = (RooRealVar*) varList.at(i)->Clone() ;
00810       var->setConstant(isConst) ;
00811       var->setVal(xval) ;
00812       if (xlo<xhi) {
00813         var->setRange(xlo,xhi) ;
00814       }
00815 
00816       if (varName.CompareTo(var->GetName())) {
00817         oocoutI((TObject*)0,Eval)  << "RooMinimizer::lastMinuitFit: fit parameter '" << varName
00818                                    << "' stored in variable '" << var->GetName() << "'" << endl ;
00819       }
00820 
00821     }
00822 
00823     if (isConst) {
00824       constPars.addOwned(*var) ;
00825     } else {
00826       var->setError(xerr) ;
00827       floatPars.addOwned(*var) ;
00828     }
00829   }
00830 
00831   res->setConstParList(constPars) ;
00832   res->setInitParList(floatPars) ;
00833   res->setFinalParList(floatPars) ;
00834   res->setMinNLL(_theFitter->Result().MinFcnValue()) ;
00835   res->setEDM(_theFitter->Result().Edm()) ;
00836   res->setCovQual(_theFitter->GetMinimizer()->CovMatrixStatus()) ;
00837   res->setStatus(_theFitter->Result().Status()) ;
00838   std::vector<double> globalCC;
00839   TMatrixDSym corrs(_theFitter->Result().Parameters().size()) ;
00840   TMatrixDSym covs(_theFitter->Result().Parameters().size()) ;
00841   for (UInt_t ic=0; ic<_theFitter->Result().Parameters().size(); ic++) {
00842     globalCC.push_back(_theFitter->Result().GlobalCC(ic));
00843     for (UInt_t ii=0; ii<_theFitter->Result().Parameters().size(); ii++) {
00844       corrs(ic,ii) = _theFitter->Result().Correlation(ic,ii);
00845       covs(ic,ii) = _theFitter->Result().CovMatrix(ic,ii);
00846     }
00847   }
00848   res->fillCorrMatrix(globalCC,corrs,covs) ;
00849 
00850   return res;
00851 
00852 }
00853 
00854 #endif

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