corr2pt_4spinor.cpp

Go to the documentation of this file.

#include "corr2pt_4spinor.h"

const std::string Corr2pt_4spinor::class_name = "Corr2pt_4spinor";

//====================================================================
void Corr2pt_4spinor::setup()
{
  assert(CommonParameters::Nc() == 3);

  int Nc = 3;
  int n;

  m_epsilon_index.resize(Nc * 6);

  n = 0;
  m_epsilon_index[Nc * n]     = 0;
  m_epsilon_index[1 + Nc * n] = 1;
  m_epsilon_index[2 + Nc * n] = 2;
  n = 1;
  m_epsilon_index[Nc * n]     = 1;
  m_epsilon_index[1 + Nc * n] = 2;
  m_epsilon_index[2 + Nc * n] = 0;
  n = 2;
  m_epsilon_index[Nc * n]     = 2;
  m_epsilon_index[1 + Nc * n] = 0;
  m_epsilon_index[2 + Nc * n] = 1;
  n = 3;
  m_epsilon_index[Nc * n]     = 2;
  m_epsilon_index[1 + Nc * n] = 1;
  m_epsilon_index[2 + Nc * n] = 0;
  n = 4;
  m_epsilon_index[Nc * n]     = 1;
  m_epsilon_index[1 + Nc * n] = 0;
  m_epsilon_index[2 + Nc * n] = 2;
  n = 5;
  m_epsilon_index[Nc * n]     = 0;
  m_epsilon_index[1 + Nc * n] = 2;
  m_epsilon_index[2 + Nc * n] = 1;
}


//====================================================================
double Corr2pt_4spinor::meson_all(const std::vector<Field_F>& sq1,
                                  const std::vector<Field_F>& sq2)
{
  int Lt = CommonParameters::Lt();

  std::vector<dcomplex> mcorr(Lt);
  GammaMatrix           qn_src, qn_sink;

  vout.general(m_vl, "PS <-- PS correlator:\n");
  qn_src  = m_gmset->get_GM(m_gmset->GAMMA5);
  qn_sink = m_gmset->get_GM(m_gmset->GAMMA5);
  meson_corr(mcorr, qn_sink, qn_src, sq1, sq2);
  for (int t = 0; t < mcorr.size(); ++t) {
    vout.general(m_vl, "  %4d  %20.12e  %20.12e\n",
                 t, real(mcorr[t]), imag(mcorr[t]));
  }
  double result = real(mcorr[0]);

  vout.general(m_vl, "V1 <-- V1 correlator:\n");
  qn_src  = m_gmset->get_GM(m_gmset->GAMMA1);
  qn_sink = m_gmset->get_GM(m_gmset->GAMMA1);
  meson_corr(mcorr, qn_sink, qn_src, sq1, sq2);
  for (int t = 0; t < mcorr.size(); ++t) {
    vout.general(m_vl, "  %4d  %20.12e  %20.12e\n",
                 t, real(mcorr[t]), imag(mcorr[t]));
  }

  vout.general(m_vl, "V2 <-- V2 correlator:\n");
  qn_src  = m_gmset->get_GM(m_gmset->GAMMA2);
  qn_sink = m_gmset->get_GM(m_gmset->GAMMA2);
  meson_corr(mcorr, qn_sink, qn_src, sq1, sq2);
  for (int t = 0; t < mcorr.size(); ++t) {
    vout.general(m_vl, "  %4d  %20.12e  %20.12e\n",
                 t, real(mcorr[t]), imag(mcorr[t]));
  }

  vout.general(m_vl, "V3 <-- V3 correlator:\n");
  qn_src  = m_gmset->get_GM(m_gmset->GAMMA3);
  qn_sink = m_gmset->get_GM(m_gmset->GAMMA3);
  meson_corr(mcorr, qn_sink, qn_src, sq1, sq2);
  for (int t = 0; t < mcorr.size(); ++t) {
    vout.general(m_vl, "  %4d  %20.12e  %20.12e\n",
                 t, real(mcorr[t]), imag(mcorr[t]));
  }

  return result;
}


