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