From 48fcc34d724386134b38cf57d280722fd348bd9f Mon Sep 17 00:00:00 2001
From: Antonin Portelli <antonin.portelli@me.com>
Date: Sun, 1 May 2016 18:31:40 -0700
Subject: [PATCH] CMeson: first implementation, still need proper output

---
 programs/Hadrons/CMeson.cc  | 31 +++++++++++++++++++++++++++++++
 programs/Hadrons/CMeson.hpp |  6 ++++++
 2 files changed, 37 insertions(+)

diff --git a/programs/Hadrons/CMeson.cc b/programs/Hadrons/CMeson.cc
index 537f157f..c171fb95 100644
--- a/programs/Hadrons/CMeson.cc
+++ b/programs/Hadrons/CMeson.cc
@@ -64,4 +64,35 @@ void CMeson::execute(Environment &env)
 {
     LOG(Message) << "computing meson contraction '" << getName() << "'"
                  << std::endl;
+    
+    XmlWriter             writer(par_.output);
+    LatticePropagator     &q1 = *env.getProp(par_.q1);
+    LatticePropagator     &q2 = *env.getProp(par_.q2);
+    LatticeComplex        c(env.getGrid());
+    SpinMatrix            g[Ns*Ns], g5;
+    std::vector<TComplex> buf;
+    Result                result;
+    unsigned int          nt = env.getGrid()->GlobalDimensions()[Tp];
+    
+    g5 = makeGammaProd(Ns*Ns - 1);
+    result.corr.resize(Ns*Ns);
+    for (unsigned int i = 0; i < Ns*Ns; ++i)
+    {
+        g[i] = makeGammaProd(i);
+    }
+    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);
 }
diff --git a/programs/Hadrons/CMeson.hpp b/programs/Hadrons/CMeson.hpp
index 8725d938..5e18cfbf 100644
--- a/programs/Hadrons/CMeson.hpp
+++ b/programs/Hadrons/CMeson.hpp
@@ -48,6 +48,12 @@ public:
                                         std::string, q2,
                                         std::string, output);
     };
+    class Result: Serializable
+    {
+    public:
+        GRID_SERIALIZABLE_CLASS_MEMBERS(Result,
+                        std::vector<std::vector<std::vector<Complex>>>, corr);
+    };
 public:
     // constructor
     CMeson(const std::string name);