Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
luxtot.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file luxtot.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE luxtot(KFL,ECM,XTOT)
5 
6 C...Purpose: to calculate total cross-section, including initial
7 C...state radiation effects.
8  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
9  SAVE /ludat1/
10  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
11  SAVE /ludat2/
12 
13 C...Status, (optimized) Q^2 scale, alpha_strong.
14  parj(151)=ecm
15  mstj(119)=10*mstj(102)+kfl
16  IF(mstj(111).EQ.0) THEN
17  q2r=ecm**2
18  ELSEIF(mstu(111).EQ.0) THEN
19  parj(168)=min(1.,max(parj(128),exp(-12.*paru(1)/
20  & ((33.-2.*mstu(112))*paru(111)))))
21  q2r=parj(168)*ecm**2
22  ELSE
23  parj(168)=min(1.,max(parj(128),paru(112)/ecm,
24  & (2.*paru(112)/ecm)**2))
25  q2r=parj(168)*ecm**2
26  ENDIF
27  alspi=ulalps(q2r)/paru(1)
28 
29 C...QCD corrections factor in R.
30  IF(mstj(101).EQ.0.OR.mstj(109).EQ.1) THEN
31  rqcd=1.
32  ELSEIF(iabs(mstj(101)).EQ.1.AND.mstj(109).EQ.0) THEN
33  rqcd=1.+alspi
34  ELSEIF(mstj(109).EQ.0) THEN
35  rqcd=1.+alspi+(1.986-0.115*mstu(118))*alspi**2
36  IF(mstj(111).EQ.1) rqcd=max(1.,rqcd+(33.-2.*mstu(112))/12.*
37  & log(parj(168))*alspi**2)
38  ELSEIF(iabs(mstj(101)).EQ.1) THEN
39  rqcd=1.+(3./4.)*alspi
40  ELSE
41  rqcd=1.+(3./4.)*alspi-(3./32.+0.519*mstu(118))*alspi**2
42  ENDIF
43 
44 C...Calculate Z0 width if default value not acceptable.
45  IF(mstj(102).GE.3) THEN
46  rva=3.*(3.+(4.*paru(102)-1.)**2)+6.*rqcd*(2.+(1.-8.*paru(102)/
47  & 3.)**2+(4.*paru(102)/3.-1.)**2)
48  DO 100 kflc=5,6
49  vq=1.
50  IF(mod(mstj(103),2).EQ.1) vq=sqrt(max(0.,1.-(2.*ulmass(kflc)/
51  & ecm)**2))
52  IF(kflc.EQ.5) vf=4.*paru(102)/3.-1.
53  IF(kflc.EQ.6) vf=1.-8.*paru(102)/3.
54  100 rva=rva+3.*rqcd*(0.5*vq*(3.-vq**2)*vf**2+vq**3)
55  parj(124)=paru(101)*parj(123)*rva/(48.*paru(102)*(1.-paru(102)))
56  ENDIF
57 
58 C...Calculate propagator and related constants for QFD case.
59  poll=1.-parj(131)*parj(132)
60  IF(mstj(102).GE.2) THEN
61  sff=1./(16.*paru(102)*(1.-paru(102)))
62  sfw=ecm**4/((ecm**2-parj(123)**2)**2+(parj(123)*parj(124))**2)
63  sfi=sfw*(1.-(parj(123)/ecm)**2)
64  ve=4.*paru(102)-1.
65  sf1i=sff*(ve*poll+parj(132)-parj(131))
66  sf1w=sff**2*((ve**2+1.)*poll+2.*ve*(parj(132)-parj(131)))
67  hf1i=sfi*sf1i
68  hf1w=sfw*sf1w
69  ENDIF
70 
71 C...Loop over different flavours: charge, velocity.
72  rtot=0.
73  rqq=0.
74  rqv=0.
75  rva=0.
76  DO 110 kflc=1,max(mstj(104),kfl)
77  IF(kfl.GT.0.AND.kflc.NE.kfl) goto 110
78  mstj(93)=1
79  pmq=ulmass(kflc)
80  IF(ecm.LT.2.*pmq+parj(127)) goto 110
81  qf=kchg(kflc,1)/3.
82  vq=1.
83  IF(mod(mstj(103),2).EQ.1) vq=sqrt(1.-(2.*pmq/ecm)**2)
84 
85 C...Calculate R and sum of charges for QED or QFD case.
86  rqq=rqq+3.*qf**2*poll
87  IF(mstj(102).LE.1) THEN
88  rtot=rtot+3.*0.5*vq*(3.-vq**2)*qf**2*poll
89  ELSE
90  vf=sign(1.,qf)-4.*qf*paru(102)
91  rqv=rqv-6.*qf*vf*sf1i
92  rva=rva+3.*(vf**2+1.)*sf1w
93  rtot=rtot+3.*(0.5*vq*(3.-vq**2)*(qf**2*poll-2.*qf*vf*hf1i+
94  & vf**2*hf1w)+vq**3*hf1w)
95  ENDIF
96  110 CONTINUE
97  rsum=rqq
98  IF(mstj(102).GE.2) rsum=rqq+sfi*rqv+sfw*rva
99 
100 C...Calculate cross-section, including QCD corrections.
101  parj(141)=rqq
102  parj(142)=rtot
103  parj(143)=rtot*rqcd
104  parj(144)=parj(143)
105  parj(145)=parj(141)*86.8/ecm**2
106  parj(146)=parj(142)*86.8/ecm**2
107  parj(147)=parj(143)*86.8/ecm**2
108  parj(148)=parj(147)
109  parj(157)=rsum*rqcd
110  parj(158)=0.
111  parj(159)=0.
112  xtot=parj(147)
113  IF(mstj(107).LE.0) RETURN
114 
115 C...Virtual cross-section.
116  xkl=parj(135)
117  xku=min(parj(136),1.-(2.*parj(127)/ecm)**2)
118  ale=2.*log(ecm/ulmass(11))-1.
119  sigv=ale/3.+2.*log(ecm**2/(ulmass(13)*ulmass(15)))/3.-4./3.+
120  &1.526*log(ecm**2/0.932)
121 
122 C...Soft and hard radiative cross-section in QED case.
123  IF(mstj(102).LE.1) THEN
124  sigv=1.5*ale-0.5+paru(1)**2/3.+2.*sigv
125  sigs=ale*(2.*log(xkl)-log(1.-xkl)-xkl)
126  sigh=ale*(2.*log(xku/xkl)-log((1.-xku)/(1.-xkl))-(xku-xkl))
127 
128 C...Soft and hard radiative cross-section in QFD case.
129  ELSE
130  szm=1.-(parj(123)/ecm)**2
131  szw=parj(123)*parj(124)/ecm**2
132  parj(161)=-rqq/rsum
133  parj(162)=-(rqq+rqv+rva)/rsum
134  parj(163)=(rqv*(1.-0.5*szm-sfi)+rva*(1.5-szm-sfw))/rsum
135  parj(164)=(rqv*szw**2*(1.-2.*sfw)+rva*(2.*sfi+szw**2-4.+3.*szm-
136  & szm**2))/(szw*rsum)
137  sigv=1.5*ale-0.5+paru(1)**2/3.+((2.*rqq+sfi*rqv)/rsum)*sigv+
138  & (szw*sfw*rqv/rsum)*paru(1)*20./9.
139  sigs=ale*(2.*log(xkl)+parj(161)*log(1.-xkl)+parj(162)*xkl+
140  & parj(163)*log(((xkl-szm)**2+szw**2)/(szm**2+szw**2))+
141  & parj(164)*(atan((xkl-szm)/szw)-atan(-szm/szw)))
142  sigh=ale*(2.*log(xku/xkl)+parj(161)*log((1.-xku)/(1.-xkl))+
143  & parj(162)*(xku-xkl)+parj(163)*log(((xku-szm)**2+szw**2)/
144  & ((xkl-szm)**2+szw**2))+parj(164)*(atan((xku-szm)/szw)-
145  & atan((xkl-szm)/szw)))
146  ENDIF
147 
148 C...Total cross-section and fraction of hard photon events.
149  parj(160)=sigh/(paru(1)/paru(101)+sigv+sigs+sigh)
150  parj(157)=rsum*(1.+(paru(101)/paru(1))*(sigv+sigs+sigh))*rqcd
151  parj(144)=parj(157)
152  parj(148)=parj(144)*86.8/ecm**2
153  xtot=parj(148)
154 
155  RETURN
156  END