/*****************************************************************************
 * Project: RooFit                                                           *
 * Package: RooFitCore                                                       *
 * @(#)root/roofitcore:$Id: RooMinuit.cxx 26410 2008-11-24 13:03:17Z wouter $
 * Authors:                                                                  *
 *   WV, Wouter Verkerke, UC Santa Barbara, verkerke@slac.stanford.edu       *
 *   DK, David Kirkby,    UC Irvine,         dkirkby@uci.edu                 *
 *                                                                           *
 * Copyright (c) 2000-2005, Regents of the University of California          *
 *                          and Stanford University. All rights reserved.    *
 *                                                                           *
 * Redistribution and use in source and binary forms,                        *
 * with or without modification, are permitted according to the terms        *
 * listed in LICENSE (http://roofit.sourceforge.net/license.txt)             *
 *****************************************************************************/

//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
// RooMinuit is a wrapper class around TFitter/TMinuit that
// provides a seamless interface between the MINUIT functionality
// and the native RooFit interface.
// <p>
// RooMinuit can minimize any RooAbsReal function with respect to
// its parameters. Usual choices for minimization are RooNLLVar
// and RooChi2Var
// <p>
// RooMinuit has methods corresponding to MINUIT functions like
// hesse(), migrad(), minos() etc. In each of these function calls
// the state of the MINUIT engine is synchronized with the state
// of the RooFit variables: any change in variables, change
// in the constant status etc is forwarded to MINUIT prior to
// execution of the MINUIT call. Afterwards the RooFit objects
// are resynchronized with the output state of MINUIT: changes
// parameter values, errors are propagated.
// <p>
// Various methods are available to control verbosity, profiling,
// automatic PDF optimization.
// END_HTML
//

#include "RooFit.h"
#include "Riostream.h"

#include "TClass.h"
#include <fstream>
#include <iomanip>
#include "TH1.h"
#include "TH2.h"
#include "TMarker.h"
#include "TGraph.h"
#include "TStopwatch.h"
#include "TFitter.h"
#include "TMinuit.h"
#include "TDirectory.h"
#include "RooMinuit.h"
#include "RooArgSet.h"
#include "RooArgList.h"
#include "RooAbsReal.h"
#include "RooAbsRealLValue.h"
#include "RooRealVar.h"
#include "RooFitResult.h"
#include "RooAbsPdf.h"
#include "RooSentinel.h"
#include "RooMsgService.h"
#include "RooPlot.h"



#if (__GNUC__==3&&__GNUC_MINOR__==2&&__GNUC_PATCHLEVEL__==3)
char* operator+( streampos&, char* );
#endif

ClassImp(RooMinuit) 
;

TVirtualFitter *RooMinuit::_theFitter = 0 ;



//_____________________________________________________________________________
void RooMinuit::cleanup()
{
  // Cleanup method called by atexit handler installed by RooSentinel
  // to delete all global heap objects when the program is terminated
  if (_theFitter) {
    delete _theFitter ;
    _theFitter =0 ;
  }
}



