I was given a translated version of a FORTRAN/PAW macro (kumac) in ROOT/C++ (which all my code is in). The macro itself first reads in data, this is then passed to a function (TF1), which is then plotted onto a canvas.
The problem is that the resulting plots (see attachments) look very different. Using my root version, I seem to get a very sharp kink at x=10^-2 and can’t figure out why. Odd behaviour can also be seen right at the top, at y=0.8, where the central (bright red) curve seems to bend into the axis.
Is there a difference in the way that root and paw deal with plotting functions which could explain this? Do I have to change a “smoothness”/number of points setting (SetNumberFitPoints doesn’t seem to work)?
The translated ROOT/C++ function:
#ifndef _fcalc_h
#define _fcalc_h 1
double fcalc(double* x, double* pars)
{
double Q2 = pars[0];
double f2Type = pars[1];
double Q2G1, Q2G2, XG1, XG2;
int IQ2L, IQ2H, IXL, IXH;
double TQ, TX, F2H, F2L, F2;
if((Q2<0.39) || (x[0]<1E-6))
{
return 1.0;
}
// calculate grid points
for(int i=0; i<160; i++)
{
if(i<=71)
{
Q2G1 = pow(10,(5.0*i/120.0))-0.7;
Q2G2 = pow(10,(5.0*(i+1)/120.0))-0.7;
}
else if(i==72)
{
Q2G1 = pow(10,(5.0*i/120.0))-0.7;
Q2G2 = pow(10,(((2.0*(i-71.0)/88.0)+3.0)));
}
else
{
Q2G1 = pow(10,(((2.0*(i-72.0)/88.0)+3.0)));
Q2G2 = pow(10,(((2.0*(i-71.0)/88.0)+3.0)));
}
if((Q2>=Q2G1) && (Q2<Q2G2))
{
IQ2L = i+1;
IQ2H = i+2;
TQ = (Q2-Q2G1)/(Q2G2-Q2G1);
}
}
for(int i=0; i<160; i++)
{
if(i<=79)
{
XG1 = pow(10,((6.0*i/120.0)-6.0));
XG2 = pow(10,((6.0*(i+1)/120.0)-6.0));
}
else if(i == 80)
{
XG1 = pow(10,((6.0*i/120.0)-6.0));
XG2 = pow(10,((2.0/80.0*(i-79))-2.0));
}
else
{
XG1 = pow(10,((2.0/80.0*(i-80))-2.0));
XG2 = pow(10,((2.0/80.0*(i-79))-2.0));
}
if((x[0]>=XG1) && (x[0]<XG2))
{
IXL = i+1;
IXH = i+2;
TX = (x[0]-XG1)/(XG2-XG1);
}
}
double f2High, f2Low;
if(f2Type == 0)
{
F2H = (1-TX)*f2p[IXL][IQ2H]+TX*f2p[IXH][IQ2H];
F2L = (1-TX)*f2p[IXL][IQ2L]+TX*f2p[IXH][IQ2L];
}
if(f2Type == 1)
{
F2H = (1-TX)*f2up[IXL][IQ2H]+TX*f2up[IXH][IQ2H];
F2L = (1-TX)*f2up[IXL][IQ2L]+TX*f2up[IXH][IQ2L];
}
if(f2Type == 2)
{
F2H = (1-TX)*f2dn[IXL][IQ2H]+TX*f2dn[IXH][IQ2H];
F2L = (1-TX)*f2dn[IXL][IQ2L]+TX*f2dn[IXH][IQ2L];
}
F2 = (1-TQ)*F2L+TQ*F2H;
if(F2>0)
{
return F2;
}
else
{
return 1E-4;
}
};
#endifThe original FORTRAN function:
REAL FUNCTION FCALC(X)
VECTOR F2V(161,161)
VECTOR Q2PASS(1)
VECTOR RESULT(1)
Q2=Q2PASS(1)
IF (Q2.LT.0.39.OR.X.LT.1E-6) THEN
FCALC=1.
RESULT(1)=1.
RETURN
ENDIF
c calculate grid points
DO I=0,159
IF(I.LE.71)THEN
Q2G1=10**(5.*I/120.)-0.7
Q2G2=10**(5.*(I+1)/120.)-0.7
ELSEIF(I.EQ.72)THEN
Q2G1=10**(5*I/120)-0.7
Q2G2=10**(2*(I-71)/88 +3.)
ELSE
Q2G1=10**(2.*(I-72)/88. +3.)
Q2G2=10**(2.*(I-71)/88. +3.)
ENDIF
IF (Q2.GE.Q2G1.AND.Q2.LT.Q2G2) THEN
IQ2L=I+1
IQ2H=I+2
TQ=(Q2-Q2G1)/(Q2G2-Q2G1)
ENDIF
ENDDO
DO I=0,159
IF(I.LE.79)THEN
XG1=10**(6.*I/120.-6.)
XG2=10**(6.*(I+1)/120.-6.)
if(I.EQ.0)then
end if
ELSEIF(I.EQ.80)THEN
XG1=10**(6*I/120-6.)
XG2=10**(2./80.*(I-79)-2.)
ELSE
XG1=10**(2./80.*(I-80)-2.)
XG2=10**(2./80.*(I-79)-2.)
ENDIF
IF (X.GE.XG1.AND.X.LT.XG2) THEN
IXL=I+1
IXH=I+2
TX=(X-XG1)/(XG2-XG1)
ENDIF
ENDDO
F2H=(1-TX)*F2V(IXL,IQ2H)+TX*F2V(IXH,IQ2H)
F2L=(1-TX)*F2V(IXL,IQ2L)+TX*F2V(IXH,IQ2L)
FCALC=(1-TQ)*F2L+TQ*F2H
IF (FCALC.GT.0) THEN
RESULT(1)=FCALC
ELSE
RESULT(1)=1E-4
ENDIF
RETURN
ENDMy root macro (I’ve attached the relevant data file):
#include "fcalc.h"
double f2p[162][162];
double f2up[162][162];
double f2dn[162][162];
int fplot()
{
// Read in the data
TString filename = "ff_tot.dat";
ifstream datafile(filename);
if(!datafile.is_open())
{
cerr<<endl<<"Error: can not open file called "<<filename<<endl;
return -1;
}
int i=0, j=0;
while(datafile.good())
{
datafile>>f2p[i][j]>>f2up[i][j]>>f2dn[i][j];
i++;
if(i==161)
{
i=0;
j++;
}
}
datafile.close();
// Plot the functions using fcalc
TCanvas* c1 = new TCanvas("c1","c1",100,100,1000,1000);
TF1* funcP;
TF1* funcUp;
TF1* funcDn;
c1->cd();
double xmin = 0.0001;
double xmax = 0.99;
funP = new TF1("funP",fcalc,xmin,xmax,2);
funP->SetParameter(0,100.0);
funP->SetParameter(1,0);
funP->SetLineColor(50);
funP->SetLineWidth(1.0);
funP->SetMinimum(0.0);
funP->SetMaximum(0.8);
funP->SetNumberFitPoints(2000);
funP->Draw("C");
funUp = new TF1("funUp",fcalc,xmin,xmax,2);
funUp->SetParameter(0,100.0);
funUp->SetParameter(1,1);
funUp->SetLineColor(48);
funUp->SetLineWidth(1.0);
funUp->Draw("C same");
funDn = new TF1("funDn",fcalc,xmin,xmax,2);
funDn->SetParameter(0,100.0);
funDn->SetParameter(1,2);
funDn->SetLineColor(45);
funDn->SetLineWidth(1.0);
funDn->Draw("C same");
gPad->SetLogx();
c1->Print("c1.eps");
return 0;
}c1.eps.gz (3.31 KB)
paw1.eps.gz (3.66 KB)
ff_tot.txt (987 KB)