More volume overlaps corrected
[u/mrichter/AliRoot.git] / ISAJET / code / frgmnt.F
1 #include "isajet/pilot.h"
2       SUBROUTINE FRGMNT
3 C
4 C          Control jet fragmentation.  Boost to frames defined in
5 C          EVOLVE and call JETGEN.
6 C
7 C          EVOLVE initializes /PJSET/ as follows--
8 C          1      - 2            = PINITS (except for E+E-)
9 C          N0W    - N0W          = QWJET  (for DRELLYAN, NJET=3)
10 C          N0JETS - N0JETS+NJET  = PJETS
11 C          N0PAIR - N0PAIR+NPAIR = PPAIR  (for WPAIR)
12 C
13 #if defined(CERNLIB_IMPNONE)
14       IMPLICIT NONE
15 #endif
16 #include "isajet/itapes.inc"
17 #include "isajet/primar.inc"
18 #include "isajet/jetpar.inc"
19 #include "isajet/pjets.inc"
20 #include "isajet/pinits.inc"
21 #include "isajet/partcl.inc"
22 #include "isajet/const.inc"
23 #include "isajet/jetset.inc"
24 #include "isajet/jwork.inc"
25 #include "isajet/keys.inc"
26 #include "isajet/q1q2.inc"
27 #include "isajet/frame.inc"
28 C
29       REAL PSUM(5),PALLJ(5),P12(5),PIN(5,2),PWREST(5),PADD(5)
30       REAL POLD(5),PNEW(5)
31       REAL PINPL,PINMN,BP,PT2AVE,PTADD,RANF,PHIADD,PALLPL,PALLMN
32       REAL PALLX,PALLY
33       INTEGER K,J,JJET,NZERO,IB,NPTCL1,NPTCL2,IFAIL,JET,NPJET1,NPLV1
34       INTEGER NPJET3,IP,NP1,NP2,NFIRST,IP1,IFR,NLJ
35       DOUBLE PRECISION DSUM(5),DPASS(5)
36 C
37 C          Initialize
38       DO 100 K=1,5
39 100   DSUM(K)=0.
40       NLJ=NJET
41       IF(KEYS(3)) NLJ=NJET+1
42       DO 101 J=1,NLJ
43         JJET=N0JETS+J-1
44         IF(JJET.EQ.N0W) GOTO 101
45         CALL DBLVEC(PJSET(1,JJET),DPASS)
46         DO 102 K=1,4
47 102     DSUM(K)=DSUM(K)+DPASS(K)
48 101   CONTINUE
49       DSUM(5)=DSQRT(DSUM(4)**2-DSUM(1)**2-DSUM(2)**2-DSUM(3)**2)
50       DO 103 K=1,5
51 103   PALLJ(K)=DSUM(K)
52 C
53       NZERO=NPTCL+1
54 C
55 C          Fragment partons from initial state shower
56 C
57       IF(.NOT.KEYS(2)) THEN
58         DO 110 J=1,NJSET
59           IF(JDCAY(J).EQ.JPACK*J+J) THEN
60             IB=JORIG(J)/JPACK-10
61             DO 120 K=1,5
62 120         PIN(K,IB)=PJSET(K,J)
63           ENDIF
64 110     CONTINUE
65 C
66         CALL FRGJET(11)
67         CALL FRGJET(12)
68 C
69         NPTCL1=NPTCL+1
70         NPTCL2=NPTCL1+1
71         IF(NPTCL1.GT.MXPTCL) GO TO 9999
72         PINPL=.5*(PIN(4,1)+PIN(3,1)+PIN(4,2)+PIN(3,2))
73         PINMN=.5*(PIN(4,1)-PIN(3,1)+PIN(4,2)-PIN(3,2))
74         PPTCL(1,NPTCL1)=0.
75         PPTCL(2,NPTCL1)=0.
76         PPTCL(3,NPTCL1)=HALFE-PINPL
77         PPTCL(4,NPTCL1)=HALFE-PINPL
78         PPTCL(5,NPTCL1)=0.
79         PPTCL(1,NPTCL2)=0.
80         PPTCL(2,NPTCL2)=0.
81         PPTCL(3,NPTCL2)=-(HALFE-PINMN)
82         PPTCL(4,NPTCL2)=HALFE-PINMN
83         PPTCL(5,NPTCL2)=0.
84         DO 130 K=1,4
85 130     PSUM(K)=-PALLJ(K)
86         PSUM(4)=PSUM(4)+ECM
87         PSUM(5)=PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2
88         IF(PSUM(5).GE.0.) THEN
89           PSUM(5)=SQRT(PSUM(5))
90           CALL RESCAL(NZERO,NPTCL2,PSUM,IFAIL)
91         ENDIF
92 C
93         DO 140 K=1,4
94 140     PBEAMS(K)=PPTCL(K,NPTCL1)+PPTCL(K,NPTCL2)
95         PBEAMS(5)=SQRT(PBEAMS(4)**2-PBEAMS(1)**2-PBEAMS(2)**2
96      $  -PBEAMS(3)**2)
97       ENDIF
98 C
99 C          Boost partons from final jets with -FRAME
100 C
101 200   DO 210 J=1,NJSET
102         JET=JORIG(J)/JPACK
103         IF ( JET.EQ.0 ) THEN
104           IFR=1
105         ELSE
106           IF(JET.GT.10) GO TO 210
107           IF(KEYS(6)) THEN
108             IF(IDJETS(JET).EQ.10) GO TO 210
109           ENDIF
110           IFR=IFRAME(JET)
111         ENDIF
112 C
113 C          Do this boost in double precision for 32-bit machines
114         CALL DBOOST(-1,FRAME(1,IFR),PJSET(1,J))
115 210   CONTINUE
116 C
117 C          Fragment partons from final jets
118 C
119       NPJET1=NPTCL+1
120       DO 220 K=1,4
121 220   PSUM(K)=0
122 C
123 C          Conserve mass of 1+2 for DRELLYAN (automatic for WPAIR)
124 C
125       IF(KEYS(3)) THEN
126         CALL FRGJET(1)
127         CALL FRGJET(2)
128         IF(STDDY) THEN
129           DO 230 K=1,4
130             PSUM(K)=PJSET(K,3)+PJSET(K,4)
131   230     CONTINUE
132         ELSE
133           DO 240 K=1,4
134             PSUM(K)=PJSET(K,N0W+1)+PJSET(K,N0W+2)
135   240     CONTINUE
136         ENDIF
137         PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2)
138         NPLV1=NPTCL
139         CALL RESCAL(NPJET1,NPLV1,PSUM,IFAIL)
140 C          EXTRADIM has only jet3 + graviton
141       ELSEIF(KEYS(11)) THEN
142         CALL FRGJET(3)
143         CALL FRGJET(0)
144         NPLV1=NPTCL
145         DO 241 K=1,4
146           PSUM(K)=PJSET(K,3)+PJSET(K,4)
147 241     CONTINUE
148         PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2)
149         CALL RESCAL (NPJET1,NPLV1,PSUM,IFAIL)
150       ELSE
151 C          All other processes
152         DO 242 J=1,NJET
153           JJET=N0JETS+J-1
154           CALL FRGJET(J)
155           DO 243 K=1,4
156   243     PSUM(K)=PSUM(K)+PJSET(K,JJET)
157   242   CONTINUE
158         PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2)
159         NPLV1=NPTCL
160         CALL RESCAL(NPJET1,NPLV1,PSUM,IFAIL)
161       ENDIF
162 C
163 C           Add extra jets for DRELLYAN
164       IF(KEYS(3).AND..NOT.STDDY) THEN
165         NPJET3=NPTCL+1
166         DO 245 J=3,NJET
167 245     CALL FRGJET(J)
168         NPTCL1=NPTCL+1
169         IF(NPTCL1.GT.MXPTCL) GO TO 9999
170         DO 250 K=1,4
171           PPTCL(K,NPTCL1)=PJSET(K,N0W)
172 250     PSUM(K)=PJSET(K,N0W)
173         PPTCL(5,NPTCL1)=PJSET(5,N0W)
174         DO 246 J=3,NJET
175           JJET=N0JETS+J-3
176           DO 246 K=1,4
177             PSUM(K)=PSUM(K)+PJSET(K,JJET)
178   246   CONTINUE    
179         PSUM(5)=SQRT(PSUM(4)**2-PSUM(1)**2-PSUM(2)**2-PSUM(3)**2)
180         CALL RESCAL(NPJET3,NPTCL1,PSUM,IFAIL)
181         DO 260 K=1,5
182 260     PWREST(K)=PPTCL(K,NPTCL1)
183       ENDIF
184 C
185 C          Boost partons back to PP COM
186 C
187       DO 300 J=1,NJSET
188         JET=JORIG(J)/JPACK
189         IF ( JET.EQ.0 ) THEN
190           IFR=1
191         ELSE
192           IF(JET.GT.10) GO TO 300
193           IF(KEYS(6)) THEN
194             IF(IDJETS(JET).EQ.10) GO TO 300
195           ENDIF
196           IFR=IFRAME(JET)
197         ENDIF
198         BP=0.
199         DO 305 K=1,3
200 305     BP=BP+FRAME(K,IFR)*PJSET(K,J)
201         BP=BP/FRAME(5,IFR)
202         DO 310 K=1,3
203 310     PJSET(K,J)=PJSET(K,J)+FRAME(K,IFR)*PJSET(4,J)/FRAME(5,IFR)
204      $  +FRAME(K,IFR)*BP/(FRAME(4,IFR)+FRAME(5,IFR))
205         PJSET(4,J)=FRAME(4,IFR)*PJSET(4,J)/FRAME(5,IFR)+BP
206 300   CONTINUE
207 C
208 C          Reset FRAME to boost hadrons to PP COM
209 C
210       IF(KEYS(1).OR.KEYS(2).OR.(KEYS(3).AND.NJET.EQ.2).OR.KEYS(5)
211      $.OR.(KEYS(7).AND.NPAIR.EQ.0).OR.KEYS(8)) THEN
212         DO 410 K=1,5
213           FRAME(K,1)=PALLJ(K)
214 410     CONTINUE
215       ELSEIF(KEYS(3).AND.NJET.GT.2) THEN
216         DO 420 K=1,5
217 420     FRAME(K,1)=PALLJ(K)
218         BP=0.
219         DO 430 K=1,3
220 430     BP=BP+FRAME(K,1)*PWREST(K)
221         BP=BP/FRAME(5,1)
222         DO 440 K=1,3
223           FRAME(K,2)=PWREST(K)+FRAME(K,1)*PWREST(4)/FRAME(5,1)
224      $    +FRAME(K,1)*BP/(FRAME(4,1)+FRAME(5,1))
225 440     CONTINUE
226         FRAME(4,2)=FRAME(4,1)*PWREST(4)/FRAME(5,1)+BP
227       ENDIF
228 C
229 C          Boost hadrons back to PP COM
230 C
231       DO 500 IP=NZERO,NPTCL
232         JET=IABS(IORIG(IP))/IPACK
233         IF(JET.GT.10) GO TO 500
234         IF(KEYS(6)) THEN
235           IF(IDJETS(JET).EQ.10) GO TO 500
236         ENDIF
237         IF(JET.EQ.0) THEN
238           IFR=1
239         ELSE
240           IFR=IFRAME(JET)
241         ENDIF
242         BP=0.
243         DO 510 K=1,3
244 510     BP=BP+FRAME(K,IFR)*PPTCL(K,IP)
245         BP=BP/FRAME(5,IFR)
246         DO 520 K=1,3
247 520     PPTCL(K,IP)=PPTCL(K,IP)+FRAME(K,IFR)*PPTCL(4,IP)/FRAME(5,IFR)
248      $  +FRAME(K,IFR)*BP/(FRAME(4,IFR)+FRAME(5,IFR))
249         PPTCL(4,IP)=FRAME(4,IFR)*PPTCL(4,IP)/FRAME(5,IFR)+BP
250 500   CONTINUE
251 C
252 C          Add intrinsic PT
253 C
254       IF(.NOT.KEYS(2)) THEN
255         PT2AVE=.1*SQRT(QSQ)
256         PTADD=SQRT(-PT2AVE*ALOG(RANF()))
257         PHIADD=2.*PI*RANF()
258         PADD(1)=2.*PTADD*COS(PHIADD)
259         PADD(2)=2.*PTADD*SIN(PHIADD)
260 C          Must use large and small components carefully to calculate
261 C          mass on 32-bit machines.
262         PALLPL=0.
263         PALLMN=0.
264         PALLX=0.
265         PALLY=0.
266         DO 525 IP=NZERO,NPTCL
267           PALLX=PALLX+PPTCL(1,IP)
268           PALLY=PALLY+PPTCL(2,IP)
269           IF(PPTCL(3,IP).GT.0.) THEN
270             PALLPL=PALLPL+(PPTCL(4,IP)+PPTCL(3,IP))
271             PALLMN=PALLMN+(PPTCL(1,IP)**2+PPTCL(2,IP)**2+PPTCL(5,IP)**2)
272      $      /(PPTCL(4,IP)+PPTCL(3,IP))
273           ELSE
274             PALLMN=PALLMN+(PPTCL(4,IP)-PPTCL(3,IP))
275             PALLPL=PALLPL+(PPTCL(1,IP)**2+PPTCL(2,IP)**2+PPTCL(5,IP)**2)
276      $      /(PPTCL(4,IP)-PPTCL(3,IP))
277           ENDIF
278 525     CONTINUE
279         POLD(1)=PALLX
280         POLD(2)=PALLY
281         POLD(3)=.5*(PALLPL-PALLMN)
282         POLD(4)=.5*(PALLPL+PALLMN)
283         POLD(5)=SQRT(PALLPL*PALLMN-PALLX**2-PALLY**2)
284         PNEW(1)=PADD(1)+POLD(1)
285         PNEW(2)=PADD(2)+POLD(2)
286         PNEW(3)=POLD(3)
287         PNEW(4)=SQRT(PNEW(1)**2+PNEW(2)**2+PNEW(3)**2+POLD(5)**2)
288         PNEW(5)=POLD(5)
289 C
290         DO 530 IP=NZERO,NPTCL
291           BP=0.
292           DO 531 K=1,3
293 531       BP=BP+POLD(K)*PPTCL(K,IP)
294           BP=BP/POLD(5)
295           DO 532 K=1,3
296 532       PPTCL(K,IP)=PPTCL(K,IP)-POLD(K)*PPTCL(4,IP)/POLD(5)
297      $    +POLD(K)*BP/(POLD(4)+POLD(5))
298           PPTCL(4,IP)=PPTCL(4,IP)*POLD(4)/POLD(5)-BP
299 C
300           BP=0.
301           DO 533 K=1,3
302 533       BP=BP+PNEW(K)*PPTCL(K,IP)
303           BP=BP/PNEW(5)
304           DO 534 K=1,3
305 534       PPTCL(K,IP)=PPTCL(K,IP)+PNEW(K)*PPTCL(4,IP)/PNEW(5)
306      $    +PNEW(K)*BP/(PNEW(4)+PNEW(5))
307           PPTCL(4,IP)=PPTCL(4,IP)*PNEW(4)/PNEW(5)+BP
308 530     CONTINUE
309 C
310 C            Add opposite PT to beam jets
311         DO 541 K=1,4
312 541     PBEAMS(K)=-PNEW(K)
313         PBEAMS(4)=PBEAMS(4)+ECM
314         PBEAMS(5)=PBEAMS(4)**2-PBEAMS(1)**2-PBEAMS(2)**2 -PBEAMS(3)**2
315         IF ( PBEAMS(5).GT.0 ) THEN
316           PBEAMS(5)=SQRT(PBEAMS(5))
317         ELSE
318           PBEAMS(4)=SQRT(PBEAMS(4)**2-PBEAMS(5)+4.)
319           PBEAMS(5)=2.
320         ENDIF
321       ENDIF
322 C
323 C          Decay hadrons
324 C
325       NP1=NZERO
326 600   NP2=NPTCL
327       DO 610 IP=NP1,NP2
328         NFIRST=NPTCL+1
329         JET=IABS(IORIG(IP))/IPACK
330         CALL DECAY(IP)
331         DO 620 IP1=NFIRST,NPTCL
332 620     IORIG(IP1)=ISIGN(IABS(IORIG(IP1))+IPACK*JET,IORIG(IP1))
333 610   CONTINUE
334       NP1=NP2+1
335       IF(NP1.LE.NPTCL) GO TO 600
336       RETURN
337 C
338 C          Error
339 C
340 9999  CALL PRTEVT(0)
341       WRITE(ITLIS,9998) NPTCL
342 9998  FORMAT(//' ERROR IN FRGMNT ... NPTCL > ',I6)
343       RETURN
344       END