TGo4FitMinuit.cxx

Go to the documentation of this file.

// $Id$
//-----------------------------------------------------------------------
//       The GSI Online Offline Object Oriented (Go4) Project
//         Experiment Data Processing at EE department, GSI
//-----------------------------------------------------------------------
// Copyright (C) 2000- GSI Helmholtzzentrum fuer Schwerionenforschung GmbH
//                     Planckstr. 1, 64291 Darmstadt, Germany
// Contact:            http://go4.gsi.de
//-----------------------------------------------------------------------
// This software can be used under the license agreements as stated
// in Go4License.txt file which is part of the distribution.
//-----------------------------------------------------------------------

#include "TGo4FitMinuit.h"

#include <iostream>

#include "TMinuit.h"
#include "TObjString.h"

#include "TGo4FitterAbstract.h"
#include "TGo4FitMinuitResult.h"

class TMinuitEx : public TMinuit {
public:
   TMinuitEx(Int_t NumPars, TGo4FitterAbstract *fitter);
   virtual ~TMinuitEx();

   Int_t Eval(Int_t npar, Double_t *grad, Double_t &fval, Double_t *pars, Int_t iflag) override;

protected:
   TGo4FitterAbstract *fxFitter{nullptr}; 
};

TMinuitEx::TMinuitEx(Int_t NumPars, TGo4FitterAbstract *fitter) : TMinuit(NumPars), fxFitter(fitter) {}

TMinuitEx::~TMinuitEx() {}

Int_t TMinuitEx::Eval(Int_t npar, Double_t *grad, Double_t &fval, Double_t *pars, Int_t iflag)
{
   if (fxFitter)
      fval = fxFitter->CalculateFitFunction(pars);
   else
      fval = 0.;
   return 0;
}

// *******************************************************************************

TGo4FitMinuit::TGo4FitMinuit() : TGo4FitterAction(), fxCommands(), fxResults() {}

TGo4FitMinuit::TGo4FitMinuit(const char *Name)
   : TGo4FitterAction(Name, "Fitter minimization using TMinuit object"), fxCommands(), fxResults()
{
   fxCommands.SetOwner(kTRUE);
   fxResults.SetOwner(kTRUE);
}

TGo4FitMinuit::~TGo4FitMinuit() {}

void TGo4FitMinuit::AddCommand(const char *iCommand)
{
   fxCommands.Add(new TObjString(iCommand));
}

const char *TGo4FitMinuit::GetCommand(Int_t n) const
{
   return fxCommands[n]->GetName();
}

void TGo4FitMinuit::DoAction(TGo4FitterAbstract *Fitter)
{

   TMinuitEx fMinuit(Fitter->NumPars(), Fitter);
   fMinuit.SetPrintLevel(-1);

   for (Int_t n = 0; n < Fitter->NumPars(); n++) {
      const char *FullName = Fitter->GetParFullName(n);
      Int_t ierflg = 0;
      Double_t epsilon = 0, RangeMin = 0, RangeMax = 0;
      if (!Fitter->GetParEpsilon(FullName, epsilon))
         epsilon = 1.;

      if (!Fitter->GetParRange(FullName, RangeMin, RangeMax)) {
         RangeMin = 0;
         RangeMax = 0;
      }

      fMinuit.mnparm(n, FullName, Fitter->GetParValue(FullName), epsilon, RangeMin, RangeMax, ierflg);

      if (Fitter->GetParFixed(FullName))
         fMinuit.FixParameter(n);
   }

   if (fxCommands.GetLast() < 0)
      fMinuit.Command("MIGRAD 500 1");
   else
      for (Int_t n = 0; n <= fxCommands.GetLast(); n++) {
         TString cmd = ((TObjString *)fxCommands[n])->GetString();
         if (cmd[0] == 'r') {
            if (cmd.Index("result", 6, 0, TString::kIgnoreCase) == 0) {
               cmd.Remove(0, 6);
               while ((cmd.Length() > 0) && (cmd[0] == ' '))
                  cmd.Remove(0, 1);
               if (cmd.IsNull())
                  cmd = "1000";
               if (cmd.Length() < 4) {
                  std::cerr << "invalid result command syntax" << std::endl;
                  break;
               }
               Bool_t getpar = (cmd[0] == '1');
               Bool_t geterr = (cmd[1] == '1');
               Bool_t getmatr = (cmd[2] == '1');
               Bool_t getcontr = (cmd[3] == '1');
               cmd.Remove(0, 4);
               while ((cmd.Length() > 0) && (cmd[0] == ' '))
                  cmd.Remove(0, 1);
               if (cmd.IsNull())
                  cmd = "Result";

               TGo4FitMinuitResult *res = new TGo4FitMinuitResult(cmd, "TMinuit result object");
               fxResults.Add(res);
               res->CallMNSTAT(&fMinuit);
               if (getpar)
                  res->CallMNPOUT(&fMinuit, Fitter->NumPars());
               if (geterr)
                  res->CallMNERRS(&fMinuit, Fitter->NumPars());
               if (getmatr)
                  res->CallMNEMAT(&fMinuit, Fitter->NumPars(), kTRUE);
               if (getcontr)
                  res->GetContourPlot(&fMinuit);
            }

         } else
            fMinuit.Command(cmd);
      }

   for (Int_t n = 0; n < Fitter->NumPars(); n++) {
      Double_t value, error;
      fMinuit.GetParameter(n, value, error);
      Fitter->SetParValue(Fitter->GetParFullName(n), value);
      Fitter->SetParError(Fitter->GetParFullName(n), error);
   }
}

TGo4FitMinuitResult *TGo4FitMinuit::GetResult(Int_t indx) const
{
   return (indx >= 0) && (indx <= fxResults.GetLast()) ? (TGo4FitMinuitResult *)fxResults.At(indx) : nullptr;
}

TGo4FitMinuitResult *TGo4FitMinuit::FindResult(const char *ResName)
{
   return (TGo4FitMinuitResult *)fxResults.FindObject(ResName);
}

void TGo4FitMinuit::AddResult(TGo4FitMinuitResult *res)
{
   fxResults.Add(res);
}

void TGo4FitMinuit::RemoveResult(TGo4FitMinuitResult *res)
{
   fxResults.Remove(res);
   fxResults.Compress();
}

void TGo4FitMinuit::Print(Option_t *option) const
{
   TGo4FitterAction::Print(option);
   if (fxCommands.GetLast() >= 0)
      std::cout << "List of commands:" << std::endl;
   for (Int_t n = 0; n <= fxCommands.GetLast(); n++)
      std::cout << "   " << fxCommands[n]->GetName() << std::endl;
   if (fxResults.GetLast() >= 0) {
      std::cout << "List of stored results:" << std::endl;
      fxResults.Print(option);
   }
}