diff --git a/projects/generate_confs/parser.ypp b/projects/generate_confs/parser.ypp
index 0b547a382..42f17f1ba 100644
--- a/projects/generate_confs/parser.ypp
+++ b/projects/generate_confs/parser.ypp
@@ -237,6 +237,10 @@
 %token TK_GAUSS_NITER_SNK
 %type <double_numb> gauss_kappa
 %type <int_numb> gauss_niter_src gauss_niter_snk
+%token TK_APE_ALPHA
+%token TK_APE_NITER
+%type <double_numb> ape_alpha
+%type <int_numb> ape_niter
 //putpourri
 %type <fermionic_putpourri_meas> fermionic_putpourri_meas
 %token TK_MEAS_PUTPOURRI
@@ -609,6 +613,8 @@ meson_corr_meas: TK_MEAS_MESON_CORRS {$$=new meson_corr_meas_pars_t();}
                | meson_corr_meas gauss_kappa {$$->gauss_kappa=$2;}
                | meson_corr_meas gauss_niter_src {$$->gauss_niter_src=$2;}
                | meson_corr_meas gauss_niter_snk {$$->gauss_niter_snk=$2;}
+               | meson_corr_meas ape_alpha {$$->ape_alpha=$2;}
+               | meson_corr_meas ape_niter {$$->ape_niter=$2;}
 ;
 
 dir: TK_DIR '=' int_numb {$$=$3;}
@@ -623,6 +629,12 @@ gauss_niter_src: TK_GAUSS_NITER_SRC '=' int_numb {$$=$3;}
 gauss_niter_snk: TK_GAUSS_NITER_SNK '=' int_numb {$$=$3;}
 ;
 
