forked from alisw/POWHEG
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathpwhg_init.f
208 lines (206 loc) · 7.24 KB
/
pwhg_init.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
subroutine pwhginit
implicit none
include 'LesHouches.h'
include 'nlegborn.h'
include 'pwhg_flst.h'
include 'pwhg_kn.h'
include 'pwhg_st.h'
include 'pwhg_em.h'
include 'pwhg_pdf.h'
include 'pwhg_rad.h'
include 'pwhg_flg.h'
real * 8 powheginput
external powheginput
integer i1,n1,n2
call init_flsttag
flg_debug=.false.
c default for scheme in Virtual
c flg_drscheme=.false. Traditional MSbar
c flg_drscheme=.true. Dimensional reduction
flg_drscheme=.false.
c On if there are electromagnetic corrections
flg_with_em=.false.
c default is 6; if em is included, it is 22 to include photons
pdf_nparton=6
if(powheginput("#flg_debug").eq.1) flg_debug=.true.
c whether to output negative weights or not
flg_withnegweights=.true.
if(powheginput("#withnegweights").eq.0)flg_withnegweights=.false.
c See if we have weighted events
flg_weightedev=.false.
if(powheginput("#bornsuppfact").gt.0) flg_weightedev=.true.
if(powheginput("#ptsupp").gt.0) then
write(*,*) ' ******** WARNING: ptsupp is deprecated'
write(*,*) ' ******** Replace it with bornsuppfact'
call flush(6)
call exit(-1)
endif
c if true, suppress emitted particle energy with respect to the emitter
c in all cases (i.e. not only for g g) in FSR. Thus, consider also
c the emitter-emitted pairs qbar q, g q.
flg_doublefsr=.false.
if(powheginput("#doublefsr").gt.0) flg_doublefsr=.true.
c Set to true to remember and use identical values of the computed
c amplitudes, for Born, real and virtual contributions
flg_smartsig=.true.
if(powheginput("#smartsig").eq.0) flg_smartsig=.false.
c if true also counterterm are output in NLO tests (default)
flg_withsubtr=.true.
if(powheginput("#withsubtr").eq.0) flg_withsubtr=.false.
c The following turns on mechanism to deal with Born zeroes
flg_bornzerodamp=.false.
c If the damp function has to be called, the following must be on
flg_withdamp=.false.
c Traditional mechanism: either withdamp 1 or hfact >0 turn on
c remnant calculation. This mechanism is prone to error, since one may
c want hfact without bornzerodamp. We keep it for backward compatibility.
c In the future, we should use the flags listed below
if(powheginput("#hfact").gt.0d0) then
flg_withdamp=.true.
flg_bornzerodamp=.true.
endif
if(powheginput("#withdamp").eq.1) then
flg_withdamp=.true.
flg_bornzerodamp=.true.
endif
c Modern mechanism;
if(powheginput("#bornzerodamp").eq.1) then
flg_withdamp=.true.
flg_bornzerodamp=.true.
endif
c hdamp will replace the old hfact in new implementations
if(powheginput("#hdamp").gt.0) then
flg_withdamp=.true.
endif
c If set do only the Born term
flg_bornonly=.false.
if (powheginput("#bornonly").eq.1) flg_bornonly=.true.
flg_LOevents=.false.
if (powheginput("#LOevents").eq.1) then
flg_LOevents=.true.
flg_bornonly=.true.
endif
flg_novirtual=.false.
if (powheginput("#novirtual").eq.1) then
flg_novirtual=.true.
endif
c initialize random number sequence
i1=powheginput('#iseed')
if (i1.lt.0) i1=0
n1=powheginput('#rand1')
if (n1.lt.0) n1=0
n2=powheginput('#rand2')
if (n2.lt.0) n2=0
c select which upper bounding function form
rad_iupperisr=powheginput("#iupperisr")
if(rad_iupperisr.lt.0) rad_iupperisr=1
rad_iupperfsr=powheginput("#iupperfsr")
if(rad_iupperfsr.lt.0) rad_iupperfsr=2
c info on pdf for each event
flg_pdfreweight=.false.
if(powheginput("#pdfreweight").eq.1)flg_pdfreweight=.true.
c infos to write LH file with infos to reweight events
flg_reweight=.not.(powheginput('#storeinfo_rwgt').eq.0)
c infos to read LH file with infos to reweight events,
c and produce a new LH file with updated weights
flg_newweight=powheginput('#compute_rwgt').eq.1
c Is MiNLO required?
flg_minlo=.false.
st_bornorder=0
if(powheginput("#minlo").eq.1) then
flg_minlo=.true.
c set st_bornorder to an impossible value; later on we check
c if the user has set it. This is only used by minlo
st_bornorder=-1000
endif
c whether to correct for upper bound violations in the generation
c of btilde and remnant events
flg_ubexcess_correct = powheginput("#ubexcess_correct") ==1
call setrandom(i1,n1,n2)
c assign a default id for the process at hand
c if the user want to assign different id's
c to each subprocess, he/she can reassign lprup(1)
c inside the user-defined subroutine init_processes
lprup(1)=10001
c initialize physical parameters
call init_phys
if(flg_minlo.and.st_bornorder.eq.-1000) then
write(*,*) '************************************'
write(*,*) ' ERROR!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*) ' st_bornorder should be set by the user'
write(*,*) ' it should equal the powers of alpha_S'
write(*,*) ' present at the Born level'
write(*,*) ' Exiting ...'
call pwhg_exit(-1)
endif
c ID of beam particles 1 and 2 (pdg convention)
c proton:
ebmup(1)=kn_beams(0,1)
ebmup(2)=kn_beams(0,2)
if(abs(pdf_ih1).eq.1) then
c proton and antiproton
idbmup(1) = 2212*pdf_ih1
elseif(abs(pdf_ih1).eq.2) then
c neutron and antineutron
idbmup(1) = 2112*pdf_ih1/abs(pdf_ih1)
else
c can be used to pass directly the pdg id of the projectilr
idbmup(1)=pdf_ih1
endif
if(abs(pdf_ih2).eq.1) then
idbmup(2) = 2212*pdf_ih2
elseif(abs(pdf_ih2).eq.2) then
idbmup(2) = 2112*pdf_ih2/abs(pdf_ih2)
else
c can be used to pass directly the pdg id of the projectilr
idbmup(2)=pdf_ih2
endif
c pdf group; negative to use internal herwig pdf's for showering
pdfgup(1)=-1
pdfgup(2)=-1
c pdf set
pdfsup(1)=-1
pdfsup(2)=-1
c If either weighted events or event with negative weights are
c required, use idwtup=-4. Thuse, in both these cases, the average of the
c event weight is the total cross section.
c Otherwise the weight is set to 1 for each event.
c User processes may override these choices. In particular, using -4 in
c all cases is recommended for new processes. Using 3 is left here for
c compatibility with older implementations.
if(flg_withnegweights.or.flg_weightedev) then
idwtup = -4
else
idwtup = 3
endif
c number of user subprocesses
c Irrelevant if idwtup=+-3,+-4
nprup = 1
call bbinit
c now the cross section is available
if(flg_weightedev) then
xsecup(1)=-1
xerrup(1)=-1
else
xsecup(1)=rad_tot *rad_branching
xerrup(1)=rad_etot *rad_branching
endif
xmaxup(1)=1
end
subroutine init_flsttag
include 'nlegborn.h'
include 'pwhg_flst.h'
integer j,l
do j=1,maxprocreal
do l=1,nlegreal
flst_realtags(l,j)=0
flst_realres(l,j)=0
enddo
enddo
do j=1,maxprocborn
do l=1,nlegborn
flst_borntags(l,j)=0
flst_bornres(l,j)=0
enddo
enddo
end