Merge pull request #13 from amanmdesai/test-complex-propagator

update cross section method

pymcabc/constants.py

@@ -2,14 +2,12 @@
 
 import math
 
-pi = math.pi
-delta = 2.0  # delta theta (radian)
 convert = 3.894e8  # GeV^-2 to pb
-g = 1.0
+g = 1
+#g = math.sqrt(math.sqrt(20/3))
 
 
 def outgoing_p(Ecm, m3, m4):
-    E = (Ecm**2 + m3**2 - m4**2) / (2 * Ecm)
+    E = (Ecm**2 - m3**2 - m4**2) / (2 * Ecm)
     E2 = E**2
-    p = E2 - m3**2
-    return math.sqrt(p)
+    return math.sqrt(E2)

pymcabc/cross_section.py

@@ -17,21 +17,19 @@ def __init__(self):
         self.mx = library["mx"][0]
         self.Ecm = library["Ecm"][0]
         self.g = pymcabc.constants.g
-        self.pi = pymcabc.constants.pi
-        self.delta = pymcabc.constants.delta
-        self.p_i = library["pi"][0]  # math.sqrt((self.Ecm / 2) ** 2 - (self.m1) ** 2)
+        self.p_f = library["outgoing_p"][0]
+        self.bw = library["bw"][0]
+        self.p_i = library["pi"][0]  
+        # math.sqrt((self.Ecm / 2) ** 2 - (self.m1) ** 2)        
+        # #self.E1 + self.E2
+        #self.E1 = library["E1"][0]
+        #self.E2 = library["E2"][0]
 
     def s_channel(self):
-        """definition for s channel"""
-        # deno = self.Ecm**2 - self.mx**2
-        deno = math.sqrt(self.p_i**2 + self.m1**2) + math.sqrt(self.p_i**2 + self.m2**2)
-        deno = deno**2 - self.mx**2
-        #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
-
+        #deno = math.sqrt(self.p_i**2 + self.m1**2) + math.sqrt(self.p_i**2 + self.m2**2)
+        deno = self.Ecm**2 - self.mx**2  + complex(0,1) * self.mx * self.bw
+        return self.g**2 / deno
+
     def t_channel(self, costh, pf):
         """definition for t channel"""
         deno = (
@@ -45,10 +43,8 @@ def t_channel(self, costh, pf):
             )
             + (2 * self.p_i * pf * costh)
         )
-        if abs(deno) <= 0.09:
-            return (self.g**2) / (deno + 100)
-        else:
-            return (self.g**2) / deno
+        deno = deno + complex(0,1) * self.mx * self.bw
+        return self.g**2 / deno
 
     def u_channel(self, costh, pf):
         """definition for u channel"""
@@ -63,10 +59,8 @@ def u_channel(self, costh, pf):
             )
             - (2 * self.p_i * pf * costh)
         )
-        if abs(deno) <= 0.09:
-            return (self.g**2) / (deno + 100)
-        else:
-            return (self.g**2) / deno
+        deno = deno + complex(0,1) * self.mx * self.bw
+        return self.g**2 / deno
 
 
 class CrossSection:
