Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pydisg.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pydisg.f
1 
2 C*********************************************************************
3 
4 C...PYDISG
5 C...Set up a DIS process as gamma* + f -> f, with beam remnant
6 C...and showering added consecutively. Photon flux by the PYGAGA
7 C...routine (if at all).
8 
9  SUBROUTINE pydisg
10 
11 C...Double precision and integer declarations.
12  IMPLICIT DOUBLE PRECISION(a-h, o-z)
13  IMPLICIT INTEGER(i-n)
14  INTEGER pyk,pychge,pycomp
15 C...Parameter statement to help give large particle numbers.
16  parameter(ksusy1=1000000,ksusy2=2000000,ktechn=3000000,
17  &kexcit=4000000,kdimen=5000000)
18 C...Commonblocks.
19  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
20  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
21  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
22  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
23  common/pypars/mstp(200),parp(200),msti(200),pari(200)
24  common/pyint1/mint(400),vint(400)
25  SAVE /pyjets/,/pydat1/,/pydat2/,/pysubs/,/pypars/,/pyint1/
26 C...Local arrays.
27  dimension pms(4)
28 
29 C...Choice of subprocess, number of documentation lines
30  idoc=7
31  mint(3)=idoc-6
32  mint(4)=idoc
33  ipu1=mint(84)+1
34  ipu2=mint(84)+2
35  ipu3=mint(84)+3
36  iside=1
37  IF(mint(107).EQ.4) iside=2
38 
39 C...Reset K, P and V vectors. Store incoming particles
40  DO 110 jt=1,mstp(126)+20
41  i=mint(83)+jt
42  DO 100 j=1,5
43  k(i,j)=0
44  p(i,j)=0d0
45  v(i,j)=0d0
46  100 CONTINUE
47  110 CONTINUE
48  DO 130 jt=1,2
49  i=mint(83)+jt
50  k(i,1)=21
51  k(i,2)=mint(10+jt)
52  DO 120 j=1,5
53  p(i,j)=vint(285+5*jt+j)
54  120 CONTINUE
55  130 CONTINUE
56  mint(6)=2
57 
58 C...Store incoming partons in hadronic CM-frame
59  DO 140 jt=1,2
60  i=mint(84)+jt
61  k(i,1)=14
62  k(i,2)=mint(14+jt)
63  k(i,3)=mint(83)+2+jt
64  140 CONTINUE
65  IF(mint(15).EQ.22) THEN
66  p(mint(84)+1,3)=0.5d0*(vint(1)+vint(307)/vint(1))
67  p(mint(84)+1,4)=0.5d0*(vint(1)-vint(307)/vint(1))
68  p(mint(84)+1,5)=-sqrt(vint(307))
69  p(mint(84)+2,3)=-0.5d0*vint(307)/vint(1)
70  p(mint(84)+2,4)=0.5d0*vint(307)/vint(1)
71  kfres=mint(16)
72  iside=2
73  ELSE
74  p(mint(84)+1,3)=0.5d0*vint(308)/vint(1)
75  p(mint(84)+1,4)=0.5d0*vint(308)/vint(1)
76  p(mint(84)+2,3)=-0.5d0*(vint(1)+vint(308)/vint(1))
77  p(mint(84)+2,4)=0.5d0*(vint(1)-vint(308)/vint(1))
78  p(mint(84)+1,5)=-sqrt(vint(308))
79  kfres=mint(15)
80  iside=1
81  ENDIF
82  sidesg=(-1d0)**(iside-1)
83 
84 C...Copy incoming partons to documentation lines.
85  DO 170 jt=1,2
86  i1=mint(83)+4+jt
87  i2=mint(84)+jt
88  k(i1,1)=21
89  k(i1,2)=k(i2,2)
90  k(i1,3)=i1-2
91  DO 150 j=1,5
92  p(i1,j)=p(i2,j)
93  150 CONTINUE
94 
95 C...Second copy for partons before ISR shower, since no such.
96  i1=mint(83)+2+jt
97  k(i1,1)=21
98  k(i1,2)=k(i2,2)
99  k(i1,3)=i1-2
100  DO 160 j=1,5
101  p(i1,j)=p(i2,j)
102  160 CONTINUE
103  170 CONTINUE
104 
105 C...Define initial partons.
106  ntry=0
107  180 ntry=ntry+1
108  IF(ntry.GT.100) THEN
109  mint(51)=1
110  RETURN
111  ENDIF
112 
113 C...Scattered quark in hadronic CM frame.
114  i=mint(83)+7
115  k(ipu3,1)=3
116  k(ipu3,2)=kfres
117  k(ipu3,3)=i
118  p(ipu3,5)=pymass(kfres)
119  p(ipu3,3)=p(ipu1,3)+p(ipu2,3)
120  p(ipu3,4)=p(ipu1,4)+p(ipu2,4)
121  p(ipu3,5)=0d0
122  k(i,1)=21
123  k(i,2)=kfres
124  k(i,3)=mint(83)+4+iside
125  p(i,3)=p(ipu3,3)
126  p(i,4)=p(ipu3,4)
127  p(i,5)=p(ipu3,5)
128  n=ipu3
129  mint(21)=kfres
130  mint(22)=0
131 
132 C...No primordial kT, or chosen according to truncated Gaussian or
133 C...exponential, or (for photon) predetermined or power law.
134  190 IF(mint(40+iside).EQ.2.AND.mint(10+iside).NE.22) THEN
135  IF(mstp(91).LE.0) THEN
136  pt=0d0
137  ELSEIF(mstp(91).EQ.1) THEN
138  pt=parp(91)*sqrt(-log(pyr(0)))
139  ELSE
140  rpt1=pyr(0)
141  rpt2=pyr(0)
142  pt=-parp(92)*log(rpt1*rpt2)
143  ENDIF
144  IF(pt.GT.parp(93)) goto 190
145  ELSEIF(mint(106+iside).EQ.3) THEN
146  pta=sqrt(vint(282+iside))
147  ptb=0d0
148  IF(mstp(66).EQ.5.AND.mstp(93).EQ.1) THEN
149  ptb=parp(99)*sqrt(-log(pyr(0)))
150  ELSEIF(mstp(66).EQ.5.AND.mstp(93).EQ.2) THEN
151  rpt1=pyr(0)
152  rpt2=pyr(0)
153  ptb=-parp(99)*log(rpt1*rpt2)
154  ENDIF
155  IF(ptb.GT.parp(100)) goto 190
156  pt=sqrt(pta**2+ptb**2+2d0*pta*ptb*cos(paru(2)*pyr(0)))
157  IF(ntry.GT.10) pt=pt*0.8d0**(ntry-10)
158  ELSEIF(iabs(mint(14+iside)).LE.8.OR.mint(14+iside).EQ.21) THEN
159  IF(mstp(93).LE.0) THEN
160  pt=0d0
161  ELSEIF(mstp(93).EQ.1) THEN
162  pt=parp(99)*sqrt(-log(pyr(0)))
163  ELSEIF(mstp(93).EQ.2) THEN
164  rpt1=pyr(0)
165  rpt2=pyr(0)
166  pt=-parp(99)*log(rpt1*rpt2)
167  ELSEIF(mstp(93).EQ.3) THEN
168  ha=parp(99)**2
169  hb=parp(100)**2
170  pt=sqrt(max(0d0,ha*(ha+hb)/(ha+hb-pyr(0)*hb)-ha))
171  ELSE
172  ha=parp(99)**2
173  hb=parp(100)**2
174  IF(mstp(93).EQ.5) hb=min(vint(48),parp(100)**2)
175  pt=sqrt(max(0d0,ha*((ha+hb)/ha)**pyr(0)-ha))
176  ENDIF
177  IF(pt.GT.parp(100)) goto 190
178  ELSE
179  pt=0d0
180  ENDIF
181  vint(156+iside)=pt
182  phi=paru(2)*pyr(0)
183  p(ipu3,1)=pt*cos(phi)
184  p(ipu3,2)=pt*sin(phi)
185  p(ipu3,4)=sqrt(p(ipu3,5)**2+pt**2+p(ipu3,3)**2)
186  pms(3-iside)=p(ipu3,5)**2+p(ipu3,1)**2+p(ipu3,2)**2
187  pcp=p(ipu3,4)+abs(p(ipu3,3))
188 
189 C...Find one or two beam remnants.
190  mint(105)=mint(102+iside)
191  mint(109)=mint(106+iside)
192  CALL pyspli(mint(10+iside),mint(12+iside),kflch,kflsp)
193  IF(mint(51).NE.0) THEN
194  mint(51)=0
195  goto 180
196  ENDIF
197 
198 C...Store first remnant parton, with colour info and kinematics.
199  i=n+1
200  k(i,1)=1
201  k(i,2)=kflsp
202  k(i,3)=mint(83)+iside
203  p(i,5)=pymass(k(i,2))
204  kcol=kchg(pycomp(kflsp),2)
205  IF(kcol.NE.0) THEN
206  k(i,1)=3
207  kfls=(3-kcol*isign(1,kflsp))/2
208  k(i,kfls+3)=mstu(5)*ipu3
209  k(ipu3,6-kfls)=mstu(5)*i
210  icolr=i
211  ENDIF
212  IF(kflch.EQ.0) THEN
213  p(i,1)=-p(ipu3,1)
214  p(i,2)=-p(ipu3,2)
215  pms(iside)=p(i,5)**2+p(i,1)**2+p(i,2)**2
216  p(i,3)=-p(ipu3,3)
217  p(i,4)=sqrt(pms(iside)+p(i,3)**2)
218  prp=p(i,4)+abs(p(i,3))
219 
220 C...When extra remnant parton or hadron: store extra remnant.
221  ELSE
222  i=i+1
223  k(i,1)=1
224  k(i,2)=kflch
225  k(i,3)=mint(83)+iside
226  p(i,5)=pymass(k(i,2))
227  kcol=kchg(pycomp(kflch),2)
228  IF(kcol.NE.0) THEN
229  k(i,1)=3
230  kfls=(3-kcol*isign(1,kflch))/2
231  k(i,kfls+3)=mstu(5)*ipu3
232  k(ipu3,6-kfls)=mstu(5)*i
233  icolr=i
234  ENDIF
235 
236 C...Relative transverse momentum when two remnants.
237  loop=0
238  200 loop=loop+1
239  CALL pyptdi(1,p(i-1,1),p(i-1,2))
240  p(i-1,1)=p(i-1,1)-0.5d0*p(ipu3,1)
241  p(i-1,2)=p(i-1,2)-0.5d0*p(ipu3,2)
242  pms(3)=p(i-1,5)**2+p(i-1,1)**2+p(i-1,2)**2
243  p(i,1)=-p(ipu3,1)-p(i-1,1)
244  p(i,2)=-p(ipu3,2)-p(i-1,2)
245  pms(4)=p(i,5)**2+p(i,1)**2+p(i,2)**2
246 
247 C...Relative distribution of energy for particle into jet plus particle.
248  imb=1
249  IF(mod(mint(10+iside)/1000,10).NE.0) imb=2
250  IF(mstp(94).LE.1) THEN
251  IF(imb.EQ.1) chi=pyr(0)
252  IF(imb.EQ.2) chi=1d0-sqrt(pyr(0))
253  IF(mod(kflch/1000,10).NE.0) chi=1d0-chi
254  ELSEIF(mstp(94).EQ.2) THEN
255  chi=1d0-pyr(0)**(1d0/(1d0+parp(93+2*imb)))
256  IF(mod(kflch/1000,10).NE.0) chi=1d0-chi
257  ELSEIF(mstp(94).EQ.3) THEN
258  CALL pyzdis(1,0,pms(4),zz)
259  chi=zz
260  ELSE
261  CALL pyzdis(1000,0,pms(4),zz)
262  chi=zz
263  ENDIF
264 
265 C...Construct total transverse mass; reject if too large.
266  chi=max(1d-8,min(1d0-1d-8,chi))
267  pms(iside)=pms(4)/chi+pms(3)/(1d0-chi)
268  IF(pms(iside).GT.p(ipu3,4)**2) THEN
269  IF(loop.LT.10) goto 200
270  goto 180
271  ENDIF
272  vint(158+iside)=chi
273 
274 C...Subdivide longitudinal momentum according to value selected above.
275  prp=sqrt(pms(iside)+p(ipu3,3)**2)+abs(p(ipu3,3))
276  pw1=(1d0-chi)*prp
277  p(i-1,4)=0.5d0*(pw1+pms(3)/pw1)
278  p(i-1,3)=0.5d0*(pw1-pms(3)/pw1)*sidesg
279  pw2=chi*prp
280  p(i,4)=0.5d0*(pw2+pms(4)/pw2)
281  p(i,3)=0.5d0*(pw2-pms(4)/pw2)*sidesg
282  ENDIF
283  n=i
284 
285 C...Boost current and remnant systems to correct frame.
286  IF(sqrt(pms(1))+sqrt(pms(2)).GT.0.99d0*vint(1)) goto 180
287  dsqlam=sqrt(max(0d0,(vint(2)-pms(1)-pms(2))**2-4d0*pms(1)*pms(2)))
288  drkc=(vint(2)+pms(3-iside)-pms(iside)+dsqlam)/
289  &(2d0*vint(1)*pcp)
290  drkr=(vint(2)+pms(iside)-pms(3-iside)+dsqlam)/
291  &(2d0*vint(1)*prp)
292  dbec=-sidesg*(drkc**2-1d0)/(drkc**2+1d0)
293  dber=sidesg*(drkr**2-1d0)/(drkr**2+1d0)
294  CALL pyrobo(ipu3,ipu3,0d0,0d0,0d0,0d0,dbec)
295  CALL pyrobo(ipu3+1,n,0d0,0d0,0d0,0d0,dber)
296 
297 C...Let current quark shower; recoil but no showering by colour partner.
298  qmax=2d0*sqrt(vint(309-iside))
299  mstj48=mstj(48)
300  mstj(48)=1
301  parj86=parj(86)
302  parj(86)=0d0
303  IF(mstp(71).EQ.1) CALL pyshow(ipu3,icolr,qmax)
304  mstj(48)=mstj48
305  parj(86)=parj86
306 
307  RETURN
308  END