I’m trying to plot a Bethe-Bloch function with ROOT.
This is my code so far:
import numpy as np
import ROOT
c1 = ROOT.TCanvas()
pi = 3.1415
e0 = 8.854*10**(-12)
e = 1.602*10**(-19)
z = -e
me = 9.109*10**(-31)
c = 299792458
NA = 6.022*10**23
Z = 6
A = 12
p = 1032
I = 16*Z**0.9*1.602*10**(-19)
# set ROOT function
def betheFunct(x):
#BetheF = -(1/(4*pi*e0**2))*(e**4*z**2/(me*c**2))*(NA*Z/A)*p*(c**2/x**2)*(ROOT.TMath.ln(2*me*x**2/(I**2*(1-x**2/c**2)))-x**2/c**2)
BetheFTMath = -(1/(4*pi*ROOT.TMath.Power(e0,2)))*(ROOT.TMath.Power(e,4)*ROOT.TMath.Power(z,2)/(me*ROOT.TMath.Power(c,2)))*(NA*Z/A)*p*(ROOT.TMath.Power(c,2)/ROOT.TMath.Power(x,2))*(ROOT.TMath.ln(2*me*ROOT.TMath.Power(x,2)/(ROOT.TMath.Power(I,2)*(1-ROOT.TMath.Power(x,2)/ROOT.TMath.Power(c,2))))-ROOT.TMath.Power(x,2)/ROOT.TMath.Power(c,2))
return BetheFTMath
ROOTBethe = ROOT.TF1('Bethe', betheFunct, float(299792000), float(299792300))
ROOTBethe.Draw()
c1.SaveAs("Bethe.pdf","pdf")
At first I tried do define my power functions with the standard ** python command (The ‘BetheF’ that is commented out), but it gives the following error:
TypeError: unsupported operand type(s) for ** or pow(): ‘ROOT.PyDoubleBuffer’ and ‘int’
Then I tried to replace the ** with the ROOT Power Function but this also nets an error:
TypeError: none of the 2 overloaded methods succeeded. Full details:
Double_t TMath::Power(Double_t x, Int_t y) =>
could not convert argument 1 (a float is required)
Double_t TMath::Power(Double_t x, Double_t y) =>
could not convert argument 1 (a float is required)
Now I guess there should be an easy solution, but I just don’t see where the float is missing for the first argument? I think it is related to x not being set as a float, but I don’t know why.
Thanks a lot!
the python function receives its argument as an array, so you’ll want to use ‘x[0]’ rather than ‘x’. (The function pointer expected by the TF1 from the C++ side is “(*fcn)(Double_t *, Double_t *)”, the second array being an array of parameters.)
Unrelated to your original question, but I would recommend doing
pow = ROOT.TMath.Power
ln = ROOT.TMath.ln
at the beginning of your function. Then in your formula you could juse use “pow” and “ln” instead of spelling out the whole name, to make it more readable.
Also I’m not sure, but is defining your constants like “6.022*10**23” equivalent to “6.022e23”, or does it actually do the exponentiation with integers and then multiplication? In python it probably works out, but for other languages that could be a bug due to overflows and such.
actually, I’d recommend using module math in the first place, as it’d be much faster than going through the pyroot layer (I just timed it, and the difference is ~5x). The constant is fine, as ‘**’ takes precedence over ‘*’. It is slower, though, as the expression needs to be evaluated, so I wouldn’t use that construct.
Interesting, is there any difference in the actual computation between using ROOT.TMath and the math module? They probably both use cmath.h under the hoods…could you make ROOT.TMath just point at math to make it run faster?
I’m going to be search-and-replacing some ROOT.TMath in my code now that I know there is a difference in speed.
nothing to do with the actual Power() function, but with the layering in PyROOT. E.g. Power() is 5x overloaded and resolving overloads happens dynamically, otoh math.pow simply assumes you want to always work in double.
Aside from performance, there are other reasons why you may want to prefer math.pow, e.g. it will give a proper python exception for overflow, whereas Power() wiil give ‘inf’.
Btw., I just pulled and see that the recent changes from Axel (we’ve been going over performance) has dropped the difference down to 4x. I also still have to add the use of the lower level call interface (I’ve only tested that for longs, for CHEP), that should give another ~2x or so. In all, when everything is said and done, the final difference should be ~2x, not ~5x. Not that that helps you today, as this is ROOT6 master, but hey.