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