EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
quench.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file quench.f
1 C
2 C
3 C
4 C
5  SUBROUTINE quench(JPJT,NTP)
6  dimension rdp(300),lqp(300),rdt(300),lqt(300)
7  common/hijcrdn/yp(3,300),yt(3,300)
8  SAVE /hijcrdn/
9  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
10  SAVE /hiparnt/
11 C
12  common/hijjet1/npj(300),kfpj(300,500),pjpx(300,500),
13  & pjpy(300,500),pjpz(300,500),pjpe(300,500),
14  & pjpm(300,500),ntj(300),kftj(300,500),
15  & pjtx(300,500),pjty(300,500),pjtz(300,500),
16  & pjte(300,500),pjtm(300,500)
17  SAVE /hijjet1/
18  common/hijjet2/nsg,njsg(900),iasg(900,3),k1sg(900,100),
19  & k2sg(900,100),pxsg(900,100),pysg(900,100),
20  & pzsg(900,100),pesg(900,100),pmsg(900,100)
21  SAVE /hijjet2/
22  common/histrng/nfp(300,15),pp(300,15),nft(300,15),pt(300,15)
23  SAVE /histrng/
24  common/ranseed/nseed
25  SAVE /ranseed/
26 C
27  bb=hint1(19) ! Uzhi
28  phi=hint1(20) ! Uzhi
29  bbx=bb*cos(phi) ! Uzhi
30  bby=bb*sin(phi) ! Uzhi
31 c
32  IF(ntp.EQ.2) go to 400
33  IF(ntp.EQ.3) go to 2000
34 C*******************************************************
35 C Jet interaction for proj jet in the direction PHIP
36 C******************************************************
37 C
38  IF(nfp(jpjt,7).NE.1) RETURN
39 
40  jp=jpjt
41  DO 290 i=1,npj(jp)
42  ptjet0=sqrt(pjpx(jp,i)**2+pjpy(jp,i)**2)
43  IF(ptjet0.LE.hipr1(11)) go to 290
44  ptot=sqrt(ptjet0*ptjet0+pjpz(jp,i)**2)
45  IF(ptot.LT.hipr1(8)) go to 290
46  phip=ulangl(pjpx(jp,i),pjpy(jp,i))
47 C******* find the wounded proj which can interact with jet***
48  kp=0
49  DO 100 i2=1,ihnt2(1)
50  IF(nfp(i2,5).NE.3 .OR. i2.EQ.jp) go to 100
51  dx=yp(1,i2)-yp(1,jp)
52  dy=yp(2,i2)-yp(2,jp)
53  phi=ulangl(dx,dy)
54  dphi=abs(phi-phip)
55  IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi ! Uzhi
56  IF(dphi.GE.hipr1(40)/2.0) go to 100
57  rd0=sqrt(dx*dx+dy*dy)
58  IF(rd0*sin(dphi).GT.hipr1(12)) go to 100
59  kp=kp+1
60  lqp(kp)=i2
61  rdp(kp)=cos(dphi)*rd0
62  100 CONTINUE
63 C******* rearrange according decending rd************
64  DO 110 i2=1,kp-1
65  DO 110 j2=i2+1,kp
66  IF(rdp(i2).LT.rdp(j2)) go to 110
67  rd=rdp(i2)
68  lq=lqp(i2)
69  rdp(i2)=rdp(j2)
70  lqp(i2)=lqp(j2)
71  rdp(j2)=rd
72  lqp(j2)=lq
73  110 CONTINUE
74 C****** find wounded targ which can interact with jet********
75  kt=0
76  DO 120 i2=1,ihnt2(3)
77  IF(nft(i2,5).NE.3) go to 120
78  dx=yt(1,i2)-yp(1,jp)-bbx
79  dy=yt(2,i2)-yp(2,jp)-bby
80  phi=ulangl(dx,dy)
81  dphi=abs(phi-phip)
82  IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi ! Uzhi
83  IF(dphi.GT.hipr1(40)/2.0) go to 120
84  rd0=sqrt(dx*dx+dy*dy)
85  IF(rd0*sin(dphi).GT.hipr1(12)) go to 120
86  kt=kt+1
87  lqt(kt)=i2
88  rdt(kt)=cos(dphi)*rd0
89  120 CONTINUE
90 C******* rearrange according decending rd************
91  DO 130 i2=1,kt-1
92  DO 130 j2=i2+1,kt
93  IF(rdt(i2).LT.rdt(j2)) go to 130
94  rd=rdt(i2)
95  lq=lqt(i2)
96  rdt(i2)=rdt(j2)
97  lqt(i2)=lqt(j2)
98  rdt(j2)=rd
99  lqt(j2)=lq
100  130 CONTINUE
101 
102  mp=0
103  mt=0
104  r0=0.0
105  nq=0
106  dp=0.0
107  ptot=sqrt(pjpx(jp,i)**2+pjpy(jp,i)**2+pjpz(jp,i)**2)
108  v1=pjpx(jp,i)/ptot
109  v2=pjpy(jp,i)/ptot
110  v3=pjpz(jp,i)/ptot
111 
112  200 rn=atl_ran(nseed)
113  210 IF(mt.GE.kt .AND. mp.GE.kp) go to 290
114  IF(mt.GE.kt) go to 220
115  IF(mp.GE.kp) go to 240
116  IF(rdp(mp+1).GT.rdt(mt+1)) go to 240
117  220 mp=mp+1
118  drr=rdp(mp)-r0
119  IF(rn.GE.1.0-exp(-drr/hipr1(13))) go to 210
120  dp=drr*hipr1(14)
121  IF(kfpj(jp,i).NE.21) dp=0.5*dp
122 C ********string tension of quark jet is 0.5 of gluon's
123  IF(dp.LE.0.2) go to 210
124  IF(ptot.LE.0.4) go to 290
125  IF(ptot.LE.dp) dp=ptot-0.2
126  de=dp
127 
128  IF(kfpj(jp,i).NE.21) THEN
129  prshu=pp(lqp(mp),1)**2+pp(lqp(mp),2)**2
130  & +pp(lqp(mp),3)**2
131  de=sqrt(pjpm(jp,i)**2+ptot**2)
132  & -sqrt(pjpm(jp,i)**2+(ptot-dp)**2)
133  ershu=(pp(lqp(mp),4)+de-dp)**2
134  amshu=ershu-prshu
135  IF(amshu.LT.hipr1(1)*hipr1(1)) go to 210
136  pp(lqp(mp),4)=sqrt(ershu)
137  pp(lqp(mp),5)=sqrt(amshu)
138  ENDIF
139 C ********reshuffle the energy when jet has mass
140  r0=rdp(mp)
141  dp1=dp*v1
142  dp2=dp*v2
143  dp3=dp*v3
144 C ********momentum and energy transfer from jet
145 
146  npj(lqp(mp))=npj(lqp(mp))+1
147  kfpj(lqp(mp),npj(lqp(mp)))=21
148  pjpx(lqp(mp),npj(lqp(mp)))=dp1
149  pjpy(lqp(mp),npj(lqp(mp)))=dp2
150  pjpz(lqp(mp),npj(lqp(mp)))=dp3
151  pjpe(lqp(mp),npj(lqp(mp)))=dp
152  pjpm(lqp(mp),npj(lqp(mp)))=0.0
153  go to 260
154 
155  240 mt=mt+1
156  drr=rdt(mt)-r0
157  IF(rn.GE.1.0-exp(-drr/hipr1(13))) go to 210
158  dp=drr*hipr1(14)
159  IF(dp.LE.0.2) go to 210
160  IF(ptot.LE.0.4) go to 290
161  IF(ptot.LE.dp) dp=ptot-0.2
162  de=dp
163 
164  IF(kfpj(jp,i).NE.21) THEN
165  prshu=pt(lqt(mt),1)**2+pt(lqt(mt),2)**2
166  & +pt(lqt(mt),3)**2
167  de=sqrt(pjpm(jp,i)**2+ptot**2)
168  & -sqrt(pjpm(jp,i)**2+(ptot-dp)**2)
169  ershu=(pt(lqt(mt),4)+de-dp)**2
170  amshu=ershu-prshu
171  IF(amshu.LT.hipr1(1)*hipr1(1)) go to 210
172  pt(lqt(mt),4)=sqrt(ershu)
173  pt(lqt(mt),5)=sqrt(amshu)
174  ENDIF
175 C ********reshuffle the energy when jet has mass
176 
177  r0=rdt(mt)
178  dp1=dp*v1
179  dp2=dp*v2
180  dp3=dp*v3
181 C ********momentum and energy transfer from jet
182  ntj(lqt(mt))=ntj(lqt(mt))+1
183  kftj(lqt(mt),ntj(lqt(mt)))=21
184  pjtx(lqt(mt),ntj(lqt(mt)))=dp1
185  pjty(lqt(mt),ntj(lqt(mt)))=dp2
186  pjtz(lqt(mt),ntj(lqt(mt)))=dp3
187  pjte(lqt(mt),ntj(lqt(mt)))=dp
188  pjtm(lqt(mt),ntj(lqt(mt)))=0.0
189 
190  260 pjpx(jp,i)=(ptot-dp)*v1
191  pjpy(jp,i)=(ptot-dp)*v2
192  pjpz(jp,i)=(ptot-dp)*v3
193  pjpe(jp,i)=pjpe(jp,i)-de
194 
195  ptot=ptot-dp
196  nq=nq+1
197  go to 200
198  290 CONTINUE
199 
200  RETURN
201 
202 C*******************************************************
203 C Jet interaction for target jet in the direction PHIT
204 C******************************************************
205 C
206 C******* find the wounded proj which can interact with jet***
207 
208  400 IF(nft(jpjt,7).NE.1) RETURN
209  jt=jpjt
210  DO 690 i=1,ntj(jt)
211  ptjet0=sqrt(pjtx(jt,i)**2+pjty(jt,i)**2)
212  IF(ptjet0.LE.hipr1(11)) go to 690
213  ptot=sqrt(ptjet0*ptjet0+pjtz(jt,i)**2)
214  IF(ptot.LT.hipr1(8)) go to 690
215  phit=ulangl(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(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(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=atl_ran(nseed)
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*hipr1(14)
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*hipr1(14)
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(ptot.LT.max(hipr1(1),hipr1(8))) go to 2690
385  phiq=ulangl(pxsg(isg,i),pysg(isg,i))
386  kp=0
387  DO 2500 i2=1,ihnt2(1)
388  IF(nfp(i2,5).NE.3.OR.i2.EQ.jp) go to 2500
389  dx=yp(1,i2)+bbx-xj
390  dy=yp(2,i2)+bby-yj
391  phi=ulangl(dx,dy)
392  dphi=abs(phi-phiq)
393  IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi ! Uzhi
394  IF(dphi.GT.hipr1(40)/2.0) go to 2500
395  rd0=sqrt(dx*dx+dy*dy)
396  IF(rd0*sin(dphi).GT.hipr1(12)) go to 2500
397  kp=kp+1
398  lqp(kp)=i2
399  rdp(kp)=cos(dphi)*rd0
400  2500 CONTINUE
401 C******* rearrange according to decending rd************
402  DO 2510 i2=1,kp-1
403  DO 2510 j2=i2+1,kp
404  IF(rdp(i2).LT.rdp(j2)) go to 2510
405  rd=rdp(i2)
406  lq=lqp(i2)
407  rdp(i2)=rdp(j2)
408  lqp(i2)=lqp(j2)
409  rdp(j2)=rd
410  lqp(j2)=lq
411  2510 CONTINUE
412 C****** find wounded targ which can interact with jet********
413  kt=0
414  DO 2520 i2=1,ihnt2(3)
415  IF(nft(i2,5).NE.3 .OR. i2.EQ.jt) go to 2520
416  dx=yt(1,i2)-xj
417  dy=yt(2,i2)-yj
418  phi=ulangl(dx,dy)
419  dphi=abs(phi-phiq)
420  IF(dphi.GE.hipr1(40)) dphi=2.*hipr1(40)-dphi ! Uzhi
421  IF(dphi.GT.hipr1(40)/2.0) go to 2520
422  rd0=sqrt(dx*dx+dy*dy)
423  IF(rd0*sin(dphi).GT.hipr1(12)) go to 2520
424  kt=kt+1
425  lqt(kt)=i2
426  rdt(kt)=cos(dphi)*rd0
427  2520 CONTINUE
428 C******* rearrange according to decending rd************
429  DO 2530 i2=1,kt-1
430  DO 2530 j2=i2+1,kt
431  IF(rdt(i2).LT.rdt(j2)) go to 2530
432  rd=rdt(i2)
433  lq=lqt(i2)
434  rdt(i2)=rdt(j2)
435  lqt(i2)=lqt(j2)
436  rdt(j2)=rd
437  lqt(j2)=lq
438  2530 CONTINUE
439 
440  mp=0
441  mt=0
442  nq=0
443  dp=0.0
444  r0=0.0
445  ptot=sqrt(pxsg(isg,i)**2+pysg(isg,i)**2
446  & +pzsg(isg,i)**2)
447  v1=pxsg(isg,i)/ptot
448  v2=pysg(isg,i)/ptot
449  v3=pzsg(isg,i)/ptot
450 
451  2600 rn=atl_ran(nseed)
452  2610 IF(mt.GE.kt .AND. mp.GE.kp) go to 2690
453  IF(mt.GE.kt) go to 2620
454  IF(mp.GE.kp) go to 2640
455  IF(rdp(mp+1).GT.rdt(mt+1)) go to 2640
456  2620 mp=mp+1
457  drr=rdp(mp)-r0
458  IF(rn.GE.1.0-exp(-drr/hipr1(13))) go to 2610
459  dp=drr*hipr1(14)/2.0
460  IF(dp.LE.0.2) go to 2610
461  IF(ptot.LE.0.4) go to 2690
462  IF(ptot.LE.dp) dp=ptot-0.2
463  de=dp
464 C
465  IF(k2sg(isg,i).NE.21) THEN
466  IF(ptot.LT.dp+hipr1(1)) go to 2690
467  prshu=pp(lqp(mp),1)**2+pp(lqp(mp),2)**2
468  & +pp(lqp(mp),3)**2
469  de=sqrt(pmsg(isg,i)**2+ptot**2)
470  & -sqrt(pmsg(isg,i)**2+(ptot-dp)**2)
471  ershu=(pp(lqp(mp),4)+de-dp)**2
472  amshu=ershu-prshu
473  IF(amshu.LT.hipr1(1)*hipr1(1)) go to 2610
474  pp(lqp(mp),4)=sqrt(ershu)
475  pp(lqp(mp),5)=sqrt(amshu)
476  ENDIF
477 C ********reshuffle the energy when jet has mass
478 C
479  r0=rdp(mp)
480  dp1=dp*v1
481  dp2=dp*v2
482  dp3=dp*v3
483 C ********momentum and energy transfer from jet
484  npj(lqp(mp))=npj(lqp(mp))+1
485  kfpj(lqp(mp),npj(lqp(mp)))=21
486  pjpx(lqp(mp),npj(lqp(mp)))=dp1
487  pjpy(lqp(mp),npj(lqp(mp)))=dp2
488  pjpz(lqp(mp),npj(lqp(mp)))=dp3
489  pjpe(lqp(mp),npj(lqp(mp)))=dp
490  pjpm(lqp(mp),npj(lqp(mp)))=0.0
491 
492  go to 2660
493 
494  2640 mt=mt+1
495  drr=rdt(mt)-r0
496  IF(rn.GE.1.0-exp(-drr/hipr1(13))) go to 2610
497  dp=drr*hipr1(14)
498  IF(dp.LE.0.2) go to 2610
499  IF(ptot.LE.0.4) go to 2690
500  IF(ptot.LE.dp) dp=ptot-0.2
501  de=dp
502 
503  IF(k2sg(isg,i).NE.21) THEN
504  IF(ptot.LT.dp+hipr1(1)) go to 2690
505  prshu=pt(lqt(mt),1)**2+pt(lqt(mt),2)**2
506  & +pt(lqt(mt),3)**2
507  de=sqrt(pmsg(isg,i)**2+ptot**2)
508  & -sqrt(pmsg(isg,i)**2+(ptot-dp)**2)
509  ershu=(pt(lqt(mt),4)+de-dp)**2
510  amshu=ershu-prshu
511  IF(amshu.LT.hipr1(1)*hipr1(1)) go to 2610
512  pt(lqt(mt),4)=sqrt(ershu)
513  pt(lqt(mt),5)=sqrt(amshu)
514  ENDIF
515 C ********reshuffle the energy when jet has mass
516 
517  r0=rdt(mt)
518  dp1=dp*v1
519  dp2=dp*v2
520  dp3=dp*v3
521 C ********momentum and energy transfer from jet
522  ntj(lqt(mt))=ntj(lqt(mt))+1
523  kftj(lqt(mt),ntj(lqt(mt)))=21
524  pjtx(lqt(mt),ntj(lqt(mt)))=dp1
525  pjty(lqt(mt),ntj(lqt(mt)))=dp2
526  pjtz(lqt(mt),ntj(lqt(mt)))=dp3
527  pjte(lqt(mt),ntj(lqt(mt)))=dp
528  pjtm(lqt(mt),ntj(lqt(mt)))=0.0
529 
530  2660 pxsg(isg,i)=(ptot-dp)*v1
531  pysg(isg,i)=(ptot-dp)*v2
532  pzsg(isg,i)=(ptot-dp)*v3
533  pesg(isg,i)=pesg(isg,i)-de
534 
535  ptot=ptot-dp
536  nq=nq+1
537  go to 2600
538  2690 CONTINUE
539  RETURN
540  END