Include corrections by V.Uzhinskii JINR, Dubna, Russia from V1.37
[u/mrichter/AliRoot.git] / HIJING / hijing1_36 / quench.F
1 * $Id$
2 C
3 C
4 C
5 C
6         SUBROUTINE QUENCH(JPJT,NTP)
7         DIMENSION RDP(300),LQP(300),RDT(300),LQT(300)
8 #include "hijcrdn.inc"
9 #include "hiparnt.inc"
10 C
11 #include "hijjet1.inc"
12 #include "hijjet2.inc"
13 #include "histrng.inc"
14 C
15         SAVE
16 C
17         BB=HINT1(19)                                            ! Uzhi
18         PHI=HINT1(20)                                           ! Uzhi
19         BBX=BB*COS(PHI)                                         ! Uzhi
20         BBY=BB*SIN(PHI)                                         ! Uzhi
21 c
22         IF(NTP.EQ.2) GO TO 400
23         IF(NTP.EQ.3) GO TO 2000 
24 C*******************************************************
25 C Jet interaction for proj jet in the direction PHIP
26 C******************************************************
27 C
28         IF(NFP(JPJT,7).NE.1) RETURN
29
30         JP=JPJT
31         DO 290 I=1,NPJ(JP)
32            PTJET0=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2)
33            IF(PTJET0.LE.HIPR1(11)) GO TO 290
34            PTOT=SQRT(PTJET0*PTJET0+PJPZ(JP,I)**2)
35            IF(PTOT.LT.HIPR1(8)) GO TO 290
36            PHIP=ULANGL_HIJING(PJPX(JP,I),PJPY(JP,I))
37 C******* find the wounded proj which can interact with jet***
38            KP=0
39            DO 100 I2=1,IHNT2(1)
40               IF(NFP(I2,5).NE.3 .OR. I2.EQ.JP) GO TO 100
41               DX=YP(1,I2)-YP(1,JP)
42               DY=YP(2,I2)-YP(2,JP)
43               PHI=ULANGL_HIJING(DX,DY)
44               DPHI=ABS(PHI-PHIP)
45               IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI      ! Uzhi
46               IF(DPHI.GE.HIPR1(40)/2.0) GO TO 100
47               RD0=SQRT(DX*DX+DY*DY)
48               IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 100
49               KP=KP+1
50               LQP(KP)=I2
51               RDP(KP)=COS(DPHI)*RD0
52  100       CONTINUE
53 C*******        rearrange according decending rd************
54            DO 110 I2=1,KP-1
55               DO 110 J2=I2+1,KP
56                  IF(RDP(I2).LT.RDP(J2)) GO TO 110
57                  RD=RDP(I2)
58                  LQ=LQP(I2)
59                  RDP(I2)=RDP(J2)
60                  LQP(I2)=LQP(J2)
61                  RDP(J2)=RD
62                  LQP(J2)=LQ
63  110          CONTINUE
64 C****** find wounded targ which can interact with jet********
65               KT=0
66               DO 120 I2=1,IHNT2(3)
67                  IF(NFT(I2,5).NE.3) GO TO 120
68                  DX=YT(1,I2)-YP(1,JP)
69                  DY=YT(2,I2)-YP(2,JP)
70                  PHI=ULANGL_HIJING(DX,DY)
71                  DPHI=ABS(PHI-PHIP)
72                  IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
73                  IF(DPHI.GT.HIPR1(40)/2.0) GO TO 120
74                  RD0=SQRT(DX*DX+DY*DY)
75                  IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 120
76                  KT=KT+1
77                  LQT(KT)=I2
78                  RDT(KT)=COS(DPHI)*RD0
79  120          CONTINUE
80 C*******        rearrange according decending rd************
81               DO 130 I2=1,KT-1
82                  DO 130 J2=I2+1,KT
83                     IF(RDT(I2).LT.RDT(J2)) GO TO 130
84                     RD=RDT(I2)
85                     LQ=LQT(I2)
86                     RDT(I2)=RDT(J2)
87                     LQT(I2)=LQT(J2)
88                     RDT(J2)=RD
89                     LQT(J2)=LQ
90  130             CONTINUE
91                 
92                  MP=0
93                  MT=0
94                  R0=0.0
95                  NQ=0
96                  DP=0.0
97                  PTOT=SQRT(PJPX(JP,I)**2+PJPY(JP,I)**2+PJPZ(JP,I)**2)
98                  V1=PJPX(JP,I)/PTOT
99                  V2=PJPY(JP,I)/PTOT
100                  V3=PJPZ(JP,I)/PTOT
101
102  200             RN=RLU_HIJING(0)
103  210             IF(MT.GE.KT .AND. MP.GE.KP) GO TO 290
104                  IF(MT.GE.KT) GO TO 220
105                  IF(MP.GE.KP) GO TO 240
106                  IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 240
107  220             MP=MP+1
108                  DRR=RDP(MP)-R0
109                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
110                  DP=DRR*HIPR1(14)
111                  IF(KFPJ(JP,I).NE.21) DP=0.5*DP
112 C       ********string tension of quark jet is 0.5 of gluon's 
113                  IF(DP.LE.0.2) GO TO 210
114                  IF(PTOT.LE.0.4) GO TO 290
115                  IF(PTOT.LE.DP) DP=PTOT-0.2
116                  DE=DP
117
118                  IF(KFPJ(JP,I).NE.21) THEN
119                     PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
120      &                   +PP(LQP(MP),3)**2
121                     DE=SQRT(PJPM(JP,I)**2+PTOT**2)
122      &                  -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
123                     ERSHU=(PP(LQP(MP),4)+DE-DP)**2
124                     AMSHU=ERSHU-PRSHU
125                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
126                     PP(LQP(MP),4)=SQRT(ERSHU)
127                     PP(LQP(MP),5)=SQRT(AMSHU)
128                  ENDIF
129 C               ********reshuffle the energy when jet has mass
130                  R0=RDP(MP)
131                  DP1=DP*V1
132                  DP2=DP*V2
133                  DP3=DP*V3
134 C               ********momentum and energy transfer from jet
135                  
136                  NPJ(LQP(MP))=NPJ(LQP(MP))+1
137                  KFPJ(LQP(MP),NPJ(LQP(MP)))=21
138                  PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
139                  PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
140                  PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
141                  PJPE(LQP(MP),NPJ(LQP(MP)))=DP
142                  PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
143                  GO TO 260
144
145  240             MT=MT+1
146                  DRR=RDT(MT)-R0
147                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 210
148                  DP=DRR*HIPR1(14)
149                  IF(DP.LE.0.2) GO TO 210
150                  IF(PTOT.LE.0.4) GO TO 290
151                  IF(PTOT.LE.DP) DP=PTOT-0.2
152                  DE=DP
153
154                  IF(KFPJ(JP,I).NE.21) THEN
155                     PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
156      &                   +PT(LQT(MT),3)**2
157                     DE=SQRT(PJPM(JP,I)**2+PTOT**2)
158      &                  -SQRT(PJPM(JP,I)**2+(PTOT-DP)**2)
159                     ERSHU=(PT(LQT(MT),4)+DE-DP)**2
160                     AMSHU=ERSHU-PRSHU
161                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 210
162                     PT(LQT(MT),4)=SQRT(ERSHU)
163                     PT(LQT(MT),5)=SQRT(AMSHU)
164                  ENDIF
165 C               ********reshuffle the energy when jet has mass
166
167                  R0=RDT(MT)
168                  DP1=DP*V1
169                  DP2=DP*V2
170                  DP3=DP*V3
171 C               ********momentum and energy transfer from jet
172                  NTJ(LQT(MT))=NTJ(LQT(MT))+1
173                  KFTJ(LQT(MT),NTJ(LQT(MT)))=21
174                  PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
175                  PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
176                  PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
177                  PJTE(LQT(MT),NTJ(LQT(MT)))=DP
178                  PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
179
180  260             PJPX(JP,I)=(PTOT-DP)*V1
181                  PJPY(JP,I)=(PTOT-DP)*V2
182                  PJPZ(JP,I)=(PTOT-DP)*V3
183                  PJPE(JP,I)=PJPE(JP,I)-DE
184
185                  PTOT=PTOT-DP
186                  NQ=NQ+1
187                  GO TO 200
188  290          CONTINUE
189
190               RETURN
191
192 C*******************************************************
193 C Jet interaction for target jet in the direction PHIT
194 C******************************************************
195 C
196 C******* find the wounded proj which can interact with jet***
197
198  400          IF(NFT(JPJT,7).NE.1) RETURN
199               JT=JPJT
200               DO 690 I=1,NTJ(JT)
201                  PTJET0=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2)
202                  IF(PTJET0.LE.HIPR1(11)) GO TO 690
203                  PTOT=SQRT(PTJET0*PTJET0+PJTZ(JT,I)**2)
204                  IF(PTOT.LT.HIPR1(8)) GO TO 690
205                  PHIT=ULANGL_HIJING(PJTX(JT,I),PJTY(JT,I))
206                  KP=0
207                  DO 500 I2=1,IHNT2(1)
208                     IF(NFP(I2,5).NE.3) GO TO 500
209                     DX=YP(1,I2)-YT(1,JT)
210                     DY=YP(2,I2)-YT(2,JT)
211                     PHI=ULANGL_HIJING(DX,DY)
212                     DPHI=ABS(PHI-PHIT)
213                     IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
214                     IF(DPHI.GT.HIPR1(40)/2.0) GO TO 500
215                     RD0=SQRT(DX*DX+DY*DY)
216                     IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 500
217                     KP=KP+1
218                     LQP(KP)=I2
219                     RDP(KP)=COS(DPHI)*RD0
220  500             CONTINUE
221 C*******        rearrange according to decending rd************
222                  DO 510 I2=1,KP-1
223                     DO 510 J2=I2+1,KP
224                        IF(RDP(I2).LT.RDP(J2)) GO TO 510
225                        RD=RDP(I2)
226                        LQ=LQP(I2)
227                        RDP(I2)=RDP(J2)
228                        LQP(I2)=LQP(J2)
229                        RDP(J2)=RD
230                        LQP(J2)=LQ
231  510                CONTINUE
232 C****** find wounded targ which can interact with jet********
233                     KT=0
234                     DO 520 I2=1,IHNT2(3)
235                        IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 520
236                        DX=YT(1,I2)-YT(1,JT)
237                        DY=YT(2,I2)-YT(2,JT)
238                        PHI=ULANGL_HIJING(DX,DY)
239                        DPHI=ABS(PHI-PHIT)
240                        IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
241                        IF(DPHI.GT.HIPR1(40)/2.0) GO TO 520
242                        RD0=SQRT(DX*DX+DY*DY)
243                        IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 520
244                        KT=KT+1
245                        LQT(KT)=I2
246                        RDT(KT)=COS(DPHI)*RD0
247  520                CONTINUE
248 C*******        rearrange according to decending rd************
249                     DO 530 I2=1,KT-1
250                        DO 530 J2=I2+1,KT
251                           IF(RDT(I2).LT.RDT(J2)) GO TO 530
252                           RD=RDT(I2)
253                           LQ=LQT(I2)
254                           RDT(I2)=RDT(J2)
255                           LQT(I2)=LQT(J2)
256                           RDT(J2)=RD
257                           LQT(J2)=LQ
258  530                   CONTINUE
259                        
260                        MP=0
261                        MT=0
262                        NQ=0
263                        DP=0.0
264                        R0=0.0
265                 PTOT=SQRT(PJTX(JT,I)**2+PJTY(JT,I)**2+PJTZ(JT,I)**2)
266                 V1=PJTX(JT,I)/PTOT
267                 V2=PJTY(JT,I)/PTOT
268                 V3=PJTZ(JT,I)/PTOT
269
270  600            RN=RLU_HIJING(0)
271  610            IF(MT.GE.KT .AND. MP.GE.KP) GO TO 690
272                 IF(MT.GE.KT) GO TO 620
273                 IF(MP.GE.KP) GO TO 640
274                 IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 640
275 620             MP=MP+1
276                 DRR=RDP(MP)-R0
277                 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
278                 DP=DRR*HIPR1(14)
279                 IF(KFTJ(JT,I).NE.21) DP=0.5*DP
280 C       ********string tension of quark jet is 0.5 of gluon's 
281                 IF(DP.LE.0.2) GO TO 610
282                 IF(PTOT.LE.0.4) GO TO 690
283                 IF(PTOT.LE.DP) DP=PTOT-0.2
284                 DE=DP
285 C
286                 IF(KFTJ(JT,I).NE.21) THEN
287                    PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
288      &                   +PP(LQP(MP),3)**2
289                    DE=SQRT(PJTM(JT,I)**2+PTOT**2)
290      &               -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
291                    ERSHU=(PP(LQP(MP),4)+DE-DP)**2
292                    AMSHU=ERSHU-PRSHU
293                    IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
294                    PP(LQP(MP),4)=SQRT(ERSHU)
295                    PP(LQP(MP),5)=SQRT(AMSHU)
296                 ENDIF
297 C               ********reshuffle the energy when jet has mass
298 C
299                 R0=RDP(MP)
300                 DP1=DP*V1
301                 DP2=DP*V2
302                 DP3=DP*V3
303 C               ********momentum and energy transfer from jet
304                 NPJ(LQP(MP))=NPJ(LQP(MP))+1
305                 KFPJ(LQP(MP),NPJ(LQP(MP)))=21
306                 PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
307                 PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
308                 PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
309                 PJPE(LQP(MP),NPJ(LQP(MP)))=DP
310                 PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
311
312                 GO TO 660
313
314 640             MT=MT+1
315                 DRR=RDT(MT)-R0
316                 IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 610
317                 DP=DRR*HIPR1(14)
318                 IF(DP.LE.0.2) GO TO 610
319                 IF(PTOT.LE.0.4) GO TO 690
320                 IF(PTOT.LE.DP) DP=PTOT-0.2
321                 DE=DP
322
323                 IF(KFTJ(JT,I).NE.21) THEN
324                    PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
325      &                   +PT(LQT(MT),3)**2
326                    DE=SQRT(PJTM(JT,I)**2+PTOT**2)
327      &               -SQRT(PJTM(JT,I)**2+(PTOT-DP)**2)
328                    ERSHU=(PT(LQT(MT),4)+DE-DP)**2
329                    AMSHU=ERSHU-PRSHU
330                    IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 610
331                    PT(LQT(MT),4)=SQRT(ERSHU)
332                    PT(LQT(MT),5)=SQRT(AMSHU)
333                 ENDIF
334 C               ********reshuffle the energy when jet has mass
335
336                 R0=RDT(MT)
337                 DP1=DP*V1
338                 DP2=DP*V2
339                 DP3=DP*V3
340 C               ********momentum and energy transfer from jet
341                 NTJ(LQT(MT))=NTJ(LQT(MT))+1
342                 KFTJ(LQT(MT),NTJ(LQT(MT)))=21
343                 PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
344                 PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
345                 PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
346                 PJTE(LQT(MT),NTJ(LQT(MT)))=DP
347                 PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
348
349 660             PJTX(JT,I)=(PTOT-DP)*V1
350                 PJTY(JT,I)=(PTOT-DP)*V2
351                 PJTZ(JT,I)=(PTOT-DP)*V3
352                 PJTE(JT,I)=PJTE(JT,I)-DE
353
354                 PTOT=PTOT-DP
355                 NQ=NQ+1
356                 GO TO 600
357 690     CONTINUE
358         RETURN
359 C********************************************************
360 C       Q-QBAR jet interaction
361 C********************************************************
362 2000    ISG=JPJT
363         IF(IASG(ISG,3).NE.1) RETURN
364 C
365         JP=IASG(ISG,1)
366         JT=IASG(ISG,2)
367         XJ=(YP(1,JP)+YT(1,JT))/2.0
368         YJ=(YP(2,JP)+YT(2,JT))/2.0
369         DO 2690 I=1,NJSG(ISG)
370            PTJET0=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2)
371            IF(PTJET0.LE.HIPR1(11).OR.PESG(ISG,I).LT.HIPR1(1))
372      &            GO TO 2690
373            PTOT=SQRT(PTJET0*PTJET0+PZSG(ISG,I)**2)
374            IF(PTOT.LT.MAX(HIPR1(1),HIPR1(8))) GO TO 2690
375            PHIQ=ULANGL_HIJING(PXSG(ISG,I),PYSG(ISG,I))
376            KP=0
377            DO 2500 I2=1,IHNT2(1)
378               IF(NFP(I2,5).NE.3.OR.I2.EQ.JP) GO TO 2500
379               DX=YP(1,I2)-XJ
380               DY=YP(2,I2)-YJ
381               PHI=ULANGL_HIJING(DX,DY)
382               DPHI=ABS(PHI-PHIQ)
383               IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI      ! Uzhi
384               IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2500
385               RD0=SQRT(DX*DX+DY*DY)
386               IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2500
387               KP=KP+1
388               LQP(KP)=I2
389               RDP(KP)=COS(DPHI)*RD0
390  2500      CONTINUE
391 C*******        rearrange according to decending rd************
392            DO 2510 I2=1,KP-1
393               DO 2510 J2=I2+1,KP
394                  IF(RDP(I2).LT.RDP(J2)) GO TO 2510
395                  RD=RDP(I2)
396                  LQ=LQP(I2)
397                  RDP(I2)=RDP(J2)
398                  LQP(I2)=LQP(J2)
399                  RDP(J2)=RD
400                  LQP(J2)=LQ
401  2510         CONTINUE
402 C****** find wounded targ which can interact with jet********
403               KT=0
404               DO 2520 I2=1,IHNT2(3)
405                  IF(NFT(I2,5).NE.3 .OR. I2.EQ.JT) GO TO 2520
406                  DX=YT(1,I2)-XJ
407                  DY=YT(2,I2)-YJ
408                  PHI=ULANGL_HIJING(DX,DY)
409                  DPHI=ABS(PHI-PHIQ)
410                  IF(DPHI.GE.HIPR1(40)) DPHI=2.*HIPR1(40)-DPHI ! Uzhi
411                  IF(DPHI.GT.HIPR1(40)/2.0) GO TO 2520
412                  RD0=SQRT(DX*DX+DY*DY)
413                  IF(RD0*SIN(DPHI).GT.HIPR1(12)) GO TO 2520
414                  KT=KT+1
415                  LQT(KT)=I2
416                  RDT(KT)=COS(DPHI)*RD0
417  2520         CONTINUE
418 C*******        rearrange according to decending rd************
419               DO 2530 I2=1,KT-1
420                  DO 2530 J2=I2+1,KT
421                     IF(RDT(I2).LT.RDT(J2)) GO TO 2530
422                     RD=RDT(I2)
423                     LQ=LQT(I2)
424                     RDT(I2)=RDT(J2)
425                     LQT(I2)=LQT(J2)
426                     RDT(J2)=RD
427                     LQT(J2)=LQ
428  2530            CONTINUE
429                 
430                  MP=0
431                  MT=0
432                  NQ=0
433                  DP=0.0
434                  R0=0.0
435                  PTOT=SQRT(PXSG(ISG,I)**2+PYSG(ISG,I)**2
436      &                +PZSG(ISG,I)**2)
437                  V1=PXSG(ISG,I)/PTOT
438                  V2=PYSG(ISG,I)/PTOT
439                  V3=PZSG(ISG,I)/PTOT
440
441  2600            RN=RLU_HIJING(0)
442  2610            IF(MT.GE.KT .AND. MP.GE.KP) GO TO 2690
443                  IF(MT.GE.KT) GO TO 2620
444                  IF(MP.GE.KP) GO TO 2640
445                  IF(RDP(MP+1).GT.RDT(MT+1)) GO TO 2640
446  2620            MP=MP+1
447                  DRR=RDP(MP)-R0
448                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
449                  DP=DRR*HIPR1(14)/2.0
450                  IF(DP.LE.0.2) GO TO 2610
451                  IF(PTOT.LE.0.4) GO TO 2690
452                  IF(PTOT.LE.DP) DP=PTOT-0.2
453                  DE=DP
454 C
455                  IF(K2SG(ISG,I).NE.21) THEN
456                     IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
457                     PRSHU=PP(LQP(MP),1)**2+PP(LQP(MP),2)**2
458      &                    +PP(LQP(MP),3)**2
459                     DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
460      &                 -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
461                     ERSHU=(PP(LQP(MP),4)+DE-DP)**2
462                     AMSHU=ERSHU-PRSHU
463                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
464                     PP(LQP(MP),4)=SQRT(ERSHU)
465                     PP(LQP(MP),5)=SQRT(AMSHU)
466                  ENDIF
467 C               ********reshuffle the energy when jet has mass
468 C
469                  R0=RDP(MP)
470                  DP1=DP*V1
471                  DP2=DP*V2
472                  DP3=DP*V3
473 C               ********momentum and energy transfer from jet
474                  NPJ(LQP(MP))=NPJ(LQP(MP))+1
475                  KFPJ(LQP(MP),NPJ(LQP(MP)))=21
476                  PJPX(LQP(MP),NPJ(LQP(MP)))=DP1
477                  PJPY(LQP(MP),NPJ(LQP(MP)))=DP2
478                  PJPZ(LQP(MP),NPJ(LQP(MP)))=DP3
479                  PJPE(LQP(MP),NPJ(LQP(MP)))=DP
480                  PJPM(LQP(MP),NPJ(LQP(MP)))=0.0
481
482                  GO TO 2660
483
484  2640            MT=MT+1
485                  DRR=RDT(MT)-R0
486                  IF(RN.GE.1.0-EXP(-DRR/HIPR1(13))) GO TO 2610
487                  DP=DRR*HIPR1(14)
488                  IF(DP.LE.0.2) GO TO 2610
489                  IF(PTOT.LE.0.4) GO TO 2690
490                  IF(PTOT.LE.DP) DP=PTOT-0.2
491                  DE=DP
492
493                  IF(K2SG(ISG,I).NE.21) THEN
494                     IF(PTOT.LT.DP+HIPR1(1)) GO TO 2690
495                     PRSHU=PT(LQT(MT),1)**2+PT(LQT(MT),2)**2
496      &                    +PT(LQT(MT),3)**2
497                     DE=SQRT(PMSG(ISG,I)**2+PTOT**2)
498      &                 -SQRT(PMSG(ISG,I)**2+(PTOT-DP)**2)
499                     ERSHU=(PT(LQT(MT),4)+DE-DP)**2
500                     AMSHU=ERSHU-PRSHU
501                     IF(AMSHU.LT.HIPR1(1)*HIPR1(1)) GO TO 2610
502                     PT(LQT(MT),4)=SQRT(ERSHU)
503                     PT(LQT(MT),5)=SQRT(AMSHU)
504                  ENDIF
505 C               ********reshuffle the energy when jet has mass
506
507                  R0=RDT(MT)
508                  DP1=DP*V1
509                  DP2=DP*V2
510                  DP3=DP*V3
511 C               ********momentum and energy transfer from jet
512                  NTJ(LQT(MT))=NTJ(LQT(MT))+1
513                  KFTJ(LQT(MT),NTJ(LQT(MT)))=21
514                  PJTX(LQT(MT),NTJ(LQT(MT)))=DP1
515                  PJTY(LQT(MT),NTJ(LQT(MT)))=DP2
516                  PJTZ(LQT(MT),NTJ(LQT(MT)))=DP3
517                  PJTE(LQT(MT),NTJ(LQT(MT)))=DP
518                  PJTM(LQT(MT),NTJ(LQT(MT)))=0.0
519
520  2660            PXSG(ISG,I)=(PTOT-DP)*V1
521                  PYSG(ISG,I)=(PTOT-DP)*V2
522                  PZSG(ISG,I)=(PTOT-DP)*V3
523                  PESG(ISG,I)=PESG(ISG,I)-DE
524
525                  PTOT=PTOT-DP
526                  NQ=NQ+1
527                  GO TO 2600
528  2690   CONTINUE
529         RETURN
530         END