Bridge++  Version 1.5.4
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
corr2pt_4spinor.cpp
Go to the documentation of this file.
1 
14 #include "corr2pt_4spinor.h"
15 
16 #include "Tools/epsilonTensor.h"
17 
18 const std::string Corr2pt_4spinor::class_name = "Corr2pt_4spinor";
19 
20 //====================================================================
21 void Corr2pt_4spinor::set_parameters(const Parameters& params)
22 {
23  m_filename_output = params.get_string("filename_output");
24  if (m_filename_output.empty()) {
25  m_filename_output = "stdout";
26  }
27 
28  const string str_vlevel = params.get_string("verbose_level");
29  m_vl = vout.set_verbose_level(str_vlevel);
30 }
31 
32 
33 //====================================================================
34 void Corr2pt_4spinor::init()
35 {
36  assert(CommonParameters::Nc() == 3);
37 
38  m_filename_output = "stdout";
39 }
40 
41 
42 //====================================================================
43 double Corr2pt_4spinor::meson_all(const std::vector<Field_F>& sq1,
44  const std::vector<Field_F>& sq2)
45 {
46  const int Lt = CommonParameters::Lt();
47 
48  std::ostream& log_file_previous = vout.getStream();
49  std::ofstream log_file;
50 
51  if (m_filename_output != "stdout") {
52  log_file.open(m_filename_output.c_str(), std::ios::app);
53  vout.init(log_file);
54  }
55 
56  std::vector<dcomplex> corr(Lt);
57 
58  GammaMatrix gm_src = m_gmset->get_GM(m_gmset->GAMMA5);
59  GammaMatrix gm_sink = m_gmset->get_GM(m_gmset->GAMMA5);
60  vout.general(m_vl, "PS <-- PS correlator:\n");
61  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
62  for (int t = 0; t < corr.size(); ++t) {
63  vout.general(m_vl, " %4d %20.12e %20.12e\n",
64  t, real(corr[t]), imag(corr[t]));
65  }
66  const double result = real(corr[0]);
67 
68  gm_src = m_gmset->get_GM(m_gmset->GAMMA1);
69  gm_sink = m_gmset->get_GM(m_gmset->GAMMA1);
70  vout.general(m_vl, "V1 <-- V1 correlator:\n");
71  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
72  for (int t = 0; t < corr.size(); ++t) {
73  vout.general(m_vl, " %4d %20.12e %20.12e\n",
74  t, real(corr[t]), imag(corr[t]));
75  }
76 
77  gm_src = m_gmset->get_GM(m_gmset->GAMMA2);
78  gm_sink = m_gmset->get_GM(m_gmset->GAMMA2);
79  vout.general(m_vl, "V2 <-- V2 correlator:\n");
80  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
81  for (int t = 0; t < corr.size(); ++t) {
82  vout.general(m_vl, " %4d %20.12e %20.12e\n",
83  t, real(corr[t]), imag(corr[t]));
84  }
85 
86  gm_src = m_gmset->get_GM(m_gmset->GAMMA3);
87  gm_sink = m_gmset->get_GM(m_gmset->GAMMA3);
88  vout.general(m_vl, "V3 <-- V3 correlator:\n");
89  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
90  for (int t = 0; t < corr.size(); ++t) {
91  vout.general(m_vl, " %4d %20.12e %20.12e\n",
92  t, real(corr[t]), imag(corr[t]));
93  }
94 
95  gm_src = m_gmset->get_GM(m_gmset->GAMMA5);
96  gm_sink = m_gmset->get_GM(m_gmset->GAMMA45);
97  vout.general(m_vl, "A4 <-- PS correlator:\n");
98  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
99  for (int t = 0; t < corr.size(); ++t) {
100  vout.general(m_vl, " %4d %20.12e %20.12e\n",
101  t, real(corr[t]), imag(corr[t]));
102  }
103 
104  gm_src = m_gmset->get_GM(m_gmset->GAMMA54);
105  gm_sink = m_gmset->get_GM(m_gmset->GAMMA5);
106  vout.general(m_vl, "PS <-- A4 correlator:\n");
107  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
108  for (int t = 0; t < corr.size(); ++t) {
109  vout.general(m_vl, " %4d %20.12e %20.12e\n",
110  t, real(corr[t]), imag(corr[t]));
111  }
112 
113  gm_src = m_gmset->get_GM(m_gmset->UNITY);
114  gm_sink = m_gmset->get_GM(m_gmset->UNITY);
115  vout.general(m_vl, "S <-- S correlator:\n");
116  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
117  for (int t = 0; t < corr.size(); ++t) {
118  vout.general(m_vl, " %4d %20.12e %20.12e\n",
119  t, real(corr[t]), imag(corr[t]));
120  }
121 
122  gm_src = m_gmset->get_GM(m_gmset->GAMMA51);
123  gm_sink = m_gmset->get_GM(m_gmset->GAMMA51);
124  vout.general(m_vl, "GAMMA51 <-- GAMMA51 correlator:\n");
125  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
126  for (int t = 0; t < corr.size(); ++t) {
127  vout.general(m_vl, " %4d %20.12e %20.12e\n",
128  t, real(corr[t]), imag(corr[t]));
129  }
130 
131  gm_src = m_gmset->get_GM(m_gmset->GAMMA52);
132  gm_sink = m_gmset->get_GM(m_gmset->GAMMA52);
133  vout.general(m_vl, "GAMMA52 <-- GAMMA52 correlator:\n");
134  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
135  for (int t = 0; t < corr.size(); ++t) {
136  vout.general(m_vl, " %4d %20.12e %20.12e\n",
137  t, real(corr[t]), imag(corr[t]));
138  }
139 
140  gm_src = m_gmset->get_GM(m_gmset->GAMMA53);
141  gm_sink = m_gmset->get_GM(m_gmset->GAMMA53);
142  vout.general(m_vl, "GAMMA53 <-- GAMMA53 correlator:\n");
143  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
144  for (int t = 0; t < corr.size(); ++t) {
145  vout.general(m_vl, " %4d %20.12e %20.12e\n",
146  t, real(corr[t]), imag(corr[t]));
147  }
148 
149  gm_src = m_gmset->get_GM(m_gmset->SIGMA12);
150  gm_sink = m_gmset->get_GM(m_gmset->SIGMA12);
151  vout.general(m_vl, "SIGMA12 <-- SIGMA12 correlator:\n");
152  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
153  for (int t = 0; t < corr.size(); ++t) {
154  vout.general(m_vl, " %4d %20.12e %20.12e\n",
155  t, real(corr[t]), imag(corr[t]));
156  }
157 
158  gm_src = m_gmset->get_GM(m_gmset->SIGMA23);
159  gm_sink = m_gmset->get_GM(m_gmset->SIGMA23);
160  vout.general(m_vl, "SIGMA23 <-- SIGMA23 correlator:\n");
161  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
162  for (int t = 0; t < corr.size(); ++t) {
163  vout.general(m_vl, " %4d %20.12e %20.12e\n",
164  t, real(corr[t]), imag(corr[t]));
165  }
166 
167  gm_src = m_gmset->get_GM(m_gmset->SIGMA31);
168  gm_sink = m_gmset->get_GM(m_gmset->SIGMA31);
169  vout.general(m_vl, "SIGMA31 <-- SIGMA31 correlator:\n");
170  meson_correlator(corr, gm_sink, gm_src, sq1, sq2);
171  for (int t = 0; t < corr.size(); ++t) {
172  vout.general(m_vl, " %4d %20.12e %20.12e\n",
173  t, real(corr[t]), imag(corr[t]));
174  }
175 
176 
177  if (m_filename_output != "stdout") {
178  log_file.close();
179  vout.init(log_file_previous);
180  }
181 
182  return result;
183 }
184 
185 
186 //====================================================================
187 void Corr2pt_4spinor::meson_correlator(std::vector<dcomplex>& corr_global,
188  const GammaMatrix& gm_sink,
189  const GammaMatrix& gm_src,
190  const std::vector<Field_F>& sq1,
191  const std::vector<Field_F>& sq2)
192 {
193  const int Nc = CommonParameters::Nc();
194  const int Nd = CommonParameters::Nd();
195  const int Lt = CommonParameters::Lt();
196  const int Nt = CommonParameters::Nt();
197 
198  assert(corr_global.size() == Lt);
199 
200  const GammaMatrix gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
201  const GammaMatrix gm_gm5_src = gm_src.mult(gm5);
202  const GammaMatrix gm5_gm_sink = gm5.mult(gm_sink);
203 
204  std::vector<dcomplex> corr_local(Nt, cmplx(0.0, 0.0));
205  for (int c0 = 0; c0 < Nc; ++c0) {
206  for (int d0 = 0; d0 < Nd; ++d0) {
207  int d1 = gm_gm5_src.index(d0);
208 
209  for (int t = 0; t < Nt; ++t) {
210  dcomplex corr_t;
211 
212  contract_at_t(corr_t, gm5_gm_sink,
213  sq1[c0 + Nc * d0], sq2[c0 + Nc * d1], t);
214 
215  corr_local[t] += gm_gm5_src.value(d0) * corr_t;
216  }
217  }
218  }
219  global_corr_t(corr_global, corr_local);
220 }
221 
222 
223 //====================================================================
224 double Corr2pt_4spinor::meson_momentum_all(const std::vector<Field_F>& sq1,
225  const std::vector<Field_F>& sq2,
226  const std::vector<int>& source_position)
227 {
228  const int Ndim = CommonParameters::Ndim();
229  const int Lt = CommonParameters::Lt();
230 
231  const int N_momentum = 10;
232 
233  typedef std::vector<int> MomentumSet;
234  std::vector<MomentumSet> momentum_sink(N_momentum);
235  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
236  momentum_sink[i_momentum].resize(Ndim - 1);
237  }
238 
239  //- momentum_sink[0] = (1,0,0)
240  int i_momentum = 0;
241  momentum_sink[i_momentum][0] = 1;
242  momentum_sink[i_momentum][1] = 0;
243  momentum_sink[i_momentum][2] = 0;
244 
245  //- momentum_sink[1] = (0,1,0)
246  i_momentum = 1;
247  momentum_sink[i_momentum][0] = 0;
248  momentum_sink[i_momentum][1] = 1;
249  momentum_sink[i_momentum][2] = 0;
250 
251  //- momentum_sink[2] = (0,0,1)
252  i_momentum = 2;
253  momentum_sink[i_momentum][0] = 0;
254  momentum_sink[i_momentum][1] = 0;
255  momentum_sink[i_momentum][2] = 1;
256 
257  //- momentum_sink[3] = (1,1,0)
258  i_momentum = 3;
259  momentum_sink[i_momentum][0] = 1;
260  momentum_sink[i_momentum][1] = 1;
261  momentum_sink[i_momentum][2] = 0;
262 
263  //- momentum_sink[4] = (0,1,1)
264  i_momentum = 4;
265  momentum_sink[i_momentum][0] = 0;
266  momentum_sink[i_momentum][1] = 1;
267  momentum_sink[i_momentum][2] = 1;
268 
269  //- momentum_sink[5] = (1,0,1)
270  i_momentum = 5;
271  momentum_sink[i_momentum][0] = 1;
272  momentum_sink[i_momentum][1] = 0;
273  momentum_sink[i_momentum][2] = 1;
274 
275  //- momentum_sink[6] = (1,1,1)
276  i_momentum = 6;
277  momentum_sink[i_momentum][0] = 1;
278  momentum_sink[i_momentum][1] = 1;
279  momentum_sink[i_momentum][2] = 1;
280 
281  //- momentum_sink[7] = (2,0,0)
282  i_momentum = 7;
283  momentum_sink[i_momentum][0] = 2;
284  momentum_sink[i_momentum][1] = 0;
285  momentum_sink[i_momentum][2] = 0;
286 
287  //- momentum_sink[8] = (0,2,0)
288  i_momentum = 8;
289  momentum_sink[i_momentum][0] = 0;
290  momentum_sink[i_momentum][1] = 2;
291  momentum_sink[i_momentum][2] = 0;
292 
293  //- momentum_sink[9] = (0,0,2)
294  i_momentum = 9;
295  momentum_sink[i_momentum][0] = 0;
296  momentum_sink[i_momentum][1] = 0;
297  momentum_sink[i_momentum][2] = 2;
298 
299 
300  std::ostream& log_file_previous = vout.getStream();
301  std::ofstream log_file;
302 
303  if (m_filename_output != "stdout") {
304  log_file.open(m_filename_output.c_str(), std::ios::app);
305  vout.init(log_file);
306  }
307 
308  std::vector<dcomplex> corr(Lt);
309 
310  GammaMatrix gm_src = m_gmset->get_GM(m_gmset->GAMMA5);
311  GammaMatrix gm_sink = m_gmset->get_GM(m_gmset->GAMMA5);
312  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
313  vout.general(m_vl, "PS_momentum(%d %d %d) <-- PS correlator:\n",
314  momentum_sink[i_momentum][0],
315  momentum_sink[i_momentum][1],
316  momentum_sink[i_momentum][2]);
317  meson_momentum_correlator(corr, momentum_sink[i_momentum], gm_sink, gm_src,
318  sq1, sq2, source_position);
319  for (int t = 0; t < corr.size(); ++t) {
320  vout.general(m_vl, " %4d %20.12e %20.12e\n",
321  t, real(corr[t]), imag(corr[t]));
322  }
323  }
324 
325  gm_src = m_gmset->get_GM(m_gmset->GAMMA1);
326  gm_sink = m_gmset->get_GM(m_gmset->GAMMA1);
327  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
328  vout.general(m_vl, "V1_momentum(%d %d %d) <-- V1 correlator:\n",
329  momentum_sink[i_momentum][0],
330  momentum_sink[i_momentum][1],
331  momentum_sink[i_momentum][2]);
332  meson_momentum_correlator(corr, momentum_sink[i_momentum], gm_sink, gm_src,
333  sq1, sq2, source_position);
334  for (int t = 0; t < corr.size(); ++t) {
335  vout.general(m_vl, " %4d %20.12e %20.12e\n",
336  t, real(corr[t]), imag(corr[t]));
337  }
338  }
339 
340  gm_src = m_gmset->get_GM(m_gmset->GAMMA2);
341  gm_sink = m_gmset->get_GM(m_gmset->GAMMA2);
342  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
343  vout.general(m_vl, "V2_momentum(%d %d %d) <-- V2 correlator:\n",
344  momentum_sink[i_momentum][0],
345  momentum_sink[i_momentum][1],
346  momentum_sink[i_momentum][2]);
347  meson_momentum_correlator(corr, momentum_sink[i_momentum], gm_sink, gm_src,
348  sq1, sq2, source_position);
349  for (int t = 0; t < corr.size(); ++t) {
350  vout.general(m_vl, " %4d %20.12e %20.12e\n",
351  t, real(corr[t]), imag(corr[t]));
352  }
353  }
354 
355  gm_src = m_gmset->get_GM(m_gmset->GAMMA3);
356  gm_sink = m_gmset->get_GM(m_gmset->GAMMA3);
357  for (int i_momentum = 0; i_momentum < N_momentum; i_momentum++) {
358  vout.general(m_vl, "V3_momentum(%d %d %d) <-- V3 correlator:\n",
359  momentum_sink[i_momentum][0],
360  momentum_sink[i_momentum][1],
361  momentum_sink[i_momentum][2]);
362  meson_momentum_correlator(corr, momentum_sink[i_momentum], gm_sink, gm_src,
363  sq1, sq2, source_position);
364  for (int t = 0; t < corr.size(); ++t) {
365  vout.general(m_vl, " %4d %20.12e %20.12e\n",
366  t, real(corr[t]), imag(corr[t]));
367  }
368  }
369 
370 
371  if (m_filename_output != "stdout") {
372  log_file.close();
373  vout.init(log_file_previous);
374  }
375 
376  return EXIT_SUCCESS;
377 }
378 
379 
380 //====================================================================
381 void Corr2pt_4spinor::meson_momentum_correlator(std::vector<dcomplex>& corr_global,
382  const std::vector<int>& momentum_sink,
383  const GammaMatrix& gm_sink,
384  const GammaMatrix& gm_src,
385  const std::vector<Field_F>& sq1,
386  const std::vector<Field_F>& sq2,
387  const std::vector<int>& source_position)
388 {
389  const int Nc = CommonParameters::Nc();
390  const int Nd = CommonParameters::Nd();
391  const int Lt = CommonParameters::Lt();
392  const int Nt = CommonParameters::Nt();
393 
394  assert(corr_global.size() == Lt);
395 
396  const GammaMatrix gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
397  const GammaMatrix gm_gm5_src = gm_src.mult(gm5);
398  const GammaMatrix gm5_gm_sink = gm5.mult(gm_sink);
399 
400  std::vector<dcomplex> corr_local(Nt, cmplx(0.0, 0.0));
401  for (int c0 = 0; c0 < Nc; ++c0) {
402  for (int d0 = 0; d0 < Nd; ++d0) {
403  int d1 = gm_gm5_src.index(d0);
404 
405  for (int t = 0; t < Nt; ++t) {
406  dcomplex corr_t;
407 
408  contract_at_t(corr_t, momentum_sink, gm5_gm_sink, source_position,
409  sq1[c0 + Nc * d0], sq2[c0 + Nc * d1], t);
410 
411  corr_local[t] += gm_gm5_src.value(d0) * corr_t;
412  }
413  }
414  }
415  global_corr_t(corr_global, corr_local);
416 }
417 
418 
419 //====================================================================
420 double Corr2pt_4spinor::proton_test(const std::vector<Field_F>& sq_u,
421  const std::vector<Field_F>& sq_d)
422 {
423  const int Lt = CommonParameters::Lt();
424 
425  std::ostream& log_file_previous = vout.getStream();
426  std::ofstream log_file;
427 
428  if (m_filename_output != "stdout") {
429  log_file.open(m_filename_output.c_str(), std::ios::app);
430  vout.init(log_file);
431  }
432 
433 
434  vout.general(m_vl, "proton <-- proton correlator(UNITY):\n");
435 
436  const GammaMatrix gm_unit = m_gmset->get_GM(m_gmset->UNITY);
437  const GammaMatrix gm_gamma0 = m_gmset->get_GM(m_gmset->GAMMA4);
438 
439  std::vector<dcomplex> p_corr_unity(Lt);
440  proton_correlator(p_corr_unity, gm_unit, sq_u, sq_d);
441 
442  for (int t = 0; t < p_corr_unity.size(); t++) {
443  vout.general(m_vl, " %4d %20.12e %20.12e\n",
444  t, real(p_corr_unity[t]), imag(p_corr_unity[t]));
445  }
446 
447  vout.general(m_vl, "proton <-- proton correlator(UPPER):\n");
448 
449  std::vector<dcomplex> p_corr_gamma0(Lt);
450  proton_correlator(p_corr_gamma0, gm_gamma0, sq_u, sq_d);
451 
452  std::vector<dcomplex> p_corr_upper(Lt);
453  for (int t = 0; t < p_corr_upper.size(); t++) {
454  p_corr_upper[t] = (p_corr_unity[t] + p_corr_gamma0[t]) * 0.5;
455  vout.general(m_vl, " %4d %20.12e %20.12e\n",
456  t, real(p_corr_upper[t]), imag(p_corr_upper[t]));
457  }
458 
459  vout.general(m_vl, "proton <-- proton correlator(GAMMA0):\n");
460 
461  for (int t = 0; t < p_corr_gamma0.size(); t++) {
462  vout.general(m_vl, " %4d %20.12e %20.12e\n",
463  t, real(p_corr_gamma0[t]), imag(p_corr_gamma0[t]));
464  }
465 
466 
467  if (m_filename_output != "stdout") {
468  log_file.close();
469  vout.init(log_file_previous);
470  }
471 
472  const double result = real(p_corr_gamma0[0]);
473 
474  return result;
475 }
476 
477 
478 //====================================================================
479 void Corr2pt_4spinor::proton_correlator(std::vector<dcomplex>& corr_global,
480  const GammaMatrix& gm,
481  const std::vector<Field_F>& sq_u,
482  const std::vector<Field_F>& sq_d)
483 {
484  const int Nc = CommonParameters::Nc();
485  const int Nd = CommonParameters::Nd();
486  const int Lt = CommonParameters::Lt();
487  const int Nt = CommonParameters::Nt();
488 
489  assert(Nc == 3);
490  assert(corr_global.size() == Lt);
491 
492  const GammaMatrix gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
493  const GammaMatrix c = m_gmset->get_GM(m_gmset->CHARGECONJG);
494  const GammaMatrix cg5 = c.mult(gm5);
495 
496 #ifdef DEBUG
497  vout.general(m_vl, "i:\tgm5\t\t\t\tc\t\t\t\tcg5\t\t\t\tgm\n");
498  for (int i = 0; i < Nd; i++) {
499  vout.general(m_vl, "%d:\t %d %e %e \t %d %e %e \t %d %e %e \t %d %e %e \n",
500  i,
501  gm5.index(i), real(gm5.value(i)), imag(gm5.value(i)),
502  c.index(i), real(c.value(i)), imag(c.value(i)),
503  cg5.index(i), real(cg5.value(i)), imag(cg5.value(i)),
504  gm.index(i), real(gm.value(i)), imag(gm.value(i))
505  );
506  }
507 #endif
508 
509  const int FactNc = 6;
510  // This is valid only when Nc =3, which was already asserted.
511 
512  std::vector<dcomplex> corr_local(Nt);
513 
514  for (int t = 0; t < Nt; t++) {
515  vout.paranoiac(m_vl, "# t= %d\n", t);
516 
517  dcomplex sum = cmplx(0.0, 0.0);
518  for (int i_alpha = 0; i_alpha < Nd; i_alpha++) {
519  int i_alphaP = gm.index(i_alpha);
520  int i_alpha3 = i_alpha;
521  int i_alpha3P = i_alphaP;
522 
523  for (int i_alpha1P = 0; i_alpha1P < Nd; i_alpha1P++) {
524  int i_alpha2P = cg5.index(i_alpha1P);
525 
526  for (int ic123P = 0; ic123P < FactNc; ic123P++) {
527  EpsilonTensor epsilon_tensor;
528 
529  int ic1P = epsilon_tensor.epsilon_3_index(ic123P, 0);
530  int ic2P = epsilon_tensor.epsilon_3_index(ic123P, 1);
531  int ic3P = epsilon_tensor.epsilon_3_index(ic123P, 2);
532  dcomplex factor = gm.value(i_alpha)
533  * cg5.value(i_alpha1P)
534  * static_cast<double>(epsilon_tensor.epsilon_3_value(ic123P));
535 
536  dcomplex sum1;
537  contract_at_t(sum1, cg5, i_alpha3,
538  sq_u[ic1P + Nc * i_alpha1P],
539  sq_d[ic2P + Nc * i_alpha2P],
540  sq_u[ic3P + Nc * i_alpha3P], t);
541 
542  dcomplex sum2;
543  contract_at_t(sum2, cg5, i_alpha3,
544  sq_u[ic3P + Nc * i_alpha3P],
545  sq_d[ic2P + Nc * i_alpha2P],
546  sq_u[ic1P + Nc * i_alpha1P], t);
547 
548  sum += factor * (sum1 - sum2);
549  }
550  }
551  }
552 
553  corr_local[t] = sum;
554  } // t loop end.
555 
556  global_corr_t(corr_global, corr_local);
557 }
558 
559 
560 //====================================================================
561 // meson =tr(gm5.qn_sink.sq1.qn_src.gm5.(sq2)^\dagger);
562 void Corr2pt_4spinor::meson_correlator_x(std::vector<dcomplex>& meson,
563  const GammaMatrix& qn_sink,
564  const GammaMatrix& qn_src,
565  const std::vector<Field_F>& sq1,
566  const std::vector<Field_F>& sq2)
567 {
568  int Nc = CommonParameters::Nc();
569  int Nd = CommonParameters::Nd();
570  int Lx = CommonParameters::Lx();
571  int Nx = CommonParameters::Nx();
572 
573  assert(meson.size() == Lx);
574 
575  GammaMatrix gm_src, gm_sink, gm5;
576  gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
577  gm_src = qn_src.mult(gm5);
578  gm_sink = gm5.mult(qn_sink);
579 
580  std::vector<dcomplex> corr_local(Nx);
581  for (int i = 0; i < Nx; ++i) {
582  corr_local[i] = 0.0;
583  }
584 
585  for (int c0 = 0; c0 < Nc; ++c0) {
586  for (int d0 = 0; d0 < Nd; ++d0) {
587  int d1 = gm_src.index(d0);
588 
589  for (int x = 0; x < Nx; ++x) {
590  dcomplex corr_x;
591  contract_at_x(corr_x, gm_sink,
592  sq1[c0 + Nc * d0], sq2[c0 + Nc * d1], x);
593 
594  corr_local[x] += gm_src.value(d0) * corr_x;
595  }
596  }
597  }
598 
599  global_corr_x(meson, corr_local);
600 }
601 
602 
603 //====================================================================
604 void Corr2pt_4spinor::meson_momentum_correlator_x(std::vector<dcomplex>& corr_global,
605  const std::vector<int>& momentum_sink,
606  const GammaMatrix& gm_sink,
607  const GammaMatrix& gm_src,
608  const std::vector<Field_F>& sq1,
609  const std::vector<Field_F>& sq2,
610  const std::vector<int>& source_position)
611 {
612  int Nc = CommonParameters::Nc();
613  int Nd = CommonParameters::Nd();
614  int Lx = CommonParameters::Lx();
615  int Nx = CommonParameters::Nx();
616 
617  assert(corr_global.size() == Lx);
618 
619  GammaMatrix gm_gm5_src, gm5_gm_sink, gm5;
620  gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
621  gm_gm5_src = gm_src.mult(gm5);
622  gm5_gm_sink = gm5.mult(gm_sink);
623 
624  std::vector<dcomplex> corr_local(Nx);
625 
626  for (int c0 = 0; c0 < Nc; ++c0) {
627  for (int d0 = 0; d0 < Nd; ++d0) {
628  int d1 = gm_gm5_src.index(d0);
629 
630  for (int x = 0; x < Nx; ++x) {
631  dcomplex corr_x;
632 
633  contract_at_x(corr_x, momentum_sink, gm5_gm_sink, source_position,
634  sq1[c0 + Nc * d0], sq2[c0 + Nc * d1], x);
635 
636  corr_local[x] += gm_gm5_src.value(d0) * corr_x;
637  }
638  }
639  }
640 
641  global_corr_x(corr_global, corr_local);
642 }
643 
644 
645 //====================================================================
646 void Corr2pt_4spinor::proton_correlator_x(std::vector<dcomplex>& proton,
647  const GammaMatrix& gm,
648  const std::vector<Field_F>& squ,
649  const std::vector<Field_F>& sqd)
650 {
651  int Nc = CommonParameters::Nc();
652  int Nd = CommonParameters::Nd();
653  int Lx = CommonParameters::Lx();
654  int Nx = CommonParameters::Nx();
655 
656  assert(Nc == 3);
657  assert(proton.size() == Lx);
658 
659  EpsilonTensor epsilon_tensor;
660 
661  GammaMatrix cg5, c, gm5;
662  gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
663  c = m_gmset->get_GM(m_gmset->CHARGECONJG);
664  cg5 = c.mult(gm5);
665 
666  int FactNc = 6;
667  // This is valid only when Nc =3, which was already asserted.
668 
669  std::vector<dcomplex> corr_local(Nx);
670 
671  for (int ix = 0; ix < Nx; ix++) {
672  dcomplex sum = 0.0;
673  dcomplex sum1, sum2;
674 
675  for (int ialph = 0; ialph < Nd; ialph++) {
676  int ialphP = gm.index(ialph);
677  int ialph3 = ialph;
678  int ialph3P = ialphP;
679 
680  for (int ialph1P = 0; ialph1P < Nd; ialph1P++) {
681  int ialph2P = cg5.index(ialph1P);
682 
683  for (int ic123P = 0; ic123P < FactNc; ic123P++) {
684  int ic1P = epsilon_tensor.epsilon_3_index(ic123P, 0);
685  int ic2P = epsilon_tensor.epsilon_3_index(ic123P, 1);
686  int ic3P = epsilon_tensor.epsilon_3_index(ic123P, 2);
687  dcomplex factor = gm.value(ialph)
688  * cg5.value(ialph1P)
689  * static_cast<double>(epsilon_tensor.epsilon_3_value(ic123P));
690 
691  contract_at_x(sum1, cg5, ialph3,
692  squ[ic1P + Nc * ialph1P],
693  sqd[ic2P + Nc * ialph2P],
694  squ[ic3P + Nc * ialph3P], ix);
695  contract_at_x(sum2, cg5, ialph3,
696  squ[ic3P + Nc * ialph3P],
697  sqd[ic2P + Nc * ialph2P],
698  squ[ic1P + Nc * ialph1P], ix);
699  sum += factor * (sum1 - sum2);
700  }
701  }
702  }
703 
704  corr_local[ix] = sum;
705  } // it loop end.
706 
707  global_corr_x(proton, corr_local);
708 }
709 
710 
711 /* moved by tanigchi
712 //====================================================================
713 void Corr2pt_4spinor::global_corr_t(std::vector<dcomplex>& corr_global,
714  const std::vector<dcomplex>& corr_local)
715 {
716  const int Lt = CommonParameters::Lt();
717  const int Nt = CommonParameters::Nt();
718 
719  assert(corr_global.size() == Lt);
720  assert(corr_local.size() == Nt);
721 
722  const int ipe_t = Communicator::ipe(3);
723 
724  std::vector<dcomplex> corr_tmp(Lt, cmplx(0.0, 0.0));
725 
726  for (int t = 0; t < Nt; ++t) {
727  int t_global = t + ipe_t * Nt;
728  corr_tmp[t_global] = corr_local[t];
729  }
730 
731  for (int t_global = 0; t_global < Lt; ++t_global) {
732  double cr_r = Communicator::reduce_sum(real(corr_tmp[t_global]));
733  double cr_i = Communicator::reduce_sum(imag(corr_tmp[t_global]));
734 
735  corr_global[t_global] = cmplx(cr_r, cr_i);
736  }
737 }
738 */
739 
740 //====================================================================
741 //============================================================END=====
Bridge::vout
BridgeIO vout
Definition: bridgeIO.cpp:503
Corr2pt_4spinor::meson_all
double meson_all(const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2)
Definition: corr2pt_4spinor.cpp:43
contract_at_t
void contract_at_t(dcomplex &corr, const GammaMatrix &gm_sink, const Field_F &v1, const Field_F &v2, const int time)
Contraction of hadron for 4-spinor fermion.
Definition: contract_4spinor.cpp:37
Bridge::BridgeIO::general
void general(const char *format,...)
Definition: bridgeIO.cpp:197
GammaMatrixSet::get_GM
GammaMatrix get_GM(GMspecies spec)
Definition: gammaMatrixSet.h:76
Corr2pt_4spinor::m_vl
Bridge::VerboseLevel m_vl
Definition: corr2pt_4spinor.h:48
Corr2pt_4spinor::proton_correlator
void proton_correlator(std::vector< dcomplex > &corr_global, const GammaMatrix &gm, const std::vector< Field_F > &sq_u, const std::vector< Field_F > &sq_d)
Definition: corr2pt_4spinor.cpp:479
Corr2pt_4spinor::m_gmset
GammaMatrixSet * m_gmset
Definition: corr2pt_4spinor.h:53
global_corr_x
void global_corr_x(std::vector< dcomplex > &corr_global, std::vector< dcomplex > &corr_local)
transform node-local correlator in x to global.
Definition: contract_4spinor.cpp:1632
Bridge::BridgeIO::init
void init(const std::string &filename)
Definition: bridgeIO.cpp:51
Parameters
Class for parameters.
Definition: parameters.h:46
Corr2pt_4spinor::proton_test
double proton_test(const std::vector< Field_F > &sq_u, const std::vector< Field_F > &sq_d)
Definition: corr2pt_4spinor.cpp:420
Corr2pt_4spinor::proton_correlator_x
void proton_correlator_x(std::vector< dcomplex > &proton, const GammaMatrix &gm, const std::vector< Field_F > &squ, const std::vector< Field_F > &sqd)
Definition: corr2pt_4spinor.cpp:646
GammaMatrix
Gamma Matrix class.
Definition: gammaMatrix.h:44
global_corr_t
void global_corr_t(std::vector< dcomplex > &corr_global, std::vector< dcomplex > &corr_local)
transform node-local correlator in t to global.
Definition: contract_4spinor.cpp:1713
corr2pt_4spinor.h
Corr2pt_4spinor::class_name
static const std::string class_name
Definition: corr2pt_4spinor.h:45
Corr2pt_4spinor::meson_correlator
void meson_correlator(std::vector< dcomplex > &corr_global, const GammaMatrix &gm_sink, const GammaMatrix &gm_src, const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2)
corr_global=(sq2)_{ab}(0,x) (gm_sink)_{bc} (sq1)_{cd}(x,0) (gm_src)_{da}=(sq2^*)_{ba}(x,0) (gamma_5 gm_sink)_{bc} (sq1)_{cd}(x,0) (gm_src gamma_5)_{da} , where sq1 and sq2 are quark propagators.
Definition: corr2pt_4spinor.cpp:187
EpsilonTensor
Epsilon tensor utility.
Definition: epsilonTensor.h:51
Corr2pt_4spinor::meson_momentum_correlator_x
void meson_momentum_correlator_x(std::vector< dcomplex > &corr_global, const std::vector< int > &momentum_sink, const GammaMatrix &gm_sink, const GammaMatrix &gm_src, const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2, const std::vector< int > &source_position)
Definition: corr2pt_4spinor.cpp:604
Bridge::BridgeIO::paranoiac
void paranoiac(const char *format,...)
Definition: bridgeIO.cpp:235
Corr2pt_4spinor::meson_correlator_x
void meson_correlator_x(std::vector< dcomplex > &meson, const GammaMatrix &gm_sink, const GammaMatrix &gm_src, const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2)
Definition: corr2pt_4spinor.cpp:562
GammaMatrix::mult
GammaMatrix mult(GammaMatrix) const
Definition: gammaMatrix.cpp:39
contract_at_x
void contract_at_x(dcomplex &corr, const GammaMatrix &gm_sink, const Field_F &v1, const Field_F &v2, int x)
Definition: contract_4spinor.cpp:477
Corr2pt_4spinor::meson_momentum_all
double meson_momentum_all(const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2, const std::vector< int > &source_position)
Definition: corr2pt_4spinor.cpp:224
Bridge::BridgeIO::getStream
std::ostream & getStream()
Definition: bridgeIO.cpp:391
Corr2pt_4spinor::m_filename_output
std::string m_filename_output
Definition: corr2pt_4spinor.h:51
GammaMatrix::value
dcomplex value(int row) const
Definition: gammaMatrix.h:88
Corr2pt_4spinor::set_parameters
virtual void set_parameters(const Parameters &params)
Definition: corr2pt_4spinor.cpp:21
Parameters::get_string
string get_string(const string &key) const
Definition: parameters.cpp:221
EpsilonTensor::epsilon_3_index
int epsilon_3_index(const int n, const int i) const
Definition: epsilonTensor.cpp:64
EpsilonTensor::epsilon_3_value
int epsilon_3_value(const int n) const
Definition: epsilonTensor.cpp:78
Corr2pt_4spinor::meson_momentum_correlator
void meson_momentum_correlator(std::vector< dcomplex > &corr_global, const std::vector< int > &momentum_sink, const GammaMatrix &gm_sink, const GammaMatrix &gm_src, const std::vector< Field_F > &sq1, const std::vector< Field_F > &sq2, const std::vector< int > &source_position)
Definition: corr2pt_4spinor.cpp:381
GammaMatrix::index
int index(int row) const
Definition: gammaMatrix.h:83
Bridge::BridgeIO::set_verbose_level
static VerboseLevel set_verbose_level(const std::string &str)
Definition: bridgeIO.cpp:131
Corr2pt_4spinor::init
void init()
Definition: corr2pt_4spinor.cpp:34