diff --git a/CHresonance_neutral_scalar_UFO.txt b/CHresonance_neutral_scalar_UFO.txt
new file mode 100644
index 0000000000000000000000000000000000000000..5cce9ef11644dd51f02aa90554ea08b24aa7e960
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO.txt
@@ -0,0 +1,4 @@
+Requestor: Zhen Wang
+Content: Singlet-catalyzed electroweak phase transitions from Ashutosh V. Kotwal
+Paper: https://journals.aps.org/prd/abstract/10.1103/PhysRevD.94.035022
+JIRA: https://its.cern.ch/jira/browse/ATLMCPROD-6413 (old ZZ sample request)
diff --git a/CHresonance_neutral_scalar_UFO/CHresonance_neutral_scalar_UFO.log b/CHresonance_neutral_scalar_UFO/CHresonance_neutral_scalar_UFO.log
new file mode 100644
index 0000000000000000000000000000000000000000..b92dd567abc51e6c6c5adbba0e35a9c2c146e1ba
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/CHresonance_neutral_scalar_UFO.log
@@ -0,0 +1,89 @@
+# This file was automatically created by FeynRules 2.0.2
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (February 23, 2011)
+# Date: Wed 8 Oct 2014 20:38:51
+
+
+# 
+# This is the logfile for the model CHresonance_neutral_scalar
+
+# Authors: N. Christensen, C. Duhr, B. Fuks
+# Model version: 1.4.5
+# Checking the Quantum numbers
+   * Electric charge defined.
+# Checking the Lagrangians
+   * All Lagrangians are ok.
+#
+# Particle definitions
+#
+
+   * No particles removed. All particles correspond to GenInt setup.
+
+# Automatically assigned PDG numbers
+   * Assigned PDG number 9000001 to particle ghA
+   * Assigned PDG number 9000002 to particle ghZ
+   * Assigned PDG number 9000003 to particle ghWp
+   * Assigned PDG number 9000004 to particle ghWm
+   * Assigned PDG number 9000005 to particle ghG
+   * Assigned PDG number 9000006 to particle s0
+
+
+# Compulsory PDG codes:
+   * Class SM leptons complete.
+   * Class SM neutrinos complete.
+   * Class SM quarks complete.
+   * Class SM gauge bosons complete.
+#
+# Parameter definitions
+#
+
+   * All parameters are ok.
+
+
+# Vertices
+   * Calling FeynmanRules for 2 Lagrangians.
+   * Number of classes vertices: 96
+   * Number of flavored vertices: 122
+   * Saved vertices in InterfaceRun[ 1 ].
+   * Checked QNumber conservation.
+      - Quantum number GhostNumber conserved in all vertices.
+      - Quantum number LeptonNumber conserved in all vertices.
+      - Quantum number Q conserved in all vertices.
+      - Quantum number Y conserved in all vertices.
+   * particles.py written.
+   * parameters.py written.
+#
+# Vertex definitions
+#
+
+   * 122 vertices written.
+   * vertices.py written.
+#
+# Lorentz structure definitions
+#
+
+   * 23 lorentz structures written.
+   * lorentz.py written.
+#
+# Coupling definitions
+#
+
+   * 69 couplings written.
+   * couplings.py written.
+#
+# Coupling order definitions
+#
+
+   * 0 couplings orders written.
+   * coupling_orders.py written.
+#
+# Decay definitions
+#
+
+   * 7 decays written.
+   * decay.py not written
+#
+# CTCoupling definitions
+#
+
+   * 0 CTcouplings written.
+   * CT_couplings.py written.
diff --git a/CHresonance_neutral_scalar_UFO/CT_couplings.py b/CHresonance_neutral_scalar_UFO/CT_couplings.py
new file mode 100644
index 0000000000000000000000000000000000000000..c8efdce6b30728dbbf31f922ab0dc0d4171e750e
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/CT_couplings.py
@@ -0,0 +1,11 @@
+# This file was automatically created by FeynRules 2.0.2
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (February 23, 2011)
+# Date: Wed 8 Oct 2014 20:39:09
+
+
+from object_library import all_couplings, Coupling
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec, cot
+
+
+
diff --git a/CHresonance_neutral_scalar_UFO/__init__.py b/CHresonance_neutral_scalar_UFO/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..76c53d16a666eb2c78dd2016063aafaf167c094a
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/__init__.py
@@ -0,0 +1,50 @@
+import function_library 
+import object_library 
+
+import particles
+import couplings
+import lorentz
+import parameters
+import vertices
+import coupling_orders
+import write_param_card
+import propagators
+
+
+all_particles = particles.all_particles
+all_vertices = vertices.all_vertices
+all_couplings = couplings.all_couplings
+all_lorentz = lorentz.all_lorentz
+all_parameters = parameters.all_parameters
+all_orders = coupling_orders.all_orders
+all_functions = function_library.all_functions
+all_propagators = propagators.all_propagators
+
+try:
+   import decays
+except ImportError:
+   pass
+else:
+   all_decays = decays.all_decays
+
+try:
+   import form_factors
+except ImportError:
+   pass
+else:
+   all_form_factors = form_factors.all_form_factors
+
+try:
+   import CT_vertices
+except ImportError:
+   pass
+else:
+   all_CTvertices = CT_vertices.all_CTvertices
+
+
+gauge = [0, 1]
+
+
+__author__ = "N. Christensen, C. Duhr, B. Fuks"
+__date__ = "21. 11. 2012"
+__version__= "1.4.5"
diff --git a/CHresonance_neutral_scalar_UFO/coupling_orders.py b/CHresonance_neutral_scalar_UFO/coupling_orders.py
new file mode 100644
index 0000000000000000000000000000000000000000..0178e7ece44ad44f2939dda6ea2e9a2e05f79b24
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/coupling_orders.py
@@ -0,0 +1,20 @@
+# This file was automatically created by FeynRules 2.0.2
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (February 23, 2011)
+# Date: Wed 8 Oct 2014 20:39:09
+
+
+from object_library import all_orders, CouplingOrder
+
+
+QCD = CouplingOrder(name = 'QCD',
+                    expansion_order = 99,
+                    hierarchy = 1)
+
+QED = CouplingOrder(name = 'QED',
+                    expansion_order = 99,
+                    hierarchy = 2)
+
+NP = CouplingOrder(name = 'NP',
+                   expansion_order = 99,
+                   hierarchy = 1)
+
diff --git a/CHresonance_neutral_scalar_UFO/couplings.py b/CHresonance_neutral_scalar_UFO/couplings.py
new file mode 100644
index 0000000000000000000000000000000000000000..0f053c4357a2c952647600c2b6274ad7390a8125
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/couplings.py
@@ -0,0 +1,287 @@
+# This file was automatically created by FeynRules 2.0.2
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (February 23, 2011)
+# Date: Wed 8 Oct 2014 20:39:09
+
+
+from object_library import all_couplings, Coupling
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec, cot
+
+
+
+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 = 'ee*complex(0,1)',
+                order = {'QED':1})
+
+GC_5 = Coupling(name = 'GC_5',
+                value = 'ee**2*complex(0,1)',
+                order = {'QED':2})
+
+GC_6 = Coupling(name = 'GC_6',
+                value = '2*ee**2*complex(0,1)',
+                order = {'QED':2})
+
+GC_7 = Coupling(name = 'GC_7',
+                value = '-ee**2/(2.*cw)',
+                order = {'QED':2})
+
+GC_8 = Coupling(name = 'GC_8',
+                value = '(ee**2*complex(0,1))/(2.*cw)',
+                order = {'QED':2})
+
+GC_9 = Coupling(name = 'GC_9',
+                value = 'ee**2/(2.*cw)',
+                order = {'QED':2})
+
+GC_10 = Coupling(name = 'GC_10',
+                 value = '-G',
+                 order = {'QCD':1})
+
+GC_11 = Coupling(name = 'GC_11',
+                 value = 'complex(0,1)*G',
+                 order = {'QCD':1})
+
+GC_12 = Coupling(name = 'GC_12',
+                 value = 'complex(0,1)*G**2',
+                 order = {'QCD':2})
+
+GC_13 = Coupling(name = 'GC_13',
+                 value = 'I1b33',
+                 order = {'QED':1})
+
+GC_14 = Coupling(name = 'GC_14',
+                 value = '-I2b33',
+                 order = {'QED':1})
+
+GC_15 = Coupling(name = 'GC_15',
+                 value = 'I3b33',
+                 order = {'QED':1})
+
+GC_16 = Coupling(name = 'GC_16',
+                 value = '-I4b33',
+                 order = {'QED':1})
+
+GC_17 = Coupling(name = 'GC_17',
+                 value = '-2*complex(0,1)*lam',
+                 order = {'QED':2})
+
+GC_18 = Coupling(name = 'GC_18',
+                 value = '-4*complex(0,1)*lam',
+                 order = {'QED':2})
+
+GC_19 = Coupling(name = 'GC_19',
+                 value = '-6*complex(0,1)*lam',
+                 order = {'QED':2})
+
+GC_20 = Coupling(name = 'GC_20',
+                 value = 'complex(0,1)*gweak*keta*MW',
+                 order = {'NP':1})
+
+GC_21 = Coupling(name = 'GC_21',
+                 value = '(complex(0,1)*gweak*keta*MZ)/cw',
+                 order = {'NP':1})
+
+GC_22 = Coupling(name = 'GC_22',
+                 value = '(ee**2*complex(0,1))/(2.*sw**2)',
+                 order = {'QED':2})
+
+GC_23 = Coupling(name = 'GC_23',
+                 value = '-((ee**2*complex(0,1))/sw**2)',
+                 order = {'QED':2})
+
+GC_24 = Coupling(name = 'GC_24',
+                 value = '(cw**2*ee**2*complex(0,1))/sw**2',
+                 order = {'QED':2})
+
+GC_25 = Coupling(name = 'GC_25',
+                 value = '-ee/(2.*sw)',
+                 order = {'QED':1})
+
+GC_26 = Coupling(name = 'GC_26',
+                 value = '-(ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_27 = Coupling(name = 'GC_27',
+                 value = '(ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_28 = Coupling(name = 'GC_28',
+                 value = '(ee*complex(0,1))/(sw*cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_29 = Coupling(name = 'GC_29',
+                 value = '-(cw*ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_30 = Coupling(name = 'GC_30',
+                 value = '(cw*ee*complex(0,1))/(2.*sw)',
+                 order = {'QED':1})
+
+GC_31 = Coupling(name = 'GC_31',
+                 value = '-((cw*ee*complex(0,1))/sw)',
+                 order = {'QED':1})
+
+GC_32 = Coupling(name = 'GC_32',
+                 value = '(cw*ee*complex(0,1))/sw',
+                 order = {'QED':1})
+
+GC_33 = Coupling(name = 'GC_33',
+                 value = '-ee**2/(2.*sw)',
+                 order = {'QED':2})
+
+GC_34 = Coupling(name = 'GC_34',
+                 value = '-(ee**2*complex(0,1))/(2.*sw)',
+                 order = {'QED':2})
+
+GC_35 = Coupling(name = 'GC_35',
+                 value = 'ee**2/(2.*sw)',
+                 order = {'QED':2})
+
+GC_36 = Coupling(name = 'GC_36',
+                 value = '(-2*cw*ee**2*complex(0,1))/sw',
+                 order = {'QED':2})
+
+GC_37 = Coupling(name = 'GC_37',
+                 value = '-(ee*complex(0,1)*sw)/(6.*cw)',
+                 order = {'QED':1})
+
+GC_38 = Coupling(name = 'GC_38',
+                 value = '(ee*complex(0,1)*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_39 = Coupling(name = 'GC_39',
+                 value = '-(cw*ee)/(2.*sw) - (ee*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_40 = Coupling(name = 'GC_40',
+                 value = '-(cw*ee*complex(0,1))/(2.*sw) + (ee*complex(0,1)*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_41 = Coupling(name = 'GC_41',
+                 value = '(cw*ee*complex(0,1))/(2.*sw) + (ee*complex(0,1)*sw)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_42 = Coupling(name = 'GC_42',
+                 value = '(cw*ee**2*complex(0,1))/sw - (ee**2*complex(0,1)*sw)/cw',
+                 order = {'QED':2})
+
+GC_43 = Coupling(name = 'GC_43',
+                 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_44 = Coupling(name = 'GC_44',
+                 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_45 = Coupling(name = 'GC_45',
+                 value = '-(ee**2*vev)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_46 = Coupling(name = 'GC_46',
+                 value = '(ee**2*vev)/(2.*cw)',
+                 order = {'QED':1})
+
+GC_47 = Coupling(name = 'GC_47',
+                 value = '-2*complex(0,1)*lam*vev',
+                 order = {'QED':1})
+
+GC_48 = Coupling(name = 'GC_48',
+                 value = '-6*complex(0,1)*lam*vev',
+                 order = {'QED':1})
+
+GC_49 = Coupling(name = 'GC_49',
+                 value = '-(ee**2*vev)/(4.*sw**2)',
+                 order = {'QED':1})
+
+GC_50 = Coupling(name = 'GC_50',
+                 value = '-(ee**2*complex(0,1)*vev)/(4.*sw**2)',
+                 order = {'QED':1})
+
+GC_51 = Coupling(name = 'GC_51',
+                 value = '(ee**2*complex(0,1)*vev)/(2.*sw**2)',
+                 order = {'QED':1})
+
+GC_52 = Coupling(name = 'GC_52',
+                 value = '(ee**2*vev)/(4.*sw**2)',
+                 order = {'QED':1})
+
+GC_53 = Coupling(name = 'GC_53',
+                 value = '-(ee**2*vev)/(2.*sw)',
+                 order = {'QED':1})
+
+GC_54 = Coupling(name = 'GC_54',
+                 value = '(ee**2*vev)/(2.*sw)',
+                 order = {'QED':1})
+
+GC_55 = Coupling(name = 'GC_55',
+                 value = '-(ee**2*vev)/(4.*cw) - (cw*ee**2*vev)/(4.*sw**2)',
+                 order = {'QED':1})
+
+GC_56 = Coupling(name = 'GC_56',
+                 value = '(ee**2*vev)/(4.*cw) - (cw*ee**2*vev)/(4.*sw**2)',
+                 order = {'QED':1})
+
+GC_57 = Coupling(name = 'GC_57',
+                 value = '-(ee**2*vev)/(4.*cw) + (cw*ee**2*vev)/(4.*sw**2)',
+                 order = {'QED':1})
+
+GC_58 = Coupling(name = 'GC_58',
+                 value = '(ee**2*vev)/(4.*cw) + (cw*ee**2*vev)/(4.*sw**2)',
+                 order = {'QED':1})
+
+GC_59 = Coupling(name = 'GC_59',
+                 value = '-(ee**2*complex(0,1)*vev)/2. - (cw**2*ee**2*complex(0,1)*vev)/(4.*sw**2) - (ee**2*complex(0,1)*sw**2*vev)/(4.*cw**2)',
+                 order = {'QED':1})
+
+GC_60 = Coupling(name = 'GC_60',
+                 value = 'ee**2*complex(0,1)*vev + (cw**2*ee**2*complex(0,1)*vev)/(2.*sw**2) + (ee**2*complex(0,1)*sw**2*vev)/(2.*cw**2)',
+                 order = {'QED':1})
+
+GC_61 = Coupling(name = 'GC_61',
+                 value = '(-2*complex(0,1)*keta)/vweak',
+                 order = {'NP':1})
+
+GC_62 = Coupling(name = 'GC_62',
+                 value = '-(yb/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_63 = Coupling(name = 'GC_63',
+                 value = '-((complex(0,1)*yb)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_64 = Coupling(name = 'GC_64',
+                 value = '-((complex(0,1)*yt)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_65 = Coupling(name = 'GC_65',
+                 value = 'yt/cmath.sqrt(2)',
+                 order = {'QED':1})
+
+GC_66 = Coupling(name = 'GC_66',
+                 value = '-ytau',
+                 order = {'QED':1})
+
+GC_67 = Coupling(name = 'GC_67',
+                 value = 'ytau',
+                 order = {'QED':1})
+
+GC_68 = Coupling(name = 'GC_68',
+                 value = '-(ytau/cmath.sqrt(2))',
+                 order = {'QED':1})
+
+GC_69 = Coupling(name = 'GC_69',
+                 value = '-((complex(0,1)*ytau)/cmath.sqrt(2))',
+                 order = {'QED':1})
+
diff --git a/CHresonance_neutral_scalar_UFO/decays.py b/CHresonance_neutral_scalar_UFO/decays.py
new file mode 100644
index 0000000000000000000000000000000000000000..8ebc6c5304ff6c62c07eaa7303e728e8033d577f
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/decays.py
@@ -0,0 +1,60 @@
+# This file was automatically created by FeynRules 2.0.2
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (February 23, 2011)
+# Date: Wed 8 Oct 2014 20:39:09
+
+
+from object_library import all_decays, Decay
+import particles as P
+
+
+Decay_H = Decay(name = 'Decay_H',
+                particle = P.H,
+                partial_widths = {(P.b,P.b__tilde__):'((-12*MB**2*yb**2 + 3*MH**2*yb**2)*cmath.sqrt(-4*MB**2*MH**2 + MH**4))/(16.*cmath.pi*abs(MH)**3)',
+                                  (P.ta__minus__,P.ta__plus__):'((MH**2*ytau**2 - 4*MTA**2*ytau**2)*cmath.sqrt(MH**4 - 4*MH**2*MTA**2))/(16.*cmath.pi*abs(MH)**3)',
+                                  (P.t,P.t__tilde__):'((3*MH**2*yt**2 - 12*MT**2*yt**2)*cmath.sqrt(MH**4 - 4*MH**2*MT**2))/(16.*cmath.pi*abs(MH)**3)',
+                                  (P.W__minus__,P.W__plus__):'(((3*ee**4*vev**2)/(4.*sw**4) + (ee**4*MH**4*vev**2)/(16.*MW**4*sw**4) - (ee**4*MH**2*vev**2)/(4.*MW**2*sw**4))*cmath.sqrt(MH**4 - 4*MH**2*MW**2))/(16.*cmath.pi*abs(MH)**3)',
+                                  (P.Z,P.Z):'(((9*ee**4*vev**2)/2. + (3*ee**4*MH**4*vev**2)/(8.*MZ**4) - (3*ee**4*MH**2*vev**2)/(2.*MZ**2) + (3*cw**4*ee**4*vev**2)/(4.*sw**4) + (cw**4*ee**4*MH**4*vev**2)/(16.*MZ**4*sw**4) - (cw**4*ee**4*MH**2*vev**2)/(4.*MZ**2*sw**4) + (3*cw**2*ee**4*vev**2)/sw**2 + (cw**2*ee**4*MH**4*vev**2)/(4.*MZ**4*sw**2) - (cw**2*ee**4*MH**2*vev**2)/(MZ**2*sw**2) + (3*ee**4*sw**2*vev**2)/cw**2 + (ee**4*MH**4*sw**2*vev**2)/(4.*cw**2*MZ**4) - (ee**4*MH**2*sw**2*vev**2)/(cw**2*MZ**2) + (3*ee**4*sw**4*vev**2)/(4.*cw**4) + (ee**4*MH**4*sw**4*vev**2)/(16.*cw**4*MZ**4) - (ee**4*MH**2*sw**4*vev**2)/(4.*cw**4*MZ**2))*cmath.sqrt(MH**4 - 4*MH**2*MZ**2))/(32.*cmath.pi*abs(MH)**3)'})
+
+Decay_Z = Decay(name = 'Decay_Z',
+                particle = P.Z,
+                partial_widths = {(P.W__minus__,P.W__plus__):'(((-12*cw**2*ee**2*MW**2)/sw**2 - (17*cw**2*ee**2*MZ**2)/sw**2 + (4*cw**2*ee**2*MZ**4)/(MW**2*sw**2) + (cw**2*ee**2*MZ**6)/(4.*MW**4*sw**2))*cmath.sqrt(-4*MW**2*MZ**2 + MZ**4))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.ve,P.ve__tilde__):'(MZ**2*(ee**2*MZ**2 + (cw**2*ee**2*MZ**2)/(2.*sw**2) + (ee**2*MZ**2*sw**2)/(2.*cw**2)))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.vm,P.vm__tilde__):'(MZ**2*(ee**2*MZ**2 + (cw**2*ee**2*MZ**2)/(2.*sw**2) + (ee**2*MZ**2*sw**2)/(2.*cw**2)))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.vt,P.vt__tilde__):'(MZ**2*(ee**2*MZ**2 + (cw**2*ee**2*MZ**2)/(2.*sw**2) + (ee**2*MZ**2*sw**2)/(2.*cw**2)))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.e__minus__,P.e__plus__):'(MZ**2*(-(ee**2*MZ**2) + (cw**2*ee**2*MZ**2)/(2.*sw**2) + (5*ee**2*MZ**2*sw**2)/(2.*cw**2)))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.mu__minus__,P.mu__plus__):'(MZ**2*(-(ee**2*MZ**2) + (cw**2*ee**2*MZ**2)/(2.*sw**2) + (5*ee**2*MZ**2*sw**2)/(2.*cw**2)))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.ta__minus__,P.ta__plus__):'((-5*ee**2*MTA**2 - ee**2*MZ**2 - (cw**2*ee**2*MTA**2)/(2.*sw**2) + (cw**2*ee**2*MZ**2)/(2.*sw**2) + (7*ee**2*MTA**2*sw**2)/(2.*cw**2) + (5*ee**2*MZ**2*sw**2)/(2.*cw**2))*cmath.sqrt(-4*MTA**2*MZ**2 + MZ**4))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.u,P.u__tilde__):'(MZ**2*(-(ee**2*MZ**2) + (3*cw**2*ee**2*MZ**2)/(2.*sw**2) + (17*ee**2*MZ**2*sw**2)/(6.*cw**2)))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.c,P.c__tilde__):'(MZ**2*(-(ee**2*MZ**2) + (3*cw**2*ee**2*MZ**2)/(2.*sw**2) + (17*ee**2*MZ**2*sw**2)/(6.*cw**2)))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.t,P.t__tilde__):'((-11*ee**2*MT**2 - ee**2*MZ**2 - (3*cw**2*ee**2*MT**2)/(2.*sw**2) + (3*cw**2*ee**2*MZ**2)/(2.*sw**2) + (7*ee**2*MT**2*sw**2)/(6.*cw**2) + (17*ee**2*MZ**2*sw**2)/(6.*cw**2))*cmath.sqrt(-4*MT**2*MZ**2 + MZ**4))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.d,P.d__tilde__):'(MZ**2*(ee**2*MZ**2 + (3*cw**2*ee**2*MZ**2)/(2.*sw**2) + (5*ee**2*MZ**2*sw**2)/(6.*cw**2)))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.s,P.s__tilde__):'(MZ**2*(ee**2*MZ**2 + (3*cw**2*ee**2*MZ**2)/(2.*sw**2) + (5*ee**2*MZ**2*sw**2)/(6.*cw**2)))/(48.*cmath.pi*abs(MZ)**3)',
+                                  (P.b,P.b__tilde__):'((-7*ee**2*MB**2 + ee**2*MZ**2 - (3*cw**2*ee**2*MB**2)/(2.*sw**2) + (3*cw**2*ee**2*MZ**2)/(2.*sw**2) - (17*ee**2*MB**2*sw**2)/(6.*cw**2) + (5*ee**2*MZ**2*sw**2)/(6.*cw**2))*cmath.sqrt(-4*MB**2*MZ**2 + MZ**4))/(48.*cmath.pi*abs(MZ)**3)'})
+
+Decay_W__plus__ = Decay(name = 'Decay_W__plus__',
+                        particle = P.W__plus__,
+                        partial_widths = {(P.ve,P.e__plus__):'(ee**2*MW**4)/(48.*cmath.pi*sw**2*abs(MW)**3)',
+                                          (P.vm,P.mu__plus__):'(ee**2*MW**4)/(48.*cmath.pi*sw**2*abs(MW)**3)',
+                                          (P.vt,P.ta__plus__):'((-MTA**2 + MW**2)*(-(ee**2*MTA**2)/(2.*sw**2) - (ee**2*MTA**4)/(2.*MW**2*sw**2) + (ee**2*MW**2)/sw**2))/(48.*cmath.pi*abs(MW)**3)',
+                                          (P.u,P.d__tilde__):'(ee**2*MW**4)/(16.*cmath.pi*sw**2*abs(MW)**3)',
+                                          (P.c,P.s__tilde__):'(ee**2*MW**4)/(16.*cmath.pi*sw**2*abs(MW)**3)',
+                                          (P.t,P.b__tilde__):'(((-3*ee**2*MB**2)/(2.*sw**2) - (3*ee**2*MT**2)/(2.*sw**2) - (3*ee**2*MB**4)/(2.*MW**2*sw**2) + (3*ee**2*MB**2*MT**2)/(MW**2*sw**2) - (3*ee**2*MT**4)/(2.*MW**2*sw**2) + (3*ee**2*MW**2)/sw**2)*cmath.sqrt(MB**4 - 2*MB**2*MT**2 + MT**4 - 2*MB**2*MW**2 - 2*MT**2*MW**2 + MW**4))/(48.*cmath.pi*abs(MW)**3)'})
+
+Decay_ta__minus__ = Decay(name = 'Decay_ta__minus__',
+                          particle = P.ta__minus__,
+                          partial_widths = {(P.W__minus__,P.vt):'((MTA**2 - MW**2)*((ee**2*MTA**2)/(2.*sw**2) + (ee**2*MTA**4)/(2.*MW**2*sw**2) - (ee**2*MW**2)/sw**2))/(32.*cmath.pi*abs(MTA)**3)'})
+
+Decay_b = Decay(name = 'Decay_b',
+                particle = P.b,
+                partial_widths = {(P.W__minus__,P.t):'(((3*ee**2*MB**2)/(2.*sw**2) + (3*ee**2*MT**2)/(2.*sw**2) + (3*ee**2*MB**4)/(2.*MW**2*sw**2) - (3*ee**2*MB**2*MT**2)/(MW**2*sw**2) + (3*ee**2*MT**4)/(2.*MW**2*sw**2) - (3*ee**2*MW**2)/sw**2)*cmath.sqrt(MB**4 - 2*MB**2*MT**2 + MT**4 - 2*MB**2*MW**2 - 2*MT**2*MW**2 + MW**4))/(96.*cmath.pi*abs(MB)**3)'})
+
+Decay_t = Decay(name = 'Decay_t',
+                particle = P.t,
+                partial_widths = {(P.W__plus__,P.b):'(((3*ee**2*MB**2)/(2.*sw**2) + (3*ee**2*MT**2)/(2.*sw**2) + (3*ee**2*MB**4)/(2.*MW**2*sw**2) - (3*ee**2*MB**2*MT**2)/(MW**2*sw**2) + (3*ee**2*MT**4)/(2.*MW**2*sw**2) - (3*ee**2*MW**2)/sw**2)*cmath.sqrt(MB**4 - 2*MB**2*MT**2 + MT**4 - 2*MB**2*MW**2 - 2*MT**2*MW**2 + MW**4))/(96.*cmath.pi*abs(MT)**3)'})
+
+Decay_s0 = Decay(name = 'Decay_s0',
+                 particle = P.s0,
+                 partial_widths = {(P.H,P.H):'(((4*keta**2*MH**4)/vweak**2 - (4*keta**2*MH**2*Ms0**2)/vweak**2 + (keta**2*Ms0**4)/vweak**2)*cmath.sqrt(-4*MH**2*Ms0**2 + Ms0**4))/(32.*cmath.pi*abs(Ms0)**3)',
+                                   (P.W__minus__,P.W__plus__):'((-(gweak**2*keta**2*Ms0**2) + (gweak**2*keta**2*Ms0**4)/(4.*MW**2) + 3*gweak**2*keta**2*MW**2)*cmath.sqrt(Ms0**4 - 4*Ms0**2*MW**2))/(16.*cmath.pi*abs(Ms0)**3)',
+                                   (P.Z,P.Z):'((-((gweak**2*keta**2*Ms0**2)/cw**2) + (gweak**2*keta**2*Ms0**4)/(4.*cw**2*MZ**2) + (3*gweak**2*keta**2*MZ**2)/cw**2)*cmath.sqrt(Ms0**4 - 4*Ms0**2*MZ**2))/(32.*cmath.pi*abs(Ms0)**3)'})
+
diff --git a/CHresonance_neutral_scalar_UFO/feynrules/Cabibbo.rst b/CHresonance_neutral_scalar_UFO/feynrules/Cabibbo.rst
new file mode 100644
index 0000000000000000000000000000000000000000..540206470c394bb07f17437e64fae6fd092e0cc9
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/feynrules/Cabibbo.rst
@@ -0,0 +1,16 @@
+(******************************************************************)
+(*     Restriction file for SM.fr                                 *)
+(*                                                                *)
+(*     All Yukawa couplings for light fermions (e,m,u,d,s)        *)
+(*     are put to zero.                                           *)
+(*                                                                *)
+(*     Only Cabibbo mixing is taken into account in CKM           *)
+(******************************************************************)
+
+M$Restrictions = {yl[i_] :> 0 /; (i === 1) || (i === 2),
+            yu[1] -> 0,
+            yd[i_] :> 0 /; (i === 1) || (i === 2),
+            CKM[3,3] -> 1,
+            CKM[i_,3] :> 0 /; (i === 1) || (i ===2),
+            CKM[3,i_] :> 0 /; (i === 1) || (i ===2)
+}
\ No newline at end of file
diff --git a/CHresonance_neutral_scalar_UFO/feynrules/DiagonalCKM.rst b/CHresonance_neutral_scalar_UFO/feynrules/DiagonalCKM.rst
new file mode 100644
index 0000000000000000000000000000000000000000..00e8816aa8f3a3ed2431a09170920a9d2d629d64
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/feynrules/DiagonalCKM.rst
@@ -0,0 +1,11 @@
+(******************************************************************)
+(*     Restriction file for SM.fr                                 *)
+(*                                                                *)                                            
+(*     The CKM matrix is diagonal                                 *)
+(******************************************************************)
+
+M$Restrictions = {
+            CKM[i_,i_] -> 1,
+            CKM[i_?NumericQ, j_?NumericQ] :> 0 /; (i =!= j),
+            cabi -> 0
+}
diff --git a/CHresonance_neutral_scalar_UFO/feynrules/Massless.rst b/CHresonance_neutral_scalar_UFO/feynrules/Massless.rst
new file mode 100644
index 0000000000000000000000000000000000000000..e8eeb32fa5513cb1fb4ba6ee8e5463db095d84de
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/feynrules/Massless.rst
@@ -0,0 +1,26 @@
+(******************************************************************)
+(*     Restriction file for SM.fr                                                     *)
+(*                                                                                                *)                                            
+(*     Only the t and b quarks, and the tau lepton are massive    *)
+(******************************************************************)
+
+M$Restrictions = {
+          Me -> 0,
+	  yme->0,
+	  ye->0,
+          MMU -> 0,
+	  ymm->0,
+	  ym->0,
+          MU -> 0,
+	  ymup->0,
+	  yup->0,
+          MD -> 0,
+	  ymdo->0,
+	  ydo->0,
+          MS -> 0,
+	  yms->0,
+	  ys->0,
+          MC -> 0,
+          ymc-> 0,
+          yc -> 0
+}
diff --git a/CHresonance_neutral_scalar_UFO/feynrules/Notebook.nb b/CHresonance_neutral_scalar_UFO/feynrules/Notebook.nb
new file mode 100644
index 0000000000000000000000000000000000000000..9e9ddd68a042edd0ffd025a223e617911bc828d4
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/feynrules/Notebook.nb
@@ -0,0 +1,1137 @@
+(* Content-type: application/vnd.wolfram.mathematica *)
+
+(*** Wolfram Notebook File ***)
+(* http://www.wolfram.com/nb *)
+
+(* CreatedBy='Mathematica 8.0' *)
+
+(*CacheID: 234*)
+(* Internal cache information:
+NotebookFileLineBreakTest
+NotebookFileLineBreakTest
+NotebookDataPosition[       157,          7]
+NotebookDataLength[     48082,       1128]
+NotebookOptionsPosition[     43974,        996]
+NotebookOutlinePosition[     44329,       1012]
+CellTagsIndexPosition[     44286,       1009]
+WindowFrame->Normal*)
+
+(* Beginning of Notebook Content *)
+Notebook[{
+
+Cell[CellGroupData[{
+Cell["\<\
+Feynrules - Composite Higgs Resonances\
+\>", "Title",
+ CellChangeTimes->{{3.605618752716743*^9, 3.605618757950574*^9}, {
+  3.6058757722119417`*^9, 3.605875777268808*^9}},
+ FontSize->24,
+ FontColor->GrayLevel[1],
+ Background->GrayLevel[0]],
+
+Cell[BoxData[
+ RowBox[{"Quit", "[", "]"}]], "Input",
+ CellChangeTimes->{{3.605618861320875*^9, 3.605618880207635*^9}, {
+  3.605635005168662*^9, 3.6056350058380632`*^9}}],
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"FR$Parallel", " ", "=", " ", "False"}], ";"}]], "Input",
+ CellChangeTimes->{{3.605619001829976*^9, 3.6056190561681433`*^9}, {
+  3.605619339401984*^9, 3.605619342273584*^9}, {3.60563473959864*^9, 
+  3.60563474030869*^9}}],
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"$FeynRulesPath", " ", "=", " ", 
+   RowBox[{
+   "SetDirectory", "[", "\"\</Users/mattlow/Software/FeynRules2/\>\"", 
+    "]"}]}], ";"}]], "Input"],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{"<<", " ", "FeynRules`"}]], "Input",
+ CellChangeTimes->{{3.605619060724291*^9, 3.6056190658268213`*^9}}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData["\<\" - FeynRules - \"\>"], "Print",
+ CellChangeTimes->{3.605875836985937*^9, 3.605875894272051*^9, 
+  3.605882343522704*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"\<\"Version: \"\>", "\[InvisibleSpace]", "\<\"2.0.2\"\>", 
+   "\[InvisibleSpace]", 
+   RowBox[{"\<\" (\"\>", " ", "\<\"21 October 2013\"\>"}], 
+   "\[InvisibleSpace]", "\<\").\"\>"}],
+  SequenceForm["Version: ", "2.0.2", " (" "21 October 2013", ")."],
+  Editable->False]], "Print",
+ CellChangeTimes->{3.605875836985937*^9, 3.605875894272051*^9, 
+  3.605882343527142*^9}],
+
+Cell[BoxData["\<\"Authors: A. Alloul, N. Christensen, C. Degrande, C. Duhr, \
+B. Fuks\"\>"], "Print",
+ CellChangeTimes->{3.605875836985937*^9, 3.605875894272051*^9, 
+  3.6058823435312967`*^9}],
+
+Cell[BoxData["\<\" \"\>"], "Print",
+ CellChangeTimes->{3.605875836985937*^9, 3.605875894272051*^9, 
+  3.605882343535253*^9}],
+
+Cell[BoxData["\<\"Please cite:\"\>"], "Print",
+ CellChangeTimes->{3.605875836985937*^9, 3.605875894272051*^9, 
+  3.605882343539241*^9}],
+
+Cell[BoxData["\<\"    - arXiv:1310.1921;\"\>"], "Print",
+ CellChangeTimes->{3.605875836985937*^9, 3.605875894272051*^9, 
+  3.6058823435431843`*^9}],
+
+Cell[BoxData["\<\"    - Comput.Phys.Commun.180:1614-1641,2009 \
+(arXiv:0806.4194).\"\>"], "Print",
+ CellChangeTimes->{3.605875836985937*^9, 3.605875894272051*^9, 
+  3.605882343547357*^9}],
+
+Cell[BoxData["\<\" \"\>"], "Print",
+ CellChangeTimes->{3.605875836985937*^9, 3.605875894272051*^9, 
+  3.6058823435516768`*^9}],
+
+Cell[BoxData["\<\"http://feynrules.phys.ucl.ac.be\"\>"], "Print",
+ CellChangeTimes->{3.605875836985937*^9, 3.605875894272051*^9, 
+  3.605882343556036*^9}],
+
+Cell[BoxData["\<\" \"\>"], "Print",
+ CellChangeTimes->{3.605875836985937*^9, 3.605875894272051*^9, 
+  3.6058823435600643`*^9}],
+
+Cell[BoxData["\<\"The FeynRules palette can be opened using the command \
+FRPalette[].\"\>"], "Print",
+ CellChangeTimes->{3.605875836985937*^9, 3.605875894272051*^9, 
+  3.60588234356409*^9}]
+}, Open  ]],
+
+Cell[BoxData[
+ RowBox[{
+  StyleBox[
+   RowBox[{"Get", "::", "noopen"}], "MessageName"], 
+  RowBox[{
+  ":", " "}], "\<\"Cannot open \[NoBreak]\\!\\(\\\"FAToFR.m\\\"\\)\[NoBreak]. \
+\\!\\(\\*ButtonBox[\\\"\[RightSkeleton]\\\", ButtonStyle->\\\"Link\\\", \
+ButtonFrame->None, ButtonData:>\\\"paclet:ref/message/General/noopen\\\", \
+ButtonNote -> \\\"Get::noopen\\\"]\\)\"\>"}]], "Message", "MSG",
+ CellChangeTimes->{{3.605619585572277*^9, 3.6056196077975883`*^9}, 
+   3.60561985205474*^9, 3.605622418597084*^9, 3.605622470939652*^9, 
+   3.6056225643872128`*^9, 3.605622678178093*^9, 3.605622772301466*^9, 
+   3.605622914978488*^9, 3.605632947062608*^9, 3.60563319057358*^9, 
+   3.605633274960883*^9, 3.605633321355266*^9, 3.605633438275936*^9, 
+   3.6056347456852093`*^9, {3.6056350124966793`*^9, 3.605635039259606*^9}, 
+   3.605637108779278*^9, {3.605637195679698*^9, 3.605637225028183*^9}, 
+   3.605637287992268*^9, 3.605637342560418*^9, 3.605637481033449*^9, 
+   3.605637539495367*^9, 3.605637610404883*^9, 3.6056376582399263`*^9, 
+   3.605637702718152*^9, 3.60565295699905*^9, 3.6056531223195457`*^9, 
+   3.6056532842527037`*^9, 3.6056539820680227`*^9, 3.605654206239273*^9, {
+   3.605654290239663*^9, 3.6056542944258842`*^9}, 3.605656362606439*^9, 
+   3.6056582925990067`*^9, 3.605658508408112*^9, {3.605658683656948*^9, 
+   3.605658709791214*^9}, 3.605658784466889*^9, 3.6056595237947407`*^9, 
+   3.6056598816239147`*^9, 3.605661357905642*^9, 3.605665013304759*^9, 
+   3.605665247749601*^9, 3.605665442368574*^9, {3.605875816573057*^9, 
+   3.605875839906644*^9}, 3.605875895188623*^9, 3.605882344549234*^9}]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell["Load Model ", "Section",
+ CellChangeTimes->{{3.6056192692699947`*^9, 3.605619290461377*^9}, {
+   3.605635085707984*^9, 3.605635087980074*^9}, 3.605875784340205*^9}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[{
+ RowBox[{"SetDirectory", "[", 
+  "\"\</Users/mattlow/Software/FeynRules2/Models/Model_CHresonance\>\"", 
+  "]"}], "\[IndentingNewLine]", 
+ RowBox[{"LoadModel", "[", 
+  RowBox[{"\"\<SM.fr\>\"", ",", "\"\<neutral_scalar.fr\>\""}], 
+  "]"}], "\[IndentingNewLine]", 
+ RowBox[{"LoadRestriction", "[", 
+  RowBox[{"\"\<DiagonalCKM.rst\>\"", ",", "\"\<Massless.rst\>\""}], 
+  "]"}]}], "Input",
+ CellChangeTimes->{{3.605619424691307*^9, 3.6056194362970877`*^9}, {
+   3.6056195317958803`*^9, 3.605619568447425*^9}, {3.6056224948230658`*^9, 
+   3.605622541801384*^9}, {3.60562265570363*^9, 3.605622665292665*^9}, 
+   3.6056332001464043`*^9, {3.6058758508413973`*^9, 3.605875859894517*^9}}],
+
+Cell[BoxData["\<\"/Users/mattlow/Software/FeynRules2/Models/Model_CHresonance\
+\"\>"], "Output",
+ CellChangeTimes->{{3.6056196162034683`*^9, 3.605619643073296*^9}, 
+   3.605619861571529*^9, 3.6056224285835543`*^9, {3.605622526112595*^9, 
+   3.605622568823061*^9}, 3.6056226814832*^9, 3.605622775100486*^9, 
+   3.60562291710483*^9, 3.60563294983956*^9, 3.6056331940543613`*^9, 
+   3.605633276844449*^9, 3.605633402448938*^9, 3.605633439866679*^9, 
+   3.605634859906128*^9, {3.605635015189953*^9, 3.605635042654422*^9}, 
+   3.605637112245749*^9, {3.605637198372223*^9, 3.605637227469658*^9}, 
+   3.605656364074188*^9, 3.605658686433711*^9, {3.6058758605426702`*^9, 
+   3.6058758982809267`*^9}, 3.60588234818725*^9}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData["\<\"Merging model-files...\"\>"], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.605882348191485*^9}],
+
+Cell[BoxData["\<\"This model implementation was created by\"\>"], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.605882348238291*^9}],
+
+Cell[BoxData["\<\"N. Christensen\"\>"], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.6058823482424088`*^9}],
+
+Cell[BoxData["\<\"C. Duhr\"\>"], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.6058823482463627`*^9}],
+
+Cell[BoxData["\<\"B. Fuks\"\>"], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.605882348250326*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"\<\"Model Version: \"\>", "\[InvisibleSpace]", "\<\"1.4.5\"\>"}],
+  SequenceForm["Model Version: ", "1.4.5"],
+  Editable->False]], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.6058823482544327`*^9}],
+
+Cell[BoxData["\<\"http://feynrules.phys.ucl.ac.be/view/Main/StandardModel\"\>\
+"], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.6058823482585087`*^9}],
+
+Cell[BoxData["\<\"For more information, type ModelInformation[].\"\>"], \
+"Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.6058823482624826`*^9}],
+
+Cell[BoxData["\<\"\"\>"], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.605882348266357*^9}],
+
+Cell[BoxData["\<\"   - Loading particle classes.\"\>"], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.605882348270246*^9}],
+
+Cell[BoxData["\<\"   - Loading gauge group classes.\"\>"], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.605882348475081*^9}],
+
+Cell[BoxData["\<\"   - Loading parameter classes.\"\>"], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.605882348508407*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"\<\"\\nModel \"\>", 
+   "\[InvisibleSpace]", "\<\"CHresonance_neutral_scalar\"\>", 
+   "\[InvisibleSpace]", "\<\" loaded.\"\>"}],
+  SequenceForm["\nModel ", "CHresonance_neutral_scalar", " loaded."],
+  Editable->False]], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.6058823486639442`*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"\<\"Loading restrictions from \"\>", 
+   "\[InvisibleSpace]", "\<\"DiagonalCKM.rst\"\>", 
+   "\[InvisibleSpace]", "\<\" : \"\>", "\[InvisibleSpace]", 
+   DynamicBox[ToBoxes[PRIVATE`FR$restrictionCounter, StandardForm],
+    ImageSizeCache->{14., {1., 8.}}], "\[InvisibleSpace]", "\<\" / \"\>", 
+   "\[InvisibleSpace]", "3"}],
+  SequenceForm["Loading restrictions from ", "DiagonalCKM.rst", " : ", 
+   Dynamic[PRIVATE`FR$restrictionCounter], " / ", 3],
+  Editable->False]], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.60588234873158*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"\<\"Loading restrictions from \"\>", 
+   "\[InvisibleSpace]", "\<\"Massless.rst\"\>", 
+   "\[InvisibleSpace]", "\<\" : \"\>", "\[InvisibleSpace]", 
+   DynamicBox[ToBoxes[PRIVATE`FR$restrictionCounter, StandardForm],
+    ImageSizeCache->{14., {1., 8.}}], "\[InvisibleSpace]", "\<\" / \"\>", 
+   "\[InvisibleSpace]", "18"}],
+  SequenceForm["Loading restrictions from ", "Massless.rst", " : ", 
+   Dynamic[PRIVATE`FR$restrictionCounter], " / ", 18],
+  Editable->False]], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.6058823488115664`*^9}],
+
+Cell[BoxData["\<\"Restrictions loaded.\"\>"], "Print",
+ CellChangeTimes->{
+  3.605619616741845*^9, 3.6056196504421997`*^9, 3.605619861575075*^9, 
+   3.6056224286336527`*^9, {3.605622526115447*^9, 3.605622568856142*^9}, 
+   3.6056226814865217`*^9, 3.605622775103384*^9, 3.605622917108451*^9, 
+   3.605632949842614*^9, 3.6056331940577917`*^9, 3.605633276847851*^9, 
+   3.605633402452071*^9, 3.605633439870932*^9, 3.605634859909452*^9, {
+   3.6056350151932497`*^9, 3.6056350426578293`*^9}, 3.605637112250169*^9, {
+   3.60563719837819*^9, 3.6056372274730253`*^9}, 3.605656364078393*^9, 
+   3.6056586864378366`*^9, {3.605875860546157*^9, 3.605875898284733*^9}, 
+   3.605882348815343*^9}]
+}, Open  ]]
+}, Open  ]]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell["Check Model", "Section",
+ CellChangeTimes->{{3.605623086167871*^9, 3.6056231049293633`*^9}, 
+   3.605635092783983*^9, 3.605875903677682*^9}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[{
+ RowBox[{"CheckMassSpectrum", "[", "Lscalar", "]"}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"CheckHermiticity", "[", "Lscalar", "]"}], ";"}]}], "Input",
+ CellChangeTimes->{{3.605619622855464*^9, 3.6056196232911158`*^9}, {
+   3.605619982287343*^9, 3.605619987388823*^9}, {3.605620054327195*^9, 
+   3.605620085036086*^9}, {3.605622638451754*^9, 3.605622646849626*^9}, {
+   3.60562268863967*^9, 3.605622718368464*^9}, {3.6056227788700457`*^9, 
+   3.60562278929569*^9}, {3.605622851705021*^9, 3.6056228520710707`*^9}, {
+   3.605622922349184*^9, 3.605622928774396*^9}, {3.605622958918388*^9, 
+   3.605622964678646*^9}, {3.605622998944463*^9, 3.60562300546445*^9}, {
+   3.6056329595409737`*^9, 3.60563296304221*^9}, {3.605656367479803*^9, 
+   3.6056563806305923`*^9}, 3.605658690688363*^9, {3.605875911163392*^9, 
+   3.6058759147629223`*^9}}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"\<\"Neglecting all terms with more than \"\>", 
+   "\[InvisibleSpace]", "\<\"2\"\>", 
+   "\[InvisibleSpace]", "\<\" particles.\"\>"}],
+  SequenceForm["Neglecting all terms with more than ", "2", " particles."],
+  Editable->False]], "Print",
+ CellChangeTimes->{
+  3.6056229291812897`*^9, 3.605622965269404*^9, {3.6056229995035*^9, 
+   3.605623006396628*^9}, {3.6056329541402884`*^9, 3.6056329635693913`*^9}, 
+   3.605633203461426*^9, 3.605633279819338*^9, 3.605633405717062*^9, 
+   3.605633445043635*^9, 3.60563486324809*^9, 3.605635045897534*^9, {
+   3.605637201873412*^9, 3.605637230937489*^9}, 3.605656382532639*^9, 
+   3.605658691840357*^9, 3.6058759152676086`*^9, 3.605882351286916*^9}],
+
+Cell[BoxData["\<\"All mass terms are diagonal.\"\>"], "Print",
+ CellChangeTimes->{
+  3.6056229291812897`*^9, 3.605622965269404*^9, {3.6056229995035*^9, 
+   3.605623006396628*^9}, {3.6056329541402884`*^9, 3.6056329635693913`*^9}, 
+   3.605633203461426*^9, 3.605633279819338*^9, 3.605633405717062*^9, 
+   3.605633445043635*^9, 3.60563486324809*^9, 3.605635045897534*^9, {
+   3.605637201873412*^9, 3.605637230937489*^9}, 3.605656382532639*^9, 
+   3.605658691840357*^9, 3.6058759152676086`*^9, 3.6058823512904463`*^9}],
+
+Cell[BoxData["\<\"Getting mass spectrum.\"\>"], "Print",
+ CellChangeTimes->{
+  3.6056229291812897`*^9, 3.605622965269404*^9, {3.6056229995035*^9, 
+   3.605623006396628*^9}, {3.6056329541402884`*^9, 3.6056329635693913`*^9}, 
+   3.605633203461426*^9, 3.605633279819338*^9, 3.605633405717062*^9, 
+   3.605633445043635*^9, 3.60563486324809*^9, 3.605635045897534*^9, {
+   3.605637201873412*^9, 3.605637230937489*^9}, 3.605656382532639*^9, 
+   3.605658691840357*^9, 3.6058759152676086`*^9, 3.605882351322145*^9}],
+
+Cell[BoxData["\<\"Checking for less then 0.1% agreement with model file \
+values.\"\>"], "Print",
+ CellChangeTimes->{
+  3.6056229291812897`*^9, 3.605622965269404*^9, {3.6056229995035*^9, 
+   3.605623006396628*^9}, {3.6056329541402884`*^9, 3.6056329635693913`*^9}, 
+   3.605633203461426*^9, 3.605633279819338*^9, 3.605633405717062*^9, 
+   3.605633445043635*^9, 3.60563486324809*^9, 3.605635045897534*^9, {
+   3.605637201873412*^9, 3.605637230937489*^9}, 3.605656382532639*^9, 
+   3.605658691840357*^9, 3.6058759152676086`*^9, 3.605882351325905*^9}]
+}, Open  ]],
+
+Cell[BoxData[
+ TagBox[GridBox[{
+    {"\<\"Particle\"\>", "\<\"Analytic value\"\>", "\<\"Numerical value\"\>", \
+"\<\"Model-file value\"\>"},
+    {"s0", 
+     SqrtBox[
+      SuperscriptBox["Ms0", "2"]], "1000.`", "1000.`"}
+   },
+   GridBoxAlignment->{
+    "Columns" -> {{Left}}, "ColumnsIndexed" -> {}, "Rows" -> {{Baseline}}, 
+     "RowsIndexed" -> {}},
+   GridBoxSpacings->{"Columns" -> {
+       Offset[0.27999999999999997`], {
+        Offset[2.0999999999999996`]}, 
+       Offset[0.27999999999999997`]}, "ColumnsIndexed" -> {}, "Rows" -> {
+       Offset[0.2], {
+        Offset[0.4]}, 
+       Offset[0.2]}, "RowsIndexed" -> {}}],
+  Function[BoxForm`e$, 
+   TableForm[BoxForm`e$]]]], "Output",
+ CellChangeTimes->{
+  3.605623006536682*^9, {3.605632954179488*^9, 3.605632963736952*^9}, 
+   3.605633203613308*^9, 3.605633279970613*^9, 3.605633405866434*^9, 
+   3.605633445195354*^9, 3.6056348633652163`*^9, 3.605635046055908*^9, {
+   3.605637202024962*^9, 3.6056372311093893`*^9}, 3.605656382681074*^9, 
+   3.605658691987759*^9, 3.60587591541451*^9, 3.605882351391329*^9}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData["\<\"Checking for hermiticity by calculating the Feynman rules \
+contained in L-HC[L].\"\>"], "Print",
+ CellChangeTimes->{
+  3.6056229291812897`*^9, 3.605622965269404*^9, {3.6056229995035*^9, 
+   3.605623006396628*^9}, {3.6056329541402884`*^9, 3.6056329635693913`*^9}, 
+   3.605633203461426*^9, 3.605633279819338*^9, 3.605633405717062*^9, 
+   3.605633445043635*^9, 3.60563486324809*^9, 3.605635045897534*^9, {
+   3.605637201873412*^9, 3.605637230937489*^9}, 3.605656382532639*^9, 
+   3.605658691840357*^9, 3.6058759152676086`*^9, 3.605882351394699*^9}],
+
+Cell[BoxData["\<\"If the lagrangian is hermitian, then the number of vertices \
+should be zero.\"\>"], "Print",
+ CellChangeTimes->{
+  3.6056229291812897`*^9, 3.605622965269404*^9, {3.6056229995035*^9, 
+   3.605623006396628*^9}, {3.6056329541402884`*^9, 3.6056329635693913`*^9}, 
+   3.605633203461426*^9, 3.605633279819338*^9, 3.605633405717062*^9, 
+   3.605633445043635*^9, 3.60563486324809*^9, 3.605635045897534*^9, {
+   3.605637201873412*^9, 3.605637230937489*^9}, 3.605656382532639*^9, 
+   3.605658691840357*^9, 3.6058759152676086`*^9, 3.605882351397386*^9}],
+
+Cell[BoxData[
+ StyleBox["\<\"Starting Feynman rule calculation.\"\>",
+  StripOnInput->False,
+  LineColor->RGBColor[1, 0.5, 0],
+  FrontFaceColor->RGBColor[1, 0.5, 0],
+  BackFaceColor->RGBColor[1, 0.5, 0],
+  GraphicsColor->RGBColor[1, 0.5, 0],
+  FontWeight->Bold,
+  FontColor->RGBColor[1, 0.5, 0]]], "Print",
+ CellChangeTimes->{
+  3.6056229291812897`*^9, 3.605622965269404*^9, {3.6056229995035*^9, 
+   3.605623006396628*^9}, {3.6056329541402884`*^9, 3.6056329635693913`*^9}, 
+   3.605633203461426*^9, 3.605633279819338*^9, 3.605633405717062*^9, 
+   3.605633445043635*^9, 3.60563486324809*^9, 3.605635045897534*^9, {
+   3.605637201873412*^9, 3.605637230937489*^9}, 3.605656382532639*^9, 
+   3.605658691840357*^9, 3.6058759152676086`*^9, 3.60588235148717*^9}],
+
+Cell[BoxData["\<\"Expanding the Lagrangian...\"\>"], "Print",
+ CellChangeTimes->{
+  3.6056229291812897`*^9, 3.605622965269404*^9, {3.6056229995035*^9, 
+   3.605623006396628*^9}, {3.6056329541402884`*^9, 3.6056329635693913`*^9}, 
+   3.605633203461426*^9, 3.605633279819338*^9, 3.605633405717062*^9, 
+   3.605633445043635*^9, 3.60563486324809*^9, 3.605635045897534*^9, {
+   3.605637201873412*^9, 3.605637230937489*^9}, 3.605656382532639*^9, 
+   3.605658691840357*^9, 3.6058759152676086`*^9, 3.605882351523753*^9}],
+
+Cell[BoxData["\<\"No vertices found.\"\>"], "Print",
+ CellChangeTimes->{
+  3.6056229291812897`*^9, 3.605622965269404*^9, {3.6056229995035*^9, 
+   3.605623006396628*^9}, {3.6056329541402884`*^9, 3.6056329635693913`*^9}, 
+   3.605633203461426*^9, 3.605633279819338*^9, 3.605633405717062*^9, 
+   3.605633445043635*^9, 3.60563486324809*^9, 3.605635045897534*^9, {
+   3.605637201873412*^9, 3.605637230937489*^9}, 3.605656382532639*^9, 
+   3.605658691840357*^9, 3.6058759152676086`*^9, 3.605882351528597*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"0", "\[InvisibleSpace]", "\<\" vertices obtained.\"\>"}],
+  SequenceForm[0, " vertices obtained."],
+  Editable->False]], "Print",
+ CellChangeTimes->{
+  3.6056229291812897`*^9, 3.605622965269404*^9, {3.6056229995035*^9, 
+   3.605623006396628*^9}, {3.6056329541402884`*^9, 3.6056329635693913`*^9}, 
+   3.605633203461426*^9, 3.605633279819338*^9, 3.605633405717062*^9, 
+   3.605633445043635*^9, 3.60563486324809*^9, 3.605635045897534*^9, {
+   3.605637201873412*^9, 3.605637230937489*^9}, 3.605656382532639*^9, 
+   3.605658691840357*^9, 3.6058759152676086`*^9, 3.605882351533533*^9}],
+
+Cell[BoxData["\<\"The lagrangian is hermitian.\"\>"], "Print",
+ CellChangeTimes->{
+  3.6056229291812897`*^9, 3.605622965269404*^9, {3.6056229995035*^9, 
+   3.605623006396628*^9}, {3.6056329541402884`*^9, 3.6056329635693913`*^9}, 
+   3.605633203461426*^9, 3.605633279819338*^9, 3.605633405717062*^9, 
+   3.605633445043635*^9, 3.60563486324809*^9, 3.605635045897534*^9, {
+   3.605637201873412*^9, 3.605637230937489*^9}, 3.605656382532639*^9, 
+   3.605658691840357*^9, 3.6058759152676086`*^9, 3.605882351537565*^9}]
+}, Open  ]]
+}, Open  ]]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell["\<\
+Feynman Rules (Fermion Partner)\
+\>", "Section",
+ CellChangeTimes->{{3.605623116445725*^9, 3.605623132577363*^9}, 
+   3.6056350952868013`*^9}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[{
+ RowBox[{"vertices", " ", "=", " ", 
+  RowBox[{"FeynmanRules", "[", "Lscalar", "]"}]}], "\[IndentingNewLine]", 
+ RowBox[{
+  RowBox[{"decays", " ", "=", " ", 
+   RowBox[{"ComputeWidths", "[", "vertices", "]"}]}], ";"}]}], "Input",
+ CellChangeTimes->{{3.605623135518034*^9, 3.6056231403631487`*^9}, {
+   3.6056329889705153`*^9, 3.6056329893374767`*^9}, {3.60563366954454*^9, 
+   3.6056336908373327`*^9}, {3.605633868933378*^9, 3.6056338702728558`*^9}, 
+   3.605633908964724*^9, {3.605633941603565*^9, 3.605633953097578*^9}, 
+   3.605634867413897*^9, {3.6056349686727057`*^9, 3.60563497715202*^9}, {
+   3.605637056408601*^9, 3.6056370590299*^9}, {3.605875923923554*^9, 
+   3.6058759561163816`*^9}}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ StyleBox["\<\"Starting Feynman rule calculation.\"\>",
+  StripOnInput->False,
+  LineColor->RGBColor[1, 0.5, 0],
+  FrontFaceColor->RGBColor[1, 0.5, 0],
+  BackFaceColor->RGBColor[1, 0.5, 0],
+  GraphicsColor->RGBColor[1, 0.5, 0],
+  FontWeight->Bold,
+  FontColor->RGBColor[1, 0.5, 0]]], "Print",
+ CellChangeTimes->{
+  3.605623141036067*^9, 3.605632989814084*^9, 3.6056332084743757`*^9, 
+   3.605633283183111*^9, 3.605633408314115*^9, 3.6056334475492563`*^9, 
+   3.605633691403208*^9, 3.6056338708193607`*^9, 3.605633909533334*^9, 
+   3.6056339535259933`*^9, 3.605634867962195*^9, 3.605635050541889*^9, 
+   3.605637238858077*^9, {3.605875937769837*^9, 3.605875956666177*^9}, 
+   3.6058823543890657`*^9}],
+
+Cell[BoxData["\<\"Expanding the Lagrangian...\"\>"], "Print",
+ CellChangeTimes->{
+  3.605623141036067*^9, 3.605632989814084*^9, 3.6056332084743757`*^9, 
+   3.605633283183111*^9, 3.605633408314115*^9, 3.6056334475492563`*^9, 
+   3.605633691403208*^9, 3.6056338708193607`*^9, 3.605633909533334*^9, 
+   3.6056339535259933`*^9, 3.605634867962195*^9, 3.605635050541889*^9, 
+   3.605637238858077*^9, {3.605875937769837*^9, 3.605875956666177*^9}, 
+   3.6058823543926353`*^9}],
+
+Cell[BoxData["\<\"Collecting the different structures that enter the \
+vertex.\"\>"], "Print",
+ CellChangeTimes->{
+  3.605623141036067*^9, 3.605632989814084*^9, 3.6056332084743757`*^9, 
+   3.605633283183111*^9, 3.605633408314115*^9, 3.6056334475492563`*^9, 
+   3.605633691403208*^9, 3.6056338708193607`*^9, 3.605633909533334*^9, 
+   3.6056339535259933`*^9, 3.605634867962195*^9, 3.605635050541889*^9, 
+   3.605637238858077*^9, {3.605875937769837*^9, 3.605875956666177*^9}, 
+   3.605882354442403*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{
+  "2", "\[InvisibleSpace]", "\<\" possible non-zero vertices have been found \
+-> starting the computation: \"\>", "\[InvisibleSpace]", 
+   DynamicBox[ToBoxes[FeynRules`FR$FeynmanRules, StandardForm],
+    ImageSizeCache->{7., {0., 8.}}], "\[InvisibleSpace]", "\<\" / \"\>", 
+   "\[InvisibleSpace]", "2", "\[InvisibleSpace]", "\<\".\"\>"}],
+  SequenceForm[
+  2, " possible non-zero vertices have been found -> starting the \
+computation: ", 
+   Dynamic[FeynRules`FR$FeynmanRules], " / ", 2, "."],
+  Editable->False]], "Print",
+ CellChangeTimes->{
+  3.605623141036067*^9, 3.605632989814084*^9, 3.6056332084743757`*^9, 
+   3.605633283183111*^9, 3.605633408314115*^9, 3.6056334475492563`*^9, 
+   3.605633691403208*^9, 3.6056338708193607`*^9, 3.605633909533334*^9, 
+   3.6056339535259933`*^9, 3.605634867962195*^9, 3.605635050541889*^9, 
+   3.605637238858077*^9, {3.605875937769837*^9, 3.605875956666177*^9}, 
+   3.605882354448645*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"2", "\[InvisibleSpace]", "\<\" vertices obtained.\"\>"}],
+  SequenceForm[2, " vertices obtained."],
+  Editable->False]], "Print",
+ CellChangeTimes->{
+  3.605623141036067*^9, 3.605632989814084*^9, 3.6056332084743757`*^9, 
+   3.605633283183111*^9, 3.605633408314115*^9, 3.6056334475492563`*^9, 
+   3.605633691403208*^9, 3.6056338708193607`*^9, 3.605633909533334*^9, 
+   3.6056339535259933`*^9, 3.605634867962195*^9, 3.605635050541889*^9, 
+   3.605637238858077*^9, {3.605875937769837*^9, 3.605875956666177*^9}, 
+   3.605882354452979*^9}]
+}, Open  ]],
+
+Cell[BoxData[
+ RowBox[{"{", 
+  RowBox[{
+   RowBox[{"{", 
+    RowBox[{
+     RowBox[{"{", 
+      RowBox[{
+       RowBox[{"{", 
+        RowBox[{"s0", ",", "1"}], "}"}], ",", 
+       RowBox[{"{", 
+        RowBox[{"W", ",", "2"}], "}"}], ",", 
+       RowBox[{"{", 
+        RowBox[{
+         SuperscriptBox["W", "\[Dagger]"], ",", "3"}], "}"}]}], "}"}], ",", 
+     FractionBox[
+      RowBox[{"\[ImaginaryI]", " ", "e", " ", "keta", " ", 
+       SubscriptBox["M", "W"], " ", 
+       SubscriptBox["\[Eta]", 
+        RowBox[{
+         SubscriptBox["\<\"\[Mu]\"\>", "2"], ",", 
+         SubscriptBox["\<\"\[Mu]\"\>", "3"]}]]}], 
+      SubscriptBox["s", "w"]]}], "}"}], ",", 
+   RowBox[{"{", 
+    RowBox[{
+     RowBox[{"{", 
+      RowBox[{
+       RowBox[{"{", 
+        RowBox[{"s0", ",", "1"}], "}"}], ",", 
+       RowBox[{"{", 
+        RowBox[{"Z", ",", "2"}], "}"}], ",", 
+       RowBox[{"{", 
+        RowBox[{"Z", ",", "3"}], "}"}]}], "}"}], ",", 
+     FractionBox[
+      RowBox[{"\[ImaginaryI]", " ", "e", " ", "keta", " ", "MZ", " ", 
+       SubscriptBox["\[Eta]", 
+        RowBox[{
+         SubscriptBox["\<\"\[Mu]\"\>", "2"], ",", 
+         SubscriptBox["\<\"\[Mu]\"\>", "3"]}]]}], 
+      RowBox[{
+       SubscriptBox["c", "w"], " ", 
+       SubscriptBox["s", "w"]}]]}], "}"}]}], "}"}]], "Output",
+ CellChangeTimes->{
+  3.605633909797886*^9, 3.6056339538109617`*^9, 3.605634868137611*^9, 
+   3.605635050886277*^9, 3.605637239169332*^9, {3.6058759379425573`*^9, 
+   3.605875956729602*^9}, 3.6058823544581203`*^9}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"\<\"Flavor expansion of the vertices: \"\>", "\[InvisibleSpace]", 
+   DynamicBox[ToBoxes[FeynRules`FR$Count1, StandardForm],
+    ImageSizeCache->{14., {1., 8.}}], "\[InvisibleSpace]", "\<\" / \"\>", 
+   "\[InvisibleSpace]", "2"}],
+  SequenceForm["Flavor expansion of the vertices: ", 
+   Dynamic[FeynRules`FR$Count1], " / ", 2],
+  Editable->False]], "Print",
+ CellChangeTimes->{
+  3.605623141036067*^9, 3.605632989814084*^9, 3.6056332084743757`*^9, 
+   3.605633283183111*^9, 3.605633408314115*^9, 3.6056334475492563`*^9, 
+   3.605633691403208*^9, 3.6056338708193607`*^9, 3.605633909533334*^9, 
+   3.6056339535259933`*^9, 3.605634867962195*^9, 3.605635050541889*^9, 
+   3.605637238858077*^9, {3.605875937769837*^9, 3.605875956666177*^9}, 
+   3.60588235451572*^9}],
+
+Cell[BoxData[
+ StyleBox["\<\"Computing the squared matrix elements relevant for the 1->2 \
+decays: \"\>",
+  StripOnInput->False,
+  LineColor->RGBColor[1, 0.5, 0],
+  FrontFaceColor->RGBColor[1, 0.5, 0],
+  BackFaceColor->RGBColor[1, 0.5, 0],
+  GraphicsColor->RGBColor[1, 0.5, 0],
+  FontWeight->Bold,
+  FontColor->RGBColor[1, 0.5, 0]]], "Print",
+ CellChangeTimes->{
+  3.605623141036067*^9, 3.605632989814084*^9, 3.6056332084743757`*^9, 
+   3.605633283183111*^9, 3.605633408314115*^9, 3.6056334475492563`*^9, 
+   3.605633691403208*^9, 3.6056338708193607`*^9, 3.605633909533334*^9, 
+   3.6056339535259933`*^9, 3.605634867962195*^9, 3.605635050541889*^9, 
+   3.605637238858077*^9, {3.605875937769837*^9, 3.605875956666177*^9}, 
+   3.605882354559538*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{
+   DynamicBox[ToBoxes[FeynRules`FR$DecayCounter, StandardForm],
+    ImageSizeCache->{14., {0., 8.}}], "\[InvisibleSpace]", "\<\" / \"\>", 
+   "\[InvisibleSpace]", "2"}],
+  SequenceForm[
+   Dynamic[FeynRules`FR$DecayCounter], " / ", 2],
+  Editable->False]], "Print",
+ CellChangeTimes->{
+  3.605623141036067*^9, 3.605632989814084*^9, 3.6056332084743757`*^9, 
+   3.605633283183111*^9, 3.605633408314115*^9, 3.6056334475492563`*^9, 
+   3.605633691403208*^9, 3.6056338708193607`*^9, 3.605633909533334*^9, 
+   3.6056339535259933`*^9, 3.605634867962195*^9, 3.605635050541889*^9, 
+   3.605637238858077*^9, {3.605875937769837*^9, 3.605875956666177*^9}, 
+   3.605882354596119*^9}]
+}, Open  ]]
+}, Open  ]],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ RowBox[{
+  RowBox[{"WriteUFO", "[", 
+   RowBox[{"LSM", ",", "Lscalar"}], "]"}], ";"}]], "Input",
+ CellChangeTimes->{{3.605623414495763*^9, 3.605623422610941*^9}, {
+  3.605635073368265*^9, 3.6056350736934977`*^9}, {3.605875960611538*^9, 
+  3.605875962361438*^9}}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData["\<\" --- Universal FeynRules Output (UFO) v 1.1 ---\"\>"], \
+"Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.6058823609681177`*^9}],
+
+Cell[BoxData[
+ StyleBox["\<\"Starting Feynman rule calculation.\"\>",
+  StripOnInput->False,
+  LineColor->RGBColor[1, 0.5, 0],
+  FrontFaceColor->RGBColor[1, 0.5, 0],
+  BackFaceColor->RGBColor[1, 0.5, 0],
+  GraphicsColor->RGBColor[1, 0.5, 0],
+  FontWeight->Bold,
+  FontColor->RGBColor[1, 0.5, 0]]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.60588236156849*^9}],
+
+Cell[BoxData["\<\"Expanding the Lagrangian...\"\>"], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882361577166*^9}],
+
+Cell[BoxData["\<\"Collecting the different structures that enter the \
+vertex.\"\>"], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.6058823669100437`*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{
+  "98", "\[InvisibleSpace]", "\<\" possible non-zero vertices have been found \
+-> starting the computation: \"\>", "\[InvisibleSpace]", 
+   DynamicBox[ToBoxes[FeynRules`FR$FeynmanRules, StandardForm],
+    ImageSizeCache->{7., {0., 8.}}], "\[InvisibleSpace]", "\<\" / \"\>", 
+   "\[InvisibleSpace]", "98", "\[InvisibleSpace]", "\<\".\"\>"}],
+  SequenceForm[
+  98, " possible non-zero vertices have been found -> starting the \
+computation: ", 
+   Dynamic[FeynRules`FR$FeynmanRules], " / ", 98, "."],
+  Editable->False]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.6058823670079393`*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"93", "\[InvisibleSpace]", "\<\" vertices obtained.\"\>"}],
+  SequenceForm[93, " vertices obtained."],
+  Editable->False]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882372273816*^9}],
+
+Cell[BoxData[
+ StyleBox["\<\"Starting Feynman rule calculation.\"\>",
+  StripOnInput->False,
+  LineColor->RGBColor[1, 0.5, 0],
+  FrontFaceColor->RGBColor[1, 0.5, 0],
+  BackFaceColor->RGBColor[1, 0.5, 0],
+  GraphicsColor->RGBColor[1, 0.5, 0],
+  FontWeight->Bold,
+  FontColor->RGBColor[1, 0.5, 0]]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882372275631*^9}],
+
+Cell[BoxData["\<\"Expanding the Lagrangian...\"\>"], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.6058823722780933`*^9}],
+
+Cell[BoxData["\<\"Collecting the different structures that enter the \
+vertex.\"\>"], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882372330551*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{
+  "2", "\[InvisibleSpace]", "\<\" possible non-zero vertices have been found \
+-> starting the computation: \"\>", "\[InvisibleSpace]", 
+   DynamicBox[ToBoxes[FeynRules`FR$FeynmanRules, StandardForm],
+    ImageSizeCache->{7., {0., 8.}}], "\[InvisibleSpace]", "\<\" / \"\>", 
+   "\[InvisibleSpace]", "2", "\[InvisibleSpace]", "\<\".\"\>"}],
+  SequenceForm[
+  2, " possible non-zero vertices have been found -> starting the \
+computation: ", 
+   Dynamic[FeynRules`FR$FeynmanRules], " / ", 2, "."],
+  Editable->False]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882372333057*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"2", "\[InvisibleSpace]", "\<\" vertices obtained.\"\>"}],
+  SequenceForm[2, " vertices obtained."],
+  Editable->False]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882372394369*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"\<\"Flavor expansion of the vertices: \"\>", "\[InvisibleSpace]", 
+   DynamicBox[ToBoxes[FeynRules`FR$Count1, StandardForm],
+    ImageSizeCache->{14., {1., 8.}}], "\[InvisibleSpace]", "\<\" / \"\>", 
+   "\[InvisibleSpace]", "95"}],
+  SequenceForm["Flavor expansion of the vertices: ", 
+   Dynamic[FeynRules`FR$Count1], " / ", 95],
+  Editable->False]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.6058823763363256`*^9}],
+
+Cell[BoxData["\<\"   - Saved vertices in InterfaceRun[ 1 ].\"\>"], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882379280229*^9}],
+
+Cell[BoxData[
+ StyleBox["\<\"Computing the squared matrix elements relevant for the 1->2 \
+decays: \"\>",
+  StripOnInput->False,
+  LineColor->RGBColor[1, 0.5, 0],
+  FrontFaceColor->RGBColor[1, 0.5, 0],
+  BackFaceColor->RGBColor[1, 0.5, 0],
+  GraphicsColor->RGBColor[1, 0.5, 0],
+  FontWeight->Bold,
+  FontColor->RGBColor[1, 0.5, 0]]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882379342313*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{
+   DynamicBox[ToBoxes[FeynRules`FR$DecayCounter, StandardForm],
+    ImageSizeCache->{14., {0., 8.}}], "\[InvisibleSpace]", "\<\" / \"\>", 
+   "\[InvisibleSpace]", "42"}],
+  SequenceForm[
+   Dynamic[FeynRules`FR$DecayCounter], " / ", 42],
+  Editable->False]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882379352541*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"\<\"Squared matrix elent compute in \"\>", "\[InvisibleSpace]", 
+   "1.580134000000001`", "\[InvisibleSpace]", "\<\" seconds.\"\>"}],
+  SequenceForm[
+  "Squared matrix elent compute in ", 1.580134000000001, " seconds."],
+  Editable->False]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.6058823814476423`*^9}]
+}, Open  ]],
+
+Cell[BoxData[
+ RowBox[{
+  StyleBox[
+   RowBox[{"Decay", "::", "Overwrite"}], "MessageName"], 
+  RowBox[{
+  ":", " "}], "\<\"Warning: Previous results will be overwritten.\"\>"}]], \
+"Message", "MSG",
+ CellChangeTimes->{3.605623441917754*^9, 3.605637266795199*^9, 
+  3.605875986706703*^9, 3.605882381659711*^9}],
+
+Cell[CellGroupData[{
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"\<\"Decay widths computed in \"\>", "\[InvisibleSpace]", 
+   "0.025485999999997233`", "\[InvisibleSpace]", "\<\" seconds.\"\>"}],
+  SequenceForm[
+  "Decay widths computed in ", 0.025485999999997233`, " seconds."],
+  Editable->False]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.60588238166367*^9}],
+
+Cell[BoxData["\<\"Preparing Python output.\"\>"], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.6058823816654043`*^9}],
+
+Cell[BoxData["\<\"    - Splitting vertices into building blocks.\"\>"], \
+"Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.6058823818620367`*^9}],
+
+Cell[BoxData[
+ InterpretationBox[
+  RowBox[{"\<\"    - Optimizing: \"\>", "\[InvisibleSpace]", 
+   DynamicBox[ToBoxes[PRIVATE`PY$SplitVertexCounter, StandardForm],
+    ImageSizeCache->{22., {0., 8.}}], "\[InvisibleSpace]", "\<\"/\"\>", 
+   "\[InvisibleSpace]", "121", "\[InvisibleSpace]", "\<\" .\"\>"}],
+  SequenceForm["    - Optimizing: ", 
+   Dynamic[PRIVATE`PY$SplitVertexCounter], "/", 121, " ."],
+  Editable->False]], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882382105023*^9}],
+
+Cell[BoxData["\<\"    - Writing files.\"\>"], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882382296308*^9}],
+
+Cell[BoxData["\<\"Done!\"\>"], "Print",
+ CellChangeTimes->{3.605623426254333*^9, 3.605637251108576*^9, 
+  3.605875969708352*^9, 3.605882382577735*^9}]
+}, Open  ]]
+}, Open  ]]
+}, Open  ]]
+}, Open  ]]
+},
+WindowSize->{740, 656},
+WindowMargins->{{Automatic, 8}, {Automatic, 0}},
+FrontEndVersion->"8.0 for Mac OS X x86 (32-bit, 64-bit Kernel) (February 23, \
+2011)",
+StyleDefinitions->"Default.nb"
+]
+(* End of Notebook Content *)
+
+(* Internal cache information *)
+(*CellTagsOutline
+CellTagsIndex->{}
+*)
+(*CellTagsIndex
+CellTagsIndex->{}
+*)
+(*NotebookFileOutline
+Notebook[{
+Cell[CellGroupData[{
+Cell[579, 22, 248, 7, 80, "Title"],
+Cell[830, 31, 169, 3, 27, "Input"],
+Cell[1002, 36, 255, 5, 27, "Input"],
+Cell[1260, 43, 181, 5, 27, "Input"],
+Cell[CellGroupData[{
+Cell[1466, 52, 127, 2, 27, "Input"],
+Cell[CellGroupData[{
+Cell[1618, 58, 138, 2, 20, "Print"],
+Cell[1759, 62, 415, 9, 20, "Print"],
+Cell[2177, 73, 192, 3, 20, "Print"],
+Cell[2372, 78, 124, 2, 20, "Print"],
+Cell[2499, 82, 135, 2, 20, "Print"],
+Cell[2637, 86, 147, 2, 20, "Print"],
+Cell[2787, 90, 187, 3, 20, "Print"],
+Cell[2977, 95, 126, 2, 20, "Print"],
+Cell[3106, 99, 154, 2, 20, "Print"],
+Cell[3263, 103, 126, 2, 20, "Print"],
+Cell[3392, 107, 190, 3, 20, "Print"]
+}, Open  ]],
+Cell[3597, 113, 1611, 25, 23, "Message"]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[5245, 143, 170, 2, 67, "Section"],
+Cell[CellGroupData[{
+Cell[5440, 149, 693, 13, 58, "Input"],
+Cell[6136, 164, 713, 10, 27, "Output"],
+Cell[CellGroupData[{
+Cell[6874, 178, 684, 10, 20, "Print"],
+Cell[7561, 190, 702, 10, 20, "Print"],
+Cell[8266, 202, 678, 10, 20, "Print"],
+Cell[8947, 214, 671, 10, 20, "Print"],
+Cell[9621, 226, 669, 10, 20, "Print"],
+Cell[10293, 238, 814, 14, 20, "Print"],
+Cell[11110, 254, 721, 11, 20, "Print"],
+Cell[11834, 267, 712, 11, 20, "Print"],
+Cell[12549, 280, 662, 10, 20, "Print"],
+Cell[13214, 292, 692, 10, 20, "Print"],
+Cell[13909, 304, 695, 10, 20, "Print"],
+Cell[14607, 316, 693, 10, 20, "Print"],
+Cell[15303, 328, 904, 16, 36, "Print"],
+Cell[16210, 346, 1152, 20, 20, "Print"],
+Cell[17365, 368, 1151, 20, 20, "Print"],
+Cell[18519, 390, 682, 10, 20, "Print"]
+}, Open  ]]
+}, Open  ]]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[19262, 407, 146, 2, 67, "Section"],
+Cell[CellGroupData[{
+Cell[19433, 413, 856, 13, 43, "Input"],
+Cell[CellGroupData[{
+Cell[20314, 430, 735, 13, 20, "Print"],
+Cell[21052, 445, 514, 7, 20, "Print"],
+Cell[21569, 454, 506, 7, 20, "Print"],
+Cell[22078, 463, 547, 8, 20, "Print"]
+}, Open  ]],
+Cell[22640, 474, 1069, 25, 67, "Output"],
+Cell[CellGroupData[{
+Cell[23734, 503, 565, 8, 20, "Print"],
+Cell[24302, 513, 561, 8, 20, "Print"],
+Cell[24866, 523, 755, 15, 20, "Print"],
+Cell[25624, 540, 511, 7, 20, "Print"],
+Cell[26138, 549, 502, 7, 20, "Print"],
+Cell[26643, 558, 624, 11, 20, "Print"],
+Cell[27270, 571, 512, 7, 20, "Print"]
+}, Open  ]]
+}, Open  ]]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[27843, 585, 152, 4, 67, "Section"],
+Cell[CellGroupData[{
+Cell[28020, 593, 710, 12, 43, "Input"],
+Cell[CellGroupData[{
+Cell[28755, 609, 713, 15, 20, "Print"],
+Cell[29471, 626, 468, 7, 20, "Print"],
+Cell[29942, 635, 499, 8, 20, "Print"],
+Cell[30444, 645, 975, 19, 20, "Print"],
+Cell[31422, 666, 579, 11, 20, "Print"]
+}, Open  ]],
+Cell[32016, 680, 1508, 44, 95, "Output"],
+Cell[CellGroupData[{
+Cell[33549, 728, 807, 15, 20, "Print"],
+Cell[34359, 745, 747, 16, 20, "Print"],
+Cell[35109, 763, 715, 15, 20, "Print"]
+}, Open  ]]
+}, Open  ]],
+Cell[CellGroupData[{
+Cell[35873, 784, 277, 6, 27, "Input"],
+Cell[CellGroupData[{
+Cell[36175, 794, 196, 3, 20, "Print"],
+Cell[36374, 799, 416, 10, 20, "Print"],
+Cell[36793, 811, 172, 2, 20, "Print"],
+Cell[36968, 815, 207, 3, 20, "Print"],
+Cell[37178, 820, 687, 14, 20, "Print"],
+Cell[37868, 836, 287, 6, 20, "Print"],
+Cell[38158, 844, 417, 10, 20, "Print"],
+Cell[38578, 856, 174, 2, 20, "Print"],
+Cell[38755, 860, 205, 3, 20, "Print"],
+Cell[38963, 865, 681, 14, 20, "Print"],
+Cell[39647, 881, 285, 6, 20, "Print"],
+Cell[39935, 889, 518, 10, 20, "Print"],
+Cell[40456, 901, 186, 2, 20, "Print"],
+Cell[40645, 905, 453, 11, 20, "Print"],
+Cell[41101, 918, 423, 10, 20, "Print"],
+Cell[41527, 930, 408, 8, 20, "Print"]
+}, Open  ]],
+Cell[41950, 941, 310, 8, 23, "Message"],
+Cell[CellGroupData[{
+Cell[42285, 953, 398, 8, 20, "Print"],
+Cell[42686, 963, 171, 2, 20, "Print"],
+Cell[42860, 967, 195, 3, 20, "Print"],
+Cell[43058, 972, 543, 10, 20, "Print"],
+Cell[43604, 984, 165, 2, 20, "Print"],
+Cell[43772, 988, 150, 2, 20, "Print"]
+}, Open  ]]
+}, Open  ]]
+}, Open  ]]
+}, Open  ]]
+}
+]
+*)
+
+(* End of internal cache information *)
diff --git a/CHresonance_neutral_scalar_UFO/feynrules/SM.fr b/CHresonance_neutral_scalar_UFO/feynrules/SM.fr
new file mode 100644
index 0000000000000000000000000000000000000000..504b0d71d77cc88c9b5b396239536609f3172477
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/feynrules/SM.fr
@@ -0,0 +1,729 @@
+(***************************************************************************************************************)
+(******                       This is the FeynRules mod-file for the Standard model                       ******)
+(******                                                                                                   ******)
+(******     Authors: N. Christensen, C. Duhr, B. Fuks                                                     ******)
+(******                                                                                                   ******)
+(****** Choose whether Feynman gauge is desired.                                                          ******)
+(****** If set to False, unitary gauge is assumed.                                                          ****)
+(****** Feynman gauge is especially useful for CalcHEP/CompHEP where the calculation is 10-100 times faster. ***)
+(****** Feynman gauge is not supported in MadGraph and Sherpa.                                              ****)
+(***************************************************************************************************************)
+
+(* ************************** *)
+(* *****  Information   ***** *)
+(* ************************** *)
+M$ModelName = "Standard Model";
+
+M$Information = {
+  Authors      -> {"N. Christensen", "C. Duhr", "B. Fuks"}, 
+  Version      -> "1.4.5",
+  Date         -> "21. 11. 2012",
+  Institutions -> {"Michigan State University", "Universite catholique de Louvain (CP3)", "IPHC Strasbourg / University of Strasbourg"},
+  Emails       -> {"neil@pa.msu.edu", "claude.duhr@uclouvain.be", "benjamin.fuks@cnrs.in2p3.fr"},
+  URLs         -> "http://feynrules.phys.ucl.ac.be/view/Main/StandardModel"
+};
+
+FeynmanGauge = True;
+
+(* ************************** *)
+(* *****  Change  log   ***** *)
+(* ************************** *)
+
+(* v1.4.5: Added widths for ghosts.                                          *)
+(* v1.4.4: Changed widths of goldstone bosons to be the same as for the W and Z bosons *)
+(* v1.4.3: Updated conventions for the symmetric structure constants of SU3. *)
+(* v1.4.2: Set FeynmanGauge=True as default again.                           *)
+(* v1.4: Added SU(2) representation.                                         *)
+(*       -> Modification in the field declarations (doublets are added)      *)
+(*       -> Modification in the Lagrangian (much simpler).                   *)
+(* v1.3: Added yukawa couplings for all fermions for gauge invariance.       *)
+(*       Added yukawa couplings for 1st generation fermions to Massless.rst. *)
+(*       Updated parameters to PDG 2010.                                     *)
+(* v1.2: Set FeynmanGauge=True as default.                                   *)
+(*       Set Gluonic ghosts to be included in both gauges.                   *)
+(* v1.1: Fixed yukawa couplings in Feynman gauge.                            *)
+(*       Changed yd[n] CKM[n,m] to yd[m] CKM[n,m].                           *)
+(*       Changed yu[n] Conjugate[CKM[m,n]] to yu[m] Conjugate[CKM[m,n]].     *)
+
+(* ************************** *)
+(* *****      vevs      ***** *)
+(* ************************** *)
+M$vevs = { {Phi[2],vev} };
+
+(* ************************** *)
+(* *****  Gauge groups  ***** *)
+(* ************************** *)
+M$GaugeGroups = {
+  U1Y  == { 
+    Abelian          -> True,  
+    CouplingConstant -> g1, 
+    GaugeBoson       -> B, 
+    Charge           -> Y
+  },
+  SU2L == { 
+    Abelian           -> False, 
+    CouplingConstant  -> gw, 
+    GaugeBoson        -> Wi, 
+    StructureConstant -> Eps, 
+    Representations   -> {Ta,SU2D}, 
+    Definitions       -> {Ta[a_,b_,c_]->PauliSigma[a,b,c]/2, FSU2L[i_,j_,k_]:> I Eps[i,j,k]}
+  },
+  SU3C == { 
+    Abelian           -> False, 
+    CouplingConstant  -> gs, 
+    GaugeBoson        -> G,
+    StructureConstant -> f, 
+    Representations   -> {T,Colour}, 
+    SymmetricTensor   -> dSUN
+  } 
+};
+
+
+(* ************************** *)
+(* *****    Indices     ***** *)
+(* ************************** *)
+
+IndexRange[Index[SU2W      ]] = Unfold[Range[3]]; 
+IndexRange[Index[SU2D      ]] = Unfold[Range[2]];
+IndexRange[Index[Gluon     ]] = NoUnfold[Range[8]];
+IndexRange[Index[Colour    ]] = NoUnfold[Range[3]]; 
+IndexRange[Index[Generation]] = Range[3];
+
+IndexStyle[SU2W,       j];
+IndexStyle[SU2D,       k];
+IndexStyle[Gluon,      a];
+IndexStyle[Colour,     m];
+IndexStyle[Generation, f];
+
+
+(* ************************** *)
+(* *** Interaction orders *** *)
+(* ***  (as used by mg5)  *** *)
+(* ************************** *)
+
+M$InteractionOrderHierarchy = {
+  {QCD, 1},
+  {QED, 2}
+};
+
+
+(* ************************** *)
+(* **** Particle classes **** *)
+(* ************************** *)
+M$ClassesDescription = {
+
+(* Gauge bosons: physical vector fields *)
+  V[1] == { 
+    ClassName       -> A, 
+    SelfConjugate   -> True,  
+    Mass            -> 0,  
+    Width           -> 0,  
+    ParticleName    -> "a", 
+    PDG             -> 22, 
+    PropagatorLabel -> "a", 
+    PropagatorType  -> W, 
+    PropagatorArrow -> None,
+    FullName        -> "Photon"
+  },
+  V[2] == { 
+    ClassName       -> Z, 
+    SelfConjugate   -> True,
+    Mass            -> {MZ, 91.1876},
+    Width           -> {WZ, 2.4952},
+    ParticleName    -> "Z", 
+    PDG             -> 23, 
+    PropagatorLabel -> "Z",
+    PropagatorType  -> Sine,
+    PropagatorArrow -> None,
+    FullName        -> "Z"
+  },
+  V[3] == {
+    ClassName        -> W,
+    SelfConjugate    -> False,
+    Mass             -> {MW, Internal},
+    Width            -> {WW, 2.085},
+    ParticleName     -> "W+",
+    AntiParticleName -> "W-",
+    QuantumNumbers   -> {Q -> 1},
+    PDG              -> 24, 
+    PropagatorLabel  -> "W",
+    PropagatorType   -> Sine,
+    PropagatorArrow  -> Forward,
+    FullName         -> "W"
+  },
+  V[4] == {
+    ClassName        -> G,
+    SelfConjugate    -> True,
+    Indices          -> {Index[Gluon]},
+    Mass             -> 0,
+    Width            -> 0,
+    ParticleName     -> "g", 
+    PDG              -> 21,
+    PropagatorLabel  -> "G",
+    PropagatorType   -> C,
+    PropagatorArrow  -> None,
+    FullName         -> "G"
+  },
+
+(* Ghosts: related to physical gauge bosons *)
+  U[1] == { 
+    ClassName       -> ghA, 
+    SelfConjugate   -> False,
+    Ghost           -> A,
+    QuantumNumbers  -> {GhostNumber -> 1},
+    Mass            -> 0,
+    Width	    -> 0,
+    PropagatorLabel -> "uA",
+    PropagatorType  -> GhostDash,
+    PropagatorArrow -> Forward
+  },
+  U[2] == {
+    ClassName       -> ghZ,
+    SelfConjugate   -> False,
+    Ghost           -> Z,
+    QuantumNumbers  -> {GhostNumber -> 1},
+    Mass            -> {MZ,91.1876},  
+    Width	    -> {WZ, 2.4952},
+    PropagatorLabel -> "uZ",
+    PropagatorType  -> GhostDash,
+    PropagatorArrow -> Forward
+  },
+  U[31] == { 
+    ClassName       -> ghWp,
+    SelfConjugate   -> False, 
+    Ghost           -> W,
+    QuantumNumbers  -> {GhostNumber -> 1, Q -> 1},
+    Mass            -> {MW,Internal}, 
+    Width           -> {WW, 2.085}, 
+    PropagatorLabel -> "uWp",
+    PropagatorType  -> GhostDash, 
+    PropagatorArrow -> Forward
+  },
+  U[32] == { 
+    ClassName       -> ghWm,
+    SelfConjugate   -> False, 
+    Ghost           -> Wbar,
+    QuantumNumbers  -> {GhostNumber -> 1, Q -> -1},
+    Mass            -> {MW,Internal}, 
+    Width           -> {WW, 2.085},
+    PropagatorLabel -> "uWm",
+    PropagatorType  -> GhostDash, 
+    PropagatorArrow -> Forward
+  },
+  U[4] == { 
+    ClassName       -> ghG, 
+    SelfConjugate   -> False,
+    Indices         -> {Index[Gluon]},
+    Ghost           -> G,
+    QuantumNumbers  ->{GhostNumber -> 1}, 
+    Mass            -> 0,
+    Width	    -> 0,
+    PropagatorLabel -> "uG",
+    PropagatorType  -> GhostDash,
+    PropagatorArrow -> Forward
+  },
+
+(* Gauge bosons: unphysical vector fields *)
+  V[11] == { 
+    ClassName     -> B, 
+    Unphysical    -> True, 
+    SelfConjugate -> True, 
+    Definitions   -> { B[mu_] -> -sw Z[mu]+cw A[mu]} 
+  },
+  V[12] == { 
+    ClassName     -> Wi,
+    Unphysical    -> True,
+    SelfConjugate -> True, 
+    Indices       -> {Index[SU2W]},
+    FlavorIndex   -> SU2W,
+    Definitions   -> { Wi[mu_,1] -> (Wbar[mu]+W[mu])/Sqrt[2], Wi[mu_,2] -> (Wbar[mu]-W[mu])/(I*Sqrt[2]), Wi[mu_,3] -> cw Z[mu] + sw A[mu]}
+  },
+
+(* Ghosts: related to unphysical gauge bosons *)
+  U[11] == {
+    ClassName     -> ghB, 
+    Unphysical    -> True,
+    SelfConjugate -> False,
+    Ghost         -> B, 
+    Definitions   -> { ghB -> -sw ghZ + cw ghA}
+  },
+  U[12] == {
+    ClassName     -> ghWi,
+    Unphysical    -> True,
+    SelfConjugate -> False,
+    Ghost         -> Wi,
+    Indices       -> {Index[SU2W]},
+    FlavorIndex   -> SU2W,
+    Definitions   -> { ghWi[1] -> (ghWp+ghWm)/Sqrt[2], ghWi[2] -> (ghWm-ghWp)/(I*Sqrt[2]), ghWi[3] -> cw ghZ+sw ghA}
+  } ,
+
+(* Fermions: physical fields *)
+  F[1] == {
+    ClassName        -> vl,
+    ClassMembers     -> {ve,vm,vt},
+    Indices          -> {Index[Generation]},
+    FlavorIndex      -> Generation,
+    SelfConjugate    -> False,
+    Mass             -> 0,
+    Width            -> 0,
+    QuantumNumbers   -> {LeptonNumber -> 1},
+    PropagatorLabel  -> {"v", "ve", "vm", "vt"} ,
+    PropagatorType   -> S,
+    PropagatorArrow  -> Forward,
+    PDG              -> {12,14,16},
+    ParticleName     -> {"ve","vm","vt"},
+    AntiParticleName -> {"ve~","vm~","vt~"},
+    FullName         -> {"Electron-neutrino", "Mu-neutrino", "Tau-neutrino"}
+  },
+  F[2] == {
+    ClassName        -> l,
+    ClassMembers     -> {e, mu, ta},
+    Indices          -> {Index[Generation]},
+    FlavorIndex      -> Generation,
+    SelfConjugate    -> False,
+    Mass             -> {Ml, {Me,5.11*^-4}, {MMU,0.10566}, {MTA,1.777}},
+    Width            -> 0,
+    QuantumNumbers   -> {Q -> -1, LeptonNumber -> 1},
+    PropagatorLabel  -> {"l", "e", "mu", "ta"},
+    PropagatorType   -> Straight,
+    PropagatorArrow  -> Forward,
+    PDG              -> {11, 13, 15},
+    ParticleName     -> {"e-", "mu-", "ta-"},
+    AntiParticleName -> {"e+", "mu+", "ta+"},
+    FullName         -> {"Electron", "Muon", "Tau"} 
+  },
+  F[3] == {
+    ClassName        -> uq,
+    ClassMembers     -> {u, c, t},
+    Indices          -> {Index[Generation], Index[Colour]},
+    FlavorIndex      -> Generation,
+    SelfConjugate    -> False,
+    Mass             -> {Mu, {MU, 2.55*^-3}, {MC,1.27}, {MT,172}},
+    Width            -> {0, 0, {WT,1.50833649}},
+    QuantumNumbers   -> {Q -> 2/3},
+    PropagatorLabel  -> {"uq", "u", "c", "t"},
+    PropagatorType   -> Straight,
+    PropagatorArrow  -> Forward,
+    PDG              -> {2, 4, 6}, 
+    ParticleName     -> {"u",  "c",  "t" },
+    AntiParticleName -> {"u~", "c~", "t~"},
+    FullName         -> {"u-quark", "c-quark", "t-quark"}
+  },
+  F[4] == {
+    ClassName        -> dq,
+    ClassMembers     -> {d, s, b},
+    Indices          -> {Index[Generation], Index[Colour]},
+    FlavorIndex      -> Generation,
+    SelfConjugate    -> False,
+    Mass             -> {Md, {MD,5.04*^-3}, {MS,0.101}, {MB,4.7}},
+    Width            -> 0,
+    QuantumNumbers   -> {Q -> -1/3},
+    PropagatorLabel  -> {"dq", "d", "s", "b"},
+    PropagatorType   -> Straight, 
+    PropagatorArrow  -> Forward,
+    PDG              -> {1,3,5},
+    ParticleName     -> {"d",  "s",  "b" },
+    AntiParticleName -> {"d~", "s~", "b~"},
+    FullName         -> {"d-quark", "s-quark", "b-quark"}
+  },
+
+(* Fermions: unphysical fields *)
+  F[11] == { 
+    ClassName      -> LL, 
+    Unphysical     -> True, 
+    Indices        -> {Index[SU2D], Index[Generation]},
+    FlavorIndex    -> SU2D,
+    SelfConjugate  -> False,
+    QuantumNumbers -> {Y -> -1/2},
+    Definitions    -> { LL[sp1_,1,ff_] :> Module[{sp2}, ProjM[sp1,sp2] vl[sp2,ff]], LL[sp1_,2,ff_] :> Module[{sp2}, ProjM[sp1,sp2] l[sp2,ff]] }
+  },
+  F[12] == { 
+    ClassName      -> lR, 
+    Unphysical     -> True, 
+    Indices        -> {Index[Generation]},
+    FlavorIndex    -> Generation,
+    SelfConjugate  -> False,
+    QuantumNumbers -> {Y -> -1},
+    Definitions    -> { lR[sp1_,ff_] :> Module[{sp2}, ProjP[sp1,sp2] l[sp2,ff]] }
+  },
+  F[13] == { 
+    ClassName      -> QL, 
+    Unphysical     -> True, 
+    Indices        -> {Index[SU2D], Index[Generation], Index[Colour]},
+    FlavorIndex    -> SU2D,
+    SelfConjugate  -> False,
+    QuantumNumbers -> {Y -> 1/6},
+    Definitions    -> { 
+      QL[sp1_,1,ff_,cc_] :> Module[{sp2}, ProjM[sp1,sp2] uq[sp2,ff,cc]], 
+      QL[sp1_,2,ff_,cc_] :> Module[{sp2,ff2}, CKM[ff,ff2] ProjM[sp1,sp2] dq[sp2,ff2,cc]] }
+  },
+  F[14] == { 
+    ClassName      -> uR, 
+    Unphysical     -> True, 
+    Indices        -> {Index[Generation], Index[Colour]},
+    FlavorIndex    -> Generation,
+    SelfConjugate  -> False,
+    QuantumNumbers -> {Y -> 2/3},
+    Definitions    -> { uR[sp1_,ff_,cc_] :> Module[{sp2}, ProjP[sp1,sp2] uq[sp2,ff,cc]] }
+  },
+  F[15] == { 
+    ClassName      -> dR, 
+    Unphysical     -> True, 
+    Indices        -> {Index[Generation], Index[Colour]},
+    FlavorIndex    -> Generation,
+    SelfConjugate  -> False,
+    QuantumNumbers -> {Y -> -1/3},
+    Definitions    -> { dR[sp1_,ff_,cc_] :> Module[{sp2}, ProjP[sp1,sp2] dq[sp2,ff,cc]] }
+  },
+
+(* Higgs: physical scalars  *)
+  S[1] == {
+    ClassName       -> H,
+    SelfConjugate   -> True,
+    Mass            -> {MH,125},
+    Width           -> {WH,0.00407},
+    PropagatorLabel -> "H",
+    PropagatorType  -> D,
+    PropagatorArrow -> None,
+    PDG             -> 25,
+    ParticleName    -> "H",
+    FullName        -> "H"
+  },
+
+(* Higgs: physical scalars  *)
+  S[2] == {
+    ClassName       -> G0,
+    SelfConjugate   -> True,
+    Goldstone       -> Z,
+    Mass            -> {MZ, 91.1876},
+    Width           -> {WZ, 2.4952},
+    PropagatorLabel -> "Go",
+    PropagatorType  -> D,
+    PropagatorArrow -> None,
+    PDG             -> 250,
+    ParticleName    -> "G0",
+    FullName        -> "G0"
+  },
+  S[3] == {
+    ClassName        -> GP,
+    SelfConjugate    -> False,
+    Goldstone        -> W,
+    Mass             -> {MW, Internal},
+    QuantumNumbers   -> {Q -> 1},
+    Width            -> {WW, 2.085},
+    PropagatorLabel  -> "GP",
+    PropagatorType   -> D,
+    PropagatorArrow  -> None,
+    PDG              -> 251,
+    ParticleName     -> "G+",
+    AntiParticleName -> "G-",
+    FullName         -> "GP"
+  },
+
+(* Higgs: unphysical scalars  *)
+  S[11] == { 
+    ClassName      -> Phi, 
+    Unphysical     -> True, 
+    Indices        -> {Index[SU2D]},
+    FlavorIndex    -> SU2D,
+    SelfConjugate  -> False,
+    QuantumNumbers -> {Y -> 1/2},
+    Definitions    -> { Phi[1] -> -I GP, Phi[2] -> (vev + H + I G0)/Sqrt[2]  }
+  }
+};
+
+
+(* ************************** *)
+(* *****     Gauge      ***** *)
+(* *****   Parameters   ***** *)
+(* *****   (FeynArts)   ***** *)
+(* ************************** *)
+
+GaugeXi[ V[1]  ] = GaugeXi[A];
+GaugeXi[ V[2]  ] = GaugeXi[Z];
+GaugeXi[ V[3]  ] = GaugeXi[W];
+GaugeXi[ V[4]  ] = GaugeXi[G];
+GaugeXi[ S[1]  ] = 1;
+GaugeXi[ S[2]  ] = GaugeXi[Z];
+GaugeXi[ S[3]  ] = GaugeXi[W];
+GaugeXi[ U[1]  ] = GaugeXi[A];
+GaugeXi[ U[2]  ] = GaugeXi[Z];
+GaugeXi[ U[31] ] = GaugeXi[W];
+GaugeXi[ U[32] ] = GaugeXi[W];
+GaugeXi[ U[4]  ] = GaugeXi[G];
+
+
+(* ************************** *)
+(* *****   Parameters   ***** *)
+(* ************************** *)
+M$Parameters = {
+
+  (* External parameters *)
+  aEWM1 == { 
+    ParameterType    -> External, 
+    BlockName        -> SMINPUTS, 
+    OrderBlock       -> 1, 
+    Value            -> 127.9,
+    InteractionOrder -> {QED,-2},
+    Description      -> "Inverse of the EW coupling constant at the Z pole"
+  },
+  Gf == {
+    ParameterType    -> External,
+    BlockName        -> SMINPUTS,
+    OrderBlock       -> 2,
+    Value            -> 1.16637*^-5, 
+    InteractionOrder -> {QED,2},
+    TeX              -> Subscript[G,f],
+    Description      -> "Fermi constant"
+  },
+  aS    == { 
+    ParameterType    -> External,
+    BlockName        -> SMINPUTS,
+    OrderBlock       -> 3,
+    Value            -> 0.1184, 
+    InteractionOrder -> {QCD,2},
+    TeX              -> Subscript[\[Alpha],s],
+    Description      -> "Strong coupling constant at the Z pole"
+  },
+  ymdo == {
+    ParameterType -> External,
+    BlockName     -> YUKAWA,
+    OrderBlock    -> 1,
+    Value         -> 5.04*^-3,
+    Description   -> "Down Yukawa mass"
+  },
+  ymup == {
+    ParameterType -> External,
+    BlockName     -> YUKAWA,
+    OrderBlock    -> 2,
+    Value         -> 2.55*^-3,
+    Description   -> "Up Yukawa mass"
+  },
+  yms == {
+    ParameterType -> External,
+    BlockName     -> YUKAWA,
+    OrderBlock    -> 3,
+    Value         -> 0.101,
+    Description   -> "Strange Yukawa mass"
+  },
+  ymc == {
+    ParameterType -> External,
+    BlockName     -> YUKAWA,
+    OrderBlock    -> 4,
+    Value         -> 1.27,
+    Description   -> "Charm Yukawa mass"
+  },
+  ymb == {
+    ParameterType -> External,
+    BlockName     -> YUKAWA,
+    OrderBlock    -> 5,
+    Value         -> 4.7,
+    Description   -> "Bottom Yukawa mass"
+  },
+  ymt == {
+    ParameterType -> External,
+    BlockName     -> YUKAWA,
+    OrderBlock    -> 6,
+    Value         -> 172,
+    Description   -> "Top Yukawa mass"
+  },
+  yme == {
+    ParameterType -> External,
+    BlockName     -> YUKAWA,
+    OrderBlock    -> 11,
+    Value         -> 5.11*^-4,
+    Description   -> "Electron Yukawa mass"
+  },
+  ymm == {
+    ParameterType -> External,
+    BlockName     -> YUKAWA,
+    OrderBlock    -> 13,
+    Value         -> 0.10566,
+    Description   -> "Muon Yukawa mass"
+  },
+  ymtau == {
+    ParameterType -> External,
+    BlockName     -> YUKAWA,
+    OrderBlock    -> 15,
+    Value         -> 1.777,
+    Description   -> "Tau Yukawa mass"
+  },
+  cabi == {
+    ParameterType -> External,
+    BlockName     -> CKMBLOCK,
+    OrderBlock    -> 1,
+    Value         -> 0.227736,
+    TeX           -> Subscript[\[Theta], c],
+    Description   -> "Cabibbo angle"
+  },
+
+  (* Internal Parameters *)
+  aEW == {
+    ParameterType    -> Internal,
+    Value            -> 1/aEWM1,
+    InteractionOrder -> {QED,2},
+    TeX              -> Subscript[\[Alpha], EW],
+    Description      -> "Electroweak coupling contant"
+  },
+  MW == { 
+    ParameterType -> Internal, 
+    Value         -> Sqrt[MZ^2/2+Sqrt[MZ^4/4-Pi/Sqrt[2]*aEW/Gf*MZ^2]], 
+    TeX           -> Subscript[M,W], 
+    Description   -> "W mass"
+  },
+  sw2 == { 
+    ParameterType -> Internal, 
+    Value         -> 1-(MW/MZ)^2, 
+    Description   -> "Squared Sin of the Weinberg angle"
+  },
+  ee == { 
+    ParameterType    -> Internal, 
+    Value            -> Sqrt[4 Pi aEW], 
+    InteractionOrder -> {QED,1}, 
+    TeX              -> e,  
+    Description      -> "Electric coupling constant"
+  },
+  cw == { 
+    ParameterType -> Internal, 
+    Value         -> Sqrt[1-sw2], 
+    TeX           -> Subscript[c,w], 
+    Description   -> "Cosine of the Weinberg angle"
+  },
+  sw == { 
+    ParameterType -> Internal, 
+    Value         -> Sqrt[sw2], 
+    TeX           -> Subscript[s,w], 
+    Description   -> "Sine of the Weinberg angle"
+  },
+  gw == { 
+    ParameterType    -> Internal, 
+    Definitions      -> {gw->ee/sw}, 
+    InteractionOrder -> {QED,1},  
+    TeX              -> Subscript[g,w], 
+    Description      -> "Weak coupling constant at the Z pole"
+  },
+  g1 == { 
+    ParameterType    -> Internal, 
+    Definitions      -> {g1->ee/cw}, 
+    InteractionOrder -> {QED,1},  
+    TeX              -> Subscript[g,1], 
+    Description      -> "U(1)Y coupling constant at the Z pole"
+  },
+  gs == { 
+    ParameterType    -> Internal, 
+    Value            -> Sqrt[4 Pi aS],
+    InteractionOrder -> {QCD,1},  
+    TeX              -> Subscript[g,s], 
+    ParameterName    -> G,
+    Description      -> "Strong coupling constant at the Z pole"
+  },
+  vev == {
+    ParameterType    -> Internal,
+    Value            -> 2*MW*sw/ee, 
+    InteractionOrder -> {QED,-1},
+    Description      -> "Higgs vacuum expectation value"
+  },
+  lam == {
+    ParameterType    -> Internal,
+    Value            -> MH^2/(2*vev^2),
+    InteractionOrder -> {QED, 2},
+    Description      -> "Higgs quartic coupling"
+  },
+  muH == {
+    ParameterType -> Internal,
+    Value         -> Sqrt[vev^2 lam],
+    TeX           -> \[Mu],
+    Description   -> "Coefficient of the quadratic piece of the Higgs potential"
+  },
+  yl == {
+    ParameterType    -> Internal,
+    Indices          -> {Index[Generation], Index[Generation]},
+    Definitions      -> {yl[i_?NumericQ, j_?NumericQ] :> 0  /; (i =!= j)},
+    Value            -> {yl[1,1] -> Sqrt[2] yme / vev, yl[2,2] -> Sqrt[2] ymm / vev, yl[3,3] -> Sqrt[2] ymtau / vev},
+    InteractionOrder -> {QED, 1},
+    ParameterName    -> {yl[1,1] -> ye, yl[2,2] -> ym, yl[3,3] -> ytau},
+    TeX              -> Superscript[y, l],
+    Description      -> "Lepton Yukawa couplings"
+  },
+  yu == {
+    ParameterType    -> Internal,
+    Indices          -> {Index[Generation], Index[Generation]},
+    Definitions      -> {yu[i_?NumericQ, j_?NumericQ] :> 0  /; (i =!= j)},
+    Value            -> {yu[1,1] -> Sqrt[2] ymup/vev, yu[2,2] -> Sqrt[2] ymc/vev, yu[3,3] -> Sqrt[2] ymt/vev},
+    InteractionOrder -> {QED, 1},
+    ParameterName    -> {yu[1,1] -> yup, yu[2,2] -> yc, yu[3,3] -> yt},
+    TeX              -> Superscript[y, u],
+    Description      -> "Up-type Yukawa couplings"
+  },
+  yd == {
+    ParameterType    -> Internal,
+    Indices          -> {Index[Generation], Index[Generation]},
+    Definitions      -> {yd[i_?NumericQ, j_?NumericQ] :> 0  /; (i =!= j)},
+    Value            -> {yd[1,1] -> Sqrt[2] ymdo/vev, yd[2,2] -> Sqrt[2] yms/vev, yd[3,3] -> Sqrt[2] ymb/vev},
+    InteractionOrder -> {QED, 1},
+    ParameterName    -> {yd[1,1] -> ydo, yd[2,2] -> ys, yd[3,3] -> yb},
+    TeX              -> Superscript[y, d],
+    Description      -> "Down-type Yukawa couplings"
+  },
+(* N. B. : only Cabibbo mixing! *)
+  CKM == { 
+    ParameterType -> Internal,
+    Indices       -> {Index[Generation], Index[Generation]},
+    Unitary       -> True,
+    Value         -> {CKM[1,1] -> Cos[cabi],  CKM[1,2] -> Sin[cabi], CKM[1,3] -> 0,
+                      CKM[2,1] -> -Sin[cabi], CKM[2,2] -> Cos[cabi], CKM[2,3] -> 0,
+                      CKM[3,1] -> 0,          CKM[3,2] -> 0,         CKM[3,3] -> 1},
+    TeX         -> Superscript[V,CKM],
+    Description -> "CKM-Matrix"}
+};
+
+(* ************************** *)
+(* *****   Lagrangian   ***** *)
+(* ************************** *)
+
+LGauge := Block[{mu,nu,ii,aa}, 
+  ExpandIndices[-1/4 FS[B,mu,nu] FS[B,mu,nu] - 1/4 FS[Wi,mu,nu,ii] FS[Wi,mu,nu,ii] - 1/4 FS[G,mu,nu,aa] FS[G,mu,nu,aa], FlavorExpand->SU2W]];
+
+LFermions := Block[{mu}, 
+  ExpandIndices[I*(
+    QLbar.Ga[mu].DC[QL, mu] + LLbar.Ga[mu].DC[LL, mu] + uRbar.Ga[mu].DC[uR, mu] + dRbar.Ga[mu].DC[dR, mu] + lRbar.Ga[mu].DC[lR, mu]), 
+  FlavorExpand->{SU2W,SU2D}]/.{CKM[a_,b_] Conjugate[CKM[a_,c_]]->IndexDelta[b,c], CKM[b_,a_] Conjugate[CKM[c_,a_]]->IndexDelta[b,c]}];
+
+LHiggs := Block[{ii,mu, feynmangaugerules},
+  feynmangaugerules = If[Not[FeynmanGauge], {G0|GP|GPbar ->0}, {}];
+ 
+  ExpandIndices[DC[Phibar[ii],mu] DC[Phi[ii],mu] + muH^2 Phibar[ii] Phi[ii] - lam Phibar[ii] Phi[ii] Phibar[jj] Phi[jj], FlavorExpand->{SU2D,SU2W}]/.feynmangaugerules
+ ];
+
+LYukawa := Block[{sp,ii,jj,cc,ff1,ff2,ff3,yuk,feynmangaugerules},
+  feynmangaugerules = If[Not[FeynmanGauge], {G0|GP|GPbar ->0}, {}];
+ 
+  yuk = ExpandIndices[
+   -yd[ff2, ff3] CKM[ff1, ff2] QLbar[sp, ii, ff1, cc].dR [sp, ff3, cc] Phi[ii] - 
+    yl[ff1, ff3] LLbar[sp, ii, ff1].lR [sp, ff3] Phi[ii] - 
+    yu[ff1, ff2] QLbar[sp, ii, ff1, cc].uR [sp, ff2, cc] Phibar[jj] Eps[ii, jj], FlavorExpand -> SU2D];
+  yuk = yuk /. { CKM[a_, b_] Conjugate[CKM[a_, c_]] -> IndexDelta[b, c], CKM[b_, a_] Conjugate[CKM[c_, a_]] -> IndexDelta[b, c]};
+  yuk+HC[yuk]/.feynmangaugerules
+ ];
+
+LGhost := Block[{LGh1,LGhw,LGhs,LGhphi,mu, generators,gh,ghbar,Vectorize,phi1,phi2,togoldstones,doublet,doublet0},
+  (* Pure gauge piece *) 	
+  LGh1 = -ghBbar.del[DC[ghB,mu],mu];
+  LGhw = -ghWibar.del[DC[ghWi,mu],mu];
+  LGhs = -ghGbar.del[DC[ghG,mu],mu];
+
+  (* Scalar pieces: see Peskin pages 739-742 *)
+  (* phi1 and phi2 are the real degrees of freedom of GP *)
+  (* Vectorize transforms a doublet in a vector in the phi-basis, i.e. the basis of real degrees of freedom *)
+  gh    = {ghB, ghWi[1], ghWi[2], ghWi[3]};
+  ghbar = {ghBbar, ghWibar[1], ghWibar[2], ghWibar[3]};
+  generators = {-I/2 g1 IdentityMatrix[2], -I/2 gw PauliSigma[1], -I/2 gw PauliSigma[2], -I/2 gw PauliSigma[3]};
+  doublet = Expand[{(-I phi1 - phi2)/Sqrt[2], Phi[2]} /. MR$Definitions /. vev -> 0]; 
+  doublet0 = {0, vev/Sqrt[2]};
+  Vectorize[{a_, b_}]:= Simplify[{Sqrt[2] Re[Expand[a]], Sqrt[2] Im[Expand[a]], Sqrt[2] Re[Expand[b]], Sqrt[2] Im[Expand[b]]}/.{Im[_]->0, Re[num_]->num}];
+  togoldstones := {phi1 -> (GP + GPbar)/Sqrt[2], phi2 -> (-GP + GPbar)/(I Sqrt[2])};
+  LGhphi=Plus@@Flatten[Table[-ghbar[[kkk]].gh[[lll]] Vectorize[generators[[kkk]].doublet0].Vectorize[generators[[lll]].(doublet+doublet0)],{kkk,4},{lll,4}]] /.togoldstones;
+
+ExpandIndices[ LGhs + If[FeynmanGauge, LGh1 + LGhw + LGhphi,0], FlavorExpand->SU2W]];
+
+LSM:= LGauge + LFermions + LHiggs + LYukawa + LGhost;
diff --git a/CHresonance_neutral_scalar_UFO/feynrules/neutral_scalar.fr b/CHresonance_neutral_scalar_UFO/feynrules/neutral_scalar.fr
new file mode 100644
index 0000000000000000000000000000000000000000..4883f4e3b09f7f99351a7712b0393f36d0780132
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/feynrules/neutral_scalar.fr
@@ -0,0 +1,40 @@
+(* ************************** *)
+(* *****  Information   ***** *)
+(* ************************** *)
+
+M$ModelName = "CHresonance_neutral_scalar";
+
+(* M$InteractionOrderHierarchy = { *)
+(*   {NP , 1},                     *)
+(*   {QCD, 1},                     *)
+(*   {QED, 2}                      *)
+(* };                              *)
+
+M$ClassesDescription = {
+  S[4] == {
+    ClassName     -> s0,
+    SelfConjugate -> True,
+    Mass          -> {Ms0, 1000.0},
+    Width         -> {Ws0, 1.0},
+    ParticleName  -> "s0"}
+}
+
+M$Parameters = {
+  keta == {
+    ParameterType    -> External,
+    Value            -> 1.0,
+    InteractionOrder -> {NP,1},
+    Description      -> "Scalar cubic coupling times (v/f)"},
+  gweak == {
+    ParameterType    -> Internal,
+    Value            -> gw,
+    InteractionOrder -> {QED,0},
+    Description      -> "gw at QED=0 order"},
+  vweak == {
+    ParameterType    -> Internal,
+    Value            -> vev,
+    InteractionOrder -> {QED,0},
+    Description      -> "vev at QED=0 order"}
+}
+
+Lscalar := 1/2 del[s0, mu] del[s0, mu] - 1/2 Ms0^2 s0^2 + keta gweak MW W[mu] Wbar[mu] s0 + keta gweak (1/2) (MZ/cw) Z[mu] Z[mu] s0 + keta (1/vweak) del[H, mu] del[H, mu] s0
diff --git a/CHresonance_neutral_scalar_UFO/feynrules/neutral_scalar.tar.gz b/CHresonance_neutral_scalar_UFO/feynrules/neutral_scalar.tar.gz
new file mode 100644
index 0000000000000000000000000000000000000000..e2e183f1cb05eeb1484cd679d7e8c31d0b47b252
Binary files /dev/null and b/CHresonance_neutral_scalar_UFO/feynrules/neutral_scalar.tar.gz differ
diff --git a/CHresonance_neutral_scalar_UFO/function_library.py b/CHresonance_neutral_scalar_UFO/function_library.py
new file mode 100644
index 0000000000000000000000000000000000000000..c4cfad7f0e64cbe3eab938f53dc84aeda4d1e26a
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/function_library.py
@@ -0,0 +1,71 @@
+# 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')
+
+# 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)')
+
+cot = Function(name = 'cot',
+               arguments = ('z',),
+               expression = '1./cmath.tan(z)')
+
+# Heaviside theta function
+
+theta_function = Function(name = 'theta_function',
+             arguments = ('x','y','z'),
+             expression = 'y if x else z')
+
+# 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))')
+
diff --git a/CHresonance_neutral_scalar_UFO/lorentz.py b/CHresonance_neutral_scalar_UFO/lorentz.py
new file mode 100644
index 0000000000000000000000000000000000000000..8f782603acc504dd24a722911e2f6bed4c664c21
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/lorentz.py
@@ -0,0 +1,106 @@
+# This file was automatically created by FeynRules 2.0.2
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (February 23, 2011)
+# Date: Wed 8 Oct 2014 20:39:09
+
+
+from object_library import all_lorentz, Lorentz
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec, cot
+try:
+   import form_factors as ForFac 
+except ImportError:
+   pass
+
+
+UUS1 = Lorentz(name = 'UUS1',
+               spins = [ -1, -1, 1 ],
+               structure = '1')
+
+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')
+
+SSS2 = Lorentz(name = 'SSS2',
+               spins = [ 1, 1, 1 ],
+               structure = 'P(-1,1)*P(-1,2)')
+
+FFS1 = Lorentz(name = 'FFS1',
+               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)')
+
+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)')
+
+VSS1 = Lorentz(name = 'VSS1',
+               spins = [ 3, 1, 1 ],
+               structure = 'P(1,2) - P(1,3)')
+
+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)')
+
+SSSS1 = Lorentz(name = 'SSSS1',
+                spins = [ 1, 1, 1, 1 ],
+                structure = '1')
+
+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.')
+
diff --git a/CHresonance_neutral_scalar_UFO/object_library.py b/CHresonance_neutral_scalar_UFO/object_library.py
new file mode 100644
index 0000000000000000000000000000000000000000..12278b669a6d95215c03269045a0e58fb70335c6
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/object_library.py
@@ -0,0 +1,377 @@
+##
+##
+## 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','counterterm','charge', 'line', 'propagating', 'goldstoneboson', 'propagator']
+
+    def __init__(self, pdg_code, name, antiname, spin, color, mass, width, texname,
+                 antitexname, charge , line=None, propagating=True, counterterm=None, goldstoneboson=False, 
+                 propagator=None, **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 not line:                                                                                                                                                                                   
+            self.line = self.find_line_type()
+        else:
+            self.line = line
+
+        if propagator:
+            if isinstance(propagator, dict):
+                self.propagator = propagator
+            else:
+                self.propagator = {0: propagator, 1: propagator}
+             
+    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):
+        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):
+
+        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, **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'
+
+        CTparam=None
+        for param in all_CTparameters:
+           pattern=re.compile(r"(?P<first>\A|\*|\+|\-|\()(?P<name>"+param.name+r")(?P<second>\Z|\*|\+|\-|\))")
+           numberOfMatches=len(pattern.findall(self.value))
+           if numberOfMatches==1:
+               if not CTparam:
+                   CTparam=param
+               else:
+                   raise UFOError("UFO does not support yet more than one occurence of CTParameters in the couplings values.")
+           elif numberOfMatches>1:
+               raise UFOError("UFO does not support yet more than one occurence of CTParameters in the couplings values.")
+
+        if not CTparam:
+            if x==0:
+                return self.value
+            else:
+                return 'ZERO'
+        else:
+            if CTparam.pole(x)=='ZERO':
+                return 'ZERO'
+            else:
+                def substitution(matchedObj):
+                    return matchedObj.group('first')+"("+CTparam.pole(x)+")"+matchedObj.group('second')
+                pattern=re.compile(r"(?P<first>\A|\*|\+|\-|\()(?P<name>"+CTparam.name+r")(?P<second>\Z|\*|\+|\-|\))")
+                return pattern.sub(substitution,self.value)
+
+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
+
+all_decays = []
+
+class Decay(UFOBaseClass):
+    require_args = ['particle','partial_widths']
+
+    def __init__(self, particle, partial_widths, **opt):
+        args = (particle, partial_widths)
+        UFOBaseClass.__init__(self, *args, **opt)
+
+        global all_decays
+        all_decays.append(self)
+    
+        # Add the information directly to the particle
+        particle.partial_widths = partial_widths
+
+all_form_factors = []
+
+class FormFactor(UFOBaseClass):
+    require_args = ['name','type','value']
+
+    def __init__(self, name, type, value, **opt):
+        args = (name, type, value)
+        UFOBaseClass.__init__(self, *args, **opt)
+
+        global all_form_factors
+        all_form_factors.append(self)
+
+        
+all_propagators = []
+
+class Propagator(UFOBaseClass):
+    
+    require_args = ['name','numerator','denominator']
+
+    def __init__(self, name, numerator, denominator=None, **opt):
+        args = (name, numerator, denominator)
+        UFOBaseClass.__init__(self, *args, **opt)
+
+        global all_propagators
+        all_propagators.append(self)
diff --git a/CHresonance_neutral_scalar_UFO/parameters.py b/CHresonance_neutral_scalar_UFO/parameters.py
new file mode 100644
index 0000000000000000000000000000000000000000..509a4dd2414cb9d8858893c7e338bf1d40ef4bea
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/parameters.py
@@ -0,0 +1,290 @@
+# This file was automatically created by FeynRules 2.0.2
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (February 23, 2011)
+# Date: Wed 8 Oct 2014 20:39:09
+
+
+
+from object_library import all_parameters, Parameter
+
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec, cot
+
+# This is a default parameter object representing 0.
+ZERO = Parameter(name = 'ZERO',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '0.0',
+                 texname = '0')
+
+# User-defined parameters.
+aEWM1 = Parameter(name = 'aEWM1',
+                  nature = 'external',
+                  type = 'real',
+                  value = 127.9,
+                  texname = '\\text{aEWM1}',
+                  lhablock = 'SMINPUTS',
+                  lhacode = [ 1 ])
+
+Gf = Parameter(name = 'Gf',
+               nature = 'external',
+               type = 'real',
+               value = 0.0000116637,
+               texname = 'G_f',
+               lhablock = 'SMINPUTS',
+               lhacode = [ 2 ])
+
+aS = Parameter(name = 'aS',
+               nature = 'external',
+               type = 'real',
+               value = 0.1184,
+               texname = '\\alpha _s',
+               lhablock = 'SMINPUTS',
+               lhacode = [ 3 ])
+
+ymb = Parameter(name = 'ymb',
+                nature = 'external',
+                type = 'real',
+                value = 4.7,
+                texname = '\\text{ymb}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 5 ])
+
+ymt = Parameter(name = 'ymt',
+                nature = 'external',
+                type = 'real',
+                value = 172,
+                texname = '\\text{ymt}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 6 ])
+
+ymtau = Parameter(name = 'ymtau',
+                  nature = 'external',
+                  type = 'real',
+                  value = 1.777,
+                  texname = '\\text{ymtau}',
+                  lhablock = 'YUKAWA',
+                  lhacode = [ 15 ])
+
+keta = Parameter(name = 'keta',
+                 nature = 'external',
+                 type = 'real',
+                 value = 0.1474,
+                 texname = '\\text{keta}',
+                 lhablock = 'FRBlock',
+                 lhacode = [ 1 ])
+
+MZ = Parameter(name = 'MZ',
+               nature = 'external',
+               type = 'real',
+               value = 91.1876,
+               texname = '\\text{MZ}',
+               lhablock = 'MASS',
+               lhacode = [ 23 ])
+
+MTA = Parameter(name = 'MTA',
+                nature = 'external',
+                type = 'real',
+                value = 1.777,
+                texname = '\\text{MTA}',
+                lhablock = 'MASS',
+                lhacode = [ 15 ])
+
+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 ])
+
+MH = Parameter(name = 'MH',
+               nature = 'external',
+               type = 'real',
+               value = 125,
+               texname = '\\text{MH}',
+               lhablock = 'MASS',
+               lhacode = [ 25 ])
+
+Ms0 = Parameter(name = 'Ms0',
+                nature = 'external',
+                type = 'real',
+                value = 7000.0,
+                #value = 1000.0,
+                texname = '\\text{Ms0}',
+                lhablock = 'MASS',
+                lhacode = [ 9000006 ])
+
+WZ = Parameter(name = 'WZ',
+               nature = 'external',
+               type = 'real',
+               value = 2.4952,
+               texname = '\\text{WZ}',
+               lhablock = 'DECAY',
+               lhacode = [ 23 ])
+
+WW = Parameter(name = 'WW',
+               nature = 'external',
+               type = 'real',
+               value = 2.085,
+               texname = '\\text{WW}',
+               lhablock = 'DECAY',
+               lhacode = [ 24 ])
+
+WT = Parameter(name = 'WT',
+               nature = 'external',
+               type = 'real',
+               value = 1.50833649,
+               texname = '\\text{WT}',
+               lhablock = 'DECAY',
+               lhacode = [ 6 ])
+
+WH = Parameter(name = 'WH',
+               nature = 'external',
+               type = 'real',
+               value = 0.00407,
+               texname = '\\text{WH}',
+               lhablock = 'DECAY',
+               lhacode = [ 25 ])
+
+Ws0 = Parameter(name = 'Ws0',
+                nature = 'external',
+                type = 'real',
+                value = 1100.0,
+                texname = '\\text{Ws0}',
+                lhablock = 'DECAY',
+                lhacode = [ 9000006 ])
+
+aEW = Parameter(name = 'aEW',
+                nature = 'internal',
+                type = 'real',
+                value = '1/aEWM1',
+                texname = '\\alpha _{\\text{EW}}')
+
+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')
+
+vev = Parameter(name = 'vev',
+                nature = 'internal',
+                type = 'real',
+                value = '(2*MW*sw)/ee',
+                texname = '\\text{vev}')
+
+gweak = Parameter(name = 'gweak',
+                  nature = 'internal',
+                  type = 'real',
+                  value = 'gw',
+                  texname = '\\text{gweak}')
+
+lam = Parameter(name = 'lam',
+                nature = 'internal',
+                type = 'real',
+                value = 'MH**2/(2.*vev**2)',
+                texname = '\\text{lam}')
+
+vweak = Parameter(name = 'vweak',
+                  nature = 'internal',
+                  type = 'real',
+                  value = 'vev',
+                  texname = '\\text{vweak}')
+
+yb = Parameter(name = 'yb',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymb*cmath.sqrt(2))/vev',
+               texname = '\\text{yb}')
+
+yt = Parameter(name = 'yt',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymt*cmath.sqrt(2))/vev',
+               texname = '\\text{yt}')
+
+ytau = Parameter(name = 'ytau',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '(ymtau*cmath.sqrt(2))/vev',
+                 texname = '\\text{ytau}')
+
+muH = Parameter(name = 'muH',
+                nature = 'internal',
+                type = 'real',
+                value = 'cmath.sqrt(lam*vev**2)',
+                texname = '\\mu')
+
+I1b33 = Parameter(name = 'I1b33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yb',
+                  texname = '\\text{I1b33}')
+
+I2b33 = Parameter(name = 'I2b33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yt',
+                  texname = '\\text{I2b33}')
+
+I3b33 = Parameter(name = 'I3b33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yt',
+                  texname = '\\text{I3b33}')
+
+I4b33 = Parameter(name = 'I4b33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yb',
+                  texname = '\\text{I4b33}')
+
diff --git a/CHresonance_neutral_scalar_UFO/parametersTemplate.py b/CHresonance_neutral_scalar_UFO/parametersTemplate.py
new file mode 100644
index 0000000000000000000000000000000000000000..d9ffa1265efd85ef18d8ba718c46d1d2571aca2c
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/parametersTemplate.py
@@ -0,0 +1,289 @@
+# This file was automatically created by FeynRules 2.0.2
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (February 23, 2011)
+# Date: Wed 8 Oct 2014 20:39:09
+
+
+
+from object_library import all_parameters, Parameter
+
+
+from function_library import complexconjugate, re, im, csc, sec, acsc, asec, cot
+
+# This is a default parameter object representing 0.
+ZERO = Parameter(name = 'ZERO',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '0.0',
+                 texname = '0')
+
+# User-defined parameters.
+aEWM1 = Parameter(name = 'aEWM1',
+                  nature = 'external',
+                  type = 'real',
+                  value = 127.9,
+                  texname = '\\text{aEWM1}',
+                  lhablock = 'SMINPUTS',
+                  lhacode = [ 1 ])
+
+Gf = Parameter(name = 'Gf',
+               nature = 'external',
+               type = 'real',
+               value = 0.0000116637,
+               texname = 'G_f',
+               lhablock = 'SMINPUTS',
+               lhacode = [ 2 ])
+
+aS = Parameter(name = 'aS',
+               nature = 'external',
+               type = 'real',
+               value = 0.1184,
+               texname = '\\alpha _s',
+               lhablock = 'SMINPUTS',
+               lhacode = [ 3 ])
+
+ymb = Parameter(name = 'ymb',
+                nature = 'external',
+                type = 'real',
+                value = 4.7,
+                texname = '\\text{ymb}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 5 ])
+
+ymt = Parameter(name = 'ymt',
+                nature = 'external',
+                type = 'real',
+                value = 172,
+                texname = '\\text{ymt}',
+                lhablock = 'YUKAWA',
+                lhacode = [ 6 ])
+
+ymtau = Parameter(name = 'ymtau',
+                  nature = 'external',
+                  type = 'real',
+                  value = 1.777,
+                  texname = '\\text{ymtau}',
+                  lhablock = 'YUKAWA',
+                  lhacode = [ 15 ])
+
+keta = Parameter(name = 'keta',
+                 nature = 'external',
+                 type = 'real',
+                 value = 0.1234,
+                 texname = '\\text{keta}',
+                 lhablock = 'FRBlock',
+                 lhacode = [ 1 ])
+
+MZ = Parameter(name = 'MZ',
+               nature = 'external',
+               type = 'real',
+               value = 91.1876,
+               texname = '\\text{MZ}',
+               lhablock = 'MASS',
+               lhacode = [ 23 ])
+
+MTA = Parameter(name = 'MTA',
+                nature = 'external',
+                type = 'real',
+                value = 1.777,
+                texname = '\\text{MTA}',
+                lhablock = 'MASS',
+                lhacode = [ 15 ])
+
+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 ])
+
+MH = Parameter(name = 'MH',
+               nature = 'external',
+               type = 'real',
+               value = 125,
+               texname = '\\text{MH}',
+               lhablock = 'MASS',
+               lhacode = [ 25 ])
+
+Ms0 = Parameter(name = 'Ms0',
+                nature = 'external',
+                type = 'real',
+                value = 1567.0,
+                texname = '\\text{Ms0}',
+                lhablock = 'MASS',
+                lhacode = [ 9000006 ])
+
+WZ = Parameter(name = 'WZ',
+               nature = 'external',
+               type = 'real',
+               value = 2.4952,
+               texname = '\\text{WZ}',
+               lhablock = 'DECAY',
+               lhacode = [ 23 ])
+
+WW = Parameter(name = 'WW',
+               nature = 'external',
+               type = 'real',
+               value = 2.085,
+               texname = '\\text{WW}',
+               lhablock = 'DECAY',
+               lhacode = [ 24 ])
+
+WT = Parameter(name = 'WT',
+               nature = 'external',
+               type = 'real',
+               value = 1.50833649,
+               texname = '\\text{WT}',
+               lhablock = 'DECAY',
+               lhacode = [ 6 ])
+
+WH = Parameter(name = 'WH',
+               nature = 'external',
+               type = 'real',
+               value = 0.00407,
+               texname = '\\text{WH}',
+               lhablock = 'DECAY',
+               lhacode = [ 25 ])
+
+Ws0 = Parameter(name = 'Ws0',
+                nature = 'external',
+                type = 'real',
+                value = 1100.0,
+                texname = '\\text{Ws0}',
+                lhablock = 'DECAY',
+                lhacode = [ 9000006 ])
+
+aEW = Parameter(name = 'aEW',
+                nature = 'internal',
+                type = 'real',
+                value = '1/aEWM1',
+                texname = '\\alpha _{\\text{EW}}')
+
+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')
+
+vev = Parameter(name = 'vev',
+                nature = 'internal',
+                type = 'real',
+                value = '(2*MW*sw)/ee',
+                texname = '\\text{vev}')
+
+gweak = Parameter(name = 'gweak',
+                  nature = 'internal',
+                  type = 'real',
+                  value = 'gw',
+                  texname = '\\text{gweak}')
+
+lam = Parameter(name = 'lam',
+                nature = 'internal',
+                type = 'real',
+                value = 'MH**2/(2.*vev**2)',
+                texname = '\\text{lam}')
+
+vweak = Parameter(name = 'vweak',
+                  nature = 'internal',
+                  type = 'real',
+                  value = 'vev',
+                  texname = '\\text{vweak}')
+
+yb = Parameter(name = 'yb',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymb*cmath.sqrt(2))/vev',
+               texname = '\\text{yb}')
+
+yt = Parameter(name = 'yt',
+               nature = 'internal',
+               type = 'real',
+               value = '(ymt*cmath.sqrt(2))/vev',
+               texname = '\\text{yt}')
+
+ytau = Parameter(name = 'ytau',
+                 nature = 'internal',
+                 type = 'real',
+                 value = '(ymtau*cmath.sqrt(2))/vev',
+                 texname = '\\text{ytau}')
+
+muH = Parameter(name = 'muH',
+                nature = 'internal',
+                type = 'real',
+                value = 'cmath.sqrt(lam*vev**2)',
+                texname = '\\mu')
+
+I1b33 = Parameter(name = 'I1b33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yb',
+                  texname = '\\text{I1b33}')
+
+I2b33 = Parameter(name = 'I2b33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yt',
+                  texname = '\\text{I2b33}')
+
+I3b33 = Parameter(name = 'I3b33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yt',
+                  texname = '\\text{I3b33}')
+
+I4b33 = Parameter(name = 'I4b33',
+                  nature = 'internal',
+                  type = 'complex',
+                  value = 'yb',
+                  texname = '\\text{I4b33}')
+
diff --git a/CHresonance_neutral_scalar_UFO/particles.py b/CHresonance_neutral_scalar_UFO/particles.py
new file mode 100644
index 0000000000000000000000000000000000000000..5f253245b40cefbf1af2201b24f39adfcab59a33
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/particles.py
@@ -0,0 +1,401 @@
+# This file was automatically created by FeynRules 2.0.2
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (February 23, 2011)
+# Date: Wed 8 Oct 2014 20:39:09
+
+
+from __future__ import division
+from object_library import all_particles, Particle
+import parameters as Param
+
+import propagators as Prop
+
+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,
+             GhostNumber = 0,
+             LeptonNumber = 0,
+             Y = 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,
+             GhostNumber = 0,
+             LeptonNumber = 0,
+             Y = 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,
+                     GhostNumber = 0,
+                     LeptonNumber = 0,
+                     Y = 0)
+
+W__minus__ = W__plus__.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,
+             GhostNumber = 0,
+             LeptonNumber = 0,
+             Y = 0)
+
+ghA = Particle(pdg_code = 9000001,
+               name = 'ghA',
+               antiname = 'ghA~',
+               spin = -1,
+               color = 1,
+               mass = Param.ZERO,
+               width = Param.ZERO,
+               texname = 'ghA',
+               antitexname = 'ghA~',
+               charge = 0,
+               GhostNumber = 1,
+               LeptonNumber = 0,
+               Y = 0)
+
+ghA__tilde__ = ghA.anti()
+
+ghZ = Particle(pdg_code = 9000002,
+               name = 'ghZ',
+               antiname = 'ghZ~',
+               spin = -1,
+               color = 1,
+               mass = Param.MZ,
+               width = Param.WZ,
+               texname = 'ghZ',
+               antitexname = 'ghZ~',
+               charge = 0,
+               GhostNumber = 1,
+               LeptonNumber = 0,
+               Y = 0)
+
+ghZ__tilde__ = ghZ.anti()
+
+ghWp = Particle(pdg_code = 9000003,
+                name = 'ghWp',
+                antiname = 'ghWp~',
+                spin = -1,
+                color = 1,
+                mass = Param.MW,
+                width = Param.WW,
+                texname = 'ghWp',
+                antitexname = 'ghWp~',
+                charge = 1,
+                GhostNumber = 1,
+                LeptonNumber = 0,
+                Y = 0)
+
+ghWp__tilde__ = ghWp.anti()
+
+ghWm = Particle(pdg_code = 9000004,
+                name = 'ghWm',
+                antiname = 'ghWm~',
+                spin = -1,
+                color = 1,
+                mass = Param.MW,
+                width = Param.WW,
+                texname = 'ghWm',
+                antitexname = 'ghWm~',
+                charge = -1,
+                GhostNumber = 1,
+                LeptonNumber = 0,
+                Y = 0)
+
+ghWm__tilde__ = ghWm.anti()
+
+ghG = Particle(pdg_code = 9000005,
+               name = 'ghG',
+               antiname = 'ghG~',
+               spin = -1,
+               color = 8,
+               mass = Param.ZERO,
+               width = Param.ZERO,
+               texname = 'ghG',
+               antitexname = 'ghG~',
+               charge = 0,
+               GhostNumber = 1,
+               LeptonNumber = 0,
+               Y = 0)
+
+ghG__tilde__ = ghG.anti()
+
+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,
+              GhostNumber = 0,
+              LeptonNumber = 1,
+              Y = 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,
+              GhostNumber = 0,
+              LeptonNumber = 1,
+              Y = 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,
+              GhostNumber = 0,
+              LeptonNumber = 1,
+              Y = 0)
+
+vt__tilde__ = vt.anti()
+
+e__minus__ = Particle(pdg_code = 11,
+                      name = 'e-',
+                      antiname = 'e+',
+                      spin = 2,
+                      color = 1,
+                      mass = Param.ZERO,
+                      width = Param.ZERO,
+                      texname = 'e-',
+                      antitexname = 'e+',
+                      charge = -1,
+                      GhostNumber = 0,
+                      LeptonNumber = 1,
+                      Y = 0)
+
+e__plus__ = e__minus__.anti()
+
+mu__minus__ = Particle(pdg_code = 13,
+                       name = 'mu-',
+                       antiname = 'mu+',
+                       spin = 2,
+                       color = 1,
+                       mass = Param.ZERO,
+                       width = Param.ZERO,
+                       texname = 'mu-',
+                       antitexname = 'mu+',
+                       charge = -1,
+                       GhostNumber = 0,
+                       LeptonNumber = 1,
+                       Y = 0)
+
+mu__plus__ = mu__minus__.anti()
+
+ta__minus__ = Particle(pdg_code = 15,
+                       name = 'ta-',
+                       antiname = 'ta+',
+                       spin = 2,
+                       color = 1,
+                       mass = Param.MTA,
+                       width = Param.ZERO,
+                       texname = 'ta-',
+                       antitexname = 'ta+',
+                       charge = -1,
+                       GhostNumber = 0,
+                       LeptonNumber = 1,
+                       Y = 0)
+
+ta__plus__ = ta__minus__.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,
+             GhostNumber = 0,
+             LeptonNumber = 0,
+             Y = 0)
+
+u__tilde__ = u.anti()
+
+c = Particle(pdg_code = 4,
+             name = 'c',
+             antiname = 'c~',
+             spin = 2,
+             color = 3,
+             mass = Param.ZERO,
+             width = Param.ZERO,
+             texname = 'c',
+             antitexname = 'c~',
+             charge = 2/3,
+             GhostNumber = 0,
+             LeptonNumber = 0,
+             Y = 0)
+
+c__tilde__ = c.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,
+             GhostNumber = 0,
+             LeptonNumber = 0,
+             Y = 0)
+
+t__tilde__ = t.anti()
+
+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,
+             GhostNumber = 0,
+             LeptonNumber = 0,
+             Y = 0)
+
+d__tilde__ = d.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,
+             GhostNumber = 0,
+             LeptonNumber = 0,
+             Y = 0)
+
+s__tilde__ = s.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,
+             GhostNumber = 0,
+             LeptonNumber = 0,
+             Y = 0)
+
+b__tilde__ = b.anti()
+
+H = Particle(pdg_code = 25,
+             name = 'H',
+             antiname = 'H',
+             spin = 1,
+             color = 1,
+             mass = Param.MH,
+             width = Param.WH,
+             texname = 'H',
+             antitexname = 'H',
+             charge = 0,
+             GhostNumber = 0,
+             LeptonNumber = 0,
+             Y = 0)
+
+G0 = Particle(pdg_code = 250,
+              name = 'G0',
+              antiname = 'G0',
+              spin = 1,
+              color = 1,
+              mass = Param.MZ,
+              width = Param.WZ,
+              texname = 'G0',
+              antitexname = 'G0',
+              goldstone = True,
+              charge = 0,
+              GhostNumber = 0,
+              LeptonNumber = 0,
+              Y = 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-',
+                     goldstone = True,
+                     charge = 1,
+                     GhostNumber = 0,
+                     LeptonNumber = 0,
+                     Y = 0)
+
+G__minus__ = G__plus__.anti()
+
+s0 = Particle(pdg_code = 9000006,
+              name = 's0',
+              antiname = 's0',
+              spin = 1,
+              color = 1,
+              mass = Param.Ms0,
+              width = Param.Ws0,
+              texname = 's0',
+              antitexname = 's0',
+              charge = 0,
+              GhostNumber = 0,
+              LeptonNumber = 0,
+              Y = 0)
+
diff --git a/CHresonance_neutral_scalar_UFO/propagators.py b/CHresonance_neutral_scalar_UFO/propagators.py
new file mode 100644
index 0000000000000000000000000000000000000000..8d421595a33106e2a7ff44435b5d55a556e2a062
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/propagators.py
@@ -0,0 +1,35 @@
+# This file was automatically created by FeynRules 1.7.69
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (November 6, 2010)
+# Date: Mon 1 Oct 2012 14:58:26
+
+from object_library import all_propagators, Propagator
+
+
+# define only once the denominator since this is always the same
+denominator = "P('mu', id) * P('mu', id) - Mass(id) * Mass(id) + complex(0,1) * Mass(id) * Width(id)"
+
+# propagator for the scalar
+S = Propagator(name = "S",
+               numerator = "complex(0,1)",
+               denominator = denominator
+               )
+
+# propagator for the incoming fermion # the one for the outcomming is computed on the flight
+F = Propagator(name = "F",
+                numerator = "complex(0,1) * (Gamma('mu', s1, s2) * P('mu', id) + Mass(id) * Identity(s1, s2))",
+                denominator = denominator
+               )
+
+# massive vector in the unitary gauge, can't be use for massless particles
+V1 = Propagator(name = "V1",
+                numerator = "complex(0,1) * (-1 * Metric(l1, l2) + Metric(l1,'mu')* P('mu', id) * P(l2, id) / Mass(id)**2 ",
+                denominator = denominator
+               )
+
+# massless vector and massive vector in unitary gauge
+V2 = Propagator(name = "V2",
+                numerator = "complex(0,-1) * Metric(l1, l2)",
+                denominator =  "P('mu', id) * P('mu', id)"
+               )
+
+
diff --git a/CHresonance_neutral_scalar_UFO/vertices.py b/CHresonance_neutral_scalar_UFO/vertices.py
new file mode 100644
index 0000000000000000000000000000000000000000..232cd4b090af587c28b3f17c69e903d236bf8e6b
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/vertices.py
@@ -0,0 +1,743 @@
+# This file was automatically created by FeynRules 2.0.2
+# Mathematica version: 8.0 for Mac OS X x86 (64-bit) (February 23, 2011)
+# Date: Wed 8 Oct 2014 20:39:09
+
+
+from object_library import all_vertices, Vertex
+import particles as P
+import couplings as C
+import lorentz as L
+
+
+V_1 = Vertex(name = 'V_1',
+             particles = [ P.G0, P.G0, P.G0, P.G0 ],
+             color = [ '1' ],
+             lorentz = [ L.SSSS1 ],
+             couplings = {(0,0):C.GC_19})
+
+V_2 = Vertex(name = 'V_2',
+             particles = [ P.G0, P.G0, P.G__minus__, P.G__plus__ ],
+             color = [ '1' ],
+             lorentz = [ L.SSSS1 ],
+             couplings = {(0,0):C.GC_17})
+
+V_3 = Vertex(name = 'V_3',
+             particles = [ P.G__minus__, P.G__minus__, P.G__plus__, P.G__plus__ ],
+             color = [ '1' ],
+             lorentz = [ L.SSSS1 ],
+             couplings = {(0,0):C.GC_18})
+
+V_4 = Vertex(name = 'V_4',
+             particles = [ P.G0, P.G0, P.H, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSSS1 ],
+             couplings = {(0,0):C.GC_17})
+
+V_5 = Vertex(name = 'V_5',
+             particles = [ P.G__minus__, P.G__plus__, P.H, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSSS1 ],
+             couplings = {(0,0):C.GC_17})
+
+V_6 = Vertex(name = 'V_6',
+             particles = [ P.H, P.H, P.H, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSSS1 ],
+             couplings = {(0,0):C.GC_19})
+
+V_7 = Vertex(name = 'V_7',
+             particles = [ P.G0, P.G0, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSS1 ],
+             couplings = {(0,0):C.GC_47})
+
+V_8 = Vertex(name = 'V_8',
+             particles = [ P.G__minus__, P.G__plus__, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSS1 ],
+             couplings = {(0,0):C.GC_47})
+
+V_9 = Vertex(name = 'V_9',
+             particles = [ P.H, P.H, P.H ],
+             color = [ '1' ],
+             lorentz = [ L.SSS1 ],
+             couplings = {(0,0):C.GC_48})
+
+V_10 = Vertex(name = 'V_10',
+              particles = [ P.a, P.a, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_6})
+
+V_11 = Vertex(name = 'V_11',
+              particles = [ P.a, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_12 = Vertex(name = 'V_12',
+              particles = [ P.ghA, P.ghWm__tilde__, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_13 = Vertex(name = 'V_13',
+              particles = [ P.ghA, P.ghWp__tilde__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_4})
+
+V_14 = Vertex(name = 'V_14',
+              particles = [ P.ghWm, P.ghA__tilde__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUS1 ],
+              couplings = {(0,0):C.GC_54})
+
+V_15 = Vertex(name = 'V_15',
+              particles = [ P.ghWm, P.ghA__tilde__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_16 = Vertex(name = 'V_16',
+              particles = [ P.ghWm, P.ghWm__tilde__, P.G0 ],
+              color = [ '1' ],
+              lorentz = [ L.UUS1 ],
+              couplings = {(0,0):C.GC_49})
+
+V_17 = Vertex(name = 'V_17',
+              particles = [ P.ghWm, P.ghWm__tilde__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.UUS1 ],
+              couplings = {(0,0):C.GC_50})
+
+V_18 = Vertex(name = 'V_18',
+              particles = [ P.ghWm, P.ghWm__tilde__, P.a ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_4})
+
+V_19 = Vertex(name = 'V_19',
+              particles = [ P.ghWm, P.ghWm__tilde__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_32})
+
+V_20 = Vertex(name = 'V_20',
+              particles = [ P.ghWm, P.ghZ__tilde__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUS1 ],
+              couplings = {(0,0):C.GC_57})
+
+V_21 = Vertex(name = 'V_21',
+              particles = [ P.ghWm, P.ghZ__tilde__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_31})
+
+V_22 = Vertex(name = 'V_22',
+              particles = [ P.ghWp, P.ghA__tilde__, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUS1 ],
+              couplings = {(0,0):C.GC_53})
+
+V_23 = Vertex(name = 'V_23',
+              particles = [ P.ghWp, P.ghA__tilde__, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_4})
+
+V_24 = Vertex(name = 'V_24',
+              particles = [ P.ghWp, P.ghWp__tilde__, P.G0 ],
+              color = [ '1' ],
+              lorentz = [ L.UUS1 ],
+              couplings = {(0,0):C.GC_52})
+
+V_25 = Vertex(name = 'V_25',
+              particles = [ P.ghWp, P.ghWp__tilde__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.UUS1 ],
+              couplings = {(0,0):C.GC_50})
+
+V_26 = Vertex(name = 'V_26',
+              particles = [ P.ghWp, P.ghWp__tilde__, P.a ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_27 = Vertex(name = 'V_27',
+              particles = [ P.ghWp, P.ghWp__tilde__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_31})
+
+V_28 = Vertex(name = 'V_28',
+              particles = [ P.ghWp, P.ghZ__tilde__, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUS1 ],
+              couplings = {(0,0):C.GC_56})
+
+V_29 = Vertex(name = 'V_29',
+              particles = [ P.ghWp, P.ghZ__tilde__, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_32})
+
+V_30 = Vertex(name = 'V_30',
+              particles = [ P.ghZ, P.ghWm__tilde__, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUS1 ],
+              couplings = {(0,0):C.GC_58})
+
+V_31 = Vertex(name = 'V_31',
+              particles = [ P.ghZ, P.ghWm__tilde__, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_31})
+
+V_32 = Vertex(name = 'V_32',
+              particles = [ P.ghZ, P.ghWp__tilde__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUS1 ],
+              couplings = {(0,0):C.GC_55})
+
+V_33 = Vertex(name = 'V_33',
+              particles = [ P.ghZ, P.ghWp__tilde__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_32})
+
+V_34 = Vertex(name = 'V_34',
+              particles = [ P.ghZ, P.ghZ__tilde__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.UUS1 ],
+              couplings = {(0,0):C.GC_59})
+
+V_35 = Vertex(name = 'V_35',
+              particles = [ P.ghG, P.ghG__tilde__, P.g ],
+              color = [ 'f(1,2,3)' ],
+              lorentz = [ L.UUV1 ],
+              couplings = {(0,0):C.GC_10})
+
+V_36 = Vertex(name = 'V_36',
+              particles = [ P.g, P.g, P.g ],
+              color = [ 'f(1,2,3)' ],
+              lorentz = [ L.VVV1 ],
+              couplings = {(0,0):C.GC_10})
+
+V_37 = Vertex(name = 'V_37',
+              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_12,(0,0):C.GC_12,(2,2):C.GC_12})
+
+V_38 = Vertex(name = 'V_38',
+              particles = [ P.t__tilde__, P.b, P.G__plus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1, L.FFS3 ],
+              couplings = {(0,0):C.GC_13,(0,1):C.GC_14})
+
+V_39 = Vertex(name = 'V_39',
+              particles = [ P.b__tilde__, P.b, P.G0 ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS2 ],
+              couplings = {(0,0):C.GC_62})
+
+V_40 = Vertex(name = 'V_40',
+              particles = [ P.b__tilde__, P.b, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS4 ],
+              couplings = {(0,0):C.GC_63})
+
+V_41 = Vertex(name = 'V_41',
+              particles = [ P.vt__tilde__, P.ta__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFS1 ],
+              couplings = {(0,0):C.GC_67})
+
+V_42 = Vertex(name = 'V_42',
+              particles = [ P.ta__plus__, P.ta__minus__, P.G0 ],
+              color = [ '1' ],
+              lorentz = [ L.FFS2 ],
+              couplings = {(0,0):C.GC_68})
+
+V_43 = Vertex(name = 'V_43',
+              particles = [ P.ta__plus__, P.ta__minus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.FFS4 ],
+              couplings = {(0,0):C.GC_69})
+
+V_44 = Vertex(name = 'V_44',
+              particles = [ P.t__tilde__, P.t, P.G0 ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS2 ],
+              couplings = {(0,0):C.GC_65})
+
+V_45 = Vertex(name = 'V_45',
+              particles = [ P.t__tilde__, P.t, P.H ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS4 ],
+              couplings = {(0,0):C.GC_64})
+
+V_46 = Vertex(name = 'V_46',
+              particles = [ P.b__tilde__, P.t, P.G__minus__ ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFS1, L.FFS3 ],
+              couplings = {(0,0):C.GC_15,(0,1):C.GC_16})
+
+V_47 = Vertex(name = 'V_47',
+              particles = [ P.a, P.W__minus__, P.G0, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_34})
+
+V_48 = Vertex(name = 'V_48',
+              particles = [ P.a, P.W__minus__, P.G__plus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_33})
+
+V_49 = Vertex(name = 'V_49',
+              particles = [ P.a, P.W__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_53})
+
+V_50 = Vertex(name = 'V_50',
+              particles = [ P.W__minus__, P.G0, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_27})
+
+V_51 = Vertex(name = 'V_51',
+              particles = [ P.W__minus__, P.G__plus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_25})
+
+V_52 = Vertex(name = 'V_52',
+              particles = [ P.a, P.W__minus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVV1 ],
+              couplings = {(0,0):C.GC_4})
+
+V_53 = Vertex(name = 'V_53',
+              particles = [ P.a, P.W__plus__, P.G0, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_34})
+
+V_54 = Vertex(name = 'V_54',
+              particles = [ P.a, P.W__plus__, P.G__minus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_35})
+
+V_55 = Vertex(name = 'V_55',
+              particles = [ P.a, P.W__plus__, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_54})
+
+V_56 = Vertex(name = 'V_56',
+              particles = [ P.W__plus__, P.G0, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_26})
+
+V_57 = Vertex(name = 'V_57',
+              particles = [ P.W__plus__, P.G__minus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_25})
+
+V_58 = Vertex(name = 'V_58',
+              particles = [ P.W__minus__, P.W__plus__, P.G0, P.G0 ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_22})
+
+V_59 = Vertex(name = 'V_59',
+              particles = [ P.W__minus__, P.W__plus__, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_22})
+
+V_60 = Vertex(name = 'V_60',
+              particles = [ P.W__minus__, P.W__plus__, P.H, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_22})
+
+V_61 = Vertex(name = 'V_61',
+              particles = [ P.W__minus__, P.W__plus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_51})
+
+V_62 = Vertex(name = 'V_62',
+              particles = [ P.a, P.a, P.W__minus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVVV2 ],
+              couplings = {(0,0):C.GC_5})
+
+V_63 = Vertex(name = 'V_63',
+              particles = [ P.W__minus__, P.W__plus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.VVV1 ],
+              couplings = {(0,0):C.GC_32})
+
+V_64 = Vertex(name = 'V_64',
+              particles = [ P.W__minus__, P.W__minus__, P.W__plus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVVV2 ],
+              couplings = {(0,0):C.GC_23})
+
+V_65 = Vertex(name = 'V_65',
+              particles = [ P.ta__plus__, P.vt, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFS3 ],
+              couplings = {(0,0):C.GC_66})
+
+V_66 = Vertex(name = 'V_66',
+              particles = [ P.a, P.Z, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_42})
+
+V_67 = Vertex(name = 'V_67',
+              particles = [ P.Z, P.G0, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_39})
+
+V_68 = Vertex(name = 'V_68',
+              particles = [ P.Z, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VSS1 ],
+              couplings = {(0,0):C.GC_40})
+
+V_69 = Vertex(name = 'V_69',
+              particles = [ P.W__minus__, P.Z, P.G0, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_8})
+
+V_70 = Vertex(name = 'V_70',
+              particles = [ P.W__minus__, P.Z, P.G__plus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_9})
+
+V_71 = Vertex(name = 'V_71',
+              particles = [ P.W__minus__, P.Z, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_46})
+
+V_72 = Vertex(name = 'V_72',
+              particles = [ P.W__plus__, P.Z, P.G0, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_8})
+
+V_73 = Vertex(name = 'V_73',
+              particles = [ P.W__plus__, P.Z, P.G__minus__, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_7})
+
+V_74 = Vertex(name = 'V_74',
+              particles = [ P.W__plus__, P.Z, P.G__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_45})
+
+V_75 = Vertex(name = 'V_75',
+              particles = [ P.a, P.W__minus__, P.W__plus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.VVVV5 ],
+              couplings = {(0,0):C.GC_36})
+
+V_76 = Vertex(name = 'V_76',
+              particles = [ P.Z, P.Z, P.G0, P.G0 ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_44})
+
+V_77 = Vertex(name = 'V_77',
+              particles = [ P.Z, P.Z, P.G__minus__, P.G__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_43})
+
+V_78 = Vertex(name = 'V_78',
+              particles = [ P.Z, P.Z, P.H, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVSS1 ],
+              couplings = {(0,0):C.GC_44})
+
+V_79 = Vertex(name = 'V_79',
+              particles = [ P.Z, P.Z, P.H ],
+              color = [ '1' ],
+              lorentz = [ L.VVS1 ],
+              couplings = {(0,0):C.GC_60})
+
+V_80 = Vertex(name = 'V_80',
+              particles = [ P.W__minus__, P.W__plus__, P.Z, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.VVVV2 ],
+              couplings = {(0,0):C.GC_24})
+
+V_81 = Vertex(name = 'V_81',
+              particles = [ P.e__plus__, P.ve, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_28})
+
+V_82 = Vertex(name = 'V_82',
+              particles = [ P.mu__plus__, P.vm, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_28})
+
+V_83 = Vertex(name = 'V_83',
+              particles = [ P.ta__plus__, P.vt, P.W__minus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_28})
+
+V_84 = Vertex(name = 'V_84',
+              particles = [ P.ve__tilde__, P.ve, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_41})
+
+V_85 = Vertex(name = 'V_85',
+              particles = [ P.vm__tilde__, P.vm, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_41})
+
+V_86 = Vertex(name = 'V_86',
+              particles = [ P.vt__tilde__, P.vt, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_41})
+
+V_87 = Vertex(name = 'V_87',
+              particles = [ P.ve__tilde__, P.e__minus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_28})
+
+V_88 = Vertex(name = 'V_88',
+              particles = [ P.vm__tilde__, P.mu__minus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_28})
+
+V_89 = Vertex(name = 'V_89',
+              particles = [ P.vt__tilde__, P.ta__minus__, P.W__plus__ ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2 ],
+              couplings = {(0,0):C.GC_28})
+
+V_90 = Vertex(name = 'V_90',
+              particles = [ P.e__plus__, P.e__minus__, P.a ],
+              color = [ '1' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_91 = Vertex(name = 'V_91',
+              particles = [ P.mu__plus__, P.mu__minus__, P.a ],
+              color = [ '1' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_92 = Vertex(name = 'V_92',
+              particles = [ P.ta__plus__, P.ta__minus__, P.a ],
+              color = [ '1' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_3})
+
+V_93 = Vertex(name = 'V_93',
+              particles = [ P.e__plus__, P.e__minus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2, L.FFV4 ],
+              couplings = {(0,0):C.GC_29,(0,1):C.GC_38})
+
+V_94 = Vertex(name = 'V_94',
+              particles = [ P.mu__plus__, P.mu__minus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2, L.FFV4 ],
+              couplings = {(0,0):C.GC_29,(0,1):C.GC_38})
+
+V_95 = Vertex(name = 'V_95',
+              particles = [ P.ta__plus__, P.ta__minus__, P.Z ],
+              color = [ '1' ],
+              lorentz = [ L.FFV2, L.FFV4 ],
+              couplings = {(0,0):C.GC_29,(0,1):C.GC_38})
+
+V_96 = Vertex(name = 'V_96',
+              particles = [ P.u__tilde__, P.u, P.a ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_2})
+
+V_97 = Vertex(name = 'V_97',
+              particles = [ P.c__tilde__, P.c, P.a ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_2})
+
+V_98 = Vertex(name = 'V_98',
+              particles = [ P.t__tilde__, P.t, P.a ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV1 ],
+              couplings = {(0,0):C.GC_2})
+
+V_99 = Vertex(name = 'V_99',
+              particles = [ P.u__tilde__, P.u, P.Z ],
+              color = [ 'Identity(1,2)' ],
+              lorentz = [ L.FFV2, L.FFV5 ],
+              couplings = {(0,0):C.GC_30,(0,1):C.GC_37})
+
+V_100 = Vertex(name = 'V_100',
+               particles = [ P.c__tilde__, P.c, P.Z ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV2, L.FFV5 ],
+               couplings = {(0,0):C.GC_30,(0,1):C.GC_37})
+
+V_101 = Vertex(name = 'V_101',
+               particles = [ P.t__tilde__, P.t, P.Z ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV2, L.FFV5 ],
+               couplings = {(0,0):C.GC_30,(0,1):C.GC_37})
+
+V_102 = Vertex(name = 'V_102',
+               particles = [ P.u__tilde__, P.u, P.g ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_11})
+
+V_103 = Vertex(name = 'V_103',
+               particles = [ P.c__tilde__, P.c, P.g ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_11})
+
+V_104 = Vertex(name = 'V_104',
+               particles = [ P.t__tilde__, P.t, P.g ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_11})
+
+V_105 = Vertex(name = 'V_105',
+               particles = [ P.d__tilde__, P.u, P.W__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV2 ],
+               couplings = {(0,0):C.GC_28})
+
+V_106 = Vertex(name = 'V_106',
+               particles = [ P.s__tilde__, P.c, P.W__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV2 ],
+               couplings = {(0,0):C.GC_28})
+
+V_107 = Vertex(name = 'V_107',
+               particles = [ P.b__tilde__, P.t, P.W__minus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV2 ],
+               couplings = {(0,0):C.GC_28})
+
+V_108 = Vertex(name = 'V_108',
+               particles = [ P.u__tilde__, P.d, P.W__plus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV2 ],
+               couplings = {(0,0):C.GC_28})
+
+V_109 = Vertex(name = 'V_109',
+               particles = [ P.c__tilde__, P.s, P.W__plus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV2 ],
+               couplings = {(0,0):C.GC_28})
+
+V_110 = Vertex(name = 'V_110',
+               particles = [ P.t__tilde__, P.b, P.W__plus__ ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV2 ],
+               couplings = {(0,0):C.GC_28})
+
+V_111 = Vertex(name = 'V_111',
+               particles = [ P.d__tilde__, P.d, P.a ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_1})
+
+V_112 = Vertex(name = 'V_112',
+               particles = [ P.s__tilde__, P.s, P.a ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_1})
+
+V_113 = Vertex(name = 'V_113',
+               particles = [ P.b__tilde__, P.b, P.a ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_1})
+
+V_114 = Vertex(name = 'V_114',
+               particles = [ P.d__tilde__, P.d, P.Z ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV2, L.FFV3 ],
+               couplings = {(0,0):C.GC_29,(0,1):C.GC_37})
+
+V_115 = Vertex(name = 'V_115',
+               particles = [ P.s__tilde__, P.s, P.Z ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV2, L.FFV3 ],
+               couplings = {(0,0):C.GC_29,(0,1):C.GC_37})
+
+V_116 = Vertex(name = 'V_116',
+               particles = [ P.b__tilde__, P.b, P.Z ],
+               color = [ 'Identity(1,2)' ],
+               lorentz = [ L.FFV2, L.FFV3 ],
+               couplings = {(0,0):C.GC_29,(0,1):C.GC_37})
+
+V_117 = Vertex(name = 'V_117',
+               particles = [ P.d__tilde__, P.d, P.g ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_11})
+
+V_118 = Vertex(name = 'V_118',
+               particles = [ P.s__tilde__, P.s, P.g ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_11})
+
+V_119 = Vertex(name = 'V_119',
+               particles = [ P.b__tilde__, P.b, P.g ],
+               color = [ 'T(3,2,1)' ],
+               lorentz = [ L.FFV1 ],
+               couplings = {(0,0):C.GC_11})
+
+V_120 = Vertex(name = 'V_120',
+               particles = [ P.H, P.H, P.s0 ],
+               color = [ '1' ],
+               lorentz = [ L.SSS2 ],
+               couplings = {(0,0):C.GC_61})
+
+V_121 = Vertex(name = 'V_121',
+               particles = [ P.W__minus__, P.W__plus__, P.s0 ],
+               color = [ '1' ],
+               lorentz = [ L.VVS1 ],
+               couplings = {(0,0):C.GC_20})
+
+V_122 = Vertex(name = 'V_122',
+               particles = [ P.Z, P.Z, P.s0 ],
+               color = [ '1' ],
+               lorentz = [ L.VVS1 ],
+               couplings = {(0,0):C.GC_21})
+
diff --git a/CHresonance_neutral_scalar_UFO/write_param_card.py b/CHresonance_neutral_scalar_UFO/write_param_card.py
new file mode 100644
index 0000000000000000000000000000000000000000..57a85b061440f522b264b837e691f5b33accc96b
--- /dev/null
+++ b/CHresonance_neutral_scalar_UFO/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 f4f397c550fc3d707d880114162d07bec8ed484c..da24eeafa50311d2c1dcef491dda5e58ca2115c2 100644
--- a/model_list.txt
+++ b/model_list.txt
@@ -160,6 +160,12 @@ Contents: Axion-like particle, coupling to up-type quarks. ALP lifetime set by c
 Provided by Adrian Carmona [adrian@ugr.es], author of the ALPs from the Top paper
 Paper ALPs from the Top: https://arxiv.org/abs/2202.09371 
 
+CHresonance_neutral_scalar_UFO
+Requestor: Zhen Wang
+Content: Singlet-catalyzed electroweak phase transitions from Ashutosh V. Kotwal
+Paper: https://journals.aps.org/prd/abstract/10.1103/PhysRevD.94.035022
+JIRA: https://its.cern.ch/jira/browse/ATLMCPROD-6413 (old ZZ sample request)
+
 clfvLO
 Requestor: Carlo Gottardo
 Content: Charged lepton flavour violation, with effective 4-fermion operators