Merge branch 'feature/hadrons' of https://github.com/paboyle/Grid into feature/rare_kaon

# Conflicts:
#	extras/Hadrons/Global.hpp
#	tests/hadrons/Test_hadrons_rarekaon.cc

extras/Hadrons/Modules/MContraction/Baryon.hpp

@@ -27,8 +27,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 *************************************************************************************/
 /*  END LEGAL */
 
-#ifndef Hadrons_Baryon_hpp_
-#define Hadrons_Baryon_hpp_
+#ifndef Hadrons_MContraction_Baryon_hpp_
+#define Hadrons_MContraction_Baryon_hpp_
 
 #include <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
@@ -55,9 +55,9 @@ template <typename FImpl1, typename FImpl2, typename FImpl3>
 class TBaryon: public Module<BaryonPar>
 {
 public:
-    TYPE_ALIASES(FImpl1, 1);
-    TYPE_ALIASES(FImpl2, 2);
-    TYPE_ALIASES(FImpl3, 3);
+    FERM_TYPE_ALIASES(FImpl1, 1);
+    FERM_TYPE_ALIASES(FImpl2, 2);
+    FERM_TYPE_ALIASES(FImpl3, 3);
     class Result: Serializable
     {
     public:
@@ -121,11 +121,11 @@ void TBaryon<FImpl1, FImpl2, FImpl3>::execute(void)
     
     // FIXME: do contractions
     
-    write(writer, "meson", result);
+    // write(writer, "meson", result);
 }
 
 END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_Baryon_hpp_
+#endif // Hadrons_MContraction_Baryon_hpp_

extras/Hadrons/Modules/MContraction/DiscLoop.hpp

@@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 *************************************************************************************/
 /*  END LEGAL */
 
-#ifndef Hadrons_DiscLoop_hpp_
-#define Hadrons_DiscLoop_hpp_
+#ifndef Hadrons_MContraction_DiscLoop_hpp_
+#define Hadrons_MContraction_DiscLoop_hpp_
 
 #include <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
@@ -52,7 +52,7 @@ public:
 template <typename FImpl>
 class TDiscLoop: public Module<DiscLoopPar>
 {
-    TYPE_ALIASES(FImpl,);
+    FERM_TYPE_ALIASES(FImpl,);
     class Result: Serializable
     {
     public:
@@ -141,4 +141,4 @@ END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_DiscLoop_hpp_
+#endif // Hadrons_MContraction_DiscLoop_hpp_

extras/Hadrons/Modules/MContraction/Gamma3pt.hpp

@@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 *************************************************************************************/
 /*  END LEGAL */
 
-#ifndef Hadrons_Gamma3pt_hpp_
-#define Hadrons_Gamma3pt_hpp_
+#ifndef Hadrons_MContraction_Gamma3pt_hpp_
+#define Hadrons_MContraction_Gamma3pt_hpp_
 
 #include <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
@@ -72,9 +72,9 @@ public:
 template <typename FImpl1, typename FImpl2, typename FImpl3>
 class TGamma3pt: public Module<Gamma3ptPar>
 {
-    TYPE_ALIASES(FImpl1, 1);
-    TYPE_ALIASES(FImpl2, 2);
-    TYPE_ALIASES(FImpl3, 3);
+    FERM_TYPE_ALIASES(FImpl1, 1);
+    FERM_TYPE_ALIASES(FImpl2, 2);
+    FERM_TYPE_ALIASES(FImpl3, 3);
     class Result: Serializable
     {
     public:
@@ -167,4 +167,4 @@ END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_Gamma3pt_hpp_
+#endif // Hadrons_MContraction_Gamma3pt_hpp_

extras/Hadrons/Modules/MContraction/Meson.hpp

@@ -29,8 +29,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 *************************************************************************************/
 /*  END LEGAL */
 
-#ifndef Hadrons_Meson_hpp_
-#define Hadrons_Meson_hpp_
+#ifndef Hadrons_MContraction_Meson_hpp_
+#define Hadrons_MContraction_Meson_hpp_
 
 #include <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
@@ -69,7 +69,7 @@ public:
                                     std::string, q1,
                                     std::string, q2,
                                     std::string, gammas,
-                                    std::string, mom,
+                                    std::string, sink,
                                     std::string, output);
 };
 
@@ -77,8 +77,10 @@ template <typename FImpl1, typename FImpl2>
 class TMeson: public Module<MesonPar>
 {
 public:
-    TYPE_ALIASES(FImpl1, 1);
-    TYPE_ALIASES(FImpl2, 2);
+    FERM_TYPE_ALIASES(FImpl1, 1);
+    FERM_TYPE_ALIASES(FImpl2, 2);
+    FERM_TYPE_ALIASES(ScalarImplCR, Scalar);
+    SINK_TYPE_ALIASES(Scalar);
     class Result: Serializable
     {
     public:
@@ -115,7 +117,7 @@ TMeson<FImpl1, FImpl2>::TMeson(const std::string name)
 template <typename FImpl1, typename FImpl2>
 std::vector<std::string> TMeson<FImpl1, FImpl2>::getInput(void)
 {
-    std::vector<std::string> input = {par().q1, par().q2};
+    std::vector<std::string> input = {par().q1, par().q2, par().sink};
     
     return input;
 }
@@ -131,12 +133,11 @@ std::vector<std::string> TMeson<FImpl1, FImpl2>::getOutput(void)
 template <typename FImpl1, typename FImpl2>
 void TMeson<FImpl1, FImpl2>::parseGammaString(std::vector<GammaPair> &gammaList)
 {
+    gammaList.clear();
     // Determine gamma matrices to insert at source/sink.
     if (par().gammas.compare("all") == 0)
     {
         // Do all contractions.
-        unsigned int n_gam = Ns * Ns;
-        gammaList.resize(n_gam*n_gam);
         for (unsigned int i = 1; i < Gamma::nGamma; i += 2)
         {
             for (unsigned int j = 1; j < Gamma::nGamma; j += 2)
@@ -155,6 +156,9 @@ void TMeson<FImpl1, FImpl2>::parseGammaString(std::vector<GammaPair> &gammaList)
 
 
 // execution ///////////////////////////////////////////////////////////////////
+#define mesonConnected(q1, q2, gSnk, gSrc) \
+(g5*(gSnk))*(q1)*(adj(gSrc)*g5)*adj(q2)
+
 template <typename FImpl1, typename FImpl2>
 void TMeson<FImpl1, FImpl2>::execute(void)
 {
@@ -162,43 +166,72 @@ void TMeson<FImpl1, FImpl2>::execute(void)
                  << " quarks '" << par().q1 << "' and '" << par().q2 << "'"
                  << std::endl;
     
-    CorrWriter              writer(par().output);
-    PropagatorField1       &q1 = *env().template getObject<PropagatorField1>(par().q1);
-    PropagatorField2       &q2 = *env().template getObject<PropagatorField2>(par().q2);
-    LatticeComplex         c(env().getGrid());
-    Gamma                  g5(Gamma::Algebra::Gamma5);
-    std::vector<GammaPair> gammaList;
+    CorrWriter             writer(par().output);
     std::vector<TComplex>  buf;
     std::vector<Result>    result;
-    std::vector<Real>      p;
-
-    p  = strToVec<Real>(par().mom);
-    LatticeComplex         ph(env().getGrid()), coor(env().getGrid());
-    Complex                i(0.0,1.0);
-    ph = zero;
-    for(unsigned int mu = 0; mu < env().getNd(); mu++)
-    {
-        LatticeCoordinate(coor, mu);
-        ph = ph + p[mu]*coor*((1./(env().getGrid()->_fdimensions[mu])));
-    }
-    ph = exp((Real)(2*M_PI)*i*ph);
+    Gamma                  g5(Gamma::Algebra::Gamma5);
+    std::vector<GammaPair> gammaList;
+    int                    nt = env().getDim(Tp);
     
     parseGammaString(gammaList);
-
     result.resize(gammaList.size());
     for (unsigned int i = 0; i < result.size(); ++i)
     {
-        Gamma gSnk(gammaList[i].first);
-        Gamma gSrc(gammaList[i].second);
-        c = trace((g5*gSnk)*q1*(adj(gSrc)*g5)*adj(q2))*ph;
-        sliceSum(c, buf, Tp);
-
         result[i].gamma_snk = gammaList[i].first;
         result[i].gamma_src = gammaList[i].second;
-        result[i].corr.resize(buf.size());
-        for (unsigned int t = 0; t < buf.size(); ++t)
+        result[i].corr.resize(nt);
+    }
+    if (env().template isObjectOfType<SlicedPropagator1>(par().q1) and
+        env().template isObjectOfType<SlicedPropagator2>(par().q2))
+    {
+        SlicedPropagator1 &q1 = *env().template getObject<SlicedPropagator1>(par().q1);
+        SlicedPropagator2 &q2 = *env().template getObject<SlicedPropagator2>(par().q2);
+        
+        LOG(Message) << "(propagator already sinked)" << std::endl;
+        for (unsigned int i = 0; i < result.size(); ++i)
         {
-            result[i].corr[t] = TensorRemove(buf[t]);
+            Gamma gSnk(gammaList[i].first);
+            Gamma gSrc(gammaList[i].second);
+            
+            for (unsigned int t = 0; t < buf.size(); ++t)
+            {
+                result[i].corr[t] = TensorRemove(trace(mesonConnected(q1[t], q2[t], gSnk, gSrc)));
+            }
+        }
+    }
+    else
+    {
+        PropagatorField1 &q1   = *env().template getObject<PropagatorField1>(par().q1);
+        PropagatorField2 &q2   = *env().template getObject<PropagatorField2>(par().q2);
+        LatticeComplex   c(env().getGrid());
+        
+        LOG(Message) << "(using sink '" << par().sink << "')" << std::endl;
+        for (unsigned int i = 0; i < result.size(); ++i)
+        {
+            Gamma       gSnk(gammaList[i].first);
+            Gamma       gSrc(gammaList[i].second);
+            std::string ns;
+                
+            ns = env().getModuleNamespace(env().getObjectModule(par().sink));
+            if (ns == "MSource")
+            {
+                PropagatorField1 &sink =
+                    *env().template getObject<PropagatorField1>(par().sink);
+                
+                c = trace(mesonConnected(q1, q2, gSnk, gSrc)*sink);
+                sliceSum(c, buf, Tp);
+            }
+            else if (ns == "MSink")
+            {
+                SinkFnScalar &sink = *env().template getObject<SinkFnScalar>(par().sink);
+                
+                c   = trace(mesonConnected(q1, q2, gSnk, gSrc));
+                buf = sink(c);
+            }
+            for (unsigned int t = 0; t < buf.size(); ++t)
+            {
+                result[i].corr[t] = TensorRemove(buf[t]);
+            }
         }
     }
     write(writer, "meson", result);
@@ -208,4 +241,4 @@ END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_Meson_hpp_
+#endif // Hadrons_MContraction_Meson_hpp_

extras/Hadrons/Modules/MContraction/WeakHamiltonian.hpp

@@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 *************************************************************************************/
 /*  END LEGAL */
 
-#ifndef Hadrons_WeakHamiltonian_hpp_
-#define Hadrons_WeakHamiltonian_hpp_
+#ifndef Hadrons_MContraction_WeakHamiltonian_hpp_
+#define Hadrons_MContraction_WeakHamiltonian_hpp_
 
 #include <Grid/Hadrons/Global.hpp>
 #include <Grid/Hadrons/Module.hpp>
@@ -83,7 +83,7 @@ public:
 class T##modname: public Module<WeakHamiltonianPar>\
 {\
 public:\
-    TYPE_ALIASES(FIMPL,)\
+    FERM_TYPE_ALIASES(FIMPL,)\
     class Result: Serializable\
     {\
     public:\
@@ -111,4 +111,4 @@ END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_WeakHamiltonian_hpp_
+#endif // Hadrons_MContraction_WeakHamiltonian_hpp_

extras/Hadrons/Modules/MContraction/WeakHamiltonianEye.hpp

@@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 *************************************************************************************/
 /*  END LEGAL */
 
-#ifndef Hadrons_WeakHamiltonianEye_hpp_
-#define Hadrons_WeakHamiltonianEye_hpp_
+#ifndef Hadrons_MContraction_WeakHamiltonianEye_hpp_
+#define Hadrons_MContraction_WeakHamiltonianEye_hpp_
 
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
 
@@ -55,4 +55,4 @@ END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_WeakHamiltonianEye_hpp_
+#endif // Hadrons_MContraction_WeakHamiltonianEye_hpp_

extras/Hadrons/Modules/MContraction/WeakHamiltonianNonEye.hpp

@@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 *************************************************************************************/
 /*  END LEGAL */
 
-#ifndef Hadrons_WeakHamiltonianNonEye_hpp_
-#define Hadrons_WeakHamiltonianNonEye_hpp_
+#ifndef Hadrons_MContraction_WeakHamiltonianNonEye_hpp_
+#define Hadrons_MContraction_WeakHamiltonianNonEye_hpp_
 
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
 
@@ -54,4 +54,4 @@ END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_WeakHamiltonianNonEye_hpp_
+#endif // Hadrons_MContraction_WeakHamiltonianNonEye_hpp_

extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp

@@ -26,8 +26,8 @@ See the full license in the file "LICENSE" in the top level distribution directo
 *************************************************************************************/
 /*  END LEGAL */
 
-#ifndef Hadrons_WeakNeutral4ptDisc_hpp_
-#define Hadrons_WeakNeutral4ptDisc_hpp_
+#ifndef Hadrons_MContraction_WeakNeutral4ptDisc_hpp_
+#define Hadrons_MContraction_WeakNeutral4ptDisc_hpp_
 
 #include <Grid/Hadrons/Modules/MContraction/WeakHamiltonian.hpp>
 
@@ -56,4 +56,4 @@ END_MODULE_NAMESPACE
 
 END_HADRONS_NAMESPACE
 
-#endif // Hadrons_WeakNeutral4ptDisc_hpp_
+#endif // Hadrons_MContraction_WeakNeutral4ptDisc_hpp_