Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
luonia.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file luonia.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE luonia(KFL,ECM)
5 
6 C...Purpose: to generate Upsilon and toponium decays into three
7 C...gluons or two gluons and a photon.
8  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
9  SAVE /lujets/
10  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
11  SAVE /ludat1/
12  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
13  SAVE /ludat2/
14 
15 C...Printout. Check input parameters.
16  IF(mstu(12).GE.1) CALL lulist(0)
17  IF(kfl.LT.0.OR.kfl.GT.8) THEN
18  CALL luerrm(16,'(LUONIA:) called with unknown flavour code')
19  IF(mstu(21).GE.1) RETURN
20  ENDIF
21  IF(ecm.LT.parj(127)+2.02*parf(101)) THEN
22  CALL luerrm(16,'(LUONIA:) called with too small CM energy')
23  IF(mstu(21).GE.1) RETURN
24  ENDIF
25 
26 C...Initial e+e- and onium state (optional).
27  nc=0
28  IF(mstj(115).GE.2) THEN
29  nc=nc+2
30  CALL lu1ent(nc-1,11,0.5*ecm,0.,0.)
31  k(nc-1,1)=21
32  CALL lu1ent(nc,-11,0.5*ecm,paru(1),0.)
33  k(nc,1)=21
34  ENDIF
35  kflc=iabs(kfl)
36  IF(mstj(115).GE.3.AND.kflc.GE.5) THEN
37  nc=nc+1
38  kf=110*kflc+3
39  mstu10=mstu(10)
40  mstu(10)=1
41  p(nc,5)=ecm
42  CALL lu1ent(nc,kf,ecm,0.,0.)
43  k(nc,1)=21
44  k(nc,3)=1
45  mstu(10)=mstu10
46  ENDIF
47 
48 C...Choose x1 and x2 according to matrix element.
49  ntry=0
50  100 x1=rlu(0)
51  x2=rlu(0)
52  x3=2.-x1-x2
53  IF(x3.GE.1..OR.((1.-x1)/(x2*x3))**2+((1.-x2)/(x1*x3))**2+
54  &((1.-x3)/(x1*x2))**2.LE.2.*rlu(0)) goto 100
55  ntry=ntry+1
56  njet=3
57  IF(mstj(101).LE.4) CALL lu3ent(nc+1,21,21,21,ecm,x1,x3)
58  IF(mstj(101).GE.5) CALL lu3ent(-(nc+1),21,21,21,ecm,x1,x3)
59 
60 C...Photon-gluon-gluon events. Small system modifications. Jet origin.
61  mstu(111)=mstj(108)
62  IF(mstj(108).EQ.2.AND.(mstj(101).EQ.0.OR.mstj(101).EQ.1))
63  &mstu(111)=1
64  paru(112)=parj(121)
65  IF(mstu(111).EQ.2) paru(112)=parj(122)
66  qf=0.
67  IF(kflc.NE.0) qf=kchg(kflc,1)/3.
68  rgam=7.2*qf**2*paru(101)/ulalps(ecm**2)
69  mk=0
70  ecmc=ecm
71  IF(rlu(0).GT.rgam/(1.+rgam)) THEN
72  IF(1.-max(x1,x2,x3).LE.max((parj(126)/ecm)**2,parj(125)))
73  & njet=2
74  IF(njet.EQ.2.AND.mstj(101).LE.4) CALL lu2ent(nc+1,21,21,ecm)
75  IF(njet.EQ.2.AND.mstj(101).GE.5) CALL lu2ent(-(nc+1),21,21,ecm)
76  ELSE
77  mk=1
78  ecmc=sqrt(1.-x1)*ecm
79  IF(ecmc.LT.2.*parj(127)) goto 100
80  k(nc+1,1)=1
81  k(nc+1,2)=22
82  k(nc+1,4)=0
83  k(nc+1,5)=0
84  IF(mstj(101).GE.5) k(nc+2,4)=mstu(5)*(nc+3)
85  IF(mstj(101).GE.5) k(nc+2,5)=mstu(5)*(nc+3)
86  IF(mstj(101).GE.5) k(nc+3,4)=mstu(5)*(nc+2)
87  IF(mstj(101).GE.5) k(nc+3,5)=mstu(5)*(nc+2)
88  njet=2
89  IF(ecmc.LT.4.*parj(127)) THEN
90  mstu10=mstu(10)
91  mstu(10)=1
92  p(nc+2,5)=ecmc
93  CALL lu1ent(nc+2,83,0.5*(x2+x3)*ecm,paru(1),0.)
94  mstu(10)=mstu10
95  njet=0
96  ENDIF
97  ENDIF
98  DO 110 ip=nc+1,n
99  110 k(ip,3)=k(ip,3)+(mstj(115)/2)+(kflc/5)*(mstj(115)/3)*(nc-1)
100 
101 C...Differential cross-sections. Upper limit for cross-section.
102  IF(mstj(106).EQ.1) THEN
103  sq2=sqrt(2.)
104  hf1=1.-parj(131)*parj(132)
105  hf3=parj(133)**2
106  ct13=(x1*x3-2.*x1-2.*x3+2.)/(x1*x3)
107  st13=sqrt(1.-ct13**2)
108  sigl=0.5*x3**2*((1.-x2)**2+(1.-x3)**2)*st13**2
109  sigu=(x1*(1.-x1))**2+(x2*(1.-x2))**2+(x3*(1.-x3))**2-sigl
110  sigt=0.5*sigl
111  sigi=(sigl*ct13/st13+0.5*x1*x3*(1.-x2)**2*st13)/sq2
112  sigmax=(2.*hf1+hf3)*abs(sigu)+2.*(hf1+hf3)*abs(sigl)+2.*(hf1+
113  & 2.*hf3)*abs(sigt)+2.*sq2*(hf1+2.*hf3)*abs(sigi)
114 
115 C...Angular orientation of event.
116  120 chi=paru(2)*rlu(0)
117  cthe=2.*rlu(0)-1.
118  phi=paru(2)*rlu(0)
119  cchi=cos(chi)
120  schi=sin(chi)
121  c2chi=cos(2.*chi)
122  s2chi=sin(2.*chi)
123  the=acos(cthe)
124  sthe=sin(the)
125  c2phi=cos(2.*(phi-parj(134)))
126  s2phi=sin(2.*(phi-parj(134)))
127  sig=((1.+cthe**2)*hf1+sthe**2*c2phi*hf3)*sigu+2.*(sthe**2*hf1-
128  & sthe**2*c2phi*hf3)*sigl+2.*(sthe**2*c2chi*hf1+((1.+cthe**2)*
129  & c2chi*c2phi-2.*cthe*s2chi*s2phi)*hf3)*sigt-2.*sq2*(2.*sthe*cthe*
130  & cchi*hf1-2.*sthe*(cthe*cchi*c2phi-schi*s2phi)*hf3)*sigi
131  IF(sig.LT.sigmax*rlu(0)) goto 120
132  CALL ludbrb(nc+1,n,0.,chi,0d0,0d0,0d0)
133  CALL ludbrb(nc+1,n,the,phi,0d0,0d0,0d0)
134  ENDIF
135 
136 C...Generate parton shower. Rearrange along strings and check.
137  IF(mstj(101).GE.5.AND.njet.GE.2) THEN
138  CALL lushow(nc+mk+1,-njet,ecmc)
139  mstj14=mstj(14)
140  IF(mstj(105).EQ.-1) mstj(14)=0
141  IF(mstj(105).GE.0) mstu(28)=0
142  CALL luprep(0)
143  mstj(14)=mstj14
144  IF(mstj(105).GE.0.AND.mstu(28).NE.0) goto 100
145  ENDIF
146 
147 C...Generate fragmentation. Information for LUTABU:
148  IF(mstj(105).EQ.1) CALL luexec
149  mstu(161)=110*kflc+3
150  mstu(162)=0
151 
152  RETURN
153  END