Merge pull request #10 from amanmdesai/beam-energy

use initial momentum as parameter

README.md

@@ -33,7 +33,7 @@ import pymcabc
 
 1. Define the process, for example:
 ```python
-pymcabc.DefineProcess('A A > B B',mA=4,mB=10,mC=1,Ecm=30)
+pymcabc.DefineProcess('A A > B B',mA=4,mB=10,mC=1,pi=30)
 ```
 
 2. Calculate the total cross section of the process (in barn):

notebooks/example.ipynb

@@ -17,7 +17,7 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "pymcabc.DefineProcess('A A > B B',mA=4,mB=10,mC=1,Ecm=30)"
+    "pymcabc.DefineProcess('A A > B B',mA=4,mB=10,mC=1,pi=30)"
    ]
   },
   {

pymcabc/constants.py

@@ -9,7 +9,7 @@
 
 
 def outgoing_p(Ecm, m3, m4):
-    p = (Ecm**2 + m3**2 - m4**2) / (2 * Ecm)
-    p = p**2
-    p = p - m3**2
+    E = (Ecm**2 + m3**2 - m4**2) / (2 * Ecm)
+    E2 = E**2
+    p = E2 - m3**2
     return math.sqrt(p)

pymcabc/cross_section.py

@@ -19,12 +19,15 @@ def __init__(self):
         self.g = pymcabc.constants.g
         self.pi = pymcabc.constants.pi
         self.delta = pymcabc.constants.delta
-        self.p_i = math.sqrt((self.Ecm / 2) ** 2 - (self.m1) ** 2)
+        self.p_i = library["pi"][0]  # math.sqrt((self.Ecm / 2) ** 2 - (self.m1) ** 2)
 
     def s_channel(self):
         """definition for s channel"""
-        deno = self.Ecm**2 - self.mx**2
-        if abs(deno) <= 0.01:
+        # deno = self.Ecm**2 - self.mx**2
+        deno = (self.p_i**2 + self.m1**2) * (self.p_i**2 + self.m2**2)
+        deno = math.sqrt(deno)
+        deno = deno + self.m1**2 + self.m2**2
+        if abs(deno) <= 0.09:
             return (self.g**2) / (deno + 100)
         else:
             return (self.g**2) / deno
@@ -35,7 +38,11 @@ def t_channel(self, costh, pf):
             self.m1**2
             + self.m3**2
             - self.mx**2
-            - (self.Ecm * math.sqrt(pf**2 + self.m3**2))
+            - (
+                2
+                * math.sqrt(self.p_i**2 + self.m1**2)
+                * math.sqrt(pf**2 + self.m3**2)
+            )
             + (2 * self.p_i * pf * costh)
         )
         if abs(deno) <= 0.09:
@@ -49,7 +56,11 @@ def u_channel(self, costh, pf):
             self.m1**2
             + self.m4**2
             - self.mx**2
-            - (self.Ecm * math.sqrt(pf**2 + self.m4**2))
+            - (
+                2
+                * math.sqrt(self.p_i**2 + self.m1**2)
+                * math.sqrt(pf**2 + self.m4**2)
+            )
             - (2 * self.p_i * pf * costh)
         )
         if abs(deno) <= 0.09:
@@ -74,7 +85,7 @@ def __init__(self):
         self.m4 = library["m4"][0]
         self.process = library["process_type"][0]
         self.p_f = pymcabc.constants.outgoing_p(self.Ecm, self.m3, self.m4)
-        self.p_i = math.sqrt((self.Ecm / 2) ** 2 - (self.m1) ** 2)
+        self.p_i = library["pi"][0]  # math.sqrt((self.Ecm / 2) ** 2 - (self.m1) ** 2)
         self.channel = library["channel"][0]
 
     def dsigma_st(self, costh):
@@ -137,8 +148,8 @@ def calc_xsection(self, N: int = 40000):
         w_sum = library["w_sum"][0]
         w_square = library["w_square"][0]
         w_max = library["w_max"][0]
-        sigma_x = w_sum * pymcabc.constants.convert / (N * 1e12)
-        variance = math.sqrt(abs(w_square / N - (w_sum / N) ** 2))  # barn unit
+        sigma_x = w_sum * pymcabc.constants.convert / (N * 1e12)  # result in barn unit
+        variance = math.sqrt(abs((w_square / N) - (w_sum / N) ** 2))  # barn unit
         error = (
             variance * pymcabc.constants.convert / (math.sqrt(N) * 1e12)
         )  # barn unit

pymcabc/detector.py

@@ -46,10 +46,11 @@ def gauss_smear(self, particle: Particle):
                 output_py[i] = random.gauss(particle.py[i], self.sigma)
                 output_pz[i] = random.gauss(particle.pz[i], self.sigma)
                 output_E[i] = random.gauss(particle.E[i], self.sigma)
