Skip to content

Commit 91638a3

Browse files
committed
Added PhaseSpace examples
Signed-off-by: Deepanshu <deepanshu2017@gmail.com>
1 parent e7a9fe9 commit 91638a3

6 files changed

Lines changed: 286 additions & 1 deletion

File tree

examples/PhaseSpace/constants.py

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,7 @@
1+
RED = "\033[1;31m"
2+
BLUE = "\033[1;34m"
3+
CYAN = "\033[1;36m"
4+
GREEN = "\033[0;32m"
5+
RESET = "\033[0;0m"
6+
BOLD = "\033[;1m"
7+
REVERSE = "\033[;7m"
Lines changed: 65 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,65 @@
1+
import HydraPython as hypy
2+
import sys
3+
import time
4+
import math
5+
import constants
6+
7+
8+
def functor(*four_vector):
9+
p1 = four_vector[0]
10+
p2 = four_vector[1]
11+
p3 = four_vector[2]
12+
13+
p = p1 + p2 + p3
14+
q = p2 + p3
15+
16+
pd = p * p2
17+
pq = p * q
18+
qd = q * p2
19+
mp2 = p.mass2()
20+
mq2 = q.mass2()
21+
md2 = p2.mass2()
22+
23+
return (pd * mq2 - pq * qd) / math.sqrt((pq * pq - mq2 * mp2) * (qd * qd - mq2 * md2))
24+
25+
26+
def main():
27+
"""
28+
This example shows how to use the Hydra's
29+
phase space Monte Carlo algorithms to calculate the
30+
average value and corresponding variance of a functor
31+
over the phase space of the decay B0 -> J/psi K pi.
32+
"""
33+
nentries = 1000000
34+
B0_mass = 5.27955
35+
Jpsi_mass = 3.0969
36+
K_mass = 0.493677
37+
pi_mass = 0.13957061
38+
if len(sys.argv) > 1:
39+
nentries = int(sys.argv[1])
40+
41+
B0 = hypy.Vector4R(B0_mass, 0.0, 0.0, 0.0)
42+
masses = [Jpsi_mass, K_mass, pi_mass]
43+
phsp = hypy.PhaseSpace3(B0_mass, masses)
44+
45+
# Device #
46+
47+
print()
48+
print(constants.RED, "ERROR: device is not supported with custom functors in current HydraPython version.", sep='')
49+
print(constants.RESET)
50+
51+
# Host #
52+
start = time.time()
53+
pair = phsp.AverageOnhost(B0, functor, nentries)
54+
end = time.time()
55+
56+
print('\n' * 2)
57+
print("----------------- Host ----------------")
58+
print("|< cos(theta_K) >(B0 -> J/psi K pi): ", pair[0], " +- ", pair[1])
59+
print("| Number of events :", nentries)
60+
print("| Time (ms) :", end - start)
61+
print("-----------------------------------------")
62+
63+
64+
if __name__ == '__main__':
65+
main()

examples/PhaseSpace/phsp_basic.py

