forked from alisw/POWHEG
-
Notifications
You must be signed in to change notification settings - Fork 0
/
LesHouches.f
479 lines (464 loc) · 13.7 KB
/
LesHouches.f
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
subroutine gen_leshouches
implicit none
include 'pwhg_math.h'
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_flg.h'
include 'pwhg_kn.h'
include 'pwhg_rad.h'
include 'LesHouches.h'
integer alr,em,rad,flem,flrad,coluborn(2)
if(kn_csi.eq.0d0) then
c it is a Born event or a LO event
c first assign flavours and colour to the underlying Born process
call momenta_lh(kn_pborn,nlegborn)
call born_lh
if (flg_LOevents) then
c set default scalup to the sqrt of the center-of-mass energy of the event.
c For processes that need different values, the scalup value can be
c overridden in the finalize_lh user subroutine.
scalup=2*kn_cmpborn(0,1)
else
scalup=sqrt(rad_ptsqmin)
endif
else
c it is an event with radiation
call momenta_lh(kn_preal,nlegreal)
call born_lh
scalup=sqrt(rad_pt2max)
nup=nlegreal
istup(nup)=1
spinup(nup)=9
vtimup(nup)=0
mothup(1,nup)=1
mothup(2,nup)=2
alr=rad_realalr
em=flst_emitter(alr)
rad=nlegreal
if(em.eq.0) then
if(kn_y.gt.0) then
em=1
else
em=2
endif
endif
flem=flst_alr(em,alr)
flrad=flst_alr(rad,alr)
coluborn(1)=icolup(1,em)
coluborn(2)=icolup(2,em)
if(em.le.2) then
c setcolour_rad works with all incoming (or all outgoing) flavours and colours;
c for ISR make everything outgoing:
c thus in input change the sign of the emitter flavour, to make it outgoing
c and conjugate its colour in output, to make it incoming.
call setcolour_rad(coluborn,-flem,flrad,
# icolup(1,em),icolup(1,rad))
call colour_conj(icolup(1,em))
else
c For ISR make everything incoming;
c Change the sign of the emitter and radiated flavours in input,
c to make it incoming;
c conjugate their colours in the output, to make them outgoing.
call setcolour_rad(coluborn,-flem,-flrad,
# icolup(1,em),icolup(1,rad))
call colour_conj(icolup(1,rad))
call colour_conj(icolup(1,em))
endif
if(flem.eq.0) flem=21
if(flrad.eq.0) flrad=21
idup(em)=flem
idup(rad)=flrad
endif
c if there are resonances, setup up appropriate LH pointers
call lh_resonances
c add resonances, perform decays, put particles on shell, etc.(or nothing!)
call finalize_lh
end
c i1<i2
subroutine add_resonance(idpdg,i1,i2)
implicit none
include 'LesHouches.h'
real * 8 iduptmp(maxnup),istuptmp(maxnup),mothuptmp(2,maxnup),
# icoluptmp(2,maxnup),puptmp(5,maxnup),vtimuptmp(maxnup),
# spinuptmp(maxnup)
integer idpdg,i1,i2,i,j
if (i1.ge.i2) then
write(*,*) 'wrong sequence in add_resonance'
stop
endif
do i=i1,nup
iduptmp(i) = idup(i)
istuptmp(i) = istup(i)
mothuptmp(1,i) = mothup(1,i)
mothuptmp(2,i) = mothup(2,i)
icoluptmp(1,i) = icolup(1,i)
icoluptmp(2,i) = icolup(2,i)
do j=1,5
puptmp(j,i)=pup(j,i)
enddo
vtimuptmp(i)=vtimup(i)
spinuptmp(i)=spinup(i)
enddo
idup(i1)=idpdg
istup(i1)=+2
mothup(1,i1)=1
mothup(2,i1)=2
icolup(1,i1)=0
icolup(2,i1)=0
do j=1,4
pup(j,i1)=puptmp(j,i1)+puptmp(j,i2)
enddo
pup(5,i1)=sqrt(pup(4,i1)**2-
# (pup(1,i1)**2+pup(2,i1)**2+pup(3,i1)**2))
vtimup(i1)=0
spinup(i1)=9
c change mothers of decaying particles
mothuptmp(1,i1)=i1
mothuptmp(2,i1)=i1
mothuptmp(1,i2)=i1
mothuptmp(2,i2)=i1
nup=nup+1
do i=i1+1, nup
idup(i) = iduptmp(i-1)
istup(i) = istuptmp(i-1)
mothup(1,i) = mothuptmp(1,i-1)
mothup(2,i) = mothuptmp(2,i-1)
icolup(1,i) = icoluptmp(1,i-1)
icolup(2,i) = icoluptmp(2,i-1)
do j=1,5
pup(j,i)=puptmp(j,i-1)
enddo
vtimup(i)=vtimuptmp(i-1)
spinup(i)=spinuptmp(i-1)
enddo
end
subroutine colour_conj(icol)
implicit none
integer icol(2)
integer itmp
itmp=icol(2)
icol(2)=icol(1)
icol(1)=itmp
end
subroutine momenta_lh(p,n)
implicit none
integer n
real * 8 p(0:3,n)
integer k,mu
include 'LesHouches.h'
do k=1,n
do mu=1,3
pup(mu,k)=p(mu,k)
enddo
pup(4,k)=p(0,k)
pup(5,k)=sqrt(abs(p(0,k)**2-p(1,k)**2
# -p(2,k)**2-p(3,k)**2))
enddo
end
subroutine born_lh
implicit none
include 'pwhg_math.h'
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_kn.h'
include 'pwhg_rad.h'
include 'LesHouches.h'
c include 'PhysPars.h'
include 'pwhg_st.h'
integer k,fl
nup=nlegborn
c id of the event
idprup=lprup(1)
ccccc CAVEAT: aqedup must not be set here
ccccc since otherwise a dependence on the
ccccc process-dependent PhysPars.h must be
ccccc introduced
aqedup=-1 ! ph_alphaem
cccccccccccccccccccccccccccccccccccccccccccccccccccccc
aqcdup=st_alpha
do k=1,nlegborn
fl=flst_born(k,rad_ubornidx)
c gluons are numbered 21 in pdg
if(fl.eq.0) fl=21
idup(k)=fl
if(k.gt.2) then
istup(k)=1
mothup(1,k)=1
mothup(2,k)=2
else
istup(k)=-1
mothup(1,k)=0
mothup(2,k)=0
endif
spinup(k)=9
vtimup(k)=0
enddo
call borncolour_lh
end
subroutine lh_resonances
implicit none
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_rad.h'
include 'LesHouches.h'
integer moth,j,iub,em
iub = rad_ubornidx
do j=1,nlegborn
moth=flst_bornres(j,iub)
if(moth.ne.0) then
mothup(1,j)=moth
mothup(2,j)=moth
istup(moth)=2
endif
enddo
if(nup.eq.nlegreal) then
em = flst_emitter(rad_realalr)
if(em.ne.0) then
moth=flst_bornres(flst_emitter(rad_realalr),iub)
if(moth.ne.0) then
mothup(1,nlegreal)=moth
mothup(2,nlegreal)=moth
endif
endif
endif
end
subroutine getnewcolor(newcolor)
implicit none
integer newcolor
include 'LesHouches.h'
integer j,k
newcolor=511
1 continue
do j=1,nup
do k=1,2
if(icolup(k,j).eq.newcolor) then
newcolor=newcolor+1
goto 1
endif
enddo
enddo
end
c Given a vertex with 3 partons, the first parton with incoming colour
c colin, the second and third partons with incoming flavours fl2, fl3,
c it assigns the incoming colours of the second and third partons
c colout2, colout3. The colour assignment is unambiguous in all cases, but in
c the one where both fl2 and fl3 are gluons, where a 50% random choice is made
c over the two possible colour connections.
subroutine setcolour_rad(colin,fl2,fl3,colout2,colout3)
implicit none
integer colin(2),fl2,fl3,colout2(2),colout3(2)
integer newcolor
real * 8 random
logical is_coloured
external random,is_coloured
call getnewcolor(newcolor)
c First handle photon case
if(fl2.eq.22) then
colout2=0
colout3(1)=colin(2)
colout3(2)=colin(1)
elseif(fl3.eq.22) then
colout3=0
colout2(1)=colin(2)
colout2(2)=colin(1)
elseif(colin(1).eq.0.and.colin(2).eq.0) then
if(is_coloured(fl2)) then
if(fl2.gt.0) then
colout2(1)=newcolor
colout2(2)=0
colout3(2)=newcolor
colout3(1)=0
elseif(fl3.gt.0) then
colout2(1)=0
colout2(2)=newcolor
colout3(2)=0
colout3(1)=newcolor
else
goto 998
endif
else
colout2=0
colout3=0
endif
else
c now handle all coloured particles
if(colin(1).ne.0.and.colin(2).ne.0) then
if(fl2.eq.0.and.fl3.eq.0) then
if(random().gt.0.5d0) then
colout2(1)=colin(2)
colout2(2)=newcolor
colout3(1)=newcolor
colout3(2)=colin(1)
else
colout3(1)=colin(2)
colout3(2)=newcolor
colout2(1)=newcolor
colout2(2)=colin(1)
endif
elseif(fl2.gt.0.and.fl3.lt.0) then
colout2(2)=0
colout2(1)=colin(2)
colout3(1)=0
colout3(2)=colin(1)
elseif(fl2.lt.0.and.fl3.gt.0) then
colout2(1)=0
colout2(2)=colin(1)
colout3(2)=0
colout3(1)=colin(2)
else
goto 998
endif
elseif(colin(2).eq.0) then
if(fl2.eq.0) then
colout3(1)=0
colout3(2)=newcolor
colout2(1)=newcolor
colout2(2)=colin(1)
elseif(fl3.eq.0) then
colout2(1)=0
colout2(2)=newcolor
colout3(1)=newcolor
colout3(2)=colin(1)
else
goto 998
endif
elseif(colin(1).eq.0) then
if(fl2.eq.0) then
colout3(1)=newcolor
colout3(2)=0
colout2(2)=newcolor
colout2(1)=colin(2)
elseif(fl3.eq.0) then
colout2(1)=newcolor
colout2(2)=0
colout3(1)=colin(2)
colout3(2)=newcolor
else
goto 998
endif
else
goto 998
endif
endif
return
998 continue
write(*,*)
# ' Error in setcolour_rad: inconsistent colour flavour input'
write(*,*) 'colin[1:2] ',colin(1),' ',colin(2)
write(*,*) 'fl2, fl3 ',fl2,' ',fl3
write(*,*) 'newcolor ',newcolor
call pwhg_exit(1)
end
subroutine displeshouches
include 'LesHouches.h'
integer mark(200),icol(200),j,k,ilist,icur
write(*,*) 'incoming beams:', idbmup(1), idbmup(2)
write(*,*) 'number of partons in subprocess',nup
write(*,'(a,10(i3,1x))') 'ist:',(istup(j),j=1,nup)
write(*,'(a,10(i3,1x))') 'ids:',(idup(j),j=1,nup)
do j=1,nup
write(*,'(5(d10.4,1x))') (pup(k,j),k=1,5)
enddo
c check that quarks have zero anticolour, and antiquarks
c have zero colour
do j=1,nup
if(idup(j).ge.1.and.idup(j).le.6.and.icolup(2,j).ne.0) then
write(*,*) ' error: anticolor in quark'
elseif(idup(j).ge.-6.and.idup(j).le.-1.and.icolup(1,j).ne.0)
# then
write(*,*) ' error: color in quark'
endif
enddo
c Print color linked lists
c first conjugate incoming colors
call colour_conj(icolup(1,1))
call colour_conj(icolup(1,2))
do j=1,nup
mark(j)=1
enddo
imarked=0
do j=1,nup
if(istup(j).eq.2.or.
# (icolup(1,j).eq.0.and.icolup(2,j).eq.0)) then
mark(j)=0
imarked=imarked+1
endif
if(istup(j).ne.2.and.idup(j).gt.0.and.idup(j).lt.6) then
icur=j
endif
enddo
mark(icur)=0
ilist=1
imarked=imarked+1
icol(ilist)=icur
12 continue
do j=1,nup
if(istup(j).ne.2.and.
# (icolup(1,j).ne.0.or.icolup(2,j).ne.0)) then
if(mark(j).ne.0.and.icolup(1,icur).eq.icolup(2,j)) then
ilist=ilist+1
imarked=imarked+1
icol(ilist)=j
mark(j)=0
icur=j
if(imarked.eq.nup) goto 22
if(icolup(1,icur).eq.0) then
do k=1,nup
if(mark(k).ne.0.and.icolup(2,k).eq.0) then
ilist=ilist+2
icol(ilist)=k
icol(ilist-1)=0
imarked=imarked+1
icur=k
mark(k)=0
goto 12
endif
enddo
write(*,*) ' inconsistent colors!'
endif
goto 12
endif
endif
enddo
write(*,*) ' inconsistent colors!'
22 continue
call colour_conj(icolup(1,1))
call colour_conj(icolup(1,2))
write(*,*) (icol(j),j=1,ilist)
end
subroutine pdfreweightinfo(id1,id2,x1,x2,xmufact,xf1,xf2)
implicit none
integer id1,id2
real * 8 x1,x2,xf1,xf2,xmufact
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_kn.h'
include 'pwhg_rad.h'
include 'pwhg_pdf.h'
include 'pwhg_st.h'
include 'LesHouches.h'
real * 8 pdf(-pdf_nparton:pdf_nparton)
if(rad_type.eq.1) then
c Btilde event: pass x1 and x2, id1, id2 etc. of the current underlying Born
id1=flst_born(1,rad_ubornidx)
id2=flst_born(2,rad_ubornidx)
x1=kn_xb1
x2=kn_xb2
call setscalesbtilde
elseif(rad_type.eq.2.or.rad_type.eq.3) then
id1=idup(1)
id2=idup(2)
if(id1.eq.21) id1=0
if(id2.eq.21) id2=0
x1=kn_x1
x2=kn_x2
call setscalesbtilde
endif
xmufact=sqrt(st_mufact2)
call pdfcall(1,x1,pdf)
xf1=x1*pdf(id1)
call pdfcall(2,x2,pdf)
xf2=x2*pdf(id2)
if(id1.eq.0) id1=21
if(id2.eq.0) id2=21
end