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