Index: /trunk/MagicSoft/Mars/Changelog
===================================================================
--- /trunk/MagicSoft/Mars/Changelog	(revision 2728)
+++ /trunk/MagicSoft/Mars/Changelog	(revision 2729)
@@ -82,4 +82,10 @@
      - set class verseion to 0 -- not ment for writing at the
        moment
+
+   * mtools/MTMinui.[h,cc]:
+     - added (will replace MMinuitInterface soon)
+
+   * mtools/ToolsLinkDef.h, mtools/Makefile:
+     - added MTMinuit
 
 
Index: /trunk/MagicSoft/Mars/mtools/MTMinuit.cc
===================================================================
--- /trunk/MagicSoft/Mars/mtools/MTMinuit.cc	(revision 2729)
+++ /trunk/MagicSoft/Mars/mtools/MTMinuit.cc	(revision 2729)
@@ -0,0 +1,329 @@
+/* ======================================================================== *\
+!
+! *
+! * 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 <mailto:wittek@mppmu.mpg.de>
+!   Author(s): Thomas Bretz, 12/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
+!
+!   Copyright: MAGIC Software Development, 2000-2003
+!
+!
+\* ======================================================================== */
+
+/////////////////////////////////////////////////////////////////////////////
+//
+// MTMinuit
+//
+// Class for interfacing with Minuit
+//
+/////////////////////////////////////////////////////////////////////////////
+#include "MTMinuit.h"
+
+#include <TMinuit.h>
+
+#include "MLog.h"
+#include "MLogManip.h"
+
+ClassImp(MTMinuit);
+
+using namespace std;
+
+// --------------------------------------------------------------------------
+//
+// Default constructor.
+//
+MTMinuit::MTMinuit(const char *name, const char *title)
+    : fFcn(NULL), fNames(NULL), fMinuit(NULL)
+{
+    fName  = name  ? name  : "MTMinuit";
+    fTitle = title ? title : "Interface for Minuit";
+}
+
+MTMinuit::~MTMinuit()
+{
+    if (fMinuit)
+        delete fMinuit;
+}
+
+TMinuit *MTMinuit::GetMinuit() const
+{
+    return fMinuit;
+}
+
+void MTMinuit::SetFcn(void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t))
+{
+    fFcn = fcn;
+}
+
+void MTMinuit::InitParameters(const TArrayD &vinit, const TArrayD *step, const TString *name)
+{
+    const int n = fVinit.GetSize();
+
+    fVinit = vinit;
+
+    fStep.Set(n);
+    memset(fStep.GetArray(), 0, n*sizeof(Double_t));
+
+    if (step && step->GetSize()!=n)
+    {
+        *fLog << warn << "MTMinuit::SetParameters: number of steps doesn't match number of parameters... ignored." << endl;
+        step = NULL;
+    }
+
+    for (int i=0; i<n; i++)
+        fStep[i] = step ? (*step)[i] : TMath::Abs(fVinit[i]/10.0);
+
+    if (fNames)
+        delete fNames;
+
+    fNames = new TString[n];
+
+    // FIXME: Sanity check for TString array size!
+    for (int i=0; i<n ; i++)
+        fNames[i] = name ? (const char*)name[i] : Form("par%02d", i);
+
+    // ---------------------------------------------------------------
+
+    fLimLo.Set(n);
+    fLimUp.Set(n);
+    memset(fLimLo.GetArray(), 0, n*sizeof(Double_t));
+    memset(fLimUp.GetArray(), 0, n*sizeof(Double_t));
+
+    fFix.Set(n);
+    memset(fFix.GetArray(), 0, n*sizeof(Char_t));
+}
+
+void MTMinuit::SetLimits(const TArrayD &limlo, const TArrayD &limup)
+{
+    if (fVinit.GetSize()==limlo.GetSize())
+        fLimLo = limlo;
+    else
+        *fLog << warn << "MTMinuit::SetLimits: size of limlo doesn't match number of parameters... ignored." << endl;
+
+    if (fVinit.GetSize()==limup.GetSize())
+        fLimUp = limup;
+    else
+        *fLog << warn << "MTMinuit::SetLimits: size of limup doesn't match number of parameters... ignored." << endl;
+}
+
+void MTMinuit::SetFixedParameters(const TArrayC &fix)
+{
+    if (fVinit.GetSize()==fix.GetSize())
+        fFix = fix;
+    else
+        *fLog << warn << "MTMinuit::SetFixedParameters: size of fix doesn't match number of parameters... ignored." << endl;
+}
+
+// -----------------------------------------------------------------------
+//
+// Interface to MINUIT
+//
+Bool_t MTMinuit::CallMinuit(TObject *objectfit , const TString &method, Bool_t nulloutput)
+{
+    if (!fFcn(NULL))
+    {
+        *fLog << err << "CallMinuit: ERROR - fFcn not set... abort." << endl;
+        return kFALSE;
+    }
+    if (!fNames)
+    {
+        *fLog << err << "CallMinuit: ERROR - Parameters not set... abort." << endl;
+        return kFALSE;
+    }
+
+    // Make a copy of fStep
+    TArrayD step(fStep);
+
+    // Save gMinuit
+    TMinuit *const save = gMinuit;
+
+    // Set the maximum number of parameters
+    if (fMinuit)
+        delete fMinuit;
+    fMinuit = new TMinuit(fVinit.GetSize());
+
+    //
+    // 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
+    //
+    fMinuit->SetPrintLevel(-1);
+
+    //..............................................
+    // Printout for warnings
+    //    SET WAR      print warnings
+    //    SET NOW      suppress warnings
+    Int_t errWarn;
+    Double_t tmpwar = 0;
+    fMinuit->mnexcm("SET NOW", &tmpwar, 0, errWarn);
+    //fMinuit->mnexcm("SET WAR", &tmpwar, 0, errWarn);
+
+    //..............................................
+    // Set the address of the minimization function
+    fMinuit->SetFCN(fFcn);
+
+    //..............................................
+    // Store address of object to be used in fcn
+    fMinuit->SetObjectFit(objectfit);
+
+    //..............................................
+    // Set starting values and step sizes for parameters
+    for (Int_t i=0; i<fVinit.GetSize(); i++)
+    {
+        if (!fMinuit->DefineParameter(i, fNames[i], fVinit[i], fStep[i], fLimLo[i], fLimUp[i]))
+            continue;
+
+        *fLog << err << "CallMinuit: Error in defining parameter " << fNames[i] << endl;
+        return kFALSE;
+    }
+
+    //
+    // 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
+    //
+    fMinuit->SetErrorDef(1.0);
+
+    // Int_t errMigrad;
+    // Double_t tmp = 0;
+    // fMinuit->mnexcm("MIGRAD", &tmp, 0, errMigrad);
+
+    // fix a parameter
+    for (Int_t i=0; i<fVinit.GetSize(); i++)
+    {
+        if (!fFix[i])
+            continue;
+
+        fMinuit->FixParameter(i);
+    }
+
+    //..............................................
+    // This doesn't seem to have any effect
+    // Set maximum number of iterations (default = 500)
+    //Int_t maxiter = 100000;
+    //fMinuit->SetMaxIterations(maxiter);
+
+    // minimization by the method of Migrad
+    Int_t fErrMinimize=0;
+    if (method.Contains("Migrad", TString::kIgnoreCase))
+    {
+        if (nulloutput)
+            fLog->SetNullOutput(kTRUE);
+        Double_t tmp = 0;
+        fMinuit->mnexcm("MIGRAD", &tmp, 0, fErrMinimize);
+        if (nulloutput)
+            fLog->SetNullOutput(kFALSE);
+    }
+
+    //..............................................
+    // same minimization as by Migrad
+    // but switches to the SIMPLEX method if MIGRAD fails to converge
+    if (method.Contains("Minimize", TString::kIgnoreCase))
+    {
+        Double_t tmp = 0;
+        fMinuit->mnexcm("MINIMIZE", &tmp, 0, fErrMinimize);
+    }
+
+    //..............................................
+    // minimization by the SIMPLEX method
+    if (method.Contains("Simplex", TString::kIgnoreCase))
+    {
+        Double_t tmp[2];
+        tmp[0] = 3000; // maxcalls;
+        tmp[1] = 0.1;  // tolerance;
+
+        if (nulloutput)
+            fLog->SetNullOutput(kTRUE);
+        fMinuit->mnexcm("SIMPLEX", &tmp[0], 2, fErrMinimize);
+        if (nulloutput)
+            fLog->SetNullOutput(kFALSE);
+    }
+
+    //..............................................
+    // 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)
+    Double_t fMin,   fEdm,   fErrdef;
+    Int_t    fNpari, fNparx, fIstat;
+    fMinuit->mnstat(fMin, fEdm, fErrdef, fNpari, fNparx, fIstat);
+
+    if (fErrMinimize != 0)
+    {
+        *fLog << err << "CallMinuit: Minimization failed with:" << endl;
+        *fLog << "  fMin = " << fMin << endl;
+        *fLog << "  fEdm = "  << fEdm << endl;
+        *fLog << "  fErrdef = "  << fErrdef << endl;
+        *fLog << "  fIstat = " << fIstat << endl;
+        *fLog << "  fErrMinimize = " << fErrMinimize << endl;
+        return kFALSE;
+    }
+
+    //..............................................
+    // minimization by the method of Migrad
+    if (method.Contains("Hesse", TString::kIgnoreCase))
+    {
+        Double_t tmp = 0;
+        fMinuit->mnexcm("HESSE", &tmp, 0, fErrMinimize);
+    }
+
+    //..............................................
+    // Minos error analysis
+    if (method.Contains("Minos", TString::kIgnoreCase))
+    {
+        Double_t tmp = 0;
+        fMinuit->mnexcm("MINOS", &tmp, 0, fErrMinimize);
+    }
+
+    //..............................................
+    // 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;
+    //fMinuit->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;
+    fMinuit->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 
+    gMinuit = save;
+
+    return kTRUE;
+}
Index: /trunk/MagicSoft/Mars/mtools/MTMinuit.h
===================================================================
--- /trunk/MagicSoft/Mars/mtools/MTMinuit.h	(revision 2729)
+++ /trunk/MagicSoft/Mars/mtools/MTMinuit.h	(revision 2729)
@@ -0,0 +1,51 @@
+#ifndef MARS_MTMinuit
+#define MARS_MTMinuit
+
+#ifndef MARS_MParContainer
+#include "MParContainer.h"
+#endif
+
+#ifndef ROOT_TArrayC
+#include <TArrayC.h>
+#endif
+
+#ifndef ROOT_TArrayD
+#include <TArrayD.h>
+#endif
+
+class TMinuit;
+
+class MTMinuit : public MParContainer
+{
+private:
+    // FIXME: Maybe we can use a TMinuit Object to store all this?
+    void (*fFcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t);
+
+    TString *fNames;
+    TMinuit *fMinuit;
+
+    TArrayD fVinit;
+    TArrayD fStep;
+    TArrayD fLimLo;
+    TArrayD fLimUp;
+    TArrayC fFix;
+
+public:
+    // FIXME: Use FCN as first argument...
+    MTMinuit(const char *name=NULL, const char *title=NULL);
+    ~MTMinuit();
+
+    Bool_t CallMinuit(TObject *fObjectFit, const TString &method, Bool_t nulloutput);
+
+    TMinuit *GetMinuit() const;
+
+    void SetFcn(void (*fcn)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t));
+    void InitParameters(const TArrayD &vinit, const TArrayD *step=0, const TString *name=0);
+
+    void SetLimits(const TArrayD &limlo, const TArrayD &limup);
+    void SetFixedParameters(const TArrayC &fix);
+
+    ClassDef(MTMinuit, 0) // Class for interfacing with Minuit
+};
+
+#endif
Index: /trunk/MagicSoft/Mars/mtools/Makefile
===================================================================
--- /trunk/MagicSoft/Mars/mtools/Makefile	(revision 2728)
+++ /trunk/MagicSoft/Mars/mtools/Makefile	(revision 2729)
@@ -31,4 +31,5 @@
 SRCFILES = MChisqEval.cc \
 	   MTFillMatrix.cc \
+           MTMinuit.cc \
 	   MagicReversi.cc \
 	   MagicSnake.cc \
Index: /trunk/MagicSoft/Mars/mtools/ToolsLinkDef.h
===================================================================
--- /trunk/MagicSoft/Mars/mtools/ToolsLinkDef.h	(revision 2728)
+++ /trunk/MagicSoft/Mars/mtools/ToolsLinkDef.h	(revision 2729)
@@ -8,4 +8,5 @@
 
 #pragma link C++ class MTFillMatrix+;
+#pragma link C++ class MTMinuit+;
 
 #pragma link C++ class MineSweeper+;
