diff --git a/Grid/qcd/utils/BaryonUtils.h b/Grid/qcd/utils/BaryonUtils.h index 25c91342..d65b9176 100644 --- a/Grid/qcd/utils/BaryonUtils.h +++ b/Grid/qcd/utils/BaryonUtils.h @@ -41,7 +41,10 @@ public: typedef typename FImpl::SitePropagator pobj; typedef typename ComplexField::vector_object vobj; - + + typedef Lattice> SpinMatrixField; + typedef typename SpinMatrixField::vector_object sobj; + static const int epsilon[6][3] ; static const Complex epsilon_sgn[6]; @@ -81,6 +84,69 @@ public: const char * quarks_right, const int parity, robj &result); + private: + template + static void Sigma_to_Nucleon_Q1_Eye_site(const mobj &Dq_loop, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); + template + static void Sigma_to_Nucleon_Q1_NonEye_site(const mobj &Du_ti, + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); + + + template + static void Sigma_to_Nucleon_Q2_Eye_site(const mobj &Dq_loop, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); + template + static void Sigma_to_Nucleon_Q2_NonEye_site(const mobj &Du_ti, + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result); + public: + template + static void Sigma_to_Nucleon_Eye(const PropagatorField &qq_loop, + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr); + template + static void Sigma_to_Nucleon_NonEye(const PropagatorField &qq_ti, + const PropagatorField &qq_tf, + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr); }; template @@ -254,4 +320,305 @@ void BaryonUtils::ContractBaryons_Sliced(const mobj &D1, result=Zero(); baryon_site(D1,D2,D3,GammaA_left,GammaB_left,GammaA_right,GammaB_right,parity,wick_contraction,result); } + +/*********************************************************************** + * End of Baryon 2pt-function code. * + * * + * The following code is for Sigma -> N rare hypeon decays * + **********************************************************************/ + +/* Dq_loop is a quark line from t_H to t_H + * Du_spec is a quark line from t_i to t_f + * Dd_tf is a quark line from t_f to t_H + * Ds_ti is a quark line from t_i to t_H */ +template +template +void BaryonUtils::Sigma_to_Nucleon_Q1_Eye_site(const mobj &Dq_loop, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) +{ + + Gamma g5(Gamma::Algebra::Gamma5); + + auto DuG = Du_spec * GammaB_nucl; + // Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) + auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5; + // Dq_loop * \gamma_\mu^L + auto DqG = Dq_loop * Gamma_H; + + for (int ie_n=0; ie_n < 6 ; ie_n++){ + int a_n = epsilon[ie_n][0]; //a + int b_n = epsilon[ie_n][1]; //b + int c_n = epsilon[ie_n][2]; //c + for (int ie_s=0; ie_s < 6 ; ie_s++){ + int a_s = epsilon[ie_s][0]; //a' + int b_s = epsilon[ie_s][1]; //b' + int c_s = epsilon[ie_s][2]; //c' + for (int alpha_s=0; alpha_s +template +void BaryonUtils::Sigma_to_Nucleon_Q1_NonEye_site(const mobj &Du_ti, + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) +{ + + Gamma g5(Gamma::Algebra::Gamma5); + + auto DuG = Du_spec * GammaB_nucl; + auto adjDu = g5 * adj(Du_tf) * g5; + auto adjDuG = adjDu * GammaB_nucl; + // Gamma^B * Ds * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) + auto GDsGDd = GammaB_sigma * Ds_ti * Gamma_H * g5 * adj(Dd_tf) * g5; + // Dq_loop * \gamma_\mu^L + auto DuGH = Du_ti * Gamma_H; + + for (int ie_n=0; ie_n < 6 ; ie_n++){ + int a_n = epsilon[ie_n][0]; //a + int b_n = epsilon[ie_n][1]; //b + int c_n = epsilon[ie_n][2]; //c + for (int ie_s=0; ie_s < 6 ; ie_s++){ + int a_s = epsilon[ie_s][0]; //a' + int b_s = epsilon[ie_s][1]; //b' + int c_s = epsilon[ie_s][2]; //c' + for (int alpha_s=0; alpha_s +template +void BaryonUtils::Sigma_to_Nucleon_Q2_Eye_site(const mobj &Dq_loop, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) +{ + + Gamma g5(Gamma::Algebra::Gamma5); + + auto DuG = Du_spec * GammaB_nucl; + // Gamma^B * Ds * \gamma_\mu^L + auto GDsG = GammaB_sigma * Ds_ti * Gamma_H; + // Dq_loop * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) + auto DqGDd = Dq_loop * Gamma_H * g5 * adj(Dd_tf) * g5; + + for (int ie_n=0; ie_n < 6 ; ie_n++){ + int a_n = epsilon[ie_n][0]; //a + int b_n = epsilon[ie_n][1]; //b + int c_n = epsilon[ie_n][2]; //c + for (int ie_s=0; ie_s < 6 ; ie_s++){ + int a_s = epsilon[ie_s][0]; //a' + int b_s = epsilon[ie_s][1]; //b' + int c_s = epsilon[ie_s][2]; //c' + for (int alpha_s=0; alpha_s +template +void BaryonUtils::Sigma_to_Nucleon_Q2_NonEye_site(const mobj &Du_ti, + const mobj &Du_tf, + const mobj2 &Du_spec, + const mobj &Dd_tf, + const mobj &Ds_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + robj &result) +{ + + Gamma g5(Gamma::Algebra::Gamma5); + + auto DuG = Du_spec * GammaB_nucl; + auto adjDu = g5 * adj(Du_tf) * g5; + auto adjDuG = adjDu * GammaB_nucl; + // Gamma^B * Ds * \gamma_\mu^L + auto GDsG = GammaB_sigma * Ds_ti * Gamma_H; + // Du * \gamma_\mu^L * (\gamma_5 * Dd^\dagger * \gamma_5) + auto DuGDd = Du_ti * Gamma_H * g5 * adj(Dd_tf) * g5; + + for (int ie_n=0; ie_n < 6 ; ie_n++){ + int a_n = epsilon[ie_n][0]; //a + int b_n = epsilon[ie_n][1]; //b + int c_n = epsilon[ie_n][2]; //c + for (int ie_s=0; ie_s < 6 ; ie_s++){ + int a_s = epsilon[ie_s][0]; //a' + int b_s = epsilon[ie_s][1]; //b' + int c_s = epsilon[ie_s][2]; //c' + for (int alpha_s=0; alpha_s +template +void BaryonUtils::Sigma_to_Nucleon_Eye(const PropagatorField &qq_loop, + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr) +{ + GridBase *grid = qs_ti.Grid(); + + auto vcorr= stn_corr.View(); + auto vq_loop = qq_loop.View(); + auto vd_tf = qd_tf.View(); + auto vs_ti = qs_ti.View(); + + // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + thread_for(ss,grid->oSites(),{ + auto Dq_loop = vq_loop[ss]; + auto Dd_tf = vd_tf[ss]; + auto Ds_ti = vs_ti[ss]; + sobj result=Zero(); + if(op == "Q1"){ + Sigma_to_Nucleon_Q1_Eye_site(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + } else if(op == "Q2"){ + Sigma_to_Nucleon_Q2_Eye_site(Dq_loop,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + } else { + assert(0 && "Weak Operator not correctly specified"); + } + vcorr[ss] = result; + } );//end loop over lattice sites +} + +template +template +void BaryonUtils::Sigma_to_Nucleon_NonEye(const PropagatorField &qq_ti, + const PropagatorField &qq_tf, + const mobj &Du_spec, + const PropagatorField &qd_tf, + const PropagatorField &qs_ti, + const Gamma Gamma_H, + const Gamma GammaB_sigma, + const Gamma GammaB_nucl, + const std::string op, + SpinMatrixField &stn_corr) +{ + GridBase *grid = qs_ti.Grid(); + + auto vcorr= stn_corr.View(); + auto vq_ti = qq_ti.View(); + auto vq_tf = qq_tf.View(); + auto vd_tf = qd_tf.View(); + auto vs_ti = qs_ti.View(); + + // accelerator_for(ss, grid->oSites(), grid->Nsimd(), { + thread_for(ss,grid->oSites(),{ + auto Dq_ti = vq_ti[ss]; + auto Dq_tf = vq_tf[ss]; + auto Dd_tf = vd_tf[ss]; + auto Ds_ti = vs_ti[ss]; + sobj result=Zero(); + if(op == "Q1"){ + Sigma_to_Nucleon_Q1_NonEye_site(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + } else if(op == "Q2"){ + Sigma_to_Nucleon_Q2_NonEye_site(Dq_ti,Dq_tf,Du_spec,Dd_tf,Ds_ti,Gamma_H,GammaB_sigma,GammaB_nucl,result); + } else { + assert(0 && "Weak Operator not correctly specified"); + } + vcorr[ss] = result; + } );//end loop over lattice sites +} + NAMESPACE_END(Grid); diff --git a/Hadrons/Modules.hpp b/Hadrons/Modules.hpp index 90698f91..3c30f206 100644 --- a/Hadrons/Modules.hpp +++ b/Hadrons/Modules.hpp @@ -12,6 +12,8 @@ #include #include #include +#include +#include #include #include #include diff --git a/Hadrons/Modules/MContraction/SigmaToNucleonEye.cc b/Hadrons/Modules/MContraction/SigmaToNucleonEye.cc new file mode 100644 index 00000000..ca4bb341 --- /dev/null +++ b/Hadrons/Modules/MContraction/SigmaToNucleonEye.cc @@ -0,0 +1,7 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MContraction; + +template class Grid::Hadrons::MContraction::TSigmaToNucleonEye; diff --git a/Hadrons/Modules/MContraction/SigmaToNucleonEye.hpp b/Hadrons/Modules/MContraction/SigmaToNucleonEye.hpp new file mode 100644 index 00000000..b65d0689 --- /dev/null +++ b/Hadrons/Modules/MContraction/SigmaToNucleonEye.hpp @@ -0,0 +1,218 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MContraction/SigmaToNucleonEye.hpp + +Copyright (C) 2015-2019 + +Author: Antonin Portelli +Author: Felix Erben + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_MContraction_SigmaToNucleonEye_hpp_ +#define Hadrons_MContraction_SigmaToNucleonEye_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * SigmaToNucleonEye * + ******************************************************************************/ +/* + * Sigma-to-nucleon 3-pt diagrams, eye topologies. + * + * Schematics: qqLoop | + * /->-¬ | + * / \ | qsTi G qdTf + * \ / | /---->------*------>----¬ + * qsTi \ / qdTf | / /-*-¬ \ + * /----->-----* *----->----¬ | / / G \ \ + * * G G * | * \ / qqLoop * + * |\ /| | |\ \-<-/ /| + * | \ / | | | \ / | + * | \---------->---------/ | | | \----------->----------/ | + * \ quSpec / | \ quSpec / + * \ / | \ / + * \---------->---------/ | \----------->----------/ + * quSpec | quSpec + * + * analogously to the rare-kaon naming, the left diagram is named 'one-trace' and + * the diagram on the right 'two-trace' + * + * Propagators: + * * qqLoop + * * quSpec, source at ti + * * qdTf, source at tf + * * qsTi, source at ti + */ +BEGIN_MODULE_NAMESPACE(MContraction) + +class SigmaToNucleonEyePar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(SigmaToNucleonEyePar, + std::string, qqLoop, + std::string, quSpec, + std::string, qdTf, + std::string, qsTi, + unsigned int, tf, + std::string, sink, + std::string, output); +}; + +template +class TSigmaToNucleonEye: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + BASIC_TYPE_ALIASES(ScalarImplCR, Scalar); + SINK_TYPE_ALIASES(Scalar); + typedef typename SpinMatrixField::vector_object::scalar_object SpinMatrix; + class Metadata: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Metadata, + Gamma::Algebra, gammaH, + Gamma::Algebra, gammaASigma, + Gamma::Algebra, gammaBSigma, + Gamma::Algebra, gammaANucl, + Gamma::Algebra, gammaBNucl, + int, trace); + }; + typedef Correlator Result; +public: + // constructor + TSigmaToNucleonEye(const std::string name); + // destructor + virtual ~TSigmaToNucleonEye(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); +protected: + // setup + virtual void setup(void); + // execution + virtual void execute(void); + // Which gamma algebra was specified + Gamma::Algebra al; +}; + +MODULE_REGISTER_TMP(SigmaToNucleonEye, ARG(TSigmaToNucleonEye), MContraction); + +/****************************************************************************** + * TSigmaToNucleonEye implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TSigmaToNucleonEye::TSigmaToNucleonEye(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TSigmaToNucleonEye::getInput(void) +{ + std::vector input = {par().qqLoop, par().quSpec, par().qdTf, par().qsTi, par().sink}; + + return input; +} + +template +std::vector TSigmaToNucleonEye::getOutput(void) +{ + std::vector out = {}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TSigmaToNucleonEye::setup(void) +{ + envTmpLat(SpinMatrixField, "c"); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TSigmaToNucleonEye::execute(void) +{ + const Gamma GammaB(Gamma::Algebra::SigmaXZ); // C*gamma_5 + const Gamma Id(Gamma::Algebra::Identity); // C*gamma_5 + + LOG(Message) << "Computing sigma-to-nucleon contractions '" << getName() << "'" << std::endl; + LOG(Message) << "' with (Gamma^A,Gamma^B)_sigma = ( Identity, C*gamma_5 ) and (Gamma^A,Gamma^B)_nucl = ( Identity, C*gamma_5 )" << std::endl; + LOG(Message) << " using sink " << par().sink << "." << std::endl; + + envGetTmp(SpinMatrixField, c); + std::vector buf; + + std::vector result; + Result r; + r.info.gammaASigma = Id.g; + r.info.gammaBSigma = GammaB.g; + r.info.gammaANucl = Id.g; + r.info.gammaBNucl = GammaB.g; + + auto &qqLoop = envGet(PropagatorField, par().qqLoop); + auto &quSpec = envGet(SlicedPropagator, par().quSpec); + auto &qdTf = envGet(PropagatorField, par().qdTf); + auto &qsTi = envGet(PropagatorField, par().qsTi); + auto qut = quSpec[par().tf]; + for (auto &G: Gamma::gall) + { + r.info.gammaH = G.g; + //Operator Q1, equivalent to the two-trace case in the rare-kaons module + c=Zero(); + BaryonUtils::Sigma_to_Nucleon_Eye(qqLoop,qut,qdTf,qsTi,G,GammaB,GammaB,"Q1",c); + sliceSum(c,buf,Tp); + r.corr.clear(); + for (unsigned int t = 0; t < buf.size(); ++t) + { + r.corr.push_back(buf[t]); + } + r.info.trace = 2; + result.push_back(r); + //Operator Q2, equivalent to the one-trace case in the rare-kaons module + c=Zero(); + BaryonUtils::Sigma_to_Nucleon_Eye(qqLoop,qut,qdTf,qsTi,G,GammaB,GammaB,"Q2",c); + sliceSum(c,buf,Tp); + r.corr.clear(); + for (unsigned int t = 0; t < buf.size(); ++t) + { + r.corr.push_back(buf[t]); + } + r.info.trace = 1; + result.push_back(r); + } + + saveResult(par().output, "stnEye", result); + +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MContraction_SigmaToNucleonEye_hpp_ diff --git a/Hadrons/Modules/MContraction/SigmaToNucleonNonEye.cc b/Hadrons/Modules/MContraction/SigmaToNucleonNonEye.cc new file mode 100644 index 00000000..117541a3 --- /dev/null +++ b/Hadrons/Modules/MContraction/SigmaToNucleonNonEye.cc @@ -0,0 +1,7 @@ +#include + +using namespace Grid; +using namespace Hadrons; +using namespace MContraction; + +template class Grid::Hadrons::MContraction::TSigmaToNucleonNonEye; diff --git a/Hadrons/Modules/MContraction/SigmaToNucleonNonEye.hpp b/Hadrons/Modules/MContraction/SigmaToNucleonNonEye.hpp new file mode 100644 index 00000000..b0c86e2e --- /dev/null +++ b/Hadrons/Modules/MContraction/SigmaToNucleonNonEye.hpp @@ -0,0 +1,224 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Hadrons/Modules/MContraction/SigmaToNucleonNonEye.hpp + +Copyright (C) 2015-2019 + +Author: Antonin Portelli +Author: Felix Erben + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#ifndef Hadrons_MContraction_SigmaToNucleonNonEye_hpp_ +#define Hadrons_MContraction_SigmaToNucleonNonEye_hpp_ + +#include +#include +#include +#include + +BEGIN_HADRONS_NAMESPACE + +/****************************************************************************** + * SigmaToNucleonNonEye * + ******************************************************************************/ +/* + * Sigma-to-Nucleon 3-pt diagrams, non-eye topologies. + * + * Schematic: + * qsTi quTf | qsTi qdTf + * /-->--¬ /-->--¬ | /-->--¬ /-->--¬ + * / \ / \ | / \ / \ + * / \ / \ | / \ / \ + * / \ / \ | / \ / \ + * * * G * | * G * * G * + * |\ * G | | |\ / \ /| + * | \ / \ /| | | \ / \ / | + * | \ / \ / | | | \ / \ / | + * | \ / \ / | | | \-->--/ \-->--/ | + * \ \-->--/ \-->--/ / | \ quTi quTf / + * \ quTi qdTf / | \ / + * \ / | \ / + * \--------->----------/ | \--------->-----------/ + * quSpec | quSpec + * + * + * analogously to the rare-kaon naming, the left diagram is named 'one-trace' and + * the diagram on the right 'two-trace' + * + * Propagators: + * * quTi, source at ti + * * quTf, source at tf + * * quSpec, source at ti + * * qdTf, source at tf + * * qsTi, source at ti + */ +BEGIN_MODULE_NAMESPACE(MContraction) + +class SigmaToNucleonNonEyePar: Serializable +{ +public: + GRID_SERIALIZABLE_CLASS_MEMBERS(SigmaToNucleonNonEyePar, + std::string, quTi, + std::string, quTf, + std::string, quSpec, + std::string, qdTf, + std::string, qsTi, + unsigned int, tf, + std::string, sink, + std::string, output); +}; + +template +class TSigmaToNucleonNonEye: public Module +{ +public: + FERM_TYPE_ALIASES(FImpl,); + BASIC_TYPE_ALIASES(ScalarImplCR, Scalar); + SINK_TYPE_ALIASES(Scalar); + typedef typename SpinMatrixField::vector_object::scalar_object SpinMatrix; + class Metadata: Serializable + { + public: + GRID_SERIALIZABLE_CLASS_MEMBERS(Metadata, + Gamma::Algebra, gammaH, + Gamma::Algebra, gammaASigma, + Gamma::Algebra, gammaBSigma, + Gamma::Algebra, gammaANucl, + Gamma::Algebra, gammaBNucl, + int, trace); + }; + typedef Correlator Result; +public: + // constructor + TSigmaToNucleonNonEye(const std::string name); + // destructor + virtual ~TSigmaToNucleonNonEye(void) {}; + // dependency relation + virtual std::vector getInput(void); + virtual std::vector getOutput(void); +protected: + // setup + virtual void setup(void); + // execution + virtual void execute(void); + // Which gamma algebra was specified + Gamma::Algebra al; +}; + +MODULE_REGISTER_TMP(SigmaToNucleonNonEye, ARG(TSigmaToNucleonNonEye), MContraction); + +/****************************************************************************** + * TSigmaToNucleonNonEye implementation * + ******************************************************************************/ +// constructor ///////////////////////////////////////////////////////////////// +template +TSigmaToNucleonNonEye::TSigmaToNucleonNonEye(const std::string name) +: Module(name) +{} + +// dependencies/products /////////////////////////////////////////////////////// +template +std::vector TSigmaToNucleonNonEye::getInput(void) +{ + std::vector input = {par().quTi, par().quTf, par().quSpec, par().qdTf, par().qsTi, par().sink}; + + return input; +} + +template +std::vector TSigmaToNucleonNonEye::getOutput(void) +{ + std::vector out = {}; + + return out; +} + +// setup /////////////////////////////////////////////////////////////////////// +template +void TSigmaToNucleonNonEye::setup(void) +{ + envTmpLat(SpinMatrixField, "c"); +} + +// execution /////////////////////////////////////////////////////////////////// +template +void TSigmaToNucleonNonEye::execute(void) +{ + const Gamma GammaB(Gamma::Algebra::SigmaXZ); // C*gamma_5 + const Gamma Id(Gamma::Algebra::Identity); // C*gamma_5 + + LOG(Message) << "Computing sigma-to-nucleon contractions '" << getName() << "'" << std::endl; + LOG(Message) << "' with (Gamma^A,Gamma^B)_sigma = ( Identity, C*gamma_5 ) and (Gamma^A,Gamma^B)_nucl = ( Identity, C*gamma_5 )" << std::endl; + LOG(Message) << " using sink " << par().sink << "." << std::endl; + + envGetTmp(SpinMatrixField, c); + std::vector buf; + + std::vector result; + Result r; + r.info.gammaASigma = Id.g; + r.info.gammaBSigma = GammaB.g; + r.info.gammaANucl = Id.g; + r.info.gammaBNucl = GammaB.g; + + auto &quTi = envGet(PropagatorField, par().quTi); + auto &quTf = envGet(PropagatorField, par().quTf); + auto &quSpec = envGet(SlicedPropagator, par().quSpec); + auto &qdTf = envGet(PropagatorField, par().qdTf); + auto &qsTi = envGet(PropagatorField, par().qsTi); + auto qut = quSpec[par().tf]; + for (auto &G: Gamma::gall) + { + r.info.gammaH = G.g; + //Operator Q1, equivalent to the two-trace case in the rare-kaons module + c=Zero(); + BaryonUtils::Sigma_to_Nucleon_NonEye(quTi,quTf,qut,qdTf,qsTi,G,GammaB,GammaB,"Q1",c); + sliceSum(c,buf,Tp); + r.corr.clear(); + for (unsigned int t = 0; t < buf.size(); ++t) + { + r.corr.push_back(buf[t]); + } + r.info.trace = 2; + result.push_back(r); + //Operator Q2, equivalent to the one-trace case in the rare-kaons module + c=Zero(); + BaryonUtils::Sigma_to_Nucleon_NonEye(quTi,quTf,qut,qdTf,qsTi,G,GammaB,GammaB,"Q2",c); + sliceSum(c,buf,Tp); + r.corr.clear(); + for (unsigned int t = 0; t < buf.size(); ++t) + { + r.corr.push_back(buf[t]); + } + r.info.trace = 1; + result.push_back(r); + } + + saveResult(par().output, "stnNonEye", result); + +} + +END_MODULE_NAMESPACE + +END_HADRONS_NAMESPACE + +#endif // Hadrons_MContraction_SigmaToNucleonNonEye_hpp_ diff --git a/Hadrons/modules.inc b/Hadrons/modules.inc index 4fd3600e..2eb1923d 100644 --- a/Hadrons/modules.inc +++ b/Hadrons/modules.inc @@ -13,6 +13,8 @@ modules_cc =\ Modules/MContraction/DiscLoop.cc \ Modules/MContraction/Gamma3pt.cc \ Modules/MContraction/Meson.cc \ + Modules/MContraction/SigmaToNucleonEye.cc \ + Modules/MContraction/SigmaToNucleonNonEye.cc \ Modules/MContraction/WeakEye3pt.cc \ Modules/MContraction/WeakMesonDecayKl2.cc \ Modules/MContraction/WeakNonEye3pt.cc \ @@ -95,6 +97,8 @@ modules_hpp =\ Modules/MContraction/DiscLoop.hpp \ Modules/MContraction/Gamma3pt.hpp \ Modules/MContraction/Meson.hpp \ + Modules/MContraction/SigmaToNucleonEye.hpp \ + Modules/MContraction/SigmaToNucleonNonEye.hpp \ Modules/MContraction/WeakEye3pt.hpp \ Modules/MContraction/WeakMesonDecayKl2.hpp \ Modules/MContraction/WeakNonEye3pt.hpp \ diff --git a/tests/hadrons/Test_sigma_to_nucleon.cc b/tests/hadrons/Test_sigma_to_nucleon.cc new file mode 100644 index 00000000..b4e79649 --- /dev/null +++ b/tests/hadrons/Test_sigma_to_nucleon.cc @@ -0,0 +1,154 @@ +/************************************************************************************* + +Grid physics library, www.github.com/paboyle/Grid + +Source file: Tests/Hadrons/Test_hadrons_spectrum.cc + +Copyright (C) 2015-2018 + + Author: Antonin Portelli + +This program is free software; you can redistribute it and/or modify +it under the terms of the GNU General Public License as published by +the Free Software Foundation; either version 2 of the License, or +(at your option) any later version. + +This program is distributed in the hope that it will be useful, +but WITHOUT ANY WARRANTY; without even the implied warranty of +MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +GNU General Public License for more details. + +You should have received a copy of the GNU General Public License along +with this program; if not, write to the Free Software Foundation, Inc., +51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + +See the full license in the file "LICENSE" in the top level distribution directory +*************************************************************************************/ +/* END LEGAL */ + +#include +#include + +using namespace Grid; +using namespace Hadrons; + +int main(int argc, char *argv[]) +{ + // initialization ////////////////////////////////////////////////////////// + Grid_init(&argc, &argv); + HadronsLogError.Active(GridLogError.isActive()); + HadronsLogWarning.Active(GridLogWarning.isActive()); + HadronsLogMessage.Active(GridLogMessage.isActive()); + HadronsLogIterative.Active(GridLogIterative.isActive()); + HadronsLogDebug.Active(GridLogDebug.isActive()); + LOG(Message) << "Grid initialized" << std::endl; + + // run setup /////////////////////////////////////////////////////////////// + Application application; + std::vector flavour = {"l", "s", "c"}; + std::vector mass = {.01, .04, .2 }; + + // global parameters + Application::GlobalPar globalPar; + globalPar.trajCounter.start = 1500; + globalPar.trajCounter.end = 1520; + globalPar.trajCounter.step = 20; + globalPar.runId = "test"; + application.setPar(globalPar); + // gauge field + application.createModule("gauge"); + // sources + MSource::Point::Par ptPar; + ptPar.position = "0 0 0 0"; + application.createModule("pt_0", ptPar); + ptPar.position = "0 0 0 4"; + application.createModule("pt_4", ptPar); + // sink + MSink::Point::Par sinkPar; + sinkPar.mom = "0 0 0"; + application.createModule("sink", sinkPar); + application.createModule("sink_spec", sinkPar); + + // set fermion boundary conditions to be periodic space, antiperiodic time. + std::string boundary = "1 1 1 -1"; + std::string twist = "0. 0. 0. 0."; + + for (unsigned int i = 0; i < flavour.size(); ++i) + { + // actions + MAction::DWF::Par actionPar; + actionPar.gauge = "gauge"; + actionPar.Ls = 12; + actionPar.M5 = 1.8; + actionPar.mass = mass[i]; + actionPar.boundary = boundary; + actionPar.twist = twist; + application.createModule("DWF_" + flavour[i], actionPar); + + // solvers + MSolver::RBPrecCG::Par solverPar; + solverPar.action = "DWF_" + flavour[i]; + solverPar.residual = 1.0e-8; + solverPar.maxIteration = 10000; + application.createModule("CG_" + flavour[i], + solverPar); + + } + + // propagators + MFermion::GaugeProp::Par quarkPar; + quarkPar.solver = "CG_l"; + quarkPar.source = "pt_0"; + application.createModule("Qpt_l_0", quarkPar); + quarkPar.source = "pt_4"; + application.createModule("Qpt_l_4", quarkPar); + quarkPar.solver = "CG_s"; + quarkPar.source = "pt_0"; + application.createModule("Qpt_s_0", quarkPar); + //This should be a loop - how do I make this? + quarkPar.solver = "CG_c"; + quarkPar.source = "pt_0"; + application.createModule("Qpt_c_loop", quarkPar); + quarkPar.solver = "CG_l"; + quarkPar.source = "pt_0"; + application.createModule("Qpt_l_loop", quarkPar); + + MSink::Smear::Par smearPar; + smearPar.q="Qpt_l_0"; + smearPar.sink = "sink_spec"; + application.createModule("Qpt_u_spec",smearPar); + + + MContraction::SigmaToNucleonEye::Par EyePar; + EyePar.output = "SigmaToNucleon/Eye_u"; + EyePar.qqLoop = "Qpt_l_loop"; + EyePar.quSpec = "Qpt_u_spec"; + EyePar.qdTf = "Qpt_l_4"; + EyePar.qsTi = "Qpt_s_0"; + EyePar.tf = 4; + EyePar.sink = "sink"; + application.createModule("SigmaToNucleonEye_u", EyePar); + EyePar.output = "SigmaToNucleon/Eye_c"; + EyePar.qqLoop = "Qpt_c_loop"; + application.createModule("SigmaToNucleonEye_c", EyePar); + MContraction::SigmaToNucleonNonEye::Par NonEyePar; + NonEyePar.output = "SigmaToNucleon/NonEye"; + NonEyePar.quTi = "Qpt_l_0"; + NonEyePar.quTf = "Qpt_l_4"; + NonEyePar.quSpec = "Qpt_u_spec"; + NonEyePar.qdTf = "Qpt_l_4"; + NonEyePar.qsTi = "Qpt_s_0"; + NonEyePar.tf = 4; + NonEyePar.sink = "sink"; + application.createModule("SigmaToNucleonNonEye", NonEyePar); + + // execution + application.saveParameterFile("stn.xml"); + application.run(); + + // epilogue + LOG(Message) << "Grid is finalizing now" << std::endl; + Grid_finalize(); + + return EXIT_SUCCESS; +}