+ape_alpha: TK_APE_ALPHA '=' double_numb {$$=$3;}
+;
+
+ape_niter: TK_APE_NITER '=' int_numb {$$=$3;}
+;
+
 operators: TK_OPERATORS '=' int_pair_list {
     for(size_t i=0;i<$3->size();i++)
 	{
diff --git a/projects/generate_confs/tokenizer.lpp b/projects/generate_confs/tokenizer.lpp
index 0560964bb..4274696a7 100644
--- a/projects/generate_confs/tokenizer.lpp
+++ b/projects/generate_confs/tokenizer.lpp
@@ -134,6 +134,10 @@ GaussKappa DEBUG_PRINTF("Found GaussKappa\n");return TK_GAUSS_KAPPA;
 GaussNiterSrc DEBUG_PRINTF("Found GaussNiterSrc\n");return TK_GAUSS_NITER_SRC;
 GaussNiterSnk DEBUG_PRINTF("Found GaussNiterSnk\n");return TK_GAUSS_NITER_SNK;
 
+ /* ape smearing (meson corrs) */
+ApeAlpha DEBUG_PRINTF("Found ApeAlpha\n");return TK_APE_ALPHA;
+ApeNiter DEBUG_PRINTF("Found ApeNiter\n");return TK_APE_NITER;
+
  /* rendens pars */
 MaxOrder DEBUG_PRINTF("Found MaxOrder\n");return TK_MAX_ORDER;
 
diff --git a/src/measures/fermions/mesons.cpp b/src/measures/fermions/mesons.cpp
index e97b2808a..4a45b1e1d 100644
--- a/src/measures/fermions/mesons.cpp
+++ b/src/measures/fermions/mesons.cpp
@@ -16,6 +16,7 @@
 #include "routines/mpi_routines.hpp"
 #include "stag.hpp"
 #include "operations/smearing/gaussian.hpp"
+#include "operations/smearing/APE.hpp"
 
 namespace nissa
 {
@@ -26,9 +27,9 @@ namespace nissa
     int nop;
     int ncombo;
     int nflavs;
-    
-    /// Compute the index where to store
-    inline int icombo(const int& iflav_so,
+  
+  /// Compute the index where to store
+  inline int icombo(const int& iflav_so,
 		      const int& iflav_si,
 		      const int& iop_so,
 		      const int& iop_si,
@@ -43,26 +44,30 @@ namespace nissa
     }
   }
   
-  // Reuse gaussian_smearing using lx as buffer
+
+
+  // reuse gaussian_smearing using lx as buffer
   template <typename T>
-  void gaussian_smear_eo(EoField<T>& out,
-			 const EoField<T>& in,
-			 LxField<T>& lx_buf,
-			 LxField<T>& lx_temp,
-			 LxField<T>& lx_H,
-			 const LxField<quad_su3>& conf_lx,
-			 const double gauss_kappa,
-			 const int gauss_niter)
-  {
-    if(gauss_niter<=0)
-      out=in;
-    else
-      {    
-	paste_eo_parts_into_lx_vector(lx_buf,in);
-	gaussian_smearing(lx_buf,lx_buf,conf_lx,gauss_kappa,gauss_niter,&lx_temp,&lx_H);
-	split_lx_vector_into_eo_parts(out,lx_buf);
-      }
-  }
+	void gaussian_smear_eo(EoField<T>& out,
+                       const EoField<T>& in,
+                       LxField<T>& lx_buf,
+                       LxField<T>& lx_temp,
+                       LxField<T>& lx_H,
+                       const LxField<quad_su3>& conf_lx,
+                       const double gauss_kappa,
+                       const int gauss_niter,
+					   const int corr_dir)
+    {
+      if(gauss_niter<=0)
+		out = in; 
+	  else
+		{
+	  paste_eo_parts_into_lx_vector(lx_buf,in);
+      gaussian_smearing(lx_buf,lx_buf,conf_lx,gauss_kappa,gauss_niter,corr_dir,&lx_temp,&lx_H);
+	  split_lx_vector_into_eo_parts(out,lx_buf);
+    	}
+	}
+  
   
   void compute_meson_corr(complex* corr,
 			  const EoField<quad_su3>& conf,
@@ -71,9 +76,12 @@ namespace nissa
   {
     const int& dir=
       meas_pars.dir;
-    const double gauss_kappa=meas_pars.gauss_kappa;
-    const int gauss_niter_src=meas_pars.gauss_niter_src;
-    const int gauss_niter_snk=meas_pars.gauss_niter_snk;
+	const double gauss_kappa=meas_pars.gauss_kappa;
+	const int gauss_niter_src=meas_pars.gauss_niter_src;
+	const int gauss_niter_snk=meas_pars.gauss_niter_snk;
+	const double ape_alpha=meas_pars.ape_alpha;
+	const int ape_niter=meas_pars.ape_niter;
+
     
     //allocate
     EoField<color> ori_source("ori_source",WITH_HALO);
@@ -83,21 +91,21 @@ namespace nissa
     std::vector<EoField<color>> quark(nflavs*nop,{"quark",WITH_HALO});
     std::vector<EoField<color>> quark0s(nflavs*nop,{"quark0s",WITH_HALO});
     
-    LxField<quad_su3> conf_lx("conf_lx",WITH_HALO);
-    paste_eo_parts_into_lx_vector(conf_lx,conf);
-    
+	LxField<quad_su3> conf_lx("conf_lx",WITH_HALO_EDGES);
+	LxField<quad_su3> conf_lx_ape("conf_lx_ape",WITH_HALO_EDGES);
+	paste_eo_parts_into_lx_vector(conf_lx,conf);
+	ape_orth_dir_smear_conf(conf_lx_ape, conf_lx, ape_alpha, ape_niter, dir);
+
     LxField<color> gauss_lx_buf("gauss_lx_buf",WITH_HALO);
     LxField<color> gauss_lx_temp("gauss_lx_temp",WITH_HALO);
     LxField<color> gauss_lx_H("gauss_lx_H",WITH_HALO);
-    
-    auto gSmear=
-      [&](EoField<color>& out,
-	  const int& nIter,
-	  const EoField<color>& in)
+
+	auto gSmear =
+      [&](EoField<color>& out, const int& nIter, const EoField<color>& in)
       {
-	gaussian_smear_eo(out,in,gauss_lx_buf,gauss_lx_temp,gauss_lx_H,conf_lx,gauss_kappa,nIter);
+        gaussian_smear_eo(out, in, gauss_lx_buf, gauss_lx_temp, gauss_lx_H, conf_lx_ape, gauss_kappa, nIter, dir);
       };
-    
+
     //form the masks
     int mask[nop],shift[nop];
     for(int iop=0;iop<nop;iop++)
@@ -128,24 +136,23 @@ namespace nissa
 	generate_fully_undiluted_eo_source(ori_source,meas_pars.rnd_type,source_coord,dir);
 	
 	EoField<color> smeared_ori_source("smeared_ori_source",WITH_HALO);
-	gSmear(smeared_ori_source,gauss_niter_src,ori_source);
-	
+    gSmear(smeared_ori_source, gauss_niter_src, ori_source);
+
 	for(int iflav=0;iflav<nflavs;iflav++)
 	  {
 	    for(int iop=0;iop<nop;iop++)
 	      {
-		const int idx=iflav*nop+iop;
-		apply_shift_op(source,temp[0],temp[1],conf,tp.backfield[iflav],shift[iop],smeared_ori_source);
-		put_stag_phases(source,mask[iop]);
-		
-		mult_Minv(quark[idx],conf,tp,iflav,meas_pars.residue,source);
-		
-		//smear also the sink prop before contraction
-		gSmear(quark[idx],gauss_niter_snk,quark[idx]);
+			const int idx=iflav*nop+iop;
+			apply_shift_op(source,temp[0],temp[1],conf,tp.backfield[iflav],shift[iop],smeared_ori_source);
+			put_stag_phases(source,mask[iop]);
+
+			mult_Minv(quark[idx],conf,tp,iflav,meas_pars.residue,source);
+
+			//smear also the sink prop before contraction
+			gSmear(quark[idx], gauss_niter_snk, quark[idx]);
 	      }
 	  }
-	
-	/// Sink
+	    /// Sink
 	for(int iflav=0;iflav<nflavs;iflav++)
 	  for(int iop=0;iop<nop;iop++)
 	    {
diff --git a/src/measures/fermions/mesons.hpp b/src/measures/fermions/mesons.hpp
index ae6d28799..c4b884ac1 100644
--- a/src/measures/fermions/mesons.hpp
+++ b/src/measures/fermions/mesons.hpp
@@ -12,6 +12,10 @@ namespace nissa
     std::vector<std::pair<int,int>> mesons;
     
     int dir;
+
+    double ape_alpha;
+
+    int ape_niter;
     
     double gauss_kappa;
     
@@ -24,7 +28,18 @@ namespace nissa
     {
       return 0;
     }
-    
+
+    //ape smearing default disabled
+    double def_ape_alpha() const 
+    {
+       return 0.5;
+    }
+
+    int def_ape_niter() const 
+    { 
+      return 0;
+    }
+
     // Smearing defaults: disabled
     double def_gauss_kappa() const
     {
@@ -72,15 +87,22 @@ namespace nissa
 	    }
 	  os<<"}\n";
 	}
+
+    if(ape_alpha!=def_ape_alpha() or full)
+  os<<" ApeAlpha\t=\t"<<ape_alpha<<"\n";
+
+    if(ape_niter!=def_ape_niter() or full)
+  os<<" ApeNiter\t=\t"<<ape_niter<<"\n";
+
+    if(gauss_kappa!=def_gauss_kappa() or full)
+  os<<" GaussKappa\t=\t"<<gauss_kappa<<"\n";
+
+    if(gauss_niter_src!=def_gauss_niter_src() or full)
+  os<<" GaussNiterSrc\t=\t"<<gauss_niter_src<<"\n";
+
+    if(gauss_niter_snk!=def_gauss_niter_snk() or full)
+  os<<" GaussNiterSnk\t=\t"<<gauss_niter_snk<<"\n";
       
-      if(gauss_kappa!=def_gauss_kappa() or full)
-	os<<" GaussKappa\t=\t"<<gauss_kappa<<"\n";
-      
-      if(gauss_niter_src!=def_gauss_niter_src() or full)
-	os<<" GaussNiterSrc\t=\t"<<gauss_niter_src<<"\n";
-      
-      if(gauss_niter_snk!=def_gauss_niter_snk() or full)
-	os<<" GaussNiterSnk\t=\t"<<gauss_niter_snk<<"\n";
       
       return os.str();
     }
@@ -91,15 +113,19 @@ namespace nissa
 	base_fermionic_meas_t::is_nonstandard() or
 	mesons.size() or
 	dir!=def_dir() or
-	gauss_kappa!=def_gauss_kappa() or
-	gauss_niter_src!=def_gauss_niter_src() or
-	gauss_niter_snk!=def_gauss_niter_snk() or
+  gauss_kappa!=def_gauss_kappa() or
+  gauss_niter_src!=def_gauss_niter_src() or
+  gauss_niter_snk!=def_gauss_niter_snk() or
+  ape_alpha!=def_ape_alpha() or
+  ape_niter!=def_ape_niter() or
 	path!=def_path();
     }
     
     meson_corr_meas_pars_t() :
       base_fermionic_meas_t(),
       dir(def_dir()),
+      ape_alpha(def_ape_alpha()),
+      ape_niter(def_ape_niter()),
       gauss_kappa(def_gauss_kappa()),
       gauss_niter_src(def_gauss_niter_src()),
       gauss_niter_snk(def_gauss_niter_snk())
diff --git a/src/operations/smearing/APE.cpp b/src/operations/smearing/APE.cpp
index 9504d9bc6..8c252a55f 100644
--- a/src/operations/smearing/APE.cpp
+++ b/src/operations/smearing/APE.cpp
@@ -11,7 +11,8 @@ namespace nissa
 		      const double& alpha,
 		      const int& nSteps,
 		      const WhichDirs& dirs,
-		      const int& minStapleDir)
+		      const int& minStapleDir,
+			  const int& skipStapleDir)
   {
     char listed_dirs[21]="";
     for(int mu=0;mu<NDIM;mu++)
@@ -32,7 +33,7 @@ namespace nissa
 	conf.updateEdges();
 	
 	PAR(0,locVol,
-	    CAPTURE(alpha,dirs,minStapleDir,
+	    CAPTURE(alpha,dirs,minStapleDir,skipStapleDir,
 		    TO_READ(conf),
 		    TO_WRITE(smearConf)),
 	    ivol,
@@ -44,10 +45,10 @@ namespace nissa
 		    su3 stap,temp1,temp2;
 		    su3_put_to_zero(stap);
 		    
-		    for(int nu=minStapleDir;nu<NDIM;nu++)     //  E---F---C
-		      if(nu!=mu)                              //  |   |   | mu
-			{                                     //  D---A---B
-			  const int A=ivol;                   //   nu
+		    for(int nu=minStapleDir;nu<NDIM;nu++)     							//  E---F---C
+		      if(nu!=mu and (skipStapleDir<0 || nu!=skipStapleDir))             //  |   |   | mu
+			{                                     							//  D---A---B
+			  const int A=ivol;                  							//   nu
 			  const int B=loclxNeighup[A][nu];
 			  const int F=loclxNeighup[A][mu];
 			  unsafe_su3_prod_su3(temp1,conf[A][nu],conf[B][mu]);
diff --git a/src/operations/smearing/APE.hpp b/src/operations/smearing/APE.hpp
index 103eb1155..4231a7f1d 100644
--- a/src/operations/smearing/APE.hpp
+++ b/src/operations/smearing/APE.hpp
@@ -67,7 +67,8 @@ namespace nissa
 		      const double& alpha,
 		      const int& nSteps,
 		      const WhichDirs& dirs=allDirs,
-		      const int& minStapleDir=0);
+		      const int& minStapleDir=0,
+          const int& skipStapleDir=-1);
   
   inline void ape_single_dir_smear_conf(LxField<quad_su3>& smear_conf,
 					const LxField<quad_su3>& origi_conf,
@@ -88,6 +89,17 @@ namespace nissa
   {
     ape_smear_conf(smear_conf,origi_conf,alpha,nstep,allOtherDirs[mu]);
   }
+
+  // - smear only mu != corr_dir
+  // - staples use only nu != corr_dir
+  inline void ape_orth_dir_smear_conf(LxField<quad_su3>& smear_conf,
+                                      const LxField<quad_su3>& origi_conf,
+                                      const double& alpha,
+                                      const int& nstep,
+                                      const int& corr_dir)
+  {
+    ape_smear_conf(smear_conf,origi_conf,alpha,nstep,allOtherDirs[corr_dir],0,corr_dir);
+  }
   
   inline void ape_temporal_smear_conf(LxField<quad_su3>& smear_conf,
 				      const LxField<quad_su3>& origi_conf,
diff --git a/src/operations/smearing/gaussian.hpp b/src/operations/smearing/gaussian.hpp
index 8d91870fa..3bd0fe84e 100644
--- a/src/operations/smearing/gaussian.hpp
+++ b/src/operations/smearing/gaussian.hpp
@@ -67,9 +67,9 @@ namespace nissa
 	if(ext_H==nullptr) _H=new LxField<TYPE>("H",WITH_HALO);	\
 	LxField<TYPE>& H=*_H;						\
   									\
-	const double norm_fact=1/(1+2*(kappa[1]+kappa[2]+kappa[3]));	\
+	const double norm_fact=1/(1+2*(kappa[0]+kappa[1]+kappa[2]+kappa[3]));	\
 									\
-	VERBOSITY_LV2_MASTER_PRINTF("GAUSSIAN smearing with kappa={%g,%g,%g}, %d iterations\n",kappa[1],kappa[2],kappa[3],niter); \
+	VERBOSITY_LV2_MASTER_PRINTF("GAUSSIAN smearing with kappa={%g,%g,%g,%g}, %d iterations\n",kappa[0],kappa[1],kappa[2],kappa[3],niter); \
 									\
 	/*iter 0*/							\
 	temp=origi_sc;							\
@@ -77,7 +77,7 @@ namespace nissa
 	/*loop over gaussian iterations*/				\
 	for(int iter=0;iter<niter;iter++)				\
 	  {								\
-	    VERBOSITY_LV3_MASTER_PRINTF("GAUSSIAN smearing with kappa={%g,%g,%g} iteration %d of %d\n",kappa[1],kappa[2],kappa[3],iter,niter); \
+	    VERBOSITY_LV3_MASTER_PRINTF("GAUSSIAN smearing with kappa={%g,%g,%g,%g} iteration %d of %d\n",kappa[0],kappa[1],kappa[2],kappa[3],iter,niter); \
 									\
 	    /*apply kappa*H*/						\
 	    NAME2(gaussian_smearing_apply_kappa_H,TYPE)(H,kappa,conf,temp); \
@@ -168,6 +168,35 @@ namespace nissa
     
     gaussian_smearing(smear_sc,origi_sc,conf,kappa,niter,ext_temp,ext_H);
   }
+
+  template <typename T>
+  void gaussian_smearing(LxField<T>& smear_sc,
+                         const LxField<T>& origi_sc,
+                         const LxField<quad_su3>& conf,
+                         double kappa_val,
+                         int niter,
+                         int corr_dir, // dir along which the correlator is measured
+                         LxField<T>* ext_temp=nullptr,
+                         LxField<T>* ext_H=nullptr)
+  {
+    Momentum kappa = {0.0, 0.0, 0.0, 0.0};
+    for(int mu=0; mu<NDIM; mu++)
+      if(mu != corr_dir)
+        kappa[mu] = kappa_val;
+
+    /* I've just added kappa[0] in the computation of norm above, it should work 
+    double norm_sum = 0.0;
+    for(int mu=0; mu<NDIM; mu++)
+      if(mu != corr_dir)
+        norm_sum += kappa[mu];
+    double norm_fact = 1.0 / (1.0 + 2.0 * norm_sum);
+	*/
+    
+    gaussian_smearing(smear_sc, origi_sc, conf, kappa, niter, ext_temp, ext_H);
+  }
+
+ 
+
 }
 
 #endif