How to plot a function expressed in an integral format?

ROOT
1

Hello dear rooters,
I am trying to plot a function expressed in the form of an integral. Specifically, my integrand is a function of two parameters (x,p). My lower limit is a constant but my upper limit itself is a function of p. The main variable to be integrated is taken to be x from 0 up to some function of p and then I need to plot this integral in terms of p. I am having difficulty to come up with a nice sample in root documentation as it seems most of the cases we have floating point values for the limits of integration. However, I have written some code in scipy but it is not responding. It gives me some understandable error but I cannot fix it. If you know python, I can upload it so that you can see it. If not, it would be nice if you know some root documentation being able to do it.

Thanks,

2

Hi,

yes, some code would be nice, as I do not even remotely understand what you’re describing. However, highly likely this is not a python issue, and so you’ll have more chance of good feedback in the ROOT support sub-forum.

Cheers,
Wim

3

Thanks for letting me know, Here is the code in scipy:

[code]# import the needed modules
import numpy as np
import matplotlib.pyplot as plt
from pylab import *
from scipy.integrate import quad

define the integrand function

def F(p, x):
return 4445.935np.power(143.753 + 19 * np.log(p/x), -1, dtype=np.int16) * np.power(x, 0.2354424, dtype=np.int16) * np.power(1-np.power(x, 0.375, dtype=np.int16), 5.55373, dtype=np.int16)
def Zmomentlow1(p, x):
return 2np.power(x, 2)F(p, x)
def Zmomentlow2(p, x):
return 2np.power(x, 1.7)F(p, x)
def Zmomenthigh(p, x):
return 2np.power(x, 2)*F(p, x)

def Zmoment§:
for p, ep in enumerate§: #make Zmoment(enegrgy) values
if ep>0 and ep<=c:
Zmoment[p] = Zmomentlow1 + Zmomentlow2
elif ep>c:
Zmoment[p] = Zmomenthigh

set the parameters

c = 5.0E6 # the energy corresponding to the knee of the cosmic ray flux
a = 0 # lower limit of first integration of low-energy domain
b = 1 # upper limit of second integration of low-energy domain
N = 1000 # number of alpha values to use in integrations

def alpha§:
return p/c

p_ary = linspace(1.0E2,1.0E11,N) # create an array of values of energy
value = zeros(N) # create an array to store results of integrations

do the integrations low1, low2 and high

idx = 0
for p in p_ary:
value[idx],errorlow = quad(Zmomentlow1, a, alpha§, args = §) + quad(Zmomentlow2, alpha§, b, args = §)
value[idx],errorhigh = quad(Zmomenthigh, a, b, args = §)
idx = idx + 0.001

plot the result

plt.figure(1)
plt.clf()
plt.xscale(‘log’)
plt.plot(p_ary,value)
plt.ylim([0, 0.007])

note the new syntax for the xlabel() and title() commands

that allows the use of LaTeX commands

xlabel(‘energy [GeV]’)
ylabel(‘Zmoment’)
title(‘Plot of Zmoment vs. Energy’)
plt.show()
plt.savefig(‘Figure1’)[/code]

Here is the error message I receive when running the code:

[quote]Traceback (most recent call last):
File “main4.py”, line 46, in
value[idx],errorlow = quad(Zmomentlow1, a, alpha§, args = §) + quad(Zmomentlow2, alpha§, b, args = §)
ValueError: too many values to unpack[/quote]

4

Hi,

nothing to do with (Py)ROOT then.

The ValueError is simply telling you that the return tuple on the right side of the assignment has more entries than there are assignable variables on the left (the infodict in this case).

Aside, it may be that for args you did not intend to write §, but (p,) rather.

Cheers,
Wim

5

p is actually the variable against which the integral must be plotter. So, I tried (p,) too but the error message is still the same. So, my understanding is that quad function is defined to have constant floating point lower and upper limits but I am giving it non-constant function namely alpha§ which I don’t know how to avoid. So, now I am trying to do it in C++ so that I can use known numerical methods (turning an integral into a series). Hopefully, I can figure that out.

Thanks though,

6

Yes, as that has nothing to do with the error, that was just an aside …

Again: you are missing a variable on the left side to unpack the ones that are created from the quad() call. Left 2, right 3, i.e. one short. Per the docs of quad, the third is “infodict”. Not that I know what that is, just that you need to receive it on the left.

Cheers,
Wim

7

Sorry, I am not getting the point.
Why then the next line is not giving the same error as I am using it exactly the same?
Could you please rewrite those two lines (46 and 47) based on the missing “infodict”?

8

[quote=“ashur”]Why then the next line is not giving the same error as I am using it exactly the same?[/quote]A priori b/c it’s not executed until it is reached, and the current line stops execution. Beyond that, they’re not exactly the same: one sums the results from quad, the other doesn’t.

Actually, reading the quad docs more carefully, I see that the extra infodict return variable is only provided if the full_output variable is set, which you don’t. So it’s really much simpler: you’re summing two tuples (in line 46), which does not yield the sum of the elements, but a concatenation of the tuple.

[quote=“ashur”]Could you please rewrite those two lines (46 and 47) based on the missing “infodict”?[/quote]Something like this is most likely what you intend to do:

And I really think it should be (I haven’t tried any of this as I have none of the required packages installed):value[idx],errorlow = map(sum(zip(quad(Zmomentlow1, a, alpha(p), args = (p,)), quad(Zmomentlow2, alpha(p), b, args = (p,)))))
Cheers,
Wim

9

Thanks,
I appreciate your patience and helpfulness.

10

New Error Message:

[quote]Traceback (most recent call last):
File “main4.py”, line 46, in
value[idx],errorlow = map(sum(zip(quad(Zmomentlow1, a, alpha§, args = (p,)), quad(Zmomentlow2, alpha§, b, args = (p,)))))
TypeError: map() requires at least two args[/quote]

And removing map, leaving sum(…+…) is giving the following error message:

[quote]Traceback (most recent call last):
File “main4.py”, line 46, in
value[idx], errorlow = sum(zip(quad(Zmomentlow1, a, alpha§, args = (p,)),quad(Zmomentlow2, alpha§, b, args = (p,))))
TypeError: ‘numpy.float64’ object is not iterable[/quote]

11

Hi,

misspelled a comma with a parens. Like I said, I don’t have any of the packages installed, so I didn’t actually try it.

The only point is that you want to sum each of the elements, right? So you can also just code it out:q1 = quad(...) q2 = quad(...) value[idx] = q1[0]+q2[0] errorlow = q1[1]+q2[1]
What ‘zip’ does is reorder the elements to create pairs, i.e. zip(q1, q2) results in [(q1[0], q2[0]), (q1[1], q2[1])]. What sum does, is it sums all elements in the list given. Finally, map maps the function it is given on each of the elements. Hence map(sum, zip(q1, q2)) does the same as the 4 lines above.

Cheers,
Wim

12

I am learning great lessons from you. Thank You,
Interestingly now there is no error message and it seems it is running. But, the resultant plot is “empty”. Now, I need to try to see what is the reason for not having a plot while its window is run and opened!

13

Here is the plot I am expecting to get:

snapshot1.png789×508 6.38 KB