@@ -75,29 +69,35 @@ class for cross section calculation
     """
 
     def __init__(self):
-        self.pi = pymcabc.constants.pi
-        self.delta = pymcabc.constants.delta
         with open("library.json", "r") as f:
             library = json.load(f)
         self.Ecm = library["Ecm"][0]
         self.m1 = library["m1"][0]
+        self.m2 = library["m2"][0]
         self.m3 = library["m3"][0]
         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 = library["pi"][0]  # math.sqrt((self.Ecm / 2) ** 2 - (self.m1) ** 2)
+        self.p_f = library["outgoing_p"][0]
+        self.p_i = library["pi"][0]  
         self.channel = library["channel"][0]
+        # math.sqrt((self.Ecm / 2) ** 2 - (self.m1) ** 2)
+        #self.E1 = library["E1"][0]
+        #self.E2 = library["E2"][0]
+        #self.E1 + self.E2
+        #self.p1 = math.sqrt(self.E1**2 - self.m1**2) 
+        #self.p2 = math.sqrt(self.E2**2 - self.m2**2) 
+        #self.phase_factor = math.sqrt( (self.E1*self.E2 + self.p1*self.p2)**2 - (self.m1*self.m2)**2)
+        #self.p_f = pymcabc.constants.outgoing_p(self.Ecm, self.m3, self.m4)
 
     def dsigma_st(self, costh):
         if self.channel == "s":
             ME = MatrixElement().s_channel()
         elif self.channel == "t":
             ME = MatrixElement().t_channel(costh, self.p_f)
         else:
-            ME = MatrixElement().s_channel() + MatrixElement().t_channel(
-                costh, self.p_f
-            )
-        dsigma_st = 1 / ((8 * self.Ecm * self.pi) ** 2)
+            ME = MatrixElement().s_channel() + MatrixElement().t_channel(costh, self.p_f)
+        ME = abs(ME)
+        dsigma_st = 1 / ((self.Ecm*8  * math.pi) ** 2)
         dsigma_st = dsigma_st * abs(self.p_f / self.p_i) * ME**2
         return dsigma_st
 
@@ -107,24 +107,23 @@ def dsigma_tu(self, costh):
         elif self.channel == "u":
             ME = MatrixElement().u_channel(costh, self.p_f)
         else:
-            ME = MatrixElement().t_channel(costh, self.p_f) + MatrixElement().u_channel(
-                costh, self.p_f
-            )
-        dsigma_tu = 0.5 / ((self.Ecm * 8 * self.pi) ** 2)
+            ME = MatrixElement().t_channel(costh, self.p_f) + MatrixElement().u_channel(costh, self.p_f)
+        ME = abs(ME)
+        dsigma_tu = 0.5 / ((self.Ecm* 8 * math.pi) ** 2)
         dsigma_tu = dsigma_tu * abs(self.p_f / self.p_i) * ME**2
         return dsigma_tu
 
     def xsection(self, w_max):
-        costh = -1 + random.random() * self.delta
+        costh = -1 + random.random() * 2
         if self.process == "st":
-            w_i = CrossSection().dsigma_st(costh) * self.delta
+            w_i = CrossSection().dsigma_st(costh) * 2 * 2 * math.pi  
         elif self.process == "tu":
-            w_i = CrossSection().dsigma_tu(costh) * self.delta
+            w_i = CrossSection().dsigma_tu(costh) * 2 * 2 * math.pi 
         if w_max < w_i:
             w_max = w_i
         return w_i, w_max
 
-    def integrate_xsec(self, N=40000):
+    def integrate_xsec(self, N=20000):
         w_sum = 0
         w_max = 0
         w_square = 0
@@ -141,14 +140,14 @@ def integrate_xsec(self, N=40000):
             json.dump(library, f)
         return None
 
-    def calc_xsection(self, N: int = 40000):
+    def calc_xsection(self, N: int = 20000):
         self.integrate_xsec(N)
         with open("library.json", "r") as f:
             library = json.load(f)
         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)  # result in 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)

pymcabc/decay_particle.py

@@ -20,7 +20,7 @@ def __init__(self):
         self.decay1_mass = library["decay1_mass"][0]
         self.decay2_mass = library["decay2_mass"][0]
         self.massive = library["massive_mass"][0]
-        self.delta = pymcabc.constants.delta
+        self.delta = 2
 
     def rotate(self, pdecay: Particle):
         """rotate particle"""

pymcabc/generate_event.py

@@ -2,7 +2,6 @@
 import random
 import json
 import numpy as np
-import pymcabc.constants
 from pymcabc.cross_section import CrossSection
 
 
@@ -12,16 +11,16 @@ class GENEvents:
     """
 
     def __init__(self, Nevent: int):
-        self.delta = pymcabc.constants.delta
+        self.delta = 2
         self.Nevent = Nevent
         with open("library.json", "r") as f:
             library = json.load(f)
         self.w_max = library["w_max"][0]
-        self.Ecm = library["Ecm"][0]
+        self.Ecm = library["Ecm"][0] #+ library["E2"][0]
         self.m3 = library["m3"][0]
         self.m4 = library["m4"][0]
         self.process = library["process_type"][0]
-        self.out_p = pymcabc.constants.outgoing_p(self.Ecm, self.m3, self.m4)
+        self.out_p = library["outgoing_p"][0]
 
     def gen_events(self):
         i = 0
@@ -36,12 +35,12 @@ def gen_events(self):
         p2_pz = np.zeros(self.Nevent)
         p2_e = np.zeros(self.Nevent)
         while i < self.Nevent:
-            costh = -1 + (random.random() * self.delta)
+            costh = -1 + (random.random() * 2)
             phi = 2 * math.pi * random.random()  # * self.delta
             if self.process == "st":
-                w_ii = CrossSection().dsigma_st(costh) * self.delta
+                w_ii = CrossSection().dsigma_st(costh) *  2 * 2 * math.pi 
             elif self.process == "tu":
-                w_ii = CrossSection().dsigma_tu(costh) * self.delta
+                w_ii = CrossSection().dsigma_tu(costh) *  2 * 2 * math.pi 
             prob = w_ii / self.w_max
             random_point = random.random()
             if random_point < prob:

pymcabc/identify_process.py

@@ -1,5 +1,6 @@
 import json
 import math
