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