diff --git a/Higgs/loop_sm_c3d4.txt b/Higgs/loop_sm_c3d4.txt
new file mode 100644
index 0000000000000000000000000000000000000000..be11ebbb7545d59e04a4c390cb234e3a5c2ed385
--- /dev/null
+++ b/Higgs/loop_sm_c3d4.txt
@@ -0,0 +1,5 @@
+Requestor: Nicholas Kyriacou. Model from Marko Stamenkovic in collaboration with others (Andreas Papaefstathiou et. al.)
+Contents: Model to Generate HHH Production. Allows to independently vary individual Higgs Coupling Parameters utitlized for Kappa-Scans. 
+Paper: https://arxiv.org/abs/2312.13562
+Source: https://github.com/mstamenk/triple-h-mc-gen/tree/main
+JIRA: https://its.cern.ch/jira/browse/ATLMCPROD-11023
diff --git a/Higgs/loop_sm_c3d4/CT_couplings.py b/Higgs/loop_sm_c3d4/CT_couplings.py
new file mode 100755
index 0000000000000000000000000000000000000000..056057dcd1d20cd31605336c19dcd59512cbd412
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/CT_couplings.py
@@ -0,0 +1,718 @@
+# This file was automatically created by FeynRules $Revision: 535 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Fri 18 Mar 2011 18:40:51
+
+from object_library import all_couplings, Coupling
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec
+
+################
+# R2 couplings #
+################
+
+# ========= #
+# Pure QCD  #
+# ========= #
+
+R2_3Gq = Coupling(name = 'R2_3Gq',
+                value = '2.0*G**3/(48.0*cmath.pi**2)',
+                order = {'QCD':3})
+
+R2_3Gg = Coupling(name = 'R2_3Gg',
+                 value = 'Ncol*G**3/(48.0*cmath.pi**2)*(7.0/4.0+lhv)',
+                 order = {'QCD':3})
+
+#=============================================================================================
+#  4-gluon R2 couplings
+#=============================================================================================
+
+# Gluon contribution to it
+
+GC_4GR2_Gluon_delta5 = Coupling(name = 'GC_4GR2_Gluon_delta5',
+                 value = '-4.0*complex(0,1)*RGR2*(2.0*lhv+5.0)',
+                 order = {'QCD':4})
+
+GC_4GR2_Gluon_delta7 = Coupling(name = 'GC_4GR2_Gluon_delta7',
+                 value = '2.0*complex(0,1)*RGR2*(2.0*lhv+7.0)',
+                 order = {'QCD':4})
+
+GC_4GR2_2Struct = Coupling(name = 'GC_4GR2_2Struct',
+                 value = '2.0*complex(0,1)*RGR2*Ncol*(lhv+3.0)',
+                 order = {'QCD':4})
+
+GC_4GR2_4Struct = Coupling(name = 'GC_4GR2_4Struct',
+                 value = '-complex(0,1)*RGR2*Ncol*(4.0*lhv+11.0)',
+                 order = {'QCD':4})
+
+# Fermion contribution to it
+
+GC_4GR2_Fermion_delta5 = Coupling(name = 'GC_4GR2_Fermion_delta5',
+                 value = '(2.0/Ncol)*5.0*complex(0,1)*RGR2',
+                 order = {'QCD':4})
+
+GC_4GR2_Fermion_delta11 = Coupling(name = 'GC_4GR2_Fermion_delta11',
+                 value = '-(2.0/Ncol)*11.0*complex(0,1)*RGR2',
+                 order = {'QCD':4})
+
+GC_4GR2_5Struct = Coupling(name = 'GC_4GR2_5Struct',
+                 value = '5.0*complex(0,1)*RGR2',
+                 order = {'QCD':4})
+
+GC_4GR2_11Struct = Coupling(name = 'GC_4GR2_11Struct',
+                 value = '-11.0*complex(0,1)*RGR2',
+                 order = {'QCD':4})
+
+# From the auto UFO from FR
+
+R2GC_137_43 = Coupling(name = 'R2GC_137_43',
+                       value = '(complex(0,1)*G**4)/(192.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_137_44 = Coupling(name = 'R2GC_137_44',
+                       value = '-(complex(0,1)*G**4)/(64.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_138_45 = Coupling(name = 'R2GC_138_45',
+                       value = '-(complex(0,1)*G**4)/(192.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_138_46 = Coupling(name = 'R2GC_138_46',
+                       value = '(complex(0,1)*G**4)/(64.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_139_47 = Coupling(name = 'R2GC_139_47',
+                       value = '-(complex(0,1)*G**4)/(48.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_140_48 = Coupling(name = 'R2GC_140_48',
+                       value = '(complex(0,1)*G**4)/(288.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_140_49 = Coupling(name = 'R2GC_140_49',
+                       value = '-(complex(0,1)*G**4)/(32.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_141_50 = Coupling(name = 'R2GC_141_50',
+                       value = '-(complex(0,1)*G**4)/(16.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_141_51 = Coupling(name = 'R2GC_141_51',
+                       value = '(-7*complex(0,1)*G**4)/(32.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_142_52 = Coupling(name = 'R2GC_142_52',
+                       value = '(complex(0,1)*G**4)/(16.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_142_53 = Coupling(name = 'R2GC_142_53',
+                       value = '(7*complex(0,1)*G**4)/(32.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_143_54 = Coupling(name = 'R2GC_143_54',
+                       value = '(11*complex(0,1)*G**4)/(192.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_143_55 = Coupling(name = 'R2GC_143_55',
+                       value = '(15*complex(0,1)*G**4)/(64.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_144_56 = Coupling(name = 'R2GC_144_56',
+                       value = '(-3*complex(0,1)*G**4)/(64.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_144_57 = Coupling(name = 'R2GC_144_57',
+                       value = '(-17*complex(0,1)*G**4)/(64.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_145_58 = Coupling(name = 'R2GC_145_58',
+                       value = '-G**4/(192.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+R2GC_145_59 = Coupling(name = 'R2GC_145_59',
+                       value = 'G**4/(64.*cmath.pi**2)',
+                       order = {'QCD':4})
+
+#=============================================================================================
+
+R2_GQQ = Coupling(name = 'R2_GQQ',
+                value = '-complex(0,1)*G**3/(16.0*cmath.pi**2)*((Ncol**2-1)/(2.0*Ncol))*(1.0+lhv)',
+                order = {'QCD':3})
+
+R2_GGq = Coupling(name = 'R2_GGq',
+                 value = '(2.0)*complex(0,1)*G**2/(48.0*cmath.pi**2)',
+                 order = {'QCD':2})
+
+R2_GGb = Coupling(name = 'R2_GGb',
+                 value = '(2.0)*complex(0,1)*G**2*(-6.0*MB**2)/(48.0*cmath.pi**2)',
+                 order = {'QCD':2})
+
+R2_GGc = Coupling(name = 'R2_GGc',
+                 value = '(2.0)*complex(0,1)*G**2*(-6.0*MC**2)/(48.0*cmath.pi**2)',
+                 order = {'QCD':2})
+
+R2_GGt = Coupling(name = 'R2_GGt',
+                 value = '(2.0)*complex(0,1)*G**2*(-6.0*MT**2)/(48.0*cmath.pi**2)',
+                 order = {'QCD':2})
+
+R2_GGg_1 = Coupling(name = 'R2_GGg_1',
+                 value = '(2.0)*complex(0,1)*G**2*Ncol/(48.0*cmath.pi**2)*(1.0/2.0+lhv)',
+                 order = {'QCD':2})
+
+R2_GGg_2 = Coupling(name = 'R2_GGg_2',
+                 value = '-(2.0)*complex(0,1)*G**2*Ncol/(48.0*cmath.pi**2)*lhv',
+                 order = {'QCD':2})
+
+R2_QQq = Coupling(name = 'R2_QQq',
+                 value =  'lhv*complex(0,1)*G**2*(Ncol**2-1)/(32.0*cmath.pi**2*Ncol)',
+                 order = {'QCD':2})
+
+R2_QQc = Coupling(name = 'R2_QQc',
+                 value =  'lhv*complex(0,1)*G**2*(Ncol**2-1)*(2.0*MC)/(32.0*cmath.pi**2*Ncol)',
+                 order = {'QCD':2})
+
+R2_QQb = Coupling(name = 'R2_QQb',
+                 value =  'lhv*complex(0,1)*G**2*(Ncol**2-1)*(2.0*MB)/(32.0*cmath.pi**2*Ncol)',
+                 order = {'QCD':2})
+
+R2_QQt = Coupling(name = 'R2_QQt',
+                 value =  'lhv*complex(0,1)*G**2*(Ncol**2-1)*(2.0*MT)/(32.0*cmath.pi**2*Ncol)',
+                 order = {'QCD':2})
+
+# ============== #
+# Mixed QCD-QED  #
+# ============== #
+
+# Quark couplings to A and Z
+
+R2_DDA = Coupling(name = 'R2_DDA',
+                value = '(-(ee*complex(0,1))/3.0)*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_UUA = Coupling(name = 'R2_UUA',
+                value = '(2.0*(ee*complex(0,1))/3.0)*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_DDZ_V2 = Coupling(name = 'R2_DDZ_V2',
+                value = '(-(cw*ee*complex(0,1))/(2.0*sw))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_DDZ_V3 = Coupling(name = 'R2_DDZ_V3',
+                value = '(-(ee*complex(0,1)*sw)/(6.0*cw))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_UUZ_V2 = Coupling(name = 'R2_UUZ_V2',
+                value = '((cw*ee*complex(0,1))/(2.*sw))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_UUZ_V5 = Coupling(name = 'R2_UUZ_V5',
+                value = '(-(ee*complex(0,1)*sw)/(6.*cw))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+# Quark couplings to W, most general CKM
+
+R2_dxuW = Coupling(name = 'R2_dxuW',
+                value = '((CKM11*ee*complex(0,1))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_dxcW = Coupling(name = 'R2_dxcW',
+                value = '((CKM21*ee*complex(0,1))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_dxtW = Coupling(name = 'R2_dxtW',
+                value = '((CKM31*ee*complex(0,1))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_sxuW = Coupling(name = 'R2_sxuW',
+                value = '((CKM12*ee*complex(0,1))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_sxcW = Coupling(name = 'R2_sxcW',
+                value = '((CKM22*ee*complex(0,1))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_sxtW = Coupling(name = 'R2_sxtW',
+                value = '((CKM32*ee*complex(0,1))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_bxuW = Coupling(name = 'R2_bxuW',
+                value = '((CKM13*ee*complex(0,1))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_bxcW = Coupling(name = 'R2_bxcW',
+                value = '((CKM23*ee*complex(0,1))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_bxtW = Coupling(name = 'R2_bxtW',
+                value = '((CKM33*ee*complex(0,1))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_uxdW = Coupling(name = 'R2_uxdW',
+                value = '((ee*complex(0,1)*complexconjugate(CKM11))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_cxdW = Coupling(name = 'R2_cxdW',
+                value = '((ee*complex(0,1)*complexconjugate(CKM21))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_txdW = Coupling(name = 'R2_txdW',
+                value = '((ee*complex(0,1)*complexconjugate(CKM31))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_uxsW = Coupling(name = 'R2_uxsW',
+                value = '((ee*complex(0,1)*complexconjugate(CKM12))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_cxsW = Coupling(name = 'R2_cxsW',
+                value = '((ee*complex(0,1)*complexconjugate(CKM22))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_txsW = Coupling(name = 'R2_txsW',
+                value = '((ee*complex(0,1)*complexconjugate(CKM32))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_uxbW = Coupling(name = 'R2_uxbW',
+                value = '((ee*complex(0,1)*complexconjugate(CKM13))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_cxbW = Coupling(name = 'R2_cxbW',
+                value = '((ee*complex(0,1)*complexconjugate(CKM23))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+R2_txbW = Coupling(name = 'R2_txbW',
+                value = '((ee*complex(0,1)*complexconjugate(CKM33))/(sw*cmath.sqrt(2)))*R2MixedFactor',
+                order = {'QCD':2,'QED':1})
+
+# R2 for the Higgs couplings to massive quarks 
+
+R2_bbH = Coupling(name = 'R2_bbH',
+                value = '(-((complex(0,1)*yb)/cmath.sqrt(2)))*(2.0*R2MixedFactor)',
+                order = {'QCD':2,'QED':1})
+
+R2_ttH = Coupling(name = 'R2_ttH',
+                value = '(-((complex(0,1)*yt)/cmath.sqrt(2)))*(2.0*R2MixedFactor)',
+                order = {'QCD':2,'QED':1})
+
+R2_ccH = Coupling(name = 'R2_ccH',
+                value = '(-((complex(0,1)*yc)/cmath.sqrt (2)))*(2.0*R2MixedFactor)',
+                order = {'QCD':2,'QED':1})
+
+# R2 interactions non proportional to the SM
+
+# R2 for the Higgs interactions
+R2_GGHc = Coupling(name = 'R2_GGHc',
+                value = '4.0*(-((complex(0,1)*yc)/cmath.sqrt(2)))*(1.0/2.0)*(G**2/(8.0*cmath.pi**2))*MC',
+                order = {'QCD':2,'QED':1})
+
+R2_GGHb = Coupling(name = 'R2_GGHb',
+                value = '4.0*(-((complex(0,1)*yb)/cmath.sqrt(2)))*(1.0/2.0)*(G**2/(8.0*cmath.pi**2))*MB',
+                order = {'QCD':2,'QED':1})
+
+R2_GGHt = Coupling(name = 'R2_GGHt',
+                value = '4.0*(-((complex(0,1)*yt)/cmath.sqrt(2)))*(1.0/2.0)*(G**2/(8.0*cmath.pi**2))*MT',
+                order = {'QCD':2,'QED':1})
+
+# R2 for the weak vector bosons interaction with gluons
+
+R2_GGZup = Coupling(name = 'R2_GGZup',
+                    value = {0:'2.0*AxialZUp*(G**2/(12.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':1})
+
+R2_GGZdown = Coupling(name = 'R2_GGZdown',
+                    value = {0:'2.0*AxialZDown*(G**2/(12.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':1})
+
+# EDIT VH
+# There is a factor four added here to all the GGVV R2 couplings below without proper understanding of it
+# I have not investigate too much about it though.
+# END EDIT VH
+
+R2_GGZAup = Coupling(name = 'R2_GGZAup',
+                    value = {0:'4.0*(-VectorAUp*VectorZUp)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGZAdown = Coupling(name = 'R2_GGZAdown',
+                    value = {0:'4.0*(-VectorADown*VectorZDown)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGZZdown = Coupling(name = 'R2_GGZZdown',
+                    value = {0:'4.0*(-VectorZDown*VectorZDown-AxialZDown*AxialZDown)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGZZup = Coupling(name = 'R2_GGZZup',
+                    value = {0:'4.0*(-VectorZUp*VectorZUp-AxialZUp*AxialZUp)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGAAdown = Coupling(name = 'R2_GGAAdown',
+                    value = {0:'4.0*(-VectorADown*VectorADown)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGAAup = Coupling(name = 'R2_GGAAup',
+                    value = {0:'4.0*(-VectorAUp*VectorAUp)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGWWud = Coupling(name = 'R2_GGWWud',
+                    value = {0:'4.0*(CKM11*complexconjugate(CKM11))*(-VectorWmDxU*VectorWpUxD-AxialWmDxU*AxialWpUxD)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGWWus = Coupling(name = 'R2_GGWWus',
+                    value = {0:'4.0*(CKM12*complexconjugate(CKM12))*(-VectorWmDxU*VectorWpUxD-AxialWmDxU*AxialWpUxD)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGWWub = Coupling(name = 'R2_GGWWub',
+                    value = {0:'4.0*(CKM13*complexconjugate(CKM13))*(-VectorWmDxU*VectorWpUxD-AxialWmDxU*AxialWpUxD)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGWWcd = Coupling(name = 'R2_GGWWcd',
+                    value = {0:'4.0*(CKM21*complexconjugate(CKM21))*(-VectorWmDxU*VectorWpUxD-AxialWmDxU*AxialWpUxD)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGWWcs = Coupling(name = 'R2_GGWWcs',
+                    value = {0:'4.0*(CKM22*complexconjugate(CKM22))*(-VectorWmDxU*VectorWpUxD-AxialWmDxU*AxialWpUxD)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGWWcb = Coupling(name = 'R2_GGWWcb',
+                    value = {0:'4.0*(CKM23*complexconjugate(CKM23))*(-VectorWmDxU*VectorWpUxD-AxialWmDxU*AxialWpUxD)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGWWtd = Coupling(name = 'R2_GGWWtd',
+                    value = {0:'4.0*(CKM31*complexconjugate(CKM31))*(-VectorWmDxU*VectorWpUxD-AxialWmDxU*AxialWpUxD)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGWWts = Coupling(name = 'R2_GGWWts',
+                    value = {0:'4.0*(CKM32*complexconjugate(CKM32))*(-VectorWmDxU*VectorWpUxD-AxialWmDxU*AxialWpUxD)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGWWtb = Coupling(name = 'R2_GGWWtb',
+                    value = {0:'4.0*(CKM33*complexconjugate(CKM33))*(-VectorWmDxU*VectorWpUxD-AxialWmDxU*AxialWpUxD)*(1.0/2.0)*(-(complex(0,1)*G**2)/(24.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGHHc = Coupling(name = 'R2_GGHHc',
+                    value = {0:'4.0*(-yc**2/2.0)*(1.0/2.0)*((complex(0,1)*G**2)/(8.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGHHb = Coupling(name = 'R2_GGHHb',
+                    value = {0:'4.0*(-yb**2/2.0)*(1.0/2.0)*((complex(0,1)*G**2)/(8.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGHHt = Coupling(name = 'R2_GGHHt',
+                    value = {0:'4.0*(-yt**2/2.0)*(1.0/2.0)*((complex(0,1)*G**2)/(8.0*cmath.pi**2))'},
+                    order = {'QCD':2,'QED':2})
+
+R2_GGGZvecUp = Coupling(name = 'R2_GGGZvecUp',
+                        value = {0:'complex(0,1)*VectorZUp*(-1.0/2.0)*(-G**3/(24.0*cmath.pi**2))'},
+                        order = {'QCD':3,'QED':1})
+
+R2_GGGZvecDown = Coupling(name = 'R2_GGGZvecDown',
+                        value = {0:'complex(0,1)*VectorZDown*(-1.0/2.0)*(-G**3/(24.0*cmath.pi**2))'},
+                        order = {'QCD':3,'QED':1})
+
+R2_GGGZaxialUp = Coupling(name = 'R2_GGGZaxialUp',
+                        value = {0:'complex(0,1)*AxialZUp*(-9.0/2.0)*(-G**3/(24.0*cmath.pi**2))'},
+                        order = {'QCD':3,'QED':1})
+
+R2_GGGZaxialDown = Coupling(name = 'R2_GGGZaxialDown',
+                        value = {0:'complex(0,1)*AxialZDown*(-9.0/2.0)*(-G**3/(24.0*cmath.pi**2))'},
+                        order = {'QCD':3,'QED':1})
+
+R2_GGGAvecUp = Coupling(name = 'R2_GGGAvecUp',
+                        value = {0:'complex(0,1)*VectorAUp*(-1.0/2.0)*(-G**3/(24.0*cmath.pi**2))'},
+                        order = {'QCD':3,'QED':1})
+R2_GGGAvecDown = Coupling(name = 'R2_GGGAvecDown',
+                        value = {0:'complex(0,1)*VectorADown*(-1.0/2.0)*(-G**3/(24.0*cmath.pi**2))'},
+                        order = {'QCD':3,'QED':1})
+
+################
+# UV couplings #
+################
+
+# ========= #
+# Pure QCD  #
+# ========= #
+
+UV_3Gg = Coupling(name = 'UV_3Gg',
+                 value = '-G_UVg*G',
+                 order = {'QCD':3})
+
+UV_3Gq = Coupling(name = 'UV_3Gq',
+                 value = '-G_UVq*G',
+                 order = {'QCD':3})
+
+UV_3Gc = Coupling(name = 'UV_3Gc',
+                 value = '-G_UVc*G',
+                 order = {'QCD':3})
+
+UV_3Gb = Coupling(name = 'UV_3Gb',
+                 value = '-G_UVb*G',
+                 order = {'QCD':3})
+
+UV_3Gt = Coupling(name = 'UV_3Gt',
+                 value = '-G_UVt*G',
+                 order = {'QCD':3})
+
+UV_4Gg = Coupling(name = 'UV_4Gg',
+                 value = '2.0*complex(0,1)*G_UVg*(G**2)',
+                 order = {'QCD':4})
+
+UV_4Gq = Coupling(name = 'UV_4Gq',
+                 value = '2.0*complex(0,1)*G_UVq*(G**2)',
+                 order = {'QCD':4})
+
+UV_4Gc = Coupling(name = 'UV_4Gc',
+                 value = '2.0*complex(0,1)*G_UVc*(G**2)',
+                 order = {'QCD':4})
+
+UV_4Gb = Coupling(name = 'UV_4Gb',
+                 value = '2.0*complex(0,1)*G_UVb*(G**2)',
+                 order = {'QCD':4})
+
+UV_4Gt = Coupling(name = 'UV_4Ggt',
+                 value = '2.0*complex(0,1)*G_UVt*(G**2)',
+                 order = {'QCD':4})
+
+UV_GQQg = Coupling(name = 'UV_GQQg',
+                 value = 'complex(0,1)*G_UVg*G',
+                 order = {'QCD':3})
+
+UV_GQQq = Coupling(name = 'UV_GQQq',
+                 value = 'complex(0,1)*G_UVq*G',
+                 order = {'QCD':3})
+
+UV_GQQc = Coupling(name = 'UV_GQQc',
+                 value = 'complex(0,1)*G_UVc*G',
+                 order = {'QCD':3})
+
+UV_GQQb = Coupling(name = 'UV_GQQb',
+                 value = 'complex(0,1)*G_UVb*G',
+                 order = {'QCD':3})
+
+UV_GQQt = Coupling(name = 'UV_GQQt',
+                 value = 'complex(0,1)*G_UVt*G',
+                 order = {'QCD':3})
+
+UV_cMass = Coupling(name = 'UV_cMass',
+                 value = 'cMass_UV',
+                 order = {'QCD':2})
+
+UV_bMass = Coupling(name = 'UV_bMass',
+                 value = 'bMass_UV',
+                 order = {'QCD':2}) 
+
+UV_tMass = Coupling(name = 'UV_tMass',
+                 value = 'tMass_UV',
+                 order = {'QCD':2})
+
+# ============== #
+# Mixed QCD-QED  #
+# ============== #
+
+UV_Hcc = Coupling(name = 'UV_Hcc',
+                 value = '-((complex(0,1)*yc)/cmath.sqrt(2))*UV_yuk_c',
+                 order = {'QED':1,'QCD':2})
+
+UV_Htt = Coupling(name = 'UV_Htt',
+                 value = '-((complex(0,1)*yt)/cmath.sqrt(2))*UV_yuk_t',
+                 order = {'QED':1,'QCD':2})
+
+UV_Hbb = Coupling(name = 'UV_Hbb',
+                 value = '-((complex(0,1)*yb)/cmath.sqrt(2))*UV_yuk_b',
+                 order = {'QED':1,'QCD':2})
+
+# ============================== #
+# Goldstone R2 CT couplings      #
+# ============================== #
+
+R2_GGGpGm_ub = Coupling(name = 'R2_GGGpGm_ub',
+                 value = 'complex(0,1)*R2_GGGpGm_factor*(Vector_ubGp*Vector_ubGm-Axial_ubGp*Axial_ubGm)',
+                 order = {'QED':2,'QCD':2})
+
+R2_GGGpGm_cd = Coupling(name = 'R2_GGGpGm_cd',
+                 value = 'complex(0,1)*R2_GGGpGm_factor*(Vector_cdGp*Vector_cdGm-Axial_cdGp*Axial_cdGm)',
+                 order = {'QED':2,'QCD':2})
+
+R2_GGGpGm_cs = Coupling(name = 'R2_GGGpGm_cs',
+                 value = 'complex(0,1)*R2_GGGpGm_factor*(Vector_csGp*Vector_csGm-Axial_csGp*Axial_csGm)',
+                 order = {'QED':2,'QCD':2})
+
+R2_GGGpGm_cb = Coupling(name = 'R2_GGGpGm_cb',
+                 value = 'complex(0,1)*R2_GGGpGm_factor*(Vector_cbGp*Vector_cbGm-Axial_cbGp*Axial_cbGm)',
+                 order = {'QED':2,'QCD':2})
+
+R2_GGGpGm_td = Coupling(name = 'R2_GGGpGm_td',
+                 value = 'complex(0,1)*R2_GGGpGm_factor*(Vector_tdGp*Vector_tdGm-Axial_tdGp*Axial_tdGm)',
+                 order = {'QED':2,'QCD':2})
+
+R2_GGGpGm_ts = Coupling(name = 'R2_GGGpGm_ts',
+                 value = 'complex(0,1)*R2_GGGpGm_factor*(Vector_tsGp*Vector_tsGm-Axial_tsGp*Axial_tsGm)',
+                 order = {'QED':2,'QCD':2})
+
+R2_GGGpGm_tb = Coupling(name = 'R2_GGGpGm_tb',
+                 value = 'complex(0,1)*R2_GGGpGm_factor*(Vector_tbGp*Vector_tbGm-Axial_tbGp*Axial_tbGm)',
+                 order = {'QED':2,'QCD':2})
+
+R2_GGG0G0_c = Coupling(name = 'R2_GGG0G0_c',
+                       value = 'complex(0,1)*R2_GGG0G0_factor*(-(1.0/2.0)*yc**2)',
+                       order = {'QED':2,'QCD':2})
+
+R2_GGG0G0_b = Coupling(name = 'R2_GGG0G0_b',
+                       value = 'complex(0,1)*R2_GGG0G0_factor*(-(1.0/2.0)*yb**2)',
+                       order = {'QED':2,'QCD':2})
+
+R2_GGG0G0_t = Coupling(name = 'R2_GGG0G0_t',
+                       value = 'complex(0,1)*R2_GGG0G0_factor*(-(1.0/2.0)*yt**2)',
+                       order = {'QED':2,'QCD':2})
+
+GC_R2_1015 = Coupling(name = 'GC_R2_1015',
+                 value = 'I1x33*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1021 = Coupling(name = 'GC_R2_1021',
+                 value = '-I2x33*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1027 = Coupling(name = 'GC_R2_1027',
+                 value = 'I3x33*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1030 = Coupling(name = 'GC_R2_1030',
+                 value = '-I4x33*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1082 = Coupling(name = 'GC_R2_1082',
+                 value = '-(yb/cmath.sqrt(2))*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1085 = Coupling(name = 'GC_R2_1085',
+                 value = '(yc/cmath.sqrt(2))*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1095 = Coupling(name = 'GC_R2_1095',
+                 value = '(yt/cmath.sqrt(2))*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1018 = Coupling(name = 'GC_R2_1018',
+                 value = '(-I2x22)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1019 = Coupling(name = 'GC_R2_1019',
+                 value = '(-I2x23)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1013 = Coupling(name = 'GC_R2_1013',
+                 value = '(I1x31)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1014 = Coupling(name = 'GC_R2_1014',
+                 value = '(I1x32)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1020 = Coupling(name = 'GC_R2_1020',
+                 value = '(-I2x32)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1028 = Coupling(name = 'GC_R2_1028',
+                 value = '(-I4x13)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1022 = Coupling(name = 'GC_R2_1022',
+                 value = '(I3x21)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1023 = Coupling(name = 'GC_R2_1023',
+                 value = '(I3x22)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1024 = Coupling(name = 'GC_R2_1024',
+                 value = '(I3x23)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1029 = Coupling(name = 'GC_R2_1029',
+                 value = '(-I4x23)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1025 = Coupling(name = 'GC_R2_1025',
+                 value = '(I3x31)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+GC_R2_1026 = Coupling(name = 'GC_R2_1026',
+                 value = '(I3x32)*(R2MixedFactor*2.0)',
+                 order = {'QED':1,'QCD':2})
+
+# ============================== #
+# Goldstone UV CT couplings      #
+# ============================== #
+
+GC_UV_1015 = Coupling(name = 'GC_UV_1015',
+                 value = 'I1x33*UV_yuk_b',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1021 = Coupling(name = 'GC_UV_1021',
+                 value = '-I2x33*UV_yuk_t',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1027 = Coupling(name = 'GC_UV_1027',
+                 value = 'I3x33*UV_yuk_t',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1030 = Coupling(name = 'GC_UV_1030',
+                 value = '-I4x33*UV_yuk_b',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1082 = Coupling(name = 'GC_UV_1082',
+                 value = '-(yb/cmath.sqrt(2))*(UV_yuk_b)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1085 = Coupling(name = 'GC_UV_1085',
+                 value = '(yc/cmath.sqrt(2))*(UV_yuk_c)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1095 = Coupling(name = 'GC_UV_1095',
+                 value = '(yt/cmath.sqrt(2))*(UV_yuk_t)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1018 = Coupling(name = 'GC_UV_1018',
+                 value = '(-I2x22)*(UV_yuk_c)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1019 = Coupling(name = 'GC_UV_1019',
+                 value = '(-I2x23)*(UV_yuk_t)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1013 = Coupling(name = 'GC_UV_1013',
+                 value = '(I1x31)*(UV_yuk_b)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1014 = Coupling(name = 'GC_UV_1014',
+                 value = '(I1x32)*(UV_yuk_b)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1020 = Coupling(name = 'GC_UV_1020',
+                 value = '(-I2x32)*(UV_yuk_c)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1028 = Coupling(name = 'GC_UV_1028',
+                 value = '(-I4x13)*(UV_yuk_b)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1022 = Coupling(name = 'GC_UV_1022',
+                 value = '(I3x21)*(UV_yuk_c)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1023 = Coupling(name = 'GC_UV_1023',
+                 value = '(I3x22)*(UV_yuk_c)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1024 = Coupling(name = 'GC_UV_1024',
+                 value = '(I3x23)*(UV_yuk_c)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1029 = Coupling(name = 'GC_UV_1029',
+                 value = '(-I4x23)*(UV_yuk_b)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1025 = Coupling(name = 'GC_UV_1025',
+                 value = '(I3x31)*(UV_yuk_t)',
+                 order = {'QED':1,'QCD':2})
+
+GC_UV_1026 = Coupling(name = 'GC_UV_1026',
+                 value = '(I3x32)*(UV_yuk_t)',
+                 order = {'QED':1,'QCD':2})
diff --git a/Higgs/loop_sm_c3d4/CT_parameters.py b/Higgs/loop_sm_c3d4/CT_parameters.py
new file mode 100755
index 0000000000000000000000000000000000000000..2aec176a8a19d455dcf88ad52a13de0be9c64fea
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/CT_parameters.py
@@ -0,0 +1,167 @@
+# This file was automatically created by FeynRules $Revision: 535 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Fri 18 Mar 2011 18:40:51
+
+from object_library import all_CTparameters, CTParameter
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec
+
+################
+# R2 vertices  #
+################
+
+# ========= #
+# Pure QCD  #
+# ========= #
+
+RGR2 = CTParameter(name = 'RGR2',
+              type = 'real',
+              value = {0:'-(3.0/2.0)*G**4/(96.0*cmath.pi**2)'},
+              texname = 'RGR2')
+
+# ============== #
+# Mixed QCD-QED  #
+# ============== #
+
+R2MixedFactor = CTParameter(name = 'R2MixedFactor',
+              type = 'real',
+              value = {0:'-(G**2*(1.0+lhv)*(Ncol**2-1.0))/(2.0*Ncol*16.0*cmath.pi**2)'},
+              texname = 'R2MixedFactor')
+
+################
+# UV vertices  #
+################
+
+# ========= #
+# Pure QCD  #
+# ========= #
+
+G_UVg = CTParameter(name = 'G_UVg',
+                    type = 'real',
+                    value = {-1:'-((G**2)/(2.0*48.0*cmath.pi**2))*11.0*CA'},
+                    texname = '\delta Gg')
+
+G_UVq = CTParameter(name = 'G_UVq',
+                    type = 'real',
+                    value = {-1:'((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF'},
+                    texname = '\delta Gq')
+
+G_UVc = CTParameter(name = 'G_UVc',
+                    type = 'real',
+                    value = {-1:'((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF',
+                              0:'cond(MC,0.0,-((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF*reglog(MC**2/MU_R**2))'},
+                    texname = '\delta Gc')
+
+G_UVb = CTParameter(name = 'G_UVb',
+                    type = 'real',
+                    value = {-1:'((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF',
+                              0:'cond(MB,0.0,-((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF*reglog(MB**2/MU_R**2))'},
+                    texname = '\delta Gb')
+
+G_UVt = CTParameter(name = 'G_UVt',
+                    type = 'real',
+                    value = {-1:'((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF',
+                              0:'cond(MT,0.0,-((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF*reglog(MT**2/MU_R**2))'},
+                    texname = '\delta Gt')
+
+GWcft_UV_c = CTParameter(name = 'GWcft_UV_c',
+                         type = 'real',
+                         value = {-1:'cond(MC,0.0,-((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF)',
+                                   0:'cond(MC,0.0,((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF*reglog(MC**2/MU_R**2))'
+                                 },
+                         texname = '\delta G_{wfct\_c}')
+
+GWcft_UV_b = CTParameter(name = 'GWcft_UV_b',
+                         type = 'real',
+                         value = {-1:'cond(MB,0.0,-((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF)',
+                                   0:'cond(MB,0.0,((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF*reglog(MB**2/MU_R**2))'
+                                 },
+                         texname = '\delta G_{wfct\_b}')
+
+GWcft_UV_t = CTParameter(name = 'GWcft_UV_t',
+                         type = 'real',
+                         value = {-1:'cond(MT,0.0,-((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF)',
+                                   0:'cond(MT,0.0,((G**2)/(2.0*48.0*cmath.pi**2))*4.0*TF*reglog(MT**2/MU_R**2))'
+                                 },
+                         texname = '\delta G_{wfct\_t}')
+
+cWcft_UV = CTParameter(name = 'cWcft_UV',
+                       type = 'real',
+                       value = {-1:'cond(MC,0.0,-((G**2)/(2.0*16.0*cmath.pi**2))*3.0*CF)',
+                                 0:'cond(MC,0.0,-((G**2)/(2.0*16.0*cmath.pi**2))*CF*(4.0-3.0*reglog(MC**2/MU_R**2)))'
+                               },
+                       texname = '\delta Z_c')
+
+bWcft_UV = CTParameter(name = 'bWcft_UV',
+                       type = 'real',
+                       value = {-1:'cond(MB,0.0,-((G**2)/(2.0*16.0*cmath.pi**2))*3.0*CF)',
+                                 0:'cond(MB,0.0,-((G**2)/(2.0*16.0*cmath.pi**2))*CF*(4.0-3.0*reglog(MB**2/MU_R**2)))'
+                               },
+                       texname = '\delta Z_b')
+
+tWcft_UV = CTParameter(name = 'tWcft_UV',
+                       type = 'real',
+                       value = {-1:'cond(MT,0.0,-((G**2)/(2.0*16.0*cmath.pi**2))*3.0*CF)',
+                                 0:'cond(MT,0.0,-((G**2)/(2.0*16.0*cmath.pi**2))*CF*(4.0-3.0*reglog(MT**2/MU_R**2)))'
+                               },
+                       texname = '\delta Z_t')
+
+bMass_UV = CTParameter(name = 'bMass_UV',
+                       type = 'complex',
+                       value = {-1:'cond(MB,0.0,complex(0,1)*((G**2)/(16.0*cmath.pi**2))*(3.0*CF)*MB)',
+                                 0:'cond(MB,0.0,complex(0,1)*((G**2)/(16.0*cmath.pi**2))*CF*(4.0-3.0*reglog(MB**2/MU_R**2))*MB)'
+                               },
+                       texname = '\delta m_b')
+
+cMass_UV = CTParameter(name = 'cMass_UV',
+                       type = 'complex',
+                       value = {-1:'cond(MC,0.0,complex(0,1)*((G**2)/(16.0*cmath.pi**2))*(3.0*CF)*MC)',
+                                 0:'cond(MC,0.0,complex(0,1)*((G**2)/(16.0*cmath.pi**2))*CF*(4.0-3.0*reglog(MC**2/MU_R**2))*MC)'
+                               },
+                       texname = '\delta m_c')
+
+tMass_UV = CTParameter(name = 'tMass_UV',
+                       type = 'complex',
+                       value = {-1:'cond(MT,0.0,complex(0,1)*((G**2)/(16.0*cmath.pi**2))*3.0*CF*MT)',
+                                 0:'cond(MT,0.0,complex(0,1)*((G**2)/(16.0*cmath.pi**2))*CF*(4.0-3.0*reglog(MT**2/MU_R**2))*MT)'
+                               },
+                       texname = '\delta m_t')
+
+# ============== #
+# Mixed QCD-QED  #
+# ============== #
+
+UV_yuk_c = CTParameter(name = 'UV_yuk_c',
+                       type = 'real',
+                       value = {-1:'-(1.0/2.0)*((G**2)/(16.0*cmath.pi**2))*3.0*CF*2.0',
+                                 0:'cond(MC,0.0,-(1.0/2.0)*((G**2)/(16.0*cmath.pi**2))*CF*(-3.0*reglog(MC**2/MU_R**2)+4.0)*2.0)'
+                               },
+                       texname = '\delta y_c')
+
+UV_yuk_b = CTParameter(name = 'UV_yuk_b',
+                       type = 'real',
+                       value = {-1:'-(1.0/2.0)*((G**2)/(16.0*cmath.pi**2))*3.0*CF*2.0',
+                                 0:'cond(MB,0.0,-(1.0/2.0)*((G**2)/(16.0*cmath.pi**2))*CF*(-3.0*reglog(MB**2/MU_R**2)+4.0)*2.0)'
+                               },
+                       texname = '\delta y_b')
+
+UV_yuk_t = CTParameter(name = 'UV_yuk_t',
+                       type = 'real',
+                       value = {-1:'-(1.0/2.0)*((G**2)/(16.0*cmath.pi**2))*3.0*CF*2.0',
+                                 0:'cond(MT,0.0,-(1.0/2.0)*((G**2)/(16.0*cmath.pi**2))*CF*(-3.0*reglog(MT**2/MU_R**2)+4.0)*2.0)'
+                               },
+                       texname = '\delta y_t')
+
+# ============================== #
+# Goldstone UV CT parameters     #
+# ============================== #
+
+R2_GGGpGm_factor = CTParameter(name = 'R2_GGGpGm_factor',
+              type = 'real',
+              value = {0:'G**2/(8.0*cmath.pi**2)'},
+              texname = 'R2_GGGpGm_factor')
+
+R2_GGG0G0_factor = CTParameter(name = 'R2_GGG0G0_factor',
+              type = 'real',
+              value = {0:'G**2/(8.0*cmath.pi**2)'},
+              texname = 'R2_GGG0G0_factor')
diff --git a/Higgs/loop_sm_c3d4/CT_vertices.py b/Higgs/loop_sm_c3d4/CT_vertices.py
new file mode 100755
index 0000000000000000000000000000000000000000..a82bade01ea88aed38f082c3194b7a2b25e7cae1
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/CT_vertices.py
@@ -0,0 +1,967 @@
+# This file was automatically created by FeynRules $Revision: 535 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Fri 18 Mar 2011 18:40:51
+
+
+from object_library import all_vertices, all_CTvertices, Vertex, CTVertex
+import particles as P
+import CT_couplings as C
+import lorentz as L
+
+################
+# R2 vertices  #
+################
+
+# ========= #
+# Pure QCD  #
+# ========= #
+
+# ggg R2
+V_R23G = CTVertex(name = 'V_R23G',
+              particles = [ P.G, P.G, P.G ],
+              color = [ 'f(1,2,3)' ],
+              lorentz = [ L.VVV1 ],
+              loop_particles = [ [[P.u], [P.d], [P.c], [P.s], [P.b], [P.t]],
+                               [[P.G]] ],
+              couplings = {(0,0,0):C.R2_3Gq, (0,0,1):C.R2_3Gg},
+              type = 'R2' )
+
+#=============================================================================================
+#  4-gluon R2 vertex
+#=============================================================================================
+
+# The CT Vertex below is the one written by hand from the expression in the QCD R2 paper
+# This implementation seems to yield correct result for g g > g g g but not for g g > g g g g
+# or anytime one of the outter gluon of the vertex is offshell.
+
+# Keep in mind that Delta8(a,b) is 2 Tr(a,b)
+
+#V_R24G = CTVertex(name = 'V_R24G',
+#              particles = [ P.G, P.G, P.G,  P.G ],
+#              color = [ 'Tr(1,2)*Tr(3,4)' , 'Tr(1,3)*Tr(2,4)' , 'Tr(1,4)*Tr(2,3)', \
+#                        'd(-1,1,2)*d(-1,3,4)' , 'd(-1,1,3)*d(-1,2,4)' , 'd(-1,1,4)*d(-1,2,3)'],
+#              lorentz = [  L.R2_4G_1234, L.R2_4G_1324, L.R2_4G_1423 ],
+#              loop_particles = [ [[P.G]], [[P.u],[P.d],[P.c],[P.s],[P.b],[P.t]] ],
+#              couplings = {(0,0,0):C.GC_4GR2_Gluon_delta5,(0,1,0):C.GC_4GR2_Gluon_delta7,(0,2,0):C.GC_4GR2_Gluon_delta7, \
+#                           (1,0,0):C.GC_4GR2_Gluon_delta7,(1,1,0):C.GC_4GR2_Gluon_delta5,(1,2,0):C.GC_4GR2_Gluon_delta7, \
+#                           (2,0,0):C.GC_4GR2_Gluon_delta7,(2,1,0):C.GC_4GR2_Gluon_delta7,(2,2,0):C.GC_4GR2_Gluon_delta5, \
+#                           (3,0,0):C.GC_4GR2_4Struct,(3,1,0):C.GC_4GR2_2Struct,(3,2,0):C.GC_4GR2_2Struct, \
+#                           (4,0,0):C.GC_4GR2_2Struct,(4,1,0):C.GC_4GR2_4Struct,(4,2,0):C.GC_4GR2_2Struct, \
+#                           (5,0,0):C.GC_4GR2_2Struct,(5,1,0):C.GC_4GR2_2Struct,(5,2,0):C.GC_4GR2_4Struct , \
+#                           (0,0,1):C.GC_4GR2_Fermion_delta11,(0,1,1):C.GC_4GR2_Fermion_delta5,(0,2,1):C.GC_4GR2_Fermion_delta5, \
+#                           (1,0,1):C.GC_4GR2_Fermion_delta5,(1,1,1):C.GC_4GR2_Fermion_delta11,(1,2,1):C.GC_4GR2_Fermion_delta5, \
+#                           (2,0,1):C.GC_4GR2_Fermion_delta5,(2,1,1):C.GC_4GR2_Fermion_delta5,(2,2,1):C.GC_4GR2_Fermion_delta11, \
+#                           (3,0,1):C.GC_4GR2_11Struct,(3,1,1):C.GC_4GR2_5Struct,(3,2,1):C.GC_4GR2_5Struct, \
+#                           (4,0,1):C.GC_4GR2_5Struct,(4,1,1):C.GC_4GR2_11Struct,(4,2,1):C.GC_4GR2_5Struct, \
+#                           (5,0,1):C.GC_4GR2_5Struct,(5,1,1):C.GC_4GR2_5Struct,(5,2,1):C.GC_4GR2_11Struct },
+#              type = 'R2')
+
+# The CT Vertex below is the one written automatically by FR
+# Gives the same result as above for g g > g g but not as soon as one of the outter gluon is offshell.
+
+V_R2RGA = CTVertex(name = 'V_R2RGA',
+               type = 'R2',
+               particles = [ P.G, P.G, P.G, P.G ],
+               color = [ 'd(-1,1,3)*d(-1,2,4)', 'd(-1,1,3)*f(-1,2,4)', 'd(-1,1,4)*d(-1,2,3)', 'd(-1,1,4)*f(-1,2,3)', 'd(-1,2,3)*f(-1,1,4)', 'd(-1,2,4)*f(-1,1,3)', 'f(-1,1,2)*f(-1,3,4)', 'f(-1,1,3)*f(-1,2,4)', 'f(-1,1,4)*f(-1,2,3)', 'Identity(1,2)*Identity(3,4)', 'Identity(1,3)*Identity(2,4)', 'Identity(1,4)*Identity(2,3)' ],
+               lorentz = [ L.R2RGA_VVVV10, L.R2RGA_VVVV2, L.R2RGA_VVVV3, L.R2RGA_VVVV5 ],
+               loop_particles = [ [ [P.b], [P.c], [P.d], [P.s], [P.t], [P.u] ], [ [P.G] ] ],
+               couplings = {(2,1,0):C.R2GC_137_43,(2,1,1):C.R2GC_137_44,(0,1,0):C.R2GC_137_43,(0,1,1):C.R2GC_137_44,(4,1,0):C.R2GC_145_58,(4,1,1):C.R2GC_145_59,(3,1,0):C.R2GC_145_58,(3,1,1):C.R2GC_145_59,(8,1,0):C.R2GC_138_45,(8,1,1):C.R2GC_138_46,(7,1,0):C.R2GC_144_56,(7,1,1):C.R2GC_144_57,(6,1,0):C.R2GC_141_50,(6,1,1):C.R2GC_141_51,(5,1,0):C.R2GC_145_58,(5,1,1):C.R2GC_145_59,(1,1,0):C.R2GC_145_58,(1,1,1):C.R2GC_145_59,(11,0,0):C.R2GC_140_48,(11,0,1):C.R2GC_140_49,(10,0,0):C.R2GC_140_48,(10,0,1):C.R2GC_140_49,(9,0,1):C.R2GC_139_47,(2,2,0):C.R2GC_137_43,(2,2,1):C.R2GC_137_44,(0,2,0):C.R2GC_137_43,(0,2,1):C.R2GC_137_44,(6,2,0):C.R2GC_142_52,(6,2,1):C.R2GC_142_53,(4,2,0):C.R2GC_145_58,(4,2,1):C.R2GC_145_59,(3,2,0):C.R2GC_145_58,(3,2,1):C.R2GC_145_59,(8,2,0):C.R2GC_144_56,(8,2,1):C.R2GC_144_57,(5,2,0):C.R2GC_145_58,(5,2,1):C.R2GC_145_59,(1,2,0):C.R2GC_145_58,(1,2,1):C.R2GC_145_59,(7,2,0):C.R2GC_138_45,(7,2,1):C.R2GC_138_46,(2,3,0):C.R2GC_137_43,(2,3,1):C.R2GC_137_44,(0,3,0):C.R2GC_137_43,(0,3,1):C.R2GC_137_44,(4,3,0):C.R2GC_145_58,(4,3,1):C.R2GC_145_59,(3,3,0):C.R2GC_145_58,(3,3,1):C.R2GC_145_59,(8,3,0):C.R2GC_143_54,(8,3,1):C.R2GC_143_55,(7,3,0):C.R2GC_143_54,(7,3,1):C.R2GC_143_55,(5,3,0):C.R2GC_145_58,(5,3,1):C.R2GC_145_59,(1,3,0):C.R2GC_145_58,(1,3,1):C.R2GC_145_59})
+
+#=============================================================================================
+
+# gdd~
+V_R2GDD = CTVertex(name = 'V_R2GDD',
+              particles = [ P.d__tilde__, P.d, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles =[[[P.d,P.G]]],                 
+              couplings = {(0,0,0):C.R2_GQQ},
+              type = 'R2')
+
+# guu~              
+V_R2GUU = CTVertex(name = 'V_R2GUU',
+               particles = [ P.u__tilde__, P.u, P.G ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               loop_particles =[[[P.u,P.G]]],
+               couplings = {(0,0,0):C.R2_GQQ},
+               type = 'R2')  
+
+# gss~
+V_R2GSS = CTVertex(name = 'V_R2GSS',
+              particles = [ P.s__tilde__, P.s, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles =[[[P.s,P.G]]],
+              couplings = {(0,0,0):C.R2_GQQ},
+              type = 'R2')
+
+# gcc~              
+V_R2GCC = CTVertex(name = 'V_R2GCC',
+               particles = [ P.c__tilde__, P.c, P.G ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               loop_particles =[[[P.c,P.G]]],               
+               couplings = {(0,0,0):C.R2_GQQ},
+               type = 'R2')  
+
+# gbb~
+V_R2GBB = CTVertex(name = 'V_R2GBB',
+              particles = [ P.b__tilde__, P.b, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles =[[[P.b,P.G]]],
+              couplings = {(0,0,0):C.R2_GQQ},
+              type = 'R2')
+
+# gtt~              
+V_R2GTT = CTVertex(name = 'V_R2GTT',
+               particles = [ P.t__tilde__, P.t, P.G ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               loop_particles =[[[P.t,P.G]]],               
+               couplings = {(0,0,0):C.R2_GQQ},
+               type = 'R2')
+
+# gg             
+V_R2GG = CTVertex(name = 'V_R2GG',
+               particles = [ P.G, P.G ],
+               color = [ 'Tr(1,2)' ],
+               lorentz = [ L.R2_GG_1, L.R2_GG_2, L.R2_GG_3],
+               loop_particles = [ [[P.u],[P.d],[P.s]],
+                                  [[P.c]],
+                                  [[P.b]],
+                                  [[P.t]],
+                                  [[P.G]] ],
+               couplings = {(0,0,0):C.R2_GGq,
+                            (0,0,1):C.R2_GGq,(0,2,1):C.R2_GGc,                            
+                            (0,0,2):C.R2_GGq,(0,2,2):C.R2_GGb,
+                            (0,0,3):C.R2_GGq,(0,2,3):C.R2_GGt,
+                            (0,0,4):C.R2_GGg_1, (0,1,4):C.R2_GGg_2},
+               type = 'R2')
+
+# d~d            
+V_R2DD = CTVertex(name = 'V_R2DD',
+               particles = [ P.d__tilde__, P.d ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.R2_QQ_1 ],
+               loop_particles = [[[P.d,P.G]]],
+               couplings = {(0,0,0):C.R2_QQq},
+               type = 'R2') 
+
+# u~u            
+V_R2UU = CTVertex(name = 'V_R2UU',
+               particles = [ P.u__tilde__, P.u ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.R2_QQ_1 ],
+               loop_particles = [[[P.u,P.G]]],            
+               couplings = {(0,0,0):C.R2_QQq},
+               type = 'R2')
+
+# s~s            
+V_R2SS = CTVertex(name = 'V_R2SS',
+               particles = [ P.s__tilde__, P.s ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.R2_QQ_1 ],
+               loop_particles = [[[P.s,P.G]]],                
+               couplings = {(0,0,0):C.R2_QQq},
+               type = 'R2')
+
+# c~c            
+V_R2CC = CTVertex(name = 'V_R2CC',
+               particles = [ P.c__tilde__, P.c ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.R2_QQ_1, L.R2_QQ_2 ],
+               loop_particles = [[[P.c,P.G]]],
+               couplings = {(0,0,0):C.R2_QQq,(0,1,0):C.R2_QQc},                
+               type = 'R2')
+
+# b~b            
+V_R2BB = CTVertex(name = 'V_R2BB',
+               particles = [ P.b__tilde__, P.b ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.R2_QQ_1, L.R2_QQ_2 ],
+               loop_particles = [[[P.b,P.G]]],
+               couplings = {(0,0,0):C.R2_QQq,(0,1,0):C.R2_QQb},                
+               type = 'R2')
+
+# t~t            
+V_R2TT = CTVertex(name = 'V_R2TT',
+               particles = [ P.t__tilde__, P.t ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.R2_QQ_1, L.R2_QQ_2 ],
+               loop_particles = [[[P.t,P.G]]],
+               couplings = {(0,0,0):C.R2_QQq,(0,1,0):C.R2_QQt},
+               type = 'R2')
+
+# ============== #
+# Mixed QCD-QED  #
+# ============== #
+
+# R2 for the A and Z couplings to the quarks
+
+V_R2ddA = CTVertex(name = 'V_R2ddA',
+              particles = [ P.d__tilde__, P.d, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.d,P.G]]],
+              couplings = {(0,0,0):C.R2_DDA},
+              type = 'R2')
+
+V_R2ssA = CTVertex(name = 'V_R2ssA',
+              particles = [ P.s__tilde__, P.s, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.s,P.G]]],
+              couplings = {(0,0,0):C.R2_DDA},
+              type = 'R2')
+
+V_R2bbA = CTVertex(name = 'V_R2bbA',
+              particles = [ P.b__tilde__, P.b, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.b,P.G]]],
+              couplings = {(0,0,0):C.R2_DDA},
+              type = 'R2')
+
+V_R2uuA = CTVertex(name = 'V_R2uuA',
+              particles = [ P.u__tilde__, P.u, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.u,P.G]]],
+              couplings = {(0,0,0):C.R2_UUA},
+              type = 'R2')
+
+V_R2ccA = CTVertex(name = 'V_R2ccA',
+              particles = [ P.c__tilde__, P.c, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.c,P.G]]],
+              couplings = {(0,0,0):C.R2_UUA},
+              type = 'R2')
+
+V_R2ttA = CTVertex(name = 'V_R2ttA',
+              particles = [ P.t__tilde__, P.t, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.t,P.G]]],
+              couplings = {(0,0,0):C.R2_UUA},
+              type = 'R2')
+
+V_R2ddZ = CTVertex(name = 'V_R2ddZ',
+              particles = [ P.d__tilde__, P.d, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV3 ],
+              loop_particles = [[[P.d,P.G]]],
+              couplings = {(0,0,0):C.R2_DDZ_V2,(0,1,0):C.R2_DDZ_V3},
+              type = 'R2')
+
+V_R2ssZ = CTVertex(name = 'V_R2ssZ',
+              particles = [ P.s__tilde__, P.s, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV3 ],
+              loop_particles = [[[P.s,P.G]]],
+              couplings = {(0,0,0):C.R2_DDZ_V2,(0,1,0):C.R2_DDZ_V3},
+              type = 'R2')
+
+V_R2bbZ = CTVertex(name = 'V_R2bbZ',
+              particles = [ P.b__tilde__, P.b, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV3 ],
+              loop_particles = [[[P.b,P.G]]],
+              couplings = {(0,0,0):C.R2_DDZ_V2,(0,1,0):C.R2_DDZ_V3},
+              type = 'R2')
+
+V_R2uuZ = CTVertex(name = 'V_R2uuZ',
+              particles = [ P.u__tilde__, P.u, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV5 ],
+              loop_particles = [[[P.u,P.G]]],
+              couplings = {(0,0,0):C.R2_UUZ_V2,(0,1,0):C.R2_UUZ_V5},
+              type = 'R2')
+
+V_R2ccZ = CTVertex(name = 'V_R2ccZ',
+              particles = [ P.c__tilde__, P.c, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV5 ],
+              loop_particles = [[[P.c,P.G]]],
+              couplings = {(0,0,0):C.R2_UUZ_V2,(0,1,0):C.R2_UUZ_V5},
+              type = 'R2')
+
+V_R2ttZ = CTVertex(name = 'V_R2ttZ',
+              particles = [ P.t__tilde__, P.t, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV5 ],
+              loop_particles = [[[P.t,P.G]]],
+              couplings = {(0,0,0):C.R2_UUZ_V2,(0,1,0):C.R2_UUZ_V5},
+              type = 'R2')
+
+# R2 for the W couplings to the quarks with most general CKM
+
+V_R2dxuW = CTVertex(name = 'V_R2dxuW',
+              particles = [ P.d__tilde__, P.u, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.d,P.u,P.G]]],                   
+              couplings = {(0,0,0):C.R2_dxuW},
+              type = 'R2')
+
+V_R2dxcW = CTVertex(name = 'V_R2dxcW',
+              particles = [ P.d__tilde__, P.c, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.d,P.c,P.G]]],                   
+              couplings = {(0,0,0):C.R2_dxcW},
+              type = 'R2')
+
+V_R2dxtW = CTVertex(name = 'V_R2dxtW',
+              particles = [ P.d__tilde__, P.t, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.d,P.t,P.G]]],                   
+              couplings = {(0,0,0):C.R2_dxtW},
+              type = 'R2')
+
+V_R2sxuW = CTVertex(name = 'V_R2sxuW',
+              particles = [ P.s__tilde__, P.u, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.s,P.u,P.G]]],                   
+              couplings = {(0,0,0):C.R2_sxuW},
+              type = 'R2')
+
+V_R2sxcW = CTVertex(name = 'V_R2sxcW',
+              particles = [ P.s__tilde__, P.c, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.s,P.c,P.G]]],                   
+              couplings = {(0,0,0):C.R2_sxcW},
+              type = 'R2')
+
+V_R2sxtW = CTVertex(name = 'V_R2sxtW',
+              particles = [ P.s__tilde__, P.t, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.s,P.t,P.G]]],                   
+              couplings = {(0,0,0):C.R2_sxtW},
+              type = 'R2')
+
+V_R2bxuW = CTVertex(name = 'V_R2bxuW',
+              particles = [ P.b__tilde__, P.u, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.b,P.u,P.G]]],                   
+              couplings = {(0,0,0):C.R2_bxuW},
+              type = 'R2')
+
+V_R2bxcW = CTVertex(name = 'V_R2bxcW',
+              particles = [ P.b__tilde__, P.c, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.b,P.c,P.G]]],                   
+              couplings = {(0,0,0):C.R2_bxcW},
+              type = 'R2')
+
+V_R2bxtW = CTVertex(name = 'V_R2bxtW',
+              particles = [ P.b__tilde__, P.t, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.b,P.t,P.G]]],                   
+              couplings = {(0,0,0):C.R2_bxtW},
+              type = 'R2')
+
+V_R2uxdW = CTVertex(name = 'V_R2uxdW',
+              particles = [ P.u__tilde__, P.d, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.u,P.d,P.G]]],                   
+              couplings = {(0,0,0):C.R2_uxdW},
+              type = 'R2')
+
+V_R2cxdW = CTVertex(name = 'V_R2cxdW',
+              particles = [ P.c__tilde__, P.d, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.c,P.d,P.G]]],                   
+              couplings = {(0,0,0):C.R2_cxdW},
+              type = 'R2')
+
+V_R2txdW = CTVertex(name = 'V_R2txdW',
+              particles = [ P.t__tilde__, P.d, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.t,P.d,P.G]]],                   
+              couplings = {(0,0,0):C.R2_txdW},
+              type = 'R2')
+
+V_R2uxsW = CTVertex(name = 'V_R2uxsW',
+              particles = [ P.u__tilde__, P.s, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.u,P.s,P.G]]],                   
+              couplings = {(0,0,0):C.R2_uxsW},
+              type = 'R2')
+
+V_R2cxsW = CTVertex(name = 'V_R2cxsW',
+              particles = [ P.c__tilde__, P.s, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.c,P.s,P.G]]],                   
+              couplings = {(0,0,0):C.R2_cxsW},
+              type = 'R2')
+
+V_R2txsW = CTVertex(name = 'V_R2txsW',
+              particles = [ P.t__tilde__, P.s, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.t,P.s,P.G]]],                   
+              couplings = {(0,0,0):C.R2_txsW},
+              type = 'R2')
+
+V_R2uxbW = CTVertex(name = 'V_R2uxbW',
+              particles = [ P.u__tilde__, P.b, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.u,P.b,P.G]]],                   
+              couplings = {(0,0,0):C.R2_uxbW},
+              type = 'R2')
+
+V_R2cxbW = CTVertex(name = 'V_R2cxbW',
+              particles = [ P.c__tilde__, P.b, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.c,P.b,P.G]]],                   
+              couplings = {(0,0,0):C.R2_cxbW},
+              type = 'R2')
+
+V_R2txbW = CTVertex(name = 'V_R2txbW',
+              particles = [ P.t__tilde__, P.b, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              loop_particles = [[[P.t,P.b,P.G]]],                   
+              couplings = {(0,0,0):C.R2_txbW},
+              type = 'R2')
+
+# R2 for the Higgs couplings to massive quarks 
+
+V_bbH = CTVertex(name = 'V_bbH',
+              particles = [ P.b__tilde__, P.b, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1 ],
+              loop_particles = [[[P.b,P.G]]],
+              couplings = {(0,0,0):C.R2_bbH},
+              type = 'R2')
+
+V_ttH = CTVertex(name = 'V_ttH',
+              particles = [ P.t__tilde__, P.t, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1 ],
+              loop_particles = [[[P.t,P.G]]],
+              couplings = {(0,0,0):C.R2_ttH},
+              type = 'R2')
+
+V_ccH = CTVertex(name = 'V_ccH',
+              particles = [ P.c__tilde__, P.c, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1 ],
+              loop_particles = [[[P.c,P.G]]],
+              couplings = {(0,0,0):C.R2_ccH},
+              type = 'R2')
+
+# R2 interactions non proportional to the SM
+
+# R2 for the Higgs interactions
+
+V_GGH = CTVertex(name = 'V_GGH',
+              particles = [ P.G, P.G, P.H ],
+              color = [ 'Tr(1,2)' ],
+              lorentz = [ L.VVS1 ],
+              loop_particles = [[[P.c]],[[P.b]],[[P.t]]],
+              couplings = {(0,0,0):C.R2_GGHc,(0,0,1):C.R2_GGHb,(0,0,2):C.R2_GGHt},
+              type = 'R2')
+
+V_GGHH = CTVertex(name = 'V_GGHH',
+              particles = [ P.G, P.G, P.H, P.H ],
+              color = [ 'Tr(1,2)' ],
+              lorentz = [ L.R2_GGHH ],
+              loop_particles = [[[P.c]],[[P.b]],[[P.t]]],
+              couplings = {(0,0,0):C.R2_GGHHc,(0,0,1):C.R2_GGHHb,(0,0,2):C.R2_GGHHt},
+              type = 'R2')
+
+# R2 for the weak vector bosons interaction with gluons
+
+V_GGZ = CTVertex(name = 'V_GGZ',
+              particles = [ P.G, P.G, P.Z ],
+              color = [ 'Tr(1,2)' ],
+              lorentz = [ L.R2_GGZ ],
+              loop_particles = [[[P.u],[P.c],[P.t]],[[P.d],[P.s],[P.b]]],
+              couplings = {(0,0,0):C.R2_GGZup,(0,0,1):C.R2_GGZdown},
+              type = 'R2')
+
+
+V_GGZZ = CTVertex(name = 'V_GGZZ',
+              particles = [ P.G, P.G, P.Z, P.Z ],
+              color = [ 'Tr(1,2)' ],
+              lorentz = [ L.R2_GGVV ],
+              loop_particles = [[[P.u],[P.c],[P.t]],[[P.d],[P.s],[P.b]]],
+              couplings = {(0,0,0):C.R2_GGZZup,(0,0,1):C.R2_GGZZdown},
+              type = 'R2')
+
+V_GGAA = CTVertex(name = 'V_GGAA',
+              particles = [ P.G, P.G, P.A, P.A ],
+              color = [ 'Tr(1,2)' ],
+              lorentz = [ L.R2_GGVV ],
+              loop_particles = [[[P.u],[P.c],[P.t]],[[P.d],[P.s],[P.b]]],
+              couplings = {(0,0,0):C.R2_GGAAup,(0,0,1):C.R2_GGAAdown},
+              type = 'R2')
+
+V_GGZA = CTVertex(name = 'V_GGZA',
+              particles = [ P.G, P.G, P.Z, P.A ],
+              color = [ 'Tr(1,2)' ],
+              lorentz = [ L.R2_GGVV ],
+              loop_particles = [[[P.u],[P.c],[P.t]],[[P.d],[P.s],[P.b]]],
+              couplings = {(0,0,0):C.R2_GGZAup,(0,0,1):C.R2_GGZAdown},
+              type = 'R2')
+
+V_GGWW = CTVertex(name = 'V_GGWW',
+              particles = [ P.G, P.G, P.W__minus__, P.W__plus__ ],
+              color = [ 'Tr(1,2)' ],
+              lorentz = [ L.R2_GGVV ],
+              loop_particles = [[[P.u,P.d]],[[P.u,P.s]],[[P.u,P.b]],
+                                [[P.c,P.d]],[[P.c,P.s]],[[P.c,P.b]],
+                                [[P.t,P.d]],[[P.t,P.s]],[[P.t,P.b]]],
+              couplings = {(0,0,0):C.R2_GGWWud,(0,0,1):C.R2_GGWWus,(0,0,2):C.R2_GGWWub,
+                           (0,0,3):C.R2_GGWWcd,(0,0,4):C.R2_GGWWcs,(0,0,5):C.R2_GGWWcb,
+                           (0,0,6):C.R2_GGWWtd,(0,0,7):C.R2_GGWWts,(0,0,8):C.R2_GGWWtb},
+              type = 'R2')
+
+V_GGGZ = CTVertex(name = 'V_GGGZ',
+              particles = [ P.G, P.G, P.G, P.Z ],
+              color = [ 'd(1,2,3)' , 'f(1,2,3)'],
+              lorentz = [ L.R2_GGVV, L.R2_GGGVa ],
+              loop_particles = [[[P.u],[P.c],[P.t]],[[P.d],[P.s],[P.b]]],
+              couplings = {(0,0,0):C.R2_GGGZvecUp,(0,0,1):C.R2_GGGZvecDown,
+                           (1,1,0):C.R2_GGGZaxialUp,(1,1,1):C.R2_GGGZaxialDown},
+              type = 'R2')
+
+V_GGGA = CTVertex(name = 'V_GGGA',
+              particles = [ P.G, P.G, P.G, P.A ],
+              color = [ 'd(1,2,3)'],
+              lorentz = [ L.R2_GGVV ],
+              loop_particles = [[[P.u],[P.c],[P.t]],[[P.d],[P.s],[P.b]]],
+              couplings = {(0,0,0):C.R2_GGGAvecUp,(0,0,1):C.R2_GGGAvecDown},
+              type = 'R2')
+
+################
+# UV vertices  #
+################
+
+# ========= #
+# Pure QCD  #
+# ========= #
+
+# There are the alpha_s renormalization vertices
+
+# ggg
+V_UV1eps3G = CTVertex(name = 'V_UV1eps3G',
+              particles = [ P.G, P.G, P.G ],
+              color = [ 'f(1,2,3)' ],
+              lorentz = [ L.VVV1 ],
+              loop_particles = [[[P.u],[P.d],[P.s]],[[P.c]],[[P.b]],[[P.t]],[[P.G]]],
+              couplings = {(0,0,0):C.UV_3Gq,(0,0,1):C.UV_3Gc,(0,0,2):C.UV_3Gb,(0,0,3):C.UV_3Gt,(0,0,4):C.UV_3Gg},
+              type = 'UV')
+
+# gggg
+V_UV4G = CTVertex(name = 'V_UV1eps4G',
+              particles = [ P.G, P.G, P.G, P.G ],
+              color = [ 'f(-1,1,2)*f(3,4,-1)', 'f(-1,1,3)*f(2,4,-1)', 'f(-1,1,4)*f(2,3,-1)' ],
+              lorentz = [ L.VVVV1, L.VVVV3, L.VVVV4 ],
+              loop_particles = [[[P.u],[P.d],[P.s]],[[P.c]],[[P.b]],[[P.t]],[[P.G]]],
+              couplings = {(0,0,0):C.UV_4Gq,(0,0,1):C.UV_4Gc,(0,0,2):C.UV_4Gb,(0,0,3):C.UV_4Gt,(0,0,4):C.UV_4Gg,
+                           (1,1,0):C.UV_4Gq,(1,1,1):C.UV_4Gc,(1,1,2):C.UV_4Gb,(1,1,3):C.UV_4Gt,(1,1,4):C.UV_4Gg,
+                           (2,2,0):C.UV_4Gq,(2,2,1):C.UV_4Gc,(2,2,2):C.UV_4Gb,(2,2,3):C.UV_4Gt,(2,2,4):C.UV_4Gg},
+              type = 'UV')
+
+# gdd~
+V_UVGDD = CTVertex(name = 'V_UVGDD',
+              particles = [ P.d__tilde__, P.d, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.u],[P.d],[P.s]],[[P.c]],[[P.b]],[[P.t]],[[P.G]]],
+              couplings = {(0,0,0):C.UV_GQQq,(0,0,1):C.UV_GQQc,(0,0,2):C.UV_GQQb,(0,0,3):C.UV_GQQt,(0,0,4):C.UV_GQQg},
+              type = 'UV')
+
+# guu~
+V_UVGUU = CTVertex(name = 'V_UVGUU',
+              particles = [ P.u__tilde__, P.u, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.u],[P.d],[P.s]],[[P.c]],[[P.b]],[[P.t]],[[P.G]]],
+              couplings = {(0,0,0):C.UV_GQQq,(0,0,1):C.UV_GQQc,(0,0,2):C.UV_GQQb,(0,0,3):C.UV_GQQt,(0,0,4):C.UV_GQQg},
+              type = 'UV')
+
+# gcc~
+V_UVGCC = CTVertex(name = 'V_UVGCC',
+              particles = [ P.c__tilde__, P.c, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.u],[P.d],[P.s]],[[P.c]],[[P.b]],[[P.t]],[[P.G]]],
+              couplings = {(0,0,0):C.UV_GQQq,(0,0,1):C.UV_GQQc,(0,0,2):C.UV_GQQb,(0,0,3):C.UV_GQQt,(0,0,4):C.UV_GQQg},
+              type = 'UV')
+
+# gss~
+V_UVGSS = CTVertex(name = 'V_UVGSS',
+              particles = [ P.s__tilde__, P.s, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.u],[P.d],[P.s]],[[P.c]],[[P.b]],[[P.t]],[[P.G]]],
+              couplings = {(0,0,0):C.UV_GQQq,(0,0,1):C.UV_GQQc,(0,0,2):C.UV_GQQb,(0,0,3):C.UV_GQQt,(0,0,4):C.UV_GQQg},
+              type = 'UV')
+
+# gbb~
+V_UVGBB = CTVertex(name = 'V_UVGBB',
+              particles = [ P.b__tilde__, P.b, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.u],[P.d],[P.s]],[[P.c]],[[P.b]],[[P.t]],[[P.G]]],
+              couplings = {(0,0,0):C.UV_GQQq,(0,0,1):C.UV_GQQc,(0,0,2):C.UV_GQQb,(0,0,3):C.UV_GQQt,(0,0,4):C.UV_GQQg},
+              type = 'UV')
+
+# gtt~
+V_UVGTT = CTVertex(name = 'V_UVGTT',
+              particles = [ P.t__tilde__, P.t, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              loop_particles = [[[P.u],[P.d],[P.s]],[[P.c]],[[P.b]],[[P.t]],[[P.G]]],
+              couplings = {(0,0,0):C.UV_GQQq,(0,0,1):C.UV_GQQc,(0,0,2):C.UV_GQQb,(0,0,3):C.UV_GQQt,(0,0,4):C.UV_GQQg},
+              type = 'UV')
+
+# These are the mass renormalization vertices.
+
+# c~c         
+V_UVcMass = CTVertex(name = 'V_UVcMass',
+               particles = [ P.c__tilde__, P.c ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.R2_QQ_2 ],
+               loop_particles = [[[P.G,P.c]]],                   
+               couplings = {(0,0,0):C.UV_cMass},
+               type = 'UVmass')
+
+# b~b         
+V_UVbMass = CTVertex(name = 'V_UVbMass',
+               particles = [ P.b__tilde__, P.b ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.R2_QQ_2 ],
+               loop_particles = [[[P.G,P.b]]],                   
+               couplings = {(0,0,0):C.UV_bMass},
+               type = 'UVmass') 
+
+# t~t         
+V_UVtMass = CTVertex(name = 'V_UVtMass',
+               particles = [ P.t__tilde__, P.t ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.R2_QQ_2 ],
+               loop_particles = [[[P.G,P.t]]],                   
+               couplings = {(0,0,0):C.UV_tMass},
+               type = 'UVmass')
+
+# ============== #
+# Mixed QCD-QED  #
+# ============== #
+
+V_UVHcc = CTVertex(name = 'V_UVHcc',
+              particles = [ P.c__tilde__, P.c, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1 ],
+              loop_particles = [[[P.G,P.c]]],                   
+              couplings = {(0,0,0):C.UV_Hcc},
+              type = 'UV')
+
+V_UVHtt = CTVertex(name = 'V_UVHtt',
+              particles = [ P.t__tilde__, P.t, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1 ],
+              loop_particles = [[[P.G,P.t]]],                   
+              couplings = {(0,0,0):C.UV_Htt},
+              type = 'UV')
+
+V_UVHbb = CTVertex(name = 'V_UVHbb',
+              particles = [ P.b__tilde__, P.b, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1 ],
+              loop_particles = [[[P.G,P.b]]],
+              couplings = {(0,0,0):C.UV_Hbb},
+              type = 'UV')
+
+
+# ============================== #
+# Goldstone R2 counter vertices  #
+# ============================== #
+
+V_R2_GGGpGm = CTVertex(name = 'V_R2_GGGpGm',
+              particles = [ P.G, P.G, P.G__plus__, P.G__minus__],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.VVSS1 ],  
+              loop_particles = [[[P.u,P.b]],
+                                [[P.c,P.d]],[[P.c,P.s]],[[P.c,P.b]],
+                                [[P.t,P.d]],[[P.t,P.s]],[[P.t,P.b]]],
+              couplings = {(0,0,0):C.R2_GGGpGm_ub,
+                           (0,0,1):C.R2_GGGpGm_cd,(0,0,2):C.R2_GGGpGm_cs,(0,0,3):C.R2_GGGpGm_cb,
+                           (0,0,4):C.R2_GGGpGm_td,(0,0,5):C.R2_GGGpGm_ts,(0,0,6):C.R2_GGGpGm_tb},
+              type = 'R2')
+
+V_R2_GGG0G0 = CTVertex(name = 'V_R2_GGG0G0',
+              particles = [ P.G, P.G, P.G0, P.G0 ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.VVSS1 ],
+              loop_particles = [[[P.c]],[[P.b]],[[P.t]]],
+              couplings = {(0,0,0):C.R2_GGG0G0_c,(0,0,1):C.R2_GGG0G0_b,(0,0,2):C.R2_GGG0G0_t},
+              type = 'R2')
+
+V_R2_txbGp = CTVertex(name = 'V_R2_txbGp',
+              particles = [ P.t__tilde__, P.b, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS8, L.FFS3 ],                    
+              loop_particles = [[[P.G,P.b,P.t]]],
+              couplings = {(0,0,0):C.GC_R2_1015,(0,1,0):C.GC_R2_1021},               
+              type = 'R2')
+
+V_R2_bxtGp = CTVertex(name = 'V_R2_bxtGp',
+              particles = [ P.b__tilde__, P.t, P.G__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS8, L.FFS3 ],                    
+              loop_particles = [[[P.G,P.b,P.t]]],
+              couplings = {(0,0,0):C.GC_R2_1027,(0,1,0):C.GC_R2_1030},                   
+              type = 'R2')
+
+V_R2_1077 = CTVertex(name = 'V_R2_1077',
+              particles = [ P.b__tilde__, P.b, P.G0 ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS2 ],
+              loop_particles = [[[P.G,P.b]]],
+              couplings = {(0,0,0):C.GC_R2_1082},
+              type = 'R2')
+
+V_R2_1138 = CTVertex(name = 'V_R2_1138',
+               particles = [ P.c__tilde__, P.c, P.G0 ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS2 ],
+               loop_particles = [[[P.G,P.c]]],
+               couplings = {(0,0,0):C.GC_R2_1085},
+               type = 'R2')
+
+V_R2_1139 = CTVertex(name = 'V_R2_1139',
+               particles = [ P.t__tilde__, P.t, P.G0 ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS2 ],
+               loop_particles = [[[P.G,P.t]]],
+               couplings = {(0,0,0):C.GC_R2_1095},
+               type = 'R2')
+
+
+V_R2_1084 = CTVertex(name = 'V_R2_1084',
+              particles = [ P.c__tilde__, P.s, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS3 ],
+              loop_particles = [[[P.G,P.c,P.s]]],
+              couplings = {(0,0,0):C.GC_R2_1018},
+              type = 'R2')
+
+
+V_R2_1085 = CTVertex(name = 'V_R2_1085',
+              particles = [ P.t__tilde__, P.s, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS3 ],
+              loop_particles = [[[P.G,P.t,P.s]]],
+              couplings = {(0,0,0):C.GC_R2_1019},
+              type = 'R2')
+
+
+V_R2_1086 = CTVertex(name = 'V_R2_1086',
+              particles = [ P.u__tilde__, P.b, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS8 ],
+              loop_particles = [[[P.G,P.u,P.b]]],
+              couplings = {(0,0,0):C.GC_R2_1013},
+              type = 'R2')
+
+V_R2_1087 = CTVertex(name = 'V_R2_1087',
+              particles = [ P.c__tilde__, P.b, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS8, L.FFS3 ],
+              loop_particles = [[[P.G,P.c,P.b]]],
+              couplings = {(0,0,0):C.GC_R2_1014,(0,1,0):C.GC_R2_1020},
+              type = 'R2')
+
+V_R2_1116 = CTVertex(name = 'V_R2_1116',
+               particles = [ P.b__tilde__, P.u, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS3 ],
+               loop_particles = [[[P.G,P.b,P.u]]],
+               couplings = {(0,0,0):C.GC_R2_1028},
+               type = 'R2')
+
+V_R2_1117 = CTVertex(name = 'V_R2_1117',
+               particles = [ P.d__tilde__, P.c, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               loop_particles = [[[P.G,P.d,P.c]]],
+               couplings = {(0,0,0):C.GC_R2_1022},
+               type = 'R2')
+
+V_R2_1118 = CTVertex(name = 'V_R2_1118',
+               particles = [ P.s__tilde__, P.c, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               loop_particles = [[[P.G,P.c,P.s]]],
+               couplings = {(0,0,0):C.GC_R2_1023},
+               type = 'R2')
+
+V_R2_1119 = CTVertex(name = 'V_R2_1119',
+               particles = [ P.b__tilde__, P.c, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8, L.FFS3 ],
+               loop_particles = [[[P.G,P.c,P.b]]],
+               couplings = {(0,0,0):C.GC_R2_1024,(0,1,0):C.GC_R2_1029},
+               type = 'R2')
+
+V_R2_1120 = CTVertex(name = 'V_R2_1120',
+               particles = [ P.d__tilde__, P.t, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               loop_particles = [[[P.G,P.d,P.t]]],
+               couplings = {(0,0,0):C.GC_R2_1025},
+               type = 'R2')
+
+V_R2_1121 = CTVertex(name = 'V_R2_1121',
+               particles = [ P.s__tilde__, P.t, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               loop_particles = [[[P.G,P.t,P.s]]],
+               couplings = {(0,0,0):C.GC_R2_1026},
+               type = 'R2')
+
+# ============================== #
+# Goldstone UV counter vertices  #
+# ============================== #
+
+V_UV_txbGp = CTVertex(name = 'V_UV_txbGp',
+              particles = [ P.t__tilde__, P.b, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS8, L.FFS3 ],
+              loop_particles = [[[P.G,P.b,P.t]]],                   
+              couplings = {(0,0,0):C.GC_UV_1015,(0,1,0):C.GC_UV_1021},
+              type = 'UV')
+
+V_UV_bxtGp = CTVertex(name = 'V_UV_bxtGp',
+               particles = [ P.b__tilde__, P.t, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8, L.FFS3 ],
+               loop_particles = [[[P.G,P.b,P.t]]],
+               couplings = {(0,0,0):C.GC_UV_1027,(0,1,0):C.GC_UV_1030},
+               type = 'UV')
+
+V_UV_1077 = CTVertex(name = 'V_UV_1077',
+              particles = [ P.b__tilde__, P.b, P.G0 ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS2 ],
+              loop_particles = [[[P.G,P.b]]],
+              couplings = {(0,0,0):C.GC_UV_1082},
+              type = 'UV')
+
+V_UV_1138 = CTVertex(name = 'V_UV_1138',
+               particles = [ P.c__tilde__, P.c, P.G0 ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS2 ],
+               loop_particles = [[[P.G,P.c]]],
+               couplings = {(0,0,0):C.GC_UV_1085},
+               type = 'UV')
+
+V_UV_1139 = CTVertex(name = 'V_UV_1139',
+               particles = [ P.t__tilde__, P.t, P.G0 ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS2 ],
+               loop_particles = [[[P.G,P.t]]],
+               couplings = {(0,0,0):C.GC_UV_1095},
+               type = 'UV')
+
+V_UV_1084 = CTVertex(name = 'V_UV_1084',
+              particles = [ P.c__tilde__, P.s, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS3 ],
+              loop_particles = [[[P.G,P.c,P.s]]],
+              couplings = {(0,0,0):C.GC_UV_1018},
+              type = 'UV')
+
+V_UV_1085 = CTVertex(name = 'V_UV_1085',
+              particles = [ P.t__tilde__, P.s, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS3 ],
+              loop_particles = [[[P.G,P.t,P.s]]],
+              couplings = {(0,0,0):C.GC_UV_1019},
+              type = 'UV')
+
+V_UV_1086 = CTVertex(name = 'V_UV_1086',
+              particles = [ P.u__tilde__, P.b, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS8 ],
+              loop_particles = [[[P.G,P.u,P.b]]],
+              couplings = {(0,0,0):C.GC_UV_1013},
+              type = 'UV')
+
+V_UV_1087 = CTVertex(name = 'V_UV_1087',
+              particles = [ P.c__tilde__, P.b, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS8, L.FFS3 ],
+              loop_particles = [[[P.G,P.c,P.b]]],
+              couplings = {(0,0,0):C.GC_UV_1014,(0,1,0):C.GC_UV_1020},
+              type = 'UV')
+
+V_UV_1116 = CTVertex(name = 'V_UV_1116',
+               particles = [ P.b__tilde__, P.u, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS3 ],
+               loop_particles = [[[P.G,P.u,P.b]]],
+               couplings = {(0,0,0):C.GC_UV_1028},
+               type = 'UV')
+
+V_UV_1117 = CTVertex(name = 'V_UV_1117',
+               particles = [ P.d__tilde__, P.c, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               loop_particles = [[[P.G,P.c,P.d]]],
+               couplings = {(0,0,0):C.GC_UV_1022},
+               type = 'UV')
+
+V_UV_1118 = CTVertex(name = 'V_UV_1118',
+               particles = [ P.s__tilde__, P.c, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               loop_particles = [[[P.G,P.c,P.s]]],
+               couplings = {(0,0,0):C.GC_UV_1023},
+               type = 'UV')
+
+V_UV_1119 = CTVertex(name = 'V_UV_1119',
+               particles = [ P.b__tilde__, P.c, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8, L.FFS3 ],
+               loop_particles = [[[P.G,P.c,P.b]]],
+               couplings = {(0,0,0):C.GC_UV_1024,(0,1,0):C.GC_UV_1029},
+               type = 'UV')
+
+V_UV_1120 = CTVertex(name = 'V_UV_1120',
+               particles = [ P.d__tilde__, P.t, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               loop_particles = [[[P.G,P.t,P.d]]],
+               couplings = {(0,0,0):C.GC_UV_1025},
+               type = 'UV')
+
+V_UV_1121 = CTVertex(name = 'V_UV_1121',
+               particles = [ P.s__tilde__, P.t, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               loop_particles = [[[P.G,P.s,P.t]]],
+               couplings = {(0,0,0):C.GC_UV_1026},
+               type = 'UV')
+
diff --git a/Higgs/loop_sm_c3d4/__init__.py b/Higgs/loop_sm_c3d4/__init__.py
new file mode 100755
index 0000000000000000000000000000000000000000..cf0df43ab5cda3285f6d93c24c85c03bf73621ca
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/__init__.py
@@ -0,0 +1,29 @@
+import function_library 
+import object_library 
+import particles
+import couplings
+import CT_couplings
+import lorentz
+import parameters
+import CT_parameters
+import vertices
+import CT_vertices
+import write_param_card
+import coupling_orders
+
+# model options
+gauge = [0, 1]
+
+all_particles = particles.all_particles
+all_vertices = vertices.all_vertices
+all_CTvertices = vertices.all_CTvertices
+all_couplings = couplings.all_couplings
+all_lorentz = lorentz.all_lorentz
+all_parameters = parameters.all_parameters
+all_CTparameters = CT_parameters.all_CTparameters
+all_functions = function_library.all_functions
+all_orders = coupling_orders.all_orders
+
+__author__ = "V. Hirschi"
+__version__ = "1.2"
+__email__ = "valentin.hirschi@gmail.com"
diff --git a/Higgs/loop_sm_c3d4/coupling_orders.py b/Higgs/loop_sm_c3d4/coupling_orders.py
new file mode 100755
index 0000000000000000000000000000000000000000..784ab5b0009abbeafc3f7f2b20868744e4bb1ecf
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/coupling_orders.py
@@ -0,0 +1,16 @@
+# This file was automatically created by FeynRules $Revision: 623 $
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (November 6, 2010)
+# Date: Thu 9 Jun 2011 17:49:34
+
+
+from object_library import all_orders, CouplingOrder
+
+QCD = CouplingOrder(name = 'QCD',
+                    hierarchy = 1,
+                    expansion_order = -1,
+                    perturbative_expansion = 1)
+
+QED = CouplingOrder(name = 'QED',
+                    hierarchy = 2,
+                    expansion_order = -1,
+                    perturbative_expansion =0)
diff --git a/Higgs/loop_sm_c3d4/couplings.py b/Higgs/loop_sm_c3d4/couplings.py
new file mode 100755
index 0000000000000000000000000000000000000000..e2579a05f3393ed4111e51e403bc5db0924525c6
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/couplings.py
@@ -0,0 +1,431 @@
+# This file was automatically created by FeynRules $Revision: 535 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Fri 18 Mar 2011 18:40:51
+
+
+from object_library import all_couplings, Coupling
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec
+
+GC_HHHH = Coupling(name = 'GC_HHHH',
+                 value = '-6*complex(0,1)*lam*(1 + d4)',
+                order = {'QED':2})
+
+GC_1 = Coupling(name = 'GC_1',
+                 value = '-(ee*complex(0,1))/3.',
+                order = {'QED':1})
+
+GC_2 = Coupling(name = 'GC_2',
+                value = '(2*ee*complex(0,1))/3.',
+                order = {'QED':1})
+
+GC_3 = Coupling(name = 'GC_3',
+                value = '-(ee*complex(0,1))',
+                order = {'QED':1})
+
+GC_4 = Coupling(name = 'GC_4',
+                value = '-G',
+                order = {'QCD':1})
+
+GC_5 = Coupling(name = 'GC_5',
+                value = 'complex(0,1)*G',
+                order = {'QCD':1})
+
+GC_6 = Coupling(name = 'GC_6',
+                value = 'complex(0,1)*G**2',
+                order = {'QCD':2})
+
+GC_7 = Coupling(name = 'GC_7',
+                value = 'cw*complex(0,1)*gw',
+                order = {'QED':1})
+
+GC_8 = Coupling(name = 'GC_8',
+                value = '-(complex(0,1)*gw**2)',
+                order = {'QED':2})
+
+GC_9 = Coupling(name = 'GC_9',
+                value = 'cw**2*complex(0,1)*gw**2',
+                order = {'QED':2})
+
+GC_10 = Coupling(name = 'GC_10',
+                 value = '(ee**2*complex(0,1))/(2.*sw**2)',
+                 order = {'QED':2})
+
+GC_11 = Coupling(name = 'GC_11',
+                 value = '(ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_12 = Coupling(name = 'GC_12',
+                 value = '(CKM11*ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_13 = Coupling(name = 'GC_13',
+                 value = '(CKM12*ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_14 = Coupling(name = 'GC_14',
+                 value = '(CKM13*ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_15 = Coupling(name = 'GC_15',
+                 value = '(CKM21*ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_16 = Coupling(name = 'GC_16',
+                 value = '(CKM22*ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_17 = Coupling(name = 'GC_17',
+                 value = '(CKM23*ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_18 = Coupling(name = 'GC_18',
+                 value = '(CKM31*ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_19 = Coupling(name = 'GC_19',
+                 value = '(CKM32*ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_20 = Coupling(name = 'GC_20',
+                 value = '(CKM33*ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_21 = Coupling(name = 'GC_21',
+                 value = '-(cw*ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_22 = Coupling(name = 'GC_22',
+                 value = '(cw*ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_23 = Coupling(name = 'GC_23',
+                 value = '-(ee*complex(0,1)*sw)/(6.*cw)',
+                 order = {'QED':1})
+
+GC_24 = Coupling(name = 'GC_24',
+                 value = '(ee*complex(0,1)*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_25 = Coupling(name = 'GC_25',
+                 value = 'complex(0,1)*gw*sw',
+                 order = {'QED':1})
+
+GC_26 = Coupling(name = 'GC_26',
+                 value = '-2*cw*complex(0,1)*gw**2*sw',
+                 order = {'QED':2})
+
+GC_27 = Coupling(name = 'GC_27',
+                 value = 'complex(0,1)*gw**2*sw**2',
+                 order = {'QED':2})
+
+GC_28 = Coupling(name = 'GC_28',
+                 value = '(cw*ee*complex(0, 1))/(2.*sw) + (ee*complex(0,1)*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_29 = Coupling(name = 'GC_29',
+                 value = 'ee**2*complex(0,1) + (cw**2*ee**2*complex(0,1))/(2.*sw**2) + (ee**2*complex(0,1)*sw**2)/(2.*cw**2)',
+                 order = {'QED':2})
+
+GC_30 = Coupling(name = 'GC_30',
+                 value = '-6*complex(0,1)*lam*v*(1 + c3)',
+                 order = {'QED':1})
+
+GC_31 = Coupling(name = 'GC_31',
+                 value = '(ee**2*complex(0,1)*v)/(2.*sw**2)',
+                 order = {'QED':1})
+
+GC_32 = Coupling(name = 'GC_32',
+                 value = 'ee**2*complex(0,1)*v + (cw**2*ee**2*complex(0,1)*v)/(2.*sw**2) + (ee**2*complex(0,1)*sw**2*v)/(2.*cw**2)',
+                 order = {'QED':1})
+
+GC_33 = Coupling(name = 'GC_33',
+                 value = '-((complex(0,1)*yb)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_34 = Coupling(name = 'GC_34',
+                 value = '-((complex(0,1)*yc)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_35 = Coupling(name = 'GC_35',
+                 value = '-((complex(0,1)*ye)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_36 = Coupling(name = 'GC_36',
+                 value = '-((complex(0,1)*ym)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_37 = Coupling(name = 'GC_37',
+                 value = '-((complex(0,1)*yt)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_38 = Coupling(name = 'GC_38',
+                 value = '-((complex(0,1)*ytau)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_39 = Coupling(name = 'GC_39',
+                 value = '(ee*complex(0,1)*complexconjugate(CKM11))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_40 = Coupling(name = 'GC_40',
+                 value = '(ee*complex(0,1)*complexconjugate(CKM12))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_41 = Coupling(name = 'GC_41',
+                 value = '(ee*complex(0,1)*complexconjugate(CKM13))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_42 = Coupling(name = 'GC_42',
+                 value = '(ee*complex(0,1)*complexconjugate(CKM21))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_43 = Coupling(name = 'GC_43',
+                 value = '(ee*complex(0,1)*complexconjugate(CKM22))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_44 = Coupling(name = 'GC_44',
+                 value = '(ee*complex(0,1)*complexconjugate(CKM23))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_45 = Coupling(name = 'GC_45',
+                 value = '(ee*complex(0,1)*complexconjugate(CKM31))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_46 = Coupling(name = 'GC_46',
+                 value = '(ee*complex(0,1)*complexconjugate(CKM32))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_47 = Coupling(name = 'GC_47',
+                 value = '(ee*complex(0,1)*complexconjugate(CKM33))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+
+# == Additional coupling for the goldstones for Feynman gauge
+
+GC_1033 = Coupling(name = 'GC_1033',
+                 value = '-6*complex(0,1)*lam',
+                 order = {'QED':2})
+
+GC_1031 = Coupling(name = 'GC_1031',
+                 value = '-2*complex(0,1)*lam',
+                 order = {'QED':2})
+
+GC_1032 = Coupling(name = 'GC_1032',
+                 value = '-4*complex(0,1)*lam',
+                 order = {'QED':2})
+
+GC_1068 = Coupling(name = 'GC_1068',
+                 value = '-2*complex(0,1)*lam*v',
+                 order = {'QED':1})
+
+GC_1006 = Coupling(name = 'GC_1006',
+                value = '2*ee**2*complex(0,1)',
+                order = {'QED':2})
+
+GC_1003 = Coupling(name = 'GC_1003',
+                value = '-(ee*complex(0,1))',
+                order = {'QED':1})
+
+GC_1055 = Coupling(name = 'GC_1055',
+                 value = '-(ee**2*complex(0,1))/(2.*sw)',
+                 order = {'QED':2})
+
+GC_1054 = Coupling(name = 'GC_1054',
+                 value = '-ee**2/(2.*sw)',
+                 order = {'QED':2})
+
+GC_1074 = Coupling(name = 'GC_1074',
+                 value = '-(ee**2*v)/(2.*sw)',
+                 order = {'QED':1})
+
+GC_1039 = Coupling(name = 'GC_1039',
+                 value = '(ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_1037 = Coupling(name = 'GC_1037',
+                 value = '-ee/(2.*sw)',
+                 order = {'QED':1})
+
+GC_1056 = Coupling(name = 'GC_1056',
+                 value = 'ee**2/(2.*sw)',
+                 order = {'QED':2})
+
+GC_1075 = Coupling(name = 'GC_1075',
+                 value = '(ee**2*v)/(2.*sw)',
+                 order = {'QED':1})
+
+GC_1038 = Coupling(name = 'GC_1038',
+                 value = '-(ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_1034 = Coupling(name = 'GC_1034',
+                 value = '(ee**2*complex(0,1))/(2.*sw**2)',
+                 order = {'QED':2})
+
+GC_1063 = Coupling(name = 'GC_1063',
+                 value = '(cw*ee**2*complex(0,1))/sw - (ee**2*complex(0,1)*sw)/cw',
+                 order = {'QED':2})
+
+GC_1060 = Coupling(name = 'GC_1060',
+                 value = '-(cw*ee)/(2.*sw) - (ee*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_1061 = Coupling(name = 'GC_1061',
+                 value = '-(cw*ee*complex(0,1))/(2.*sw) + (ee*complex(0,1)*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_1008 = Coupling(name = 'GC_1008',
+                value = '(ee**2*complex(0,1))/(2.*cw)',
+                order = {'QED':2})
+
+GC_1009 = Coupling(name = 'GC_1009',
+                value = 'ee**2/(2.*cw)',
+                order = {'QED':2})
+
+GC_1067 = Coupling(name = 'GC_1067',
+                 value = '(ee**2*v)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_1007 = Coupling(name = 'GC_1007',
+                value = '-ee**2/(2.*cw)',
+                order = {'QED':2})
+
+GC_1066 = Coupling(name = 'GC_1066',
+                 value = '-(ee**2*v)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_1065 = Coupling(name = 'GC_1065',
+                 value = 'ee**2*complex(0,1) + (cw**2*ee**2*complex(0,1))/(2.*sw**2) + (ee**2*complex(0,1)*sw**2)/(2.*cw**2)',
+                 order = {'QED':2})
+
+GC_1064 = Coupling(name = 'GC_1064',
+                 value = '-(ee**2*complex(0,1)) + (cw**2*ee**2*complex(0,1))/(2.*sw**2) + (ee**2*complex(0,1)*sw**2)/(2.*cw**2)',
+                 order = {'QED':2})
+
+GC_1082 = Coupling(name = 'GC_1082',
+                 value = '-(yb/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_1016 = Coupling(name = 'GC_1016',
+                 value = '-I2x12',
+                 order = {'QED':1})
+
+GC_1017 = Coupling(name = 'GC_1017',
+                 value = '-I2x13',
+                 order = {'QED':1})
+
+GC_1018 = Coupling(name = 'GC_1018',
+                 value = '-I2x22',
+                 order = {'QED':1})
+
+GC_1019 = Coupling(name = 'GC_1019',
+                 value = '-I2x23',
+                 order = {'QED':1})
+
+GC_1013 = Coupling(name = 'GC_1013',
+                 value = 'I1x31',
+                 order = {'QED':1})
+
+GC_1014 = Coupling(name = 'GC_1014',
+                 value = 'I1x32',
+                 order = {'QED':1})
+
+GC_1020 = Coupling(name = 'GC_1020',
+                 value = '-I2x32',
+                 order = {'QED':1})
+
+GC_1015 = Coupling(name = 'GC_1015',
+                 value = 'I1x33',
+                 order = {'QED':1})
+
+GC_1021 = Coupling(name = 'GC_1021',
+                 value = '-I2x33',
+                 order = {'QED':1})
+
+GC_1088 = Coupling(name = 'GC_1088',
+                 value = '-(ye/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_1092 = Coupling(name = 'GC_1092',
+                 value = '-(ym/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_1098 = Coupling(name = 'GC_1098',
+                 value = '-(ytau/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_1087 = Coupling(name = 'GC_1087',
+                 value = 'ye',
+                 order = {'QED':1})
+
+GC_1091 = Coupling(name = 'GC_1091',
+                 value = 'ym',
+                 order = {'QED':1})
+
+GC_1097 = Coupling(name = 'GC_1097',
+                 value = 'ytau',
+                 order = {'QED':1})
+
+GC_1028 = Coupling(name = 'GC_1028',
+                 value = '-I4x13',
+                 order = {'QED':1})
+
+GC_1022 = Coupling(name = 'GC_1022',
+                 value = 'I3x21',
+                 order = {'QED':1})
+
+GC_1023 = Coupling(name = 'GC_1023',
+                 value = 'I3x22',
+                 order = {'QED':1})
+
+GC_1024 = Coupling(name = 'GC_1024',
+                 value = 'I3x23',
+                 order = {'QED':1})
+
+GC_1025 = Coupling(name = 'GC_1025',
+                 value = 'I3x31',
+                 order = {'QED':1})
+
+GC_1029 = Coupling(name = 'GC_1029',
+                 value = '-I4x23',
+                 order = {'QED':1})
+
+GC_1026 = Coupling(name = 'GC_1026',
+                 value = 'I3x32',
+                 order = {'QED':1})
+
+GC_1027 = Coupling(name = 'GC_1027',
+                 value = 'I3x33',
+                 order = {'QED':1})
+
+GC_1030 = Coupling(name = 'GC_1030',
+                 value = '-I4x33',
+                 order = {'QED':1})
+
+GC_1085 = Coupling(name = 'GC_1085',
+                 value = 'yc/cmath.sqrt(2)',
+                 order = {'QED':1})
+
+GC_1095 = Coupling(name = 'GC_1095',
+                 value = 'yt/cmath.sqrt(2)',
+                 order = {'QED':1})
+
+GC_1086 = Coupling(name = 'GC_1086',
+                 value = '-ye',
+                 order = {'QED':1})
+
+GC_1090 = Coupling(name = 'GC_1090',
+                 value = '-ym',
+                 order = {'QED':1})
+
+GC_1096 = Coupling(name = 'GC_1096',
+                 value = '-ytau',
+                 order = {'QED':1})
+
+
+
+
+
+
+
+
diff --git a/Higgs/loop_sm_c3d4/function_library.py b/Higgs/loop_sm_c3d4/function_library.py
new file mode 100755
index 0000000000000000000000000000000000000000..375b74e41d613cf4647b00b1ea9963a189ddb158
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/function_library.py
@@ -0,0 +1,64 @@
+# This file is part of the UFO.
+#
+# This file contains definitions for functions that
+# are extensions of the cmath library, and correspond
+# either to functions that are in cmath, but inconvenient
+# to access from there (e.g. z.conjugate()),
+# or functions that are simply not defined.
+#
+#
+
+__date__ = "22 July 2010"
+__author__ = "claude.duhr@durham.ac.uk"
+
+import cmath
+from object_library import all_functions, Function
+
+#
+# shortcuts for functions from cmath
+#
+
+complexconjugate = Function(name = 'complexconjugate',
+                            arguments = ('z',),
+                            expression = 'z.conjugate()')
+
+
+re = Function(name = 're',
+              arguments = ('z',),
+              expression = 'z.real')
+
+im = Function(name = 'im',
+              arguments = ('z',),
+              expression = 'z.imag')
+
+# Auxiliary functions for NLO
+
+cond = Function(name = 'cond',
+                arguments = ('condition','ExprTrue','ExprFalse'),
+                expression = '(ExprTrue if condition==0.0 else ExprFalse)')
+
+reglog = Function(name = 'reglog',
+                arguments = ('z'),
+                expression = '(0.0 if z==0.0 else cmath.log(z))')
+
+# New functions (trigonometric)
+
+sec = Function(name = 'sec',
+             arguments = ('z',),
+             expression = '1./cmath.cos(z)')
+
+asec = Function(name = 'asec',
+             arguments = ('z',),
+             expression = 'cmath.acos(1./z)')
+
+csc = Function(name = 'csc',
+             arguments = ('z',),
+             expression = '1./cmath.sin(z)')
+
+acsc = Function(name = 'acsc',
+             arguments = ('z',),
+             expression = 'cmath.asin(1./z)')
+
+
+
+
diff --git a/Higgs/loop_sm_c3d4/lorentz.py b/Higgs/loop_sm_c3d4/lorentz.py
new file mode 100755
index 0000000000000000000000000000000000000000..2b5938d37283714deed80036354c9c464e5a8d48
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/lorentz.py
@@ -0,0 +1,185 @@
+# This file was automatically created by FeynRules $Revision: 535 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Fri 18 Mar 2011 18:40:51
+
+
+from object_library import all_lorentz, Lorentz
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec
+
+###################################
+# CounterTerms Lorentz structures #
+###################################
+
+R2_GG_1 = Lorentz(name = 'R2_GG_1',
+               spins = [ 3, 3 ],
+               structure = 'P(-1,1)*P(-1,1)*Metric(1,2)')
+
+R2_GG_2 = Lorentz(name = 'R2_GG_2',
+               spins = [ 3, 3 ],
+               structure = 'P(1,1)*P(2,1)')
+
+R2_GG_3 = Lorentz(name = 'R2_GG_3',
+               spins = [ 3, 3 ],
+               structure = 'Metric(1,2)')
+
+R2_QQ_1 = Lorentz(name = 'R2_QQ_1',
+               spins = [ 2, 2 ],
+               structure = 'P(-1,1)*Gamma(-1,2,1)')
+
+R2_QQ_2 = Lorentz(name = 'R2_QQ_2',
+               spins = [ 2, 2 ],
+               structure = 'Identity(1,2)')
+
+GHGHG = Lorentz(name = 'GHGHG',
+                 spins = [ 1, 1, 3 ],
+                structure = 'P(3,2)')
+
+#=============================================================================================
+#  4-gluon R2 vertex
+#=============================================================================================
+
+
+R2_4G_1234 = Lorentz(name = 'R2_4G_1234',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,2)*Metric(3,4)')
+
+R2_4G_1324 = Lorentz(name = 'R2_4G_1324',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,3)*Metric(2,4)')
+
+R2_4G_1423 = Lorentz(name = 'R2_4G_1423',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3)')
+
+R2RGA_VVVV10 = Lorentz(name = 'R2RGA_VVVV10',
+                 spins = [ 3, 3, 3, 3 ],
+                 structure = 'Metric(1,4)*Metric(2,3) + Metric(1,3)*Metric(2,4) + Metric(1,2)*Metric(3,4)')
+
+R2RGA_VVVV2 = Lorentz(name = 'R2RGA_VVVV2',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3)')
+
+R2RGA_VVVV3 = Lorentz(name = 'R2RGA_VVVV3',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,3)*Metric(2,4)')
+
+R2RGA_VVVV5 = Lorentz(name = 'R2RGA_VVVV5',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,2)*Metric(3,4)')
+
+#=============================================================================================
+
+R2_GGZ = Lorentz(name = 'R2_GGZ',
+                 spins = [ 3, 3, 3 ],
+                 structure = 'Epsilon(3,1,2,-1)*P(-1,2)-Epsilon(3,1,2,-1)*P(-1,1)') 
+
+R2_GGVV = Lorentz(name = 'R2_GGVV',
+                 spins = [ 3, 3, 3, 3 ],
+                 structure = 'Metric(1,2)*Metric(3,4)+Metric(1,3)*Metric(2,4)+Metric(1,4)*Metric(2,3)')
+
+R2_GGHH = Lorentz(name = 'R2_GGHH',
+                 spins = [ 3, 3, 1, 1 ],
+                 structure = 'Metric(1,2)')
+
+R2_GGGVa = Lorentz(name = 'R2_GGGVa',
+                 spins = [ 3, 3, 3, 3 ],
+                 structure = 'Epsilon(4,1,2,3)')
+
+###################
+# Base structures #
+###################
+
+UUV1 = Lorentz(name = 'UUV1',
+               spins = [ -1, -1, 3 ],
+               structure = 'P(3,2) + P(3,3)')
+
+SSS1 = Lorentz(name = 'SSS1',
+               spins = [ 1, 1, 1 ],
+               structure = '1')
+
+FFS1 = Lorentz(name = 'FFS1',
+               spins = [ 2, 2, 1 ],
+               structure = 'Identity(2,1)')
+
+FFV1 = Lorentz(name = 'FFV1',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,1)')
+
+FFV2 = Lorentz(name = 'FFV2',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,-1)*ProjM(-1,1)')
+
+FFV3 = Lorentz(name = 'FFV3',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,-1)*ProjM(-1,1) - 2*Gamma(3,2,-1)*ProjP(-1,1)')
+
+FFV4 = Lorentz(name = 'FFV4',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,-1)*ProjM(-1,1) + 2*Gamma(3,2,-1)*ProjP(-1,1)')
+
+FFV5 = Lorentz(name = 'FFV5',
+               spins = [ 2, 2, 3 ],
+               structure = 'Gamma(3,2,-1)*ProjM(-1,1) + 4*Gamma(3,2,-1)*ProjP(-1,1)')
+
+VVS1 = Lorentz(name = 'VVS1',
+               spins = [ 3, 3, 1 ],
+               structure = 'Metric(1,2)')
+
+VVV1 = Lorentz(name = 'VVV1',
+               spins = [ 3, 3, 3 ],
+               structure = 'P(3,1)*Metric(1,2) - P(3,2)*Metric(1,2) - P(2,1)*Metric(1,3) + P(2,3)*Metric(1,3) + P(1,2)*Metric(2,3) - P(1,3)*Metric(2,3)')
+
+VVSS1 = Lorentz(name = 'VVSS1',
+                spins = [ 3, 3, 1, 1 ],
+                structure = 'Metric(1,2)')
+
+VVVV1 = Lorentz(name = 'VVVV1',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) - Metric(1,3)*Metric(2,4)')
+
+VVVV2 = Lorentz(name = 'VVVV2',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) + Metric(1,3)*Metric(2,4) - 2*Metric(1,2)*Metric(3,4)')
+
+VVVV3 = Lorentz(name = 'VVVV3',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) - Metric(1,2)*Metric(3,4)')
+
+VVVV4 = Lorentz(name = 'VVVV4',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,3)*Metric(2,4) - Metric(1,2)*Metric(3,4)')
+
+VVVV5 = Lorentz(name = 'VVVV5',
+                spins = [ 3, 3, 3, 3 ],
+                structure = 'Metric(1,4)*Metric(2,3) - (Metric(1,3)*Metric(2,4))/2. - (Metric(1,2)*Metric(3,4))/2.')
+
+# == Additional lorentz structure for Fenyman goldstones
+
+UUS1 = Lorentz(name = 'UUS1',
+               spins = [ -1, -1, 1 ],
+               structure = '1')
+
+FFS8 = Lorentz(name = 'FFS8',
+               spins = [ 2, 2, 1 ],
+               structure = 'ProjM(2,1)')
+
+FFS2 = Lorentz(name = 'FFS2',
+               spins = [ 2, 2, 1 ],
+               structure = 'ProjM(2,1) - ProjP(2,1)')
+
+FFS3 = Lorentz(name = 'FFS3',
+               spins = [ 2, 2, 1 ],
+               structure = 'ProjP(2,1)')
+
+FFS4 = Lorentz(name = 'FFS4',
+               spins = [ 2, 2, 1 ],
+               structure = 'ProjM(2,1) + ProjP(2,1)')
+
+VSS1 = Lorentz(name = 'VSS1',
+               spins = [ 3, 1, 1 ],
+               structure = 'P(1,2) - P(1,3)')
+
+SSSS1 = Lorentz(name = 'SSSS1',
+                spins = [ 1, 1, 1, 1 ],
+                structure = '1')
diff --git a/Higgs/loop_sm_c3d4/object_library.py b/Higgs/loop_sm_c3d4/object_library.py
new file mode 100755
index 0000000000000000000000000000000000000000..af97b3f70e66e4827ddd091e9088a009ee8608a8
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/object_library.py
@@ -0,0 +1,312 @@
+##
+##
+## Feynrules Header
+##
+##
+##
+##
+##
+
+import cmath
+import re
+
+
+class UFOError(Exception):
+        """Exception raised if when inconsistencies are detected in the UFO model."""
+        pass
+
+class UFOBaseClass(object):
+    """The class from which all FeynRules classes are derived."""
+
+    require_args = []
+
+    def __init__(self, *args, **options):
+        assert(len(self.require_args) == len (args))
+    
+        for i, name in enumerate(self.require_args):
+            setattr(self, name, args[i])
+    
+        for (option, value) in options.items():
+            setattr(self, option, value)
+
+    def get(self, name):
+        return getattr(self, name)
+    
+    def set(self, name, value):
+        setattr(self, name, value)
+        
+    def get_all(self):
+        """Return a dictionary containing all the information of the object"""
+        return self.__dict__
+
+    def __str__(self):
+        return self.name
+
+    def nice_string(self):
+        """ return string with the full information """
+        return '\n'.join(['%s \t: %s' %(name, value) for name, value in self.__dict__.items()])
+
+    def __repr__(self):
+        replacements = [
+            ('+','__plus__'),
+            ('-','__minus__'),
+            ('@','__at__'),
+            ('!','__exclam__'),
+            ('?','__quest__'),
+            ('*','__star__'),
+            ('~','__tilde__')
+            ]
+        text = self.name
+        for orig,sub in replacements:
+            text = text.replace(orig,sub)
+        return text
+
+
+
+all_particles = []
+
+    
+
+class Particle(UFOBaseClass):
+    """A standard Particle""" 
+
+    require_args=['pdg_code', 'name', 'antiname', 'spin', 'color', 'mass', 'width', 'texname', 'antitexname', 'charge']
+
+    require_args_all = ['pdg_code', 'name', 'antiname', 'spin', 'color', 'mass', 'width', 'texname', 'antitexname', 'charge', 'loop_particles', 'counterterm','line', 'propagating', 'goldstoneboson']
+
+    def __init__(self, pdg_code, name, antiname, spin, color, mass, width, texname,
+                 antitexname, charge , loop_particles=None, counterterm=None, line=None, propagating=True, goldstoneboson=False, **options):
+
+        args= (pdg_code, name, antiname, spin, color, mass, width, texname,
+                 antitexname, float(charge))
+
+        UFOBaseClass.__init__(self, *args,  **options)
+
+        global all_particles
+        all_particles.append(self)
+
+        self.propagating = propagating
+        self.goldstoneboson= goldstoneboson
+
+        self.selfconjugate = (name == antiname)
+        if 1: #not line:
+            self.line = self.find_line_type()
+        else:
+            self.line = line
+
+    def find_line_type(self):
+        """ find how we draw a line if not defined
+        valid output: dashed/straight/wavy/curly/double/swavy/scurly
+        """
+        
+        spin = self.spin
+        color = self.color
+        
+        #use default
+        if spin == 1:
+            return 'dashed'
+        elif spin == 2:
+            if not self.selfconjugate:
+                return 'straight'
+            elif color == 1:
+                return 'swavy'
+            else:
+                return 'scurly'
+        elif spin == 3:
+            if color == 1:
+                return 'wavy'
+            
+            else:
+                return 'curly'
+        elif spin == 5:
+            return 'double'
+        elif spin == -1:
+            return 'dotted'
+        else:
+            return 'dashed' # not supported yet
+
+    def anti(self):
+        # We do not copy the UV wavefunction renormalization as it is defined for the particle only.
+        if self.selfconjugate:
+            raise Exception('%s has no anti particle.' % self.name) 
+        outdic = {}
+        for k,v in self.__dict__.items():
+            if k not in self.require_args_all:                
+                outdic[k] = -v
+        if self.color in [1,8]:
+            newcolor = self.color
+        else:
+            newcolor = -self.color
+                
+        return Particle(-self.pdg_code, self.antiname, self.name, self.spin, newcolor, self.mass, self.width,
+                        self.antitexname, self.texname, -self.charge, self.line, self.propagating, self.goldstoneboson, **outdic)
+
+
+
+all_parameters = []
+
+class Parameter(UFOBaseClass):
+
+    require_args=['name', 'nature', 'type', 'value', 'texname']
+
+    def __init__(self, name, nature, type, value, texname, lhablock=None, lhacode=None, loop_particles=None, counterterm=None):
+
+        args = (name,nature,type,value,texname)
+
+        UFOBaseClass.__init__(self, *args)
+
+        args=(name,nature,type,value,texname)
+
+        global all_parameters
+        all_parameters.append(self)
+
+        if (lhablock is None or lhacode is None)  and nature == 'external':
+            raise Exception('Need LHA information for external parameter "%s".' % name)
+        self.lhablock = lhablock
+        self.lhacode = lhacode
+
+all_CTparameters = []
+
+class CTParameter(UFOBaseClass):
+
+    require_args=['name', 'nature,', 'type', 'value', 'texname']
+
+    def __init__(self, name, type, value, texname):
+
+        args = (name,'internal',type,value,texname)
+
+        UFOBaseClass.__init__(self, *args)
+
+        args=(name,'internal',type,value,texname)
+
+        self.nature='interal'
+
+        global all_CTparameters
+        all_CTparameters.append(self)
+
+    def finite(self):
+        try:
+            return self.value[0]
+        except KeyError:
+            return 'ZERO'
+    
+    def pole(self, x):
+        try:
+            return self.value[-x]
+        except KeyError:
+            return 'ZERO'
+
+all_vertices = []
+
+class Vertex(UFOBaseClass):
+
+    require_args=['name', 'particles', 'color', 'lorentz', 'couplings']
+
+    def __init__(self, name, particles, color, lorentz, couplings, **opt):
+ 
+        args = (name, particles, color, lorentz, couplings)
+
+        UFOBaseClass.__init__(self, *args, **opt)
+
+        args=(particles,color,lorentz,couplings)
+        
+        global all_vertices
+        all_vertices.append(self)
+
+all_CTvertices = []
+
+class CTVertex(UFOBaseClass):
+
+    require_args=['name', 'particles', 'color', 'lorentz', 'couplings', 'type', 'loop_particles']
+
+    def __init__(self, name, particles, color, lorentz, couplings, type, loop_particles, **opt):
+ 
+        args = (name, particles, color, lorentz, couplings, type, loop_particles)
+
+        UFOBaseClass.__init__(self, *args, **opt)
+
+        args=(particles,color,lorentz,couplings, type, loop_particles)
+        
+        global all_CTvertices
+        all_CTvertices.append(self)
+
+all_couplings = []
+
+class Coupling(UFOBaseClass):
+
+    require_args=['name', 'value', 'order']
+
+    require_args_all=['name', 'value', 'order', 'loop_particles', 'counterterm']
+
+    def __init__(self, name, value, order, loop_particles=None, counterterm=None, **opt):
+
+        args =(name, value, order)	
+        UFOBaseClass.__init__(self, *args, **opt)
+        global all_couplings
+        all_couplings.append(self)
+   
+    def value(self):
+        return self.pole(0)
+
+    def pole(self, x):
+        """ the self.value attribute can be a dictionary directly specifying the Laurent serie using normal
+        parameter or just a string which can possibly contain CTparameter defining the Laurent serie."""
+        
+        if isinstance(self.value,dict):
+            if -x in self.value.keys():
+                return self.value[-x]
+            else:
+                return 'ZERO'
+
+        if x==0:
+            return self.value
+        else:
+            return 'ZERO'
+
+all_lorentz = []
+
+class Lorentz(UFOBaseClass):
+
+    require_args=['name','spins','structure']
+    
+    def __init__(self, name, spins, structure='external', **opt):
+        args = (name, spins, structure)
+        UFOBaseClass.__init__(self, *args, **opt)
+
+        global all_lorentz
+        all_lorentz.append(self)
+
+
+all_functions = []
+
+class Function(object):
+
+    def __init__(self, name, arguments, expression):
+
+        global all_functions
+        all_functions.append(self)
+
+        self.name = name
+        self.arguments = arguments
+        self.expr = expression
+    
+    def __call__(self, *opt):
+
+        for i, arg in enumerate(self.arguments):
+            exec('%s = %s' % (arg, opt[i] ))
+
+        return eval(self.expr)
+
+all_orders = []
+
+class CouplingOrder(object):
+
+    def __init__(self, name, expansion_order, hierarchy, perturbative_expansion = 0):
+        
+        global all_orders
+        all_orders.append(self)
+
+        self.name = name
+        self.expansion_order = expansion_order
+        self.hierarchy = hierarchy
+        self.perturbative_expansion = perturbative_expansion
diff --git a/Higgs/loop_sm_c3d4/parameters.py b/Higgs/loop_sm_c3d4/parameters.py
new file mode 100755
index 0000000000000000000000000000000000000000..84d6bf9be074002c808c1f7a0c70e93ce8c2b92e
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/parameters.py
@@ -0,0 +1,838 @@
+# This file was automatically created by FeynRules $Revision: 634 $
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (November 6, 2010)
+# Date: Wed 6 Jul 2011 14:07:37
+
+
+
+from object_library import all_parameters, Parameter
+
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec
+
+# This is a default parameter object representing 0.
+ZERO = Parameter(name = 'ZERO',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '0.0',
+                 texname = '0')
+
+# Loop related parameters
+
+lhv = Parameter(name = 'lhv',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '1.0',
+                 texname = '\lambda_{HV}')
+
+Ncol = Parameter(name = 'Ncol',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '3.0',
+                 texname = 'N_{col}')
+
+CA = Parameter(name = 'CA',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '3.0',
+                 texname = 'C_{A}')
+
+TF = Parameter(name = 'TF',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '0.5',
+                 texname = 'T_{F}')
+
+CF = Parameter(name = 'CF',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '(4.0/3.0)',
+                 texname = 'C_{F}')
+
+MU_R = Parameter(name = 'MU_R',
+              nature = 'external',
+              type = 'real',
+              value = 91.188,
+              texname = '\\text{\\mu_r}',
+              lhablock = 'LOOP',
+              lhacode = [ 1 ])
+
+# User-defined parameters.
+aEWM1 = Parameter(name = 'aEWM1',
+                  nature = 'external',
+                  type = 'real',
+                  value = 132.50698,
+                  texname = '\\text{aEWM1}',
+                  lhablock = 'SMINPUTS',
+                  lhacode = [ 1 ])
+
+Gf = Parameter(name = 'Gf',
+               nature = 'external',
+               type = 'real',
+               value = 0.0000116639,
+               texname = 'G_f',
+               lhablock = 'SMINPUTS',
+               lhacode = [ 2 ])
+
+aS = Parameter(name = 'aS',
+               nature = 'external',
+               type = 'real',
+               value = 0.118,
+               texname = '\\text{aS}',
+               lhablock = 'SMINPUTS',
+               lhacode = [ 3 ])
+
+lamWS = Parameter(name = 'lamWS',
+                  nature = 'external',
+                  type = 'real',
+                  value = 0.2253,
+                  texname = '\\text{lamWS}',
+                  lhablock = 'Wolfenstein',
+                  lhacode = [ 1 ])
+
+AWS = Parameter(name = 'AWS',
+                nature = 'external',
+                type = 'real',
+                value = 0.808,
+                texname = '\\text{AWS}',
+                lhablock = 'Wolfenstein',
+                lhacode = [ 2 ])
+
+rhoWS = Parameter(name = 'rhoWS',
+                  nature = 'external',
+                  type = 'real',
+                  value = 0.132,
+                  texname = '\\text{rhoWS}',
+                  lhablock = 'Wolfenstein',
+                  lhacode = [ 3 ])
+
+etaWS = Parameter(name = 'etaWS',
+                  nature = 'external',
+                  type = 'real',
+                  value = 0.341,
+                  texname = '\\text{etaWS}',
+                  lhablock = 'Wolfenstein',
+                  lhacode = [ 4 ])
+
+ymc = Parameter(name = 'ymc',
+                nature = 'external',
+                type = 'real',
+                value = 1.27,
+                texname = '\\text{ymc}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 4 ])
+
+ymb = Parameter(name = 'ymb',
+                nature = 'external',
+                type = 'real',
+                value = 4.2,
+                texname = '\\text{ymb}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 5 ])
+
+ymt = Parameter(name = 'ymt',
+                nature = 'external',
+                type = 'real',
+                value = 164.5,
+                texname = '\\text{ymt}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 6 ])
+
+yme = Parameter(name = 'yme',
+                nature = 'external',
+                type = 'real',
+                value = 0.0005110000000000001,
+                texname = '\\text{yme}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 11 ])
+
+ymm = Parameter(name = 'ymm',
+                nature = 'external',
+                type = 'real',
+                value = 0.10566,
+                texname = '\\text{ymm}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 13 ])
+
+ymtau = Parameter(name = 'ymtau',
+                  nature = 'external',
+                  type = 'real',
+                  value = 1.777,
+                  texname = '\\text{ymtau}',
+                  lhablock = 'YUKAWA',
+                  lhacode = [ 15 ])
+
+MC = Parameter(name = 'MC',
+               nature = 'external',
+               type = 'real',
+               value = 1.42,
+               texname = '\\text{MC}',
+               lhablock = 'MASS',
+               lhacode = [ 4 ])
+
+MT = Parameter(name = 'MT',
+               nature = 'external',
+               type = 'real',
+               value = 172.,
+               texname = '\\text{MT}',
+               lhablock = 'MASS',
+               lhacode = [ 6 ])
+
+MB = Parameter(name = 'MB',
+               nature = 'external',
+               type = 'real',
+               value = 4.7,
+               texname = '\\text{MB}',
+               lhablock = 'MASS',
+               lhacode = [ 5 ])
+
+MZ = Parameter(name = 'MZ',
+               nature = 'external',
+               type = 'real',
+               value = 91.188,
+               texname = '\\text{MZ}',
+               lhablock = 'MASS',
+               lhacode = [ 23 ])
+
+MH = Parameter(name = 'MH',
+               nature = 'external',
+               type = 'real',
+               value = 125.,
+               texname = '\\text{MH}',
+               lhablock = 'MASS',
+               lhacode = [ 25 ])
+
+Me = Parameter(name = 'Me',
+               nature = 'external',
+               type = 'real',
+               value = 0.0005110000000000001,
+               texname = '\\text{Me}',
+               lhablock = 'MASS',
+               lhacode = [ 11 ])
+
+MM = Parameter(name = 'MM',
+               nature = 'external',
+               type = 'real',
+               value = 0.10566,
+               texname = '\\text{MM}',
+               lhablock = 'MASS',
+               lhacode = [ 13 ])
+
+MTA = Parameter(name = 'MTA',
+                nature = 'external',
+                type = 'real',
+                value = 1.777,
+                texname = '\\text{MTA}',
+                lhablock = 'MASS',
+                lhacode = [ 15 ])
+
+WT = Parameter(name = 'WT',
+               nature = 'external',
+               type = 'real',
+               value = 1.50833649,
+               texname = '\\text{WT}',
+               lhablock = 'DECAY',
+               lhacode = [ 6 ])
+
+WZ = Parameter(name = 'WZ',
+               nature = 'external',
+               type = 'real',
+               value = 2.44140351,
+               texname = '\\text{WZ}',
+               lhablock = 'DECAY',
+               lhacode = [ 23 ])
+
+WW = Parameter(name = 'WW',
+               nature = 'external',
+               type = 'real',
+               value = 2.04759951,
+               texname = '\\text{WW}',
+               lhablock = 'DECAY',
+               lhacode = [ 24 ])
+
+WH = Parameter(name = 'WH',
+               nature = 'external',
+               type = 'real',
+               value = 6.38233934e-03,
+               texname = '\\text{WH}',
+               lhablock = 'DECAY',
+               lhacode = [ 25 ])
+
+WTau = Parameter(name = 'WTau',
+                 nature = 'external',
+                 type = 'real',
+                 value = 2.27e-12,
+                 texname = '\\text{WTau}',
+                 lhablock = 'DECAY',
+                 lhacode = [ 15 ])
+
+c3  = Parameter(name = 'c3',
+                nature = 'external',
+                type = 'real',
+                value = 2.0,
+                lhablock = 'TRIPCOUP',
+                lhacode = [ 4 ],
+                texname = '\\text{c3}')
+
+d4  = Parameter(name = 'd4',
+                nature = 'external',
+                type = 'real',
+                value = 3.0,
+                lhablock = 'QUARTCOUP',
+                lhacode = [ 6 ],
+                texname = '\\text{d4}')
+
+CKM11 = Parameter(name = 'CKM11',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '1 - lamWS**2/2.',
+                  texname = '\\text{CKM11}')
+
+CKM12 = Parameter(name = 'CKM12',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'lamWS',
+                  texname = '\\text{CKM12}')
+
+CKM13 = Parameter(name = 'CKM13',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'AWS*lamWS**3*(-(etaWS*complex(0,1)) + rhoWS)',
+                  texname = '\\text{CKM13}')
+
+CKM21 = Parameter(name = 'CKM21',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-lamWS',
+                  texname = '\\text{CKM21}')
+
+CKM22 = Parameter(name = 'CKM22',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '1 - lamWS**2/2.',
+                  texname = '\\text{CKM22}')
+
+CKM23 = Parameter(name = 'CKM23',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'AWS*lamWS**2',
+                  texname = '\\text{CKM23}')
+
+CKM31 = Parameter(name = 'CKM31',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'AWS*lamWS**3*(1 - etaWS*complex(0,1) - rhoWS)',
+                  texname = '\\text{CKM31}')
+
+CKM32 = Parameter(name = 'CKM32',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-(AWS*lamWS**2)',
+                  texname = '\\text{CKM32}')
+
+CKM33 = Parameter(name = 'CKM33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '1',
+                  texname = '\\text{CKM33}')
+
+aEW = Parameter(name = 'aEW',
+                nature = 'internal',
+                type = 'real',
+                value = '1/aEWM1',
+                texname = '\\text{aEW}')
+
+G = Parameter(name = 'G',
+              nature = 'internal',
+              type = 'real',
+              value = '2*cmath.sqrt(aS)*cmath.sqrt(cmath.pi)',
+              texname = 'G')
+
+MW = Parameter(name = 'MW',
+               nature = 'internal',
+               type = 'real',
+               value = 'cmath.sqrt(MZ**2/2. + cmath.sqrt(MZ**4/4. - (aEW*cmath.pi*MZ**2)/(Gf*cmath.sqrt(2))))',
+               texname = 'M_W')
+
+ee = Parameter(name = 'ee',
+               nature = 'internal',
+               type = 'real',
+               value = '2*cmath.sqrt(aEW)*cmath.sqrt(cmath.pi)',
+               texname = 'e')
+
+sw2 = Parameter(name = 'sw2',
+                nature = 'internal',
+                type = 'real',
+                value = '1 - MW**2/MZ**2',
+                texname = '\\text{sw2}')
+
+cw = Parameter(name = 'cw',
+               nature = 'internal',
+               type = 'real',
+               value = 'cmath.sqrt(1 - sw2)',
+               texname = 'c_w')
+
+sw = Parameter(name = 'sw',
+               nature = 'internal',
+               type = 'real',
+               value = 'cmath.sqrt(sw2)',
+               texname = 's_w')
+
+g1 = Parameter(name = 'g1',
+               nature = 'internal',
+               type = 'real',
+               value = 'ee/cw',
+               texname = 'g_1')
+
+gw = Parameter(name = 'gw',
+               nature = 'internal',
+               type = 'real',
+               value = 'ee/sw',
+               texname = 'g_w')
+
+v = Parameter(name = 'v',
+              nature = 'internal',
+              type = 'real',
+              value = '(2*MW*sw)/ee',
+              texname = 'v')
+
+lam = Parameter(name = 'lam',
+                nature = 'internal',
+                type = 'real',
+                value = 'MH**2/(2.*v**2)',
+                texname = '\\text{lam}')
+
+yb = Parameter(name = 'yb',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymb*cmath.sqrt(2))/v',
+               texname = '\\text{yb}')
+
+yc = Parameter(name = 'yc',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymc*cmath.sqrt(2))/v',
+               texname = '\\text{yc}')
+
+ye = Parameter(name = 'ye',
+               nature = 'internal',
+               type = 'real',
+               value = '(yme*cmath.sqrt(2))/v',
+               texname = '\\text{ye}')
+
+ym = Parameter(name = 'ym',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymm*cmath.sqrt(2))/v',
+               texname = '\\text{ym}')
+
+yt = Parameter(name = 'yt',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymt*cmath.sqrt(2))/v',
+               texname = '\\text{yt}')
+
+ytau = Parameter(name = 'ytau',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '(ymtau*cmath.sqrt(2))/v',
+                 texname = '\\text{ytau}')
+
+muH = Parameter(name = 'muH',
+                nature = 'internal',
+                type = 'real',
+                value = 'cmath.sqrt(lam*v**2)',
+                texname = '\\mu')
+
+# To facilitate R2 writings
+
+AxialZUp = Parameter(name = 'AxialZUp',
+              nature = 'internal',
+              type = 'real',
+              value = '(3.0/2.0)*(-(ee*sw)/(6.*cw))-(1.0/2.0)*((cw*ee)/(2.*sw))',
+              texname = 'AxialZUp')
+
+AxialZDown = Parameter(name = 'AxialZDown',
+              nature = 'internal',
+              type = 'real',
+              value = '(-1.0/2.0)*(-(cw*ee)/(2.*sw))+(-3.0/2.0)*(-(ee*sw)/(6.*cw))',
+              texname = 'AxialZdown')
+
+VectorZUp = Parameter(name = 'VectorZUp',
+              nature = 'internal',                      
+              type = 'real',
+              value = '(1.0/2.0)*((cw*ee)/(2.*sw))+(5.0/2.0)*(-(ee*sw)/(6.*cw))',
+              texname = 'AxialZUp')
+
+VectorZDown = Parameter(name = 'VectorZDown',
+              nature = 'internal',                        
+              type = 'real',
+              value = '(1.0/2.0)*(-(cw*ee)/(2.*sw))+(-1.0/2.0)*(-(ee*sw)/(6.*cw))',
+              texname = 'AxialZdown')
+
+VectorAUp = Parameter(name = 'VectorAUp',
+              nature = 'internal',                      
+              type = 'real',
+              value = '(2*ee)/3.',
+              texname = 'VectorAUp')
+
+VectorADown = Parameter(name = 'VectorADown',
+              nature = 'internal',                        
+              type = 'real',
+              value = '-(ee)/3.',
+              texname = 'VectorADown')
+
+VectorWmDxU = Parameter(name = 'VectorWmDxU',
+              nature = 'internal',                        
+              type = 'real',
+              value = '(1.0/2.0)*((ee)/(sw*cmath.sqrt(2)))',
+              texname = 'VectorWmDxU')
+
+AxialWmDxU = Parameter(name = 'AxialWmDxU',
+              nature = 'internal',                        
+              type = 'real',
+              value = '(-1.0/2.0)*((ee)/(sw*cmath.sqrt(2)))',
+              texname = 'AxialWmDxU')
+
+VectorWpUxD = Parameter(name = 'VectorWpUxD',
+              nature = 'internal',                        
+              type = 'real',
+              value = '(1.0/2.0)*((ee)/(sw*cmath.sqrt(2)))',
+              texname = 'VectorWpUxD')
+
+AxialWpUxD = Parameter(name = 'AxialWpUxD',
+              nature = 'internal',                        
+              type = 'real',
+              value = '-(1.0/2.0)*((ee)/(sw*cmath.sqrt(2)))',
+              texname = 'AxialWpUxD')
+
+# == Additional parameters for goldstones
+
+CKM1x1 = Parameter(name = 'CKM1x1',
+                   nature = 'internal',
+                   type = 'complex',
+                   value = '1 - lamWS**2/2.',
+                   texname = '\\text{CKM1x1}')
+
+CKM1x2 = Parameter(name = 'CKM1x2',
+                   nature = 'internal',
+                   type = 'complex',
+                   value = 'lamWS',
+                   texname = '\\text{CKM1x2}')
+
+CKM1x3 = Parameter(name = 'CKM1x3',
+                   nature = 'internal',
+                   type = 'complex',
+                   value = 'AWS*lamWS**3*(-(etaWS*complex(0,1)) + rhoWS)',
+                   texname = '\\text{CKM1x3}')
+
+CKM2x1 = Parameter(name = 'CKM2x1',
+                   nature = 'internal',
+                   type = 'complex',
+                   value = '-lamWS',
+                   texname = '\\text{CKM2x1}')
+
+CKM2x2 = Parameter(name = 'CKM2x2',
+                   nature = 'internal',
+                   type = 'complex',
+                   value = '1 - lamWS**2/2.',
+                   texname = '\\text{CKM2x2}')
+
+CKM2x3 = Parameter(name = 'CKM2x3',
+                   nature = 'internal',
+                   type = 'complex',
+                   value = 'AWS*lamWS**2',
+                   texname = '\\text{CKM2x3}')
+
+CKM3x1 = Parameter(name = 'CKM3x1',
+                   nature = 'internal',
+                   type = 'complex',
+                   value = 'AWS*lamWS**3*(1 - etaWS*complex(0,1) - rhoWS)',
+                   texname = '\\text{CKM3x1}')
+
+CKM3x2 = Parameter(name = 'CKM3x2',
+                   nature = 'internal',
+                   type = 'complex',
+                   value = '-(AWS*lamWS**2)',
+                   texname = '\\text{CKM3x2}')
+
+CKM3x3 = Parameter(name = 'CKM3x3',
+                   nature = 'internal',
+                   type = 'complex',
+                   value = '1',
+                   texname = '\\text{CKM3x3}')
+
+I1x31 = Parameter(name = 'I1x31',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yb*complexconjugate(CKM1x3)',
+                  texname = '\\text{I1x31}')
+
+I1x32 = Parameter(name = 'I1x32',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yb*complexconjugate(CKM2x3)',
+                  texname = '\\text{I1x32}')
+
+I1x33 = Parameter(name = 'I1x33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yb*complexconjugate(CKM3x3)',
+                  texname = '\\text{I1x33}')
+
+I2x12 = Parameter(name = 'I2x12',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yc*complexconjugate(CKM2x1)',
+                  texname = '\\text{I2x12}')
+
+I2x13 = Parameter(name = 'I2x13',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yt*complexconjugate(CKM3x1)',
+                  texname = '\\text{I2x13}')
+
+I2x22 = Parameter(name = 'I2x22',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yc*complexconjugate(CKM2x2)',
+                  texname = '\\text{I2x22}')
+
+I2x23 = Parameter(name = 'I2x23',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yt*complexconjugate(CKM3x2)',
+                  texname = '\\text{I2x23}')
+
+I2x32 = Parameter(name = 'I2x32',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yc*complexconjugate(CKM2x3)',
+                  texname = '\\text{I2x32}')
+
+I2x33 = Parameter(name = 'I2x33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yt*complexconjugate(CKM3x3)',
+                  texname = '\\text{I2x33}')
+
+I3x21 = Parameter(name = 'I3x21',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'CKM2x1*yc',
+                  texname = '\\text{I3x21}')
+
+I3x22 = Parameter(name = 'I3x22',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'CKM2x2*yc',
+                  texname = '\\text{I3x22}')
+
+I3x23 = Parameter(name = 'I3x23',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'CKM2x3*yc',
+                  texname = '\\text{I3x23}')
+
+I3x31 = Parameter(name = 'I3x31',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'CKM3x1*yt',
+                  texname = '\\text{I3x31}')
+
+I3x32 = Parameter(name = 'I3x32',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'CKM3x2*yt',
+                  texname = '\\text{I3x32}')
+
+I3x33 = Parameter(name = 'I3x33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'CKM3x3*yt',
+                  texname = '\\text{I3x33}')
+
+I4x13 = Parameter(name = 'I4x13',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'CKM1x3*yb',
+                  texname = '\\text{I4x13}')
+
+I4x23 = Parameter(name = 'I4x23',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'CKM2x3*yb',
+                  texname = '\\text{I4x23}')
+
+I4x33 = Parameter(name = 'I4x33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'CKM3x3*yb',
+                  texname = '\\text{I4x33}')
+
+Vector_ubGp = Parameter(name = 'Vector_ubGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'I1x31',
+                  texname = '\\text{Vector_ubGp}')
+
+Vector_cdGp = Parameter(name = 'Vector_cdGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I2x12',
+                  texname = '\\text{Vector_cdGp}')
+
+Vector_csGp = Parameter(name = 'Vector_csGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I2x22',
+                  texname = '\\text{Vector_csGp}')
+
+Vector_cbGp = Parameter(name = 'Vector_cbGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'I1x32-I2x32',
+                  texname = '\\text{Vector_cbGp}')
+
+Vector_tdGp = Parameter(name = 'Vector_tdGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I2x13',
+                  texname = '\\text{Vector_tdGp}')
+
+Vector_tsGp = Parameter(name = 'Vector_tsGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I2x23',
+                  texname = '\\text{Vector_tsGp}')
+
+Vector_tbGp = Parameter(name = 'Vector_tbGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'I1x33-I2x33',
+                  texname = '\\text{Vector_tbGp}')
+
+Axial_ubGp = Parameter(name = 'Axial_ubGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I1x31',
+                  texname = '\\text{Axial_ubGp}')
+
+Axial_cdGp = Parameter(name = 'Axial_cdGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I2x12',
+                  texname = '\\text{Axial_cdGp}')
+
+Axial_csGp = Parameter(name = 'Axial_csGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I2x22',
+                  texname = '\\text{Axial_csGp}')
+
+Axial_cbGp = Parameter(name = 'Axial_cbGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I2x32-I1x32',
+                  texname = '\\text{Axial_cbGp}')
+
+Axial_tdGp = Parameter(name = 'Axial_tdGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I2x13',
+                  texname = '\\text{Axial_tdGp}')
+
+Axial_tsGp = Parameter(name = 'Axial_tsGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I2x23',
+                  texname = '\\text{Axial_tsGp}')
+
+Axial_tbGp = Parameter(name = 'Axial_tbGp',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I2x33-I1x33',
+                  texname = '\\text{Axial_tbGp}')
+
+Vector_ubGm = Parameter(name = 'Vector_ubGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I4x13',
+                  texname = '\\text{Vector_ubGm}')
+
+Vector_cdGm = Parameter(name = 'Vector_cdGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'I3x21',
+                  texname = '\\text{Vector_cdGm}')
+
+Vector_csGm = Parameter(name = 'Vector_csGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'I3x22',
+                  texname = '\\text{Vector_csGm}')
+
+Vector_cbGm = Parameter(name = 'Vector_cbGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'I3x23-I4x23',
+                  texname = '\\text{Vector_cbGm}')
+
+Vector_tdGm = Parameter(name = 'Vector_tdGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'I3x31',
+                  texname = '\\text{Vector_tdGm}')
+
+Vector_tsGm = Parameter(name = 'Vector_tsGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'I3x32',
+                  texname = '\\text{Vector_tsGm}')
+
+Vector_tbGm = Parameter(name = 'Vector_tbGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'I3x33-I4x33',
+                  texname = '\\text{Vector_tbGm}')
+
+Axial_ubGm = Parameter(name = 'Axial_ubGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I4x13',
+                  texname = '\\text{Axial_ubGm}')
+
+Axial_cdGm = Parameter(name = 'Axial_cdGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I3x21',
+                  texname = '\\text{Axial_cdGm}')
+
+Axial_csGm = Parameter(name = 'Axial_csGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I3x22',
+                  texname = '\\text{Axial_csGm}')
+
+Axial_cbGm = Parameter(name = 'Axial_cbGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I4x23-I3x23',
+                  texname = '\\text{Axial_cbGm}')
+
+Axial_tdGm = Parameter(name = 'Axial_tdGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I3x31',
+                  texname = '\\text{Axial_tdGm}')
+
+Axial_tsGm = Parameter(name = 'Axial_tsGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I3x32',
+                  texname = '\\text{Axial_tsGm}')
+
+Axial_tbGm = Parameter(name = 'Axial_tbGm',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = '-I4x33-I3x33',
+                  texname = '\\text{Axial_tbGm}')
diff --git a/Higgs/loop_sm_c3d4/particles.py b/Higgs/loop_sm_c3d4/particles.py
new file mode 100755
index 0000000000000000000000000000000000000000..8bf2b53fa0a193024806d1c31f0843cfc5b15fff
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/particles.py
@@ -0,0 +1,342 @@
+# This file was automatically created by FeynRules $Revision: 535 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Fri 18 Mar 2011 18:40:51
+
+
+from __future__ import division
+from object_library import all_particles, Particle
+import parameters as Param
+import CT_parameters as CTParam
+
+# ======================================================================
+# QCD particles
+# ======================================================================
+
+d = Particle(pdg_code = 1,
+             name = 'd',
+             antiname = 'd~',
+             spin = 2,
+             color = 3,
+             mass = Param.ZERO,
+             width = Param.ZERO,
+             texname = 'd',
+             antitexname = 'd',
+             charge = -1/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+d__tilde__ = d.anti()
+
+u = Particle(pdg_code = 2,
+             name = 'u',
+             antiname = 'u~',
+             spin = 2,
+             color = 3,
+             mass = Param.ZERO,
+             width = Param.ZERO,
+             texname = 'u',
+             antitexname = 'u',
+             charge = 2/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+u__tilde__ = u.anti()
+
+s = Particle(pdg_code = 3,
+             name = 's',
+             antiname = 's~',
+             spin = 2,
+             color = 3,
+             mass = Param.ZERO,
+             width = Param.ZERO,
+             texname = 's',
+             antitexname = 's',
+             charge = -1/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+s__tilde__ = s.anti()
+
+c = Particle(pdg_code = 4,
+             name = 'c',
+             antiname = 'c~',
+             spin = 2,
+             color = 3,
+             mass = Param.MC,
+             width = Param.ZERO,
+             texname = 'c',
+             antitexname = 'c',
+             charge = 2/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+c__tilde__ = c.anti()
+
+b = Particle(pdg_code = 5,
+             name = 'b',
+             antiname = 'b~',
+             spin = 2,
+             color = 3,
+             mass = Param.MB,
+             width = Param.ZERO,
+             texname = 'b',
+             antitexname = 'b',
+             charge = -1/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+b__tilde__ = b.anti()
+
+t = Particle(pdg_code = 6,
+             name = 't',
+             antiname = 't~',
+             spin = 2,
+             color = 3,
+             mass = Param.MT,
+             width = Param.WT,
+             texname = 't',
+             antitexname = 't',
+             charge = 2/3,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+t__tilde__ = t.anti()
+
+G = Particle(pdg_code = 21,
+             name = 'G',
+             antiname = 'G',
+             spin = 3,
+             color = 8,
+             mass = Param.ZERO,
+             width = Param.ZERO,
+             texname = 'G',
+             antitexname = 'G',
+             charge = 0,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+gh = Particle(pdg_code = 82,
+             name = 'gh',
+             antiname = 'gh~',
+             spin = -1,
+             color = 8,
+             mass = Param.ZERO,
+             width = Param.ZERO,
+             texname = 'gh',
+             antitexname = 'gh',
+             charge = 0,
+             LeptonNumber = 0,
+             line='dashed',
+             GhostNumber = 1)
+
+gh__tilde__ = gh.anti()
+
+# ======================================================================
+# Other particles 
+# ======================================================================
+
+ve = Particle(pdg_code = 12,
+              name = 've',
+              antiname = 've~',
+              spin = 2,
+              color = 1,
+              mass = Param.ZERO,
+              width = Param.ZERO,
+              texname = 've',
+              antitexname = 've',
+              charge = 0,
+              LeptonNumber = 1,
+              GhostNumber = 0)
+
+ve__tilde__ = ve.anti()
+
+vm = Particle(pdg_code = 14,
+              name = 'vm',
+              antiname = 'vm~',
+              spin = 2,
+              color = 1,
+              mass = Param.ZERO,
+              width = Param.ZERO,
+              texname = 'vm',
+              antitexname = 'vm',
+              charge = 0,
+              LeptonNumber = 1,
+              GhostNumber = 0)
+
+vm__tilde__ = vm.anti()
+
+vt = Particle(pdg_code = 16,
+              name = 'vt',
+              antiname = 'vt~',
+              spin = 2,
+              color = 1,
+              mass = Param.ZERO,
+              width = Param.ZERO,
+              texname = 'vt',
+              antitexname = 'vt',
+              charge = 0,
+              LeptonNumber = 1,
+              GhostNumber = 0)
+
+vt__tilde__ = vt.anti()
+
+A = Particle(pdg_code = 22,
+             name = 'A',
+             antiname = 'A',
+             spin = 3,
+             color = 1,
+             mass = Param.ZERO,
+             width = Param.ZERO,
+             texname = 'A',
+             antitexname = 'A',
+             charge = 0,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+Z = Particle(pdg_code = 23,
+             name = 'Z',
+             antiname = 'Z',
+             spin = 3,
+             color = 1,
+             mass = Param.MZ,
+             width = Param.WZ,
+             texname = 'Z',
+             antitexname = 'Z',
+             charge = 0,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+W__plus__ = Particle(pdg_code = 24,
+                     name = 'W+',
+                     antiname = 'W-',
+                     spin = 3,
+                     color = 1,
+                     mass = Param.MW,
+                     width = Param.WW,
+                     texname = 'W+',
+                     antitexname = 'W+',
+                     charge = 1,
+                     LeptonNumber = 0,
+                     GhostNumber = 0)
+
+W__minus__ = W__plus__.anti()
+
+H = Particle(pdg_code = 25,
+             name = 'H',
+             antiname = 'H',
+             spin = 1,
+             color = 1,
+             mass = Param.MH,
+             width = Param.WH,
+             texname = '\\phi',
+             antitexname = '\\phi',
+             charge = 0,
+             LeptonNumber = 0,
+             GhostNumber = 0)
+
+e__minus__ = Particle(pdg_code = 11,
+                      name = 'e-',
+                      antiname = 'e+',
+                      spin = 2,
+                      color = 1,
+                      mass = Param.Me,
+                      width = Param.ZERO,
+                      texname = 'e-',
+                      antitexname = 'e-',
+                      charge = -1,
+                      LeptonNumber = 1,
+                      GhostNumber = 0)
+
+e__plus__ = e__minus__.anti()
+
+m__minus__ = Particle(pdg_code = 13,
+                      name = 'm-',
+                      antiname = 'm+',
+                      spin = 2,
+                      color = 1,
+                      mass = Param.MM,
+                      width = Param.ZERO,
+                      texname = 'm-',
+                      antitexname = 'm-',
+                      charge = -1,
+                      LeptonNumber = 1,
+                      GhostNumber = 0)
+
+m__plus__ = m__minus__.anti()
+
+tt__minus__ = Particle(pdg_code = 15,
+                       name = 'tt-',
+                       antiname = 'tt+',
+                       spin = 2,
+                       color = 1,
+                       mass = Param.MTA,
+                       width = Param.WTau,
+                       texname = 'tt-',
+                       antitexname = 'tt-',
+                       charge = -1,
+                       LeptonNumber = 1,
+                       GhostNumber = 0)
+
+tt__plus__ = tt__minus__.anti()
+
+# ======================================================================
+# Goldstones 
+# ======================================================================
+
+G0 = Particle(pdg_code = 250,
+              name = 'G0',
+              antiname = 'G0',
+              spin = 1,
+              color = 1,
+              mass = Param.MZ,
+              width = Param.WZ,
+              texname = 'G0',
+              antitexname = 'G0',
+              GoldstoneBoson = True,
+              charge = 0,
+              GhostNumber = 0,
+              LeptonNumber = 0)
+
+G__plus__ = Particle(pdg_code = 251,
+                     name = 'G+',
+                     antiname = 'G-',
+                     spin = 1,
+                     color = 1,
+                     mass = Param.MW,
+                     width = Param.WW,
+                     texname = 'G+',
+                     antitexname = 'G-',
+                     GoldstoneBoson = True,
+                     charge = 1,
+                     GhostNumber = 0,
+                     LeptonNumber = 0)
+
+G__minus__ = G__plus__.anti()
+
+# Wavefunction renormalization
+
+b.loop_particles = [[[5,21]]]
+b.counterterm = {(1,0,0):CTParam.bWcft_UV.value}
+
+c.loop_particles = [[[4,21]]]
+c.counterterm = {(1,0,0):CTParam.cWcft_UV.value}
+
+t.loop_particles = [[[6,21]]]
+t.counterterm = {(1,0,0):CTParam.tWcft_UV.value}
+
+G.loop_particles = [[[4]],[[5]],[[6]]]
+G.counterterm = {(1,0,0):CTParam.GWcft_UV_c.value,(1,0,1):CTParam.GWcft_UV_b.value,(1,0,2):CTParam.GWcft_UV_t.value}
+
+# Set counterterms values
+
+#Param.MB.loop_particles= [[[5,21]]]
+#Param.MB.counterterm = {(1,0,0):CTParam.bMass_UV.value}
+
+#Param.MC.loop_particles= [[[4,21]]]
+#Param.MC.counterterm = {(1,0,0):CTParam.cMass_UV.value}
+
+#Param.MT.loop_particles= [[[6,21]]]
+#Param.MT.counterterm = {(1,0,0):CTParam.tMass_UV.value}
+
+#Param.G.loop_particles = [[[2],[1],[3]],[[4]],[[5]],[[6]],[[21]]],
+#Param.G.counterterm = {(1,0,0):CTParam.G_UVq.value,(1,0,1):CTParam.G_UVc.value,(1,0,2):CTParam.G_UVb.value,(1,0,3):CTParam.G_UVt.value,(1,0,4):CTParam.G_UVg.value},
diff --git a/Higgs/loop_sm_c3d4/restrict_c_mass.dat b/Higgs/loop_sm_c3d4/restrict_c_mass.dat
new file mode 100644
index 0000000000000000000000000000000000000000..063853148c39057e482794a212a668e26e4a7e36
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/restrict_c_mass.dat
@@ -0,0 +1,59 @@
+######################################################################
+## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################
+######################################################################
+
+###################################
+### INFORMATION FOR LOOP
+####################################
+Block loop
+  1 9.118800e+01 # MU_R
+
+###################################
+## INFORMATION FOR SMINPUTS
+###################################
+Block SMINPUTS 
+    1 1.325070e+02 # aEWM1 
+    2 1.166390e-05 # Gf 
+    3 1.180000e-01 # aS 
+
+###################################
+## INFORMATION FOR MASS
+###################################
+Block MASS 
+    4 1.550000e+00 # MC 
+    5 4.700000e+00 # MB 
+    6 1.730000e+02 # MT 
+   11 0.000000e+00 # Me 
+   13 0.000000e+00 # MM 
+   15 1.777000e+00 # MTA 
+   23 9.118800e+01 # MZ 
+   25 1.250000e+02 # MH 
+
+###################################
+## INFORMATION FOR DECAY
+###################################
+DECAY   6 1.491500E+00
+DECAY  15 0.000000e+00
+DECAY  23 2.441404e+00 
+DECAY  24 2.047600e+00 
+DECAY  25 6.38233934e-03 
+
+###################################
+## INFORMATION FOR WOLFENSTEIN
+###################################
+Block Wolfenstein 
+    1 0.000000e+00 # lamWS 
+    2 0.000000e+00 # AWS 
+    3 0.000000e+00 # rhoWS 
+    4 0.000000e+00 # etaWS 
+
+###################################
+## INFORMATION FOR YUKAWA
+###################################
+Block YUKAWA 
+    4 1.550000e+00 # ymc 
+    5 4.700000e+00 # ymb 
+    6 1.730000e+02 # ymt 
+   11 0.000000e+00 # yme 
+   13 0.000000e+00 # ymm 
+   15 1.777000e+00 # ymtau 
diff --git a/Higgs/loop_sm_c3d4/restrict_ckm.dat b/Higgs/loop_sm_c3d4/restrict_ckm.dat
new file mode 100644
index 0000000000000000000000000000000000000000..464c9a7630a343dc5debdcb064463a2f5dbfec43
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/restrict_ckm.dat
@@ -0,0 +1,59 @@
+######################################################################
+## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################
+######################################################################
+
+###################################
+### INFORMATION FOR LOOP
+####################################
+Block loop
+  1 9.118800e+01 # MU_R
+
+###################################
+## INFORMATION FOR SMINPUTS
+###################################
+Block SMINPUTS 
+    1 1.325070e+02 # aEWM1 
+    2 1.166390e-05 # Gf 
+    3 1.180000e-01 # aS 
+
+###################################
+## INFORMATION FOR MASS
+###################################
+Block MASS 
+    4 0.000000e+00 # MC 
+    5 4.700000e+00 # MB 
+    6 1.730000e+02 # MT 
+   11 0.000000e+00 # Me 
+   13 0.000000e+00 # MM 
+   15 1.777000e+00 # MTA 
+   23 9.118800e+01 # MZ 
+   25 1.250000e+02 # MH 
+
+###################################
+## INFORMATION FOR DECAY
+###################################
+DECAY   6 1.491500E+00
+DECAY  15 0.000000e+00
+DECAY  23 2.441404e+00 
+DECAY  24 2.047600e+00 
+DECAY  25 6.38233934e-03 
+
+###################################
+## INFORMATION FOR WOLFENSTEIN
+###################################
+Block Wolfenstein 
+    1 2.253000e-01 # lamWS 
+    2 8.080000e-01 # AWS 
+    3 1.320000e-01 # rhoWS 
+    4 3.410000e-01 # etaWS 
+
+###################################
+## INFORMATION FOR YUKAWA
+###################################
+Block YUKAWA 
+    4 0.000000e+00 # ymc 
+    5 4.700000e+00 # ymb 
+    6 1.730000e+02 # ymt 
+   11 0.000000e+00 # yme 
+   13 0.000000e+00 # ymm 
+   15 1.777000e+00 # ymtau 
diff --git a/Higgs/loop_sm_c3d4/restrict_default.dat b/Higgs/loop_sm_c3d4/restrict_default.dat
new file mode 100644
index 0000000000000000000000000000000000000000..6ac783615c7e3e63336b0b598447909350e0e4e0
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/restrict_default.dat
@@ -0,0 +1,65 @@
+######################################################################
+## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################
+######################################################################
+
+###################################
+### INFORMATION FOR LOOP
+####################################
+Block loop
+  1 9.118800e+01 # MU_R
+
+###################################
+## INFORMATION FOR SMINPUTS
+###################################
+Block SMINPUTS 
+    1 1.325070e+02 # aEWM1 
+    2 1.166390e-05 # Gf 
+    3 1.180000e-01 # aS
+
+Block tripcoup
+    4 2.0000000000 # c3
+
+Block quartcoup
+    6 3.0000000000 # d4
+
+###################################
+## INFORMATION FOR MASS
+###################################
+Block MASS 
+    4 0.000000e+00 # MC 
+    5 4.700000e+00 # MB 
+    6 1.730000e+02 # MT 
+   11 0.000000e+00 # Me 
+   13 0.000000e+00 # MM 
+   15 1.777000e+00 # MTA 
+   23 9.118800e+01 # MZ 
+   25 1.250000e+02 # MH 
+
+###################################
+## INFORMATION FOR DECAY
+###################################
+DECAY   6 1.491500E+00
+DECAY  15 0.000000e+00
+DECAY  23 2.441404e+00 
+DECAY  24 2.047600e+00 
+DECAY  25 6.38233934e-03
+
+###################################
+## INFORMATION FOR WOLFENSTEIN
+###################################
+Block Wolfenstein 
+    1 0.000000e+00 # lamWS 
+    2 0.000000e+00 # AWS 
+    3 0.000000e+00 # rhoWS 
+    4 0.000000e+00 # etaWS 
+
+###################################
+## INFORMATION FOR YUKAWA
+###################################
+Block YUKAWA 
+    4 0.000000e+00 # ymc 
+    5 4.700000e+00 # ymb 
+    6 1.730000e+02 # ymt 
+   11 0.000000e+00 # yme 
+   13 0.000000e+00 # ymm 
+   15 1.777000e+00 # ymtau 
diff --git a/Higgs/loop_sm_c3d4/restrict_lepton_masses.dat b/Higgs/loop_sm_c3d4/restrict_lepton_masses.dat
new file mode 100644
index 0000000000000000000000000000000000000000..7512d061b872d8a18c4f5217540690eb87d41925
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/restrict_lepton_masses.dat
@@ -0,0 +1,59 @@
+######################################################################
+## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################
+######################################################################
+
+###################################
+### INFORMATION FOR LOOP
+####################################
+Block loop
+  1 9.118800e+01 # MU_R
+
+###################################
+## INFORMATION FOR SMINPUTS
+###################################
+Block SMINPUTS 
+    1 1.325070e+02 # aEWM1 
+    2 1.166390e-05 # Gf 
+    3 1.180000e-01 # aS 
+
+###################################
+## INFORMATION FOR MASS
+###################################
+Block MASS 
+    4 0.000000e+00 # MC 
+    5 4.700000e+00 # MB 
+    6 1.730000e+02 # MT 
+   11 5.110000e-04 # Me 
+   13 1.056600e-01 # MM 
+   15 1.777000e+00 # MTA 
+   23 9.118800e+01 # MZ 
+   25 1.250000e+02 # MH 
+
+###################################
+## INFORMATION FOR DECAY
+###################################
+DECAY   6 1.491500E+00
+DECAY  15 2.270000e-12
+DECAY  23 2.441404e+00 
+DECAY  24 2.047600e+00 
+DECAY  25 6.38233934e-03 
+
+###################################
+## INFORMATION FOR WOLFENSTEIN
+###################################
+Block Wolfenstein 
+    1 0.000000e+00 # lamWS 
+    2 0.000000e+00 # AWS 
+    3 0.000000e+00 # rhoWS 
+    4 0.000000e+00 # etaWS 
+
+###################################
+## INFORMATION FOR YUKAWA
+###################################
+Block YUKAWA 
+    4 0.000000e+00 # ymc 
+    5 4.700000e+00 # ymb
+    6 1.730000e+02 # ymt
+   11 5.110000e-04 # yme
+   13 1.056600e-01 # ymm
+   15 1.777000e+00 # ymtau 
diff --git a/Higgs/loop_sm_c3d4/restrict_no_b_mass.dat b/Higgs/loop_sm_c3d4/restrict_no_b_mass.dat
new file mode 100644
index 0000000000000000000000000000000000000000..a23cd75616e68b0921fa65477602c3e18cd2793c
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/restrict_no_b_mass.dat
@@ -0,0 +1,59 @@
+######################################################################
+## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################
+######################################################################
+
+###################################
+### INFORMATION FOR LOOP
+####################################
+Block loop
+  1 9.118800e+01 # MU_R
+
+###################################
+## INFORMATION FOR SMINPUTS
+###################################
+Block SMINPUTS 
+    1 1.325070e+02 # aEWM1 
+    2 1.166390e-05 # Gf 
+    3 1.180000e-01 # aS 
+
+###################################
+## INFORMATION FOR MASS
+###################################
+Block MASS 
+    4 0.000000e+00 # MC 
+    5 0.000000e+00 # MB 
+    6 1.730000e+02 # MT 
+   11 0.000000e+00 # Me 
+   13 0.000000e+00 # MM 
+   15 1.777000e+00 # MTA 
+   23 9.118800e+01 # MZ 
+   25 1.250000e+02 # MH 
+
+###################################
+## INFORMATION FOR DECAY
+###################################
+DECAY   6 1.491500E+00
+DECAY  15 0.000000e+00
+DECAY  23 2.441404e+00 
+DECAY  24 2.047600e+00 
+DECAY  25 6.38233934e-03 
+
+###################################
+## INFORMATION FOR WOLFENSTEIN
+###################################
+Block Wolfenstein 
+    1 0.000000e+00 # lamWS 
+    2 0.000000e+00 # AWS 
+    3 0.000000e+00 # rhoWS 
+    4 0.000000e+00 # etaWS 
+
+###################################
+## INFORMATION FOR YUKAWA
+###################################
+Block YUKAWA 
+    4 0.000000e+00 # ymc 
+    5 0.000000e+00 # ymb 
+    6 1.730000e+02 # ymt 
+   11 0.000000e+00 # yme 
+   13 0.000000e+00 # ymm 
+   15 1.777000e+00 # ymtau 
diff --git a/Higgs/loop_sm_c3d4/restrict_no_masses.dat b/Higgs/loop_sm_c3d4/restrict_no_masses.dat
new file mode 100644
index 0000000000000000000000000000000000000000..700612593eda124b4a8896e9cdef609474d12c25
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/restrict_no_masses.dat
@@ -0,0 +1,59 @@
+######################################################################
+## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################
+######################################################################
+
+###################################
+### INFORMATION FOR LOOP
+####################################
+Block loop
+  1 9.118800e+01 # MU_R
+
+###################################
+## INFORMATION FOR SMINPUTS
+###################################
+Block SMINPUTS 
+    1 1.325070e+02 # aEWM1 
+    2 1.166390e-05 # Gf 
+    3 1.180000e-01 # aS 
+
+###################################
+## INFORMATION FOR MASS
+###################################
+Block MASS 
+    4 0.000000e+00 # MC 
+    5 0.000000e+00 # MB 
+    6 1.730000e+02 # MT 
+   11 0.000000e+00 # Me 
+   13 0.000000e+00 # MM 
+   15 0.000000e+00 # MTA 
+   23 9.118800e+01 # MZ 
+   25 1.250000e+02 # MH 
+
+###################################
+## INFORMATION FOR DECAY
+###################################
+DECAY   6 1.491500E+00
+DECAY  15 0.000000e+00
+DECAY  23 2.441404e+00 
+DECAY  24 2.047600e+00 
+DECAY  25 6.38233934e-03 
+
+###################################
+## INFORMATION FOR WOLFENSTEIN
+###################################
+Block Wolfenstein 
+    1 0.000000e+00 # lamWS 
+    2 0.000000e+00 # AWS 
+    3 0.000000e+00 # rhoWS 
+    4 0.000000e+00 # etaWS 
+
+###################################
+## INFORMATION FOR YUKAWA
+###################################
+Block YUKAWA 
+    4 0.000000e+00 # ymc 
+    5 0.000000e+00 # ymb 
+    6 1.730000e+02 # ymt 
+   11 0.000000e+00 # yme 
+   13 0.000000e+00 # ymm 
+   15 0.000000e+00 # ymtau 
diff --git a/Higgs/loop_sm_c3d4/restrict_no_tau_mass.dat b/Higgs/loop_sm_c3d4/restrict_no_tau_mass.dat
new file mode 100644
index 0000000000000000000000000000000000000000..9626186db64f4bf2e30c9e55758d193aebe89298
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/restrict_no_tau_mass.dat
@@ -0,0 +1,59 @@
+######################################################################
+## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################
+######################################################################
+
+###################################
+### INFORMATION FOR LOOP
+####################################
+Block loop
+  1 9.118800e+01 # MU_R
+
+###################################
+## INFORMATION FOR SMINPUTS
+###################################
+Block SMINPUTS 
+    1 1.325070e+02 # aEWM1 
+    2 1.166390e-05 # Gf 
+    3 1.180000e-01 # aS 
+
+###################################
+## INFORMATION FOR MASS
+###################################
+Block MASS 
+    4 0.000000e+00 # MC 
+    5 4.700000e+00 # MB 
+    6 1.730000e+02 # MT 
+   11 0.000000e+00 # Me 
+   13 0.000000e+00 # MM 
+   15 0.000000e+00 # MTA 
+   23 9.118800e+01 # MZ 
+   25 1.250000e+02 # MH 
+
+###################################
+## INFORMATION FOR DECAY
+###################################
+DECAY   6 1.491500E+00
+DECAY  15 0.000000e+00
+DECAY  23 2.441404e+00 
+DECAY  24 2.047600e+00 
+DECAY  25 6.38233934e-03 
+
+###################################
+## INFORMATION FOR WOLFENSTEIN
+###################################
+Block Wolfenstein 
+    1 0.000000e+00 # lamWS 
+    2 0.000000e+00 # AWS 
+    3 0.000000e+00 # rhoWS 
+    4 0.000000e+00 # etaWS 
+
+###################################
+## INFORMATION FOR YUKAWA
+###################################
+Block YUKAWA 
+    4 0.000000e+00 # ymc 
+    5 4.700000e+00 # ymb 
+    6 1.730000e+02 # ymt 
+   11 0.000000e+00 # yme 
+   13 0.000000e+00 # ymm 
+   15 0.000000e+00 # ymtau 
diff --git a/Higgs/loop_sm_c3d4/restrict_no_widths.dat b/Higgs/loop_sm_c3d4/restrict_no_widths.dat
new file mode 100644
index 0000000000000000000000000000000000000000..1b9aa0329544d721fc36b7e6d6d208008616363e
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/restrict_no_widths.dat
@@ -0,0 +1,59 @@
+######################################################################
+## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################
+######################################################################
+
+###################################
+### INFORMATION FOR LOOP
+####################################
+Block loop
+  1 9.118800e+01 # MU_R
+
+###################################
+## INFORMATION FOR SMINPUTS
+###################################
+Block SMINPUTS 
+    1 1.325070e+02 # aEWM1 
+    2 1.166390e-05 # Gf 
+    3 1.180000e-01 # aS 
+
+###################################
+## INFORMATION FOR MASS
+###################################
+Block MASS 
+    4 0.000000e+00 # MC 
+    5 4.700000e+00 # MB 
+    6 1.730000e+02 # MT 
+   11 0.000000e+00 # Me 
+   13 0.000000e+00 # MM 
+   15 1.777000e+00 # MTA 
+   23 9.118800e+01 # MZ 
+   25 1.250000e+02 # MH 
+
+###################################
+## INFORMATION FOR DECAY
+###################################
+DECAY  6 0.000000e+00
+DECAY  15 0.000000e+00
+DECAY  23 0.000000e+00
+DECAY  24 0.000000e+00
+DECAY  25 0.000000e+00
+
+###################################
+## INFORMATION FOR WOLFENSTEIN
+###################################
+Block Wolfenstein 
+    1 0.000000e+00 # lamWS 
+    2 0.000000e+00 # AWS 
+    3 0.000000e+00 # rhoWS 
+    4 0.000000e+00 # etaWS 
+
+###################################
+## INFORMATION FOR YUKAWA
+###################################
+Block YUKAWA 
+    4 0.000000e+00 # ymc 
+    5 4.700000e+00 # ymb 
+    6 1.730000e+02 # ymt 
+   11 0.000000e+00 # yme 
+   13 0.000000e+00 # ymm 
+   15 1.777000e+00 # ymtau 
diff --git a/Higgs/loop_sm_c3d4/restrict_parallel_test.dat b/Higgs/loop_sm_c3d4/restrict_parallel_test.dat
new file mode 100644
index 0000000000000000000000000000000000000000..185029d34c995efd5ae353eadbf1a124e30f1ca9
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/restrict_parallel_test.dat
@@ -0,0 +1,59 @@
+######################################################################
+## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################
+######################################################################
+
+###################################
+### INFORMATION FOR LOOP
+####################################
+Block loop
+  1 9.118800e+01 # MU_R
+
+###################################
+## INFORMATION FOR SMINPUTS
+###################################
+Block SMINPUTS 
+    1 1.32506980e+02 # aEWM1 
+    2 1.16639000e-05 # Gf 
+    3 1.18000000e-01 # aS 
+
+###################################
+## INFORMATION FOR MASS
+###################################
+Block MASS 
+    4 0.000000e+00 # MC 
+    5 4.620000e+00 # MB 
+    6 1.743000e+02 # MT 
+   11 0.000000e+00 # Me 
+   13 0.000000e+00 # MM 
+   15 1.777000e+00 # MTA 
+   23 9.118800e+01 # MZ 
+   25 1.200000e+02 # MH 
+
+###################################
+## INFORMATION FOR DECAY
+###################################
+DECAY   6 0.000000e+00
+DECAY  15 0.000000e+00
+DECAY  23 0.000000e+00 
+DECAY  24 0.000000e+00 
+DECAY  25 0.000000e+00 
+
+###################################
+## INFORMATION FOR WOLFENSTEIN
+###################################
+Block Wolfenstein 
+    1 0.000000e+00 # lamWS 
+    2 0.000000e+00 # AWS 
+    3 0.000000e+00 # rhoWS 
+    4 0.000000e+00 # etaWS 
+
+###################################
+## INFORMATION FOR YUKAWA
+###################################
+Block YUKAWA 
+    4 0.000000e+00 # ymc 
+    5 4.6200000e+00 # ymb 
+    6 1.743000e+02 # ymt 
+   11 0.000000e+00 # yme 
+   13 0.000000e+00 # ymm 
+   15 1.777000e+00 # ymtau 
diff --git a/Higgs/loop_sm_c3d4/restrict_test.dat b/Higgs/loop_sm_c3d4/restrict_test.dat
new file mode 100644
index 0000000000000000000000000000000000000000..4f33591caceb089f53df245247ffba147f3ebe72
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/restrict_test.dat
@@ -0,0 +1,61 @@
+######################################################################
+## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################
+######################################################################
+
+###################################
+### INFORMATION FOR LOOP
+####################################
+Block loop
+  1 9.118800e+01 # MU_R
+
+###################################
+## INFORMATION FOR SMINPUTS
+###################################
+Block SMINPUTS 
+    1 1.325070e+02 # aEWM1 
+    2 1.166390e-05 # Gf 
+    3 1.180000e-01 # aS 
+
+###################################
+## INFORMATION FOR MASS
+###################################
+Block MASS 
+    4 0.000000e+00 # MC 
+    5 4.700000e+00 # MB 
+    6 1.730000e+02 # MT 
+   11 0.000000e+00 # Me 
+   13 0.000000e+00 # MM 
+   15 1.777000e+00 # MTA 
+   23 9.118800e+01 # MZ 
+   25 1.200000e+02 # MH 
+
+###################################
+## INFORMATION FOR DECAY
+###################################
+DECAY   6 0.000000E+00
+DECAY  15 0.000000e+00
+DECAY  23 0.000000e+00 
+DECAY  24 0.000000e+00 
+DECAY  25 0.000000e+00 
+
+###################################
+## INFORMATION FOR WOLFENSTEIN
+###################################
+Block Wolfenstein 
+    1 0.000000e+00 # lamWS 
+    2 0.000000e+00 # AWS 
+    3 0.000000e+00 # rhoWS 
+    4 0.000000e+00 # etaWS 
+
+###################################
+## INFORMATION FOR YUKAWA
+###################################
+Block YUKAWA 
+    4 0.000000e+00 # ymc 
+#    5 4.200000e+00 # ymb 
+#    6 1.645000e+02 # ymt 
+    5 0.000000e+00 # ymb 
+    6 0.00000e+00 # ymt 
+   11 0.000000e+00 # yme 
+   13 0.000000e+00 # ymm 
+   15 1.777000e+00 # ymtau 
diff --git a/Higgs/loop_sm_c3d4/restrict_zeromass_ckm.dat b/Higgs/loop_sm_c3d4/restrict_zeromass_ckm.dat
new file mode 100644
index 0000000000000000000000000000000000000000..6846e8975447cc11ed5f440bcc9e49cf182c11d6
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/restrict_zeromass_ckm.dat
@@ -0,0 +1,59 @@
+######################################################################
+## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################
+######################################################################
+
+###################################
+### INFORMATION FOR LOOP
+####################################
+Block loop
+  1 9.118800e+01 # MU_R
+
+###################################
+## INFORMATION FOR SMINPUTS
+###################################
+Block SMINPUTS 
+    1 1.325070e+02 # aEWM1 
+    2 1.166390e-05 # Gf 
+    3 1.180000e-01 # aS 
+
+###################################
+## INFORMATION FOR MASS
+###################################
+Block MASS 
+    4 0.000000e+00 # MC 
+    5 0.000000e+00 # MB 
+    6 1.730000e+02 # MT 
+   11 0.000000e+00 # Me 
+   13 0.000000e+00 # MM 
+   15 0.000000e+00 # MTA 
+   23 9.118800e+01 # MZ 
+   25 1.250000e+02 # MH 
+
+###################################
+## INFORMATION FOR DECAY
+###################################
+DECAY   6 1.491500E+00
+DECAY  15 0.000000e+00
+DECAY  23 2.441404e+00 
+DECAY  24 2.047600e+00 
+DECAY  25 6.38233934e-03 
+
+###################################
+## INFORMATION FOR WOLFENSTEIN
+###################################
+Block Wolfenstein 
+    1 2.253000e-01 # lamWS 
+    2 8.080000e-01 # AWS 
+    3 1.320000e-01 # rhoWS 
+    4 3.410000e-01 # etaWS 
+
+###################################
+## INFORMATION FOR YUKAWA
+###################################
+Block YUKAWA 
+    4 0.000000e+00 # ymc 
+    5 0.000000e+00 # ymb 
+    6 1.730000e+02 # ymt 
+   11 0.000000e+00 # yme 
+   13 0.000000e+00 # ymm 
+   15 0.000000e+00 # ymtau 
diff --git a/Higgs/loop_sm_c3d4/vertices.py b/Higgs/loop_sm_c3d4/vertices.py
new file mode 100755
index 0000000000000000000000000000000000000000..49494a06ad0f2e07d8e8c9b1d402df87d30bcae3
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/vertices.py
@@ -0,0 +1,802 @@
+# This file was automatically created by FeynRules $Revision: 535 $
+# Mathematica version: 7.0 for Mac OS X x86 (64-bit) (November 11, 2008)
+# Date: Fri 18 Mar 2011 18:40:51
+
+
+from object_library import all_vertices, all_CTvertices, Vertex, CTVertex
+import particles as P
+import couplings as C
+import lorentz as L
+
+# ======================================================================
+# QCD base vertices
+# ======================================================================
+
+V_3 = Vertex(name = 'V_3',
+              particles = [ P.G, P.G, P.G ],
+              color = [ 'f(1,2,3)' ],
+              lorentz = [ L.VVV1 ],
+              couplings = {(0,0):C.GC_4})
+              
+V_4 = Vertex(name = 'V_4',
+              particles = [ P.G, P.G, P.G, P.G ],
+              color = [ 'f(-1,1,2)*f(3,4,-1)', 'f(-1,1,3)*f(2,4,-1)', 'f(-1,1,4)*f(2,3,-1)' ],
+              lorentz = [ L.VVVV1, L.VVVV3, L.VVVV4 ],
+              couplings = {(1,1):C.GC_6,(0,0):C.GC_6,(2,2):C.GC_6})
+
+V_24 = Vertex(name = 'V_24',
+              particles = [ P.d__tilde__, P.d, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_5})
+
+V_49 = Vertex(name = 'V_49',
+               particles = [ P.u__tilde__, P.u, P.G ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_5})
+
+V_25 = Vertex(name = 'V_25',
+              particles = [ P.s__tilde__, P.s, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_5})
+
+V_50 = Vertex(name = 'V_50',
+               particles = [ P.c__tilde__, P.c, P.G ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_5})
+
+V_26 = Vertex(name = 'V_26',
+              particles = [ P.b__tilde__, P.b, P.G ],
+              color = [ 'T(3,2,1)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_5})
+
+V_51 = Vertex(name = 'V_51',
+               particles = [ P.t__tilde__, P.t, P.G ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_5})
+
+# QCD ghost
+V_Ggh = Vertex(name = 'V_Ggh',
+               particles = [ P.gh__tilde__, P.gh, P.G ],
+               color = [ 'f(2,3,1)' ],
+               lorentz = [ L.GHGHG ],
+               couplings = {(0,0):C.GC_4})
+
+# ======================================================================
+# Non-QCD base vertices
+# ======================================================================
+
+V_HHHH = Vertex(name = 'V_HHHH',
+             particles = [ P.H, P.H, P.H, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSSS1 ],
+             couplings = {(0,0):C.GC_HHHH})
+
+V_1 = Vertex(name = 'V_1',
+             particles = [ P.H, P.H, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSS1 ],
+             couplings = {(0,0):C.GC_30})
+
+V_5 = Vertex(name = 'V_5',
+             particles = [ P.A, P.W__minus__, P.W__plus__ ],
+             color = [ '1' ],
+             lorentz = [ L.VVV1 ],
+             couplings = {(0,0):C.GC_25})
+
+V_6 = Vertex(name = 'V_6',
+             particles = [ P.W__minus__, P.W__plus__, P.H, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.VVSS1 ],
+             couplings = {(0,0):C.GC_10})
+
+V_7 = Vertex(name = 'V_7',
+             particles = [ P.W__minus__, P.W__plus__, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.VVS1 ],
+             couplings = {(0,0):C.GC_31})
+
+V_8 = Vertex(name = 'V_8',
+             particles = [ P.A, P.A, P.W__minus__, P.W__plus__ ],
+             color = [ '1' ],
+             lorentz = [ L.VVVV2 ],
+             couplings = {(0,0):C.GC_27})
+
+V_9 = Vertex(name = 'V_9',
+             particles = [ P.W__minus__, P.W__plus__, P.Z ],
+             color = [ '1' ],
+             lorentz = [ L.VVV1 ],
+             couplings = {(0,0):C.GC_7})
+
+V_10 = Vertex(name = 'V_10',
+              particles = [ P.W__minus__, P.W__minus__, P.W__plus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVVV2 ],
+              couplings = {(0,0):C.GC_8})
+
+V_11 = Vertex(name = 'V_11',
+              particles = [ P.A, P.W__minus__, P.W__plus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.VVVV5 ],
+              couplings = {(0,0):C.GC_26})
+
+V_12 = Vertex(name = 'V_12',
+              particles = [ P.Z, P.Z, P.H, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_29})
+
+V_13 = Vertex(name = 'V_13',
+              particles = [ P.Z, P.Z, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_32})
+
+V_14 = Vertex(name = 'V_14',
+              particles = [ P.W__minus__, P.W__plus__, P.Z, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.VVVV2 ],
+              couplings = {(0,0):C.GC_9})
+
+V_15 = Vertex(name = 'V_15',
+              particles = [ P.d__tilde__, P.d, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_1})
+
+V_16 = Vertex(name = 'V_16',
+              particles = [ P.s__tilde__, P.s, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_1})
+
+V_17 = Vertex(name = 'V_17',
+              particles = [ P.b__tilde__, P.b, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_1})
+
+V_18 = Vertex(name = 'V_18',
+              particles = [ P.e__plus__, P.e__minus__, P.A ],
+              color = [ '1' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_19 = Vertex(name = 'V_19',
+              particles = [ P.m__plus__, P.m__minus__, P.A ],
+              color = [ '1' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_20 = Vertex(name = 'V_20',
+              particles = [ P.tt__plus__, P.tt__minus__, P.A ],
+              color = [ '1' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_21 = Vertex(name = 'V_21',
+              particles = [ P.u__tilde__, P.u, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_2})
+
+V_22 = Vertex(name = 'V_22',
+              particles = [ P.c__tilde__, P.c, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_2})
+
+V_23 = Vertex(name = 'V_23',
+              particles = [ P.t__tilde__, P.t, P.A ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_2})
+
+V_27 = Vertex(name = 'V_27',
+              particles = [ P.b__tilde__, P.b, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1 ],
+              couplings = {(0,0):C.GC_33})
+
+V_28 = Vertex(name = 'V_28',
+              particles = [ P.d__tilde__, P.d, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV3 ],
+              couplings = {(0,0):C.GC_21,(0,1):C.GC_23})
+
+V_29 = Vertex(name = 'V_29',
+              particles = [ P.s__tilde__, P.s, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV3 ],
+              couplings = {(0,0):C.GC_21,(0,1):C.GC_23})
+
+V_30 = Vertex(name = 'V_30',
+              particles = [ P.b__tilde__, P.b, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV3 ],
+              couplings = {(0,0):C.GC_21,(0,1):C.GC_23})
+
+V_31 = Vertex(name = 'V_31',
+              particles = [ P.d__tilde__, P.u, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_12})
+
+V_32 = Vertex(name = 'V_32',
+              particles = [ P.d__tilde__, P.c, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_15})
+
+V_33 = Vertex(name = 'V_33',
+              particles = [ P.d__tilde__, P.t, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_18})
+
+V_34 = Vertex(name = 'V_34',
+              particles = [ P.s__tilde__, P.u, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_13})
+
+V_35 = Vertex(name = 'V_35',
+              particles = [ P.s__tilde__, P.c, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_16})
+
+V_36 = Vertex(name = 'V_36',
+              particles = [ P.s__tilde__, P.t, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_19})
+
+V_37 = Vertex(name = 'V_37',
+              particles = [ P.b__tilde__, P.u, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_14})
+
+V_38 = Vertex(name = 'V_38',
+              particles = [ P.b__tilde__, P.c, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_17})
+
+V_39 = Vertex(name = 'V_39',
+              particles = [ P.b__tilde__, P.t, P.W__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_20})
+
+V_40 = Vertex(name = 'V_40',
+              particles = [ P.u__tilde__, P.d, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_39})
+
+V_41 = Vertex(name = 'V_41',
+              particles = [ P.c__tilde__, P.d, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_42})
+
+V_42 = Vertex(name = 'V_42',
+              particles = [ P.t__tilde__, P.d, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_45})
+
+V_43 = Vertex(name = 'V_43',
+              particles = [ P.u__tilde__, P.s, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_40})
+
+V_44 = Vertex(name = 'V_44',
+              particles = [ P.c__tilde__, P.s, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_43})
+
+V_45 = Vertex(name = 'V_45',
+              particles = [ P.t__tilde__, P.s, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_46})
+
+V_46 = Vertex(name = 'V_46',
+              particles = [ P.u__tilde__, P.b, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_41})
+
+V_47 = Vertex(name = 'V_47',
+              particles = [ P.c__tilde__, P.b, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_44})
+
+V_48 = Vertex(name = 'V_48',
+              particles = [ P.t__tilde__, P.b, P.W__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_47})
+
+V_52 = Vertex(name = 'V_52',
+              particles = [ P.e__plus__, P.e__minus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.FFS1 ],
+              couplings = {(0,0):C.GC_35})
+
+V_53 = Vertex(name = 'V_53',
+              particles = [ P.m__plus__, P.m__minus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.FFS1 ],
+              couplings = {(0,0):C.GC_36})
+
+V_54 = Vertex(name = 'V_54',
+              particles = [ P.tt__plus__, P.tt__minus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.FFS1 ],
+              couplings = {(0,0):C.GC_38})
+
+V_55 = Vertex(name = 'V_55',
+              particles = [ P.c__tilde__, P.c, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1 ],
+              couplings = {(0,0):C.GC_34})
+
+V_56 = Vertex(name = 'V_56',
+              particles = [ P.t__tilde__, P.t, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1 ],
+              couplings = {(0,0):C.GC_37})
+
+V_57 = Vertex(name = 'V_57',
+              particles = [ P.e__plus__, P.e__minus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2, L.FFV4 ],
+              couplings = {(0,0):C.GC_21,(0,1):C.GC_24})
+
+V_58 = Vertex(name = 'V_58',
+              particles = [ P.m__plus__, P.m__minus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2, L.FFV4 ],
+              couplings = {(0,0):C.GC_21,(0,1):C.GC_24})
+
+V_59 = Vertex(name = 'V_59',
+              particles = [ P.tt__plus__, P.tt__minus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2, L.FFV4 ],
+              couplings = {(0,0):C.GC_21,(0,1):C.GC_24})
+
+V_60 = Vertex(name = 'V_60',
+              particles = [ P.e__plus__, P.ve, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_61 = Vertex(name = 'V_61',
+              particles = [ P.m__plus__, P.vm, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_62 = Vertex(name = 'V_62',
+              particles = [ P.tt__plus__, P.vt, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_63 = Vertex(name = 'V_63',
+              particles = [ P.ve__tilde__, P.e__minus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_64 = Vertex(name = 'V_64',
+              particles = [ P.vm__tilde__, P.m__minus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_65 = Vertex(name = 'V_65',
+              particles = [ P.vt__tilde__, P.tt__minus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_11})
+
+V_66 = Vertex(name = 'V_66',
+              particles = [ P.u__tilde__, P.u, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV5 ],
+              couplings = {(0,0):C.GC_22,(0,1):C.GC_23})
+
+V_67 = Vertex(name = 'V_67',
+              particles = [ P.c__tilde__, P.c, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV5 ],
+              couplings = {(0,0):C.GC_22,(0,1):C.GC_23})
+
+V_68 = Vertex(name = 'V_68',
+              particles = [ P.t__tilde__, P.t, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV5 ],
+              couplings = {(0,0):C.GC_22,(0,1):C.GC_23})
+
+V_69 = Vertex(name = 'V_69',
+              particles = [ P.ve__tilde__, P.ve, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_28})
+
+V_70 = Vertex(name = 'V_70',
+              particles = [ P.vm__tilde__, P.vm, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_28})
+
+V_71 = Vertex(name = 'V_71',
+              particles = [ P.vt__tilde__, P.vt, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_28})
+
+# ======================
+# Goldstone interactions
+# ======================
+
+V_1001 = Vertex(name = 'V_1001',
+             particles = [ P.G0, P.G0, P.G0, P.G0 ],
+             color = [ '1' ],
+             lorentz = [ L.SSS1 ],
+             couplings = {(0,0):C.GC_1033})
+
+V_1002 = Vertex(name = 'V_1002',
+             particles = [ P.G0, P.G0, P.G__minus__, P.G__plus__ ],
+             color = [ '1' ],
+             lorentz = [ L.SSSS1 ],
+             couplings = {(0,0):C.GC_1031})
+
+V_1003 = Vertex(name = 'V_1003',
+             particles = [ P.G__minus__, P.G__minus__, P.G__plus__, P.G__plus__ ],
+             color = [ '1' ],
+             lorentz = [ L.SSSS1 ],
+             couplings = {(0,0):C.GC_1032})
+
+V_1004 = Vertex(name = 'V_1004',
+             particles = [ P.G0, P.G0, P.H, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSSS1 ],
+             couplings = {(0,0):C.GC_1031})
+
+V_1005 = Vertex(name = 'V_1005',
+             particles = [ P.G__minus__, P.G__plus__, P.H, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSSS1 ],
+             couplings = {(0,0):C.GC_1031})
+
+V_1007 = Vertex(name = 'V_1007',
+             particles = [ P.G0, P.G0, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSS1 ],
+             couplings = {(0,0):C.GC_1068})
+
+V_1008 = Vertex(name = 'V_1008',
+             particles = [ P.G__minus__, P.G__plus__, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSS1 ],
+             couplings = {(0,0):C.GC_1068})
+
+V_1010 = Vertex(name = 'V_1010',
+              particles = [ P.A, P.A, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1006})
+
+V_1011 = Vertex(name = 'V_1011',
+              particles = [ P.A, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_1003})
+
+V_1038 = Vertex(name = 'V_1038',
+              particles = [ P.A, P.W__minus__, P.G0, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1055})
+
+V_1039 = Vertex(name = 'V_1039',
+              particles = [ P.A, P.W__minus__, P.G__plus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1054})
+
+V_1040 = Vertex(name = 'V_1040',
+              particles = [ P.A, P.W__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_1074})
+
+V_1041 = Vertex(name = 'V_1041',
+              particles = [ P.W__minus__, P.G0, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_1039})
+
+V_1042 = Vertex(name = 'V_1042',
+              particles = [ P.W__minus__, P.G__plus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_1037})
+
+V_1044 = Vertex(name = 'V_1044',
+              particles = [ P.A, P.W__plus__, P.G0, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1055})
+
+V_1045 = Vertex(name = 'V_1045',
+              particles = [ P.A, P.W__plus__, P.G__minus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1056})
+
+V_1046 = Vertex(name = 'V_1046',
+              particles = [ P.A, P.W__plus__, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_1075})
+
+V_1047 = Vertex(name = 'V_1047',
+              particles = [ P.W__plus__, P.G0, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_1038})
+
+V_1048 = Vertex(name = 'V_1048',
+              particles = [ P.W__plus__, P.G__minus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_1037})
+
+V_1049 = Vertex(name = 'V_1049',
+              particles = [ P.W__minus__, P.W__plus__, P.G0, P.G0 ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1034})
+
+V_1050 = Vertex(name = 'V_1050',
+              particles = [ P.W__minus__, P.W__plus__, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1034})
+
+V_1056 = Vertex(name = 'V_1056',
+              particles = [ P.A, P.Z, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1063})
+
+V_1057 = Vertex(name = 'V_1057',
+              particles = [ P.Z, P.G0, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_1060})
+
+V_1058 = Vertex(name = 'V_1058',
+              particles = [ P.Z, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_1061})
+
+V_1059 = Vertex(name = 'V_1059',
+              particles = [ P.W__minus__, P.Z, P.G0, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1008})
+
+V_1060 = Vertex(name = 'V_1060',
+              particles = [ P.W__minus__, P.Z, P.G__plus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1009})
+
+V_1061 = Vertex(name = 'V_1061',
+              particles = [ P.W__minus__, P.Z, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_1067})
+
+V_1062 = Vertex(name = 'V_1062',
+              particles = [ P.W__plus__, P.Z, P.G0, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1008})
+
+V_1063 = Vertex(name = 'V_1063',
+              particles = [ P.W__plus__, P.Z, P.G__minus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1007})
+
+V_1064 = Vertex(name = 'V_1064',
+              particles = [ P.W__plus__, P.Z, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_1066})
+
+V_1066 = Vertex(name = 'V_1066',
+              particles = [ P.Z, P.Z, P.G0, P.G0 ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1065})
+
+V_1067 = Vertex(name = 'V_1067',
+              particles = [ P.Z, P.Z, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_1064})
+
+V_1077 = Vertex(name = 'V_1077',
+              particles = [ P.b__tilde__, P.b, P.G0 ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS2 ],
+              couplings = {(0,0):C.GC_1082})
+
+V_1082 = Vertex(name = 'V_1082',
+              particles = [ P.c__tilde__, P.d, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS3 ],
+              couplings = {(0,0):C.GC_1016})
+
+V_1083 = Vertex(name = 'V_1083',
+              particles = [ P.t__tilde__, P.d, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS3 ],
+              couplings = {(0,0):C.GC_1017})
+
+V_1084 = Vertex(name = 'V_1084',
+              particles = [ P.c__tilde__, P.s, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS3 ],
+              couplings = {(0,0):C.GC_1018})
+
+V_1085 = Vertex(name = 'V_1085',
+              particles = [ P.t__tilde__, P.s, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS3 ],
+              couplings = {(0,0):C.GC_1019})
+
+V_1086 = Vertex(name = 'V_1086',
+              particles = [ P.u__tilde__, P.b, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS8 ],
+              couplings = {(0,0):C.GC_1013})
+
+V_1087 = Vertex(name = 'V_1087',
+              particles = [ P.c__tilde__, P.b, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS8, L.FFS3 ],
+              couplings = {(0,0):C.GC_1014,(0,1):C.GC_1020})
+
+V_1088 = Vertex(name = 'V_1088',
+              particles = [ P.t__tilde__, P.b, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS8, L.FFS3 ],
+              couplings = {(0,0):C.GC_1015,(0,1):C.GC_1021})
+
+V_1101 = Vertex(name = 'V_1101',
+               particles = [ P.e__plus__, P.e__minus__, P.G0 ],
+               color = [ '1' ],
+               lorentz = [ L.FFS2 ],
+               couplings = {(0,0):C.GC_1088})
+
+V_1102 = Vertex(name = 'V_1102',
+               particles = [ P.m__plus__, P.m__minus__, P.G0 ],
+               color = [ '1' ],
+               lorentz = [ L.FFS2 ],
+               couplings = {(0,0):C.GC_1092})
+
+V_1103 = Vertex(name = 'V_1103',
+               particles = [ P.tt__plus__, P.tt__minus__, P.G0 ],
+               color = [ '1' ],
+               lorentz = [ L.FFS2 ],
+               couplings = {(0,0):C.GC_1098})
+
+V_1110 = Vertex(name = 'V_1110',
+               particles = [ P.ve__tilde__, P.e__minus__, P.G__plus__ ],
+               color = [ '1' ],
+               lorentz = [ L.FFS8 ],
+               couplings = {(0,0):C.GC_1087})
+
+V_1111 = Vertex(name = 'V_1111',
+               particles = [ P.vm__tilde__, P.m__minus__, P.G__plus__ ],
+               color = [ '1' ],
+               lorentz = [ L.FFS8 ],
+               couplings = {(0,0):C.GC_1091})
+
+V_1112 = Vertex(name = 'V_1112',
+               particles = [ P.vt__tilde__, P.tt__minus__, P.G__plus__ ],
+               color = [ '1' ],
+               lorentz = [ L.FFS8 ],
+               couplings = {(0,0):C.GC_1097})
+
+V_1116 = Vertex(name = 'V_1116',
+               particles = [ P.b__tilde__, P.u, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS3 ],
+               couplings = {(0,0):C.GC_1028})
+
+V_1117 = Vertex(name = 'V_1117',
+               particles = [ P.d__tilde__, P.c, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               couplings = {(0,0):C.GC_1022})
+
+V_1118 = Vertex(name = 'V_1118',
+               particles = [ P.s__tilde__, P.c, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               couplings = {(0,0):C.GC_1023})
+
+V_1119 = Vertex(name = 'V_1119',
+               particles = [ P.b__tilde__, P.c, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8, L.FFS3 ],
+               couplings = {(0,0):C.GC_1024,(0,1):C.GC_1029})
+
+V_1120 = Vertex(name = 'V_1120',
+               particles = [ P.d__tilde__, P.t, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               couplings = {(0,0):C.GC_1025})
+
+V_1121 = Vertex(name = 'V_1121',
+               particles = [ P.s__tilde__, P.t, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8 ],
+               couplings = {(0,0):C.GC_1026})
+
+V_1122 = Vertex(name = 'V_1122',
+               particles = [ P.b__tilde__, P.t, P.G__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS8, L.FFS3 ],
+               couplings = {(0,0):C.GC_1027,(0,1):C.GC_1030})
+
+V_1138 = Vertex(name = 'V_1138',
+               particles = [ P.c__tilde__, P.c, P.G0 ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS2 ],
+               couplings = {(0,0):C.GC_1085})
+
+V_1139 = Vertex(name = 'V_1139',
+               particles = [ P.t__tilde__, P.t, P.G0 ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFS2 ],
+               couplings = {(0,0):C.GC_1095})
+
+V_1145 = Vertex(name = 'V_1145',
+               particles = [ P.e__plus__, P.ve, P.G__minus__ ],
+               color = [ '1' ],
+               lorentz = [ L.FFS3 ],
+               couplings = {(0,0):C.GC_1086})
+
+V_1146 = Vertex(name = 'V_1146',
+               particles = [ P.m__plus__, P.vm, P.G__minus__ ],
+               color = [ '1' ],
+               lorentz = [ L.FFS3 ],
+               couplings = {(0,0):C.GC_1090})
+
+V_1147 = Vertex(name = 'V_1147',
+               particles = [ P.tt__plus__, P.vt, P.G__minus__ ],
+               color = [ '1' ],
+               lorentz = [ L.FFS3 ],
+               couplings = {(0,0):C.GC_1096})
diff --git a/Higgs/loop_sm_c3d4/write_param_card.py b/Higgs/loop_sm_c3d4/write_param_card.py
new file mode 100755
index 0000000000000000000000000000000000000000..57a85b061440f522b264b837e691f5b33accc96b
--- /dev/null
+++ b/Higgs/loop_sm_c3d4/write_param_card.py
@@ -0,0 +1,182 @@
+
+__date__ = "3 june 2010"
+__author__ = 'olivier.mattelaer@uclouvain.be'
+
+class ParamCardWriter(object):
+    
+    header = \
+    """######################################################################\n""" + \
+    """## PARAM_CARD AUTOMATICALY GENERATED BY THE UFO  #####################\n""" + \
+    """######################################################################\n"""   
+    
+    def __init__(self, filename, list_of_parameters=None, generic=False):
+        """write a valid param_card.dat"""
+        
+        if not list_of_parameters:
+            from parameters import all_parameters
+            list_of_parameters = [param for param in all_parameters if \
+                                                       param.nature=='external']
+        
+        self.generic_output = generic
+        if generic:
+            self.define_not_dep_param(list_of_parameters)
+
+        
+        self.fsock = open(filename, 'w')
+        self.fsock.write(self.header)
+        
+        self.write_card(list_of_parameters)
+    
+    def define_not_dep_param(self, list_of_parameters):
+        """define self.dep_mass and self.dep_width in case that they are 
+        requested in the param_card.dat"""
+        from particles import all_particles
+        
+        self.dep_mass = [(part, part.mass) for part in all_particles \
+                            if part.pdg_code > 0 and \
+                                            part.mass not in list_of_parameters]
+        self.dep_width = [(part, part.width) for part in all_particles\
+                             if part.pdg_code > 0 and \
+                                part.width not in list_of_parameters]        
+    
+    @staticmethod
+    def order_param(obj1, obj2):
+        """ order parameter of a given block """
+        
+        maxlen = min([len(obj1.lhacode), len(obj2.lhacode)])
+    
+        for i in range(maxlen):
+            if obj1.lhacode[i] < obj2.lhacode[i]:
+                return -1
+            elif obj1.lhacode[i] == obj2.lhacode[i]:
+                return 0
+            else:
+                return 1
+        #identical up to the first finish
+        if len(obj1.lhacode) > len(obj2.lhacode):
+            return 1
+        elif  len(obj1.lhacode) == len(obj2.lhacode):
+            return 0
+        else:
+            return -1
+        
+    def write_card(self, all_ext_param):
+        """ """
+        
+        # list all lhablock
+        all_lhablock = set([param.lhablock for param in all_ext_param])
+        
+        # ordonate lhablock alphabeticaly
+        all_lhablock = list(all_lhablock)
+        all_lhablock.sort()
+        # put at the beginning SMINPUT + MASS + DECAY
+        for name in ['DECAY', 'MASS','SMINPUTS']:
+            if name in all_lhablock:
+                all_lhablock.remove(name)
+                all_lhablock.insert(0, name)
+        
+        for lhablock in all_lhablock:
+            self.write_block(lhablock)
+            need_writing = [ param for param in all_ext_param if \
+                                                     param.lhablock == lhablock]
+            from functools import cmp_to_key
+            need_writing.sort(key=cmp_to_key(self.order_param))
+            [self.write_param(param, lhablock) for param in need_writing]
+            
+            if self.generic_output:
+                if lhablock in ['MASS', 'DECAY']:
+                    self.write_dep_param_block(lhablock)
+
+        if self.generic_output:
+            self.write_qnumber()
+                               
+    def write_block(self, name):
+        """ write a comment for a block"""
+        
+        self.fsock.writelines(
+        """\n###################################""" + \
+        """\n## INFORMATION FOR %s""" % name.upper() +\
+        """\n###################################\n"""
+         )
+        if name!='DECAY':
+            self.fsock.write("""Block %s \n""" % name)
+
+    def write_param(self, param, lhablock):
+        
+        lhacode=' '.join(['%3s' % key for key in param.lhacode])
+        if lhablock != 'DECAY':
+            text = """  %s %e # %s \n""" % (lhacode, complex(param.value).real, param.name ) 
+        else:
+            text = '''DECAY %s %e \n''' % (lhacode, complex(param.value).real)
+        self.fsock.write(text) 
+                    
+
+
+    
+    def write_dep_param_block(self, lhablock):
+        import cmath
+        from parameters import all_parameters
+        for parameter in all_parameters:
+            try:
+                exec("%s = %s" % (parameter.name, parameter.value))
+            except Exception:
+                pass
+        text = "##  Not dependent paramater.\n"
+        text += "## Those values should be edited following analytical the \n"
+        text += "## analytical expression. Some generator could simply ignore \n"
+        text += "## those values and use the analytical expression\n"
+        
+        if lhablock == 'MASS':
+            data = self.dep_mass
+            prefix = " "
+        else:
+            data = self.dep_width
+            prefix = "DECAY "
+        for part, param in data:
+            if isinstance(param.value, str):
+                value = complex(eval(param.value)).real
+            else:
+                value = param.value
+            
+            text += """%s %s %f # %s : %s \n""" %(prefix, part.pdg_code, 
+                        value, part.name, param.value)
+        self.fsock.write(text)    
+    
+    sm_pdg = [1,2,3,4,5,6,11,12,13,13,14,15,16,21,22,23,24,25]
+    data="""Block QNUMBERS %(pdg)d  # %(name)s 
+        1 %(charge)d  # 3 times electric charge
+        2 %(spin)d  # number of spin states (2S+1)
+        3 %(color)d  # colour rep (1: singlet, 3: triplet, 8: octet)
+        4 %(antipart)d  # Particle/Antiparticle distinction (0=own anti)\n"""
+    
+    def write_qnumber(self):
+        """ write qnumber """
+        from particles import all_particles
+        
+        text="""#===========================================================\n"""
+        text += """# QUANTUM NUMBERS OF NEW STATE(S) (NON SM PDG CODE)\n"""
+        text += """#===========================================================\n\n"""
+        
+        for part in all_particles:
+            if part.pdg_code in self.sm_pdg or part.pdg_code < 0:
+                continue
+            text += self.data % {'pdg': part.pdg_code,
+                                 'name': part.name,
+                                 'charge': 3 * part.charge,
+                                 'spin': 2 * part.spin + 1,
+                                 'color': part.color,
+                                 'antipart': part.name != part.antiname and 1 or 0}
+        
+        self.fsock.write(text)
+        
+        
+            
+            
+            
+            
+        
+            
+if '__main__' == __name__:
+    ParamCardWriter('./param_card.dat', generic=True)
+    print('write ./param_card.dat')
+    
diff --git a/model_list.txt b/model_list.txt
index 236e8c5f2408ade57abd4f4dd9d1e5f8a3bc4ca4..92b491db5f7dbafac272a8da7c1087fa9f207826 100644
--- a/model_list.txt
+++ b/model_list.txt
@@ -1039,6 +1039,13 @@ Requestor: Bernhard Meirose
 Content: Higgs-induced lepton-jets, replacing the usrmodv4_lj* models
 JIRA: https://its.cern.ch/jira/browse/AGENE-1105
 
+Higgs/loop_sm_c3d4
+Requestor: Nicholas Kyriacou. Model from Marko Stamenkovic in collaboration with others (Andreas Papaefstathiou et. al.)
+Contents: Model to Generate HHH Production. Allows to independently vary individual Higgs Coupling Parameters utitlized for Kappa-Scans. 
+Paper: https://arxiv.org/abs/2312.13562
+Source: https://github.com/mstamenk/triple-h-mc-gen/tree/main
+JIRA: https://its.cern.ch/jira/browse/ATLMCPROD-11023
+
 Higgs/loop_sm_twoscalar
 Requestor: Holly Pacey
 Content: Two real scalar model, with two extra higgs-like particles added. Loop-level implementation for gg-Scalar production vertices.