Add feynman diagram (#5)

* add feynman diagram

* add example notebook

* update readme and new version

* add formatting

* add dependencies

* add actions

.github/workflows/publish.yml

@@ -36,7 +36,6 @@ jobs:
       with:
         name: artifact
         path: dist
-
-    - uses: pypa/[email protected]
+    - uses: pypa/[email protected]
       with:
         password: ${{ secrets.pypi_password }}

.gitignore

@@ -127,3 +127,8 @@ dmypy.json
 
 # Pyre type checker
 .pyre/
+
+
+notebooks/*png
+*root
+*.json

README.md

@@ -26,13 +26,17 @@ pip install pymcabc
 The physics of ABC model (theory) is described in Grifiths.
 
 ## Simple script to start using the package:
+1. Import pymcabc:
+```python
+import pymcabc
+```
 
 1. Define the process, for example:
 ```python
 pymcabc.DefineProcess('A A > B B',mA=4,mB=10,mC=1,Ecm=30)
 ```
 
-2. Calculate the total cross section of the process:
+2. Calculate the total cross section of the process (in barn):
 ```python
 pymcabc.CrossSection().calc_xsection()
 ```

notebooks/README.md

@@ -1 +1 @@
-## How to use notebooks?
+Example Notebook

notebooks/example.ipynb

@@ -0,0 +1,85 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "85a8f377",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "import pymcabc"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0a1f4ab0",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pymcabc.DefineProcess('A A > B B',mA=4,mB=10,mC=1,Ecm=30)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "56f7abc3",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pymcabc.FeynmanDiagram() ## draw feynman diagrams for the process"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "2d556baf",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pymcabc.CrossSection().calc_xsection()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "f1bf311b",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pymcabc.SaveEvent(10000,boolDecay=True,boolDetector=True).to_root('name.root')"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "0d6d763d",
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "pymcabc.PlotData.file('name.root')"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.10.6"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

pymcabc/__init__.py

@@ -1,5 +1,5 @@
 """
-Copyright (c) 2022 Aman Desai. All rights reserved.
+Copyright (c) 2023 Aman Desai. All rights reserved.
 """
 
 
@@ -8,3 +8,4 @@
 from pymcabc.save_events import SaveEvent
 from pymcabc.generate_event import GENEvents
 from pymcabc.plotting import PlotData
+from pymcabc.feynman_diagram import FeynmanDiagram

pymcabc/cross_section.py

@@ -5,8 +5,8 @@
 
 
 class MatrixElement:
-    """ internal class for matrix element calculation
-    """
+    """internal class for matrix element calculation"""
+
     def __init__(self):
         with open("library.json", "r") as f:
             library = json.load(f)
@@ -58,6 +58,7 @@ class CrossSection:
     """
     class for cross section calculation
     """
+
     def __init__(self):
         self.pi = pymcabc.constants.pi
         self.delta = pymcabc.constants.delta
@@ -70,26 +71,30 @@ def __init__(self):
         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.channel =  library["channel"][0]
+        self.channel = library["channel"][0]
 
     def dsigma_st(self, costh):
-        if self.channel == 's':
+        if self.channel == "s":
             ME = MatrixElement().s_channel()
-        elif self.channel == 't':
+        elif self.channel == "t":
             ME = MatrixElement().t_channel(costh, self.p_f)
         else:
-            ME = MatrixElement().s_channel() + MatrixElement().t_channel(costh, self.p_f)
+            ME = MatrixElement().s_channel() + MatrixElement().t_channel(
+                costh, self.p_f
+            )
         dsigma_st = 1 / ((8 * self.Ecm * self.pi) ** 2)
         dsigma_st = dsigma_st * abs(self.p_f / self.p_i) * ME**2
         return dsigma_st
 
     def dsigma_tu(self, costh):
-        if self.channel == 't':
+        if self.channel == "t":
             ME = MatrixElement().t_channel(costh, self.p_f)
-        elif self.channel == 'u':
+        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)
+            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)
         dsigma_tu = dsigma_tu * abs(self.p_f / self.p_i) * ME**2
         return dsigma_tu
@@ -121,7 +126,7 @@ 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 = 40000):
         self.integrate_xsec(N)
         with open("library.json", "r") as f:
             library = json.load(f)

pymcabc/decay_particle.py

@@ -10,6 +10,7 @@ class DecayParticle:
     """
     decays particle
     """
+
     def __init__(self):
         # self.Ecm = library["Ecm"][0]
         with open("library.json", "r") as f:
@@ -23,10 +24,9 @@ def __init__(self):
         self.massive = library["massive_mass"][0]
         self.delta = pymcabc.constants.delta
 
-    def rotate(self,pdecay: Particle,size: int):
-        """rotate particle
-        """
-        costh =  (np.random.rand(size) * 2) -1
+    def rotate(self, pdecay: Particle, size: int):
+        """rotate particle"""
+        costh = (np.random.rand(size) * 2) - 1
         sinth = np.sqrt(1 - costh**2)
         phi = 2 * math.pi * np.random.rand(size)
         sinPhi = np.sin(phi)
@@ -40,34 +40,42 @@ def rotate(self,pdecay: Particle,size: int):
 
         return pdecay
 
-
     def decay(self, top: Particle):
         """decay particle"""
-        self.decay_p = 1 / (2 * top.mass()) * np.sqrt(
-            (self.mA**4+ self.mB**4+ self.mC**4)
-            - 2 * (self.mA**2 * self.mB**2 +  self.mA**2 * self.mC**2 + self.mB**2 * self.mC**2)
+        self.decay_p = (
+            1
+            / (2 * top.mass())
+            * np.sqrt(
+                (self.mA**4 + self.mB**4 + self.mC**4)
+                - 2
+                * (
+                    self.mA**2 * self.mB**2
+                    + self.mA**2 * self.mC**2
+                    + self.mB**2 * self.mC**2
+                )
+            )
         )
 
         decay1 = Particle(-9, -9, -9, -9)
         decay2 = Particle(-9, -9, -9, -9)
         # decay2.mass() = self.decay2_mass
 
-        E1 =  np.sqrt(
+        E1 = np.sqrt(
             self.decay1_mass * self.decay1_mass + self.decay_p * self.decay_p
-        )*np.ones(top.E.shape[0])
+        ) * np.ones(top.E.shape[0])
         E2 = np.sqrt(
-            self.decay2_mass * self.decay2_mass  + self.decay_p * self.decay_p
-        )*np.ones(top.E.shape[0])
+            self.decay2_mass * self.decay2_mass + self.decay_p * self.decay_p
+        ) * np.ones(top.E.shape[0])
 
         decay1.set4momenta(E1, 0, self.decay_p, 0)
         decay2.set4momenta(E2, 0, -self.decay_p, 0)
 
-        decay1 = self.rotate(decay1, size = top.E.shape[0])
-        decay2 = self.rotate(decay2, size = top.E.shape[0])
+        decay1 = self.rotate(decay1, size=top.E.shape[0])
+        decay2 = self.rotate(decay2, size=top.E.shape[0])
 
         return decay1, decay2
 
-    def nearlyequal(self,a, b):
+    def nearlyequal(self, a, b):
         if abs(a - b) < 0.001:
             return True
         else:
@@ -77,14 +85,16 @@ def prepare_decay(self, top: Particle):
         if self.decay_process != "NaN":
             # decay_process = decay_process.replace(" < "," ")
             # decay_process = decay_process.split(" ")
-            #print(top.mass()[0], self.massive)
-            if self.nearlyequal(top.mass()[0], self.massive) and top.mass()[0]>(self.decay1_mass + self.decay2_mass):
+            # print(top.mass()[0], self.massive)
+            if self.nearlyequal(top.mass()[0], self.massive) and top.mass()[0] > (
+                self.decay1_mass + self.decay2_mass
+            ):
                 decay1, decay2 = self.decay(top)
                 decay1 = decay1.boost(top)
                 decay2 = decay2.boost(top)
                 return decay1, decay2
             else:
-                output = np.ones(top.E.size)*(-9)
-                decay1 = Particle(output,output, output, output)
-                decay2 = Particle(output,output, output, output)
+                output = np.ones(top.E.size) * (-9)
+                decay1 = Particle(output, output, output, output)
+                decay2 = Particle(output, output, output, output)
                 return decay1, decay2

pymcabc/detector.py

@@ -4,19 +4,21 @@
 
 class Detector:
     """Applies gaussian smearing on E and momenta"""
+
     def identify_smear(particle: Particle, type: str = "gauss"):
         if type == "gauss":
             particle = self.gauss_smear(particle)
         else:
             print("Type Not found")
         return particle
+
     def gauss_smear(self, particle: Particle, sigma: float = 0.5):
-        size=particle.E.shape[0]
-        if particle.px[0] ==-9 and particle.py[0] == -9:
+        size = particle.E.shape[0]
+        if particle.px[0] == -9 and particle.py[0] == -9:
             return particle
         else:
-            particle.px = np.random.normal(particle.px,sigma,size)
-            particle.py = np.random.normal(particle.py,sigma,size)
-            particle.pz = np.random.normal(particle.pz,sigma,size)
-            particle.E = np.random.normal(particle.E,sigma,size)
+            particle.px = np.random.normal(particle.px, sigma, size)
+            particle.py = np.random.normal(particle.py, sigma, size)
+            particle.pz = np.random.normal(particle.pz, sigma, size)
+            particle.E = np.random.normal(particle.E, sigma, size)
         return particle