]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ISAJET/code/twojet.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Coding conventions.
[u/mrichter/AliRoot.git] / ISAJET / code / twojet.F
CommitLineData
0795afa3 1#include "isajet/pilot.h"
2 SUBROUTINE TWOJET
3C
4C Driving routine to generate initial parameters for jets,
5C assuming zero initial transverse momentum, ie PT(1)=PT(2).
6C
7C Parameters are PT,YJ,PHI with P,YJ,XJ as dependent variables,
8C where YJ=RAPIDITY, XJ=Feynman X.
9C All parameters are stored in COMMON/JETPAR/.
10C Cross section is called from NOGOOD.
11C
12#if defined(CERNLIB_IMPNONE)
13 IMPLICIT NONE
14#endif
15#include "isajet/idrun.inc"
16#include "isajet/itapes.inc"
17#include "isajet/keys.inc"
18#include "isajet/mbpar.inc"
19#include "isajet/pjets.inc"
20#include "isajet/pinits.inc"
21#include "isajet/jetlim.inc"
22#include "isajet/ptpar.inc"
23#include "isajet/jetpar.inc"
24#include "isajet/primar.inc"
25#include "isajet/partcl.inc"
26#include "isajet/const.inc"
27#include "isajet/jetsig.inc"
28#include "isajet/totals.inc"
29#include "isajet/isloop.inc"
30#include "isajet/sstype.inc"
31#include "isajet/xmssm.inc"
32C
33 REAL ACOSH,XXX,WTFCN,PPP,RANF,SIGN,SGN,AMQ1,AMASS,AMQ2
34 REAL PPLUS,PMINUS,PSUM3,PSUM4,PPL,PMN,SQ1,SQ2,ROOT,P1PL,P1MN
35 REAL P2PL,P2MN,AMI1,AMI2
36 INTEGER NREJ,I,II,IS,IFL1,IFL2
37 REAL X(2)
38 EQUIVALENCE (X(1),X1)
39 LOGICAL NOGOOD
40 LOGICAL YGENJ
41 INTEGER LISTJ(17),LISTW(4),LISTSS(85),LISTSM(30)
42C
43C SUSY IDENT codes from /SSTYPE/. (Fortran 77 allows - signs
44C in parameter statements but not data statements.)
45 INTEGER MSUPL,MSDNL,MSSTL,MSCHL,MSBT1,MSTP1,
46 $MSUPR,MSDNR,MSSTR,MSCHR,MSBT2,MSTP2,MSW1,MSW2,
47 $MSNEL,MSEL,MSNML,MSMUL,MSNTL,MSTAU1,MSER,MSMUR,MSTAU2
48 PARAMETER (MSUPL=-ISUPL)
49 PARAMETER (MSDNL=-ISDNL)
50 PARAMETER (MSSTL=-ISSTL)
51 PARAMETER (MSCHL=-ISCHL)
52 PARAMETER (MSBT1=-ISBT1)
53 PARAMETER (MSTP1=-ISTP1)
54 PARAMETER (MSUPR=-ISUPR)
55 PARAMETER (MSDNR=-ISDNR)
56 PARAMETER (MSSTR=-ISSTR)
57 PARAMETER (MSCHR=-ISCHR)
58 PARAMETER (MSBT2=-ISBT2)
59 PARAMETER (MSTP2=-ISTP2)
60 PARAMETER (MSW1=-ISW1)
61 PARAMETER (MSW2=-ISW2)
62 PARAMETER (MSNEL=-ISNEL)
63 PARAMETER (MSEL=-ISEL)
64 PARAMETER (MSNML=-ISNML)
65 PARAMETER (MSMUL=-ISMUL)
66 PARAMETER (MSNTL=-ISNTL)
67 PARAMETER (MSTAU1=-ISTAU1)
68 PARAMETER (MSER=-ISER)
69 PARAMETER (MSMUR=-ISMUR)
70 PARAMETER (MSTAU2=-ISTAU2)
71C
72 DATA LISTSS/ISGL,
73 $ISUPL,MSUPL,ISDNL,MSDNL,ISSTL,MSSTL,ISCHL,MSCHL,ISBT1,MSBT1,
74 $ISTP1,MSTP1,
75 $ISUPR,MSUPR,ISDNR,MSDNR,ISSTR,MSSTR,ISCHR,MSCHR,ISBT2,MSBT2,
76 $ISTP2,MSTP2,
77 $ISW1,MSW1,ISW2,MSW2,ISZ1,ISZ2,ISZ3,ISZ4,
78 $ISNEL,MSNEL,ISEL,MSEL,ISNML,MSNML,ISMUL,MSMUL,ISNTL,MSNTL,
79 $ISTAU1,MSTAU1,ISER,MSER,ISMUR,MSMUR,ISTAU2,MSTAU2,
80 $9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,11,-11,12,-12,13,-13,
81 $14,-14,15,-15,16,-16,10,80,-80,90,82,83,84,86,-86/
82 DATA LISTSM/9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,11,-11,12,-12,13,-13,
83 $14,-14,15,-15,16,-16,10,80,-80,90,81/
84 DATA LISTJ/9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,7,-7,8,-8/
85 DATA LISTW/10,80,-80,90/
86C Inverse hyperbolic cosine function
87 ACOSH(XXX)=ALOG(XXX+SQRT(XXX**2-1.))
88 WTFCN(PPP)=2.*PPP*PTGEN2*PTGEN3*PPP**((PTGEN3-1.)/PTGEN3)
89C
90C Initialize
91C
92 NPTCL=0
93 PHI(1)=PHIMIN(1)+(PHIMAX(1)-PHIMIN(1))*RANF()
94 PHI(2)=AMOD(PHI(1)+PI,2.*PI)
95 NREJ=-1
96 SIGMA=0.
97 WT=1.
98 IF(.NOT.FIXPT(2)) GOTO 101
99 FIXPT(1)=.TRUE.
100 PT(1)=PT(2)
101 101 CONTINUE
102 IF(FIXPT(1)) GOTO 400
103 DO 110 I=1,2
104 IF(FIXP(I)) GOTO 200
105 IF(FIXXJ(I)) GOTO 300
106 110 CONTINUE
107C
108C Genetate PT and YJ with no variables fixed
109C
110 111 NREJ=NREJ+1
111 IF(NREJ.GT.NTRIES) GO TO 910
112 SUMWT=SUMWT+SIGMA*WT/(NEVOLV*NFRGMN)
113 NKINPT=NKINPT+1
114 SIGMA=0.
115 WT=1.
116C Generate PT with a power law distribution
117 PT(1)=(PTGEN1+PTGEN2*RANF())**PTGEN3
118 PT(2)=PT(1)
119 SIGMAX=PTFUN1*PT(1)**PTFUN2
120C GENERATE FLAT IN YJ, CALCULATE CORRESPONDING TH
121 DO 115 I=1,2
122 IF(FIXYJ(I)) GOTO 115
123 IF(.NOT.YGENJ(I)) GOTO 111
124 115 CONTINUE
125 DO 116 I=1,2
126 P(I)=PT(I)/STH(I)
127 IF(P(I).LT.PMIN(I).OR.P(I).GT.PMAX(I)) GOTO 111
128 XJ(I)=P(I)*CTH(I)/HALFE
129 IF(XJ(I).LT.XJMIN(I).OR.XJ(I).GT.XJMAX(I)) GOTO 111
130 116 CONTINUE
131 WT=WT*WTFCN(PT(1))
132 IF(NOGOOD(1)) GOTO 111
133 SUMWT=SUMWT+SIGMA*WT/(NEVOLV*NFRGMN)
134 NKEEP=NKEEP+1
135 GO TO 500
136C
137C Generate PT and YJ fixing P
138C
139 200 CONTINUE
140 II=3-I
141 211 NREJ=NREJ+1
142 IF(NREJ.GT.NTRIES) GO TO 910
143 NKINPT=NKINPT+1
144 WT=0.
145 IF(FIXYJ(I)) GOTO 212
146C Generate PT with a power law distribution
147 PT(1)=(PTGEN1+PTGEN2*RANF())**PTGEN3
148 SIGMAX=PTFUN1*PT(1)**PTFUN2
149 PT(2)=PT(1)
150C Given PT, TH is fixed except for a sign
151 STH(I)=PT(I)/P(I)
152 SIGN=1.0
153 IF(RANF().GT.0.5) SIGN=-1.0
154 CTH(I)=SIGN*SQRT(1.-STH(I)**2)
155 TH(I)=ATAN2(STH(I),CTH(I))
156 YJ(I)=-ALOG(TAN(TH(I)/2.))
157 IF(YJ(I).LT.YJMIN(I).OR.YJ(I).GT.YJMAX(I)) GOTO 211
158 GOTO 213
159 212 PT(1)=P(I)*STH(I)
160 213 CONTINUE
161 XJ(I)=P(I)*CTH(I)/HALFE
162 IF(XJ(I).LT.XJMIN(I).OR.XJ(I).GT.XJMAX(I)) GOTO 211
163 IF(FIXP(II)) GOTO 220
164 IF(FIXXJ(II)) GOTO 230
165 IF(FIXYJ(II)) GOTO 215
166 IF(.NOT.YGENJ(II)) GOTO 211
167 215 CONTINUE
168 P(II)=PT(II)/STH(II)
169 IF(P(II).LT.PMIN(II).OR.P(II).GT.PMAX(II)) GOTO 211
170 XJ(II)=P(II)*CTH(II)/HALFE
171 IF(XJ(II).LT.XJMIN(II).OR.XJ(II).GT.XJMAX(II)) GOTO 211
172 GOTO 250
173220 STH(II)=PT(II)/P(II)
174 SGN=1.0
175 IF(RANF().GT.0.5) SGN=-1.0
176 CTH(II)=SGN*SQRT(1.-STH(II)**2)
177 TH(II)=ATAN2(STH(II),CTH(II))
178 YJ(II)=-ALOG(TAN(TH(II)/2.))
179 IF(YJ(II).LT.YJMIN(II).OR.YJ(II).GT.YJMAX(II)) GOTO 211
180 XJ(II)=P(II)*CTH(II)/HALFE
181 IF(XJ(II).LT.XJMIN(II).OR.XJ(II).GT.XJMAX(II)) GOTO 211
182 GOTO 250
183 230 TH(II)=ATAN2(PT(II),XJ(II)*HALFE)
184 YJ(II)=-ALOG(TAN(TH(II)/2.))
185 IF(YJ(II).LT.YJMIN(II).OR.YJ(II).GT.YJMAX(II)) GOTO 211
186 CTH(II)=COS(TH(II))
187 STH(II)=SIN(TH(II))
188 250 CONTINUE
189 IF(NOGOOD(1)) GOTO 211
190 NKEEP=NKEEP+1
191 GO TO 500
192C
193C Generate PT and YJ at fixed XJ
194C
195 300 CONTINUE
196 II=3-I
197 311 NREJ=NREJ+1
198 IF(NREJ.GT.NTRIES) GO TO 910
199 NKINPT=NKINPT+1
200 WT=0.
201C Generate PT with a power law distribution
202 PT(1)=(PTGEN1+PTGEN2*RANF())**PTGEN3
203 SIGMAX=PTFUN1*PT(1)**PTFUN2
204 PT(2)=PT(1)
205 TH(I)=ATAN2(PT(I),XJ(I)*HALFE)
206 YJ(I)=-ALOG(TAN(TH(I)/2.))
207 IF(YJ(I).LT.YJMIN(I).OR.YJ(I).GT.YJMAX(I)) GOTO 311
208 CTH(I)=COS(TH(I))
209 STH(I)=SIN(TH(I))
210 P(I)=PT(I)/STH(I)
211 IF(FIXYJ(II)) GOTO 315
212 IF(FIXP(II)) GOTO 314
213 YJ(II)=YJMIN(II)+(YJMAX(II)-YJMIN(II))*RANF()
214 TH(II)=2.*ATAN(EXP(-YJ(II)))
215 CTH(II)=COS(TH(II))
216 STH(II)=SIN(TH(II))
217 GOTO 315
218 314 CONTINUE
219 STH(II)=PT(II)/P(II)
220 CTH(II)=SQRT(1.-STH(II)**2)
221 IF(RANF().GT.0.5) CTH(II)=-CTH(II)
222 TH(II)=ATAN2(STH(II),CTH(II))
223 YJ(II)=-ALOG(TAN(TH(II)/2.))
224 315 CONTINUE
225 P(II)=PT(II)/STH(II)
226 XJ(II)=P(II)*CTH(II)/HALFE
227 IF(XJ(II).LT.XJMIN(II).OR.XJ(II).GT.XJMAX(II)) GOTO 311
228 IF(NOGOOD(1)) GOTO 311
229 NKEEP=NKEEP+1
230 GO TO 500
231C
232C Generate YJ at fixed PT
233C
234 400 CONTINUE
235 PT(2)=PT(1)
236 411 NREJ=NREJ+1
237 IF(NREJ.GT.NTRIES) GO TO 910
238 NKINPT=NKINPT+1
239 WT=0.
240 DO 415 I=1,2
241 IF(FIXYJ(I)) GOTO 415
242 IF(FIXP(I)) GOTO 413
243 IF(.NOT.YGENJ(I)) GO TO 411
244 GOTO 414
245 413 CONTINUE
246 IS=1
247 IF(RANF().GT.0.5) IS=2
248 CTH(I)=CTHS(IS,I)
249 TH(I)=THS(IS,I)
250 YJ(I)=YJS(IS,I)
251 414 CONTINUE
252 P(I)=PT(I)/STH(I)
253 XJ(I)=P(I)*CTH(I)/HALFE
254 415 CONTINUE
255 IF(NOGOOD(1)) GOTO 411
256 NKEEP=NKEEP+1
257C
258C Reset /JETPAR/
259C
260 500 CONTINUE
261 IF(KEYS(1)) THEN
262 IFL1=LISTJ(JETTYP(1))
263 IFL2=LISTJ(JETTYP(2))
264 AMQ1=AMASS(IFL1)
265 AMQ2=AMASS(IFL2)
266 AMI1=AMASS(LISTJ(INITYP(1)))
267 AMI2=AMASS(LISTJ(INITYP(2)))
268 CALL TWOKIN(AMI1,AMI2,AMQ1,AMQ2)
269 ELSEIF(KEYS(5).OR.(KEYS(10).AND.GOMSSM)) THEN
270 IFL1=LISTSS(JETTYP(1))
271 IFL2=LISTSS(JETTYP(2))
272 AMQ1=AMASS(IFL1)
273 AMQ2=AMASS(IFL2)
274 CALL TWOKIN(0.,0.,AMQ1,AMQ2)
275 ELSEIF(KEYS(6)) THEN
276 IFL1=LISTW(JETTYP(1))
277 IFL2=LISTW(JETTYP(2))
278 AMQ1=AMASS(IFL1)
279 AMQ2=AMASS(IFL2)
280 CALL TWOKIN(0.,0.,AMQ1,AMQ2)
281 ELSEIF(KEYS(8)) THEN
282 IF(JETTYP(1).LE.13) THEN
283 IFL1=LISTJ(JETTYP(1))
284 ELSE
285 IFL1=10
286 ENDIF
287 IF(JETTYP(2).LE.13) THEN
288 IFL2=LISTJ(JETTYP(2))
289 ELSE
290 IFL2=10
291 ENDIF
292 AMQ1=AMASS(IFL1)
293 AMQ2=AMASS(IFL2)
294 CALL TWOKIN(0.,0.,AMQ1,AMQ2)
295 ELSEIF(KEYS(10).AND.(.NOT.GOMSSM)) THEN
296 IFL1=LISTSM(JETTYP(1))
297 IFL2=LISTSM(JETTYP(2))
298 AMQ1=AMASS(IFL1)
299 AMQ2=AMASS(IFL2)
300 CALL TWOKIN(0.,0.,AMQ1,AMQ2)
301 ENDIF
302C
303C Set PBEAM and PJETS
304C
305 PBEAM(1)=(1.-X1)*HALFE
306 PBEAM(2)=(1.-X2)*HALFE
307 DO 501 I=1,2
308 PJETS(3,I)=P(I)*CTH(I)
309 PJETS(1,I)=PT(I)*COS(PHI(I))
310 PJETS(2,I)=PT(I)*SIN(PHI(I))
311 IF(KEYS(1)) THEN
312 IDJETS(I)=LISTJ(JETTYP(I))
313 ELSEIF(KEYS(5).OR.(KEYS(10).AND.GOMSSM)) THEN
314 IDJETS(I)=LISTSS(JETTYP(I))
315 ELSEIF(KEYS(6)) THEN
316 IDJETS(I)=LISTW(JETTYP(I))
317 ELSEIF(KEYS(8)) THEN
318 IDJETS(1)=IFL1
319 IDJETS(2)=IFL2
320 ELSEIF(KEYS(10)) THEN
321 IDJETS(I)=LISTSM(JETTYP(I))
322 ENDIF
323 PJETS(5,I)=AMASS(IDJETS(I))
324 PJETS(4,I)=SQRT(P(I)**2+PJETS(5,I)**2)
325 501 CONTINUE
326C
327C Set PINITS
328C
329 DO 600 I=1,2
330 IDINIT(I)=LISTJ(INITYP(I))
331 PINITS(5,I)=AMASS(IDINIT(I))
332 PPLUS=X(I)*ECM
333 PMINUS=PINITS(5,I)**2/PPLUS
334 PINITS(4,I)=.5*(PPLUS+PMINUS)
335 PINITS(3,I)=.5*(PPLUS-PMINUS)*(3-2*I)
336 PINITS(2,I)=0.
337 PINITS(1,I)=0.
338600 CONTINUE
339C Calculate PINITS exactly.
340 PSUM3=PJETS(3,1)+PJETS(3,2)
341 PSUM4=PJETS(4,1)+PJETS(4,2)
342 IF(PSUM3.GT.0.) THEN
343 PPL=PSUM4+PSUM3
344 PMN=SHAT/PPL
345 ELSE
346 PMN=PSUM4-PSUM3
347 PPL=SHAT/PMN
348 ENDIF
349 SQ1=PINITS(5,1)**2
350 SQ2=PINITS(5,2)**2
351 ROOT=SQRT((PPL*PMN-SQ1-SQ2)**2-4.*SQ1*SQ2)
352 P1PL=(PPL*PMN+SQ1-SQ2+ROOT)/(2.*PMN)
353 P1MN=SQ1/P1PL
354 P2MN=(PPL*PMN+SQ2-SQ1+ROOT)/(2.*PPL)
355 P2PL=SQ2/P2MN
356 PINITS(4,1)=.5*(P1PL+P1MN)
357 PINITS(3,1)=.5*(P1PL-P1MN)
358 PINITS(4,2)=.5*(P2PL+P2MN)
359 PINITS(3,2)=.5*(P2PL-P2MN)
360 RETURN
361C
362C Error
363C
364910 CALL PRTEVT(0)
365 WRITE(ITLIS,1000) NREJ
366 1000 FORMAT(//' IT IS TAKING MORE THAN',I5,' TRIES TO GENERATE AN',
367 $' EVENT. CHECK LIMITS OR INCREASE NTRIES.')
368 STOP 99
369 END