35def resample_meson(input_filename, output_filename, meson_id, meson_width ):
36 event = False
37 findW=False
38
39 with open(input_filename, 'r') as input_file:
40 with open(output_filename, 'w') as output_file:
41 for line in input_file:
42 Q2=0
43 if not event and '<event>' not in line:
44 output_file.write(line)
45 continue
46
47 if not event and '<event>' in line:
48 output_file.write(line)
49 event = True
50 findW = False
51 continue
52
53 if '<' in line or '>' in line:
54 event = False
55 if line.startswith('<mgrwt>') or line.startswith('</mgrwt>'):
56 output_file.write(line)
57 continue
58
59 if abs(int(line.split()[0])) == 24:
60
61 findW=True
62
63
64
65 if (abs(int(line.split()[0])) == abs(int(meson_id)) and abs(int(line.split()[2])) == 4) or (not findW and abs(int(line.split()[0])) == abs(int(meson_id)) and abs(int(line.split()[2]))==3 ):
66
67
68 line_parts = line.split()
69 line_parts[0] = str(int(line.split()[0]))
70
71
72
73
74 momentum_x = float(line.split()[6])
75 momentum_y = float(line.split()[7])
76 momentum_z= float(line.split()[8])
77 Energy = float(line.split()[9])
78
79 meson=build_TLorentz(Energy,momentum_x, momentum_y, momentum_z)
80
81 continue
82
83 if ((abs(int(line.split()[0])) == 13 or abs(int(line.split()[0])) ==11) and int(line.split()[2])==4) or ( not findW and (abs(int(line.split()[0])) == 13 or abs(int(line.split()[0])) ==11) and int(line.split()[2])==3):
84
85 E_l = float(line.split()[9])
86 px_l= float(line.split()[6])
87 py_l= float(line.split()[7])
88 pz_l= float(line.split()[8])
89
90 lep1=build_TLorentz(E_l,px_l, py_l, pz_l)
91
92
93 combined = meson+lep1
94
95 beta=combined.BoostVector()
96 meson.Boost(-beta)
97 lep1.Boost(-beta)
98
99
100
101 combined_B =meson+lep1
102
103 meson_massp =meson.M()
104
105
106 theta = math.acos(meson.Pz() / math.sqrt(meson.Px()**2 + meson.Py()**2 + meson.Pz()**2))
107
108 phi = math.atan2(meson.Py(), meson.Px())
109
110 Q2 = generate_virtuality(meson_massp, float(meson_width))
111
112 while (Q2 > meson.E()*meson.E()+meson.Px()**2 + meson.Py()**2 + meson.Pz()**2 or not (0.3 < math.sqrt(Q2) < 1.5)):
113
114 Q2 = generate_virtuality(meson_massp, float(meson_width))
115
116 meson_energy = combined_B.E()/2 +((Q2-lep1.M()**2)/(2*combined_B.E()))
117
118 p3_magnitude=calculate_momentum_magnitude(meson_energy, Q2)
119 four_momentum = build_four_momentum(meson_energy, p3_magnitude, theta, phi)
120 meson.SetPxPyPzE(four_momentum[1],four_momentum[2],four_momentum[3],four_momentum[0])
121
122 px_lep= combined_B.Px()-meson.Px()
123 py_lep= combined_B.Py()-meson.Py()
124 pz_lep= combined_B.Pz()-meson.Pz()
125
126 lep1_energy = math.sqrt(px_lep**2+py_lep**2+pz_lep**2+lep1.M()**2)
127 lep1.SetPxPyPzE(px_lep,py_lep,pz_lep,lep1_energy)
128
129
130 meson.Boost(beta)
131 lep1.Boost(beta)
132
133 line_parts[9] = str("{:+.10e}".format(meson.E()))
134 line_parts[6] =str("{:+.10e}".format(meson.Px()))
135 line_parts[7] =str("{:+.10e}".format(meson.Py()))
136 line_parts[8] =str("{:+.10e}".format(meson.Pz()))
137 line_parts[10] =str("{:+.10e}".format(meson.M()))
138
139 updated_line = ' '.join(line_parts) + '\n'
140 output_file.write(' ' +''.join(updated_line) )
141
142 if ((abs(int(line.split()[0])) == 13 or abs(int(line.split()[0])) ==11) and int(line.split()[2])==4) or (not findW and (abs(int(line.split()[0])) == 13 or abs(int(line.split()[0])) ==11) and int(line.split()[2])==3):
143 line_parts_l = line.split()
144 line_parts_l[9] = str("{:+.10e}".format(lep1.E()))
145 line_parts_l[6] =str("{:+.10e}".format(lep1.Px()))
146 line_parts_l[7] =str("{:+.10e}".format(lep1.Py()))
147 line_parts_l[8] =str("{:+.10e}".format(lep1.Pz()))
148 line_parts_l[10] =str("{:+.10e}".format(lep1.M()))
149
150 updated_line_l = ' '.join(line_parts_l) + '\n'
151 output_file.write(' ' +''.join(updated_line_l) + '\n')
152
153
154 else:
155 output_file.write(line)
156
157
158
1.14.0