diff --git a/current_generator/include/helspin.frm b/current_generator/include/helspin.frm
index 1c9f236..7c7b007 100644
--- a/current_generator/include/helspin.frm
+++ b/current_generator/include/helspin.frm
@@ -1,244 +1,185 @@
 */**
 *  \brief Procedures for dealing with helicity spinors
 *
 *  \authors   The HEJ collaboration (see AUTHORS for details)
 *  \date      2019
 *  \copyright GPLv2 or later
 */
 on shortstats;
 cf Conjugate,Current;
 cf angle,square;
 #define CHAINS "AngleChain,SquareChain,SpinorChain"
 ctensor `CHAINS';
 
 * internal symbols, rename?
 v p,q,k,l,p1;
 cf dum;
 ctensor ct;
 s x,y;
 i mu1,...,mu50;
 
-#procedure WriteOptimised(OUTPUT,EXPRESSION,NUMHELICITIES,?MOMENTA)
-
-   #write <`OUTPUT'> "template<%"
-   #define FIRST "1"
-   #do i=1,`NUMHELICITIES'
-      #ifdef `FIRST'
-         #undefine FIRST
-         #write <`OUTPUT'> "unsigned%"
-         #else
-         #write <`OUTPUT'> ", unsigned%"
-      #endif
-   #enddo
-   #write <`OUTPUT'> ">"
-   #write <`OUTPUT'> "std::complex<double> `EXPRESSION'("
-      #call WriteMomenta(`?MOMENTA')
-   #write <`OUTPUT'> "\n);\n"
-
-   #call WriteOptimisedHelper(`OUTPUT',`EXPRESSION',`NUMHELICITIES',,`?MOMENTA')
-
-#endprocedure
-
-#procedure WriteOptimisedHelper(OUTPUT,EXPRESSION,NUMHELICITIES,HELSTRING,?MOMENTA)
-
-   #if `NUMHELICITIES' > 0
-      #do H={+,-}
-         #call WriteOptimisedHelper(`OUTPUT',`EXPRESSION',{`NUMHELICITIES'-1},`HELSTRING'`H',`?MOMENTA')
-      #enddo
-
-      #else
-
-      #optimize [`EXPRESSION'`HELSTRING']
-      #write <`OUTPUT'> "inline template<>"
-*     the template specialization string is not syntactically correct
-*     and will be fixed later on by a small script
-      #write <`OUTPUT'> "std::complex<double> `EXPRESSION'<`HELSTRING'>("
-      #call WriteMomenta(`?MOMENTA')
-      #write <`OUTPUT'> "\n) {"
-      #write <`OUTPUT'> "%O"
-      #write <`OUTPUT'> "   return %E;" [`EXPRESSION'`HELSTRING']
-      #write <`OUTPUT'> "}" [`EXPRESSION'`HELSTRING']
-      #clearoptimize
-   #endif
-
-#endprocedure
-
-#procedure WriteMomenta(?MOMENTA)
-
-   #define FIRST "1"
-   #do P={`?MOMENTA'}
-      #ifdef `FIRST'
-         #undefine FIRST
-         #write <`OUTPUT'> "   CLHEP::HepLorentzVector const & `P'%"
-         #else
-         #write <`OUTPUT'> ",\n   CLHEP::HepLorentzVector const & `P'%"
-      #endif
-   #enddo
-
-#endprocedure
-
 #procedure ContractCurrents
 
    #call InsCurrent;
    .sort
 
    #call GetMomentaInChains($momenta)
 
    #if "`$momenta'" != ""
       #message Rewrite momenta `$momenta' into spinors
       #call ToSpinors(`$momenta')
    #endif
 
    while(match(SpinorChain(p?, mu1?, q?)*SpinorChain(k?, mu1?, l?)));
       #call Fierz(mu1?)
    endwhile;
    #call SortArgs
    id AngleChain(p?, q?) = angle(p, q);
    id SquareChain(p?, q?) = square(p, q);
    argument denom_;
       id AngleChain(p?, q?) = angle(p, q);
       id SquareChain(p?, q?) = square(p, q);
    endargument;
 
 #endprocedure
 
 #procedure GetMomentaInChains(MOMENTA)
 
    id ct?{`CHAINS'}(?a) = dum(?a)*ct(?a);
    #$tmp = 0;
    #$p = 0;
    argument dum;
       $p = term_;
       $tmp = $tmp + $p;
    endargument;
    id dum(?a) = 1;
    moduleoption local $p;
    moduleoption sum $tmp;
    .sort
    #$tmp = dum_($tmp);
    #inside $tmp
       splitarg dum_;
       argument dum_;
          dropcoefficient;
       endargument;
       repeat id dum_(?a, mu1?index_, ?b) = dum_(?a, ?b);
       id dum_(?a`MOMENTA') = 0;
    #endinside
 
 #endprocedure
 
 #procedure InsCurrent
 
    id Current(+1, ?a) = Current(-1, reverse_(?a));
 *  hack
 *  with ?a as pattern FORM has problems if the arguments are sums of momenta
    #do N=3,50
       id Current(-1, mu1?,...,mu`N'?) = SpinorChain(mu1,...,mu`N');
    #enddo
 
 #endprocedure
 
 #procedure SortArgs
 
    #do F={AngleChain,SquareChain}
       antisymmetrize `F':2;
    #enddo
 
 #endprocedure
 
 #procedure DoConjugate
 
    splitarg Conjugate;
    repeat id Conjugate(x?, y?, ?a) = Conjugate(x) + Conjugate(y, ?a);
    normalize Conjugate;
    id Conjugate(AngleChain(p?, q?)) = -SquareChain(p, q);
    id Conjugate(SquareChain(p?, q?)) = -AngleChain(p, q);
    id Conjugate(SpinorChain(?a)) = dum(reverse_(?a));
    id dum(?a) = SpinorChain(?a);
 
 #endprocedure
 
 #procedure Fierz(MU)
 
    once SpinorChain(p?, `MU', q?)*SpinorChain(k?, `MU', l?) = (
       2*AngleChain(p, k)*SquareChain(l, q)
    );
 
 #endprocedure
 
 #procedure MomentumToSpinorsInSpinorChain(P)
 
    id SpinorChain(p?, ?a, `P', ?b, q?) = dum(
       mod_(nargs_(p, ?a),2), mod_(nargs_(?b, q),2),
       SpinorChain(p, ?a, `P', ?b, q)
    );
    id dum(1, 1, SpinorChain(p?, ?a, `P', ?b, q?)) = (
       AngleChain(p, ?a, `P')*SquareChain(`P', ?b, q)
    );
    id dum(0, 0, SpinorChain(p?, ?a, `P', ?b, q?)) = (
       SpinorChain(p, ?a, `P')*SpinorChain(`P', ?b, q)
    );
    if(count(dum,1)>0);
       print "even number of arguments in spinor chain in %t";
       exit;
    endif;
 
 #endprocedure
 
 #procedure MomentumToSpinorsInAngleChain(P)
 
    id AngleChain(p?, ?a, `P', ?b, q?) = dum(
       mod_(nargs_(p, ?a),2), mod_(nargs_(?b, q),2),
       AngleChain(p, ?a, `P', ?b, q)
    );
    id dum(1, 0, AngleChain(p?, ?a, `P', ?b, q?)) = (
       AngleChain(p, ?a, `P')*Conjugate(SpinorChain(`P', ?b, q))
    );
    id dum(0, 1, AngleChain(p?, ?a, `P', ?b, q?)) = (
       SpinorChain(p, ?a, `P')*AngleChain(`P', ?b, q)
    );
    if(count(dum,1)>0);
       print "odd number of arguments in angle chain in %t";
       exit;
    endif;
 
 #endprocedure
 
 #procedure MomentumToSpinorsInSquareChain(P)
 
    id SquareChain(p?, ?a, `P', ?b, q?) = dum(
       mod_(nargs_(p, ?a),2), mod_(nargs_(?b, q),2),
       SquareChain(p, ?a, `P', ?b, q)
    );
    id dum(1, 0, SquareChain(p?, ?a, `P', ?b, q?)) = (
       SquareChain(p, ?a, `P')*SpinorChain(`P', ?b, q)
    );
    id dum(0, 1, SquareChain(p?, ?a, `P', ?b, q?)) = (
       Conjugate(SpinorChain(p, ?a, `P'))*SquareChain(`P', ?b, q)
    );
    if(count(dum,1)>0);
       print "odd number of arguments in angle chain in %t";
       exit;
    endif;
 
 #endprocedure
 
 
 #procedure MomentumToSpinors(P)
 
    repeat;
       #do CHAIN={`CHAINS'}
          id `CHAIN'(?a, `P', `P', ?b) = 0;
          #call MomentumToSpinorsIn`CHAIN'(`P')
       #enddo
    endrepeat;
 
 #endprocedure
 
 #procedure ToSpinors(?PS)
 
    #do P={`?PS'}
       #call MomentumToSpinors(`P')
    #enddo
 
 #endprocedure
