~maddevelopers/mg5amcnlo/WWW5_caching

li = [self.no, self.id, self.status,self.mo1, self.mo2, self.co1, self.co2, self.mom.px,self.mom.py, self.mom.pz, self.mom.E, self.mom.m, self.life, self.polar]

print "%2i | %9i | %4i | %4i %4i | %4i %4i | %18.10E %18.10E %18.10E %18.10E %18.10E | %1.0f. %2.0f" % tuple(li)

def writeMe(self):

li = [self.id, self.status,self.mo1, self.mo2, self.co1, self.co2, self.mom.px,self.mom.py, self.mom.pz, self.mom.E, self.mom.m, self.life, self.polar]

return "%9i %4i %4i %4i %4i %4i %18.10E %18.10E %18.10E %18.10E %18.10E %1.0f. %2.0f\n" % tuple(li)

#useful function for converting a string to variables

def parseStringToVars(input):

if input.find(".")>-1 :

return float(input)

else:

100

return int(input)

101

102

#main part of analysis

103

if len(sys.argv)!=3:

104

print "Usage: addmasses.py <infile> <outfile> "

105

print " Last modified: Fri Nov 21 10:49:14 CST 2008 "

106

sys.exit(0)

107

else:

108

print "Running addmasses.py to add masses and correct kinematics of fixed particles"

109

110

#first print out leading information

111

try:

112

f=open(sys.argv[1],'r')

113

except IOError:

114

print "need a file for reading"

115

sys.exit(0)

116

117

try:

118

g=open(sys.argv[2],'w')

119

except IOError:

120

print "need a file for writing"

121

sys.exit(0)

122

123

while 1:

124

try:

125

line=f.readline()

126

except IOError:

127

print "Problem reading from file ",sys.argv[1]

128

sys.exit(0)

129

if line.find("<event>")==-1:

130

g.write(line)

131

else:

132

break

133

134

f.close()

135

136

#let xml find the event tags

137

try:

138

xmldoc = minidom.parse(sys.argv[1])

139

except IOError:

140

print " could not open file for xml parsing ",sys.argv[1]

141

sys.exit(0)

142

143

144

reflist = xmldoc.getElementsByTagName('event')

145

146

for ref in reflist:

147

lines = ref.toxml()

148

slines = lines.split("\n")

149

next = 0

150

nup = 0

151

nlines = len(slines)

152

counter = 0

153

event = []

154

event_description=""

155

event_poundSign=""

156

157

while counter<nlines:

158

s=slines[counter]

159

if s.find("<event>")>-1:

160

next=1

161

elif s.find("</event>")>-1:

162

pass

163

elif s.find("#")>-1:

164

event_poundSign=s

165

elif next==1:

166

event_description=s

167

next=0

168

elif not s:

169

continue

170

else:

171

t=[]

172

for l in s.split(): t.append(parseStringToVars(l))

173

nup = nup+1

174

event.append(Particle(nup,t))

175

counter=counter+1

176

177

#default is to skip this

178

if motherFlag:

179

180

noMotherList=[]

181

for p in event:

182

if p.status == -1: continue

183

idabs = abs(p.id)

184

idmo = p.mo1

185

pmo = event[idmo-1]

186

if idabs>=11 and idabs<=16:

187

if p.mo1==1 or (pmo.co1 !=0 or pmo.co2 !=0):

188

noMotherList.append(p.no)

189

elif idabs<=5 and abs(pmo.id)==6:

190

if not ( p.co1==pmo.co1 and p.co2==pmo.co2):

191

noMotherList.append(p.no)

192

193

nAdded=0

194

if len(noMotherList)==0:

195

pass

196

elif len(noMotherList)%2 != 0:

197

print "single orphan; do not know how to process"

198

else:

199

ki=0

200

while ki<len(noMotherList)-1:

201

l1=event[noMotherList[ki]-1]

202

ki=ki+1

203

l2=event[noMotherList[ki]-1]

204

lq=l1.id + l2.id

205

if lq==-1 or lq==1:

206

idMom=24*lq

207

elif lq==0:

208

idMom=23

209

else:

210

break

211

lMom=l1.mom+l2.mom

212

lm=[idMom,2,l1.mo1,l1.mo2,0,0,lMom.px,lMom.py,lMom.pz,lMom.E,lMom.calcMass(),0,0]

213

nup=nup+1

214

nAdded=nAdded+1

215

event.append(Particle(nup,lm))

216

l1.mo1=l1.mo2=l2.mo1=l2.mo2=nup

217

ki=ki+1

218

219

# update number of partons if mothers added

220

if nAdded>0:

221

s1=event_description.split()[0]

222

mySub = re.compile(s1)

223

event_description = mySub.sub(str(nup),event_description)

224

225

if nAdded>0:

226

for ip in range(len(event)):

227

l=event[ip]

228

if l.mo1 > ip + 1:

229

nmo=l.mo1

230

event.insert(ip, event.pop(l.mo1-1))

231

event[ip].no = ip + 1

232

for l2 in event[ip + 1:]:

233

if l2.no > ip and l2.no < nmo + 1:

234

l2.no += 1

235

if l2.mo1 == nmo:

236

l2.mo1 = l2.mo2 = ip + 1

237

elif l2.mo1 > ip and l2.mo1 < nmo:

238

l2.mo1 = l2.mo2 = l2.mo1 + 1

239

240

# identify mothers

241

particleDict={}

242

for p in event:

243

idabs = abs(p.id)

244

if idabs>=11 and idabs<=16:

245

if p.mo1==1:

246

pass

247

else:

248

if p.mo1 in particleDict:

249

l=particleDict[p.mo1]

250

l.append(p.no)

251

else:

252

l=[p.no]

253

particleDict[p.mo1]=l

254

255

# repair kinematics

256

for k in particleDict:

257

if len(particleDict[k]) != 2: continue

258

t=particleDict[k]

259

p1=event[t[0]-1]

260

p1.mom.boost(event[k-1].mom,-1)

261

p2=event[t[1]-1]

262

p2.mom.boost(event[k-1].mom,-1)

263

rsh=event[k-1].mom.m

264

if p1.id in massDict: p1.mom.m=massDict[p1.id]

265

if p2.id in massDict: p2.mom.m=massDict[p2.id]

266

p1.mom.E = (rsh*rsh + p1.mom.m**2 - p2.mom.m**2)/(2.0*rsh)

267

pmagOld=math.sqrt(p1.mom.px**2+p1.mom.py**2+p1.mom.pz**2)

268

pmagNew=math.sqrt(p1.mom.E**2-p1.mom.m**2)

269

p1.mom.reScale(pmagOld,pmagNew)

270

p2.mom = Momentum(0,0,0,rsh,rsh) - p1.mom

271

p1.mom.boost(event[k-1].mom,1)

272

p2.mom.boost(event[k-1].mom,1)

273

274

pSum = Momentum(0,0,0,0,0)

275

for p in event:

276

if p.status== 2 :

277

pass

278

elif p.status==-1:

279

pSum = pSum - p.mom

280

elif p.status==1:

281

pSum = pSum + p.mom

282

283

if abs(pSum.px)>toler or abs(pSum.py)>toler or abs(pSum.pz)>toler or abs(pSum.E)>toler:

284

print "Event does not pass tolerance ",toler

285

pSum.printMe()

286

287

if 1:

288

g.write("<event>\n")

289

g.write(event_description+"\n")

290

for p in event:

291

g.write(p.writeMe())

292

g.write(event_poundSign+"\n")

293

g.write("</event>\n")

294

#at the end

295

g.write("</LesHouchesEvents>\n")

296

g.close()

297

298

299

Older »