Lines changed: 58 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,58 @@
1+
import HydraPython as hypy
2+
import sys
3+
import time
4+
5+
6+
def main():
7+
"""
8+
This example shows how to use the Hydra's
9+
phase space Monte Carlo algorithms to
10+
generate a sample of B0 -> J/psi K pi.
11+
"""
12+
nentries = 1000000
13+
B0_mass = 5.27955
14+
Jpsi_mass = 3.0969
15+
K_mass = 0.493677
16+
pi_mass = 0.13957061
17+
if len(sys.argv) > 1:
18+
nentries = int(sys.argv[1])
19+
20+
B0 = hypy.Vector4R(B0_mass, 0.0, 0.0, 0.0)
21+
masses = [Jpsi_mass, K_mass, pi_mass]
22+
phsp = hypy.PhaseSpace3(B0_mass, masses)
23+
24+
# Device #
25+
events_d = hypy.device_events_3(nentries)
26+
start = time.time()
27+
phsp.GenerateOndevice(B0, events_d)
28+
end = time.time()
29+
30+
print('\n' * 2)
31+
print("----------------- Device ----------------")
32+
print("| B0 -> J/psi K pi")
33+
print("| Number of events :", nentries)
34+
print("| Time (ms) :", end - start)
35+
print("-----------------------------------------")
36+
37+
for i in range(10):
38+
print(events_d[i])
39+
40+
# Host #
41+
events_h = hypy.host_events_3(nentries)
42+
start = time.time()
43+
phsp.GenerateOnhost(B0, events_h)
44+
end = time.time()
45+
46+
print('\n' * 2)
47+
print("------------------ Host -----------------")
48+
print("| B0 -> J/psi K pi")
49+
print("| Number of events :", nentries)
50+
print("| Time (ms) :", end - start)
51+
print("-----------------------------------------")
52+
53+
for i in range(10):
54+
print(events_h[i])
55+
56+
57+
if __name__ == '__main__':
58+
main()

examples/PhaseSpace/phsp_chain.py

Lines changed: 74 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,74 @@
1+
import HydraPython as hypy
2+
import sys
3+
import time
4+
5+
6+
def main():
7+
"""
8+
This example shows how to use the Hydra's
9+
phase space Monte Carlo algorithms to
10+
generate a sample of B0 -> J/psi K pi and
11+
use the multiple daughters to generate the
12+
grand-daughters.
13+
"""
14+
nentries = 1000000
15+
B0_mass = 5.27955
16+
Jpsi_mass = 3.0969
17+
K_mass = 0.493677
18+
pi_mass = 0.13957061
19+
mu_mass = 0.1056583745
20+
if len(sys.argv) > 1:
21+
nentries = int(sys.argv[1])
22+
23+
B0 = hypy.Vector4R(B0_mass, 0.0, 0.0, 0.0)
24+
masses1 = [Jpsi_mass, K_mass, pi_mass]
25+
masses2 = [mu_mass, mu_mass]
26+
27+
phsp1 = hypy.PhaseSpace3(B0_mass, masses1)
28+
phsp2 = hypy.PhaseSpace2(Jpsi_mass, masses2)
29+
30+
# Device #
31+
32+
daughters_d = hypy.host_events_3(nentries)
33+
grand_daughters_d = hypy.host_events_2(nentries)
34+
start = time.time()
35+
36+
phsp1.GenerateOnhost(B0, daughters_d)
37+
phsp2.GenerateOnhost(daughters_d.getDaughters(0), grand_daughters_d)
38+
39+
end = time.time()
40+
41+
print('\n' * 2)
42+
print("------------------ Host -----------------")
43+
print("| B0 -> J/psi K pi | J/psi -> mu+ mu-")
44+
print("| Number of events :", nentries)
45+
print("| Time (ms) :", end - start)
46+
print("-----------------------------------------")
47+
48+
for i in range(10):
49+
print(daughters_d[i], grand_daughters_d[i])
50+
51+
# Host #
52+
53+
daughters_h = hypy.host_events_3(nentries)
54+
grand_daughters_h = hypy.host_events_2(nentries)
55+
start = time.time()
56+
57+
phsp1.GenerateOnhost(B0, daughters_h)
58+
phsp2.GenerateOnhost(daughters_h.getDaughters(0), grand_daughters_h)
59+
60+
end = time.time()
61+
62+
print('\n' * 2)
63+
print("----------------- Device ----------------")
64+
print("| B0 -> J/psi K pi | J/psi -> mu+ mu-")
65+
print("| Number of events :", nentries)
66+
print("| Time (ms) :", end - start)
67+
print("-----------------------------------------")
68+
69+
for i in range(10):
70+
print(daughters_h[i], grand_daughters_h[i])
71+
72+
73+
if __name__ == '__main__':
74+
main()
Lines changed: 71 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,71 @@
1+
import HydraPython as hypy
2+
import sys
3+
import time
4+
import math
5+
import constants
6+
7+
8+
def functor(*four_vector):
9+
p1 = four_vector[0]
10+
p2 = four_vector[1]
11+
p3 = four_vector[2]
12+
13+
p = p1 + p2 + p3
14+
q = p2 + p3
15+
16+
pd = p * p2
17+
pq = p * q
18+
qd = q * p2
19+
mp2 = p.mass2()
20+
mq2 = q.mass2()
21+
md2 = p2.mass2()
22+
23+
return (pd * mq2 - pq * qd) / math.sqrt((pq * pq - mq2 * mp2) * (qd * qd - mq2 * md2))
24+
25+
26+
def main():
27+
"""
28+
This example shows how to use the Hydra's
29+
phase space Monte Carlo algorithms to
30+
generate events of B0 -> J/psi K pi on fly and
31+
evaluate a set of functors.
32+
"""
33+
nentries = 1000000
34+
B0_mass = 5.27955
35+
Jpsi_mass = 3.0969
36+
K_mass = 0.493677
37+
pi_mass = 0.13957061
38+
if len(sys.argv) > 1:
39+
nentries = int(sys.argv[1])
40+
41+
B0 = hypy.Vector4R(B0_mass, 0.0, 0.0, 0.0)
42+
masses = [Jpsi_mass, K_mass, pi_mass]
43+
phsp = hypy.PhaseSpace3(B0_mass, masses)
44+
45+
# Device #
46+
47+
print()
48+
print(constants.RED, "ERROR: device is not supported with custom functors in current HydraPython version.", sep='')
49+
print(constants.RESET)
50+
51+
# Host #
52+
result = hypy.host_vector_float2(nentries)
53+
start = time.time()
54+
phsp.EvaluateOnhost(B0, result, functor)
55+
end = time.time()
56+
57+
print('\n' * 2)
58+
print("----------------- Host ----------------")
59+
print("| B0 -> J/psi K pi")
60+
print("| Number of events :", nentries)
61+
print("| Time (ms) :", end - start)
62+
print("-----------------------------------------")
63+
64+
65+
for i in range(10):
66+
print("<", i, ">: ", result[i])
67+
68+
69+
70+
if __name__ == '__main__':
71+
main()

