EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
luxdif.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file luxdif.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE luxdif(NC,NJET,KFL,ECM,CHI,THE,PHI)
5 
6 C...Purpose: to give the angular orientation of events.
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 
14 C...Charge. Factors depending on polarization for QED case.
15  qf=kchg(kfl,1)/3.
16  poll=1.-parj(131)*parj(132)
17  pold=parj(132)-parj(131)
18  IF(mstj(102).LE.1.OR.mstj(109).EQ.1) THEN
19  hf1=poll
20  hf2=0.
21  hf3=parj(133)**2
22  hf4=0.
23 
24 C...Factors depending on flavour, energy and polarization for QFD case.
25  ELSE
26  sff=1./(16.*paru(102)*(1.-paru(102)))
27  sfw=ecm**4/((ecm**2-parj(123)**2)**2+(parj(123)*parj(124))**2)
28  sfi=sfw*(1.-(parj(123)/ecm)**2)
29  ae=-1.
30  ve=4.*paru(102)-1.
31  af=sign(1.,qf)
32  vf=af-4.*qf*paru(102)
33  hf1=qf**2*poll-2.*qf*vf*sfi*sff*(ve*poll-ae*pold)+
34  & (vf**2+af**2)*sfw*sff**2*((ve**2+ae**2)*poll-2.*ve*ae*pold)
35  hf2=-2.*qf*af*sfi*sff*(ae*poll-ve*pold)+2.*vf*af*sfw*sff**2*
36  & (2.*ve*ae*poll-(ve**2+ae**2)*pold)
37  hf3=parj(133)**2*(qf**2-2.*qf*vf*sfi*sff*ve+(vf**2+af**2)*
38  & sfw*sff**2*(ve**2-ae**2))
39  hf4=-parj(133)**2*2.*qf*vf*sfw*(parj(123)*parj(124)/ecm**2)*
40  & sff*ae
41  ENDIF
42 
43 C...Mass factor. Differential cross-sections for two-jet events.
44  sq2=sqrt(2.)
45  qme=0.
46  IF(mstj(103).GE.4.AND.iabs(mstj(101)).LE.1.AND.mstj(102).LE.1.AND.
47  &mstj(109).NE.1) qme=(2.*ulmass(kfl)/ecm)**2
48  IF(njet.EQ.2) THEN
49  sigu=4.*sqrt(1.-qme)
50  sigl=2.*qme*sqrt(1.-qme)
51  sigt=0.
52  sigi=0.
53  siga=0.
54  sigp=4.
55 
56 C...Kinematical variables. Reduce four-jet event to three-jet one.
57  ELSE
58  IF(njet.EQ.3) THEN
59  x1=2.*p(nc+1,4)/ecm
60  x2=2.*p(nc+3,4)/ecm
61  ELSE
62  ecmr=p(nc+1,4)+p(nc+4,4)+sqrt((p(nc+2,1)+p(nc+3,1))**2+
63  & (p(nc+2,2)+p(nc+3,2))**2+(p(nc+2,3)+p(nc+3,3))**2)
64  x1=2.*p(nc+1,4)/ecmr
65  x2=2.*p(nc+4,4)/ecmr
66  ENDIF
67 
68 C...Differential cross-sections for three-jet (or reduced four-jet).
69  xq=(1.-x1)/(1.-x2)
70  ct12=(x1*x2-2.*x1-2.*x2+2.+qme)/sqrt((x1**2-qme)*(x2**2-qme))
71  st12=sqrt(1.-ct12**2)
72  IF(mstj(109).NE.1) THEN
73  sigu=2.*x1**2+x2**2*(1.+ct12**2)-qme*(3.+ct12**2-x1-x2)-
74  & qme*x1/xq+0.5*qme*((x2**2-qme)*st12**2-2.*x2)*xq
75  sigl=(x2*st12)**2-qme*(3.-ct12**2-2.5*(x1+x2)+x1*x2+qme)+
76  & 0.5*qme*(x1**2-x1-qme)/xq+0.5*qme*((x2**2-qme)*ct12**2-x2)*xq
77  sigt=0.5*(x2**2-qme-0.5*qme*(x2**2-qme)/xq)*st12**2
78  sigi=((1.-0.5*qme*xq)*(x2**2-qme)*st12*ct12+qme*(1.-x1-x2+
79  & 0.5*x1*x2+0.5*qme)*st12/ct12)/sq2
80  siga=x2**2*st12/sq2
81  sigp=2.*(x1**2-x2**2*ct12)
82 
83 C...Differential cross-sect for scalar gluons (no mass or QFD effects).
84  ELSE
85  sigu=2.*(2.-x1-x2)**2-(x2*st12)**2
86  sigl=(x2*st12)**2
87  sigt=0.5*sigl
88  sigi=-(2.-x1-x2)*x2*st12/sq2
89  siga=0.
90  sigp=0.
91  ENDIF
92  ENDIF
93 
94 C...Upper bounds for differential cross-section.
95  hf1a=abs(hf1)
96  hf2a=abs(hf2)
97  hf3a=abs(hf3)
98  hf4a=abs(hf4)
99  sigmax=(2.*hf1a+hf3a+hf4a)*abs(sigu)+2.*(hf1a+hf3a+hf4a)*
100  &abs(sigl)+2.*(hf1a+2.*hf3a+2.*hf4a)*abs(sigt)+2.*sq2*
101  &(hf1a+2.*hf3a+2.*hf4a)*abs(sigi)+4.*sq2*hf2a*abs(siga)+
102  &2.*hf2a*abs(sigp)
103 
104 C...Generate angular orientation according to differential cross-sect.
105  100 chi=paru(2)*rlu(0)
106  cthe=2.*rlu(0)-1.
107  phi=paru(2)*rlu(0)
108  cchi=cos(chi)
109  schi=sin(chi)
110  c2chi=cos(2.*chi)
111  s2chi=sin(2.*chi)
112  the=acos(cthe)
113  sthe=sin(the)
114  c2phi=cos(2.*(phi-parj(134)))
115  s2phi=sin(2.*(phi-parj(134)))
116  sig=((1.+cthe**2)*hf1+sthe**2*(c2phi*hf3-s2phi*hf4))*sigu+
117  &2.*(sthe**2*hf1-sthe**2*(c2phi*hf3-s2phi*hf4))*sigl+
118  &2.*(sthe**2*c2chi*hf1+((1.+cthe**2)*c2chi*c2phi-2.*cthe*s2chi*
119  &s2phi)*hf3-((1.+cthe**2)*c2chi*s2phi+2.*cthe*s2chi*c2phi)*hf4)*
120  &sigt-2.*sq2*(2.*sthe*cthe*cchi*hf1-2.*sthe*(cthe*cchi*c2phi-
121  &schi*s2phi)*hf3+2.*sthe*(cthe*cchi*s2phi+schi*c2phi)*hf4)*sigi+
122  &4.*sq2*sthe*cchi*hf2*siga+2.*cthe*hf2*sigp
123  IF(sig.LT.sigmax*rlu(0)) goto 100
124 
125  RETURN
126  END