JHUGen MELA  v2.4.1
Matrix element calculations as used in JHUGen. MELA is an important tool that was used for the Higgs boson discovery and for precise measurements of its structure and interactions. Please see the website https://spin.pha.jhu.edu/ and papers cited there for more details, and kindly cite those papers when using this code.
testME_more.py
Go to the documentation of this file.
1 #!/usr/bin/env python
2 
3 """
4 test python and some MEs
5 """
6 
7 if __name__ == "__main__":
8  import argparse
9  parser = argparse.ArgumentParser()
10  parser.add_argument("seed", type=int)
11  parser.add_argument('unittest_args', nargs='*')
12  args = parser.parse_args()
13  seed = args.seed
14 
15 import os
16 import random
17 import re
18 import sys
19 import unittest
20 
21 from mela import Mela, TVar
22 
23 event1_ggH = """
24 <event>
25  9 0 1.0000000E+00 6.2499553E+01 7.8125000E-03 1.3828638E-01
26  21 -1 0 0 501 502 0.00000000000E+00 0.00000000000E+00 5.37086667033E+02 5.37086667033E+02 0.00000000000E+00 0.00000000000E+00 1.
27  21 -1 0 0 502 501 0.00000000000E+00 0.00000000000E+00 -7.27293062837E+00 7.27293062837E+00 0.00000000000E+00 0.00000000000E+00 1.
28  25 2 1 2 0 0 2.66453525910E-15 0.00000000000E+00 5.29813736405E+02 5.44359597662E+02 1.24999105129E+02 0.00000000000E+00 1.
29  23 2 3 3 0 0 -1.07287063341E+00 2.03413884495E+01 2.12763358296E+02 2.18095469074E+02 4.33873698414E+01 0.00000000000E+00 1.
30  23 2 3 3 0 0 1.07287063341E+00 -2.03413884495E+01 3.17050378109E+02 3.26264128588E+02 7.42456477420E+01 0.00000000000E+00 1.
31  11 1 4 4 0 0 -4.96385765951E+00 1.62933935971E+01 2.11696003997E+02 2.12380113632E+02 5.10998597706E-04 0.00000000000E+00 1.
32  -11 1 4 4 0 0 3.89098702611E+00 4.04799485238E+00 1.06735429842E+00 5.71535544142E+00 5.10999912595E-04 0.00000000000E+00 1.
33  12 1 5 5 0 0 8.04215953543E+00 -2.29213235677E+01 2.77968427622E+01 3.69151442044E+01 4.63538696652E-06 0.00000000000E+00 1.
34  -12 1 5 5 0 0 -6.96928890202E+00 2.57993511828E+00 2.89253535347E+02 2.89348984384E+02 5.39479660939E-06 0.00000000000E+00 1.
35 </event>
36 """
37 
38 event2_VBF = """
39 <event>
40 11 60 1.0000000E+00 1.2500000E+02 7.8125000E-03 1.2380607E-01
41  -4 -1 0 0 501 0 0.00000000000E+00 0.00000000000E+00 3.03951895455E+01 3.03951895455E+01 0.00000000000E+00 0.00000000000E+00 1.
42  1 -1 0 0 502 0 0.00000000000E+00 0.00000000000E+00 -7.16454865818E+02 7.16454865818E+02 0.00000000000E+00 0.00000000000E+00 1.
43  -4 1 1 2 501 0 -6.29623305954E+01 -4.67693451618E+01 -5.85204921720E+01 9.78584422772E+01 0.00000000000E+00 0.00000000000E+00 1.
44  1 1 1 2 502 0 9.29867095677E+00 2.61406963578E+01 -9.53450713739E+01 9.92999694004E+01 0.00000000000E+00 0.00000000000E+00 1.
45  25 2 1 2 0 0 5.36636596386E+01 2.06286488040E+01 -5.32194112726E+02 5.49691643686E+02 1.25000000000E+02 0.00000000000E+00 1.
46  23 2 5 5 0 0 -4.38443764422E-01 -1.86323403603E+01 -9.56577645198E+01 9.86840261342E+01 1.55168540914E+01 0.00000000000E+00 1.
47  23 2 5 5 0 0 5.41021034030E+01 3.92609891643E+01 -4.36536348206E+02 4.51007617552E+02 9.15173476546E+01 0.00000000000E+00 1.
48  -11 1 6 6 0 0 3.31752759767E+00 1.23639112123E+00 -8.43121325514E+00 9.14439771558E+00 5.10999986247E-04 0.00000000000E+00 1.
49  11 1 6 6 0 0 -3.75597136210E+00 -1.98687314816E+01 -8.72265512647E+01 8.95396284186E+01 5.11001208360E-04 0.00000000000E+00 1.
50  -13 1 7 7 0 0 6.93884009242E+01 3.00845351746E+01 -2.00352143436E+02 2.14151399675E+02 1.05660000076E-01 0.00000000000E+00 1.
51  13 1 7 7 0 0 -1.52862975212E+01 9.17645398978E+00 -2.36184204771E+02 2.36856217877E+02 1.05659999949E-01 0.00000000000E+00 1.
52 </event>
53 """
54 
55 event3_ZH = """
56 <event>
57 12 50 1.0000000E+00 1.2500000E+02 7.8125000E-03 1.2380607E-01
58  1 -1 0 0 503 0 0.00000000000E+00 0.00000000000E+00 5.98124082110E+01 5.98124082110E+01 0.00000000000E+00 0.00000000000E+00 1.
59  -1 -1 0 0 0 503 0.00000000000E+00 0.00000000000E+00 -3.44771438947E+02 3.44771438947E+02 0.00000000000E+00 0.00000000000E+00 1.
60  23 2 1 2 0 0 -1.33768582227E+01 8.88541882213E+01 -1.21353594052E+02 1.79939405232E+02 9.78646395061E+01 0.00000000000E+00 1.
61  25 2 1 2 0 0 1.33768582227E+01 -8.88541882213E+01 -1.63605436683E+02 2.24644441925E+02 1.24997517076E+02 0.00000000000E+00 1.
62  13 1 3 3 0 0 3.17738029751E+01 -1.54337455361E+00 -2.87056115717E+01 4.28483355599E+01 1.05660000003E-01 0.00000000000E+00 1.
63  -13 1 3 3 0 0 -4.51506611978E+01 9.03975627749E+01 -9.26479824805E+01 1.37091069673E+02 1.05660000012E-01 0.00000000000E+00 1.
64  23 2 4 4 0 0 4.40015216426E+01 -3.16198925284E+01 -5.72188020543E+01 7.95131988002E+01 1.06021679173E+01 0.00000000000E+00 1.
65  23 2 4 4 0 0 -3.06246634199E+01 -5.72342956929E+01 -1.06386634629E+02 1.45131243125E+02 7.43728921737E+01 0.00000000000E+00 1.
66  -11 1 8 8 0 0 1.78742924427E+01 -1.12102312034E+01 -2.03424438331E+00 2.11966451221E+01 5.10999863718E-04 0.00000000000E+00 1.
67  11 1 7 7 0 0 3.86248305415E+01 -3.12040320381E+01 -4.96421050079E+01 7.02133017289E+01 5.11000665199E-04 0.00000000000E+00 1.
68  -11 1 7 7 0 0 5.37669110106E+00 -4.15860490296E-01 -7.57669704638E+00 9.29989707133E+00 5.11000011708E-04 0.00000000000E+00 1.
69  11 1 8 8 0 0 -4.84989558626E+01 -4.60240644895E+01 -1.04352390245E+02 1.23934598003E+02 5.10996863057E-04 0.00000000000E+00 1.
70 </event>
71 """
72 
73 class TestMela(unittest.TestCase):
74  @classmethod
75  def setUpClass(cls):
76  cls.m = Mela()
77  def setUp(self):
78  random.seed(seed)
79  def tearDown(self):
80  self.m.resetInputEvent()
81  def testcontact_decay(self):
82  self.runcontact(event1_ggH, False, TVar.SelfDefine_spin0, TVar.JHUGen, TVar.ZZINDEPENDENT)
83  def testcontact_VBF(self):
84  self.runcontact(event2_VBF, True, TVar.SelfDefine_spin0, TVar.JHUGen, TVar.JJVBF)
85  def testcontact_ZH(self):
86  self.runcontact(event3_ZH, True, TVar.SelfDefine_spin0, TVar.JHUGen, TVar.Lep_ZH)
87  def runcontact(self, event, isprod, *setprocessargs):
88  """
89  correspondence between a1 L1 L1Zg and contact terms
90  """
91 
92  M_Z = 91.1876
93  Ga_Z = 2.4952
94  aL_lep = -0.53762
95  aR_lep = 0.46238
96  aL_nu = 1
97  aR_nu = 0
98 
99  sitW2 = 0.23119
100  aL_up = (-2 * sitW2 * 2./3 + 1)
101  aR_up = -2 * sitW2 * 2./3
102  aL_dn = (-2 * sitW2 * -1./3 - 1)
103  aR_dn = -2 * sitW2 * -1./3
104 
105  e = 0.8431872482432357 # = cL_lep = cR_lep from mod_Parameters
106  c_up = -2./3 * e
107  c_dn = 1./3 * e
108  L1 = 10000.
109 
110  m = self.m
111  computeP = m.computeProdP if isprod else m.computeP
112 
113  ghz1, ghz1_prime2, ghzgs1_prime2 = random.uniform(-1, 1), random.uniform(-10000, 10000), random.uniform(-10000, 10000)
114  m.setInputEvent_fromLHE_Hwithdecay(event, True)
115 
116  m.setProcess(*setprocessargs)
117  m.ghz1 = ghz1
118  m.ghz1_prime2 = ghz1_prime2
119  m.ghzgs1_prime2 = ghzgs1_prime2
120  me1 = computeP()
121 
122  m.setProcess(*setprocessargs)
123  m.ghz1 = ghz1 + 2 * ghz1_prime2 * (M_Z**2 - 1j*M_Z*Ga_Z)/L1**2
124  m.ghzzp1 = M_Z**2/L1**2
125  m.ezp_El_left = m.ezp_Mu_left = m.ezp_Ta_left = aL_lep * ghz1_prime2 + e * ghzgs1_prime2
126  m.ezp_El_right = m.ezp_Mu_right = m.ezp_Ta_right = aR_lep * ghz1_prime2 + e * ghzgs1_prime2
127  m.ezp_Up_left = m.ezp_Chm_left = m.ezp_Top_left = aL_up * ghz1_prime2 + c_up * ghzgs1_prime2
128  m.ezp_Up_right = m.ezp_Chm_right = m.ezp_Top_right = aR_up * ghz1_prime2 + c_up * ghzgs1_prime2
129  m.ezp_Dn_left = m.ezp_Str_left = m.ezp_Bot_left = aL_dn * ghz1_prime2 + c_dn * ghzgs1_prime2
130  m.ezp_Dn_right = m.ezp_Str_right = m.ezp_Bot_right = aR_dn * ghz1_prime2 + c_dn * ghzgs1_prime2
131  m.ezp_NuE_left = aL_nu * ghz1_prime2
132  m.ezp_NuE_right = aR_nu * ghz1_prime2
133  me2 = computeP()
134 
135  self.assertEquals(me1, me2)
136  self.assertNotEquals(me1, 0)
137 
138  def Vprimekwargs(self, process, usevpvp):
139  def couplings():
140  if process in (0, 50, 60):
141  for i in range(1, 5):
142  for j in range(8):
143  if j in (3, 6) and not usevpvp: continue #for ZZ', q1^2-q2^2 is defined differently
144  coupling = "ghz{}".format(i)
145  if j: coupling += "_prime"
146  if j>1: coupling += str(j)
147  yield coupling
148  if process == 0 and not usevpvp:
149  yield "ghzgs1_prime2"; yield "ghzgs2"; yield "ghzgs3"; yield "ghzgs4"
150  elif process == 2:
151  for i in range(1, 11): yield "b{}".format(i)
152  if not usevpvp:
153  for i in range(1, 5)+[8]: yield "bzgs".format(i)
154  else:
155  assert False, process
156 
157  return {coupling:
158  (random.uniform(0, 1) + random.uniform(0, 1) * 1j)
159  * (10000 if "prime" in coupling else 1)
160  for coupling in couplings()}
161 
162  def testzp_decay(self):
163  self.runVprime(event1_ggH, False, False, TVar.SelfDefine_spin0, TVar.JHUGen, TVar.ZZINDEPENDENT, **self.Vprimekwargs(0, False))
164  def testzp_spin2(self):
165  self.runVprime(event1_ggH, False, False, TVar.SelfDefine_spin2, TVar.JHUGen, TVar.ZZINDEPENDENT, **self.Vprimekwargs(2, False))
166  def testzp_VBF(self):
167  self.runVprime(event2_VBF, True, False, TVar.SelfDefine_spin0, TVar.JHUGen, TVar.JJVBF, **self.Vprimekwargs(60, False))
168  def testzp_ZH(self):
169  self.runVprime(event3_ZH, True, False, TVar.SelfDefine_spin0, TVar.JHUGen, TVar.Lep_ZH, **self.Vprimekwargs(50, False))
170  def testzpzp_decay(self):
171  self.runVprime(event1_ggH, False, True, TVar.SelfDefine_spin0, TVar.JHUGen, TVar.ZZINDEPENDENT, **self.Vprimekwargs(0, True))
172  def testzp_spin2(self):
173  self.runVprime(event1_ggH, False, True, TVar.SelfDefine_spin2, TVar.JHUGen, TVar.ZZINDEPENDENT, **self.Vprimekwargs(2, True))
174  def testzpzp_VBF(self):
175  self.runVprime(event2_VBF, True, True, TVar.SelfDefine_spin0, TVar.JHUGen, TVar.JJVBF, **self.Vprimekwargs(60, True))
176  def testzpzp_ZH(self):
177  self.runVprime(event3_ZH, True, True, TVar.SelfDefine_spin0, TVar.JHUGen, TVar.Lep_ZH, **self.Vprimekwargs(50, True))
178  def runVprime(self, event, isprod, usevpvp, *setprocessargs, **couplings):
179  """
180  test that with mZ' == mZ, we get the same behavor with Z and Z'
181  """
182 
183  M_Z = 91.1876
184  Ga_Z = 2.4952
185  M_W = 80.399
186  Ga_W = 2.085
187  sitW2 = 0.23119
188 
189  aL_lep = -0.53762
190  aR_lep = 0.46238
191  aL_nu = 1
192  aR_nu = 0
193 
194  sitW2 = 0.23119
195  aL_up = (-2 * sitW2 * 2./3 + 1)
196  aR_up = -2 * sitW2 * 2./3
197  aL_dn = (-2 * sitW2 * -1./3 - 1)
198  aR_dn = -2 * sitW2 * -1./3
199 
200  bL = (2*(1-sitW2))**.5
201 
202  m = self.m
203  computeP = m.computeProdP if isprod else m.computeP
204  m.setInputEvent_fromLHE_Hwithdecay(event, True)
205 
206  m.setProcess(*setprocessargs)
207  for coupling, value in couplings.iteritems():
208  if re.match("(ghz[1-4](_prime[2-7]?)?|b([1-9]|10))", coupling):
209  pass
210  elif re.match("(ghzgs([2-4]|1_prime2)|bzgs[12348])", coupling):
211  if usevpvp:
212  raise ValueError("Can't use gs coupling {} for vpvp test".format(coupling))
213  else:
214  raise ValueError("Invalid coupling {}".format(coupling))
215  setattr(m, coupling, value)
216 
217  me1 = computeP()
218 
219  m.setProcess(*setprocessargs)
220 
221  m.a1 = 1
222  for coupling, value in couplings.iteritems():
223  if usevpvp:
224  newcoupling = coupling.replace("ghz", "ghzpzp").replace("b", "bzpzp")
225  else:
226  newcoupling = coupling.replace("ghz", "ghzzp").replace("b", "bzzp").replace("zzpgs", "zpgs").replace("zzpzgs", "zpgs")
227  if "gs" not in coupling:
228  value /= 2
229  assert newcoupling != coupling
230  setattr(m, newcoupling, value)
231  m.M_Zprime = M_Z
232  m.Ga_Zprime = Ga_Z
233  m.M_Wprime = M_W
234  m.Ga_Wprime = Ga_W
235 
236  m.ezp_El_left = m.ezp_Mu_left = m.ezp_Ta_left = aL_lep
237  m.ezp_El_right = m.ezp_Mu_right = m.ezp_Ta_right = aR_lep
238  m.ezp_Up_left = m.ezp_Chm_left = m.ezp_Top_left = aL_up
239  m.ezp_Up_right = m.ezp_Chm_right = m.ezp_Top_right = aR_up
240  m.ezp_Dn_left = m.ezp_Str_left = m.ezp_Bot_left = aL_dn
241  m.ezp_Dn_right = m.ezp_Str_right = m.ezp_Bot_right = aR_dn
242  m.ezp_NuE_left = aL_nu
243  m.ezp_NuE_right = aR_nu
244 
245  m.ewp_El_left = m.ewp_Mu_left = m.ewp_Ta_left = m.ewp_Up_left = m.ewp_Chm_left = m.ewp_Top_left = bL
246 
247  me2 = computeP()
248 
249  self.assertEquals(me1, me2)
250  self.assertNotEquals(me1, 0)
251 
252  def testmultiplecalls_VBF(self):
253  self.runmultiplecalls(event2_VBF)
254  def testmultiplecalls_ZH(self):
255  self.runmultiplecalls(event3_ZH.replace(" 13", " 3").replace("-13", " -3"))
256 
257  def runmultiplecalls(self, event):
258  m = self.m
259  m.setInputEvent_fromLHE_Hwithdecay(event, False)
260 
261  firstrun = []
262  secondrun = []
263 
264  for lst in firstrun, secondrun:
265  m.setProcess(TVar.SelfDefine_spin0, TVar.JHUGen, TVar.JJVBF)
266  m.ghz1 = 1
267  lst.append(m.computeProdP(True))
268  lst.append(m.getConstant())
269 
270  m.setProcess(TVar.SelfDefine_spin0, TVar.JHUGen, TVar.JJQCD)
271  m.ghg2 = 1
272  lst.append(m.computeProdP(True))
273  lst.append(m.getConstant())
274 
275  m.setProcess(TVar.SelfDefine_spin0, TVar.JHUGen, TVar.Had_ZH)
276  m.ghz3 = 1
277  lst.append(m.computeProdP(True))
278  lst.append(m.getConstant())
279 
280  m.setProcess(TVar.SelfDefine_spin0, TVar.JHUGen, TVar.Had_WH)
281  m.ghz2_prime4 = 1
282  lst.append(m.computeProdP(True))
283  lst.append(m.getConstant())
284 
285  m.setProcess(TVar.SelfDefine_spin0, TVar.JHUGen, TVar.ZZINDEPENDENT)
286  m.ghzzp1 = ezp_El_left = 1
287  lst.append(m.computeP(True))
288  lst.append(m.getConstant())
289 
290  self.assertEquals(firstrun, secondrun)
291 
292 if __name__ == "__main__":
293  unittest.main(argv=[sys.argv[0]]+args.unittest_args)
testME_more.TestMela.testcontact_ZH
def testcontact_ZH(self)
Definition: testME_more.py:85
testME_more.TestMela.runcontact
def runcontact(self, event, isprod, *setprocessargs)
Definition: testME_more.py:87
testME_more.TestMela.Vprimekwargs
def Vprimekwargs(self, process, usevpvp)
Definition: testME_more.py:138
testME_more.TestMela.tearDown
def tearDown(self)
Definition: testME_more.py:79
testME_more.TestMela
Definition: testME_more.py:73
testME_more.TestMela.testmultiplecalls_VBF
def testmultiplecalls_VBF(self)
Definition: testME_more.py:252
testME_more.TestMela.testzpzp_ZH
def testzpzp_ZH(self)
Definition: testME_more.py:176
testME_more.TestMela.setUp
def setUp(self)
Definition: testME_more.py:77
mela.Mela
Definition: mela.py:49
testME_more.TestMela.runVprime
def runVprime(self, event, isprod, usevpvp, *setprocessargs, **couplings)
Definition: testME_more.py:178
testME_more.TestMela.setUpClass
def setUpClass(cls)
Definition: testME_more.py:75
testME_more.TestMela.testzp_VBF
def testzp_VBF(self)
Definition: testME_more.py:166
testME_more.TestMela.testcontact_decay
def testcontact_decay(self)
Definition: testME_more.py:81
testME_more.TestMela.testzp_ZH
def testzp_ZH(self)
Definition: testME_more.py:168
testME_more.TestMela.testzp_spin2
def testzp_spin2(self)
Definition: testME_more.py:164
testME_more.TestMela.testzpzp_VBF
def testzpzp_VBF(self)
Definition: testME_more.py:174
testME_more.TestMela.testzp_decay
def testzp_decay(self)
Definition: testME_more.py:162
testME_more.TestMela.testcontact_VBF
def testcontact_VBF(self)
Definition: testME_more.py:83
testME_more.TestMela.testmultiplecalls_ZH
def testmultiplecalls_ZH(self)
Definition: testME_more.py:254
testME_more.TestMela.runmultiplecalls
def runmultiplecalls(self, event)
Definition: testME_more.py:257
testME_more.TestMela.m
m
Definition: testME_more.py:76
testME_more.TestMela.testzpzp_decay
def testzpzp_decay(self)
Definition: testME_more.py:170
mela.couplings
tuple couplings
Definition: mela.py:791
Generated by   doxygen 1.8.17