1655
1769
return FourMomentum(mom)
1772
class OneNLOWeight(object):
1774
def __init__(self, input):
1777
if isinstance(input, str):
1780
def parse(self, text):
1781
"""parse the line and create the related object"""
1782
#0.546601845792D+00 0.000000000000D+00 0.000000000000D+00 0.119210435309D+02 0.000000000000D+00 5 -1 2 -11 12 21 0 0.24546101D-01 0.15706890D-02 0.12586055D+04 0.12586055D+04 0.12586055D+04 1 2 2 2 5 2 2 0.539995789976D+04
1783
# below comment are from Rik description email
1786
# 1. The first three doubles are, as before, the 'wgt', i.e., the overall event of this
1787
# contribution, and the ones multiplying the log[mu_R/QES] and the log[mu_F/QES]
1788
# stripped of alpha_s and the PDFs.
1789
self.pwgt = [float(f) for f in data[:3]]
1790
# 2. The next two doubles are the values of the (corresponding) Born and
1791
# real-emission matrix elements. You can either use these values to check
1792
# that the newly computed original matrix element weights are correct,
1793
# or directly use these so that you don't have to recompute the original weights.
1794
# For contributions for which the real-emission matrix elements were
1795
# not computed, the 2nd of these numbers is zero. The opposite is not true,
1796
# because each real-emission phase-space configuration has an underlying Born one
1797
# (this is not unique, but on our code we made a specific choice here).
1798
# This latter information is useful if the real-emission matrix elements
1799
# are unstable; you can then reweight with the Born instead.
1800
# (see also point 9 below, where the momentum configurations are assigned).
1801
# I don't think this instability is real problem when reweighting the real-emission
1802
# with tree-level matrix elements (as we generally would do), but is important
1803
# when reweighting with loop-squared contributions as we have been doing for gg->H.
1804
# (I'm not sure that reweighting tree-level with loop^2 is something that
1805
# we can do in general, because we don't really know what to do with the
1806
# virtual matrix elements because we cannot generate 2-loop diagrams.)
1807
self.born = float(data[3])
1808
self.real = float(data[4])
1809
# 3. integer: number of external particles of the real-emission configuration (as before)
1810
self.nexternal = int(data[5])
1811
# 4. PDG codes corresponding to the real-emission configuration (as before)
1812
self.pdgs = [int(i) for i in data[6:6+self.nexternal]]
1813
flag = 6+self.nexternal # new starting point for the position
1814
# 5. next integer is the power of g_strong in the matrix elements (as before)
1815
self.qcdpower = int(data[flag])
1816
# 6. 2 doubles: The bjorken x's used for this contribution (as before)
1817
self.bjks = [float(f) for f in data[flag+1:flag+3]]
1818
# 7. 3 doubles: The Ellis-sexton scale, the renormalisation scale and the factorisation scale, all squared, used for this contribution (as before)
1819
self.scales2 = [float(f) for f in data[flag+3:flag+6]]
1820
# 8.the value of g_strong
1821
self.gs = float(data[flag+6])
1822
# 9. 2 integers: the corresponding Born and real-emission type kinematics. (in the list of momenta)
1823
# Note that also the Born-kinematics has n+1 particles, with, in general,
1824
# one particle with zero momentum (this is not ALWAYS the case,
1825
# there could also be 2 particles with perfectly collinear momentum).
1826
# To convert this from n+1 to a n particles, you have to sum the momenta
1827
# of the two particles that 'merge', see point 12 below.
1828
self.born_related = int(data[flag+7])
1829
self.real_related = int(data[flag+8])
1830
# 10. 1 integer: the 'type'. This is the information you should use to determine
1831
# if to reweight with Born, virtual or real-emission matrix elements.
1832
# (Apart from the possible problems with complicated real-emission matrix elements
1833
# that need to be computed very close to the soft/collinear limits, see point 2 above.
1834
# I guess that for tree-level this is always okay, but when reweighting
1835
# a tree-level contribution with a one-loop squared one, as we do
1836
# for gg->Higgs, this is important).
1837
# type=1 : real-emission:
1839
# type=3 : integrated counter terms:
1840
# type=4 : soft counter-term :
1841
# type=5 : collinear counter-term :
1842
# type=6 : soft-collinear counter-term:
1843
# type=7 : O(alphaS) expansion of Sudakov factor for NNLL+NLO :
1844
# type=8 : soft counter-term (with n+1-body kin.):
1845
# type=9 : collinear counter-term (with n+1-body kin.):
1846
# type=10: soft-collinear counter-term (with n+1-body kin.):
1847
# type=11: real-emission (with n-body kin.):
1848
# type=12: MC subtraction with n-body kin.:
1849
# type=13: MC subtraction with n+1-body kin.:
1850
# type=14: virtual corrections minus approximate virtual
1851
# type=15: approximate virtual corrections:
1852
self.type = int(data[flag+9])
1853
# 11. 1 integer: The FKS configuration for this contribution (not really
1854
# relevant for anything, but is used in checking the reweighting to
1855
# get scale & PDF uncertainties).
1856
self.nfks = int(data[flag+10])
1857
# 12. 2 integers: the two particles that should be merged to form the
1858
# born contribution from the real-emission one. Remove these two particles
1859
# from the (ordered) list of PDG codes, and insert a newly created particle
1860
# at the location of the minimum of the two particles removed.
1861
# I.e., if you merge particles 2 and 4, you have to insert the new particle
1862
# as the 2nd particle. And particle 5 and above will be shifted down by one.
1863
self.to_merge_pdg = [int (f) for f in data[flag+11:flag+13]]
1864
# 13. 1 integer: the PDG code of the particle that is created after merging the two particles at point 12.
1865
self.merge_new_pdg = int(data[flag+13])
1866
# 14. 1 double: the reference number that one should be able to reconstruct
1867
# form the weights (point 1 above) and the rest of the information of this line.
1868
# This is really the contribution to this event as computed by the code
1869
# (and is passed to the integrator). It contains everything.
1870
self.ref_wgt = float(data[flag+14])
1872
#check the momenta configuration linked to the event
1873
if self.type in [1,11]:
1874
self.momenta_config = self.real_related
1876
self.momenta_config = self.born_related
1879
class NLO_PARTIALWEIGHT(object):
1881
class BasicEvent(list):
1883
def __init__(self, momenta, wgts, event):
1884
list.__init__(self, momenta)
1888
self.pdgs = list(wgts[0].pdgs)
1891
if wgts[0].momenta_config == wgts[0].born_related:
1892
# need to remove one momenta.
1893
ind1, ind2 = [ind-1 for ind in wgts[0].to_merge_pdg]
1895
ind1, ind2 = ind2, ind1
1896
if ind1 >= sum(1 for p in event if p.status==-1):
1897
new_p = self[ind1] + self[ind2]
1899
new_p = self[ind1] - self[ind2]
1901
self.insert(ind1, new_p)
1904
self.pdgs.insert(ind1, wgts[0].merge_new_pdg )
1906
# DO NOT update the pdgs of the partial weight!
1907
elif any(w.type in [1,11] for w in wgts):
1908
if any(w.type not in [1,11] for w in wgts):
1910
# check if this is too soft/colinear if so use the born
1911
ind1, ind2 = [ind-1 for ind in wgts[0].to_merge_pdg]
1913
ind1, ind2 = ind2, ind1
1914
if ind1 >= sum(1 for p in event if p.status==-1):
1915
new_p = self[ind1] + self[ind2]
1917
new_p = self[ind1] - self[ind2]
1920
ptot = FourMomentum()
1921
for i in xrange(len(self)):
1926
if ptot.mass_sqr > 1e-16:
1927
misc.sprint(ptot, ptot.mass_sqr)
1929
inv_mass = new_p.mass_sqr
1930
shat = (self[0]+self[1]).mass_sqr
1931
if (abs(inv_mass)/shat < 1e-6):
1932
# misc.sprint(abs(inv_mass)/shat)
1934
self.insert(ind1, new_p)
1937
self.pdgs.insert(ind1, wgts[0].merge_new_pdg )
1939
# DO NOT update the pdgs of the partial weight!
1942
ptot = FourMomentum()
1943
for i in xrange(len(self)):
1948
if ptot.mass_sqr > 1e-16:
1949
misc.sprint(ptot, ptot.mass_sqr)
1952
def get_pdg_code(self):
1955
def get_tag_and_order(self):
1956
""" return the tag and order for this basic event"""
1957
(initial, _), _ = self.event.get_tag_and_order()
1958
order = self.get_pdg_code()
1961
initial, out = order[:len(initial)], order[len(initial):]
1964
return (tuple(initial), tuple(out)), order
1966
def get_momenta(self, get_order, allow_reversed=True):
1967
"""return the momenta vector in the order asked for"""
1969
#avoid to modify the input
1970
order = [list(get_order[0]), list(get_order[1])]
1971
out = [''] *(len(order[0])+len(order[1]))
1972
pdgs = self.get_pdg_code()
1973
for pos, part in enumerate(self):
1974
if pos < len(get_order[0]): #initial
1976
ind = order[0].index(pdgs[pos])
1977
except ValueError, error:
1978
if not allow_reversed:
1981
order = [[-i for i in get_order[0]],[-i for i in get_order[1]]]
1983
return self.get_momenta(order, False)
1992
ind = order[1].index(pdgs[pos])
1993
except ValueError, error:
1994
if not allow_reversed:
1997
order = [[-i for i in get_order[0]],[-i for i in get_order[1]]]
1999
return self.get_momenta(order, False)
2002
position = len(order[0]) + ind
2005
out[position] = (part.E, part.px, part.py, part.pz)
2010
def get_helicity(self, *args):
2011
return [9] * len(self)
2015
return self.event.aqcd
2017
def __init__(self, input, event):
2020
if isinstance(input, str):
2024
def parse(self, text):
2025
"""create the object from the string information (see example below)"""
2026
#0.2344688900d+00 8 2 0
2027
#0.4676614699d+02 0.0000000000d+00 0.0000000000d+00 0.4676614699d+02
2028
#0.4676614699d+02 0.0000000000d+00 0.0000000000d+00 -.4676614699d+02
2029
#0.4676614699d+02 0.2256794794d+02 0.4332148227d+01 0.4073073437d+02
2030
#0.4676614699d+02 -.2256794794d+02 -.4332148227d+01 -.4073073437d+02
2031
#0.0000000000d+00 -.0000000000d+00 -.0000000000d+00 -.0000000000d+00
2032
#0.4780341163d+02 0.0000000000d+00 0.0000000000d+00 0.4780341163d+02
2033
#0.4822581633d+02 0.0000000000d+00 0.0000000000d+00 -.4822581633d+02
2034
#0.4729127470d+02 0.2347155377d+02 0.5153455534d+01 0.4073073437d+02
2035
#0.4627255267d+02 -.2167412893d+02 -.3519736379d+01 -.4073073437d+02
2036
#0.2465400591d+01 -.1797424844d+01 -.1633719155d+01 -.4224046944d+00
2037
#0.473706252575d-01 0.000000000000d+00 0.000000000000d+00 5 -3 3 -11 11 21 0 0.11849903d-02 0.43683926d-01 0.52807978d+03 0.52807978d+03 0.52807978d+03 1 2 1 0.106660059627d+03
2038
#-.101626389492d-02 0.000000000000d+00 -.181915673961d-03 5 -3 3 -11 11 21 2 0.11849903d-02 0.43683926d-01 0.52807978d+03 0.52807978d+03 0.52807978d+03 1 3 1 -.433615206719d+01
2039
#0.219583436285d-02 0.000000000000d+00 0.000000000000d+00 5 -3 3 -11 11 21 2 0.11849903d-02 0.43683926d-01 0.52807978d+03 0.52807978d+03 0.52807978d+03 1 15 1 0.936909375537d+01
2040
#0.290043597283d-03 0.000000000000d+00 0.000000000000d+00 5 -3 3 -11 11 21 2 0.12292838d-02 0.43683926d-01 0.58606724d+03 0.58606724d+03 0.58606724d+03 1 12 1 0.118841547979d+01
2041
#-.856330613460d-01 0.000000000000d+00 0.000000000000d+00 5 -3 3 -11 11 21 2 0.11849903d-02 0.43683926d-01 0.52807978d+03 0.52807978d+03 0.52807978d+03 1 4 1 -.365375546483d+03
2042
#0.854918237609d-01 0.000000000000d+00 0.000000000000d+00 5 -3 3 -11 11 21 2 0.12112732d-02 0.45047393d-01 0.58606724d+03 0.58606724d+03 0.58606724d+03 2 11 1 0.337816057347d+03
2043
#0.359257891118d-05 0.000000000000d+00 0.000000000000d+00 5 21 3 -11 11 3 2 0.12292838d-02 0.43683926d-01 0.58606724d+03 0.58606724d+03 0.58606724d+03 1 12 3 0.334254554762d+00
2044
#0.929944817736d-03 0.000000000000d+00 0.000000000000d+00 5 21 3 -11 11 3 2 0.12112732d-02 0.45047393d-01 0.58606724d+03 0.58606724d+03 0.58606724d+03 2 11 3 0.835109616010d+02
2046
text = text.lower().replace('d','e')
2047
all_line = text.split('\n')
2048
#get global information
2050
while not first_line.strip():
2051
first_line = all_line.pop(0)
2053
wgt, nb_wgt, nb_event, _ = first_line.split()
2054
nb_wgt, nb_event = int(nb_wgt), int(nb_event)
2058
for line in all_line:
2061
p = FourMomentum(data)
2064
wgt = OneNLOWeight(line)
2067
assert len(wgts) == int(nb_wgt)
2069
get_weights_for_momenta = {}
2072
if wgt.momenta_config in get_weights_for_momenta:
2073
get_weights_for_momenta[wgt.momenta_config].append(wgt)
2075
if size_momenta == 0: size_momenta = wgt.nexternal
2076
assert size_momenta == wgt.nexternal
2077
get_weights_for_momenta[wgt.momenta_config] = [wgt]
2079
assert sum(len(c) for c in get_weights_for_momenta.values()) == int(nb_wgt), "%s != %s" % (sum(len(c) for c in get_weights_for_momenta.values()), nb_wgt)
2084
for key in range(1, nb_event+1):
2085
if key in get_weights_for_momenta:
2086
wgt = get_weights_for_momenta[key]
2087
evt = self.BasicEvent(momenta[:size_momenta], get_weights_for_momenta[key], self.event)
2088
self.cevents.append(evt)
2089
momenta = momenta[size_momenta:]
2092
for cevt in self.cevents:
2093
nb_wgt_check += len(cevt.wgts)
2094
assert nb_wgt_check == int(nb_wgt)
1658
2098
if '__main__' == __name__:
1660
2100
# Example 1: adding some missing information to the event (here distance travelled)
added
removed
madgraph/various/lhe_parser.py
