Hadrons fixed for new gammas, Meson only does one contraction but this’ll change in the future

extras/Hadrons/Modules/MContraction/Meson.hpp

@@ -45,9 +45,11 @@ class MesonPar: Serializable
 {
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(MesonPar,
-                                    std::string, q1,
+                                    std::string,    q1,
-                                    std::string, q2,
+                                    std::string,    q2,
-                                    std::string, output);
+                                    std::string,    output,
+                                    Gamma::Algebra, gammaSource,
+                                    Gamma::Algebra, gammaSink);
 };
 
 template <typename FImpl1, typename FImpl2>
@@ -59,8 +61,7 @@ public:
     class Result: Serializable
     {
     public:
-        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result, std::vector<Complex>, corr);
-                                        std::vector<std::vector<std::vector<Complex>>>, corr);
     };
 public:
     // constructor
@@ -114,29 +115,17 @@ void TMeson<FImpl1, FImpl2>::execute(void)
     PropagatorField1      &q1 = *env().template getObject<PropagatorField1>(par().q1);
     PropagatorField2      &q2 = *env().template getObject<PropagatorField2>(par().q2);
     LatticeComplex        c(env().getGrid());
-    SpinMatrix            g[Ns*Ns], g5;
+    Gamma                 gSrc(par().gammaSource), gSnk(par().gammaSink);
+    Gamma                 g5(Gamma::Algebra::Gamma5);
     std::vector<TComplex> buf;
     Result                result;
     
-    //g5 = makeGammaProd(Ns*Ns - 1);
+    c = trace(gSnk*q1*adj(gSrc)*g5*adj(q2)*g5);
-    result.corr.resize(Ns*Ns);
+    sliceSum(c, buf, Tp);
-    for (unsigned int i = 0; i < Ns*Ns; ++i)
+    result.corr.resize(buf.size());
+    for (unsigned int t = 0; t < buf.size(); ++t)
     {
-        //g[i] = makeGammaProd(i);
+        result.corr[t] = TensorRemove(buf[t]);
-    }
-    for (unsigned int iSink = 0; iSink < Ns*Ns; ++iSink)
-    {
-        result.corr[iSink].resize(Ns*Ns);
-        for (unsigned int iSrc = 0; iSrc < Ns*Ns; ++iSrc)
-        {
-            c = trace(g[iSink]*q1*g[iSrc]*g5*adj(q2)*g5);
-            sliceSum(c, buf, Tp);
-            result.corr[iSink][iSrc].resize(buf.size());
-            for (unsigned int t = 0; t < buf.size(); ++t)
-            {
-                result.corr[iSink][iSrc][t] = TensorRemove(buf[t]);
-            }
-        }
     }
     write(writer, "meson", result);
 }

extras/Hadrons/Modules/MSource/SeqGamma.hpp

@@ -60,11 +60,11 @@ class SeqGammaPar: Serializable
 {
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(SeqGammaPar,
-                                    std::string,  q,
+                                    std::string,    q,
-                                    unsigned int, tA,
+                                    unsigned int,   tA,
-                                    unsigned int, tB,
+                                    unsigned int,   tB,
-                                    unsigned int, gamma,
+                                    Gamma::Algebra, gamma,
-                                    std::string,  mom);
+                                    std::string,    mom);
 };
 
 template <typename FImpl>
@@ -140,11 +140,10 @@ void TSeqGamma<FImpl>::execute(void)
     PropagatorField &q   = *env().template getObject<PropagatorField>(par().q);
     Lattice<iScalar<vInteger>> t(env().getGrid());
     LatticeComplex             ph(env().getGrid()), coor(env().getGrid());
-    SpinMatrix                 g;
+    Gamma                      g(par().gamma);
     std::vector<Real>          p;
     Complex                    i(0.0,1.0);
     
-    //g  = makeGammaProd(par().gamma);
     p  = strToVec<Real>(par().mom);
     ph = zero;
     for(unsigned int mu = 0; mu < env().getNd(); mu++)
@@ -154,7 +153,7 @@ void TSeqGamma<FImpl>::execute(void)
     }
     ph = exp(i*ph);
     LatticeCoordinate(t, Tp);
-    src = where((t >= par().tA) and (t <= par().tB), g*ph*q, 0.*q);
+    src = where((t >= par().tA) and (t <= par().tB), ph*(g*q), 0.*q);
 }
 
 END_MODULE_NAMESPACE

tests/hadrons/Test_hadrons_meson_3pt.cc

@@ -30,6 +30,14 @@
 using namespace Grid;
 using namespace Hadrons;
 
+static Gamma::Algebra gmu[4] =
+{
+    Gamma::Algebra::GammaX,
+    Gamma::Algebra::GammaY,
+    Gamma::Algebra::GammaZ,
+    Gamma::Algebra::GammaT
+}
+
 int main(int argc, char *argv[])
 {
     // initialization //////////////////////////////////////////////////////////
@@ -102,7 +110,7 @@ int main(int argc, char *argv[])
             seqName.push_back(std::vector<std::string>(Nd));
             for (unsigned int mu = 0; mu < Nd; ++mu)
             {
-                seqPar.gamma   = 0x1 << mu;
+                seqPar.gamma   = gmu[mu];
                 seqName[i][mu] = "G" + std::to_string(seqPar.gamma)
                                  + "_" + std::to_string(seqPar.tA) + "-"
                                  + qName[i];
@@ -127,9 +135,11 @@ int main(int argc, char *argv[])
         for (unsigned int i = 0; i < flavour.size(); ++i)
         for (unsigned int j = i; j < flavour.size(); ++j)
         {
-            mesPar.output = "mesons/Z2_" + flavour[i] + flavour[j];
+            mesPar.output      = "mesons/Z2_" + flavour[i] + flavour[j];
-            mesPar.q1     = qName[i];
+            mesPar.q1          = qName[i];
-            mesPar.q2     = qName[j];
+            mesPar.q2          = qName[j];
+            mesPar.gammaSource = Gamma::Algebra::Gamma5;
+            mesPar.gammaSink   = Gamma::Algebra::Gamma5;
             application.createModule<MContraction::Meson>("meson_Z2_"
                                                           + std::to_string(t)
                                                           + "_"