Merge pull request #11 from amanmdesai/update-saveevent

update save event and detector

pymcabc/cross_section.py

@@ -24,9 +24,9 @@ def __init__(self):
     def s_channel(self):
         """definition for s channel"""
         # 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
+        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:

pymcabc/decay_particle.py

@@ -45,7 +45,7 @@ def decay(self, top: Particle):
         )
 
         E2 = top.mass() - E1
-
+        """
         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)
@@ -65,7 +65,8 @@ def decay(self, top: Particle):
                 )
             )
         )
-
+
+        """
         # decay2.mass() = self.decay2_mass
 
         E1 = math.sqrt(

pymcabc/detector.py

@@ -26,6 +26,7 @@ def gauss_smear(self, particle: Particle):
             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
                 )
@@ -39,18 +40,16 @@ def gauss_smear(self, particle: Particle):
                 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_px[i] = random.gauss(particle.px[i], self.factor*self.sigma)
+                output_py[i] = random.gauss(particle.py[i], self.factor*self.sigma)
+                output_pz[i] = random.gauss(particle.pz[i], self.factor*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
-                )
-                """
 
             particle_output = Particle(output_E, output_px, output_py, output_pz)
         return particle_output
+
+#mass = output_E[i] ** 2 - (output_px[i] ** 2 + output_py[i] ** 2 + output_pz[i] ** 2)
+

pymcabc/feynman_diagram.py

@@ -35,18 +35,16 @@ def __init__(self):
                 elif p == "u":
                     self.u_chan()
                 else:
-                    print("Possible channels: s, t, and u")
                     return 0
+        elif channel =="s":
+            self.s_chan()
+        elif channel == "t":
+            self.t_chan()
+        elif channel == "u":
+            self.u_chan()
         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
+            print("Possible channels: s, t, and u")
+            return 0
 
     def s_chan(self):
         fig = plt.figure(figsize=(5.0, 5.0))

pymcabc/save_events.py

@@ -48,6 +48,10 @@ def __init__(
         input_string = input_string.split(" ")
         self.output_1 = input_string[2]
         self.output_2 = input_string[3]
+        if self.output_1 == self.output_2:
+            self.output_1 = self.output_1 + "_1"
+            self.output_2 = self.output_2 + "_2"
+
         self.decay_process = library["decay_process"]
         if self.decay_process[0] != "NaN":
             decay_split = self.decay_process[0].replace(" > ", " ")
@@ -117,14 +121,14 @@ def to_root(self, name: str = "ABC_events.root"):
             decay1, decay2 = DecayParticle().prepare_decay(top1)
             decay3, decay4 = DecayParticle().prepare_decay(top2)
             if self.boolDetector == True:
-                if decay1.px[0] == -9 and decay1.E[0] == -9:
-                    self.top1 = Detector(
-                        self.detector_sigma, self.detector_factor
-                    ).gauss_smear(self.top1)
-                if decay2.px[0] == -9 and decay2.E[0] == -9:
-                    self.top2 = Detector(
-                        self.detector_sigma, self.detector_factor
-                    ).gauss_smear(self.top2)
+                #if decay1.px[0] == -9 and decay1.E[0] == -9:
+                self.top1 = Detector(
+                    self.detector_sigma, self.detector_factor
+                ).gauss_smear(self.top1)
+                #if decay2.px[0] == -9 and decay2.E[0] == -9:
+                self.top2 = Detector(
+                    self.detector_sigma, self.detector_factor
+                ).gauss_smear(self.top2)
                 # self.top1 = Detector(self.detector_sigma,self.detector_factor).gauss_smear(self.top1)
                 # self.top2 = Detector(self.detector_sigma,self.detector_factor).gauss_smear(self.top2)
 
@@ -141,6 +145,7 @@ def to_root(self, name: str = "ABC_events.root"):
                     self.detector_sigma, self.detector_factor
                 ).gauss_smear(decay4)
 
+
             file["events"] = {
                 self.output_1 + "_E": self.top1.E,
                 self.output_1 + "_Px": self.top1.px,

setup.py

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

tests/test_plot.py

@@ -3,11 +3,21 @@
 
 
 # test if the file is prepared
-def test_plot():
+def test_plot_st():
+    pymcabc.DefineProcess("A B > A 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_st.root"
+    )
+    pymcabc.PlotData.file("test_eventGen_detector_decay_st.root")
+    assert "B_E.png" in os.listdir(), "file not created"
+
+def test_plot_tu():
     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"
+        "test_eventGen_detector_decay_tu.root"
     )
-    pymcabc.PlotData.file("test_eventGen_detector_decay.root")
-    assert "B_E_decay_A.png" in os.listdir(), "file not created"
+    pymcabc.PlotData.file("test_eventGen_detector_decay_tu.root")
+    assert "B_1_E_decay_A.png" in os.listdir(), "file not created"
+