From c7bb76198a3cb42352102cf8ee208b434bdb06cf Mon Sep 17 00:00:00 2001
From: Johannes Schumann <johannes.schumann@fau.de>
Date: Thu, 14 Oct 2021 11:35:25 +0200
Subject: [PATCH] =?UTF-8?q?Add=20Q=C2=B2=20to=20the=20output=20dataset?=
MIME-Version: 1.0
Content-Type: text/plain; charset=UTF-8
Content-Transfer-Encoding: 8bit

---
 km3buu/output.py | 41 +++++++++++++++++++++++++++++++++++------
 1 file changed, 35 insertions(+), 6 deletions(-)

diff --git a/km3buu/output.py b/km3buu/output.py
index cc975cf..d1119c0 100644
--- a/km3buu/output.py
+++ b/km3buu/output.py
@@ -346,13 +346,10 @@ class GiBUUOutput:
         return retval
 
     @staticmethod
-    def bjorken_x(roottuple_data):
+    def _q(roottuple_data):
         """
-        Calculate Bjorken x variable for the GiBUU events
-
-        Parameters
-        ----------
-            roottuple_data: awkward.highlevel.Array                
+        Calculate the difference of the four vectors from 
+        in and out going lepton (k_in - k_out)
         """
         d = roottuple_data
         k_in = np.vstack([
@@ -368,8 +365,39 @@ class GiBUUOutput:
             np.array(d.lepOut_Pz)
         ])
         q = k_in - k_out
+        return q
+
+    @staticmethod
+    def Qsquared(roottuple_data):
+        """
+        Calculate the passed energy from the neutrino interaction to the 
+        nucleus denoted by the variable Q²
+
+        Parameters
+        ----------
+            roottuple_data: awkward.highlevel.Array                
+
+        Returns
+        -------
+            
+        """
+        q = GiBUUOutput._q(roottuple_data)
         qtmp = np.power(q, 2)
         q2 = qtmp[0, :] - np.sum(qtmp[1:4, :], axis=0)
+        return q2
+
+    @staticmethod
+    def bjorken_x(roottuple_data):
+        """
+        Calculate Bjorken x variable for the GiBUU events
+
+        Parameters
+        ----------
+            roottuple_data: awkward.highlevel.Array                
+        """
+        d = roottuple_data
+        q = GiBUUOutput._q(d)
+        q2 = GiBUUOutput.Qsquared(d)
         pq = q[0, :] * d.nuc_E - q[1, :] * d.nuc_Px - q[2, :] * d.nuc_Py - q[
             3, :] * d.nuc_Pz
         x = np.divide(-q2, 2 * pq)
@@ -439,6 +467,7 @@ class GiBUUOutput:
         retval["xsec"] = self._event_xsec(retval)
         retval["Bx"] = GiBUUOutput.bjorken_x(retval)
         retval["By"] = GiBUUOutput.bjorken_y(retval)
+        retval["Q2"] = GiBUUOutput.Qsquared(retval)
         return retval
 
     @property
-- 
GitLab