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