diff --git a/docs/pyerrors/input/hadrons.html b/docs/pyerrors/input/hadrons.html index c5cef399..495c4a91 100644 --- a/docs/pyerrors/input/hadrons.html +++ b/docs/pyerrors/input/hadrons.html @@ -274,289 +274,295 @@ 172 # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case. 173 174 for diagram in diagrams: -175 real_data = np.zeros(Nt) -176 for x0 in range(Nt): -177 raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:]["re"].astype(np.double) -178 real_data += np.roll(raw_data, -x0) -179 real_data /= Nt +175 +176 if diagram == "triangle" and "Identity" not in str(identifier): +177 part = "im" +178 else: +179 part = "re" 180 -181 corr_data[diagram].append(real_data) -182 h5file.close() -183 -184 res_dict[str(identifier)] = {} -185 -186 for diagram in diagrams: -187 -188 tmp_data = np.array(corr_data[diagram]) +181 real_data = np.zeros(Nt) +182 for x0 in range(Nt): +183 raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:][part].astype(np.double) +184 real_data += np.roll(raw_data, -x0) +185 real_data /= Nt +186 +187 corr_data[diagram].append(real_data) +188 h5file.close() 189 -190 l_obs = [] -191 for c in tmp_data.T: -192 l_obs.append(Obs([c], [ens_id], idl=[idx])) +190 res_dict[str(identifier)] = {} +191 +192 for diagram in diagrams: 193 -194 corr = Corr(l_obs) -195 corr.tag = str(identifier) -196 -197 res_dict[str(identifier)][diagram] = corr -198 -199 return res_dict -200 -201 -202class Npr_matrix(np.ndarray): -203 -204 def __new__(cls, input_array, mom_in=None, mom_out=None): -205 obj = np.asarray(input_array).view(cls) -206 obj.mom_in = mom_in -207 obj.mom_out = mom_out -208 return obj +194 tmp_data = np.array(corr_data[diagram]) +195 +196 l_obs = [] +197 for c in tmp_data.T: +198 l_obs.append(Obs([c], [ens_id], idl=[idx])) +199 +200 corr = Corr(l_obs) +201 corr.tag = str(identifier) +202 +203 res_dict[str(identifier)][diagram] = corr +204 +205 return res_dict +206 +207 +208class Npr_matrix(np.ndarray): 209 -210 @property -211 def g5H(self): -212 """Gamma_5 hermitean conjugate -213 -214 Uses the fact that the propagator is gamma5 hermitean, so just the -215 in and out momenta of the propagator are exchanged. -216 """ -217 return Npr_matrix(self, -218 mom_in=self.mom_out, -219 mom_out=self.mom_in) -220 -221 def _propagate_mom(self, other, name): -222 s_mom = getattr(self, name, None) -223 o_mom = getattr(other, name, None) -224 if s_mom is not None and o_mom is not None: -225 if not np.allclose(s_mom, o_mom): -226 raise Exception(name + ' does not match.') -227 return o_mom if o_mom is not None else s_mom -228 -229 def __matmul__(self, other): -230 return self.__new__(Npr_matrix, -231 super().__matmul__(other), -232 self._propagate_mom(other, 'mom_in'), -233 self._propagate_mom(other, 'mom_out')) +210 def __new__(cls, input_array, mom_in=None, mom_out=None): +211 obj = np.asarray(input_array).view(cls) +212 obj.mom_in = mom_in +213 obj.mom_out = mom_out +214 return obj +215 +216 @property +217 def g5H(self): +218 """Gamma_5 hermitean conjugate +219 +220 Uses the fact that the propagator is gamma5 hermitean, so just the +221 in and out momenta of the propagator are exchanged. +222 """ +223 return Npr_matrix(self, +224 mom_in=self.mom_out, +225 mom_out=self.mom_in) +226 +227 def _propagate_mom(self, other, name): +228 s_mom = getattr(self, name, None) +229 o_mom = getattr(other, name, None) +230 if s_mom is not None and o_mom is not None: +231 if not np.allclose(s_mom, o_mom): +232 raise Exception(name + ' does not match.') +233 return o_mom if o_mom is not None else s_mom 234 -235 def __array_finalize__(self, obj): -236 if obj is None: -237 return -238 self.mom_in = getattr(obj, 'mom_in', None) -239 self.mom_out = getattr(obj, 'mom_out', None) +235 def __matmul__(self, other): +236 return self.__new__(Npr_matrix, +237 super().__matmul__(other), +238 self._propagate_mom(other, 'mom_in'), +239 self._propagate_mom(other, 'mom_out')) 240 -241 -242def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None): -243 """Read hadrons ExternalLeg hdf5 file and output an array of CObs -244 -245 Parameters -246 ---------- -247 path : str -248 path to the files to read -249 filestem : str -250 namestem of the files to read -251 ens_id : str -252 name of the ensemble, required for internal bookkeeping -253 idl : range -254 If specified only configurations in the given range are read in. -255 """ -256 -257 files, idx = _get_files(path, filestem, idl) -258 -259 mom = None -260 -261 corr_data = [] -262 for hd5_file in files: -263 file = h5py.File(path + '/' + hd5_file, "r") -264 raw_data = file['ExternalLeg/corr'][0][0].view('complex') -265 corr_data.append(raw_data) -266 if mom is None: -267 mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) -268 file.close() -269 corr_data = np.array(corr_data) -270 -271 rolled_array = np.rollaxis(corr_data, 0, 5) -272 -273 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) -274 for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]): -275 real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx]) -276 imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx]) -277 matrix[si, sj, ci, cj] = CObs(real, imag) +241 def __array_finalize__(self, obj): +242 if obj is None: +243 return +244 self.mom_in = getattr(obj, 'mom_in', None) +245 self.mom_out = getattr(obj, 'mom_out', None) +246 +247 +248def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None): +249 """Read hadrons ExternalLeg hdf5 file and output an array of CObs +250 +251 Parameters +252 ---------- +253 path : str +254 path to the files to read +255 filestem : str +256 namestem of the files to read +257 ens_id : str +258 name of the ensemble, required for internal bookkeeping +259 idl : range +260 If specified only configurations in the given range are read in. +261 """ +262 +263 files, idx = _get_files(path, filestem, idl) +264 +265 mom = None +266 +267 corr_data = [] +268 for hd5_file in files: +269 file = h5py.File(path + '/' + hd5_file, "r") +270 raw_data = file['ExternalLeg/corr'][0][0].view('complex') +271 corr_data.append(raw_data) +272 if mom is None: +273 mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) +274 file.close() +275 corr_data = np.array(corr_data) +276 +277 rolled_array = np.rollaxis(corr_data, 0, 5) 278 -279 return Npr_matrix(matrix, mom_in=mom) -280 -281 -282def read_Bilinear_hd5(path, filestem, ens_id, idl=None): -283 """Read hadrons Bilinear hdf5 file and output an array of CObs +279 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) +280 for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]): +281 real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx]) +282 imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx]) +283 matrix[si, sj, ci, cj] = CObs(real, imag) 284 -285 Parameters -286 ---------- -287 path : str -288 path to the files to read -289 filestem : str -290 namestem of the files to read -291 ens_id : str -292 name of the ensemble, required for internal bookkeeping -293 idl : range -294 If specified only configurations in the given range are read in. -295 """ -296 -297 files, idx = _get_files(path, filestem, idl) -298 -299 mom_in = None -300 mom_out = None -301 -302 corr_data = {} -303 for hd5_file in files: -304 file = h5py.File(path + '/' + hd5_file, "r") -305 for i in range(16): -306 name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8') -307 if name not in corr_data: -308 corr_data[name] = [] -309 raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex') -310 corr_data[name].append(raw_data) -311 if mom_in is None: -312 mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) -313 if mom_out is None: -314 mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float) -315 -316 file.close() -317 -318 result_dict = {} -319 -320 for key, data in corr_data.items(): -321 local_data = np.array(data) -322 -323 rolled_array = np.rollaxis(local_data, 0, 5) -324 -325 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) -326 for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]): -327 real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx]) -328 imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx]) -329 matrix[si, sj, ci, cj] = CObs(real, imag) +285 return Npr_matrix(matrix, mom_in=mom) +286 +287 +288def read_Bilinear_hd5(path, filestem, ens_id, idl=None): +289 """Read hadrons Bilinear hdf5 file and output an array of CObs +290 +291 Parameters +292 ---------- +293 path : str +294 path to the files to read +295 filestem : str +296 namestem of the files to read +297 ens_id : str +298 name of the ensemble, required for internal bookkeeping +299 idl : range +300 If specified only configurations in the given range are read in. +301 """ +302 +303 files, idx = _get_files(path, filestem, idl) +304 +305 mom_in = None +306 mom_out = None +307 +308 corr_data = {} +309 for hd5_file in files: +310 file = h5py.File(path + '/' + hd5_file, "r") +311 for i in range(16): +312 name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8') +313 if name not in corr_data: +314 corr_data[name] = [] +315 raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex') +316 corr_data[name].append(raw_data) +317 if mom_in is None: +318 mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) +319 if mom_out is None: +320 mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float) +321 +322 file.close() +323 +324 result_dict = {} +325 +326 for key, data in corr_data.items(): +327 local_data = np.array(data) +328 +329 rolled_array = np.rollaxis(local_data, 0, 5) 330 -331 result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out) -332 -333 return result_dict -334 -335 -336def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]): -337 """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs +331 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) +332 for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]): +333 real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx]) +334 imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx]) +335 matrix[si, sj, ci, cj] = CObs(real, imag) +336 +337 result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out) 338 -339 Parameters -340 ---------- -341 path : str -342 path to the files to read -343 filestem : str -344 namestem of the files to read -345 ens_id : str -346 name of the ensemble, required for internal bookkeeping -347 idl : range -348 If specified only configurations in the given range are read in. -349 vertices : list -350 Vertex functions to be extracted. -351 """ -352 -353 files, idx = _get_files(path, filestem, idl) -354 -355 mom_in = None -356 mom_out = None -357 -358 vertex_names = [] -359 for vertex in vertices: -360 vertex_names += _get_lorentz_names(vertex) -361 -362 corr_data = {} +339 return result_dict +340 +341 +342def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]): +343 """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs +344 +345 Parameters +346 ---------- +347 path : str +348 path to the files to read +349 filestem : str +350 namestem of the files to read +351 ens_id : str +352 name of the ensemble, required for internal bookkeeping +353 idl : range +354 If specified only configurations in the given range are read in. +355 vertices : list +356 Vertex functions to be extracted. +357 """ +358 +359 files, idx = _get_files(path, filestem, idl) +360 +361 mom_in = None +362 mom_out = None 363 -364 tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_' -365 -366 for hd5_file in files: -367 file = h5py.File(path + '/' + hd5_file, "r") -368 -369 for i in range(32): -370 name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8')) -371 if name in vertex_names: -372 if name not in corr_data: -373 corr_data[name] = [] -374 raw_data = file[tree + str(i) + '/corr'][0][0].view('complex') -375 corr_data[name].append(raw_data) -376 if mom_in is None: -377 mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) -378 if mom_out is None: -379 mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float) -380 -381 file.close() -382 -383 intermediate_dict = {} -384 -385 for vertex in vertices: -386 lorentz_names = _get_lorentz_names(vertex) -387 for v_name in lorentz_names: -388 if v_name in [('SigmaXY', 'SigmaZT'), -389 ('SigmaXT', 'SigmaYZ'), -390 ('SigmaYZ', 'SigmaXT'), -391 ('SigmaZT', 'SigmaXY')]: -392 sign = -1 -393 else: -394 sign = 1 -395 if vertex not in intermediate_dict: -396 intermediate_dict[vertex] = sign * np.array(corr_data[v_name]) -397 else: -398 intermediate_dict[vertex] += sign * np.array(corr_data[v_name]) -399 -400 result_dict = {} -401 -402 for key, data in intermediate_dict.items(): -403 -404 rolled_array = np.moveaxis(data, 0, 8) +364 vertex_names = [] +365 for vertex in vertices: +366 vertex_names += _get_lorentz_names(vertex) +367 +368 corr_data = {} +369 +370 tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_' +371 +372 for hd5_file in files: +373 file = h5py.File(path + '/' + hd5_file, "r") +374 +375 for i in range(32): +376 name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8')) +377 if name in vertex_names: +378 if name not in corr_data: +379 corr_data[name] = [] +380 raw_data = file[tree + str(i) + '/corr'][0][0].view('complex') +381 corr_data[name].append(raw_data) +382 if mom_in is None: +383 mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float) +384 if mom_out is None: +385 mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float) +386 +387 file.close() +388 +389 intermediate_dict = {} +390 +391 for vertex in vertices: +392 lorentz_names = _get_lorentz_names(vertex) +393 for v_name in lorentz_names: +394 if v_name in [('SigmaXY', 'SigmaZT'), +395 ('SigmaXT', 'SigmaYZ'), +396 ('SigmaYZ', 'SigmaXT'), +397 ('SigmaZT', 'SigmaXY')]: +398 sign = -1 +399 else: +400 sign = 1 +401 if vertex not in intermediate_dict: +402 intermediate_dict[vertex] = sign * np.array(corr_data[v_name]) +403 else: +404 intermediate_dict[vertex] += sign * np.array(corr_data[v_name]) 405 -406 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) -407 for index in np.ndindex(rolled_array.shape[:-1]): -408 real = Obs([rolled_array[index].real], [ens_id], idl=[idx]) -409 imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx]) -410 matrix[index] = CObs(real, imag) +406 result_dict = {} +407 +408 for key, data in intermediate_dict.items(): +409 +410 rolled_array = np.moveaxis(data, 0, 8) 411 -412 result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out) -413 -414 return result_dict -415 -416 -417def _get_lorentz_names(name): -418 lorentz_index = ['X', 'Y', 'Z', 'T'] +412 matrix = np.empty((rolled_array.shape[:-1]), dtype=object) +413 for index in np.ndindex(rolled_array.shape[:-1]): +414 real = Obs([rolled_array[index].real], [ens_id], idl=[idx]) +415 imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx]) +416 matrix[index] = CObs(real, imag) +417 +418 result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out) 419 -420 res = [] +420 return result_dict 421 -422 if name == "TT": -423 for i in range(4): -424 for j in range(i + 1, 4): -425 res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[i] + lorentz_index[j])) -426 return res +422 +423def _get_lorentz_names(name): +424 lorentz_index = ['X', 'Y', 'Z', 'T'] +425 +426 res = [] 427 -428 if name == "TTtilde": +428 if name == "TT": 429 for i in range(4): 430 for j in range(i + 1, 4): -431 for k in range(4): -432 for o in range(k + 1, 4): -433 fac = epsilon_tensor_rank4(i, j, k, o) -434 if not np.isclose(fac, 0.0): -435 res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[k] + lorentz_index[o])) -436 return res -437 -438 assert len(name) == 2 -439 -440 if 'S' in name or 'P' in name: -441 if not set(name) <= set(['S', 'P']): -442 raise Exception("'" + name + "' is not a Lorentz scalar") +431 res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[i] + lorentz_index[j])) +432 return res +433 +434 if name == "TTtilde": +435 for i in range(4): +436 for j in range(i + 1, 4): +437 for k in range(4): +438 for o in range(k + 1, 4): +439 fac = epsilon_tensor_rank4(i, j, k, o) +440 if not np.isclose(fac, 0.0): +441 res.append(("Sigma" + lorentz_index[i] + lorentz_index[j], "Sigma" + lorentz_index[k] + lorentz_index[o])) +442 return res 443 -444 g_names = {'S': 'Identity', -445 'P': 'Gamma5'} -446 -447 res.append((g_names[name[0]], g_names[name[1]])) -448 -449 else: -450 if not set(name) <= set(['V', 'A']): -451 raise Exception("'" + name + "' is not a Lorentz scalar") +444 assert len(name) == 2 +445 +446 if 'S' in name or 'P' in name: +447 if not set(name) <= set(['S', 'P']): +448 raise Exception("'" + name + "' is not a Lorentz scalar") +449 +450 g_names = {'S': 'Identity', +451 'P': 'Gamma5'} 452 -453 for ind in lorentz_index: -454 res.append(('Gamma' + ind + (name[0] == 'A') * 'Gamma5', -455 'Gamma' + ind + (name[1] == 'A') * 'Gamma5')) -456 -457 return res +453 res.append((g_names[name[0]], g_names[name[1]])) +454 +455 else: +456 if not set(name) <= set(['V', 'A']): +457 raise Exception("'" + name + "' is not a Lorentz scalar") +458 +459 for ind in lorentz_index: +460 res.append(('Gamma' + ind + (name[0] == 'A') * 'Gamma5', +461 'Gamma' + ind + (name[1] == 'A') * 'Gamma5')) +462 +463 return res @@ -729,31 +735,37 @@ If specified only configurations in the given range are read in. 173 # "DistillationContraction/Metadata/DmfSuffix" contains info about different quarks, irrelevant in the SU(3) case. 174 175 for diagram in diagrams: -176 real_data = np.zeros(Nt) -177 for x0 in range(Nt): -178 raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:]["re"].astype(np.double) -179 real_data += np.roll(raw_data, -x0) -180 real_data /= Nt +176 +177 if diagram == "triangle" and "Identity" not in str(identifier): +178 part = "im" +179 else: +180 part = "re" 181 -182 corr_data[diagram].append(real_data) -183 h5file.close() -184 -185 res_dict[str(identifier)] = {} -186 -187 for diagram in diagrams: -188 -189 tmp_data = np.array(corr_data[diagram]) +182 real_data = np.zeros(Nt) +183 for x0 in range(Nt): +184 raw_data = h5file["DistillationContraction/Correlators/" + diagram + "/" + str(x0)][:][part].astype(np.double) +185 real_data += np.roll(raw_data, -x0) +186 real_data /= Nt +187 +188 corr_data[diagram].append(real_data) +189 h5file.close() 190 -191 l_obs = [] -192 for c in tmp_data.T: -193 l_obs.append(Obs([c], [ens_id], idl=[idx])) +191 res_dict[str(identifier)] = {} +192 +193 for diagram in diagrams: 194 -195 corr = Corr(l_obs) -196 corr.tag = str(identifier) -197 -198 res_dict[str(identifier)][diagram] = corr -199 -200 return res_dict +195 tmp_data = np.array(corr_data[diagram]) +196 +197 l_obs = [] +198 for c in tmp_data.T: +199 l_obs.append(Obs([c], [ens_id], idl=[idx])) +200 +201 corr = Corr(l_obs) +202 corr.tag = str(identifier) +203 +204 res_dict[str(identifier)][diagram] = corr +205 +206 return res_dict @@ -786,44 +798,44 @@ If specified only configurations in the given range are read in. -
203class Npr_matrix(np.ndarray):
-204
-205    def __new__(cls, input_array, mom_in=None, mom_out=None):
-206        obj = np.asarray(input_array).view(cls)
-207        obj.mom_in = mom_in
-208        obj.mom_out = mom_out
-209        return obj
+            
209class Npr_matrix(np.ndarray):
 210
