corr2pt_Wilson_SF.cpp

Go to the documentation of this file.

#include "corr2pt_Wilson_SF.h"
#include "communicator.h"

using std::valarray;
using Bridge::vout;

//====================================================================

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

//====================================================================
double Corr2pt_Wilson_SF::fAfP(const valarray<Field_F>& sq1,
                               const valarray<Field_F>& sq2)
{
  int Lt = CommonParameters::Lt();

  valarray<dcomplex> mcorr(Lt);
  GammaMatrix        qn_src, qn_sink;

  int    Lvol = CommonParameters::Lvol();
  int    m_Nc = CommonParameters::Nc();
  double norm = 0.5 / Lvol * Lt;
  //  double norm=0.5/Lvol*Lt*(2*0.130*2*0.130);
  //  vout.general(m_vl,"norm=%lf\n",norm);

  vout.general(m_vl, "AO correlator:\n");

  qn_src  = m_gmset->get_GM(m_gmset->GAMMA5);
  qn_sink = m_gmset->get_GM(m_gmset->GAMMA54);

  meson_corr(mcorr, qn_sink, qn_src, sq1, sq1);

  for (int t = 0; t < mcorr.size(); ++t) {
    vout.general(m_vl, "fA  %4d  %20.12e  %20.12e\n",
                 t, -norm * real(mcorr[t]), -norm * imag(mcorr[t]));
  }


  vout.general(m_vl, "PO correlator:\n");

  qn_src  = m_gmset->get_GM(m_gmset->GAMMA5);
  qn_sink = m_gmset->get_GM(m_gmset->GAMMA5);
  //  qn_src  = m_gmset->get_GM(m_gmset->UNITY);
  //  qn_sink = m_gmset->get_GM(m_gmset->UNITY);

  meson_corr(mcorr, qn_sink, qn_src, sq1, sq1);

  for (int t = 0; t < mcorr.size(); ++t) {
    vout.general(m_vl, "fP  %4d  %20.12e  %20.12e\n",
                 t, norm * real(mcorr[t]), norm * imag(mcorr[t]));
  }


  vout.general(m_vl, "AO' correlator:\n");

  qn_src  = m_gmset->get_GM(m_gmset->GAMMA5);
  qn_sink = m_gmset->get_GM(m_gmset->GAMMA54);

  meson_corr(mcorr, qn_sink, qn_src, sq2, sq2);

  for (int t = 0; t < mcorr.size(); ++t) {
    vout.general(m_vl, "fA'  %4d  %20.12e  %20.12e\n",
                 t, norm * real(mcorr[Lt - 1 - t]), norm * imag(mcorr[Lt - 1 - t]));
  }


  vout.general(m_vl, "PO' correlator:\n");

  qn_src  = m_gmset->get_GM(m_gmset->GAMMA5);
  qn_sink = m_gmset->get_GM(m_gmset->GAMMA5);
  //  qn_src  = m_gmset->get_GM(m_gmset->UNITY);
  //  qn_sink = m_gmset->get_GM(m_gmset->UNITY);

  meson_corr(mcorr, qn_sink, qn_src, sq2, sq2);

  for (int t = 0; t < mcorr.size(); ++t) {
    vout.general(m_vl, "fP'  %4d  %20.12e  %20.12e\n",
                 t, norm * real(mcorr[Lt - 1 - t]), norm * imag(mcorr[Lt - 1 - t]));
  }


  //- NB. use meson correlator at t=1 for a non-trivial check.
  // double result = real(mcorr[0]);
  double result = real(mcorr[1]);

  return result;
}


//====================================================================
double Corr2pt_Wilson_SF::meson_corr(valarray<dcomplex>& meson,
                                     const GammaMatrix& qn_sink,
                                     const GammaMatrix& qn_src,
                                     const valarray<Field_F>& sq1,
                                     const valarray<Field_F>& sq2)
{
  int Nc   = CommonParameters::Nc();
  int Nd   = CommonParameters::Nd();
  int Nvol = CommonParameters::Nvol();
  int Lt   = CommonParameters::Lt();
  int Nx   = CommonParameters::Nx();
  int Ny   = CommonParameters::Ny();
  int Nz   = CommonParameters::Nz();
  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;
  valarray<dcomplex> corr_local(Nt);
  valarray<double>   corr_r(Nd), corr_i(Nd);
  valarray<int>      s2(Nd);

  Field corrF(2, Nvol, 2);

  corr_local = dcomplex(0.0);

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

      for (int t = 0; t < Nt; ++t) {
        corr_r = 0.0;
        corr_i = 0.0;
        for (int z = 0; z < Nz; ++z) {
          for (int y = 0; y < Ny; ++y) {
            for (int x = 0; x < Nx; ++x) {
              int site = index.site(x, y, z, t);

              for (int s0 = 0; s0 < Nd; ++s0) {
                int s1 = gm_sink.index(s0);
                //int s1 = qn_sink.index(s0);

                for (int c1 = 0; c1 < Nc; ++c1) {
                  corr_r[s0] += sq1[c0 + Nc * d0].cmp_r(c1, s1, site)
                                * sq2[c0 + Nc * d1].cmp_r(c1, s0, site)
                                + sq1[c0 + Nc * d0].cmp_i(c1, s1, site)
                                * sq2[c0 + Nc * d1].cmp_i(c1, s0, site);

                  corr_i[s0] += sq1[c0 + Nc * d0].cmp_r(c1, s1, site)
                                * sq2[c0 + Nc * d1].cmp_i(c1, s0, site)
                                - sq1[c0 + Nc * d0].cmp_i(c1, s1, site)
                                * sq2[c0 + Nc * d1].cmp_r(c1, s0, site);
                }
              }
            }
          }
        }
        for (int s0 = 0; s0 < Nd; ++s0) {
          dcomplex gmf = gm_src.value(d0) * gm_sink.value(s0);
          //dcomplex gmf = qn_src.value(d0)*qn_sink.value(s0);
          dcomplex corr = cmplx(corr_r[s0], corr_i[s0]);
          corr_local[t] += gmf * corr;
          //corr_local[t] += corr;
        }
        //vout.general(m_vl,"%d %lf %lf\n",t,real(corr_local[t]),imag(corr_local[t]));
      }
    }
  }

  valarray<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;
    //    vout.general(m_vl,"%d %d\n",t,t2);
    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);
    meson[t] = cmplx(crr, cri);
    //    vout.general(m_vl,"%d %lf %lf\n",t,crr,cri);
  }


  //- NB. use meson correlator at t=1 for a non-trivial check.
  // double result = real(mcorr[0]);
  double result = real(meson[1]);


  return result;
}


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