-                """
+
                 mass = output_E[i] ** 2 - (
                     output_px[i] ** 2 + output_py[i] ** 2 + output_pz[i] ** 2
                 )
-                print(mass)
+                """
+
             particle_output = Particle(output_E, output_px, output_py, output_pz)
         return particle_output

pymcabc/feynman_diagram.py

@@ -16,8 +16,16 @@ def __init__(self):
         process = library["process"]
         mediator = library["mediator"][0]
         channel = library["channel"][0]
+
+        decay_process = library["decay_process"]
+        if decay_process !="NaN":
+            decay_process = decay_process[0].replace(">","")
+            self.decay_process = decay_process.split()
+            self.decay()
+
         process = process[0].replace(">", mediator)
         self.process = process.split()
+
         if channel == "none":
             for p in library["process_type"][0]:
                 if p == "s":
@@ -120,3 +128,23 @@ def u_chan(self):
 
         diagram.plot()
         plt.savefig("uchan.pdf")
+
+    def decay(self):
+        fig = plt.figure(figsize=(10.,10.))
+        ax = fig.add_axes([0,0,1,1], frameon=False)
+
+        diagram = Diagram(ax)
+        in1 = diagram.vertex(xy=(.1,.5), marker='')
+        v1  = diagram.vertex(xy=(.5,.5))
+        out1 = diagram.vertex(xy=(.9,.75),marker='')
+        out2 = diagram.vertex(xy=(.9,.25),marker='')
+
+        a = diagram.line(in1,v1,arrow=False)
+        b = diagram.line(v1,out1,arrow=False)
+        c = diagram.line(v1,out2,arrow=False)
+        a.text(self.decay_process[0],fontsize=30,t=.1,y=.1)
+        b.text(self.decay_process[1],fontsize=30,t=-.01,y=.1)
+        c.text(self.decay_process[2],fontsize=30,t=-.1, y=-.1)
+
+        diagram.plot()
+        plt.savefig('decay.pdf')

pymcabc/identify_process.py

@@ -1,4 +1,5 @@
 import json
+import math
 
 
 def build_json():
@@ -17,6 +18,7 @@ def build_json():
         "decay1_mass": [],
         "decay2_mass": [],
         "mediator": [],
+        "pi": [],
         "Ecm": [],
         "process": [],
         "process_type": [],
