Minimization for blastwave function returning nan results

Newbie
1

Hello all,

I am a grad student working on reproducing some plots using existing data. I am very new to ROOT/C++ and have very little background in coding, so this is all a learning experience for me. The macro I’ve built with the help of a postdoc runs without errors, however it returns nan results and prints out the guesses we input for the parameters. Any insights on to what could be going wrong is much appreciated. I have attached the data files.

Here is my code:

#include "TMinuit.h"
#include "TString.h"
#include "TDatime.h"
#include "TGraph.h"
#include "TAxis.h"
#include "TCanvas.h"

    const char *f1 = "K0s_mult130-inf_pp_7000GeV.txt";   //Data File Array
    const char *f2 = "Lambda_mult130-inf_pp_7000GeV.txt";

int linenumber = 0;
Float_t x1ks[30], y1ks[30], ey1ks[30];
Float_t x1la[30], y1la[30], ey1la[30];
Float_t x2ks[30], y2ks[30], ey2ks[30];
Float_t x2la[30], y2la[30], ey2la[30];
Float_t x3ks[30], y3ks[30], ey3ks[30];
Float_t x3la[30], y3la[30], ey3la[30];


//______________________________________________________________________________
double getChi2;
int getNdf;

const double R = 1.0;
const double dr = 1e-3; // FIXME

double integral(double beta_T, double T, double n, double pt, double mt)
{
  
  double s = 0.;
  for(double r = dr/2; r < R; r += dr)
  {
    double beta_r = (beta_T*(n+2)/2) * TMath::Power((r/R),n);
    double rho = TMath::ATanH(beta_r);

    s += r * dr * TMath::BesselK1((mt*cosh(rho))/T) * TMath::BesselI0((pt*sinh(rho))/T);
  }

  return s;
}

//______________________________________________________________________________
void func(int &npar, double *gin, double &f, double *par, int flag)
{
  double aka    = par[0];
  double aka1   = par[4];
  double T      = par[2];
  double beta_T = par[1];
  double n = par[3];

  double chi2 = 0.;

  for(int i = 0; i < 14; i++)
  {
    double m = 0.497;

    double pt = x1ks[i];
    double mt = sqrt(m*m + pt*pt);

    double a = 0.;
    a = aka;
  
    double theo = a * mt * integral(beta_T, T, n, pt, mt);
    double q = (y1ks[i] - theo)/ey1ks[i];

    chi2 += q*q;
  }

  for(int i = 1; i < 13; i++){

    double m1 = 1.115;

    double pt1 = x1la[i];
    double mt1 = sqrt(m1*m1 + pt1*pt1);

    double a1 = 0.;
    a1 = aka1;

    double theo1 = a1 * mt1 * integral(beta_T,T,n,pt1,mt1);
    double q1 = (y1la[i]-theo1)/ey1la[i];

    chi2 += q1*q1;
  }

  f = chi2;

  getChi2 = f;
  getNdf  = 14 + 12 - 5 - 1;
}
//______________________________________________________________________________
void blastwave()
{

std::ifstream ifs1(f1, std::ifstream::in);           //K0s 130-inf pp Data File
    
    Double_t x1, y1, ey1;
    
    TGraphErrors *g1 = new TGraphErrors();
    g1->SetName("g1");
    
    
    while(ifs1 >> x1 >> y1 >> ey1)
    {
      //g1->SetPoint(g1->GetN(), x1, y1);
      //g1->SetPointError(linenumber, 0, ey1);
      //linenumber++ ;
        int i=0;
        x1ks[i]=x1;
        y1ks[i]=y1;
        ey1ks[i]=ey1;
        i++;
    }
    
std::ifstream ifs2(f2, std::ifstream::in);           //Lambda 130-inf pp Data File
    
    Double_t x2, y2, ey2;
    
    TGraphErrors *g2 = new TGraphErrors();
    g2->SetName("g2");
    
    while(ifs2 >> x2 >> y2 >> ey2)
    {
      //g2->SetPoint(g2->GetN(), x2, y2);
      //g2->SetPointError(linenumber, 0, ey2);
      //linenumber++ ;
        int i=0;
        x1la[i]=x2;
        y1la[i]=y2;
        ey1la[i]=ey2;
        i++;
    }
   
    double beta_T,Tkin,aka,aka1,n;
    double ebeta_T,eTkin,eaka,eaka1,en;
    
      TMinuit * gMinuit;
      gMinuit = new TMinuit();

      //set the function to minimize with minuit;
      gMinuit->SetFCN(func);

      double arglist[10];
      arglist[0] = 1;
      Int_t ierflg = 0;

      gMinuit->mnexcm("SET ERR", arglist, 1,ierflg);

      gMinuit->mnparm(0, "aka",    10,    0.1,   1,    10000, ierflg);
      gMinuit->mnparm(1, "beta_T", 0.6,   0.0001,  0.005, 1.0,   ierflg);
      gMinuit->mnparm(2, "Tkin",   0.16,  0.0001,  0.005, 1.0,   ierflg);
      gMinuit->mnparm(3, "n",      1.0,   0.01,  0.1,  10.0,  ierflg);
      gMinuit->mnparm(4, "aka1",   10,    0.1,   1,    10000, ierflg);

      gMinuit->mnexcm("MIGRAD",    arglist,  1,   ierflg);
      gMinuit->mnexcm("MINOS",     arglist,  1,   ierflg);
      gMinuit->mnexcm("CALL FCN",  arglist,  1,   ierflg);
      //gMinuit->mnexcm("HESSE",     arglist,  1,   ierflg);

      gMinuit->GetParameter(0, aka,    eaka);
      gMinuit->GetParameter(1, beta_T, ebeta_T);
      gMinuit->GetParameter(2, Tkin,   eTkin);
      gMinuit->GetParameter(3, n,      en);
      gMinuit->GetParameter(4, aka1,   eaka1);

      gMinuit->SetErrorDef(1);
      //test = (TGraph*)gMinuit->Contour(100,1,2);
      //test->SetFillColor(42);
      
      cout<<"beta_T = "<<beta_T<<"   "<<"Tkin = "<<Tkin<<endl;
      
}

Thanks in advance!
Best,
Mike Reynolds

K0s_mult130-inf_pp_7000GeV.txt (650 Bytes) Lambda_mult130-inf_pp_7000GeV.txt (399 Bytes)

2

Hi,
There is a problem with your FCN function (func).
At the first evaluation done by Minuit on the initial parameter values you have an infinite result, which makes the computed derivatives equal to a nan , since it compute them by f(x+h)-f(x). If both f(x+h) and f(x) are infinite then the result is a nan.

You need to debug and validate first your code implementing the function

Lorenzo