//_____________________________________________________________________________
RooMinuit::RooMinuit(RooAbsReal& function)
{
  // Construct MINUIT interface to given function. Function can be anything,
  // but is typically a -log(likelihood) implemented by RooNLLVar or a chi^2
  // (implemented by RooChi2Var). Other frequent use cases are a RooAddition
  // of a RooNLLVar plus a penalty or constraint term. This class propagates
  // all RooFit information (floating parameters, their values and errors)
  // to MINUIT before each MINUIT call and propagates all MINUIT information
  // back to the RooFit object at the end of each call (updated parameter
  // values, their (asymmetric errors) etc. The default MINUIT error level
  // for HESSE and MINOS error analysis is taken from the defaultErrorLevel()
  // value of the input function.

  RooSentinel::activate() ;

  // Store function reference
  _func = &function ;
  _logfile = 0 ;
  _optConst = kFALSE ;
  _verbose = kFALSE ;
  _profile = kFALSE ;
  _handleLocalErrors = kTRUE ;
  _printLevel = 1 ;
  _printEvalErrors = 10 ;
  _warnLevel = -999 ;
  _doEvalErrorWall = kTRUE ;

  // Examine parameter list
  RooArgSet* paramSet = function.getParameters(RooArgSet()) ;
  RooArgList paramList(*paramSet) ;
  delete paramSet ;

  _floatParamList = (RooArgList*) paramList.selectByAttrib("Constant",kFALSE) ; 
  if (_floatParamList->getSize()>1) {
    _floatParamList->sort() ;
  }
  _floatParamList->setName("floatParamList") ;

  _constParamList = (RooArgList*) paramList.selectByAttrib("Constant",kTRUE) ;
  if (_constParamList->getSize()>1) {
    _constParamList->sort() ;
  }
  _constParamList->setName("constParamList") ;

  // Remove all non-RooRealVar parameters from list (MINUIT cannot handle them)
  TIterator* pIter = _floatParamList->createIterator() ;
  RooAbsArg* arg ;
  while((arg=(RooAbsArg*)pIter->Next())) {
    if (!arg->IsA()->InheritsFrom(RooAbsRealLValue::Class())) {
      coutW(Minimization) << "RooMinuit::RooMinuit: removing parameter " << arg->GetName() 
			  << " from list because it is not of type RooRealVar" << endl ;
      _floatParamList->remove(*arg) ;
    }
  }
  _nPar      = _floatParamList->getSize() ;  
  delete pIter ;

  // Save snapshot of initial lists
  _initFloatParamList = (RooArgList*) _floatParamList->snapshot(kFALSE) ;
  _initConstParamList = (RooArgList*) _constParamList->snapshot(kFALSE) ;

  // Initialize MINUIT
  Int_t nPar= _floatParamList->getSize();
  if (_theFitter) delete _theFitter ;
  _theFitter = new TFitter(nPar*2) ; //WVE Kludge, nPar*2 works around TMinuit memory allocation bug
  _theFitter->SetObjectFit(this) ;

  // Shut up for now
  setPrintLevel(-1) ;
  _theFitter->Clear();
  
  // Tell MINUIT to use our global glue function
  _theFitter->SetFCN(RooMinuitGlue);

  // Use +0.5 for 1-sigma errors
  setErrorLevel(function.defaultErrorLevel()) ;

  // Declare our parameters to MINUIT
  synchronize(kFALSE) ;

  // Reset the *largest* negative log-likelihood value we have seen so far
  _maxFCN= -1e30 ;
  _numBadNLL = 0 ;

  // Now set default verbosity
  if (RooMsgService::instance().silentMode()) {
    setWarnLevel(-1) ;
    setPrintLevel(-1) ;
  } else {
    setWarnLevel(1) ;
    setPrintLevel(1) ;
  }
}



//_____________________________________________________________________________
RooMinuit::~RooMinuit() 
{
  // Destructor

  delete _floatParamList ;
  delete _initFloatParamList ;
  delete _constParamList ;
  delete _initConstParamList ;
}



//_____________________________________________________________________________
void RooMinuit::setStrategy(Int_t istrat) 
{
  // Change MINUIT strategy to istrat. Accepted codes
  // are 0,1,2 and represent MINUIT strategies for dealing
  // most efficiently with fast FCNs (0), expensive FCNs (2)
  // and 'intermediate' FCNs (1)

  Double_t stratArg(istrat) ;
  _theFitter->ExecuteCommand("SET STR",&stratArg,1) ;
}



//_____________________________________________________________________________
void RooMinuit::setErrorLevel(Double_t level)
{
  // Set the level for MINUIT error analysis to the given
  // value. This function overrides the default value
  // that is taken in the RooMinuit constructor from
  // the defaultErrorLevel() method of the input function
  _theFitter->ExecuteCommand("SET ERR",&level,1);
}



//_____________________________________________________________________________
void RooMinuit::setEps(Double_t eps)
{
  // Change MINUIT epsilon 
  _theFitter->ExecuteCommand("SET EPS",&eps,1) ;  
}



//_____________________________________________________________________________
RooFitResult* RooMinuit::fit(const char* options)
{
  // Parse traditional RooAbsPdf::fitTo driver options
  // 
  //  s - Run Hesse first to estimate initial step size
  //  m - Run Migrad only
  //  h - Run Hesse to estimate errors
  //  v - Verbose mode
  //  l - Log parameters after each Minuit steps to file
  //  t - Activate profile timer
  //  r - Save fit result
  //  0 - Run Migrad with strategy 0

  TString opts(options) ;
  opts.ToLower() ;
  
  // Initial configuration
  if (opts.Contains("v")) setVerbose(1) ;
  if (opts.Contains("t")) setProfile(1) ;
  if (opts.Contains("l")) setLogFile(Form("%s.log",_func->GetName())) ;
  if (opts.Contains("c")) optimizeConst(1) ;

  // Fitting steps
  if (opts.Contains("s")) hesse() ;
  if (opts.Contains("0")) setStrategy(0) ;
  migrad() ;
  if (opts.Contains("0")) setStrategy(1) ;
  if (opts.Contains("h")||!opts.Contains("m")) hesse() ;
  if (!opts.Contains("m")) minos() ;
  
  return (opts.Contains("r")) ? save() : 0 ; 
}



