diff --git a/Include/scattering.h b/Include/scattering.h index f1fa8adae..3a8f0d76a 100644 --- a/Include/scattering.h +++ b/Include/scattering.h @@ -68,7 +68,8 @@ struct prop_p{ * @brief Bundle of meson_observables with momentum 0. */ struct mo_0{ - meson_observable *rho[3][3]; + // YD: Go to 6 to include g0gi + meson_observable *rho[6][6]; meson_observable *pi; }; @@ -78,7 +79,8 @@ struct mo_0{ struct mo_p{ int p[3]; meson_observable *d, *r1, *r2, *r3, *r4, *pi; - meson_observable *t1[3], *t2[3], *rho[3][3]; + // YD: Go to 6 to include g0gi + meson_observable *t1[6], *t2[6], *rho[6][6]; }; @@ -105,8 +107,10 @@ void io2pt(meson_observable* mo, int pmax, int sourceno, char* path, char* name, void io4pt(meson_observable* mo, int pmax, int sourceno, char* path, char* name,char * cnfg_filename); void io2pt_logfile(meson_observable* mo, int pmax, int sourceno, char* path, char* name,char * cnfg_filename); void io4pt_logfile(meson_observable* mo, int pmax, int sourceno, char* path, char* name,char * cnfg_filename); -void IOold_0(struct mo_0* molist[], int numsources, char* path, char* cnfg_filename, int pmax); -void IOold_p(struct mo_p* molist[], int numsources, char* path, char* cnfg_filename, int pmax); +void IO_0(struct mo_0* molist[], int numsources, char* path, char* cnfg_filename); +void IO_0_axial(struct mo_0* molist[], int numsources, char* path, char* cnfg_filename); +void IO_p(struct mo_p* molist[], int numsources, char* path, char* cnfg_filename, int pmax); +void IO_p_axial(struct mo_p* molist[], int numsources, char* path, char* cnfg_filename, int pmax); void IO_json_0(struct mo_0* molist[], int numsources, char* path,char * cnfg_filename); void IO_json_p(struct mo_p* molist[], int numsources, char* path, char* cnfg_filename); void init_mo_0(struct mo_0* mo, int pmax); diff --git a/LibHR/Observables/measure_scattering_tools.c b/LibHR/Observables/measure_scattering_tools.c index 1bea8c066..459541fb6 100644 --- a/LibHR/Observables/measure_scattering_tools.c +++ b/LibHR/Observables/measure_scattering_tools.c @@ -488,14 +488,15 @@ void init_mo_0(struct mo_0 *mo, int pmax) mo->pi = (meson_observable *)malloc(sizeof(meson_observable)); lprintf("init_mo_0", 0, "pi initiated\n"); init_mo(mo->pi, "pi", Nmom * GLB_T); - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) + // YD: Go to 6 to include g0gi + for (int i = 0; i < 6; i++) + for (int j = 0; j < 6; j++) { mo->rho[i][j] = (meson_observable *)malloc(sizeof(meson_observable)); init_mo(mo->rho[i][j], "rho", Nmom * GLB_T); mo->rho[i][j]->ind1 = i + 3; mo->rho[i][j]->ind2 = j + 3; - lprintf("init_mo_0", 0, "rho %d,%d initiated\n", i, j); + lprintf("init_mo_0", 0, "rho ind1=%d, ind2=%d initiated\n", i, j); } lprintf("MAIN", 0, "Meson observable 0 initiated!\n"); } @@ -526,7 +527,8 @@ void init_mo_p(struct mo_p *mo, int px, int py, int pz, int pmax) #undef X #undef R /// \endcond - for (int i = 0; i < 3; i++) + // YD: Go to 6 to include g0gi + for (int i = 0; i < 6; i++) { mo->t1[i] = (meson_observable *)malloc(sizeof(meson_observable)); init_mo(mo->t1[i], "t1", Nmom * GLB_T); @@ -535,15 +537,17 @@ void init_mo_p(struct mo_p *mo, int px, int py, int pz, int pmax) mo->t2[i] = (meson_observable *)malloc(sizeof(meson_observable)); init_mo(mo->t2[i], "t2", Nmom * GLB_T); mo->t2[i]->ind2 = i + 3; - lprintf("init_mo_p", 0, "Allocated t gamma=%d\n", i); + lprintf("init_mo_p", 0, "Allocated t ind2 = %d\n", i+3); - for (int j = 0; j < 3; j++) + // YD: Go to 6 to include g0gi + for (int j = 0; j < 6; j++) { mo->rho[i][j] = (meson_observable *)malloc(sizeof(meson_observable)); init_mo(mo->rho[i][j], "rho", Nmom * GLB_T); mo->rho[i][j]->ind1 = i + 3; mo->rho[i][j]->ind2 = j + 3; - lprintf("init_mo_p", 0, "Allocated rho gamma=%d,%d\n", i, j); + // YD: Renamed to make it more clear + lprintf("init_mo_p", 0, "Allocated rho ind1=%d, ind2 = %d\n", i+3, j+3); } } lprintf("MAIN", 0, "Meson observable p initiated!\n"); @@ -560,8 +564,9 @@ void gen_mo_0(struct mo_0 *mo, struct prop_common *p0, struct src_common *s0, in { measure_mesons_core(p0->Q_0, p0->Q_0, s0->src_0, mo->pi, 1, tau, pmax, 0, GLB_T); do_global_sum(mo->pi, 1.0); - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) + // YD: Go to 6 to include g0gi + for (int i = 0; i < 6; i++) + for (int j = 0; j < 6; j++) { measure_mesons_core(p0->Q_0, p0->Q_0, s0->src_0, mo->rho[i][j], 1, tau, pmax, 0, GLB_T); do_global_sum(mo->rho[i][j], 1.0); @@ -584,7 +589,8 @@ void gen_mo_p(struct mo_p *mo, struct prop_common *p0, struct prop_p *pp, struct do_global_sum(mo->pi, 1.0); measure_scattering_AD_core(mo->d, pp->Q_p, p0->Q_0, p0->Q_0_eta, p0->Q_0_eta, tau, 0, pmax-1, mo->p[0], mo->p[1], mo->p[2]); do_global_sum(mo->d, 1.0); - for (int i = 0; i < 3; i++) + // YD: Go to 6 to include g0gi + for (int i = 0; i < 6; i++) { measure_mesons_core(pp->W_0_mp, p0->Q_0, s0->src_0, mo->t1[i], 1, tau, pmax, 0, GLB_T); do_global_sum(mo->t1[i], 1.0); @@ -592,7 +598,8 @@ void gen_mo_p(struct mo_p *mo, struct prop_common *p0, struct prop_p *pp, struct measure_mesons_core(p0->Q_0, pp->W_0_p, s0->src_0, mo->t2[i], 1, tau, pmax, 0, GLB_T); do_global_sum(mo->t2[i], 1.0); - for (int j = 0; j < 3; j++) + // YD: Go to 6 to include g0gi + for (int j = 0; j < 6; j++) { measure_mesons_core(p0->Q_0, pp->Q_p, s0->src_0, mo->rho[i][j], 1, tau, pmax, 0, GLB_T); do_global_sum(mo->rho[i][j], 1.0); @@ -647,6 +654,7 @@ void gen_mo_p(struct mo_p *mo, struct prop_common *p0, struct prop_p *pp, struct /** * @brief Resets the mo_p structure, setting all correlation functions to 0. */ +// YD: Decide not th change anything here as I do not see it being used void reset_mo_p(struct mo_p *mo) { /// \cond @@ -693,8 +701,9 @@ void copy_mo(meson_observable *mo_dest, meson_observable *mo_src) void free_mo_0(struct mo_0 *mo) { free_mo(mo->pi); - for (int i = 0; i < 3; i++) - for (int j = 0; j < 3; j++) + // YD: Go to 6 to include g0gi + for (int i = 0; i < 6; i++) + for (int j = 0; j < 6; j++) { free_mo(mo->rho[i][j]); } @@ -720,11 +729,12 @@ void free_mo_p(struct mo_p *mo) #undef X #undef R /// \endcond - for (int i = 0; i < 3; i++) + // YD: Go to 6 to include g0gi + for (int i = 0; i < 6; i++) { free_mo(mo->t1[i]); free_mo(mo->t2[i]); - for (int j = 0; j < 3; j++) + for (int j = 0; j < 6; j++) { free_mo(mo->rho[i][j]); } @@ -840,31 +850,78 @@ void io4pt(meson_observable *mo, int pmax, int sourceno, char *path, char *name, } /** - * @brief "Old style" IO where each correlation function is saved to separate file. Prints zero-momentum only. + * @brief Prints zero-momentum only. * @param molist an array of mo_0 objects, where each index corresponds to a different noise source * @param numsources number of noise sources * @param path path to which the file should be saved * @param cnfg_filename name of the configuration * @see IO_json_0 */ -void IOold_0(struct mo_0 *molist[], int numsources, char *path, char *cnfg_filename, int pmax) +void IO_0(struct mo_0 *molist[], int numsources, char *path, char *cnfg_filename) { for (int src = 0; src < numsources; src++) { lprintf("IOold_0", 0, "Printing pi for source %d\n", src); - io2pt_logfile(molist[src]->pi,pmax, src, path, "pi", cnfg_filename); + io2pt_logfile(molist[src]->pi, 2, src, path, "pi", cnfg_filename); for (int i = 0; i < 3; i++) { - char tmp[100]; - lprintf("IOold_0", 0, "Printing rho %d for source %d\n", i + 1, src); - sprintf(tmp, "rho_p0_g%d", i + 1); - io2pt_logfile(molist[src]->rho[i][i], pmax, src, path, tmp, cnfg_filename); + for (int j = 0; j < 3; j++) + { + char tmp[100]; + lprintf("IOold_0", 0, "Printing rho %d for source %d\n", i + 1, src); + sprintf(tmp, "rho_p0_g%d_g%d", i + 1, j + 1); + io2pt_logfile(molist[src]->rho[i][j], 2, src, path, tmp, cnfg_filename); + } } } } /** - * @brief "Old style" IO where each correlation function is saved to separate file. Prints momentum p contractions. + * @brief Prints zero-momentum only for the contractions including g0gi operators. + * @param molist an array of mo_0 objects, where each index corresponds to a different noise source + * @param numsources number of noise sources + * @param path path to which the file should be saved + * @param cnfg_filename name of the configuration + * @see IO_json_0 + */ +void IO_0_axial(struct mo_0 *molist[], int numsources, char *path, char *cnfg_filename) +{ + for (int src = 0; src < numsources; src++) + { + lprintf("IOold_0", 0, "Printing pi for source %d\n", src); + io2pt_logfile(molist[src]->pi, 2, src, path, "pi", cnfg_filename); + for (int i = 0; i < 3; i++) + { + for (int j = 3; j < 6; j++) + { + char tmp[100]; + lprintf("IOold_0", 0, "Printing rho %d for source %d\n", i + 1, src); + sprintf(tmp, "rho_p0_g%d_g0g%d", i + 1, j - 2); + io2pt_logfile(molist[src]->rho[i][j], 2, src, path, tmp, cnfg_filename); + } + } + for (int i = 3; i < 6; i++) + { + for (int j = 0; j < 3; j++) + { + char tmp[100]; + lprintf("IOold_0", 0, "Printing rho %d for source %d\n", i + 1, src); + sprintf(tmp, "rho_p0_g0g%d_g%d", i - 2, j + 1); + io2pt_logfile(molist[src]->rho[i][j], 2, src, path, tmp, cnfg_filename); + } + for (int j = 3; j < 6; j++) + { + char tmp[100]; + lprintf("IOold_0", 0, "Printing rho %d for source %d\n", i + 1, src); + sprintf(tmp, "rho_p0_g0g%d_g0g%d", i - 2, j - 2); + io2pt_logfile(molist[src]->rho[i][j], 2, src, path, tmp, cnfg_filename); + } + } + } +} + +/** + * @brief Prints momentum p contractions. * @param molist an array of mo_p objects, where each index corresponds to a different noise source * @param numsources number of noise sources * @param path path to which the file should be saved @@ -872,7 +929,7 @@ void IOold_0(struct mo_0 *molist[], int numsources, char *path, char *cnfg_filen * * @see IO_json_p */ -void IOold_p(struct mo_p *molist[], int numsources, char *path, char *cnfg_filename, int pmax) +void IO_p(struct mo_p *molist[], int numsources, char *path, char *cnfg_filename, int pmax) { char tmp[100]; for (int src = 0; src < numsources; src++) @@ -901,15 +958,61 @@ void IOold_p(struct mo_p *molist[], int numsources, char *path, char *cnfg_filen for (int i = 0; i < 3; i++) { - lprintf("IOold_p", 0, "Printing for source %d momentum (%d,%d,%d) gamma %d\n", src, px, py, pz, i + 1); + lprintf("IOold_p", 0, "Printing T's and rho's for source %d momentum (%d,%d,%d) gamma %d\n", src, px, py, pz, i + 1); sprintf(tmp, "t1_p(%d,%d,%d)_g%d", px, py, pz, i + 1); io2pt_logfile(molist[src]->t1[i], pmax, src, path, tmp, cnfg_filename); sprintf(tmp, "t2_p(%d,%d,%d)_g%d", px, py, pz, i + 1); io2pt_logfile(molist[src]->t2[i], pmax, src, path, tmp, cnfg_filename); for (int j = 0; j < 3; j++) { - sprintf(tmp, "rho_p(%d,%d,%d)_g%d%d", px, py, pz, i + 1, j + 1); - lprintf("IOold_p", 0, "Printing for source %d momentum (%d,%d,%d) gamma %d\n", src, px, py, pz, i + 1, j + 1); + sprintf(tmp, "rho_p(%d,%d,%d)_g%d_g%d", px, py, pz, i + 1, j + 1); + io2pt_logfile(molist[src]->rho[i][j], pmax, src, path, tmp, cnfg_filename); + } + } + } +} + +/** + * @brief Prints momentum p contractions for the contractions including g0gi operators. + * @param molist an array of mo_p objects, where each index corresponds to a different noise source + * @param numsources number of noise sources + * @param path path to which the file should be saved + * @param cnfg_filename name of the configuration + * + * @see IO_json_p + */ +void IO_p_axial(struct mo_p *molist[], int numsources, char *path, char *cnfg_filename, int pmax) +{ + char tmp[100]; + for (int src = 0; src < numsources; src++) + { + int px = molist[src]->p[0]; + int py = molist[src]->p[1]; + int pz = molist[src]->p[2]; + + for (int i = 0; i < 3; i++) + { + for (int j = 3; j < 6; j++) + { + sprintf(tmp, "rho_p(%d,%d,%d)_g%d_g0g%d", px, py, pz, i + 1, j - 2); + io2pt_logfile(molist[src]->rho[i][j], pmax, src, path, tmp, cnfg_filename); + } + } + for (int i = 3; i < 6; i++) + { + lprintf("IOold_p", 0, "Printing T's and rho's for source %d momentum (%d,%d,%d) gamma 0 gamma %d\n", src, px, py, pz, i + 1); + sprintf(tmp, "t1_p(%d,%d,%d)_g0g%d", px, py, pz, i - 2); + io2pt_logfile(molist[src]->t1[i], pmax, src, path, tmp, cnfg_filename); + sprintf(tmp, "t2_p(%d,%d,%d)_g0g%d", px, py, pz, i - 2); + io2pt_logfile(molist[src]->t2[i], pmax, src, path, tmp, cnfg_filename); + for (int j = 0; j < 3; j++) + { + sprintf(tmp, "rho_p(%d,%d,%d)_g0g%d_g%d", px, py, pz, i - 2, j + 1); + io2pt_logfile(molist[src]->rho[i][j], pmax, src, path, tmp, cnfg_filename); + } + for (int j = 3; j < 6; j++) + { + sprintf(tmp, "rho_p(%d,%d,%d)_g0g%d_g0g%d", px, py, pz, i - 2, j - 2); io2pt_logfile(molist[src]->rho[i][j], pmax, src, path, tmp, cnfg_filename); } } diff --git a/Scattering/scatter.c b/Scattering/scatter.c index f3b5bba95..229db65e5 100644 --- a/Scattering/scatter.c +++ b/Scattering/scatter.c @@ -57,6 +57,7 @@ typedef struct _input_scatt int nhits; int tsrc; int free_theory; + int axial; char outdir[256], bc[16], p[256], configlist[256]; /* for the reading function */ @@ -77,6 +78,7 @@ typedef struct _input_scatt {"Boundary conditions:", "mes:bc = %s", STRING_T, &(varname).bc}, \ {"Momenta:", "mes:p = %s", STRING_T, &(varname).p}, \ {"Measure free theory:", "mes:free_theory = %d", INT_T, &(varname).free_theory}, \ + {"Measure Axialvector:", "mes:axial = %d", INT_T, &(varname).axial}, \ {NULL, NULL, INT_T, NULL} \ } \ } @@ -185,6 +187,7 @@ int main(int argc, char *argv[]) lprintf("MAIN", 0, "Configuration from %s\n", cnfg_filename); represent_gauge_field(); + // YD: Going to 6 here to include g0gi struct mo_0 *mo_p0[numsources]; struct mo_p *mo_p[Nmom][numsources]; for (int i = 0; i < numsources; i++) @@ -227,12 +230,20 @@ int main(int argc, char *argv[]) } } lprintf("MAIN", 0, "num sources: %d, path: %s\n", numsources, path); - IOold_0(mo_p0, numsources, path, cnfg_filename, pmax); + IO_0(mo_p0, numsources, path, cnfg_filename); + if (mes_var.axial) + { + IO_0_axial(mo_p0, numsources, path, cnfg_filename); + } //IO_json_0(mo_p0, numsources, path,cnfg_filename); for (int i = 0; i < Nmom; i++) { - IOold_p(mo_p[i], numsources, path, cnfg_filename, pmax); + IO_p(mo_p[i], numsources, path, cnfg_filename, pmax); //IO_json_p(mo_p[i], numsources, path,cnfg_filename); + if (mes_var.axial) + { + IO_p_axial(mo_p[i], numsources, path, cnfg_filename, pmax); + } } for (int src = 0; src < numsources; src++)