+import pymcabc.constants
 
 
 def build_json():
@@ -12,6 +13,8 @@ def build_json():
         "m3": [],
         "m4": [],
         "mx": [],
+        "outgoing_p": [],
+        "bw": [],
         "massive": [],
         "massive_mass": [],
         "decay_process": [],
@@ -41,7 +44,8 @@ class DefineProcess:
         mA (float): mass of particle A
         mB (float): mass of particle B
         mC (float): mass of particle C
-        Ecm (float): center of mass energy
+        E1 (float): center of mass energy for beam 1
+        E2 (float): center of mass energy for beam 2
         channel (str): optional, use to study effect a particular channel
     """
 
@@ -52,18 +56,10 @@ def __init__(
         mB: float,
         mC: float,
         pi: float,
+        #E1: float,
+        #E2: float,
         channel: str = "none",
     ):
-        """
-        Defines the process, masses of particles and center of mass energy
-        Parameters:
-            input_string (str): Physics process. Example > 'A A > B B'
-            mA (float): mass of particle A
-            mB (float): mass of particle B
-            mC (float): mass of particle C
-            Ecm (float): center of mass energy
-            channel (str): optional, use to study effect a particular channel
-        """
 
         build_json()
         with open("library.json", "r") as f:
@@ -73,21 +69,28 @@ def __init__(
         self.mB = mB
         self.mC = mC
         self.p_i = pi
+        #self.E1 = E1
+        #self.E2 = E2
         if self.mA < 0 or self.mB < 0 or self.mC < 0:
             raise Exception("Negative masses not accepted")
+        #if self.E1 < 0:
+        #    raise Exception("Negative 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["pi"].append(pi)
+        #self.library["E2"].append(self.E2)
         self.library["channel"].append(channel)
         self.process()
         self.channel()
         self.masses()
-        self.ECM()
         self.identify_mediator()
         self.identify_decay()
+        self.ECM()
+        self.final_momenta()
+        self.bw()
 
     def process(self):
         """identify the physics process"""
@@ -143,20 +146,30 @@ def masses(self):
         with open("library.json", "w") as f:
             json.dump(self.library, f)
         return None
-
+    
     def ECM(self):
-        """center of mass energy"""
+        #center of mass energy
         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)
+        if len(self.library["Ecm"]) == 0:
+            self.library["Ecm"].append(Ecm)
         with open("library.json", "w") as f:
             json.dump(self.library, f)
-        return Ecm
+        print(
+              "\n",
+            "Energy Beam 1 : ", E1,
+              "\n",
+            "Energy Beam 2 : ", E2,
+              "\n",
+            "Energy CM : ", Ecm,
+              )
+        return E1, E2, Ecm
+
 
     def identify_mediator(self):
         """identify the mediator of the process"""
@@ -191,6 +204,18 @@ def identify_mediator(self):
         with open("library.json", "w") as f:
             json.dump(self.library, f)
         return None
+
+    def final_momenta(self):
+        p_f = pymcabc.constants.outgoing_p(self.library["Ecm"][0], self.library["m3"][0], self.library["m4"][0])
+        self.library["outgoing_p"].append(p_f)
+        with open("library.json", "w") as f:
+            json.dump(self.library, f)
+
+    def bw(self):
+        deno  = 8*math.pi*(self.library["mx"][0])**2
+        self.library["bw"].append((pymcabc.constants.g**2*self.library["outgoing_p"][0])/deno)
+        with open("library.json", "w") as f:
+            json.dump(self.library, f)
 
     def identify_decay(self):
         """identify the decay chain associated with the process"""

pymcabc/save_events.py

@@ -42,7 +42,7 @@ def __init__(
         with open("library.json", "r") as f:
             library = json.load(f)
         self.w_max = library["w_max"][0]
-        self.Ecm = library["Ecm"][0]
+        self.Ecm = library["Ecm"][0] 
         input_string = library["process"][0]
         input_string = input_string.replace(" > ", " ")
         input_string = input_string.split(" ")

setup.py

@@ -5,7 +5,7 @@
 
 setuptools.setup(
     name="pymcabc",
-    version="0.6.0",
+    version="0.7.0",
     description="Monte Carlo Event Generator for the ABC theory",
     author="Aman Desai",
     author_email="[email protected]",

tests/test_crosssection.py

@@ -5,12 +5,13 @@
 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 < 9e-11, "Sigma over estimated"
-    assert sigma > 9e-13, "Sigma over estimated"
+    print(sigma)
+    assert sigma > 7e-13, "Sigma under estimated"
+    assert sigma < 1e-11, "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"
+    assert sigma > 1e-13, "Sigma under estimated"
+    assert sigma < 7e-11, "Sigma over estimated"