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