From 581be32ed24b9a39aaf22fd05cda185f13cb884e Mon Sep 17 00:00:00 2001
From: James Harrison <jch1g10@soton.ac.uk>
Date: Thu, 14 Dec 2017 13:42:41 +0000
Subject: [PATCH] Implement infrared improvement for v=0 on-shell self-energy

---
 extras/Hadrons/Modules/MGauge/StochEm.cc  |  2 +-
 extras/Hadrons/Modules/MGauge/StochEm.hpp |  3 +-
 extras/Hadrons/Modules/MGauge/UnitEm.cc   |  2 +-
 lib/qcd/action/gauge/Photon.h             | 34 ++++++++++++++++++++---
 4 files changed, 34 insertions(+), 7 deletions(-)

diff --git a/extras/Hadrons/Modules/MGauge/StochEm.cc b/extras/Hadrons/Modules/MGauge/StochEm.cc
index c7a9fc4f..8f84fe94 100644
--- a/extras/Hadrons/Modules/MGauge/StochEm.cc
+++ b/extras/Hadrons/Modules/MGauge/StochEm.cc
@@ -67,7 +67,7 @@ void TStochEm::setup(void)
 // execution ///////////////////////////////////////////////////////////////////
 void TStochEm::execute(void)
 {
-    PhotonR photon(par().gauge, par().zmScheme);
+    PhotonR photon(par().gauge, par().zmScheme, par().improvement);
     EmField &a = *env().createLattice<EmField>(getName());
     EmComp  *w;
     
diff --git a/extras/Hadrons/Modules/MGauge/StochEm.hpp b/extras/Hadrons/Modules/MGauge/StochEm.hpp
index 12ce9fdc..b6bfa262 100644
--- a/extras/Hadrons/Modules/MGauge/StochEm.hpp
+++ b/extras/Hadrons/Modules/MGauge/StochEm.hpp
@@ -44,7 +44,8 @@ class StochEmPar: Serializable
 public:
     GRID_SERIALIZABLE_CLASS_MEMBERS(StochEmPar,
                                     PhotonR::Gauge,    gauge,
-                                    PhotonR::ZmScheme, zmScheme);
+                                    PhotonR::ZmScheme, zmScheme,
+                                    Integer, improvement);
 };
 
 class TStochEm: public Module<StochEmPar>
diff --git a/extras/Hadrons/Modules/MGauge/UnitEm.cc b/extras/Hadrons/Modules/MGauge/UnitEm.cc
index f33dfed3..747ef142 100644
--- a/extras/Hadrons/Modules/MGauge/UnitEm.cc
+++ b/extras/Hadrons/Modules/MGauge/UnitEm.cc
@@ -61,7 +61,7 @@ void TUnitEm::setup(void)
 // execution ///////////////////////////////////////////////////////////////////
 void TUnitEm::execute(void)
 {
-    PhotonR photon(0, 0); // Just chose arbitrary input values here
+    PhotonR photon(0, 0, 0); // Just chose arbitrary input values here
     EmField &a = *env().createLattice<EmField>(getName());
     LOG(Message) << "Generating unit EM potential..." << std::endl;
     photon.UnitField(a);
diff --git a/lib/qcd/action/gauge/Photon.h b/lib/qcd/action/gauge/Photon.h
index 1429c2ba..e0329bed 100644
--- a/lib/qcd/action/gauge/Photon.h
+++ b/lib/qcd/action/gauge/Photon.h
@@ -60,7 +60,7 @@ namespace QCD{
     GRID_SERIALIZABLE_ENUM(Gauge, undef, feynman, 1, coulomb, 2, landau, 3);
     GRID_SERIALIZABLE_ENUM(ZmScheme, undef, qedL, 1, qedTL, 2);
   public:
-    Photon(Gauge gauge, ZmScheme zmScheme);
+    Photon(Gauge gauge, ZmScheme zmScheme, Integer improvement);
     virtual ~Photon(void) = default;
     void FreePropagator(const GaugeField &in, GaugeField &out);
     void MomentumSpacePropagator(const GaugeField &in, GaugeField &out);
@@ -75,13 +75,14 @@ namespace QCD{
   private:
     Gauge    gauge_;
     ZmScheme zmScheme_;
+    Integer  improvement_;
   };
 
   typedef Photon<QedGimplR>  PhotonR;
   
   template<class Gimpl>
-  Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme)
-  : gauge_(gauge), zmScheme_(zmScheme)
+  Photon<Gimpl>::Photon(Gauge gauge, ZmScheme zmScheme, Integer improvement)
+  : gauge_(gauge), zmScheme_(zmScheme), improvement_(improvement)
   {}
   
   template<class Gimpl>
@@ -128,6 +129,7 @@ namespace QCD{
   {
     GridBase           *grid = out._grid;
     const unsigned int nd    = grid->_ndimension;
+    std::vector<int>   &l    = grid->_fdimensions;
     
     switch (zmScheme_)
     {
@@ -149,9 +151,33 @@ namespace QCD{
         for(int d = 0; d < grid->_ndimension - 1; d++)
         {
           LatticeCoordinate(coor,d);
+          coor = where(coor < Integer(l[d]/2), coor, coor-Integer(l[d]));
           spNrm = spNrm + coor*coor;
         }
         out = where(spNrm == Integer(0), 0.*out, out);
+
+        // IR improvement
+        switch (improvement_)
+        {
+          case 0:
+            break;
+          case 1:
+          {
+            Real f1 = sqrt(2.48560548);
+            out = where(spNrm == Integer(1), f1*out, out);
+            break;
+          }
+          case 2:
+          {
+            Real f1 = sqrt(4.93053406);
+            Real f2 = sqrt(-1.44492857);
+            out = where(spNrm == Integer(1), f1*out, out);
+            out = where(spNrm == Integer(2), f2*out, out);
+            break;
+          }
+          default:
+            break;
+        }
         
         break;
       }
@@ -159,7 +185,7 @@ namespace QCD{
         break;
     }
   }
-  
+
   template<class Gimpl>
   void Photon<Gimpl>::MomentumSpacePropagator(const GaugeField &in,
                                                GaugeField &out)