include/PyEvents.h

Lines changed: 11 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -117,6 +117,16 @@ namespace hydra_python {
117117
}, "idx"_a, \
118118
"Daughter iterator. Iterate over the all N events of given " \
119119
"particle idx.") \
120+
.def("getDaughters", \
121+
[](const hydra::Events<N, hydra::BACKEND::sys_t>& e, \
122+
hydra::GInt_t idx) { \
123+
hypy::BACKEND##_vector_float4 daughters(e.size()); \
124+
for (auto i = e.DaughtersBegin(idx); i != e.DaughtersEnd(idx); ++i) \
125+
daughters[idx] = *i; \
126+
return daughters;\
127+
}, "idx"_a, \
128+
"Daughters. Get all N daughters of given " \
129+
"particle idx.") \
120130
.def("Events", \
121131
[](const hydra::Events<N, hydra::BACKEND::sys_t>& e) { \
122132
typedef decltype(e.begin()) iter_t; \
@@ -144,7 +154,7 @@ namespace hydra_python {
144154
&hydra::Events<N, hydra::BACKEND::sys_t>::capacity) \
145155
.def("resize", &hydra::Events<N, hydra::BACKEND::sys_t>::resize, \
146156
"Resize the number of Events.") \
147-
.def("unweight", \
157+
.def("Unweight", \
148158
&hydra::Events<N, hydra::BACKEND::sys_t>::Unweight, \
149159
"Unweight all Events with seed.") \
150160
.def("size", &hydra::Events<N, hydra::BACKEND::sys_t>::size) \

0 commit comments

Comments
 (0)