ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyhidiff.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhidiff.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE pyhidiff
5 
6 C...Handles diffractive and elastic scattering.
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/pyhipars/mstp(200),parp(200),msti(200),pari(200)
12  SAVE /pyhipars/
13  common/pyhiint1/mint(400),vint(400)
14  SAVE /pyhiint1/
15 
16 C...Reset K, P and V vectors. Store incoming particles.
17  DO 100 jt=1,mstp(126)+10
18  i=mint(83)+jt
19  DO 100 j=1,5
20  k(i,j)=0
21  p(i,j)=0.
22  100 v(i,j)=0.
23  n=mint(84)
24  mint(3)=0
25  mint(21)=0
26  mint(22)=0
27  mint(23)=0
28  mint(24)=0
29  mint(4)=4
30  DO 110 jt=1,2
31  i=mint(83)+jt
32  k(i,1)=21
33  k(i,2)=mint(10+jt)
34  p(i,5)=vint(2+jt)
35  p(i,3)=vint(5)*(-1)**(jt+1)
36  110 p(i,4)=sqrt(p(i,3)**2+p(i,5)**2)
37  mint(6)=2
38 
39 C...Subprocess; kinematics.
40  isub=mint(1)
41  sqlam=(vint(2)-vint(63)-vint(64))**2-4.*vint(63)*vint(64)
42  pz=sqrt(sqlam)/(2.*vint(1))
43  DO 150 jt=1,2
44  i=mint(83)+jt
45  pe=(vint(2)+vint(62+jt)-vint(65-jt))/(2.*vint(1))
46 
47 C...Elastically scattered particle.
48  IF(mint(16+jt).LE.0) THEN
49  n=n+1
50  k(n,1)=1
51  k(n,2)=k(i,2)
52  k(n,3)=i+2
53  p(n,3)=pz*(-1)**(jt+1)
54  p(n,4)=pe
55  p(n,5)=p(i,5)
56 
57 C...Diffracted particle: valence quark kicked out.
58  ELSEIF(mstp(101).EQ.1) THEN
59  n=n+2
60  k(n-1,1)=2
61  k(n,1)=1
62  k(n-1,3)=i+2
63  k(n,3)=i+2
64  CALL pyhispli(k(i,2),21,k(n,2),k(n-1,2))
65  p(n-1,5)=ulmass(k(n-1,2))
66  p(n,5)=ulmass(k(n,2))
67  sqlam=(vint(62+jt)-p(n-1,5)**2-p(n,5)**2)**2-
68  & 4.*p(n-1,5)**2*p(n,5)**2
69  p(n-1,3)=(pe*sqrt(sqlam)+pz*(vint(62+jt)+p(n-1,5)**2-
70  & p(n,5)**2))/(2.*vint(62+jt))*(-1)**(jt+1)
71  p(n-1,4)=sqrt(p(n-1,3)**2+p(n-1,5)**2)
72  p(n,3)=pz*(-1)**(jt+1)-p(n-1,3)
73  p(n,4)=sqrt(p(n,3)**2+p(n,5)**2)
74 
75 C...Diffracted particle: gluon kicked out.
76  ELSE
77  n=n+3
78  k(n-2,1)=2
79  k(n-1,1)=2
80  k(n,1)=1
81  k(n-2,3)=i+2
82  k(n-1,3)=i+2
83  k(n,3)=i+2
84  CALL pyhispli(k(i,2),21,k(n,2),k(n-2,2))
85  k(n-1,2)=21
86  p(n-2,5)=ulmass(k(n-2,2))
87  p(n-1,5)=0.
88  p(n,5)=ulmass(k(n,2))
89 C...Energy distribution for particle into two jets.
90  120 imb=1
91  IF(mod(k(i,2)/1000,10).NE.0) imb=2
92  chik=parp(92+2*imb)
93  IF(mstp(92).LE.1) THEN
94  IF(imb.EQ.1) chi=rlu(0)
95  IF(imb.EQ.2) chi=1.-sqrt(rlu(0))
96  ELSEIF(mstp(92).EQ.2) THEN
97  chi=1.-rlu(0)**(1./(1.+chik))
98  ELSEIF(mstp(92).EQ.3) THEN
99  cut=2.*0.3/vint(1)
100  130 chi=rlu(0)**2
101  IF((chi**2/(chi**2+cut**2))**0.25*(1.-chi)**chik.LT.
102  & rlu(0)) goto 130
103  ELSE
104  cut=2.*0.3/vint(1)
105  cutr=(1.+sqrt(1.+cut**2))/cut
106  140 chir=cut*cutr**rlu(0)
107  chi=(chir**2-cut**2)/(2.*chir)
108  IF((1.-chi)**chik.LT.rlu(0)) goto 140
109  ENDIF
110  IF(chi.LT.p(n,5)**2/vint(62+jt).OR.chi.GT.1.-p(n-2,5)**2/
111  & vint(62+jt)) goto 120
112  sqm=p(n-2,5)**2/(1.-chi)+p(n,5)**2/chi
113  IF((sqrt(sqm)+parj(32))**2.GE.vint(62+jt)) goto 120
114  pzi=(pe*(vint(62+jt)-sqm)+pz*(vint(62+jt)+sqm))/
115  & (2.*vint(62+jt))
116  pei=sqrt(pzi**2+sqm)
117  pqqp=(1.-chi)*(pei+pzi)
118  p(n-2,3)=0.5*(pqqp-p(n-2,5)**2/pqqp)*(-1)**(jt+1)
119  p(n-2,4)=sqrt(p(n-2,3)**2+p(n-2,5)**2)
120  p(n-1,3)=(pz-pzi)*(-1)**(jt+1)
121  p(n-1,4)=abs(p(n-1,3))
122  p(n,3)=pzi*(-1)**(jt+1)-p(n-2,3)
123  p(n,4)=sqrt(p(n,3)**2+p(n,5)**2)
124  ENDIF
125 
126 C...Documentation lines.
127  k(i+2,1)=21
128  IF(mint(16+jt).EQ.0) k(i+2,2)=mint(10+jt)
129  IF(mint(16+jt).NE.0) k(i+2,2)=10*(mint(10+jt)/10)
130  k(i+2,3)=i
131  p(i+2,3)=pz*(-1)**(jt+1)
132  p(i+2,4)=pe
133  p(i+2,5)=sqrt(vint(62+jt))
134  150 CONTINUE
135 
136 C...Rotate outgoing partons/particles using cos(theta).
137  CALL ludbrb(mint(83)+3,n,acos(vint(23)),vint(24),0d0,0d0,0d0)
138 
139  RETURN
140  END