-211    @property
-212    def g5H(self):
-213        """Gamma_5 hermitean conjugate
-214
-215        Uses the fact that the propagator is gamma5 hermitean, so just the
-216        in and out momenta of the propagator are exchanged.
-217        """
-218        return Npr_matrix(self,
-219                          mom_in=self.mom_out,
-220                          mom_out=self.mom_in)
-221
-222    def _propagate_mom(self, other, name):
-223        s_mom = getattr(self, name, None)
-224        o_mom = getattr(other, name, None)
-225        if s_mom is not None and o_mom is not None:
-226            if not np.allclose(s_mom, o_mom):
-227                raise Exception(name + ' does not match.')
-228        return o_mom if o_mom is not None else s_mom
-229
-230    def __matmul__(self, other):
-231        return self.__new__(Npr_matrix,
-232                            super().__matmul__(other),
-233                            self._propagate_mom(other, 'mom_in'),
-234                            self._propagate_mom(other, 'mom_out'))
+211    def __new__(cls, input_array, mom_in=None, mom_out=None):
+212        obj = np.asarray(input_array).view(cls)
+213        obj.mom_in = mom_in
+214        obj.mom_out = mom_out
+215        return obj
+216
+217    @property
+218    def g5H(self):
+219        """Gamma_5 hermitean conjugate
+220
+221        Uses the fact that the propagator is gamma5 hermitean, so just the
+222        in and out momenta of the propagator are exchanged.
+223        """
+224        return Npr_matrix(self,
+225                          mom_in=self.mom_out,
+226                          mom_out=self.mom_in)
+227
+228    def _propagate_mom(self, other, name):
+229        s_mom = getattr(self, name, None)
+230        o_mom = getattr(other, name, None)
+231        if s_mom is not None and o_mom is not None:
+232            if not np.allclose(s_mom, o_mom):
+233                raise Exception(name + ' does not match.')
+234        return o_mom if o_mom is not None else s_mom
 235
