Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
luexec.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file luexec.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE luexec
5 
6 C...Purpose: to administrate the fragmentation and decay chain.
7  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
8  SAVE /lujets/
9  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
10  SAVE /ludat1/
11  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
12  SAVE /ludat2/
13  common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5)
14  SAVE /ludat3/
15  dimension ps(2,6)
16 
17 C...Initialize and reset.
18  mstu(24)=0
19  IF(mstu(12).GE.1) CALL lulist(0)
20  mstu(31)=mstu(31)+1
21  mstu(1)=0
22  mstu(2)=0
23  mstu(3)=0
24  mcons=1
25 
26 C...Sum up momentum, energy and charge for starting entries.
27  nsav=n
28  DO 100 i=1,2
29  DO 100 j=1,6
30  100 ps(i,j)=0.
31  DO 120 i=1,n
32  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 120
33  DO 110 j=1,4
34  110 ps(1,j)=ps(1,j)+p(i,j)
35  ps(1,6)=ps(1,6)+luchge(k(i,2))
36  120 CONTINUE
37  paru(21)=ps(1,4)
38 
39 C...Prepare system for subsequent fragmentation/decay.
40  CALL luprep(0)
41 
42 C...Loop through jet fragmentation and particle decays.
43  mbe=0
44  130 mbe=mbe+1
45  ip=0
46  140 ip=ip+1
47  kc=0
48  IF(k(ip,1).GT.0.AND.k(ip,1).LE.10) kc=lucomp(k(ip,2))
49  IF(kc.EQ.0) THEN
50 
51 C...Particle decay if unstable and allowed. Save long-lived particle
52 C...decays until second pass after Bose-Einstein effects.
53  ELSEIF(kchg(kc,2).EQ.0) THEN
54  IF(mstj(21).GE.1.AND.mdcy(kc,1).GE.1.AND.(mstj(51).LE.0.OR.mbe.
55  & eq.2.OR.pmas(kc,2).GE.parj(91).OR.iabs(k(ip,2)).EQ.311))
56  & CALL ludecy(ip)
57 
58 C...Decay products may develop a shower.
59  IF(mstj(92).GT.0) THEN
60  ip1=mstj(92)
61  qmax=sqrt(max(0.,(p(ip1,4)+p(ip1+1,4))**2-(p(ip1,1)+p(ip1+1,
62  & 1))**2-(p(ip1,2)+p(ip1+1,2))**2-(p(ip1,3)+p(ip1+1,3))**2))
63  CALL lushow(ip1,ip1+1,qmax)
64  CALL luprep(ip1)
65  mstj(92)=0
66  ELSEIF(mstj(92).LT.0) THEN
67  ip1=-mstj(92)
68  CALL lushow(ip1,-3,p(ip,5))
69  CALL luprep(ip1)
70  mstj(92)=0
71  ENDIF
72 
73 C...Jet fragmentation: string or independent fragmentation.
74  ELSEIF(k(ip,1).EQ.1.OR.k(ip,1).EQ.2) THEN
75  mfrag=mstj(1)
76  IF(mfrag.GE.1.AND.k(ip,1).EQ.1) mfrag=2
77  IF(mstj(21).GE.2.AND.k(ip,1).EQ.2.AND.n.GT.ip) THEN
78  IF(k(ip+1,1).EQ.1.AND.k(ip+1,3).EQ.k(ip,3).AND.
79  & k(ip,3).GT.0.AND.k(ip,3).LT.ip) THEN
80  IF(kchg(lucomp(k(k(ip,3),2)),2).EQ.0) mfrag=min(1,mfrag)
81  ENDIF
82  ENDIF
83  IF(mfrag.EQ.1) CALL lustrf(ip)
84  IF(mfrag.EQ.2) CALL luindf(ip)
85  IF(mfrag.EQ.2.AND.k(ip,1).EQ.1) mcons=0
86  IF(mfrag.EQ.2.AND.(mstj(3).LE.0.OR.mod(mstj(3),5).EQ.0)) mcons=0
87  ENDIF
88 
89 C...Loop back if enough space left in LUJETS and no error abort.
90  IF(mstu(24).NE.0.AND.mstu(21).GE.2) THEN
91  ELSEIF(ip.LT.n.AND.n.LT.mstu(4)-20-mstu(32)) THEN
92  goto 140
93  ELSEIF(ip.LT.n) THEN
94  CALL luerrm(11,'(LUEXEC:) no more memory left in LUJETS')
95  ENDIF
96 
97 C...Include simple Bose-Einstein effect parametrization if desired.
98  IF(mbe.EQ.1.AND.mstj(51).GE.1) THEN
99  CALL luboei(nsav)
100  goto 130
101  ENDIF
102 
103 C...Check that momentum, energy and charge were conserved.
104  DO 160 i=1,n
105  IF(k(i,1).LE.0.OR.k(i,1).GT.10) goto 160
106  DO 150 j=1,4
107  150 ps(2,j)=ps(2,j)+p(i,j)
108  ps(2,6)=ps(2,6)+luchge(k(i,2))
109  160 CONTINUE
110  pdev=(abs(ps(2,1)-ps(1,1))+abs(ps(2,2)-ps(1,2))+abs(ps(2,3)-
111  &ps(1,3))+abs(ps(2,4)-ps(1,4)))/(1.+abs(ps(2,4))+abs(ps(1,4)))
112  IF(mcons.EQ.1.AND.pdev.GT.paru(11)) CALL luerrm(15,
113  &'(LUEXEC:) four-momentum was not conserved')
114  IF(mcons.EQ.1.AND.abs(ps(2,6)-ps(1,6)).GT.0.1) CALL luerrm(15,
115  &'(LUEXEC:) charge was not conserved')
116 
117  RETURN
118  END