Index: trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc (revision 2311) +++ trunk/MagicSoft/Mars/manalysis/MMinuitInterface.cc (revision 2313) @@ -36,4 +36,6 @@ #include // fabs +#include +#include #include #include @@ -42,5 +44,4 @@ #include "MLogManip.h" #include "MParContainer.h" - ClassImp(MMinuitInterface); @@ -57,15 +58,5 @@ fName = name ? name : "MMinuitInterface"; fTitle = title ? title : "Interface for Minuit"; - } - -// -------------------------------------------------------------------------- -// -// Default destructor. -// -MMinuitInterface::~MMinuitInterface() -{ -} - // ----------------------------------------------------------------------- @@ -76,250 +67,224 @@ Bool_t MMinuitInterface::CallMinuit( void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t), - UInt_t npar, char name[80][100], - Double_t vinit[80], Double_t step[80], - Double_t limlo[80], Double_t limup[80], Int_t fix[80], - TObject *objectfit , TString method, Bool_t nulloutput) + const TString* name, const TArrayD &vinit, const TArrayD &step2, + const TArrayD &limlo, const TArrayD &limup, const TArrayI &fix, + TObject *objectfit , const TString &method, Bool_t nulloutput) { - TString str(method); - - // - // Be carefull: This is not thread safe - // - UInt_t maxpar = 100; - - if (npar > maxpar) + TArrayD step(step2); + + // + // Be carefull: This is not thread safe + // + if (vinit.GetSize() != step.GetSize() || + vinit.GetSize() != limlo.GetSize() || + vinit.GetSize() != limup.GetSize() || + vinit.GetSize() != fix.GetSize()) + { + *fLog << "CallMinuit: Parameter Size Mismatch" << endl; + return kFALSE; + } + + + + //.............................................. + // Set the maximum number of parameters + TMinuit *const save = gMinuit; + + TMinuit &minuit = *new TMinuit(vinit.GetSize()); + minuit.SetName("MinuitSupercuts"); + + + //.............................................. + // Set the print level + // -1 no output except SHOW comands + // 0 minimum output + // 1 normal output (default) + // 2 additional ouput giving intermediate results + // 3 maximum output, showing progress of minimizations + // + minuit.SetPrintLevel(1); + + //.............................................. + // Printout for warnings + // SET WAR print warnings + // SET NOW suppress warnings + Int_t errWarn; + Double_t tmpwar = 0; + minuit.mnexcm("SET NOW", &tmpwar, 0, errWarn); + //minuit.mnexcm("SET WAR", &tmpwar, 0, errWarn); + + //.............................................. + // Set the address of the minimization function + minuit.SetFCN(fcn); + + //.............................................. + // Store address of object to be used in fcn + minuit.SetObjectFit(objectfit ); + + //.............................................. + // Set starting values and step sizes for parameters + for (Int_t i=0; i 0) { - if (minuit.DefineParameter(i, &name[i][0], vinit[i], step[i], - limlo[i], limup[i])) - { - *fLog << "CallMinuit: Error in defining parameter " - << name[i][0] << endl; - return kFALSE; - } + minuit.FixParameter(i); + step[i] = 0.0; } - - //.............................................. - //Int_t NumPars; - //NumPars = minuit.GetNumPars(); - //*fLog << "CallMinuit : number of free parameters = " - // << NumPars << endl; - - - //.............................................. - // Error definition : - // - // for chisquare function : - // up = 1.0 means calculate 1-standard deviation error - // = 4.0 means calculate 2-standard deviation error - // - // for log(likelihood) function : - // up = 0.5 means calculate 1-standard deviation error - // = 2.0 means calculate 2-standard deviation error - Double_t up = 1.0; - minuit.SetErrorDef(up); - - - - // Int_t errMigrad; - // Double_t tmp = 0; - // minuit.mnexcm("MIGRAD", &tmp, 0, errMigrad); - - - //.............................................. - // fix a parameter - for (UInt_t i=0; i 0) - { - Int_t parNo = i; - minuit.FixParameter(parNo); - step[parNo] = 0.0; - } - } - - //.............................................. - //NumPars = minuit.GetNumPars(); - //*fLog << "CallMinuit : number of free parameters = " - // << NumPars << endl; - - //.............................................. - // Set maximum number of iterations (default = 500) - //Int_t maxiter = 100000; - Int_t maxiter = 10; - minuit.SetMaxIterations(maxiter); - - //.............................................. - // minimization by the method of Migrad - if (str.Contains("Migrad", TString::kIgnoreCase)) - { - //*fLog << "call MIGRAD" << endl; - if (nulloutput) + } + + //.............................................. + //NumPars = minuit.GetNumPars(); + //*fLog << "CallMinuit : number of free parameters = " + // << NumPars << endl; + + //.............................................. + // Set maximum number of iterations (default = 500) + //Int_t maxiter = 100000; + minuit.SetMaxIterations(10); + + //.............................................. + // minimization by the method of Migrad + if (method.Contains("Migrad", TString::kIgnoreCase)) + { + //*fLog << "call MIGRAD" << endl; + if (nulloutput) fLog->SetNullOutput(kTRUE); - Double_t tmp = 0; - minuit.mnexcm("MIGRAD", &tmp, 0, fErrMinimize); - if (nulloutput) + Double_t tmp = 0; + minuit.mnexcm("MIGRAD", &tmp, 0, fErrMinimize); + if (nulloutput) fLog->SetNullOutput(kFALSE); - //*fLog << "return from MIGRAD" << endl; - } - - //.............................................. - // same minimization as by Migrad - // but switches to the SIMPLEX method if MIGRAD fails to converge - if (str.Contains("Minimize", TString::kIgnoreCase)) - { - *fLog << "call MINIMIZE" << endl; - Double_t tmp = 0; - minuit.mnexcm("MINIMIZE", &tmp, 0, fErrMinimize); - *fLog << "return from MINIMIZE" << endl; - } - - //.............................................. - // minimization by the SIMPLEX method - if (str.Contains("Simplex", TString::kIgnoreCase)) - { - *fLog << "call SIMPLEX" << endl; - if (nulloutput) + //*fLog << "return from MIGRAD" << endl; + } + + //.............................................. + // same minimization as by Migrad + // but switches to the SIMPLEX method if MIGRAD fails to converge + if (method.Contains("Minimize", TString::kIgnoreCase)) + { + *fLog << "call MINIMIZE" << endl; + Double_t tmp = 0; + minuit.mnexcm("MINIMIZE", &tmp, 0, fErrMinimize); + *fLog << "return from MINIMIZE" << endl; + } + + //.............................................. + // minimization by the SIMPLEX method + if (method.Contains("Simplex", TString::kIgnoreCase)) + { + *fLog << "call SIMPLEX" << endl; + if (nulloutput) fLog->SetNullOutput(kTRUE); - Double_t tmp = 0; - minuit.mnexcm("SIMPLEX", &tmp, 0, fErrMinimize); - if (nulloutput) + Double_t tmp = 0; + minuit.mnexcm("SIMPLEX", &tmp, 0, fErrMinimize); + if (nulloutput) fLog->SetNullOutput(kFALSE); - //*fLog << "return from SIMPLEX" << endl; - } - - //.............................................. - // check quality of minimization - // istat = 0 covariance matrix not calculated - // 1 diagonal approximation only (not accurate) - // 2 full matrix, but forced positive-definite - // 3 full accurate covariance matrix - // (indication of normal convergence) - minuit.mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat); - - //if (fErrMinimize != 0 || fIstat < 3) - if (fErrMinimize != 0) - { - *fLog << "CallMinuit : Minimization failed" << endl; - *fLog << " fMin = " << fMin << ", fEdm = " << fEdm - << ", fErrdef = " << fErrdef << ", fIstat = " << fIstat - << ", fErrMinimize = " << fErrMinimize << endl; - return kFALSE; - } - - //*fLog << "CallMinuit : Minimization was successful" << endl; - //*fLog << " fMin = " << fMin << ", fEdm = " << fEdm - // << ", fErrdef = " << fErrdef << ", fIstat = " << fIstat - // << ", fErrMinimize = " << fErrMinimize << endl; - - - //.............................................. - // minimization by the method of Migrad - if (str.Contains("Hesse", TString::kIgnoreCase)) - { - //*fLog << "call HESSE" << endl; - Double_t tmp = 0; - minuit.mnexcm("HESSE", &tmp, 0, fErrMinimize); - //*fLog << "return from HESSE" << endl; - } - - //.............................................. - // Minos error analysis - if (str.Contains("Minos", TString::kIgnoreCase)) - { - //*fLog << "call MINOS" << endl; - Double_t tmp = 0; - minuit.mnexcm("MINOS", &tmp, 0, fErrMinimize); - //*fLog << "return from MINOS" << endl; - } - - //.............................................. - // Print current status of minimization - // if nkode = 0 only function value - // 1 parameter values, errors, limits - // 2 values, errors, step sizes, internal values - // 3 values, errors, step sizes, 1st derivatives - // 4 values, paraboloc errors, MINOS errors - - //Int_t nkode = 4; - //minuit.mnprin(nkode, fmin); - - //.............................................. - // call fcn with IFLAG = 3 (final calculation : calculate p(chi2)) - // iflag = 1 initial calculations only - // 2 calculate 1st derivatives and function - // 3 calculate function only - // 4 calculate function + final calculations - const char *command = "CALL"; - Double_t iflag = 3; - Int_t errfcn3; - minuit.mnexcm(command, &iflag, 1, errfcn3); - - //*fLog << "Exit CallMinuit" << endl; - - return kTRUE; + //*fLog << "return from SIMPLEX" << endl; + } + + //.............................................. + // check quality of minimization + // istat = 0 covariance matrix not calculated + // 1 diagonal approximation only (not accurate) + // 2 full matrix, but forced positive-definite + // 3 full accurate covariance matrix + // (indication of normal convergence) + minuit.mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat); + + //if (fErrMinimize != 0 || fIstat < 3) + if (fErrMinimize != 0) + { + *fLog << "CallMinuit : Minimization failed" << endl; + *fLog << " fMin = " << fMin << ", fEdm = " << fEdm + << ", fErrdef = " << fErrdef << ", fIstat = " << fIstat + << ", fErrMinimize = " << fErrMinimize << endl; + return kFALSE; + } + + //*fLog << "CallMinuit : Minimization was successful" << endl; + //*fLog << " fMin = " << fMin << ", fEdm = " << fEdm + // << ", fErrdef = " << fErrdef << ", fIstat = " << fIstat + // << ", fErrMinimize = " << fErrMinimize << endl; + + + //.............................................. + // minimization by the method of Migrad + if (method.Contains("Hesse", TString::kIgnoreCase)) + { + //*fLog << "call HESSE" << endl; + Double_t tmp = 0; + minuit.mnexcm("HESSE", &tmp, 0, fErrMinimize); + //*fLog << "return from HESSE" << endl; + } + + //.............................................. + // Minos error analysis + if (method.Contains("Minos", TString::kIgnoreCase)) + { + //*fLog << "call MINOS" << endl; + Double_t tmp = 0; + minuit.mnexcm("MINOS", &tmp, 0, fErrMinimize); + //*fLog << "return from MINOS" << endl; + } + + //.............................................. + // Print current status of minimization + // if nkode = 0 only function value + // 1 parameter values, errors, limits + // 2 values, errors, step sizes, internal values + // 3 values, errors, step sizes, 1st derivatives + // 4 values, paraboloc errors, MINOS errors + + //Int_t nkode = 4; + //minuit.mnprin(nkode, fmin); + + //.............................................. + // call fcn with IFLAG = 3 (final calculation : calculate p(chi2)) + // iflag = 1 initial calculations only + // 2 calculate 1st derivatives and function + // 3 calculate function only + // 4 calculate function + final calculations + Double_t iflag = 3; + Int_t errfcn3; + minuit.mnexcm("CALL", &iflag, 1, errfcn3); + + delete &minuit; + gMinuit = save; + + return kTRUE; } - - -//=========================================================================== - - - - - - - - - - - Index: trunk/MagicSoft/Mars/manalysis/MMinuitInterface.h =================================================================== --- trunk/MagicSoft/Mars/manalysis/MMinuitInterface.h (revision 2311) +++ trunk/MagicSoft/Mars/manalysis/MMinuitInterface.h (revision 2313) @@ -2,36 +2,29 @@ #define MARS_MMinuitInterface -#ifndef MARS_MTask -#include "MTask.h" +#ifndef MARS_MParContainer +#include "MParContainer.h" #endif -#ifndef ROOT_TArrayD -#include -#endif +class TArrayD; +class TArrayI; - - -class MMinuitInterface : public MTask +class MMinuitInterface : public MParContainer { private: - - UInt_t fNpar; - Double_t fMin, fEdm, fErrdef; - Int_t fNpari, fNparx, fIstat; - Int_t fErrMinimize; - + UInt_t fNpar; + Double_t fMin, fEdm, fErrdef; + Int_t fNpari, fNparx, fIstat; + Int_t fErrMinimize; public: - MMinuitInterface(const char *name=NULL, const char *title=NULL); - ~MMinuitInterface(); + MMinuitInterface(const char *name=NULL, const char *title=NULL); -Bool_t MMinuitInterface::CallMinuit( - void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t), - UInt_t npar, char name[80][100], - Double_t vinit[80], Double_t step[80], - Double_t limlo[80], Double_t limup[80], Int_t fix[80], - TObject *fObjectFit, TString method, Bool_t nulloutput); + Bool_t CallMinuit( + void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t), + const TString *name, const TArrayD &vinit, const TArrayD &step, + const TArrayD &limlo, const TArrayD &limup, const TArrayI &fix, + TObject *fObjectFit, const TString &method, Bool_t nulloutput); - ClassDef(MMinuitInterface, 1) // Class for interfacing with Minuit + ClassDef(MMinuitInterface, 0) // Class for interfacing with Minuit };