-236    def __array_finalize__(self, obj):
-237        if obj is None:
-238            return
-239        self.mom_in = getattr(obj, 'mom_in', None)
-240        self.mom_out = getattr(obj, 'mom_out', None)
+236    def __matmul__(self, other):
+237        return self.__new__(Npr_matrix,
+238                            super().__matmul__(other),
+239                            self._propagate_mom(other, 'mom_in'),
+240                            self._propagate_mom(other, 'mom_out'))
+241
+242    def __array_finalize__(self, obj):
+243        if obj is None:
+244            return
+245        self.mom_in = getattr(obj, 'mom_in', None)
+246        self.mom_out = getattr(obj, 'mom_out', None)
 
@@ -1078,44 +1090,44 @@ in and out momenta of the propagator are exchanged.

-
243def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
-244    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
-245
-246    Parameters
-247    ----------
-248    path : str
-249        path to the files to read
-250    filestem : str
-251        namestem of the files to read
-252    ens_id : str
-253        name of the ensemble, required for internal bookkeeping
-254    idl : range
-255        If specified only configurations in the given range are read in.
-256    """
-257
-258    files, idx = _get_files(path, filestem, idl)
-259
-260    mom = None
-261
-262    corr_data = []
-263    for hd5_file in files:
-264        file = h5py.File(path + '/' + hd5_file, "r")
-265        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
-266        corr_data.append(raw_data)
-267        if mom is None:
-268            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
-269        file.close()
-270    corr_data = np.array(corr_data)
-271
-272    rolled_array = np.rollaxis(corr_data, 0, 5)
-273
-274    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
-275    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
-276        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
-277        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
-278        matrix[si, sj, ci, cj] = CObs(real, imag)
+            
249def read_ExternalLeg_hd5(path, filestem, ens_id, idl=None):
+250    """Read hadrons ExternalLeg hdf5 file and output an array of CObs
+251
+252    Parameters
+253    ----------
+254    path : str
+255        path to the files to read
+256    filestem : str
+257        namestem of the files to read
+258    ens_id : str
+259        name of the ensemble, required for internal bookkeeping
+260    idl : range
+261        If specified only configurations in the given range are read in.
+262    """
+263
+264    files, idx = _get_files(path, filestem, idl)
+265
+266    mom = None
+267
+268    corr_data = []
+269    for hd5_file in files:
+270        file = h5py.File(path + '/' + hd5_file, "r")
+271        raw_data = file['ExternalLeg/corr'][0][0].view('complex')
+272        corr_data.append(raw_data)
+273        if mom is None:
+274            mom = np.array(str(file['ExternalLeg/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
+275        file.close()
+276    corr_data = np.array(corr_data)
+277
+278    rolled_array = np.rollaxis(corr_data, 0, 5)
 279
-280    return Npr_matrix(matrix, mom_in=mom)
+280    matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
+281    for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
+282        real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
+283        imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
+284        matrix[si, sj, ci, cj] = CObs(real, imag)
+285
+286    return Npr_matrix(matrix, mom_in=mom)
 
@@ -1148,58 +1160,58 @@ If specified only configurations in the given range are read in.
-
283def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
-284    """Read hadrons Bilinear hdf5 file and output an array of CObs
-285
-286    Parameters
-287    ----------
-288    path : str
-289        path to the files to read
-290    filestem : str
-291        namestem of the files to read
-292    ens_id : str
-293        name of the ensemble, required for internal bookkeeping
-294    idl : range
-295        If specified only configurations in the given range are read in.
-296    """
-297
-298    files, idx = _get_files(path, filestem, idl)
-299
-300    mom_in = None
-301    mom_out = None
-302
-303    corr_data = {}
-304    for hd5_file in files:
-305        file = h5py.File(path + '/' + hd5_file, "r")
-306        for i in range(16):
-307            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
-308            if name not in corr_data:
-309                corr_data[name] = []
-310            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
-311            corr_data[name].append(raw_data)
-312            if mom_in is None:
-313                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
-314            if mom_out is None:
-315                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
-316
-317        file.close()
-318
-319    result_dict = {}
-320
-321    for key, data in corr_data.items():
-322        local_data = np.array(data)
-323
-324        rolled_array = np.rollaxis(local_data, 0, 5)
-325
-326        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
-327        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
-328            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
-329            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
-330            matrix[si, sj, ci, cj] = CObs(real, imag)
+            
289def read_Bilinear_hd5(path, filestem, ens_id, idl=None):
+290    """Read hadrons Bilinear hdf5 file and output an array of CObs
+291
+292    Parameters
+293    ----------
+294    path : str
+295        path to the files to read
+296    filestem : str
+297        namestem of the files to read
+298    ens_id : str
+299        name of the ensemble, required for internal bookkeeping
+300    idl : range
+301        If specified only configurations in the given range are read in.
+302    """
+303
+304    files, idx = _get_files(path, filestem, idl)
+305
+306    mom_in = None
+307    mom_out = None
+308
+309    corr_data = {}
+310    for hd5_file in files:
+311        file = h5py.File(path + '/' + hd5_file, "r")
+312        for i in range(16):
+313            name = file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['gamma'][0].decode('UTF-8')
+314            if name not in corr_data:
+315                corr_data[name] = []
+316            raw_data = file['Bilinear/Bilinear_' + str(i) + '/corr'][0][0].view('complex')
+317            corr_data[name].append(raw_data)
+318            if mom_in is None:
+319                mom_in = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
+320            if mom_out is None:
+321                mom_out = np.array(str(file['Bilinear/Bilinear_' + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
+322
+323        file.close()
+324
+325    result_dict = {}
+326
+327    for key, data in corr_data.items():
+328        local_data = np.array(data)
+329
+330        rolled_array = np.rollaxis(local_data, 0, 5)
 331
-332        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
-333
-334    return result_dict
+332        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
+333        for si, sj, ci, cj in np.ndindex(rolled_array.shape[:-1]):
+334            real = Obs([rolled_array[si, sj, ci, cj].real], [ens_id], idl=[idx])
+335            imag = Obs([rolled_array[si, sj, ci, cj].imag], [ens_id], idl=[idx])
+336            matrix[si, sj, ci, cj] = CObs(real, imag)
+337
+338        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
+339
+340    return result_dict
 
@@ -1232,85 +1244,85 @@ If specified only configurations in the given range are read in.
-
337def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
-338    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
-339
-340    Parameters
-341    ----------
-342    path : str
-343        path to the files to read
-344    filestem : str
-345        namestem of the files to read
-346    ens_id : str
-347        name of the ensemble, required for internal bookkeeping
-348    idl : range
-349        If specified only configurations in the given range are read in.
-350    vertices : list
-351        Vertex functions to be extracted.
-352    """
-353
-354    files, idx = _get_files(path, filestem, idl)
-355
-356    mom_in = None
-357    mom_out = None
-358
-359    vertex_names = []
-360    for vertex in vertices:
-361        vertex_names += _get_lorentz_names(vertex)
-362
-363    corr_data = {}
+            
343def read_Fourquark_hd5(path, filestem, ens_id, idl=None, vertices=["VA", "AV"]):
+344    """Read hadrons FourquarkFullyConnected hdf5 file and output an array of CObs
+345
+346    Parameters
+347    ----------
+348    path : str
+349        path to the files to read
+350    filestem : str
+351        namestem of the files to read
+352    ens_id : str
+353        name of the ensemble, required for internal bookkeeping
+354    idl : range
+355        If specified only configurations in the given range are read in.
+356    vertices : list
+357        Vertex functions to be extracted.
+358    """
+359
+360    files, idx = _get_files(path, filestem, idl)
+361
+362    mom_in = None
+363    mom_out = None
 364