//_____________________________________________________________________________
Int_t RooMinuit::migrad() 
{
  // Execute MIGRAD. Changes in parameter values
  // and calculated errors are automatically
  // propagated back the RooRealVars representing
  // the floating parameters in the MINUIT operation

  Double_t arglist[2];
  arglist[0]= 500*_nPar; // maximum iterations
  arglist[1]= 1.0;       // tolerance

  synchronize(_verbose) ;
  profileStart() ;
  RooAbsReal::enableEvalErrorLogging(kTRUE) ;
  RooAbsReal::clearEvalErrorLog() ;  
  _status= _theFitter->ExecuteCommand("MIGRAD",arglist,2);
  RooAbsReal::enableEvalErrorLogging(kFALSE) ;
  profileStop() ;
  backProp() ;
  return _status ;
}



//_____________________________________________________________________________
Int_t RooMinuit::hesse() 
{
  // Execute HESSE. Changes in parameter values
  // and calculated errors are automatically
  // propagated back the RooRealVars representing
  // the floating parameters in the MINUIT operation


  Double_t arglist[2];
  arglist[0]= 500*_nPar; // maximum iterations

  synchronize(_verbose) ;
  profileStart() ;
  RooAbsReal::enableEvalErrorLogging(kTRUE) ;
  RooAbsReal::clearEvalErrorLog() ;  
  _status= _theFitter->ExecuteCommand("HESSE",arglist,1);
  RooAbsReal::enableEvalErrorLogging(kFALSE) ;
  profileStop() ;
  backProp() ;
  return _status ;
}



//_____________________________________________________________________________
Int_t RooMinuit::minos() 
{
  // Execute MINOS. Changes in parameter values
  // and calculated errors are automatically
  // propagated back the RooRealVars representing
  // the floating parameters in the MINUIT operation

  Double_t arglist[2];
  arglist[0]= 500*_nPar; // maximum iterations

  synchronize(_verbose) ;
  profileStart() ;
  RooAbsReal::enableEvalErrorLogging(kTRUE) ;
  RooAbsReal::clearEvalErrorLog() ;  
  _status= _theFitter->ExecuteCommand("MINOS",arglist,1);
  RooAbsReal::enableEvalErrorLogging(kFALSE) ;
  profileStop() ;
  backProp() ;
  return _status ;
}


// added FMV, 08/18/03 

//_____________________________________________________________________________
Int_t RooMinuit::minos(const RooArgSet& minosParamList) 
{
  // Execute MINOS for given list of parameters. Changes in parameter values
  // and calculated errors are automatically
  // propagated back the RooRealVars representing
  // the floating parameters in the MINUIT operation

  Int_t nMinosPar(0) ;
  Double_t* arglist = new Double_t[_nPar+1];

  if (minosParamList.getSize()>0) {
    TIterator* aIter = minosParamList.createIterator() ;
    RooAbsArg* arg ;
    while((arg=(RooAbsArg*)aIter->Next())) {
      RooAbsArg* par = _floatParamList->find(arg->GetName());
      if (par && !par->isConstant()) {
	Int_t index = _floatParamList->index(par);
	nMinosPar++;
        arglist[nMinosPar]=index+1;
      }
    }
    delete aIter ;
  }
  arglist[0]= 500*_nPar; // maximum iterations

  synchronize(_verbose) ;
  profileStart() ;
  RooAbsReal::enableEvalErrorLogging(kTRUE) ;
  RooAbsReal::clearEvalErrorLog() ;  
  _status= _theFitter->ExecuteCommand("MINOS",arglist,1+nMinosPar);
  RooAbsReal::enableEvalErrorLogging(kFALSE) ;
  profileStop() ;
  backProp() ;

  delete[] arglist ;
  return _status ;
}



//_____________________________________________________________________________
Int_t RooMinuit::seek() 
{
  // Execute SEEK. Changes in parameter values
  // and calculated errors are automatically
  // propagated back the RooRealVars representing
  // the floating parameters in the MINUIT operation

  Double_t arglist[2];
  arglist[0]= 500*_nPar; // maximum iterations

  synchronize(_verbose) ;
  profileStart() ;
  RooAbsReal::enableEvalErrorLogging(kTRUE) ;
  RooAbsReal::clearEvalErrorLog() ;  
  _status= _theFitter->ExecuteCommand("SEEK",arglist,1);
  RooAbsReal::enableEvalErrorLogging(kFALSE) ;
  profileStop() ;
  backProp() ;
  return _status ;
}



