TGo4FitModelFormula.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 "TGo4FitModelFormula.h"

#include <iostream>

#include "TFormula.h"
#include "TGo4FitParameter.h"

TGo4FitModelFormula::TGo4FitModelFormula()
   : TGo4FitModel(), fxExpression(), fxPosIndex(), fxWidthIndex(), fxFormula(nullptr)
{
}

TGo4FitModelFormula::TGo4FitModelFormula(const char *iName, const char *iExpressionStr, Int_t iNPars,
                                         Bool_t AddAmplitude)
   : TGo4FitModel(iName, "based on TFormula line shape", AddAmplitude), fxExpression(iExpressionStr), fxPosIndex(),
     fxWidthIndex(), fxFormula(nullptr)
{

   for (Int_t n = 0; n < iNPars; n++)
      NewParameter(GetExprParName(n), "formula parameter", 0.);
}

TGo4FitModelFormula::~TGo4FitModelFormula()
{
   CloseFormula();
}

Int_t TGo4FitModelFormula::GetNumberOfExprPar()
{
   Int_t num = TGo4FitParsList::NumPars();
   if (GetAmplPar())
      num--;
   return num;
}

TGo4FitParameter *TGo4FitModelFormula::GetExprPar(Int_t n)
{
   if ((n < 0) || (n >= GetNumberOfExprPar()))
      return nullptr;
   if ((GetAmplIndex() >= 0) && (n >= GetAmplIndex()))
      n++;
   return GetPar(n);
}

Bool_t TGo4FitModelFormula::SetNumberOfExprPar(Int_t num)
{
   Int_t oldnum = GetNumberOfExprPar();
   if ((num < 0) || (oldnum == num))
      return kFALSE;
   if (num < oldnum) {
      for (Int_t n = oldnum - 1; n >= num; n--)
         RemovePar(GetExprPar(n)->GetName());
   } else {
      for (Int_t n = oldnum; n < num; n++)
         NewParameter(GetExprParName(n), "formula parameter", 0.);
   }
   return kTRUE;
}

void TGo4FitModelFormula::SetPosParIndex(Int_t naxis, Int_t indx)
{
   if (naxis < 0)
      return;
   Int_t oldsize = fxPosIndex.GetSize();
   if (naxis >= oldsize) {
      fxPosIndex.Set(naxis + 1);
      for (Int_t n = oldsize; n < fxPosIndex.GetSize(); n++)
         fxPosIndex[n] = -1;
   }
   fxPosIndex[naxis] = indx;
}

void TGo4FitModelFormula::SetWidthParIndex(Int_t naxis, Int_t indx)
{
   if (naxis < 0)
      return;
   Int_t oldsize = fxWidthIndex.GetSize();
   if (naxis >= oldsize) {
      fxWidthIndex.Set(naxis + 1);
      for (Int_t n = oldsize; n < fxWidthIndex.GetSize(); n++)
         fxWidthIndex[n] = -1;
   }
   fxWidthIndex[naxis] = indx;
}

Int_t TGo4FitModelFormula::GetPosParIndex(Int_t naxis)
{
   if ((naxis < 0) || (naxis >= fxPosIndex.GetSize()))
      return -1;
   return GetParIndex(GetExprPar(fxPosIndex[naxis]));
}

Int_t TGo4FitModelFormula::GetWidthParIndex(Int_t naxis)
{
   if ((naxis < 0) || (naxis >= fxWidthIndex.GetSize()))
      return -1;
   return GetParIndex(GetExprPar(fxWidthIndex[naxis]));
}

TString TGo4FitModelFormula::GetExprParName(Int_t n)
{
   TString res;
   res.Form("Par%d", n);
   return res;
}

Bool_t TGo4FitModelFormula::Initialize(Int_t UseBuffers)
{
   if (!CompileFormula())
      return kFALSE;
   CloseFormula();
   return TGo4FitModel::Initialize(UseBuffers);
}

Bool_t TGo4FitModelFormula::BeforeEval(Int_t ndim)
{
   if (!CompileFormula())
      return kFALSE;
   Par_ndim = ndim;
   for (Int_t n = 0; n < NumPars(); n++)
      fxFormula->SetParameter(n, GetPar(n)->GetValue());
   return kTRUE;
}

Double_t TGo4FitModelFormula::EvalN(const Double_t *v)
{
   switch (Par_ndim) {
      case 0: return 0.;
      case 1: return fxFormula->Eval(v[0]);
      case 2: return fxFormula->Eval(v[0], v[1]);
      case 3: return fxFormula->Eval(v[0], v[1], v[2]);
      default: return fxFormula->EvalPar(v, nullptr);
   }
}

void TGo4FitModelFormula::AfterEval()
{
   CloseFormula();
}

void TGo4FitModelFormula::Finalize()
{
   CloseFormula();
   TGo4FitModel::Finalize();
}

void TGo4FitModelFormula::Print(Option_t *option) const
{
   TGo4FitModel::Print(option);
   std::cout << "  Expression = " << *fxExpression << std::endl;
   for (Int_t naxis = 0; naxis < fxPosIndex.GetSize(); naxis++) {
      TGo4FitParameter *par = ((TGo4FitModelFormula *)this)->GetExprPar(fxPosIndex[naxis]);
      if (par)
         std::cout << "  Position on " << naxis << " axis is " << par->GetName() << std::endl;
   }
   for (Int_t naxis = 0; naxis < fxWidthIndex.GetSize(); naxis++) {
      TGo4FitParameter *par = ((TGo4FitModelFormula *)this)->GetExprPar(fxWidthIndex[naxis]);
      if (par)
         std::cout << "  Width on " << naxis << " axis is " << par->GetName() << std::endl;
   }
}

Bool_t TGo4FitModelFormula::CompileFormula()
{
   CloseFormula();

   TString Expr(fxExpression);
   for (Int_t n = 0; n < NumPars(); n++) {
      TString str;
      str.Form("[%d]", n);
      Expr.ReplaceAll(GetParName(n), str);
   }

   fxFormula = new TFormula("Expression", Expr);

   Int_t err = fxFormula->Compile(Expr);
   if (err != 0) {
      std::cerr << "Error in formula: " << fxExpression.Data() << "     code " << err << std::endl;
      CloseFormula();
      return kFALSE;
   }
   return kTRUE;
}

void TGo4FitModelFormula::CloseFormula()
{
   if (fxFormula) {
      delete fxFormula;
      fxFormula = nullptr;
   }
}