Bridge++  Ver. 1.1.x
 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 <cassert>
15 #include "corr2pt_4spinor.h"
16 
17 using std::valarray;
18 using Bridge::vout;
19 
20 //====================================================================
21 void Corr2pt_4spinor::setup()
22 {
23  assert(CommonParameters::Nc() == 3);
24 
25  int Nc = 3;
26  int n;
27 
28  m_epsilon_index.resize(Nc * 6);
29 
30  n = 0;
31  m_epsilon_index[Nc * n] = 0;
32  m_epsilon_index[1 + Nc * n] = 1;
33  m_epsilon_index[2 + Nc * n] = 2;
34  n = 1;
35  m_epsilon_index[Nc * n] = 1;
36  m_epsilon_index[1 + Nc * n] = 2;
37  m_epsilon_index[2 + Nc * n] = 0;
38  n = 2;
39  m_epsilon_index[Nc * n] = 2;
40  m_epsilon_index[1 + Nc * n] = 0;
41  m_epsilon_index[2 + Nc * n] = 1;
42  n = 3;
43  m_epsilon_index[Nc * n] = 2;
44  m_epsilon_index[1 + Nc * n] = 1;
45  m_epsilon_index[2 + Nc * n] = 0;
46  n = 4;
47  m_epsilon_index[Nc * n] = 1;
48  m_epsilon_index[1 + Nc * n] = 0;
49  m_epsilon_index[2 + Nc * n] = 2;
50  n = 5;
51  m_epsilon_index[Nc * n] = 0;
52  m_epsilon_index[1 + Nc * n] = 2;
53  m_epsilon_index[2 + Nc * n] = 1;
54 }
55 
56 
57 //====================================================================
58 double Corr2pt_4spinor::meson_all(const valarray<Field_F>& sq1,
59  const valarray<Field_F>& sq2)
60 {
61  int Lt = CommonParameters::Lt();
62 
63  valarray<dcomplex> mcorr(Lt);
64  GammaMatrix qn_src, qn_sink;
65 
66  vout.general(m_vl, "PS <-- PS correlator:\n");
67  qn_src = m_gmset->get_GM(m_gmset->GAMMA5);
68  qn_sink = m_gmset->get_GM(m_gmset->GAMMA5);
69  meson_corr(mcorr, qn_sink, qn_src, sq1, sq2);
70  for (int t = 0; t < mcorr.size(); ++t) {
71  vout.general(m_vl, " %4d %20.12e %20.12e\n",
72  t, real(mcorr[t]), imag(mcorr[t]));
73  }
74  double result = real(mcorr[0]);
75 
76  vout.general(m_vl, "V1 <-- V1 correlator:\n");
77  qn_src = m_gmset->get_GM(m_gmset->GAMMA1);
78  qn_sink = m_gmset->get_GM(m_gmset->GAMMA1);
79  meson_corr(mcorr, qn_sink, qn_src, sq1, sq2);
80  for (int t = 0; t < mcorr.size(); ++t) {
81  vout.general(m_vl, " %4d %20.12e %20.12e\n",
82  t, real(mcorr[t]), imag(mcorr[t]));
83  }
84 
85  vout.general(m_vl, "V2 <-- V2 correlator:\n");
86  qn_src = m_gmset->get_GM(m_gmset->GAMMA2);
87  qn_sink = m_gmset->get_GM(m_gmset->GAMMA2);
88  meson_corr(mcorr, qn_sink, qn_src, sq1, sq2);
89  for (int t = 0; t < mcorr.size(); ++t) {
90  vout.general(m_vl, " %4d %20.12e %20.12e\n",
91  t, real(mcorr[t]), imag(mcorr[t]));
92  }
93 
94  vout.general(m_vl, "V3 <-- V3 correlator:\n");
95  qn_src = m_gmset->get_GM(m_gmset->GAMMA3);
96  qn_sink = m_gmset->get_GM(m_gmset->GAMMA3);
97  meson_corr(mcorr, qn_sink, qn_src, sq1, sq2);
98  for (int t = 0; t < mcorr.size(); ++t) {
99  vout.general(m_vl, " %4d %20.12e %20.12e\n",
100  t, real(mcorr[t]), imag(mcorr[t]));
101  }
102 
103  return result;
104 }
105 
106 
107 //====================================================================
108 void Corr2pt_4spinor::meson_corr(valarray<dcomplex>& meson,
109  const GammaMatrix& qn_sink,
110  const GammaMatrix& qn_src,
111  const valarray<Field_F>& sq1,
112  const valarray<Field_F>& sq2)
113 {
114  int Nc = CommonParameters::Nc();
115  int Nd = CommonParameters::Nd();
116  int Lt = CommonParameters::Lt();
117  int Nt = CommonParameters::Nt();
118 
119  assert(meson.size() == Lt);
120 
121  GammaMatrix gm_src, gm_sink, gm5;
122  gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
123  gm_src = qn_src.mult(gm5);
124  gm_sink = gm5.mult(qn_sink);
125 
126  Index_lex index;
127  valarray<dcomplex> corr_local(Nt);
128  valarray<double> corr_r(Nd), corr_i(Nd);
129  valarray<int> s2(Nd);
130 
131  for (int c0 = 0; c0 < Nc; ++c0) {
132  for (int d0 = 0; d0 < Nd; ++d0) {
133  int d1 = gm_src.index(d0);
134 
135  dcomplex corr_t;
136  for (int t = 0; t < Nt; ++t) {
137  contract_at_t(corr_t, gm_sink,
138  sq1[c0 + Nc * d0], sq2[c0 + Nc * d1], t);
139 
140  dcomplex corr = cmplx(real(corr_t), imag(corr_t));
141  corr_local[t] += gm_src.value(d0) * corr;
142  }
143  }
144  }
145 
146  global_corr_t(meson, corr_local);
147 }
148 
149 
150 //====================================================================
151 double Corr2pt_4spinor::proton_test(const valarray<Field_F>& squ,
152  const valarray<Field_F>& sqd)
153 {
154  int Lt = CommonParameters::Lt();
155 
156  valarray<dcomplex> p_corr_unity(Lt), p_corr_gamm0(Lt), p_corr_upper(Lt);
157  GammaMatrix qn_unit, qn_gamm0;
158 
159 //#define DEBUG 0
160 
161 #if (DEBUG & 1)
162  vout.general(m_vl, "proton <-- proton correlator(UNITY):\n");
163 #endif
164 
165  qn_unit = m_gmset->get_GM(m_gmset->UNITY);
166  qn_gamm0 = m_gmset->get_GM(m_gmset->GAMMA4);
167 
168  proton_corr(p_corr_unity, qn_unit, squ, sqd);
169 
170 #if (DEBUG & 1)
171  for (int it = 0; it < p_corr_unity.size(); it++) {
172  vout.general(m_vl, " %4d %20.12e %20.12e\n",
173  it, real(p_corr_unity[it]), imag(p_corr_unity[it]));
174  }
175 #endif
176 
177  vout.general(m_vl, "proton <-- proton correlator(UPPER):\n");
178 
179  proton_corr(p_corr_gamm0, qn_gamm0, squ, sqd);
180  for (int it = 0; it < p_corr_upper.size(); it++) {
181  p_corr_upper[it] = (p_corr_unity[it] + p_corr_gamm0[it]) * 0.5;
182  vout.general(m_vl, " %4d %20.12e %20.12e\n",
183  it, real(p_corr_upper[it]), imag(p_corr_upper[it]));
184  }
185 
186 #if (DEBUG & 1)
187  vout.general(m_vl, "proton <-- proton correlator(GAMMA0):\n");
188 
189  for (int it = 0; it < p_corr_gamm0.size(); it++) {
190  vout.general(m_vl, " %4d %20.12e %20.12e\n",
191  it, real(p_corr_gamm0[it]), imag(p_corr_gamm0[it]));
192  }
193 #endif
194 
195  double result = real(p_corr_gamm0[0]);
196 
197  return result;
198 }
199 
200 
201 //====================================================================
202 void Corr2pt_4spinor::proton_corr(valarray<dcomplex>& proton,
203  const GammaMatrix& gm,
204  const valarray<Field_F>& squ,
205  const valarray<Field_F>& sqd)
206 {
207  int Nc = CommonParameters::Nc();
208  int Nd = CommonParameters::Nd();
209  int Lt = CommonParameters::Lt();
210  int Nt = CommonParameters::Nt();
211 
212  assert(Nc == 3);
213  assert(proton.size() == Lt);
214 
215  GammaMatrix cg5, c, gm5;
216  gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
217  c = m_gmset->get_GM(m_gmset->CHARGECONJG);
218  cg5 = c.mult(gm5);
219 
220 #if (DEBUG & 1)
221  vout.general(m_vl, "i:\tgm5\t\t\t\tc\t\t\t\tcg5\t\t\t\tgm\n");
222  for (int i = 0; i < Nd; i++) {
223  vout.general(m_vl, "%d:\t %d %e %e \t %d %e %e \t %d %e %e \t %d %e %e \n",
224  i,
225  gm5.index(i), real(gm5.value(i)), imag(gm5.value(i)),
226  c.index(i), real(c.value(i)), imag(c.value(i)),
227  cg5.index(i), real(cg5.value(i)), imag(cg5.value(i)),
228  gm.index(i), real(gm.value(i)), imag(gm.value(i))
229  );
230  }
231 #endif
232 
233  int FactNc = 6;
234  // This is valid only when Nc =3, which was already asserted.
235 
236  valarray<dcomplex> corr_local(Nt);
237  corr_local = cmplx(0.0, 0.0);
238 
239  for (int it = 0; it < Nt; it++) {
240 #if ((DEBUG & 65535) & 0)
241  vout.general(m_vl, "# it= %d\n", it);
242 #endif
243 
244  dcomplex sum = 0.0;
245  dcomplex sum1, sum2;
246 
247  for (int ialph = 0; ialph < Nd; ialph++) {
248  int ialphP = gm.index(ialph);
249  int ialph3 = ialph;
250  int ialph3P = ialphP;
251 
252  for (int ialph1P = 0; ialph1P < Nd; ialph1P++) {
253  int ialph2P = cg5.index(ialph1P);
254 
255  for (int ic123P = 0; ic123P < FactNc; ic123P++) {
256  int ic1P = epsilon_index(0, ic123P);
257  int ic2P = epsilon_index(1, ic123P);
258  int ic3P = epsilon_index(2, ic123P);
259  dcomplex factor = gm.value(ialph)
260  * cg5.value(ialph1P) * epsilon_value(ic123P);
261 
262  contract_at_t(sum1, cg5, ialph3,
263  squ[ic1P + Nc * ialph1P],
264  sqd[ic2P + Nc * ialph2P],
265  squ[ic3P + Nc * ialph3P], it);
266  contract_at_t(sum2, cg5, ialph3,
267  squ[ic3P + Nc * ialph3P],
268  sqd[ic2P + Nc * ialph2P],
269  squ[ic1P + Nc * ialph1P], it);
270  sum += factor * (sum1 - sum2);
271  }
272  }
273  }
274 
275  corr_local[it] = sum;
276  } // it loop end.
277 
278  global_corr_t(proton, corr_local);
279 }
280 
281 
282 //====================================================================
283 void Corr2pt_4spinor::global_corr_t(valarray<dcomplex>& corr_global,
284  valarray<dcomplex>& corr_local)
285 {
286  int Lt = CommonParameters::Lt();
287  int Nt = CommonParameters::Nt();
288 
289  assert(corr_global.size() == Lt);
290  assert(corr_local.size() == Nt);
291 
292  valarray<dcomplex> corr_tmp(Lt);
293 
294  int ipet = Communicator::ipe(3);
295 
296  for (int t = 0; t < Lt; ++t) {
297  corr_tmp[t] = cmplx(0.0, 0.0);
298  }
299 
300  for (int t = 0; t < Nt; ++t) {
301  int t2 = t + ipet * Nt;
302  corr_tmp[t2] = corr_local[t];
303  }
304 
305  for (int t = 0; t < Lt; ++t) {
306  double crr = real(corr_tmp[t]);
307  double cri = imag(corr_tmp[t]);
308  crr = Communicator::reduce_sum(crr);
309  cri = Communicator::reduce_sum(cri);
310  corr_global[t] = cmplx(crr, cri);
311  }
312 }
313 
314 
315 //====================================================================
316 //============================================================END=====