//_____________________________________________________________________________
Int_t RooMinuit::simplex() 
{
  // Execute SIMPLEX. Changes in parameter values
  // and calculated errors are automatically
  // propagated back the RooRealVars representing
  // the floating parameters in the MINUIT operation

  Double_t arglist[2];
  arglist[0]= 500*_nPar; // maximum iterations
  arglist[1]= 1.0;       // tolerance

  synchronize(_verbose) ;
  profileStart() ;
  RooAbsReal::enableEvalErrorLogging(kTRUE) ;
  RooAbsReal::clearEvalErrorLog() ;  
  _status= _theFitter->ExecuteCommand("SIMPLEX",arglist,2);
  RooAbsReal::enableEvalErrorLogging(kFALSE) ;
  profileStop() ;
  backProp() ;
  return _status ;
}



//_____________________________________________________________________________
Int_t RooMinuit::improve()
{
  // Execute IMPROVE. Changes in parameter values
  // and calculated errors are automatically
  // propagated back the RooRealVars representing
  // the floating parameters in the MINUIT operation

  Double_t arglist[2];
  arglist[0]= 500*_nPar; // maximum iterations

  synchronize(_verbose) ;
  profileStart() ;
  RooAbsReal::enableEvalErrorLogging(kTRUE) ;
  RooAbsReal::clearEvalErrorLog() ;  
  _status= _theFitter->ExecuteCommand("IMPROVE",arglist,1);
  RooAbsReal::enableEvalErrorLogging(kFALSE) ;
  profileStop() ;
  backProp() ;
  return _status ;
}



//_____________________________________________________________________________
Int_t RooMinuit::setPrintLevel(Int_t newLevel) 
{
  // Change the MINUIT internal printing level
  Int_t ret = _printLevel ;
  Double_t arg(newLevel) ;
  _theFitter->ExecuteCommand("SET PRINT",&arg,1);
  _printLevel = newLevel ;
  return ret ;
}



//_____________________________________________________________________________
void RooMinuit::setNoWarn() 
{
  // Instruct MINUIT to suppress warnings

  Double_t arg(0) ;
  _theFitter->ExecuteCommand("SET NOWARNINGS",&arg,1);
  _warnLevel = -1 ;
}



//_____________________________________________________________________________
Int_t RooMinuit::setWarnLevel(Int_t newLevel) 
{
  // Set MINUIT warning level to given level

  if (newLevel==_warnLevel) {
    return _warnLevel ;
  }

  Int_t ret = _warnLevel ;
  Double_t arg(newLevel) ;

  if (newLevel>=0) {
    _theFitter->ExecuteCommand("SET WARNINGS",&arg,1);
  } else {
    Double_t arg2(0) ;
    _theFitter->ExecuteCommand("SET NOWARNINGS",&arg2,1);    
  }
  _warnLevel = newLevel ;
  
  return ret ;
}
      