diff --git a/current_generator/include/write.frm b/current_generator/include/write.frm
new file mode 100644
index 0000000..2123f26
--- /dev/null
+++ b/current_generator/include/write.frm
@@ -0,0 +1,76 @@
+#procedure WriteHeader(OUTPUT)
+
+   #write <`OUTPUT'> "#pragma once"
+   #write <`OUTPUT'> ""
+   #write <`OUTPUT'> "#include <complex>"
+   #write <`OUTPUT'> "#include \"CLHEP/Vector/LorentzVector.h\""
+   #write <`OUTPUT'> "#include \"HEJ/Tensor.hh\""
+   #write <`OUTPUT'> ""
+   #write <`OUTPUT'> "namespace HEJ {"
+
+#endprocedure
+
+#procedure WriteOptimised(OUTPUT,EXPRESSION,NUMHELICITIES,?MOMENTA)
+
+   #write <`OUTPUT'> "   template<%"
+   #define FIRST "1"
+   #do i=1,`NUMHELICITIES'
+      #ifdef `FIRST'
+         #undefine FIRST
+         #write <`OUTPUT'> "Helicity%"
+         #else
+         #write <`OUTPUT'> ", Helicity%"
+      #endif
+   #enddo
+   #write <`OUTPUT'> ">"
+   #write <`OUTPUT'> "   std::complex<double> `EXPRESSION'("
+      #call WriteMomenta(`?MOMENTA')
+   #write <`OUTPUT'> "\n   );\n"
+
+   #call WriteOptimisedHelper(`OUTPUT',`EXPRESSION',`NUMHELICITIES',,`?MOMENTA')
+
+#endprocedure
+
+#procedure WriteOptimisedHelper(OUTPUT,EXPRESSION,NUMHELICITIES,HELSTRING,?MOMENTA)
+
+   #if `NUMHELICITIES' > 0
+      #do H={+,-}
+         #call WriteOptimisedHelper(`OUTPUT',`EXPRESSION',{`NUMHELICITIES'-1},`HELSTRING'`H',`?MOMENTA')
+      #enddo
+
+      #else
+
+      #optimize [`EXPRESSION'`HELSTRING']
+      #write <`OUTPUT'> "   inline template<>"
+*     the template specialization string is not syntactically correct
+*     and will be fixed later on by a small script
+      #write <`OUTPUT'> "   std::complex<double> `EXPRESSION'<`HELSTRING'>("
+      #call WriteMomenta(`?MOMENTA')
+      #write <`OUTPUT'> "\n   ) {"
+      #write <`OUTPUT'> "   %O"
+      #write <`OUTPUT'> "      return %E;" [`EXPRESSION'`HELSTRING']
+      #write <`OUTPUT'> "   }\n" [`EXPRESSION'`HELSTRING']
+      #clearoptimize
+   #endif
+
+#endprocedure
+
+#procedure WriteMomenta(?MOMENTA)
+
+   #define FIRST "1"
+   #do P={`?MOMENTA'}
+      #ifdef `FIRST'
+         #undefine FIRST
+         #write <`OUTPUT'> "      CLHEP::HepLorentzVector const & `P'%"
+         #else
+         #write <`OUTPUT'> ",\n      CLHEP::HepLorentzVector const & `P'%"
+      #endif
+   #enddo
+
+#endprocedure
+
+#procedure WriteFooter(OUTPUT)
+
+   #write <`OUTPUT'> "}"
+
+#endprocedure
diff --git a/current_generator/jW_uno.frm b/current_generator/jW_uno.frm
index e2750fd..147af1a 100644
--- a/current_generator/jW_uno.frm
+++ b/current_generator/jW_uno.frm
@@ -1,78 +1,81 @@
 */**
 *  \brief Calculation of W unordered current
 *
 *  \authors   The HEJ collaboration (see AUTHORS for details)
 *  \date      2019
 *  \copyright GPLv2 or later
 */
 #include- include/helspin.frm