@@ -49,7 +51,7 @@ def __init__(
         mA: float,
         mB: float,
         mC: float,
-        Ecm: float,
+        pi: float,
         channel: str = "none",
     ):
         """
@@ -70,14 +72,15 @@ def __init__(
         self.mA = mA
         self.mB = mB
         self.mC = mC
-        self.Ecm = Ecm
-        if  self.mA<0 or self.mB<0 or self.mC < 0:
+        self.p_i = pi
+        if self.mA < 0 or self.mB < 0 or self.mC < 0:
             raise Exception("Negative masses not accepted")
-        if self.Ecm < 0:
-            raise Exception("Negative center of mass energy not accepted")
+        if self.p_i <= 0:
+            raise Exception("Negative or Zero absolute momentum not accepted")
         self.library["mA"].append(mA)
         self.library["mB"].append(mB)
         self.library["mC"].append(mC)
+        self.library["pi"].append(self.p_i)
         self.library["channel"].append(channel)
         self.process()
         self.channel()
@@ -143,10 +146,17 @@ def masses(self):
 
     def ECM(self):
         """center of mass energy"""
-        self.library["Ecm"].append(self.Ecm)
+        with open("library.json", "r") as f:
+            library = json.load(f)
+        m1 = library["m1"][0]
+        m2 = library["m2"][0]
+        E1 = math.sqrt(m1**2 + self.p_i**2)
+        E2 = math.sqrt(m2**2 + self.p_i**2)
+        Ecm = E1 + E2
+        self.library["Ecm"].append(Ecm)
         with open("library.json", "w") as f:
             json.dump(self.library, f)
-        return None
+        return Ecm
 
     def identify_mediator(self):
         """identify the mediator of the process"""
@@ -216,4 +226,4 @@ def identify_decay(self):
             self.library["decay_process"].append("NaN")
         with open("library.json", "w") as f:
             json.dump(self.library, f)
-        return None
+        return None

pymcabc/plotting.py

@@ -16,17 +16,12 @@ def plot(data, key):
         label = key.replace("_", " ")
         plt.xlabel(label)
         plt.ylabel("Events")
-        # plt.show()
         plt.savefig(key + ".png")
 
     def file(filename="ABC_events.root"):
-        # if ".root" in filename:
         file = uproot.open(filename)
         tree = file["events"]
         branches = tree.arrays()
         for key in tree.keys():
             PlotData.plot(branches[key], key)
-        # if ".csv" in filename:
-        # df = pd.read_csv(filename)
-        # for col in df.columns:
-        #    PlotData.plot(df[col], col)
+

setup.py

@@ -1,28 +1,30 @@
 import setuptools
 
-with open('README.md','r') as file:
+with open("README.md", "r") as file:
     long_description = file.read()
 
 setuptools.setup(
-name="pymcabc",
-version="0.3.0",
-description= "Monte Carlo Event Generator for the ABC theory",
-author="Aman Desai",
-author_email="[email protected]",
-maintainer="Aman Desai",
-maintainer_email="[email protected]",
-url = "https://github.com/amanmdesai/pymcabc",
-long_description=long_description,
-packages=["pymcabc"],
-install_requires=["numpy","uproot","matplotlib","feynman"],
-long_description_content_type="text/markdown",
-classifiers=["License :: OSI Approved :: MIT License",
-             "Programming Language :: Python :: 3",
-             "Programming Language :: Python :: 3.8",
-             "Programming Language :: Python :: 3.9",
-             "Programming Language :: Python :: 3.10",
-             "Topic :: Scientific/Engineering :: Physics",
-             "Operating System :: Unix",
-             "Operating System :: Microsoft :: Windows",
-             "Operating System :: MacOS"]
+    name="pymcabc",
+    version="0.5.0",
+    description="Monte Carlo Event Generator for the ABC theory",
+    author="Aman Desai",
+    author_email="[email protected]",
+    maintainer="Aman Desai",
+    maintainer_email="[email protected]",
+    url="https://github.com/amanmdesai/pymcabc",
+    long_description=long_description,
+    packages=["pymcabc"],
+    install_requires=["numpy", "uproot", "matplotlib", "feynman"],
+    long_description_content_type="text/markdown",
+    classifiers=[
+        "License :: OSI Approved :: MIT License",
+        "Programming Language :: Python :: 3",
+        "Programming Language :: Python :: 3.8",
+        "Programming Language :: Python :: 3.9",
+        "Programming Language :: Python :: 3.10",
+        "Topic :: Scientific/Engineering :: Physics",
+        "Operating System :: Unix",
+        "Operating System :: Microsoft :: Windows",
+        "Operating System :: MacOS",
+    ],
 )

tests/test_crosssection.py

@@ -2,8 +2,15 @@
 import os
 
 
-def test_xsec():
-    pymcabc.DefineProcess("A A > B B", mA=4, mB=10, mC=1, Ecm=30)
+def test_xsec_st():
+    pymcabc.DefineProcess("A B > A B", mA=4, mB=10, mC=1, pi=15)
     sigma, error = pymcabc.CrossSection().calc_xsection()
-    assert sigma >= 10e-14, "Sigma under estimated"
-    assert sigma <= 10e-12, "Sigma over estimated"
+    assert sigma < 9e-11, "Sigma over estimated"
+    assert sigma > 9e-13, "Sigma over estimated"
+
+
+def test_xsec_tu():
+    pymcabc.DefineProcess("A A > B B", mA=4, mB=10, mC=1, pi=15)
+    sigma, error = pymcabc.CrossSection().calc_xsection()
+    assert sigma < 4e-13, "Sigma over estimated"
+    assert sigma > 1e-14, "Sigma over estimated"

tests/test_evengen.py

@@ -4,7 +4,7 @@
 
 # test if the file is prepared
 def test_eventGen_detector_decay():
-    pymcabc.DefineProcess("A A > B B", mA=4, mB=10, mC=1, Ecm=30)
+    pymcabc.DefineProcess("A A > B B", mA=4, mB=10, mC=1, pi=30)
     pymcabc.CrossSection().calc_xsection()
     pymcabc.SaveEvent(100, boolDecay=True, boolDetector=True).to_root(
         "test_eventGen_detector_decay.root"
@@ -14,7 +14,7 @@ def test_eventGen_detector_decay():
 
 
 def test_eventGen_decay():
-    pymcabc.DefineProcess("A A > B B", mA=4, mB=10, mC=1, Ecm=30)
+    pymcabc.DefineProcess("A A > B B", mA=4, mB=10, mC=1, pi=30)
     pymcabc.CrossSection().calc_xsection()
     pymcabc.SaveEvent(100, boolDecay=True, boolDetector=False).to_root(
         "test_eventGen_decay.root"
@@ -24,7 +24,7 @@ def test_eventGen_decay():
 
 
 def test_eventGen_detector():
-    pymcabc.DefineProcess("A A > B B", mA=4, mB=10, mC=1, Ecm=30)
+    pymcabc.DefineProcess("A A > B B", mA=4, mB=10, mC=1, pi=30)
     pymcabc.CrossSection().calc_xsection()
     pymcabc.SaveEvent(100, boolDecay=False, boolDetector=True).to_root(
         "test_eventGen_detector.root"
@@ -34,7 +34,7 @@ def test_eventGen_detector():
 
 
 def test_eventGen():
-    pymcabc.DefineProcess("A A > B B", mA=4, mB=10, mC=1, Ecm=30)
+    pymcabc.DefineProcess("A A > B B", mA=4, mB=10, mC=1, pi=30)
     pymcabc.CrossSection().calc_xsection()
     pymcabc.SaveEvent(100, boolDecay=False, boolDetector=False).to_root(
         "test_eventGen.root"