//_____________________________________________________________________________
Bool_t RooMinuit::synchronize(Bool_t verbose)
{
  // Internal function to synchronize TMinuit with current
  // information in RooAbsReal function parameters
  
  Int_t oldPrint = setPrintLevel(-1) ;
  Int_t oldWarn = setWarnLevel(-1) ;
  Bool_t constValChange(kFALSE) ;
  Bool_t constStatChange(kFALSE) ;

  Int_t index(0) ;

  // Handle eventual migrations from constParamList -> floatParamList 
  for(index= 0; index < _constParamList->getSize() ; index++) {
    RooRealVar *par= dynamic_cast<RooRealVar*>(_constParamList->at(index)) ;
    if (!par) continue ;

    RooRealVar *oldpar= dynamic_cast<RooRealVar*>(_initConstParamList->at(index)) ;

    // Test if constness changed
    if (!par->isConstant()) {
      
      // Remove from constList, add to floatList
      _constParamList->remove(*par) ;
      _floatParamList->add(*par) ;
      _initFloatParamList->addClone(*oldpar) ;      
      _initConstParamList->remove(*oldpar) ;
      constStatChange=kTRUE ;
      _nPar++ ;

      if (verbose) {
	coutI(Minimization) << "RooMinuit::synchronize: parameter " << par->GetName() << " is now floating." << endl ;
      }
    } 

    // Test if value changed
    if (par->getVal()!= oldpar->getVal()) {
      constValChange=kTRUE ;      
      if (verbose) {
	coutI(Minimization) << "RooMinuit::synchronize: value of constant parameter " << par->GetName() 
			    << " changed from " << oldpar->getVal() << " to " << par->getVal() << endl ;
      }      
    }    

  }

  // Update reference list
  *_initConstParamList = *_constParamList ;


  // Synchronize MINUIT with function state
  for(index= 0; index < _nPar; index++) {
    RooRealVar *par= dynamic_cast<RooRealVar*>(_floatParamList->at(index)) ;

    Double_t pstep(0) ;
    Double_t pmin(0) ;
    Double_t pmax(0) ;

    if(!par->isConstant()) {

      // Verify that floating parameter is indeed of type RooRealVar 
      if (!par->IsA()->InheritsFrom(RooRealVar::Class())) {
	coutW(Minimization) << "RooMinuit::fit: Error, non-constant parameter " << par->GetName() 
			    << " is not of type RooRealVar, skipping" << endl ;
	continue ;
      }

      // Set the limits, if not infinite
      if (par->hasMin() && par->hasMax()) {
	pmin = par->getMin();
	pmax = par->getMax();
      }
      
      // Calculate step size
      pstep= par->getError();
      if(pstep <= 0) {
	// Floating parameter without error estitimate
	if (par->hasMin() && par->hasMax()) {
	  pstep= 0.1*(pmax-pmin);

	  // Trim default choice of error if within 2 sigma of limit
	  if (pmax - par->getVal() < 2*pstep) {
	    pstep = (pmax - par->getVal())/2 ;
	  } else if (par->getVal() - pmin < 2*pstep) {
	    pstep = (par->getVal() - pmin )/2 ;	    
	  }	  

	  // If trimming results in zero error, restore default
	  if (pstep==0) {
	    pstep= 0.1*(pmax-pmin);
	  }

	} else {
	  pstep=1 ;
	}						  
	if(_verbose) {
	  coutW(Minimization) << "RooMinuit::synchronize: WARNING: no initial error estimate available for "
			      << par->GetName() << ": using " << pstep << endl;
	}
      }       
    } else {
      pmin = par->getVal() ;
      pmax = par->getVal() ;      
    }

    // Extract previous information
    Double_t oldVar,oldVerr,oldVlo,oldVhi ;
    char oldParname[100] ;
    Int_t ierr = _theFitter->GetParameter(index,oldParname,oldVar,oldVerr,oldVlo,oldVhi)  ;

    // Determine if parameters is currently fixed in MINUIT

    Int_t ix ;
    Bool_t oldFixed(kFALSE) ;
    if (ierr>=0) {
      for (ix = 1; ix <= gMinuit->fNpfix; ++ix) {
        if (gMinuit->fIpfix[ix-1] == index+1) oldFixed=kTRUE ;
      }
    }
    
    if (par->isConstant() && !oldFixed) {

      // Parameter changes floating -> constant : update only value if necessary
      if (oldVar!=par->getVal()) {
	Double_t arglist[2] ;
	arglist[0] = index+1 ;
	arglist[1] = par->getVal() ;
	_theFitter->ExecuteCommand("SET PAR",arglist,2) ;
	if (verbose) {
	  coutI(Minimization) << "RooMinuit::synchronize: value of parameter " << par->GetName() << " changed from " << oldVar << " to " << par->getVal() << endl ;
	}
      }

      _theFitter->FixParameter(index) ;
      constStatChange=kTRUE ;
      if (verbose) {
	coutI(Minimization) << "RooMinuit::synchronize: parameter " << par->GetName() << " is now fixed." << endl ;
      }

    } else if (par->isConstant() && oldFixed) {
      
      // Parameter changes constant -> constant : update only value if necessary
      if (oldVar!=par->getVal()) {
	Double_t arglist[2] ;
	arglist[0] = index+1 ;
	arglist[1] = par->getVal() ;
	_theFitter->ExecuteCommand("SET PAR",arglist,2) ;
	constValChange=kTRUE ;

	if (verbose) {
	  coutI(Minimization) << "RooMinuit::synchronize: value of fixed parameter " << par->GetName() << " changed from " << oldVar << " to " << par->getVal() << endl ;
	}
      }

    } else {
      
      if (!par->isConstant() && oldFixed) {
	_theFitter->ReleaseParameter(index) ;
	constStatChange=kTRUE ;
	
	if (verbose) {
	  coutI(Minimization) << "RooMinuit::synchronize: parameter " << par->GetName() << " is now floating." << endl ;
	}
      } 

      // Parameter changes constant -> floating : update all if necessary
      if (oldVar!=par->getVal() || oldVlo!=pmin || oldVhi != pmax || oldVerr!=pstep) {
	_theFitter->SetParameter(index, par->GetName(), par->getVal(), pstep, pmin, pmax);
      }

      // Inform user about changes in verbose mode
      if (verbose && ierr>=0) {
	// if ierr<0, par was moved from the const list and a message was already printed

	if (oldVar!=par->getVal()) {
	  coutI(Minimization) << "RooMinuit::synchronize: value of parameter " << par->GetName() << " changed from " << oldVar << " to " << par->getVal() << endl ;
	}
	if (oldVlo!=pmin || oldVhi!=pmax) {
	  coutI(Minimization) << "RooMinuit::synchronize: limits of parameter " << par->GetName() << " changed from [" << oldVlo << "," << oldVhi 
	       << "] to [" << pmin << "," << pmax << "]" << endl ;
	}

	// If oldVerr=0, then parameter was previously fixed
	if (oldVerr!=pstep && oldVerr!=0) {
	coutI(Minimization) << "RooMinuit::synchronize: error/step size of parameter " << par->GetName() << " changed from " << oldVerr << " to " << pstep << endl ;
	}
      }      
    }
  }




  oldWarn = setWarnLevel(oldWarn) ;
  oldPrint = setPrintLevel(oldPrint) ;

  if (_optConst) {
    if (constStatChange) {

      RooAbsReal::enableEvalErrorLogging(kTRUE) ;

      coutI(Minimization) << "RooMinuit::synchronize: set of constant parameters changed, rerunning const optimizer" << endl ;
      _func->constOptimizeTestStatistic(RooAbsArg::ConfigChange) ;
    } else if (constValChange) {
      coutI(Minimization) << "RooMinuit::synchronize: constant parameter values changed, rerunning const optimizer" << endl ;
      _func->constOptimizeTestStatistic(RooAbsArg::ValueChange) ;
    }

    RooAbsReal::enableEvalErrorLogging(kFALSE) ;  

  }

  return 0 ;  
}
	



