From 69d6d5ae436075bdaa21e524819be90d24f0dccd Mon Sep 17 00:00:00 2001
From: Fabian Joswig <fabian.joswig@ed.ac.uk>
Date: Wed, 20 Oct 2021 11:25:48 +0100
Subject: [PATCH] Npr_matrix class and read_ExternalLeg improved

---
 pyerrors/input/hadrons.py |  5 +++--
 pyerrors/npr.py           | 36 ++++++++++++++++++++----------------
 2 files changed, 23 insertions(+), 18 deletions(-)

diff --git a/pyerrors/input/hadrons.py b/pyerrors/input/hadrons.py
index b903b3bb..561887cb 100644
--- a/pyerrors/input/hadrons.py
+++ b/pyerrors/input/hadrons.py
@@ -57,10 +57,10 @@ def read_meson_hd5(path, filestem, ens_id, meson='meson_0', tree='meson'):
     corr_data = []
     for hd5_file in files:
         file = h5py.File(path + '/' + hd5_file, "r")
-
         raw_data = list(file[tree + '/' + meson + '/corr'])
         real_data = [o[0] for o in raw_data]
         corr_data.append(real_data)
+        file.close()
     corr_data = np.array(corr_data)
 
     l_obs = []
@@ -90,6 +90,7 @@ def read_ExternalLeg_hd5(path, filestem, ens_id, order='C'):
         file = h5py.File(path + '/' + hd5_file, "r")
         raw_data = file['ExternalLeg/corr'][0][0].view('complex')
         corr_data.append(raw_data)
+        file.close()
     corr_data = np.array(corr_data)
 
     mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(' '), dtype=int)
@@ -103,4 +104,4 @@ def read_ExternalLeg_hd5(path, filestem, ens_id, order='C'):
         matrix[si, sj, ci, cj] = CObs(real, imag)
         matrix[si, sj, ci, cj].gamma_method()
 
-    return Npr_matrix(matrix.reshape((12,12), order=order), mom_in=mom)
+    return Npr_matrix(matrix.swapaxes(1, 2).reshape((12,12), order=order), mom_in=mom)
diff --git a/pyerrors/npr.py b/pyerrors/npr.py
index 71b222ad..889d4f8b 100644
--- a/pyerrors/npr.py
+++ b/pyerrors/npr.py
@@ -73,27 +73,31 @@ class Npr_matrix(np.ndarray):
 
     @property
     def g5H(self):
-        new_matrix = Npr_matrix.g5 @ self.conj().T @ Npr_matrix.g5
+        """Gamma_5 hermitean conjugate
+
+        Returns gamma_5 @ M.T.conj() @ gamma_5 and exchanges in and out going
+        momenta. Works only for 12x12 matrices.
+        """
+        if self.shape != (12, 12):
+            raise Exception('g5H only works for 12x12 matrices.')
+        extended_g5 = np.kron(np.eye(3, dtype=int), Npr_matrix.g5)
+        new_matrix = extended_g5 @ self.conj().T @ extended_g5
         new_matrix.mom_in = self.mom_out
         new_matrix.mom_out = self.mom_in
         return new_matrix
 
+    def _propagate_mom(self, other, name):
+        s_mom = getattr(self, name, None)
+        o_mom = getattr(other, name, None)
+        if s_mom != o_mom and s_mom and o_mom:
+                raise Exception(name + ' does not match.')
+        return o_mom if o_mom else s_mom
+
     def __matmul__(self, other):
-        if hasattr(other, 'mom_in'):
-            if self.mom_in != other.mom_in and self.mom_in and other.mom_in:
-                    raise Exception('mom_in does not match.')
-            mom_in = self.mom_in if self.mom_in else other.mom_in
-        else:
-            mom_in = self.mom_in
-
-        if hasattr(other, 'mom_out'):
-            if self.mom_out != other.mom_out and self.mom_out and other.mom_out:
-                    raise Exception('mom_out does not match.')
-            mom_out = self.mom_out if self.mom_out else other.mom_out
-        else:
-            mom_out = self.mom_out
-
-        return self.__new__(Npr_matrix, super().__matmul__(other), mom_in, mom_out)
+        return self.__new__(Npr_matrix,
+                            super().__matmul__(other),
+                            self._propagate_mom(other, 'mom_in'),
+                            self._propagate_mom(other, 'mom_out'))
 
     def __array_finalize__(self, obj):
         if obj is None: return