/* ======================================================================== *\ ! ! * ! * This file is part of MARS, the MAGIC Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appear in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Wolfgang Wittek 7/2003 ! ! Copyright: MAGIC Software Development, 2000-2003 ! ! \* ======================================================================== */ ///////////////////////////////////////////////////////////////////////////// // // // MMinuitInterface // // // // Class for interfacing with Minuit // // // // // // // ///////////////////////////////////////////////////////////////////////////// #include "MMinuitInterface.h" #include // fabs #include #include #include #include #include "MLog.h" #include "MLogManip.h" #include "MParContainer.h" ClassImp(MMinuitInterface); using namespace std; // -------------------------------------------------------------------------- // // Default constructor. // MMinuitInterface::MMinuitInterface(const char *name, const char *title) { fName = name ? name : "MMinuitInterface"; fTitle = title ? title : "Interface for Minuit"; } // ----------------------------------------------------------------------- // // Interface to MINUIT // // Bool_t MMinuitInterface::CallMinuit( void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t), 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) { 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) { minuit.FixParameter(i); step[i] = 0.0; } } //.............................................. //NumPars = minuit.GetNumPars(); //*fLog << "CallMinuit : number of free parameters = " // << NumPars << endl; //.............................................. // This doesn't seem to have any effect // Set maximum number of iterations (default = 500) //Int_t maxiter = 100000; //minuit.SetMaxIterations(maxiter); //.............................................. // 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) 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 (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); Int_t maxcalls = 3000; Double_t tolerance = 0.1; Double_t tmp[2]; tmp[0] = maxcalls; tmp[1] = tolerance; minuit.mnexcm("SIMPLEX", &tmp[0], 2, 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 (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, parabolic 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); // WW : the following statements were commented out because the // Minuit object will still be used; // this may be changed in the future //delete &minuit; //gMinuit = save; return kTRUE; }