//_____________________________________________________________________________
void RooMinuit::optimizeConst(Bool_t flag) 
{
  // If flag is true, perform constant term optimization on
  // function being minimized.

  RooAbsReal::enableEvalErrorLogging(kTRUE) ;

  if (_optConst && !flag){ 
    if (_printLevel>-1) coutI(Minimization) << "RooMinuit::optimizeConst: deactivating const optimization" << endl ;
    _func->constOptimizeTestStatistic(RooAbsArg::DeActivate) ;
    _optConst = flag ;
  } else if (!_optConst && flag) {
    if (_printLevel>-1) coutI(Minimization) << "RooMinuit::optimizeConst: activating const optimization" << endl ;
    _func->constOptimizeTestStatistic(RooAbsArg::Activate) ;
    _optConst = flag ;
  } else if (_optConst && flag) {
    if (_printLevel>-1) coutI(Minimization) << "RooMinuit::optimizeConst: const optimization already active" << endl ;
  } else {
    if (_printLevel>-1) coutI(Minimization) << "RooMinuit::optimizeConst: const optimization wasn't active" << endl ;
  }

  RooAbsReal::enableEvalErrorLogging(kFALSE) ;

}



//_____________________________________________________________________________
RooFitResult* RooMinuit::save(const char* userName, const char* userTitle) 
{
  // Save and return a RooFitResult snaphot of current minimizer status.
  // This snapshot contains the values of all constant parameters,
  // the value of all floating parameters at RooMinuit construction and
  // after the last MINUIT operation, the MINUIT status, variance quality,
  // EDM setting, number of calls with evaluation problems, the minimized
  // function value and the full correlation matrix

  TString name,title ;
  name = userName ? userName : Form(_func->GetName()) ;
  title = userTitle ? userTitle : Form(_func->GetTitle()) ;  
  RooFitResult* fitRes = new RooFitResult(name,title) ;

  // Move eventual fixed paramaters in floatList to constList
  Int_t i ;
  RooArgList saveConstList(*_constParamList) ;
  RooArgList saveFloatInitList(*_initFloatParamList) ;
  RooArgList saveFloatFinalList(*_floatParamList) ;
  for (i=0 ; i<_floatParamList->getSize() ; i++) {
    RooAbsArg* par = _floatParamList->at(i) ;
    if (par->isConstant()) {
      saveFloatInitList.remove(*saveFloatInitList.find(par->GetName()),kTRUE) ;
      saveFloatFinalList.remove(*par) ;
      saveConstList.add(*par) ;
    }
  }
  saveConstList.sort() ;
  
  fitRes->setConstParList(saveConstList) ;
  fitRes->setInitParList(saveFloatInitList) ;
  
  Double_t edm, errdef, minVal;
  Int_t nvpar, nparx;
  Int_t icode = _theFitter->GetStats(minVal, edm, errdef, nvpar, nparx);
  fitRes->setStatus(_status) ;
  fitRes->setCovQual(icode) ;
  fitRes->setMinNLL(minVal) ;
  fitRes->setNumInvalidNLL(_numBadNLL) ;
  fitRes->setEDM(edm) ;    
  fitRes->setFinalParList(saveFloatFinalList) ;
  fitRes->fillCorrMatrix() ;

  return fitRes ;
}




