/
lepto.F
337 lines (308 loc) · 9.96 KB
/
lepto.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
C **********************************************************************
SUBROUTINE LEPTO
IMPLICIT NONE
C...Administer the generation of an event.
C...Note: if error flag LST(21) is non-zero, no proper event generated.
COMMON /LINTRL/ PSAVE(3,4,5),KSAVE(4),XMIN,XMAX,YMIN,YMAX,
&Q2MIN,Q2MAX,W2MIN,W2MAX,ILEP,INU,IG,IZ
REAL PSAVE,XMIN,XMAX,YMIN,YMAX,Q2MIN,Q2MAX,W2MIN,W2MAX
INTEGER KSAVE,ILEP,INU,IG,IZ
SAVE /LINTRL/
*
* to avoid variable conflictions, a second keep element is necessary
* with the same common block name (see LPTOU2)
*
COMMON /LEPTOU/ CUT(14),LST(40),PARL(30),
& X,Y,W2,Q2,U
REAL CUT,PARL,X,Y,W2,Q2,U
INTEGER LST
SAVE /LEPTOU/
COMMON /LINTER/ PARI(50),EWQC(2,2,8),QC(8),ZL(2,4),ZQ(2,8),PQ(17)
REAL PARI,EWQC,QC,ZL,ZQ,PQ
SAVE /LINTER/
INTEGER NLUPDM,NPLBUF
PARAMETER (NLUPDM=4000,NPLBUF=5)
COMMON/LUJETS/N,K(NLUPDM,5),P(NLUPDM,NPLBUF),V(NLUPDM,5)
INTEGER N,K
REAL P,V
SAVE /LUJETS/
COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
INTEGER MSTU,MSTJ
REAL PARU,PARJ
SAVE /LUDAT1/
COMMON /LBOOST/ DBETA(2,3),STHETA(2),SPHI(2),PB(5),PHIR
DOUBLE PRECISION DBETA
REAL STHETA,SPHI,PB,PHIR
SAVE /LBOOST/
* if ARIADNE is used to simulate the parton shower evolution, the
* ARDAT1 Common block is neccessary for a proper interface.
* if ARIADNE is used to simulate the parton shower evolution, the
* ARDAT1 Common block is neccessary for a proper interface.
COMMON /ARDAT1/ PARA(40),MSTA(40)
REAL PARA
INTEGER MSTA
SAVE /ARDAT1/
INTEGER NUMMIS,NWARN,I,J,NS,L17
REAL ULMASS,ULANGL,RLU,ULALPS
REAL QG,QQB,SPQ,SRLU,PLU,PT,PHI,ENERGY,P2,ATAN
DOUBLE PRECISION DETOT,DARI29,DARI30
DIMENSION SPQ(17)
DATA NUMMIS,NWARN/0,10/,DARI29,DARI30/2*0.D0/
L17=0
1 LST(21)=0
DO 10 I=1,10
DO 10 J=1,5
K(I,J)=0
10 V(I,J)=0.
DO 15 I=1,4
K(I,1)=21
15 K(I,2)=KSAVE(I)
K(4,1)=1
N=2
IF(LST(17).NE.0.AND.LST(2).GT.0) THEN
C...Lepton and/or nucleon energy may vary from event to event,
IF(L17.EQ.0) THEN
C...Momentum vectors from P(i,j) i=1,2 j=1,2,3 on entry in LEPTO
DO 20 I=1,2
P(I,5)=ULMASS(K(I,2))
P(I,4)=SQRT(P(I,1)**2+P(I,2)**2+P(I,3)**2+P(I,5)**2)
DO 20 J=1,5
20 PSAVE(3,I,J)=P(I,J)
ELSE
C...Momentum vectors from PSAVE if new try, i.e. jump back to 1
DO 25 I=1,2
DO 25 J=1,5
25 P(I,J)=PSAVE(3,I,J)
ENDIF
L17=1
C...Transform to cms of incoming particles, lepton along +z axis.
DO 30 J=1,3
30 DBETA(1,J)=(DBLE(P(1,J))+DBLE(P(2,J)))/
& (DBLE(P(1,4))+DBLE(P(2,4)))
CALL LUDBRB(0,0,0.,0.,-DBETA(1,1),-DBETA(1,2),-DBETA(1,3))
SPHI(1)=ULANGL(P(1,1),P(1,2))
CALL LUDBRB(0,0,0.,-SPHI(1),0.D0,0.D0,0.D0)
STHETA(1)=ULANGL(P(1,3),P(1,1))
CALL LUDBRB(0,0,-STHETA(1),0.,0.D0,0.D0,0.D0)
LST(28)=2
PARL(21)=2.*(P(1,4)*P(2,4)-P(1,3)*P(2,3))
ELSE
C...Initial state momenta fixed from LINIT call.
DO 42 I=1,2
DO 40 J=1,5
40 P(I,J)=PSAVE(3,I,J)
42 IF(PSAVE(3,1,3).LT.0.) P(I,3)=-PSAVE(3,I,3)
LST(28)=3
ENDIF
CALL LEPTOX
C...Return if error or if no event to be generated.
IF(LST(21).NE.0.OR.LST(2).LE.0.OR.LST(7).EQ.-1) RETURN
IF(PARI(29).LT.0.5) THEN
C...For first call, reset double precision counters.
DARI29=0.D0
DARI30=0.D0
ENDIF
DARI29=DARI29+1.D0
PARI(29)=DARI29
C CALL GULIST(-3,2)
C...Scattered lepton and exchanged boson added to event record in LKINEM
C...Transform to lepton-nucleon cms if not made earlier
IF(LST(17).EQ.0) THEN
DO 46 I=3,4
DO 45 J=1,5
45 PSAVE(3,I,J)=P(I,J)
46 IF(PSAVE(3,1,3).LT.0.) PSAVE(3,I,3)=-P(I,3)
CALL LUDBRB(0,0,0.,0.,0.D0,0.D0,-DBETA(1,3))
LST(28)=2
ENDIF
DO 50 I=1,4
DO 50 J=1,5
50 PSAVE(2,I,J)=P(I,J)
C CALL GULIST(-2,2)
C...Prepare for parton cascade.
IF(LST(8).GE.2.AND.MOD(LST(8),10).NE.9) CALL LSHOWR(0)
C...Transform to hadronic cms, boost parameters in double precision.
DETOT=DBLE(P(1,4))-DBLE(P(4,4))+DBLE(P(2,4))
DBETA(2,1)=-DBLE(P(4,1))/DETOT
DBETA(2,2)=-DBLE(P(4,2))/DETOT
DBETA(2,3)=(DBLE(P(1,3))-DBLE(P(4,3))+DBLE(P(2,3)))/DETOT
CALL LUDBRB(0,0,0.,0.,-DBETA(2,1),-DBETA(2,2),-DBETA(2,3))
SPHI(2)=0.
STHETA(2)=ULANGL(P(3,3),P(3,1))
CALL LUDBRB(0,0,-STHETA(2),0.,0.D0,0.D0,0.D0)
LST(28)=1
DO 60 I=1,4
DO 60 J=1,5
60 PSAVE(1,I,J)=P(I,J)
C...Save momentum of exchanged boson (used in subroutine LFRAME).
DO 70 J=1,5
70 PB(J)=P(3,J)
C CALL GULIST(-1,2)
90 N=4
MSTU(1)=N+1
LST(26)=N+1
LST(27)=0
PARL(25)=ULALPS(Q2)
IF(LST(8).EQ.1.OR.LST(8)/10.EQ.1.OR.MOD(LST(8),10).EQ.9) THEN
C...Probabilities for hard, first order QCD events.
CAE...Corrected what to do when LQGEV or LQQBEV fail. Now make LQEV.
CALL LQCDPR(QG,QQB)
DO 100 I=1,17
100 SPQ(I)=PQ(I)
200 SRLU=RLU(0)
IF(SRLU.GT.QQB+QG) THEN
CALL LQEV
ELSEIF(SRLU.GT.QQB) THEN
IF(LST(8).EQ.9) THEN
DO 211 I=1,17
211 PQ(I)=SPQ(I)
CALL LQEV
ELSE
CALL LQGEV
DO 212 I=1,17
PQ(I)=SPQ(I)
212 CONTINUE
ENDIF
ELSE
CALL LQQBEV
DO 213 I=1,17
PQ(I)=SPQ(I)
213 CONTINUE
IF(LST(8).EQ.9.AND.LST(21).EQ.0) THEN
IF(PLU(5,11).LT.Q2*PARA(20)) THEN
DO 220 I=1,17
220 PQ(I)=SPQ(I)
CALL LQEVAR(K(5,2),K(7,2))
ENDIF
ENDIF
ENDIF
IF(LST(21).NE.0) THEN
230 CALL LQEV
IF(LST(21).NE.0) GOTO 230
ENDIF
ELSE
C...QPM model without QCD corrections (cascade applied later).
300 CALL LQEV
IF(LST(21).NE.0) GOTO 300
ENDIF
NS=MSTU(1)
MSTU(1)=0
C CALL GULIST(-3,2)
C WRITE(6,*) ' LST(24)=',LST(24)
CJR-- no preclustering of small systems
MSTJ(14)=-1
CJR--
IF(LST(8).LE.1.OR.MOD(LST(8),10).EQ.9) THEN
C...No parton cascade, introduce primordial kt.
IF(PARL(3).GT.1.E-03) THEN
CALL LPRIKT(PARL(3),PT,PHI)
CALL LUDBRB(NS,N,0.,-PHI,0.D0,0.D0,0.D0)
CALL LUDBRB(NS,N,ATAN(2.*PT/SQRT(W2)),PHI,0.D0,0.D0,0.D0)
ENDIF
IF(MOD(LST(8),10).NE.9) THEN
C...Check system against fragmentation cuts.
MSTU(24)=0
CALL LUPREP(0)
IF(MSTU(24).NE.0) THEN
IF(LST(3).GE.1) WRITE(6,*)'LUPREP error MSTU(24)=',MSTU(24),
& ', New event generated'
LST(21)=1
GOTO 1
ENDIF
ENDIF
ELSEIF(LST(24).EQ.1) THEN
C...Include parton cascades (+ remnant & kt) on q-event
CALL LSHOWR(1)
ELSE
C...Include parton cascades (+ remnant & kt) on qg- or qqbar-event
CALL LMEPS
ENDIF
IF(LST(21).NE.0) THEN
C IF(LST(3).GE.1)
C & WRITE(6,*)'Cascade error LST(21)= ',LST(21),
C & ', New event generated'
GOTO 1
ENDIF
CJR-- Soft colour interactions
IF(LST(34).EQ.1) CALL LSCI(PARL(7))
IF(LST(21).NE.0) GOTO 1
CJR-- take care of small systems
CALL LSMALL
IF(LST(21).NE.0) THEN
IF(LST(3).GE.1) WRITE(6,*)' LSMALL error LST(21)= ',LST(21),
& ', New event generated'
GOTO 1
ENDIF
MSTJ(14)=1
CALL LUPREP(0)
IF(MSTU(24).NE.0) THEN
IF(LST(3).GE.1) WRITE(6,*)' LUPREP error MSTU(24)= ',MSTU(24),
& ', New event generated'
LST(21)=1
ENDIF
CJR--
IF(LST(21).NE.0) GOTO 1
DO 400 I=1,N
C...Correct energy-momentum-mass mismatch for real particle
IF(P(I,5).LT.0.) GOTO 400
ENERGY=SQRT(DBLE(P(I,5))**2+DBLE(P(I,1))**2+DBLE(P(I,2))**2+
&DBLE(P(I,3))**2)
P2=DBLE(P(I,4))**2-DBLE(P(I,1))**2-DBLE(P(I,2))**2-DBLE(P(I,3))**2
IF(ABS(ENERGY-P(I,4))/(PSAVE(3,1,4)+PSAVE(3,2,4)).GT.PARU(11))THEN
NUMMIS=NUMMIS+1
C...For testing purposes
C IF(LST(3).GE.1.AND.NUMMIS.LE.NWARN) THEN
C WRITE(6,1000) I,(K(I,J),J=1,2),(P(I,J),J=1,5),
C & SIGN(SQRT(ABS(P2)),P2),ENERGY,INT(DARI29),NWARN
C IF(ABS(P2-P(I,5)**2).GT.400.) CALL LULIST(2)
C ENDIF
CAE WRITE(6,*) 'Energy mismatch',LST(24),PARL(28),PARL(29),NUMMIS
GOTO 90
ENDIF
P(I,4)=ENERGY
400 CONTINUE
DARI30=DARI30+1.D0
PARI(30)=DARI30
Ctest IF(LST(23).EQ.2) PARL(24)=PARL(24)*DARI30/DARI29
DO 500 I=1,N
DO 500 J=1,5
500 V(I,J)=0.
IF(LST(7).EQ.1) THEN
CALL LUEXEC
IF(MSTU(24).NE.0) THEN
WRITE(6,*) ' Error from JETSET, new event made'
GOTO 90
ENDIF
ENDIF
C CALL GULIST(-1,2)
C...Transform to desired frame
C LST(28)=1
LST(29)=0
PHIR=6.2832*RLU(0)
IF(LST(17).EQ.0) THEN
IF(LST(5).GE.2) CALL LFRAME(LST(5),0)
C...Restore momenta (e,p,boson,L) due to numerical errors from boosts
C...But only in frames 1-3
IF(LST(28).GT.0 .AND. LST(28).LT.4) THEN
DO 600 I=1,4
DO 600 J=1,5
600 P(I,J)=PSAVE(LST(28),I,J)
ENDIF
IF(LST(6).EQ.1.AND.LST(28).GE.2) THEN
C...Random rotation in azimuthal angle
CALL LUDBRB(0,0,0.,PHIR,0.D0,0.D0,0.D0)
LST(29)=1
ENDIF
ELSE
IF(LST(5).GE.2) CALL LFRAME(LST(5),LST(6))
ENDIF
C...Deactivate scattered lepton
IF(MOD(LST(4),10).EQ.0) K(4,1)=21
C CALL GULIST(0,2)
RETURN
1000 FORMAT(' Warning: too large numerical mismatch in ',
&'particle energy-momentum-mass',
&/,3X,'I K(I,1) ..2) P(I,1) P(I,2) P(I,3)',
&' P(I,4) P(I,5) mass energy',/,I4,2I6,7F8.3,/,
&' Event no.',I8,' regenerated. Only first',I5,' warnings printed')
END