# TFormula in Python

Hello,
i need to integrate a very complicated function. I tried with scipy but the numerical integration is incorrect. So, i decided to use ROOT for numerical integration.
The function is a two-dimensional function, very long function! I decided to use python because i need to operate with matrices and trigonometric functions, if you suppose that it is better to use C++ please tell me.

The main problem is that i can’t implement my function inside a TF2. A very simple example:

``````import ROOT
import numpy as np

def myFunc(x,y):
return y*x

f = ROOT.TF2("f","myFunc(x,y)",0,10**8,0,2*ROOT.TMath.Pi())

wf1 = ROOT.Math.WrappedMultiTF1(f)
ig.SetFunction(wf1)
ig.SetRelTolerance(0.001)

xmin = ROOT.std.vector('double')(2)
xmax = ROOT.std.vector('double')(2)

xmin[0] = 0
xmin[1] = 0
xmax[0] = 10**8
xmax[1] = ROOT.TMath.PiOver2()

print ig.Integral(np.array(xmin),np.array(xmax))
print f.Integral(0,10**8,0,ROOT.TMath.PiOver2())
``````

This gives the error:

``````Error in <TFormula::Eval>: Formula is invalid and not ready to execute. myFunc is unknown.
``````

The complete code is:

``````import ROOT
import numpy as np
from math import sin, cos
from sympy import simplify
from scipy import integrate

pi = ROOT.TMath.Pi()
e = ROOT.TMath.E()

mu0 = 4*pi*10**(-7)
e0 = 8.854*10**(-12)
c = 299792458
h = 1.05457*10**(-34)
wr = 1

wk, wm, wkr, wpl, gk, gm, gr, al, A = 0.7*wr, wr, wr, 0, 0.05*wr, 0.005*wr, 0.005*wr, 0.001, 0.001
gpl, we = 0.005*wpl, wkr

kk = lambda x,t: 1j*wk*(x*cos(t))/(-(x*cos(t))**2-wkr**2-gk*(x*cos(t)))
mu = lambda x,t: 1 + al + A*(x*cos(t))**2/(-(x*cos(t))**2-wm**2-gm*(x*cos(t)))
eps = lambda x,t: 1 - wpl**2/(-(x*cos(t))**2-gpl*(x*cos(t))) - we**2/(-(x*cos(t))**2-wr**2-gr*(x*cos(t)))

eta0 = np.sqrt(mu0/e0)
eta = lambda x,t: eta0*np.sqrt(mu(x,t)/eps(x,t))
n1 = lambda x,t: np.sqrt(eps(x,t)*mu(x,t)) + kk(x,t)
n2 = lambda x,t: np.sqrt(eps(x,t)*mu(x,t)) - kk(x,t)

gamm1 = lambda x,t: simplify((eta0**2 + eta(x,t)**2)/(2*eta(x,t)*eta0))
gamm2 = lambda x,t: simplify((eta0**2 - eta(x,t)**2)/(2*eta(x,t)*eta0))
chi1 = lambda x,t: simplify(np.sqrt((x*sin(t))**2+n1(x,t)*(x*cos(t))**2)/(n1(x,t)*x))
chi2 = lambda x,t: simplify(np.sqrt((x*sin(t))**2+n2(x,t)*(x*cos(t))**2)/(n2(x,t)*x))

delta = lambda x,t: simplify(gamm1(x,t)*(chi1(x,t)+chi2(x,t)) + chi1(x,t)*chi2(x,t) + 1)

rss = lambda x,t: simplify((-gamm2(x,t)*(chi1(x,t)+chi2(x,t)) - chi1(x,t)*chi2(x,t) + 1)/delta(x,t))
rpp = lambda x,t: simplify((gamm2(x,t)*(chi1(x,t)+chi2(x,t)) - chi1(x,t)*chi2(x,t) + 1)/delta(x,t))
rsp = lambda x,t: simplify(1j*(chi1(x,t) - chi2(x,t))/delta(x,t))
rps = lambda x,t: simplify(-1j*(chi1(x,t) - chi2(x,t))/delta(x,t))

r1 = lambda x,t: np.array([[rss(x,t),rsp(x,t)],[rps(x,t),rpp(x,t)]],dtype='float')
r2 = lambda x,t: r1(x,t)
identity = np.array([[1,0],[0,1]])

g = lambda d,x,t: np.subtract(identity,np.dot(np.dot(r1(x,t),r2(x,t)),identity*e**(-2*x*d)))

funz = lambda d,x,t: h/((2*pi)**2)*np.log(np.linalg.det(g(d,x,t)))*c*sin(t)*x**2

""" ---------------------------------------------- """

def myFunc(x,y):
return funz(10**(-8),x,y)

f = ROOT.TF2("f","myFunc(x,y)",0,10**8,0,2*ROOT.TMath.Pi())

wf1 = ROOT.Math.WrappedMultiTF1(f)
ig.SetFunction(wf1)
ig.SetRelTolerance(0.001)

xmin = ROOT.std.vector('double')(2)
xmax = ROOT.std.vector('double')(2)

xmin[0] = 0
xmin[1] = 0
xmax[0] = 10**8
xmax[1] = ROOT.TMath.PiOver2()

print ig.Integral(np.array(xmin),np.array(xmax))
print f.Integral(0,10**8,0,ROOT.TMath.PiOver2())

funz2 = lambda x,y: funz(10**(-8),x,y)

print punti
``````

Thanks so much,
Alberto.

Hi Alberto,

You can pass a Python callable to the constructor of `TF2` from Python:

``````f = ROOT.TF2("f",myFunc,0,10**8,0,2*ROOT.TMath.Pi())
``````

You can see examples of its use here for `TF1`:

Please note the callable should receive two arrays, one with the variables and one with the parameters of the function.

Enric

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.