//====================================================================
void Corr2pt_4spinor::meson_corr(std::vector<dcomplex>& meson,
                                 const GammaMatrix& qn_sink,
                                 const GammaMatrix& qn_src,
                                 const std::vector<Field_F>& sq1,
                                 const std::vector<Field_F>& sq2)
{
  int Nc = CommonParameters::Nc();
  int Nd = CommonParameters::Nd();
  int Lt = CommonParameters::Lt();
  int Nt = CommonParameters::Nt();

  assert(meson.size() == Lt);

  GammaMatrix gm_src, gm_sink, gm5;
  gm5     = m_gmset->get_GM(m_gmset->GAMMA5);
  gm_src  = qn_src.mult(gm5);
  gm_sink = gm5.mult(qn_sink);

  Index_lex             index;
  std::vector<dcomplex> corr_local(Nt);
  std::vector<double>   corr_r(Nd), corr_i(Nd);
  std::vector<int>      s2(Nd);

  for (int c0 = 0; c0 < Nc; ++c0) {
    for (int d0 = 0; d0 < Nd; ++d0) {
      int d1 = gm_src.index(d0);

      dcomplex corr_t;
      for (int t = 0; t < Nt; ++t) {
        contract_at_t(corr_t, gm_sink,
                      sq1[c0 + Nc * d0], sq2[c0 + Nc * d1], t);

        dcomplex corr = cmplx(real(corr_t), imag(corr_t));
        corr_local[t] += gm_src.value(d0) * corr;
      }
    }
  }

  global_corr_t(meson, corr_local);
}


//====================================================================
double Corr2pt_4spinor::proton_test(const std::vector<Field_F>& squ,
                                    const std::vector<Field_F>& sqd)
{
  int Lt = CommonParameters::Lt();

  std::vector<dcomplex> p_corr_unity(Lt), p_corr_gamm0(Lt), p_corr_upper(Lt);
  GammaMatrix           qn_unit, qn_gamm0;

//#define DEBUG 0

#if (DEBUG & 1)
  vout.general(m_vl, "proton <-- proton correlator(UNITY):\n");
#endif

  qn_unit  = m_gmset->get_GM(m_gmset->UNITY);
  qn_gamm0 = m_gmset->get_GM(m_gmset->GAMMA4);

  proton_corr(p_corr_unity, qn_unit, squ, sqd);

#if (DEBUG & 1)
  for (int it = 0; it < p_corr_unity.size(); it++) {
    vout.general(m_vl, "  %4d  %20.12e  %20.12e\n",
                 it, real(p_corr_unity[it]), imag(p_corr_unity[it]));
  }
#endif

  vout.general(m_vl, "proton <-- proton correlator(UPPER):\n");

  proton_corr(p_corr_gamm0, qn_gamm0, squ, sqd);
  for (int it = 0; it < p_corr_upper.size(); it++) {
    p_corr_upper[it] = (p_corr_unity[it] + p_corr_gamm0[it]) * 0.5;
    vout.general(m_vl, "  %4d  %20.12e  %20.12e\n",
                 it, real(p_corr_upper[it]), imag(p_corr_upper[it]));
  }

#if (DEBUG & 1)
  vout.general(m_vl, "proton <-- proton correlator(GAMMA0):\n");

  for (int it = 0; it < p_corr_gamm0.size(); it++) {
    vout.general(m_vl, "  %4d  %20.12e  %20.12e\n",
                 it, real(p_corr_gamm0[it]), imag(p_corr_gamm0[it]));
  }
#endif

  double result = real(p_corr_gamm0[0]);

  return result;
}


