Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pyhiremn.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pyhiremn.f
1 
2 C*********************************************************************
3 
4  SUBROUTINE pyhiremn(IPU1,IPU2)
5 
6 C...Adds on target remnants (one or two from each side) and
7 C...includes primordial kT.
8  common/hiparnt/hipr1(100),ihpr2(50),hint1(100),ihnt2(50)
9  SAVE /hiparnt/
10  common/histrng/nfp(300,15),pphi(300,15),nft(300,15),pthi(300,15)
11  SAVE /histrng/
12 C...COMMON BLOCK FROM HIJING
13  common/lujets/n,k(9000,5),p(9000,5),v(9000,5)
14  SAVE /lujets/
15  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
16  SAVE /ludat1/
17  common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4)
18  SAVE /ludat2/
19  common/pyhipars/mstp(200),parp(200),msti(200),pari(200)
20  SAVE /pyhipars/
21  common/pyhiint1/mint(400),vint(400)
22  SAVE /pyhiint1/
23  dimension kflch(2),kflsp(2),chi(2),pms(6),is(2),robo(5)
24 
25 C...Special case for lepton-lepton interaction.
26  IF(mint(43).EQ.1) THEN
27  DO 100 jt=1,2
28  i=mint(83)+jt+2
29  k(i,1)=21
30  k(i,2)=k(i-2,2)
31  k(i,3)=i-2
32  DO 100 j=1,5
33  100 p(i,j)=p(i-2,j)
34  ENDIF
35 
36 C...Find event type, set pointers.
37  IF(ipu1.EQ.0.AND.ipu2.EQ.0) RETURN
38  isub=mint(1)
39  ilep=0
40  IF(ipu1.EQ.0) ilep=1
41  IF(ipu2.EQ.0) ilep=2
42  IF(isub.EQ.95) ilep=-1
43  IF(ilep.EQ.1) iq=mint(84)+1
44  IF(ilep.EQ.2) iq=mint(84)+2
45  ip=max(ipu1,ipu2)
46  ilepr=mint(83)+5-ilep
47  ns=n
48 
49 C...Define initial partons, including primordial kT.
50  110 DO 130 jt=1,2
51  i=mint(83)+jt+2
52  IF(jt.EQ.1) ipu=ipu1
53  IF(jt.EQ.2) ipu=ipu2
54  k(i,1)=21
55  k(i,3)=i-2
56  IF(isub.EQ.95) THEN
57  k(i,2)=21
58  shs=0.
59  ELSEIF(mint(40+jt).EQ.1.AND.ipu.NE.0) THEN
60  k(i,2)=k(ipu,2)
61  p(i,5)=p(ipu,5)
62  p(i,1)=0.
63  p(i,2)=0.
64  pms(jt)=p(i,5)**2
65  ELSEIF(ipu.NE.0) THEN
66  k(i,2)=k(ipu,2)
67  p(i,5)=p(ipu,5)
68 C...No primordial kT or chosen according to truncated Gaussian or
69 C...exponential.
70 C
71 c X.N. Wang (7.22.97)
72 c
73  rpt1=0.0
74  rpt2=0.0
75  ss_w2=(pphi(ihnt2(11),4)+pthi(ihnt2(12),4))**2
76  & -(pphi(ihnt2(11),1)+pthi(ihnt2(12),1))**2
77  & -(pphi(ihnt2(11),2)+pthi(ihnt2(12),2))**2
78  & -(pphi(ihnt2(11),3)+pthi(ihnt2(12),3))**2
79 C
80 C********this is s of the current NN collision
81  IF(ss_w2.LE.4.0*parp(93)**2) goto 1211
82 c
83  IF(ihpr2(5).LE.0) THEN
84 120 IF(mstp(91).LE.0) THEN
85  pt=0.
86  ELSEIF(mstp(91).EQ.1) THEN
87  pt=parp(91)*sqrt(-log(rlu(0)))
88  ELSE
89  rpt1=rlu(0)
90  rpt2=rlu(0)
91  pt=-parp(92)*log(rpt1*rpt2)
92  ENDIF
93  IF(pt.GT.parp(93)) goto 120
94  phi=paru(2)*rlu(0)
95  rpt1=pt*cos(phi)
96  rpt2=pt*sin(phi)
97  ELSE IF(ihpr2(5).EQ.1) THEN
98  IF(jt.EQ.1) jpt=nfp(ihnt2(11),11)
99  IF(jt.EQ.2) jpt=nft(ihnt2(12),11)
100 1205 ptgs=parp(91)*sqrt(-log(rlu(0)))
101  IF(ptgs.GT.parp(93)) go to 1205
102  phi=2.0*hipr1(40)*rlu(0)
103  rpt1=ptgs*cos(phi)
104  rpt2=ptgs*sin(phi)
105  DO 1210 i_int=1,jpt-1
106  pkcsq=parp(91)*sqrt(-log(rlu(0)))
107  phi=2.0*hipr1(40)*rlu(0)
108  rpt1=rpt1+pkcsq*cos(phi)
109  rpt2=rpt2+pkcsq*sin(phi)
110 1210 CONTINUE
111  IF(rpt1**2+rpt2**2.GE.ss_w2/4.0) go to 1205
112  ENDIF
113 C X.N. Wang
114 C ********When initial interaction among soft partons is
115 C assumed the primordial pt comes from the sum of
116 C pt of JPT-1 number of initial interaction, JPT
117 C is the number of interaction including present
118 C one that nucleon hassuffered
119 1211 p(i,1)=rpt1
120  p(i,2)=rpt2
121  pms(jt)=p(i,5)**2+p(i,1)**2+p(i,2)**2
122  ELSE
123  k(i,2)=k(iq,2)
124  q2=vint(52)
125  p(i,5)=-sqrt(q2)
126  pms(jt)=-q2
127  shs=(1.-vint(43-jt))*q2/vint(43-jt)+vint(5-jt)**2
128  ENDIF
129  130 CONTINUE
130 
131 C...Kinematics construction for initial partons.
132  i1=mint(83)+3
133  i2=mint(83)+4
134  IF(ilep.EQ.0) shs=vint(141)*vint(142)*vint(2)+
135  &(p(i1,1)+p(i2,1))**2+(p(i1,2)+p(i2,2))**2
136  shr=sqrt(max(0.,shs))
137  IF(ilep.EQ.0) THEN
138  IF((shs-pms(1)-pms(2))**2-4.*pms(1)*pms(2).LE.0.) goto 110
139  p(i1,4)=0.5*(shr+(pms(1)-pms(2))/shr)
140  p(i1,3)=sqrt(max(0.,p(i1,4)**2-pms(1)))
141  p(i2,4)=shr-p(i1,4)
142  p(i2,3)=-p(i1,3)
143  ELSEIF(ilep.EQ.1) THEN
144  p(i1,4)=p(iq,4)
145  p(i1,3)=p(iq,3)
146  p(i2,4)=p(ip,4)
147  p(i2,3)=p(ip,3)
148  ELSEIF(ilep.EQ.2) THEN
149  p(i1,4)=p(ip,4)
150  p(i1,3)=p(ip,3)
151  p(i2,4)=p(iq,4)
152  p(i2,3)=p(iq,3)
153  ENDIF
154  IF(mint(43).EQ.1) RETURN
155 
156 C...Transform partons to overall CM-frame (not for leptoproduction).
157  IF(ilep.EQ.0) THEN
158  robo(3)=(p(i1,1)+p(i2,1))/shr
159  robo(4)=(p(i1,2)+p(i2,2))/shr
160  CALL ludbrb(i1,i2,0.,0.,-dble(robo(3)),-dble(robo(4)),0d0)
161  robo(2)=ulangl(p(i1,1),p(i1,2))
162  CALL ludbrb(i1,i2,0.,-robo(2),0d0,0d0,0d0)
163  robo(1)=ulangl(p(i1,3),p(i1,1))
164  CALL ludbrb(i1,i2,-robo(1),0.,0d0,0d0,0d0)
165  nmax=max(mint(52),ipu1,ipu2)
166  CALL ludbrb(i1,nmax,robo(1),robo(2),dble(robo(3)),dble(robo(4)),
167  & 0d0)
168  robo(5)=max(-0.999999,min(0.999999,(vint(141)-vint(142))/
169  & (vint(141)+vint(142))))
170  CALL ludbrb(i1,nmax,0.,0.,0d0,0d0,dble(robo(5)))
171  ENDIF
172 
173 C...Check invariant mass of remnant system:
174 C...hadronic events or leptoproduction.
175  IF(ilep.LE.0) THEN
176  IF(mstp(81).LE.0.OR.mstp(82).LE.0.OR.isub.EQ.95) THEN
177  vint(151)=0.
178  vint(152)=0.
179  ENDIF
180  peh=p(i1,4)+p(i2,4)+0.5*vint(1)*(vint(151)+vint(152))
181  pzh=p(i1,3)+p(i2,3)+0.5*vint(1)*(vint(151)-vint(152))
182  shh=(vint(1)-peh)**2-(p(i1,1)+p(i2,1))**2-(p(i1,2)+p(i2,2))**2-
183  & pzh**2
184  pmmin=p(mint(83)+1,5)+p(mint(83)+2,5)+ulmass(k(i1,2))+
185  & ulmass(k(i2,2))
186  IF(shr.GE.vint(1).OR.shh.LE.(pmmin+parp(111))**2) THEN
187  mint(51)=1
188  RETURN
189  ENDIF
190  shr=sqrt(shh+(p(i1,1)+p(i2,1))**2+(p(i1,2)+p(i2,2))**2)
191  ELSE
192  pei=p(iq,4)+p(ip,4)
193  pzi=p(iq,3)+p(ip,3)
194  pms(ilep)=max(0.,pei**2-pzi**2)
195  pmmin=p(ilepr-2,5)+ulmass(k(ilepr,2))+sqrt(pms(ilep))
196  IF(shr.LE.pmmin+parp(111)) THEN
197  mint(51)=1
198  RETURN
199  ENDIF
200  ENDIF
201 
202 C...Subdivide remnant if necessary, store first parton.
203  140 i=ns
204  DO 190 jt=1,2
205  IF(jt.EQ.ilep) goto 190
206  IF(jt.EQ.1) ipu=ipu1
207  IF(jt.EQ.2) ipu=ipu2
208  CALL pyhispli(mint(10+jt),mint(12+jt),kflch(jt),kflsp(jt))
209  i=i+1
210  is(jt)=i
211  DO 150 j=1,5
212  k(i,j)=0
213  p(i,j)=0.
214  150 v(i,j)=0.
215  k(i,1)=3
216  k(i,2)=kflsp(jt)
217  k(i,3)=mint(83)+jt
218  p(i,5)=ulmass(k(i,2))
219 
220 C...First parton colour connections and transverse mass.
221  kfls=(3-kchg(lucomp(kflsp(jt)),2)*isign(1,kflsp(jt)))/2
222  k(i,kfls+3)=ipu
223  k(ipu,6-kfls)=mod(k(ipu,6-kfls),mstu(5))+mstu(5)*i
224  IF(kflch(jt).EQ.0) THEN
225  p(i,1)=-p(mint(83)+jt+2,1)
226  p(i,2)=-p(mint(83)+jt+2,2)
227  pms(jt)=p(i,5)**2+p(i,1)**2+p(i,2)**2
228 
229 C...When extra remnant parton or hadron: find relative pT, store.
230  ELSE
231  CALL luptdi(1,p(i,1),p(i,2))
232  pms(jt+2)=p(i,5)**2+p(i,1)**2+p(i,2)**2
233  i=i+1
234  DO 160 j=1,5
235  k(i,j)=0
236  p(i,j)=0.
237  160 v(i,j)=0.
238  k(i,1)=1
239  k(i,2)=kflch(jt)
240  k(i,3)=mint(83)+jt
241  p(i,5)=ulmass(k(i,2))
242  p(i,1)=-p(mint(83)+jt+2,1)-p(i-1,1)
243  p(i,2)=-p(mint(83)+jt+2,2)-p(i-1,2)
244  pms(jt+4)=p(i,5)**2+p(i,1)**2+p(i,2)**2
245 C...Relative distribution of energy for particle into two jets.
246  imb=1
247  IF(mod(mint(10+jt)/1000,10).NE.0) imb=2
248  IF(iabs(kflch(jt)).LE.10.OR.kflch(jt).EQ.21) THEN
249  chik=parp(92+2*imb)
250  IF(mstp(92).LE.1) THEN
251  IF(imb.EQ.1) chi(jt)=rlu(0)
252  IF(imb.EQ.2) chi(jt)=1.-sqrt(rlu(0))
253  ELSEIF(mstp(92).EQ.2) THEN
254  chi(jt)=1.-rlu(0)**(1./(1.+chik))
255  ELSEIF(mstp(92).EQ.3) THEN
256  cut=2.*0.3/vint(1)
257  170 chi(jt)=rlu(0)**2
258  IF((chi(jt)**2/(chi(jt)**2+cut**2))**0.25*(1.-chi(jt))**chik
259  & .LT.rlu(0)) goto 170
260  ELSE
261  cut=2.*0.3/vint(1)
262  cutr=(1.+sqrt(1.+cut**2))/cut
263  180 chir=cut*cutr**rlu(0)
264  chi(jt)=(chir**2-cut**2)/(2.*chir)
265  IF((1.-chi(jt))**chik.LT.rlu(0)) goto 180
266  ENDIF
267 C...Relative distribution of energy for particle into jet plus particle.
268  ELSE
269  IF(mstp(92).LE.1) THEN
270  IF(imb.EQ.1) chi(jt)=rlu(0)
271  IF(imb.EQ.2) chi(jt)=1.-sqrt(rlu(0))
272  ELSE
273  chi(jt)=1.-rlu(0)**(1./(1.+parp(93+2*imb)))
274  ENDIF
275  IF(mod(kflch(jt)/1000,10).NE.0) chi(jt)=1.-chi(jt)
276  ENDIF
277  pms(jt)=pms(jt+4)/chi(jt)+pms(jt+2)/(1.-chi(jt))
278  kfls=kchg(lucomp(kflch(jt)),2)*isign(1,kflch(jt))
279  IF(kfls.NE.0) THEN
280  k(i,1)=3
281  kfls=(3-kfls)/2
282  k(i,kfls+3)=ipu
283  k(ipu,6-kfls)=mod(k(ipu,6-kfls),mstu(5))+mstu(5)*i
284  ENDIF
285  ENDIF
286  190 CONTINUE
287  IF(shr.LE.sqrt(pms(1))+sqrt(pms(2))) goto 140
288  n=i
289 
290 C...Reconstruct kinematics of remnants.
291  DO 200 jt=1,2
292  IF(jt.EQ.ilep) goto 200
293  pe=0.5*(shr+(pms(jt)-pms(3-jt))/shr)
294  pz=sqrt(pe**2-pms(jt))
295  IF(kflch(jt).EQ.0) THEN
296  p(is(jt),4)=pe
297  p(is(jt),3)=pz*(-1)**(jt-1)
298  ELSE
299  pw1=chi(jt)*(pe+pz)
300  p(is(jt)+1,4)=0.5*(pw1+pms(jt+4)/pw1)
301  p(is(jt)+1,3)=0.5*(pw1-pms(jt+4)/pw1)*(-1)**(jt-1)
302  p(is(jt),4)=pe-p(is(jt)+1,4)
303  p(is(jt),3)=pz*(-1)**(jt-1)-p(is(jt)+1,3)
304  ENDIF
305  200 CONTINUE
306 
307 C...Hadronic events: boost remnants to correct longitudinal frame.
308  IF(ilep.LE.0) THEN
309  CALL ludbrb(ns+1,n,0.,0.,0d0,0d0,-dble(pzh/(vint(1)-peh)))
310 C...Leptoproduction events: boost colliding subsystem.
311  ELSE
312  nmax=max(ip,mint(52))
313  pef=shr-pe
314  pzf=pz*(-1)**(ilep-1)
315  pt2=p(ilepr,1)**2+p(ilepr,2)**2
316  phipt=ulangl(p(ilepr,1),p(ilepr,2))
317  CALL ludbrb(mint(84)+1,nmax,0.,-phipt,0d0,0d0,0d0)
318  rqp=p(iq,3)*(pt2+pei**2)-p(iq,4)*pei*pzi
319  sinth=p(iq,4)*sqrt(pt2*(pt2+pei**2)/(rqp**2+pt2*
320  & p(iq,4)**2*pzi**2))*sign(1.,-rqp)
321  CALL ludbrb(mint(84)+1,nmax,asin(sinth),0.,0d0,0d0,0d0)
322  betax=(-pei*pzi*sinth+sqrt(pt2*(pt2+pei**2-(pzi*sinth)**2)))/
323  & (pt2+pei**2)
324  CALL ludbrb(mint(84)+1,nmax,0.,0.,dble(betax),0d0,0d0)
325  CALL ludbrb(mint(84)+1,nmax,0.,phipt,0d0,0d0,0d0)
326  pem=p(iq,4)+p(ip,4)
327  pzm=p(iq,3)+p(ip,3)
328  betaz=(-pem*pzm+pzf*sqrt(pzf**2+pem**2-pzm**2))/(pzf**2+pem**2)
329  CALL ludbrb(mint(84)+1,nmax,0.,0.,0d0,0d0,dble(betaz))
330  CALL ludbrb(i1,i2,asin(sinth),0.,dble(betax),0d0,0d0)
331  CALL ludbrb(i1,i2,0.,phipt,0d0,0d0,dble(betaz))
332  ENDIF
333 
334  RETURN
335  END