//_____________________________________________________________________________
RooPlot* RooMinuit::contour(RooRealVar& var1, RooRealVar& var2, Double_t n1, Double_t n2, Double_t n3, Double_t n4, Double_t n5, Double_t n6) 
{
  // Create and draw a TH2 with the error contours in parameters var1 and v2 at up to 6 'sigma' settings
  // where 'sigma' is calculated as n*n*errorLevel

  // Verify that both variables are floating parameters of PDF
  Int_t index1= _floatParamList->index(&var1);
  if(index1 < 0) {
    coutE(Minimization) << "RooMinuit::contour(" << GetName() 
			<< ") ERROR: " << var1.GetName() << " is not a floating parameter of " << _func->GetName() << endl ;
    return 0;
  }

  Int_t index2= _floatParamList->index(&var2);
  if(index2 < 0) {
    coutE(Minimization) << "RooMinuit::contour(" << GetName() 
			<< ") ERROR: " << var2.GetName() << " is not a floating parameter of PDF " << _func->GetName() << endl ;
    return 0;
  }
  
  // create and draw a frame
  RooPlot* frame = new RooPlot(var1,var2) ;

  // draw a point at the current parameter values
  TMarker *point= new TMarker(var1.getVal(), var2.getVal(), 8);
  frame->addObject(point) ;

  // remember our original value of ERRDEF
  Double_t errdef= gMinuit->fUp;

  Double_t n[6] ;
  n[0] = n1 ; n[1] = n2 ; n[2] = n3 ; n[3] = n4 ; n[4] = n5 ; n[5] = n6 ;
  

  for (Int_t ic = 0 ; ic<6 ; ic++) {
    if(n[ic] > 0) {
      // set the value corresponding to an n1-sigma contour
      gMinuit->SetErrorDef(n[ic]*n[ic]*errdef);
      // calculate and draw the contour
      TGraph* graph= (TGraph*)gMinuit->Contour(50, index1, index2);
      if (!graph) {
	coutE(Minimization) << "RooMinuit::contour(" << GetName() << ") ERROR: MINUIT did not return a contour graph for n=" << n[ic] << endl ;
      } else {
	graph->SetName(Form("contour_%s_n%f",_func->GetName(),n[ic])) ;
	graph->SetLineStyle(ic+1) ;
	graph->SetLineWidth(2) ;
	graph->SetLineColor(kBlue) ;
	frame->addObject(graph,"L") ;
      }
    }
  }

  // restore the original ERRDEF
  gMinuit->SetErrorDef(errdef);
  
  return frame ;
}



//_____________________________________________________________________________
Bool_t RooMinuit::setLogFile(const char* inLogfile) 
{
  // Change the file name for logging of a RooMinuit of all MINUIT steppings
  // through the parameter space. If inLogfile is null, the current log file
  // is closed and logging is stopped.

  if (_logfile) {
    coutI(Minimization) << "RooMinuit::setLogFile: closing previous log file" << endl ;
    _logfile->close() ;
    delete _logfile ;
    _logfile = 0 ;
  }
  _logfile = new ofstream(inLogfile) ;
  if (!_logfile->good()) {
    coutI(Minimization) << "RooMinuit::setLogFile: cannot open file " << inLogfile << endl ;
    _logfile->close() ;
    delete _logfile ;
    _logfile= 0;
  }  
  return kFALSE ;
}



//_____________________________________________________________________________
Double_t RooMinuit::getPdfParamVal(Int_t index)
{
  // Access PDF parameter value by ordinal index (needed by MINUIT)

  return ((RooRealVar*)_floatParamList->at(index))->getVal() ;
}



//_____________________________________________________________________________
Double_t RooMinuit::getPdfParamErr(Int_t index)
{
  // Access PDF parameter error by ordinal index (needed by MINUIT)

  return ((RooRealVar*)_floatParamList->at(index))->getError() ;  
}



//_____________________________________________________________________________
Bool_t RooMinuit::setPdfParamVal(Int_t index, Double_t value, Bool_t verbose)
{
  // Modify PDF parameter value by ordinal index (needed by MINUIT)

  RooRealVar* par = (RooRealVar*)_floatParamList->at(index) ;

  if (par->getVal()!=value) {
    if (verbose) cout << par->GetName() << "=" << value << ", " ;
    par->setVal(value) ;  
    return kTRUE ;
  }

  return kFALSE ;
}



//_____________________________________________________________________________
void RooMinuit::setPdfParamErr(Int_t index, Double_t value)
{
  // Modify PDF parameter error by ordinal index (needed by MINUIT)

  ((RooRealVar*)_floatParamList->at(index))->setError(value) ;    
}