//====================================================================
void Corr2pt_4spinor::proton_corr(std::vector<dcomplex>& proton,
                                  const GammaMatrix& gm,
                                  const std::vector<Field_F>& squ,
                                  const std::vector<Field_F>& sqd)
{
  int Nc = CommonParameters::Nc();
  int Nd = CommonParameters::Nd();
  int Lt = CommonParameters::Lt();
  int Nt = CommonParameters::Nt();

  assert(Nc == 3);
  assert(proton.size() == Lt);

  GammaMatrix cg5, c, gm5;
  gm5 = m_gmset->get_GM(m_gmset->GAMMA5);
  c   = m_gmset->get_GM(m_gmset->CHARGECONJG);
  cg5 = c.mult(gm5);

#if (DEBUG & 1)
  vout.general(m_vl, "i:\tgm5\t\t\t\tc\t\t\t\tcg5\t\t\t\tgm\n");
  for (int i = 0; i < Nd; i++) {
    vout.general(m_vl, "%d:\t %d %e %e \t %d  %e %e \t %d  %e %e \t %d  %e %e \n",
                 i,
                 gm5.index(i), real(gm5.value(i)), imag(gm5.value(i)),
                 c.index(i), real(c.value(i)), imag(c.value(i)),
                 cg5.index(i), real(cg5.value(i)), imag(cg5.value(i)),
                 gm.index(i), real(gm.value(i)), imag(gm.value(i))
                 );
  }
#endif

  int FactNc = 6;
  // This is valid only when Nc =3, which was already asserted.

  std::vector<dcomplex> corr_local(Nt);
//  corr_local = cmplx(0.0, 0.0);

  for (int it = 0; it < Nt; it++) {
#if ((DEBUG & 65535) & 0)
    vout.general(m_vl, "# it= %d\n", it);
#endif

    dcomplex sum = 0.0;
    dcomplex sum1, sum2;

    for (int ialph = 0; ialph < Nd; ialph++) {
      int ialphP  = gm.index(ialph);
      int ialph3  = ialph;
      int ialph3P = ialphP;

      for (int ialph1P = 0; ialph1P < Nd; ialph1P++) {
        int ialph2P = cg5.index(ialph1P);

        for (int ic123P = 0; ic123P < FactNc; ic123P++) {
          int      ic1P   = epsilon_index(0, ic123P);
          int      ic2P   = epsilon_index(1, ic123P);
          int      ic3P   = epsilon_index(2, ic123P);
          dcomplex factor = gm.value(ialph)
                            * cg5.value(ialph1P) * epsilon_value(ic123P);

          contract_at_t(sum1, cg5, ialph3,
                        squ[ic1P + Nc * ialph1P],
                        sqd[ic2P + Nc * ialph2P],
                        squ[ic3P + Nc * ialph3P], it);
          contract_at_t(sum2, cg5, ialph3,
                        squ[ic3P + Nc * ialph3P],
                        sqd[ic2P + Nc * ialph2P],
                        squ[ic1P + Nc * ialph1P], it);
          sum += factor * (sum1 - sum2);
        }
      }
    }

    corr_local[it] = sum;
  } // it loop end.

  global_corr_t(proton, corr_local);
}


//====================================================================
void Corr2pt_4spinor::global_corr_t(std::vector<dcomplex>& corr_global,
                                    std::vector<dcomplex>& corr_local)
{
  int Lt = CommonParameters::Lt();
  int Nt = CommonParameters::Nt();

  assert(corr_global.size() == Lt);
  assert(corr_local.size() == Nt);

  std::vector<dcomplex> corr_tmp(Lt);

  int ipet = Communicator::ipe(3);

  for (int t = 0; t < Lt; ++t) {
    corr_tmp[t] = cmplx(0.0, 0.0);
  }

  for (int t = 0; t < Nt; ++t) {
    int t2 = t + ipet * Nt;
    corr_tmp[t2] = corr_local[t];
  }

  for (int t = 0; t < Lt; ++t) {
    double crr = real(corr_tmp[t]);
    double cri = imag(corr_tmp[t]);
    crr            = Communicator::reduce_sum(crr);
    cri            = Communicator::reduce_sum(cri);
    corr_global[t] = cmplx(crr, cri);
  }
}


//====================================================================
//============================================================END=====