/*************************************************************************************

Grid physics library, www.github.com/paboyle/Grid 

Source file: extras/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.cc

Copyright (C) 2015-2018

Author: Antonin Portelli <antonin.portelli@me.com>
Author: Lanny91 <andrew.lawson@gmail.com>

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 <Grid/Hadrons/Modules/MContraction/WeakNeutral4ptDisc.hpp>

using namespace Grid;
using namespace Hadrons;
using namespace MContraction;

/*
 * Weak Hamiltonian + current contractions, disconnected topology for neutral 
 * mesons.
 * 
 * These contractions are generated by operators Q_1,...,10 of the dS=1 Weak
 * Hamiltonian in the physical basis and an additional current J (see e.g. 
 * Fig 11 of arXiv:1507.03094).
 * 
 * Schematic:
 *                        
 *           q2          q4             q3
 *       /--<--¬     /---<--¬       /---<--¬
 *     /         \ /         \     /        \
 *  i *           * H_W      |  J *          * f
 *     \         / \         /     \        /
 *      \--->---/   \-------/       \------/
 *          q1 
 * 
 * options
 * - q1: input propagator 1 (string)
 * - q2: input propagator 2 (string)
 * - q3: input propagator 3 (string), assumed to be sequential propagator 
 * - q4: input propagator 4 (string), assumed to be a loop
 * 
 * type 1: trace(q1*adj(q2)*g5*gL[mu])*trace(loop*gL[mu])*trace(q3*g5)
 * type 2: trace(q1*adj(q2)*g5*gL[mu]*loop*gL[mu])*trace(q3*g5)
 */

/*******************************************************************************
 *                  TWeakNeutral4ptDisc implementation                         *
 ******************************************************************************/
// constructor /////////////////////////////////////////////////////////////////
TWeakNeutral4ptDisc::TWeakNeutral4ptDisc(const std::string name)
: Module<WeakHamiltonianPar>(name)
{}

// dependencies/products ///////////////////////////////////////////////////////
std::vector<std::string> TWeakNeutral4ptDisc::getInput(void)
{
    std::vector<std::string> in = {par().q1, par().q2, par().q3, par().q4};
    
    return in;
}

std::vector<std::string> TWeakNeutral4ptDisc::getOutput(void)
{
    std::vector<std::string> out = {};
    
    return out;
}

// setup ///////////////////////////////////////////////////////////////////////
void TWeakNeutral4ptDisc::setup(void)
{
    unsigned int ndim = env().getNd();

    envTmpLat(LatticeComplex,  "expbuf");
    envTmpLat(PropagatorField, "tmp");
    envTmpLat(LatticeComplex,  "curr");
    envTmp(std::vector<PropagatorField>, "meson", 1, ndim, PropagatorField(env().getGrid()));
    envTmp(std::vector<PropagatorField>, "loop", 1, ndim,  PropagatorField(env().getGrid()));
}

// execution ///////////////////////////////////////////////////////////////////
void TWeakNeutral4ptDisc::execute(void)
{
    LOG(Message) << "Computing Weak Hamiltonian neutral disconnected contractions '" 
                 << getName() << "' using quarks '" << par().q1 << "', '" 
                 << par().q2 << ", '" << par().q3 << "' and '" << par().q4 
                 << "'." << std::endl;

    auto                  &q1 = envGet(PropagatorField, par().q1);
    auto                  &q2 = envGet(PropagatorField, par().q2);
    auto                  &q3 = envGet(PropagatorField, par().q3);
    auto                  &q4 = envGet(PropagatorField, par().q4);
    Gamma                 g5  = Gamma(Gamma::Algebra::Gamma5);
    std::vector<TComplex> corrbuf;
    std::vector<Result>   result(n_neut_disc_diag);
    unsigned int          ndim = env().getNd();

    envGetTmp(LatticeComplex,               expbuf); 
    envGetTmp(PropagatorField,              tmp);
    envGetTmp(LatticeComplex,               curr);
    envGetTmp(std::vector<PropagatorField>, meson);
    envGetTmp(std::vector<PropagatorField>, loop);

    // Setup for type 1 contractions.
    for (int mu = 0; mu < ndim; ++mu)
    {
        meson[mu] = MAKE_DISC_MESON(q1, q2, GammaL(Gamma::gmu[mu]));
        loop[mu] = MAKE_DISC_LOOP(q4, GammaL(Gamma::gmu[mu]));
    }
    curr = MAKE_DISC_CURR(q3, GammaL(Gamma::Algebra::Gamma5));

    // Perform type 1 contractions.    
    SUM_MU(expbuf, trace(meson[mu]*loop[mu]))
    expbuf *= curr;
    MAKE_DIAG(expbuf, corrbuf, result[neut_disc_1_diag], "HW_disc0_1")

    // Perform type 2 contractions.
    SUM_MU(expbuf, trace(meson[mu])*trace(loop[mu]))
    expbuf *= curr;
    MAKE_DIAG(expbuf, corrbuf, result[neut_disc_2_diag], "HW_disc0_2")

    // IO
    saveResult(par().output, "HW_disc0", result);
}