Merge pull request #8 from amanmdesai/modify-numpy-dependencies

Modify numpy dependencies

.gitignore

@@ -131,4 +131,6 @@ dmypy.json
 
 notebooks/*png
 *root
-*.json
+*.json
+*.png
+*.pdf

pymcabc/__init__.py

@@ -9,3 +9,4 @@
 from pymcabc.generate_event import GENEvents
 from pymcabc.plotting import PlotData
 from pymcabc.feynman_diagram import FeynmanDiagram
+from pymcabc.particle import Particle

pymcabc/decay_particle.py

@@ -1,7 +1,6 @@
 import math
 import random
 import json
-import numpy as np
 import pymcabc.constants
 from pymcabc.particle import Particle
 
@@ -12,7 +11,6 @@ class DecayParticle:
     """
 
     def __init__(self):
-        # self.Ecm = library["Ecm"][0]
         with open("library.json", "r") as f:
             library = json.load(f)
         self.mA = library["mA"][0]
@@ -24,28 +22,35 @@ def __init__(self):
         self.massive = library["massive_mass"][0]
         self.delta = pymcabc.constants.delta
 
-    def rotate(self, pdecay: Particle, size: int):
+    def rotate(self, pdecay: Particle):
         """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)
-        cosPhi = np.cos(phi)
+        costh = (2*random.random()) - 1
+        sinth = math.sqrt(1 - costh**2)
+        phi = 2 * math.pi * random.random()
+        sinPhi = math.sin(phi)
+        cosPhi = math.cos(phi)
+        pdecay_out = Particle(0,0,0,0)
+        pdecay_out.px = pdecay.pz * sinth *cosPhi
+        pdecay_out.py = pdecay.pz * sinth *sinPhi
+        pdecay_out.pz = pdecay.pz * costh
+        pdecay_out.E = pdecay.E 
+
+        return pdecay_out
 
-        pdecay.px = pdecay.px * cosPhi - pdecay.py * sinPhi
-        pdecay.py = pdecay.px * sinPhi + pdecay.py * cosPhi
+    def decay(self, top: Particle):
+        """decay particle"""
 
-        pdecay.px = pdecay.px * costh - pdecay.pz * sinth
-        pdecay.pz = pdecay.px * sinth + pdecay.pz * costh
+        E1 = (top.mass()**2 - self.decay2_mass**2 + self.decay1_mass**2)/(2*top.mass())
+
+        E2 =  top.mass() - E1
 
-        return pdecay
+        self.decay_p = math.sqrt((top.mass()**2  - (self.decay1_mass + self.decay2_mass)**2) * (top.mass()**2  - (self.decay1_mass - self.decay2_mass)**2)) / (2*top.mass())
 
-    def decay(self, top: Particle):
-        """decay particle"""
+        """ 
         self.decay_p = (
             1
             / (2 * top.mass())
-            * np.sqrt(
+            * math.sqrt(
                 (self.mA**4 + self.mB**4 + self.mC**4)
                 - 2
                 * (
@@ -56,27 +61,30 @@ def decay(self, top: Particle):
             )
         )
 
-        decay1 = Particle(-9, -9, -9, -9)
-        decay2 = Particle(-9, -9, -9, -9)
         # decay2.mass() = self.decay2_mass
 
-        E1 = np.sqrt(
+        E1 = math.sqrt(
             self.decay1_mass * self.decay1_mass + self.decay_p * self.decay_p
-        ) * np.ones(top.E.shape[0])
-        E2 = np.sqrt(
+        ) 
+        E2 = math.sqrt(
             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 = Particle(0,0,0,0)
+        decay2 = Particle(0,0,0,0)
 
-        decay1 = self.rotate(decay1, size=top.E.shape[0])
-        decay2 = self.rotate(decay2, size=top.E.shape[0])
+        decay1.set4momenta(E1, 0, 0, self.decay_p)
+        decay2.set4momenta(E2, 0, 0, -self.decay_p)
+
+        decay1 = self.rotate(decay1)
+        decay2.set4momenta(E2, -decay1.px, -decay1.py, -decay1.pz)
+        #decay2 = self.rotate(decay2)
 
         return decay1, decay2
 
     def nearlyequal(self, a, b):
-        if abs(a - b) < 0.001:
+        if abs(a - b) < 1e-3:
             return True
         else:
             return False
@@ -86,15 +94,50 @@ def prepare_decay(self, top: Particle):
             # 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
-            ):
-                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)
-                return decay1, decay2
+            decay_array1_px = [0]*len(top.px)
+            decay_array1_py = [0]*len(top.px)
+            decay_array1_pz = [0]*len(top.px)
+            decay_array1_E = [0]*len(top.px)
+
+            decay_array2_px = [0]*len(top.px)
+            decay_array2_py = [0]*len(top.px)
+            decay_array2_pz = [0]*len(top.px)
+            decay_array2_E = [0]*len(top.px)
+
+            for i in range(len(top.px)):
+                heavy_state = Particle(top.E[i],top.px[i], top.py[i],top.pz[i])
+                if self.nearlyequal(heavy_state.mass(), self.massive) and heavy_state.mass() > (
+                    self.decay1_mass + self.decay2_mass
+                ):
+
+                    decay1, decay2 = self.decay(heavy_state)
+                    decay1 = decay1.boost(heavy_state)
+                    decay2 = decay2.boost(heavy_state)
+
+                    decay_array1_px[i] = decay1.px
+                    decay_array1_py[i] = decay1.py
+                    decay_array1_pz[i] = decay1.pz
+                    decay_array1_E[i] = decay1.E
+
+                    decay_array2_px[i] = decay2.px
+                    decay_array2_py[i] = decay2.py
+                    decay_array2_pz[i] = decay2.pz
+                    decay_array2_E[i] = decay2.E
+                else:
+                    output = -9
+                    decay_array1_px[i] = output
+                    decay_array1_py[i] = output
+                    decay_array1_pz[i] = output
+                    decay_array1_E[i] = output
+
+                    decay_array2_px[i] = output
+                    decay_array2_py[i] = output
+                    decay_array2_pz[i] = output
+                    decay_array2_E[i] = output
+
+            decay1 = Particle(decay_array1_E,decay_array1_px, decay_array1_py, decay_array1_pz)
+            decay2 = Particle(decay_array2_E,decay_array2_px, decay_array2_py, decay_array2_pz)
+            return decay1, decay2
+
+
+

pymcabc/detector.py

@@ -1,24 +1,47 @@
-import numpy as np
+import random
+import math
 from pymcabc.particle import Particle
 
 
 class Detector:
     """Applies gaussian smearing on E and momenta"""
+    def __init__(self, sigma: float =1., factor: float =1.):
+        self.sigma = sigma
+        self.factor = factor
 
-    def identify_smear(particle: Particle, type: str = "gauss"):
+    def identify_smear(self, 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]
+    def gauss_smear(self, particle: Particle):
         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)
-        return particle
+            output_px = [0]*len(particle.px)
+            output_py = [0]*len(particle.px)
+            output_pz = [0]*len(particle.px)
+            output_E = [0]*len(particle.px)
+            for i in range(len(particle.px)):
+
+                momentum =  math.sqrt(particle.px[i]**2  + particle.py[i]**2  + particle.pz[i]**2 )
+                random_measure_momentum = random.gauss(momentum, self.factor*self.sigma) / momentum
+                random_measure_energy = random.gauss(particle.E[i], self.sigma) / particle.E[i]
+
+                output_px[i] = random_measure_momentum*particle.px[i]
+                output_py[i] = random_measure_momentum*particle.py[i]
+                output_pz[i] = random_measure_momentum*particle.pz[i]
+
+                output_E[i] = random_measure_energy*particle.E[i]
+
+                """output_px[i] = random.gauss(particle.px[i], self.sigma)
+                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

@@ -15,16 +15,27 @@ def __init__(self):
 
         process = library["process"]
         mediator = library["mediator"][0]
+        channel = library["channel"][0]
         process = process[0].replace(">", mediator)
         self.process = process.split()
-
-        for p in library["process_type"][0]:
-            if p == "s":
-                self.s_chan()
-            elif p == "t":
-                self.t_chan()
-            elif p == "u":
-                self.u_chan()
+        if channel == "none":
+            for p in library["process_type"][0]:
+                if p == "s":
+                    self.s_chan()
+                elif p == "t":
+                    self.t_chan()
+                elif p == "u":
+                    self.u_chan()
+                else:
+                    print("Possible channels: s, t, and u")
+                    return 0
+        else:
+            if channel == "s":
+                    self.s_chan()
+            if channel == "t":
+                    self.t_chan()
+            if channel == "u":
+                    self.u_chan()
             else:
                 print("Possible channels: s, t, and u")
                 return

pymcabc/particle.py

@@ -60,16 +60,21 @@ def pT(self):
 
     def mass(self):
         """returns particle's mass"""
-        # try:
-        #    x = math.sqrt(self.E**2 - sum([self.px**2, self.py**2, self.pz**2]))
-        #    return x
-        # except:
-        #    return 0
-        x = self.E**2 - self.px**2 - self.py**2 - self.pz**2
-        if x[0] < 0:
-            x = np.zeros(self.E.size())
-        else:
-            x = np.sqrt(x)
+        try:
+            x = self.E**2 - self.px**2 - self.py**2 - self.pz**2
+            if x < 0:
+                x = 0
+            else:
+                x = math.sqrt(x)
+        except:
+            x = [0]*len(self.px)
+            for i in range(len(x)):
+                x[i] = self.E[i]**2 - self.px[i]**2 - self.py[i]**2 - self.pz[i]**2
+                print(x[i])
+                if x[i] < 0:
+                    x[i] = 0
+                else:
+                    x[i] = math.sqrt(x)
         return x
 
     def set4momenta(self, new_E, new_px, new_py, new_pz):
@@ -83,15 +88,23 @@ def boost(self, other):
         """
         boosts a particle four momentum.
          # boost motivated from ROOT TLorentzVector class
+
+         other is used to boost
         """
-        new = Particle(-9, -9, -9, -9)
-        new_other = Particle(-9, -9, -9, -9)
+        new = Particle(0., 0., 0., 0.)
+        new_other = Particle(0., 0., 0., 0.)
+
         new_other.set4momenta(
             other.E, other.px / other.E, other.py / other.E, other.pz / other.E
         )
-        beta = new_other.p()
-        gamma = 1.0 / np.sqrt(1 - beta**2)
-        gamma_2 = (gamma - 1.0) / beta
+        beta = new_other.p2()
+        gamma = 1.0 / math.sqrt(1.0 - beta)
+
+        if beta>0:
+            gamma_2 = (gamma - 1.0) / beta
+        else:
+            gamma_2 = 0.0
+
 
         dotproduct = (
             self.px * new_other.px + self.py * new_other.py + self.pz * new_other.pz
@@ -100,4 +113,4 @@ def boost(self, other):
         new.py = self.py + (gamma_2 * dotproduct + gamma * self.E) * new_other.py
         new.pz = self.pz + (gamma_2 * dotproduct + gamma * self.E) * new_other.pz
         new.E = gamma * (self.E + dotproduct)
-        return new
+        return new

pymcabc/plotting.py

@@ -16,7 +16,7 @@ def plot(data, key):
         label = key.replace("_", " ")
         plt.xlabel(label)
         plt.ylabel("Events")
-        plt.show()
+        #plt.show()
         plt.savefig(key + ".png")
 
     def file(filename="ABC_events.root"):