+#include- include/write.frm
 
 s h,s1g,s1W,s1gW,s2g,sbg,tag,taW,taW1,taWg;
 v p,p1,p2,pa,pb,pg,pl,plbar,pW,pr,q1,q1g;
 i mu,nu,rho,sigma;
 i mu1,...,mu20;
 cf epsW,epsg,j2b,m2,sqrt;
 
 #do h1={+,-}
    #do h2={+,-}
       #do hg={+,-}
 *        eq:U1tensor in developer manual
          l [U1 `h1'`h2'`hg'] = epsW(-1, rho)*epsg(`hg'1, nu)*j2b(`h2'1, mu)*(
             + Current(`h1'1, p1, nu, p1+pg, mu, pa-pW, rho, pa)/(s1g*taW)
             + Current(`h1'1, p1, nu, p1+pg, rho, p1+pg+pW, mu, pa)/(s1g*s1gW)
             + Current(`h1'1, p1, rho, p1+pW, nu, p1+pg+pW, mu, pa)/(s1W*s1gW)
          );
 
 *        eq:U2tensor in developer manual
          l [U2 `h1'`h2'`hg'] = epsW(-1, rho)*epsg(`hg'1, nu)*j2b(`h2'1, mu)*(
             + Current(`h1'1, p1, mu, pa - pW - pg, nu, pa - pW, rho, pa)/(taW*taWg)
             + Current(`h1'1, p1, mu, pa - pW - pg, rho, pa - pg, nu, pa)/(tag*taWg)
             + Current(`h1'1, p1, rho, p1 + pW, mu, pa - pg, nu, pa)/(s1W*tag)
          );
 
 *        eq:Ltensor in developer manual
          l [L `h1'`h2'`hg'] = epsW(-1, rho)*epsg(`hg'1, nu)*j2b(`h2'1, mu)*(
             Current(`h1'1, p1, sigma, pa - pW, rho, pa)/taW +
             Current(`h1'1, p1, rho, p1 + pW, sigma, pa)/s1W
          )*(
             ((pb(nu)/sbg + p2(nu)/s2g)*m2(q1g) + 2*q1(nu) - pg(nu))*d_(mu, sigma)
             - 2*pg(mu)*d_(nu, sigma)
             + (2*pg(sigma) - q1(sigma))*d_(mu, nu)
          )/taW1;
       #enddo
    #enddo
 #enddo
 
 * choice of best reference vector (p2 or pb)
 id epsg(h?, mu?)*j2b(h?, nu?) = epsg(h, mu, p2)*j2b(h, nu);
 also epsg(h?, mu?) = epsg(h, mu, pb);
 
 id epsW(h?, mu?) = Current(h, pl, mu, plbar);
 
 id epsg(-1, mu?, pr?) = sqrt(2)/2*SpinorChain(pr, mu, pg)/AngleChain(pg,pr);
 id epsg(+1, mu?, pr?) = sqrt(2)/2*SpinorChain(pg, mu, pr)/SquareChain(pg,pr);
 id j2b(h?, mu?) = Current(h, p2, mu, pb);
 
 multiply replace_(q1g,q1-pg);
 multiply replace_(q1,pa-p1-pW);
 multiply replace_(pW,pl+plbar);
 
 .sort
 #call ContractCurrents
 multiply replace_(
    s1g,m2(p1+pg),
    s1W,m2(p1+pl+plbar),
    s1gW,m2(p1+pg+pl+plbar),
    s2g,m2(p2+pg),
    tag,m2(pa-pg),
    taW,m2(pa-pl-plbar),
    taW1,m2(pa-pl-plbar-p1),
    taW1,m2(pa-pl-plbar-pg)
 );
 .sort
 format O4;
 format c;
+#call WriteHeader(`OUTPUT')
 #call WriteOptimised(`OUTPUT',U1,3,p1,p2,pa,pb,pg,pl,plbar)
 #call WriteOptimised(`OUTPUT',U2,3,p1,p2,pa,pb,pg,pl,plbar)
 #call WriteOptimised(`OUTPUT',L,3,p1,p2,pa,pb,pg,pl,plbar)
+#call WriteFooter(`OUTPUT')
 .end