Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pydiff.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pydiff.f
1 
2 C*********************************************************************
3 
4 C...PYDIFF
5 C...Handles diffractive and elastic scattering.
6 
7  SUBROUTINE pydiff
8 
9 C...Double precision and integer declarations.
10  IMPLICIT DOUBLE PRECISION(a-h, o-z)
11  IMPLICIT INTEGER(i-n)
12  INTEGER pyk,pychge,pycomp
13 C...Commonblocks.
14  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
15  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
16  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
17  common/pypars/mstp(200),parp(200),msti(200),pari(200)
18  common/pyint1/mint(400),vint(400)
19  SAVE /pyjets/,/pydat1/,/pypars/,/pyint1/
20  include "mc_set.inc"
21 
22 C...Reset K, P and V vectors. Store incoming particles.
23  DO 110 jt=1,mstp(126)+10
24  i=mint(83)+jt
25  DO 100 j=1,5
26  k(i,j)=0
27  p(i,j)=0d0
28  v(i,j)=0d0
29  100 CONTINUE
30  110 CONTINUE
31  n=mint(84)
32  mint(3)=0
33  mint(21)=0
34  mint(22)=0
35  mint(23)=0
36  mint(24)=0
37  mint(4)=4
38  DO 130 jt=1,2
39  i=mint(83)+jt
40  k(i,1)=21
41  k(i,2)=mint(10+jt)
42  DO 120 j=1,5
43  p(i,j)=vint(285+5*jt+j)
44  120 CONTINUE
45  130 CONTINUE
46  mint(6)=2
47 
48 C...Subprocess; kinematics.
49  sqlam=(vint(2)-vint(63)-vint(64))**2-4d0*vint(63)*vint(64)
50  pz=sqrt(sqlam)/(2d0*vint(1))
51  DO 200 jt=1,2
52  i=mint(83)+jt
53  pe=(vint(2)+vint(62+jt)-vint(65-jt))/(2d0*vint(1))
54  kfh=mint(102+jt)
55 
56 C...Elastically scattered particle. (Except elastic GVMD states.)
57  IF(mint(16+jt).LE.0.AND.(mint(10+jt).NE.22.OR.
58  & mint(106+jt).NE.3)) THEN
59  n=n+1
60  k(n,1)=1
61  k(n,2)=kfh
62  k(n,3)=i+2
63  p(n,3)=pz*(-1)**(jt+1)
64  p(n,4)=pe
65  p(n,5)=sqrt(vint(62+jt))
66 
67 C...Decay rho from elastic scattering of gamma with sin**2(theta)
68 C...distribution of decay products (in rho rest frame).
69  IF(kfh.EQ.113.AND.mint(10+jt).EQ.22.AND.mstp(102).EQ.1) THEN
70  nsav=n
71  dbetaz=p(n,3)/sqrt(p(n,3)**2+p(n,5)**2)
72  p(n,3)=0d0
73  p(n,4)=p(n,5)
74  CALL pydecy(nsav)
75  IF(n.EQ.nsav+2.AND.iabs(k(nsav+1,2)).EQ.211) THEN
76  phi=pyangl(p(nsav+1,1),p(nsav+1,2))
77  CALL pyrobo(nsav+1,nsav+2,0d0,-phi,0d0,0d0,0d0)
78  the=pyangl(p(nsav+1,3),p(nsav+1,1))
79  CALL pyrobo(nsav+1,nsav+2,-the,0d0,0d0,0d0,0d0)
80  140 cthe=2d0*pyr(0)-1d0
81 
82 C... Changing parameters for R_rho with values corresponding to W<7
83 C... (measured by HERMES R_rho=1/eps * r0400/(1. - r0400)
84  pmvirt=pmas(pycomp(113),1)
85  r_rho=parp(165)*(vint(307)/(pmvirt**2))**parp(166)
86  beamas=pymass(11)
87 C new epsilon (f_L/f_T) as used in pysigh.F with proton mass
88  eps=1d0/(1d0+(vint(309)**2*(1d0-2d0*beamas**2/
89  & vint(307)))/(2d0/(1d0+vint(307)/vint(309)**2/
90  & ebeamenucl**2)*(1d0-vint(309)-
91  & (vint(307)/4d0/ebeamenucl**2))))
92  r0400=eps*r_rho / ( 1. + eps * r_rho)
93  w_ang=0.75d0*(1.d0-r0400+(3.d0*r0400-1.d0)*cthe**2.)
94  if( r0400 .le. 1.d0/3.d0 ) then
95  w_ang_max_x = 0.d0
96  else
97  w_ang_max_x = 1.d0
98  endif
99  w_ang_max= 0.75d0*(1.d0-r0400+(3.d0*r0400-1.d0)
100  $ *w_ang_max_x**2.)
101  IF(pyr(0).gt.w_ang/w_ang_max) goto 140
102 C IF(1D0-CTHE**2.LT.PYR(0)) GOTO 140
103  CALL pyrobo(nsav+1,nsav+2,acos(cthe),phi,0d0,0d0,0d0)
104  ENDIF
105  CALL pyrobo(nsav,nsav+2,0d0,0d0,0d0,0d0,dbetaz)
106  ENDIF
107 
108 C...Diffracted particle: low-mass system to two particles.
109  ELSEIF(vint(62+jt).LT.(vint(66+jt)+parp(103))**2) THEN
110  n=n+2
111  k(n-1,1)=1
112  k(n,1)=1
113  k(n-1,3)=i+2
114  k(n,3)=i+2
115  pmmas=sqrt(vint(62+jt))
116  ntry=0
117  150 ntry=ntry+1
118  IF(ntry.LT.20) THEN
119  mint(105)=mint(102+jt)
120  mint(109)=mint(106+jt)
121  CALL pyspli(kfh,21,kfl1,kfl2)
122  CALL pykfdi(kfl1,0,kfl3,kf1)
123  IF(kf1.EQ.0) goto 150
124  CALL pykfdi(kfl2,-kfl3,kfldum,kf2)
125  IF(kf2.EQ.0) goto 150
126  ELSE
127  kf1=kfh
128  kf2=111
129  ENDIF
130  pm1=pymass(kf1)
131  pm2=pymass(kf2)
132  IF(pm1+pm2+parj(64).GT.pmmas) goto 150
133  k(n-1,2)=kf1
134  k(n,2)=kf2
135  p(n-1,5)=pm1
136  p(n,5)=pm2
137  pzp=sqrt(max(0d0,(pmmas**2-pm1**2-pm2**2)**2-
138  & 4d0*pm1**2*pm2**2))/(2d0*pmmas)
139  p(n-1,3)=pzp
140  p(n,3)=-pzp
141  p(n-1,4)=sqrt(pm1**2+pzp**2)
142  p(n,4)=sqrt(pm2**2+pzp**2)
143  CALL pyrobo(n-1,n,acos(2d0*pyr(0)-1d0),paru(2)*pyr(0),
144  & 0d0,0d0,0d0)
145  dbetaz=pz*(-1)**(jt+1)/sqrt(pz**2+pmmas**2)
146  CALL pyrobo(n-1,n,0d0,0d0,0d0,0d0,dbetaz)
147 
148 C...Diffracted particle: valence quark kicked out.
149  ELSEIF(mstp(101).EQ.1.OR.(mstp(101).EQ.3.AND.pyr(0).LT.
150  & parp(101))) THEN
151  n=n+2
152  k(n-1,1)=2
153  k(n,1)=1
154  k(n-1,3)=i+2
155  k(n,3)=i+2
156  mint(105)=mint(102+jt)
157  mint(109)=mint(106+jt)
158  CALL pyspli(kfh,21,k(n,2),k(n-1,2))
159  p(n-1,5)=pymass(k(n-1,2))
160  p(n,5)=pymass(k(n,2))
161  sqlam=(vint(62+jt)-p(n-1,5)**2-p(n,5)**2)**2-
162  & 4d0*p(n-1,5)**2*p(n,5)**2
163  p(n-1,3)=(pe*sqrt(sqlam)+pz*(vint(62+jt)+p(n-1,5)**2-
164  & p(n,5)**2))/(2d0*vint(62+jt))*(-1)**(jt+1)
165  p(n-1,4)=sqrt(p(n-1,3)**2+p(n-1,5)**2)
166  p(n,3)=pz*(-1)**(jt+1)-p(n-1,3)
167  p(n,4)=sqrt(p(n,3)**2+p(n,5)**2)
168 
169 C...Diffracted particle: gluon kicked out.
170  ELSE
171  n=n+3
172  k(n-2,1)=2
173  k(n-1,1)=2
174  k(n,1)=1
175  k(n-2,3)=i+2
176  k(n-1,3)=i+2
177  k(n,3)=i+2
178  mint(105)=mint(102+jt)
179  mint(109)=mint(106+jt)
180  CALL pyspli(kfh,21,k(n,2),k(n-2,2))
181  k(n-1,2)=21
182  p(n-2,5)=pymass(k(n-2,2))
183  p(n-1,5)=0d0
184  p(n,5)=pymass(k(n,2))
185 C...Energy distribution for particle into two jets.
186  160 imb=1
187  IF(mod(kfh/1000,10).NE.0) imb=2
188  chik=parp(92+2*imb)
189  IF(mstp(92).LE.1) THEN
190  IF(imb.EQ.1) chi=pyr(0)
191  IF(imb.EQ.2) chi=1d0-sqrt(pyr(0))
192  ELSEIF(mstp(92).EQ.2) THEN
193  chi=1d0-pyr(0)**(1d0/(1d0+chik))
194  ELSEIF(mstp(92).EQ.3) THEN
195  cut=2d0*0.3d0/vint(1)
196  170 chi=pyr(0)**2
197  IF((chi**2/(chi**2+cut**2))**0.25d0*(1d0-chi)**chik.LT.
198  & pyr(0)) goto 170
199  ELSEIF(mstp(92).EQ.4) THEN
200  cut=2d0*0.3d0/vint(1)
201  cutr=(1d0+sqrt(1d0+cut**2))/cut
202  180 chir=cut*cutr**pyr(0)
203  chi=(chir**2-cut**2)/(2d0*chir)
204  IF((1d0-chi)**chik.LT.pyr(0)) goto 180
205  ELSE
206  cut=2d0*0.3d0/vint(1)
207  cuta=cut**(1d0-parp(98))
208  cutb=(1d0+cut)**(1d0-parp(98))
209  190 chi=(cuta+pyr(0)*(cutb-cuta))**(1d0/(1d0-parp(98)))
210  IF(((chi+cut)**2/(2d0*(chi**2+cut**2)))**
211  & (0.5d0*parp(98))*(1d0-chi)**chik.LT.pyr(0)) goto 190
212  ENDIF
213  IF(chi.LT.p(n,5)**2/vint(62+jt).OR.chi.GT.1d0-p(n-2,5)**2/
214  & vint(62+jt)) goto 160
215  sqm=p(n-2,5)**2/(1d0-chi)+p(n,5)**2/chi
216  pzi=(pe*(vint(62+jt)-sqm)+pz*(vint(62+jt)+sqm))/
217  & (2d0*vint(62+jt))
218  pei=sqrt(pzi**2+sqm)
219  pqqp=(1d0-chi)*(pei+pzi)
220  p(n-2,3)=0.5d0*(pqqp-p(n-2,5)**2/pqqp)*(-1)**(jt+1)
221  p(n-2,4)=sqrt(p(n-2,3)**2+p(n-2,5)**2)
222  p(n-1,4)=0.5d0*(vint(62+jt)-sqm)/(pei+pzi)
223  p(n-1,3)=p(n-1,4)*(-1)**jt
224  p(n,3)=pzi*(-1)**(jt+1)-p(n-2,3)
225  p(n,4)=sqrt(p(n,3)**2+p(n,5)**2)
226  ENDIF
227 
228 C...Documentation lines.
229  k(i+2,1)=21
230  IF(mint(16+jt).EQ.0) k(i+2,2)=kfh
231  IF(mint(16+jt).NE.0.OR.(mint(10+jt).EQ.22.AND.
232  & mint(106+jt).EQ.3)) k(i+2,2)=isign(9900000,kfh)+10*(kfh/10)
233  k(i+2,3)=i
234  p(i+2,3)=pz*(-1)**(jt+1)
235  p(i+2,4)=pe
236  p(i+2,5)=sqrt(vint(62+jt))
237  200 CONTINUE
238 
239 C...Rotate outgoing partons/particles using cos(theta).
240  IF(vint(23).LT.0.9d0) THEN
241  CALL pyrobo(mint(83)+3,n,acos(vint(23)),vint(24),0d0,0d0,0d0)
242  ELSE
243  CALL pyrobo(mint(83)+3,n,asin(vint(59)),vint(24),0d0,0d0,0d0)
244  ENDIF
245 
246  RETURN
247  END