EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
lyremn.F
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file lyremn.F
1 
2 C **********************************************************************
3 
4  SUBROUTINE lyremn(IPU1,IPU2)
5 
6  IMPLICIT NONE
7 
8 C...ADDS ON TARGET REMNANTS (ONE OR TWO FROM EACH SIDE) AND
9 C...INCLUDES PRIMORDIAL KT.
10 *
11 * to avoid variable conflictions, a second keep element is necessary
12 * with the same common block name (see LEPTO2)
13 *
14 
15  COMMON /leptou/ cut(14),lst(40),parl(30),
16  & xlp,ylp,w2lp,q2lp,ulp
17  REAL cut,parl,xlp,ylp,w2lp,q2lp,ulp
18  INTEGER lst
19  SAVE /leptou/
20 
21  INTEGER nlupdm,nplbuf
22  parameter(nlupdm=4000,nplbuf=5)
23  common/lujets/n,k(nlupdm,5),p(nlupdm,nplbuf),v(nlupdm,5)
24  INTEGER n,k
25  REAL p,v
26  SAVE /lujets/
27 
28  common/ludat1/mstu(200),paru(200),mstj(200),parj(200)
29  INTEGER mstu,mstj
30  REAL paru,parj
31  SAVE /ludat1/
32 
33  COMMON /pypara/ ipy(80),pypar(80),pyvar(80)
34  REAL pypar,pyvar
35  INTEGER ipy
36  SAVE /pypara/
37 
38 *
39 * to avoid variable conflictions, a second keep element is necessary
40 * with the same common block name (see LYPRO2)
41 *
42  COMMON /lyproc/ isub,kfl(3,2),x(2),sh,th,uh,q2,xsec(0:40)
43  REAL x,sh,th,uh,q2,xsec
44  INTEGER isub,kfl
45  SAVE /lyproc/
46 
47 
48  INTEGER ipu1,ipu2,kflch,kflsp,ipu,iq,ilep,ip,ns,i,j,
49  +jt,ifls,imin,imax,is,kfi2,ntry
50  INTEGER lucomp
51  REAL chi,pms,robo,pei,pe,pzi,pz,shs,pzh,peh,ptspl,
52  +phispl,shr,pmmin,pw1,pef,pzf,pt2,shh,phipt,rqp,sinth,
53  +pem,pzm
54  REAL ulmass,ulangl,amk42,amk32
55  dimension kflch(2),kflsp(2),chi(2),pms(6),is(2),robo(5)
56  DOUBLE PRECISION dbetax,dbetaz,drobo(5)
57  DATA ipu,iq/0,0/,pei,pe,pzi,pz,shs,pzh,peh/7*0./
58 
59 C...FIND EVENT TYPE, SET POINTERS
60  IF(ipu1.EQ.0.AND.ipu2.EQ.0) RETURN
61  ilep=0
62  IF(ipu1.EQ.0) ilep=1
63  IF(ipu2.EQ.0) ilep=2
64  IF(isub.EQ.7) ilep=-1
65  IF(ilep.EQ.1) iq=21
66  IF(ilep.EQ.2) iq=23
67  ip=max(ipu1,ipu2)
68  ns=n
69 
70 C...DEFINE INITIAL PARTONS, INCLUDING PRIMORDIAL KT
71  100 DO 120 i=3,4
72  IF(i.EQ.3) ipu=ipu1
73  IF(i.EQ.4) ipu=ipu2
74  k(i,1)=21
75  k(i,3)=i-2
76  DO 110 j=1,5
77  110 p(i,j)=0.
78  IF(isub.EQ.7) THEN
79  k(i,2)=21
80  shs=0.
81  ELSEIF(ipu.NE.0) THEN
82  k(i,2)=k(ipu,2)
83  p(i,5)=p(ipu,5)
84  CALL lprikt(parl(3),ptspl,phispl)
85  p(i,1)=ptspl*cos(phispl)
86  p(i,2)=ptspl*sin(phispl)
87  pms(i-2)=p(i,5)**2+p(i,1)**2+p(i,2)**2
88  ELSE
89  k(i,2)=k(iq,2)
90  p(i,5)=-sqrt(q2)
91  pms(i-2)=-q2
92  shs=(1.-x(5-i))*q2/x(5-i)+pyvar(7-i)**2
93  ENDIF
94  120 CONTINUE
95 
96 C...KINEMATICS CONSTRUCTION FOR INITIAL PARTONS
97  IF(ilep.EQ.0) shs=pyvar(31)*pyvar(32)*pyvar(2)+
98  &(p(3,1)+p(4,1))**2+(p(3,2)+p(4,2))**2
99  shr=sqrt(max(0.,shs))
100  IF(ilep.EQ.0) THEN
101  IF((shs-pms(1)-pms(2))**2-4.*pms(1)*pms(2).LE.0.) goto 100
102  p(3,4)=0.5*(shr+(pms(1)-pms(2))/shr)
103  p(3,3)=sqrt(max(0.,p(3,4)**2-pms(1)))
104  p(4,4)=shr-p(3,4)
105  p(4,3)=-p(3,3)
106  ELSEIF(ilep.EQ.1) THEN
107  p(3,4)=p(iq,4)
108  p(3,3)=p(iq,3)
109  p(4,4)=p(ip,4)
110  p(4,3)=p(ip,3)
111  ELSEIF(ilep.EQ.2) THEN
112  p(3,4)=p(ip,4)
113  p(3,3)=p(ip,3)
114  p(4,4)=p(iq,4)
115  p(4,3)=p(iq,3)
116  ENDIF
117 
118 C...TRANSFORM PARTONS TO OVERALL CM-FRAME (NOT FOR LEPTOPRODUCTION)
119  IF(ilep.EQ.0) THEN
120  mstu(1)=3
121  mstu(2)=4
122  drobo(3)=(p(3,1)+p(4,1))/shr
123  drobo(4)=(p(3,2)+p(4,2))/shr
124  CALL ludbrb(mstu(1),mstu(2),0.,0.,-drobo(3),-drobo(4),0.d0)
125  robo(2)=ulangl(p(3,1),p(3,2))
126  CALL ludbrb(mstu(1),mstu(2),0.,-robo(2),0.d0,0.d0,0.d0)
127  robo(1)=ulangl(p(3,3),p(3,1))
128  CALL ludbrb(mstu(1),mstu(2),-robo(1),0.,0.d0,0.d0,0.d0)
129  mstu(2)=max(ipy(47),ipu1,ipu2)
130  CALL ludbrb(mstu(1),mstu(2),
131  & robo(1),robo(2),drobo(3),drobo(4),0.d0)
132  drobo(5)=max(-0.999999,min(0.999999,(pyvar(31)-pyvar(32))/
133  & (pyvar(31)+pyvar(32))))
134  CALL ludbrb(mstu(1),mstu(2),0.,0.,0.d0,0.d0,drobo(5))
135  mstu(1)=0
136  mstu(2)=0
137  ENDIF
138 
139 C...CHECK INVARIANT MASS OF REMNANT SYSTEM:
140 C...HADRONIC EVENTS OR LEPTOPRODUCTION
141  IF(ilep.LE.0) THEN
142  IF(ipy(12).LE.0.OR.isub.EQ.7) pyvar(33)=0.
143  IF(ipy(12).LE.0.OR.isub.EQ.7) pyvar(34)=0.
144  peh=p(3,4)+p(4,4)+0.5*pyvar(1)*(pyvar(33)+pyvar(34))
145  pzh=p(3,3)+p(4,3)+0.5*pyvar(1)*(pyvar(33)-pyvar(34))
146  shh=(pyvar(1)-peh)**2-(p(3,1)+p(4,1))**2-(p(3,2)+p(4,2))**2-
147  & pzh**2
148  mstj(93)=1
149  amk32=ulmass(k(3,2))
150  mstj(93)=1
151  amk42=ulmass(k(4,2))
152  pmmin=p(1,5)+p(2,5)+amk32+amk42
153  IF(shr.GE.pyvar(1).OR.shh.LE.(pmmin+pypar(12))**2) THEN
154  ipy(48)=1
155  RETURN
156  ENDIF
157  shr=sqrt(shh+(p(3,1)+p(4,1))**2+(p(3,2)+p(4,2))**2)
158  ELSE
159  pei=p(iq,4)+p(ip,4)
160  pzi=p(iq,3)+p(ip,3)
161  pms(ilep)=max(0.,pei**2-pzi**2+p(5-ilep,1)**2+p(5-ilep,2)**2)
162  mstj(93)=1
163  pmmin=p(3-ilep,5)+ulmass(k(5-ilep,2))+sqrt(pms(ilep))
164  IF(shr.LE.pmmin+pypar(12)) THEN
165  ipy(48)=1
166  RETURN
167  ENDIF
168  ENDIF
169 
170 C... SUBDIVIDE REMNANT IF NECESSARY, STORE FIRST PARTON
171 CJR-- begin
172 CJR-- try to find kinematically allowed solution
173 CJR-- no more than 100 times
174  ntry=0
175  130 CONTINUE
176  ntry=ntry+1
177  IF(ntry.GT.100) THEN
178  ipy(48)=1
179  RETURN
180  ENDIF
181 CJR-- end
182  i=ns-1
183  DO 160 jt=1,2
184  IF(jt.EQ.ilep) goto 160
185  IF(jt.EQ.1) ipu=ipu1
186  IF(jt.EQ.2) ipu=ipu2
187  CALL lyspli(ipy(40+jt),kfl(1,jt),kflch(jt),kflsp(jt))
188  i=i+2
189  is(jt)=i
190  k(i,1)=3
191  k(i,2)=kflsp(jt)
192  k(i,3)=jt
193  mstj(93)=1
194  p(i,5)=ulmass(k(i,2))
195 CJR--
196  kfi2=lucomp(k(i,2))
197  IF (kfi2.EQ.90) THEN
198  p(i,5)=p(i,5)-2.*parl(20)
199  ELSEIF (1.LE.kfi2 .AND. kfi2.LE.6) THEN
200  p(i,5)=p(i,5)-parl(20)
201  ENDIF
202 CJR--
203 C...FIRST PARTON COLOUR CONNECTIONS AND TRANSVERSE MASS
204  k(i+1,1)=-1
205  k(i+1,3)=i
206  k(i+1,2)=1000
207  IF(ipy(34).GE.1) k(i+1,2)=1000+jt
208  DO 140 j=1,5
209  140 p(i+1,j)=0.
210  IF(kflsp(jt).EQ.21) THEN
211  p(i+1,3)=ipu
212  p(i+1,4)=ipu
213  p(ipu+1,1)=i
214  p(ipu+1,2)=i
215  k(i,4)=ipu+ipu*mstu(5)
216  k(i,5)=ipu+ipu*mstu(5)
217  k(ipu,4)=mod(k(ipu,4),mstu(5))+i*mstu(5)
218  k(ipu,5)=mod(k(ipu,5),mstu(5))+i*mstu(5)
219  ELSE
220  ifls=(3-isign(1,kflsp(jt)*(1102-iabs(kflsp(jt)))))/2
221  p(i+1,ifls+2)=ipu
222  p(ipu+1,3-ifls)=i
223  k(i,ifls+3)=ipu
224  k(ipu,6-ifls)=mod(k(ipu,6-ifls),mstu(5))+i*mstu(5)
225  ENDIF
226  IF(kflch(jt).EQ.0) THEN
227  p(i,1)=-p(jt+2,1)
228  p(i,2)=-p(jt+2,2)
229  pms(jt)=p(i,5)**2+p(i,1)**2+p(i,2)**2
230  ELSE
231 C...WHEN EXTRA REMNANT PARTON OR HADRON: FIND RELATIVE PT, STORE
232 C...PRIMORDIAL KT SPLIT SHARED BETWEEN REMNANTS
233  CALL lprikt(parl(14),ptspl,phispl)
234 C...RELATIVE DISTRIBUTION OF ENERGY; EXTRA PARTON COLOUR CONNECTION
235  CALL lremh(0,ptspl,kflsp(jt),kflch(jt),chi(jt))
236  p(i,1)=-p(jt+2,1)*(1.-chi(jt))+ptspl*cos(phispl)
237  p(i,2)=-p(jt+2,2)*(1.-chi(jt))+ptspl*sin(phispl)
238  pms(jt+2)=p(i,5)**2+p(i,1)**2+p(i,2)**2
239  i=i+2
240  DO 150 j=1,5
241  k(i,j)=0
242  k(i+1,j)=0
243  p(i,j)=0.
244  150 p(i+1,j)=0.
245  k(i,1)=1
246  k(i,2)=kflch(jt)
247  k(i,3)=jt
248  mstj(93)=1
249  p(i,5)=ulmass(k(i,2))
250 CJR--
251  kfi2=lucomp(k(i,2))
252  IF (kfi2.EQ.90) THEN
253  p(i,5)=p(i,5)-2.*parl(20)
254  ELSEIF (1.LE.kfi2 .AND. kfi2.LE.6) THEN
255  p(i,5)=p(i,5)-parl(20)
256  ENDIF
257 CJR--
258  p(i,1)=-p(jt+2,1)*chi(jt)-ptspl*cos(phispl)
259  p(i,2)=-p(jt+2,2)*chi(jt)-ptspl*sin(phispl)
260  pms(jt+4)=p(i,5)**2+p(i,1)**2+p(i,2)**2
261 C...end of update
262  pms(jt)=pms(jt+4)/chi(jt)+pms(jt+2)/(1.-chi(jt))
263  k(i+1,1)=-1
264  k(i+1,3)=i
265  k(i+1,2)=1000
266  IF(ipy(34).GE.1) k(i+1,2)=1000+jt
267  IF((iabs(kflch(jt)).GE.1.AND.iabs(kflch(jt)).LE.8).OR.
268  & iabs(kflch(jt)).EQ.21.OR.lucomp(iabs(kflch(jt))).EQ.90) THEN
269  ifls=(3-isign(1,kflch(jt)*(1102-iabs(kflch(jt)))))/2
270  p(i+1,ifls+2)=ipu
271  p(ipu+1,3-ifls)=i
272  k(i,1)=3
273  k(i,ifls+3)=ipu
274  k(ipu,6-ifls)=mod(k(ipu,6-ifls),mstu(5))+i*mstu(5)
275  ELSE
276  IF(ipy(34).GE.1) THEN
277  k(i,1)=1
278  k(i,3)=jt
279  ENDIF
280  ENDIF
281  ENDIF
282  160 CONTINUE
283  IF(shr.LE.sqrt(pms(1))+sqrt(pms(2))) goto 130
284  n=i+1
285 
286 C...RECONSTRUCT KINEMATICS OF REMNANTS
287  DO 170 jt=1,2
288  IF(jt.EQ.ilep) goto 170
289  pe=0.5*(shr+(pms(jt)-pms(3-jt))/shr)
290  pz=sqrt(pe**2-pms(jt))
291  IF(kflch(jt).EQ.0) THEN
292  p(is(jt),4)=pe
293  p(is(jt),3)=pz*(-1)**(jt-1)
294  ELSE
295  pw1=chi(jt)*(pe+pz)
296  p(is(jt)+2,4)=0.5*(pw1+pms(jt+4)/pw1)
297  p(is(jt)+2,3)=0.5*(pw1-pms(jt+4)/pw1)*(-1)**(jt-1)
298  p(is(jt),4)=pe-p(is(jt)+2,4)
299  p(is(jt),3)=pz*(-1)**(jt-1)-p(is(jt)+2,3)
300  ENDIF
301  170 CONTINUE
302 
303 C CALL GULIST(31,2)
304 C...HADRONIC EVENTS: BOOST REMNANTS TO CORRECT LONGITUDINAL FRAME
305  IF(ilep.LE.0) THEN
306  mstu(1)=ns+1
307  CALL ludbrb(mstu(1),mstu(2),
308  & 0.,0.,0.d0,0.d0,-dble(pzh)/(dble(pyvar(1))-dble(peh)))
309  mstu(1)=0
310 C...LEPTOPRODUCTION EVENTS: BOOST COLLIDING SUBSYSTEM
311  ELSE
312  imin=21
313  imax=max(ip,ipy(47))
314  pef=shr-pe
315  pzf=pz*(-1)**(ilep-1)
316  pt2=p(5-ilep,1)**2+p(5-ilep,2)**2
317  phipt=ulangl(p(5-ilep,1),p(5-ilep,2))
318  CALL ludbrb(imin,imax,0.,-phipt,0.d0,0.d0,0.d0)
319  rqp=p(iq,3)*(pt2+pei**2)-p(iq,4)*pei*pzi
320  sinth=p(iq,4)*sqrt(pt2*(pt2+pei**2)/(rqp**2+pt2*
321  & p(iq,4)**2*pzi**2))*sign(1.,-rqp)
322  CALL ludbrb(imin,imax,asin(sinth),0.,0.d0,0.d0,0.d0)
323  dbetax=(-dble(pei)*pzi*sinth+
324  & sqrt(dble(pt2)*(pt2+pei**2-(pzi*sinth)**2)))/
325  & (dble(pt2)+pei**2)
326  CALL ludbrb(imin,imax,0.,0.,dbetax,0.d0,0.d0)
327  CALL ludbrb(imin,imax,0.,phipt,0.d0,0.d0,0.d0)
328  pem=p(iq,4)+p(ip,4)
329  pzm=p(iq,3)+p(ip,3)
330  dbetaz=(-dble(pem)*pzm+
331  & pzf*sqrt(dble(pzf)**2+pem**2-pzm**2))/(dble(pzf)**2+pem**2)
332  CALL ludbrb(imin,imax,0.,0.,0.d0,0.d0,dbetaz)
333 C...Avoid double application of kt
334  p(4,1)=0.
335  p(4,2)=0.
336  CALL ludbrb(3,4,asin(sinth),0.,dbetax,0.d0,0.d0)
337  CALL ludbrb(3,4,0.,phipt,0.d0,0.d0,dbetaz)
338  ENDIF
339 
340  RETURN
341  END