//_____________________________________________________________________________
void RooMinuit::clearPdfParamAsymErr(Int_t index) 
{
  // Modify PDF parameter error by ordinal index (needed by MINUIT)

  ((RooRealVar*)_floatParamList->at(index))->removeAsymError() ;      
}


//_____________________________________________________________________________
void RooMinuit::setPdfParamErr(Int_t index, Double_t loVal, Double_t hiVal) 
{
  // Modify PDF parameter error by ordinal index (needed by MINUIT)

  ((RooRealVar*)_floatParamList->at(index))->setAsymError(loVal,hiVal) ;    
}



//_____________________________________________________________________________
void RooMinuit::profileStart() 
{
  // Start profiling timer
  if (_profile) {
    _timer.Start() ;
    _cumulTimer.Start(kFALSE) ;
  }
}




//_____________________________________________________________________________
void RooMinuit::profileStop() 
{
  // Stop profiling timer and report results of last session
  if (_profile) {
    _timer.Stop() ;
    _cumulTimer.Stop() ;
    coutI(Minimization) << "Command timer: " ; _timer.Print() ;
    coutI(Minimization) << "Session timer: " ; _cumulTimer.Print() ;
  }
}





//_____________________________________________________________________________
void RooMinuit::backProp() 
{
  // Transfer MINUIT fit results back into RooFit objects

  Double_t val,err,vlo,vhi, eplus, eminus, eparab, globcc;
  char buffer[10240];
  Int_t index ;
  for(index= 0; index < _nPar; index++) {
    _theFitter->GetParameter(index, buffer, val, err, vlo, vhi);
    setPdfParamVal(index, val);
    _theFitter->GetErrors(index, eplus, eminus, eparab, globcc);

    // Set the parabolic error
    setPdfParamErr(index, err);

    if(eplus > 0 || eminus < 0) {
      // Store the asymmetric error, if it is available
      setPdfParamErr(index, eminus,eplus);
    } else {
      // Clear the asymmetric error
      clearPdfParamAsymErr(index) ;
    }
  }
}



void RooMinuitGlue(Int_t& /*np*/, Double_t* /*gin*/,
		   Double_t &f, Double_t *par, Int_t /*flag*/)
{
  // Static function that interfaces minuit with RooMinuit

  // Retrieve fit context and its components
  RooMinuit* context = (RooMinuit*) RooMinuit::_theFitter->GetObjectFit() ;
  ofstream* logf   = context->logfile() ;
  Double_t& maxFCN = context->maxFCN() ;
  Bool_t verbose   = context->_verbose ;

  // Set the parameter values for this iteration
  Int_t nPar= context->getNPar();
  for(Int_t index= 0; index < nPar; index++) {
    if (logf) (*logf) << par[index] << " " ;
    context->setPdfParamVal(index, par[index],verbose);
  }

  // Calculate the function for these parameters
  f= context->_func->getVal() ;
  if ( RooAbsPdf::evalError() || RooAbsReal::numEvalErrors()>0 ) {

    if (context->_printEvalErrors>=0) {

      if (context->_doEvalErrorWall) {
	oocoutW(context,Minimization) << "RooFitGlue: Minimized function has error status." << endl 
				      << "Returning maximum FCN so far (" << maxFCN 
				      << ") to force MIGRAD to back out of this region. Error log follows" << endl ;
      } else {
	oocoutW(context,Minimization) << "RooFitGlue: Minimized function has error status but is ignored" << endl ;
      }	

      TIterator* iter = context->_floatParamList->createIterator() ;
      RooRealVar* var ;
      Bool_t first(kTRUE) ;
      ooccoutW(context,Minimization) << "Parameter values: " ;
      while((var=(RooRealVar*)iter->Next())) {
	if (first) { first = kFALSE ; } else ooccoutW(context,Minimization) << ", " ;
	ooccoutW(context,Minimization) << var->GetName() << "=" << var->getVal() ;
      }
      delete iter ;
      ooccoutW(context,Minimization) << endl ;
      
      RooAbsReal::printEvalErrors(ooccoutW(context,Minimization),context->_printEvalErrors) ;
      ooccoutW(context,Minimization) << endl ;
    } 

    if (context->_doEvalErrorWall) {
      f = maxFCN ;
    }

    RooAbsPdf::clearEvalError() ;
    RooAbsReal::clearEvalErrorLog() ;
    context->_numBadNLL++ ;
  } else if (f>maxFCN) {
    maxFCN = f ;
  }

  // Optional logging
  if (logf) (*logf) << setprecision(15) << f << setprecision(4) << endl;
  if (verbose) {
    cout << "\nprevFCN = " << setprecision(10) << f << setprecision(4) << "  " ;
    cout.flush() ;
  }
}