-365    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
-366
-367    for hd5_file in files:
-368        file = h5py.File(path + '/' + hd5_file, "r")
-369
-370        for i in range(32):
-371            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
-372            if name in vertex_names:
-373                if name not in corr_data:
-374                    corr_data[name] = []
-375                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
-376                corr_data[name].append(raw_data)
-377                if mom_in is None:
-378                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
-379                if mom_out is None:
-380                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
-381
-382        file.close()
-383
-384    intermediate_dict = {}
-385
-386    for vertex in vertices:
-387        lorentz_names = _get_lorentz_names(vertex)
-388        for v_name in lorentz_names:
-389            if v_name in [('SigmaXY', 'SigmaZT'),
-390                          ('SigmaXT', 'SigmaYZ'),
-391                          ('SigmaYZ', 'SigmaXT'),
-392                          ('SigmaZT', 'SigmaXY')]:
-393                sign = -1
-394            else:
-395                sign = 1
-396            if vertex not in intermediate_dict:
-397                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
-398            else:
-399                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
-400
-401    result_dict = {}
-402
-403    for key, data in intermediate_dict.items():
-404
-405        rolled_array = np.moveaxis(data, 0, 8)
+365    vertex_names = []
+366    for vertex in vertices:
+367        vertex_names += _get_lorentz_names(vertex)
+368
+369    corr_data = {}
+370
+371    tree = 'FourQuarkFullyConnected/FourQuarkFullyConnected_'
+372
+373    for hd5_file in files:
+374        file = h5py.File(path + '/' + hd5_file, "r")
+375
+376        for i in range(32):
+377            name = (file[tree + str(i) + '/info'].attrs['gammaA'][0].decode('UTF-8'), file[tree + str(i) + '/info'].attrs['gammaB'][0].decode('UTF-8'))
+378            if name in vertex_names:
+379                if name not in corr_data:
+380                    corr_data[name] = []
+381                raw_data = file[tree + str(i) + '/corr'][0][0].view('complex')
+382                corr_data[name].append(raw_data)
+383                if mom_in is None:
+384                    mom_in = np.array(str(file[tree + str(i) + '/info'].attrs['pIn'])[3:-2].strip().split(), dtype=float)
+385                if mom_out is None:
+386                    mom_out = np.array(str(file[tree + str(i) + '/info'].attrs['pOut'])[3:-2].strip().split(), dtype=float)
+387
+388        file.close()
+389
+390    intermediate_dict = {}
+391
+392    for vertex in vertices:
+393        lorentz_names = _get_lorentz_names(vertex)
+394        for v_name in lorentz_names:
+395            if v_name in [('SigmaXY', 'SigmaZT'),
+396                          ('SigmaXT', 'SigmaYZ'),
+397                          ('SigmaYZ', 'SigmaXT'),
+398                          ('SigmaZT', 'SigmaXY')]:
+399                sign = -1
+400            else:
+401                sign = 1
+402            if vertex not in intermediate_dict:
+403                intermediate_dict[vertex] = sign * np.array(corr_data[v_name])
+404            else:
+405                intermediate_dict[vertex] += sign * np.array(corr_data[v_name])
 406
-407        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
-408        for index in np.ndindex(rolled_array.shape[:-1]):
-409            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
-410            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
-411            matrix[index] = CObs(real, imag)
+407    result_dict = {}
+408
+409    for key, data in intermediate_dict.items():
+410
+411        rolled_array = np.moveaxis(data, 0, 8)
 412
-413        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
-414
-415    return result_dict
+413        matrix = np.empty((rolled_array.shape[:-1]), dtype=object)
+414        for index in np.ndindex(rolled_array.shape[:-1]):
+415            real = Obs([rolled_array[index].real], [ens_id], idl=[idx])
+416            imag = Obs([rolled_array[index].imag], [ens_id], idl=[idx])
+417            matrix[index] = CObs(real, imag)
+418
+419        result_dict[key] = Npr_matrix(matrix, mom_in=mom_in, mom_out=mom_out)
+420
+421    return result_dict