Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
pygaga.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file pygaga.f
1 C*********************************************************************
2 
3 C...PYGAGA
4 C...For lepton beams it gives photon-hadron or photon-photon systems
5 C...to be treated with the ordinary machinery and combines this with a
6 C...description of the lepton -> lepton + photon branching.
7 
8  SUBROUTINE pygaga(IGAGA,WTGAGA)
9 
10 C...Double precision and integer declarations.
11  IMPLICIT DOUBLE PRECISION(a-h, o-z)
12  IMPLICIT INTEGER(i-n)
13  INTEGER pyk,pychge,pycomp
14 C...Commonblocks.
15  common/pyjets/n,npad,k(4000,5),p(4000,5),v(4000,5)
16  common/pydat1/mstu(200),paru(200),mstj(200),parj(200)
17  common/pydat2/kchg(500,4),pmas(500,4),parf(2000),vckm(4,4)
18  common/pysubs/msel,mselpd,msub(500),kfin(2,-40:40),ckin(200)
19  common/pypars/mstp(200),parp(200),msti(200),pari(200)
20  common/pyint1/mint(400),vint(400)
21  common/pyint5/ngenpd,ngen(0:500,3),xsec(0:500,3)
22  DOUBLE PRECISION minq2, rccorr, pyth_xsec, temp, etheta
23  DOUBLE PRECISION ppt, ppx, ppy
24  include "mcRadCor.inc"
25  include "mc_set.inc"
26  include "radgen.inc"
27  include "phiout.inc"
28  include "py6strf.inc"
29 
30  SAVE /pyjets/,/pydat1/,/pydat2/,/pysubs/,/pypars/,/pyint1/,
31  &/pyint5/
32 C...Local variables and data statement.
33  dimension pms(2),xmin(2),xmax(2),q2min(2),q2max(2),pmc(3),
34  &x(2),q2(2),y(2),theta(2),phi(2),pt(2),beta(3)
35  SAVE pms,xmin,xmax,q2min,q2max,pmc,x,q2,theta,phi,pt,w2min,
36  & ymin,ymax
37  DATA eps/1d-4/
38 
39 C...Initialize generation of photons inside leptons.
40  IF(igaga.EQ.1) THEN
41 
42 C...Save quantities on incoming lepton system.
43  vint(301)=vint(1)
44  vint(302)=vint(2)
45  pms(1)=vint(303)**2
46  IF(mint(141).EQ.0) pms(1)=sign(vint(3)**2,vint(3))
47  pms(2)=vint(304)**2
48  IF(mint(142).EQ.0) pms(2)=sign(vint(4)**2,vint(4))
49  pmc(3)=vint(302)-pms(1)-pms(2)
50  w2min=max(ckin(77),2d0*ckin(3),2d0*ckin(5))**2
51 
52 C...Calculate range of x and Q2 values allowed in generation.
53  DO 100 i=1,2
54  pmc(i)=vint(302)+pms(i)-pms(3-i)
55  IF(mint(140+i).NE.0) THEN
56  xmin(i)=max(ckin(59+2*i),eps)
57  xmax(i)=min(ckin(60+2*i),1d0-2d0*vint(301)*sqrt(pms(i))/
58  & pmc(i),1d0-eps)
59  ymin=max(ckin(71+2*i),eps)
60  ymax=min(ckin(72+2*i),1d0-eps)
61  IF(ckin(64+2*i).GT.0d0) xmin(i)=max(xmin(i),
62  & (ymin*pmc(3)-ckin(64+2*i))/pmc(i))
63  xmax(i)=min(xmax(i),(ymax*pmc(3)-ckin(63+2*i))/pmc(i))
64  themin=max(ckin(67+2*i),0d0)
65  themax=min(ckin(68+2*i),paru(1))
66  IF(ckin(68+2*i).LT.0d0) themax=paru(1)
67  q2min(i)=max(ckin(63+2*i),xmin(i)**2*pms(i)/(1d0-xmin(i))+
68  & ((1d0-xmax(i))*(vint(302)-2d0*pms(3-i))-
69  & 2d0*pms(i)/(1d0-xmax(i)))*sin(themin/2d0)**2,0d0)
70  q2max(i)=xmax(i)**2*pms(i)/(1d0-xmax(i))+
71  & ((1d0-xmin(i))*(vint(302)-2d0*pms(3-i))-
72  & 2d0*pms(i)/(1d0-xmin(i)))*sin(themax/2d0)**2
73  IF(ckin(64+2*i).GT.0d0) q2max(i)=min(ckin(64+2*i),q2max(i))
74 C...W limits when lepton on one side only.
75  IF(mint(143-i).EQ.0) THEN
76  xmin(i)=max(xmin(i),(w2min-pms(3-i))/pmc(i))
77  IF(ckin(78).GT.0d0) xmax(i)=min(xmax(i),
78  & (ckin(78)**2-pms(3-i))/pmc(i))
79  ENDIF
80  ENDIF
81  100 CONTINUE
82 
83 C...W limits when lepton on both sides.
84  IF(mint(141).NE.0.AND.mint(142).NE.0) THEN
85  IF(ckin(78).GT.0d0) xmax(1)=min(xmax(1),
86  & (ckin(78)**2+pmc(3)-pmc(2)*xmin(2))/pmc(1))
87  IF(ckin(78).GT.0d0) xmax(2)=min(xmax(2),
88  & (ckin(78)**2+pmc(3)-pmc(1)*xmin(1))/pmc(2))
89  IF(iabs(mint(141)).NE.iabs(mint(142))) THEN
90  xmin(1)=max(xmin(1),(pms(1)-pms(2)+vint(302)*(w2min-
91  & pms(1)-pms(2))/(pmc(2)*xmax(2)+pms(1)-pms(2)))/pmc(1))
92  xmin(2)=max(xmin(2),(pms(2)-pms(1)+vint(302)*(w2min-
93  & pms(1)-pms(2))/(pmc(1)*xmax(1)+pms(2)-pms(1)))/pmc(2))
94  ELSE
95  xmin(1)=max(xmin(1),w2min/(vint(302)*xmax(2)))
96  xmin(2)=max(xmin(2),w2min/(vint(302)*xmax(1)))
97  ENDIF
98  ENDIF
99 
100 C...Q2 and W values and photon flux weight factors for initialization.
101  ELSEIF(igaga.EQ.2) THEN
102  isub=mint(1)
103  mint(15)=0
104  mint(16)=0
105 
106 C...W value for photon on one or both sides, and for processes
107 C...with gamma-gamma cross section peaked at small shat.
108  IF(mint(141).NE.0.AND.mint(142).EQ.0) THEN
109  vint(2)=vint(302)+pms(1)-pmc(1)*(1d0-xmax(1))
110  ELSEIF(mint(141).EQ.0.AND.mint(142).NE.0) THEN
111  vint(2)=vint(302)+pms(2)-pmc(2)*(1d0-xmax(2))
112  ELSEIF(isub.GE.137.AND.isub.LE.140) THEN
113  vint(2)=max(ckin(77)**2,12d0*max(ckin(3),ckin(5))**2)
114  IF(ckin(78).GT.0d0) vint(2)=min(vint(2),ckin(78)**2)
115  ELSE
116  vint(2)=xmax(1)*xmax(2)*vint(302)
117  IF(ckin(78).GT.0d0) vint(2)=min(vint(2),ckin(78)**2)
118  ENDIF
119  vint(1)=sqrt(max(0d0,vint(2)))
120 
121 C...Upper estimate of photon flux weight factor.
122 C...Initialization Q2 scale. Flag incoming unresolved photon.
123  wtgaga=1d0
124  DO 110 i=1,2
125  IF(mint(140+i).NE.0) THEN
126  IF(mstp(199).EQ.1) then
127  wtgaga=wtgaga*2d0*(paru(101)/paru(2))*
128  & (log(real(mcset_ymax/mcset_ymin)))*
129  & (log(real(mcset_q2max/mcset_q2min)))
130  ELSE
131  wtgaga=wtgaga*2d0*(paru(101)/paru(2))*
132  & log(xmax(i)/xmin(i))*log(q2max(i)/q2min(i))
133  ENDIF
134  IF(isub.EQ.99.AND.mint(106+i).EQ.4.AND.mint(109-i).EQ.3)
135  & THEN
136  q2init=5d0+q2min(3-i)
137  ELSEIF(isub.EQ.99.AND.mint(106+i).EQ.4) THEN
138  q2init=pmas(pycomp(113),1)**2+q2min(3-i)
139  ELSEIF(isub.EQ.132.OR.isub.EQ.134.OR.isub.EQ.136) THEN
140  q2init=max(ckin(1),2d0*ckin(3),2d0*ckin(5))**2/3d0
141  ELSEIF((isub.EQ.138.AND.i.EQ.2).OR.
142  & (isub.EQ.139.AND.i.EQ.1)) THEN
143  q2init=vint(2)/3d0
144  ELSEIF(isub.EQ.140) THEN
145  q2init=vint(2)/2d0
146  ELSE
147  q2init=q2min(i)
148  ENDIF
149  vint(2+i)=-sqrt(max(q2min(i),min(q2max(i),q2init)))
150  IF(mstp(14).EQ.0.OR.(isub.GE.131.AND.isub.LE.140))
151  & mint(14+i)=22
152  vint(306+i)=vint(2+i)**2
153  ENDIF
154  110 CONTINUE
155  vint(320)=wtgaga
156 
157 C...Update pTmin and cross section information.
158  IF(mstp(82).LE.1) THEN
159  ptmn=parp(81)*(vint(1)/parp(89))**parp(90)
160  ELSE
161  ptmn=parp(82)*(vint(1)/parp(89))**parp(90)
162  ENDIF
163  vint(149)=4d0*ptmn**2/vint(2)
164  vint(154)=ptmn
165  CALL pyxtot
166  vint(318)=vint(317)
167 
168 C...Generate photons inside leptons and
169 C...calculate photon flux weight factors.
170  ELSEIF(igaga.EQ.3) THEN
171  isub=mint(1)
172  mint(15)=0
173  mint(16)=0
174 
175 C...Generate phase space point and check against cuts.
176  loop=0
177  120 loop=loop+1
178  DO 130 i=1,2
179  IF(mint(140+i).NE.0) THEN
180 C...Pick x and Q2
181  x(i)=xmin(i)*(xmax(i)/xmin(i))**pyr(0)
182  q2(i)=q2min(i)*(q2max(i)/q2min(i))**pyr(0)
183 C...Cuts on internal consistency in x and Q2.
184  IF(q2(i).LT.x(i)**2*pms(i)/(1d0-x(i))) goto 120
185  IF(q2(i).GT.(1d0-x(i))*(vint(302)-2d0*pms(3-i))-
186  & (2d0-x(i)**2)*pms(i)/(1d0-x(i))) goto 120
187 C...Cuts on y and theta.
188  y(i)=(pmc(i)*x(i)+q2(i))/pmc(3)
189  IF(y(i).LT.ckin(71+2*i).OR.y(i).GT.ckin(72+2*i)) goto 120
190  rat=((1d0-x(i))*q2(i)-x(i)**2*pms(i))/
191  & ((1d0-x(i))**2*(vint(302)-2d0*pms(3-i)-2d0*pms(i)))
192  theta(i)=2d0*asin(sqrt(max(0d0,min(1d0,rat))))
193  IF(theta(i).LT.ckin(67+2*i)) goto 120
194  IF(ckin(68+2*i).GT.0d0.AND.theta(i).GT.ckin(68+2*i))
195  & goto 120
196 
197 C...Phi angle isotropic. Reconstruct pT.
198  phi(i)=paru(2)*pyr(0)
199  pt(i)=sqrt(((1d0-x(i))*pmc(i))**2/(4d0*vint(302))-
200  & pms(i))*sin(theta(i))
201 
202 C...Store info on variables selected, for documentation purposes.
203  vint(2+i)=-sqrt(q2(i))
204  vint(304+i)=x(i)
205  vint(306+i)=q2(i)
206  vint(308+i)=y(i)
207  vint(310+i)=theta(i)
208  vint(312+i)=phi(i)
209  ELSE
210  vint(304+i)=1d0
211  vint(306+i)=0d0
212  vint(308+i)=1d0
213  vint(310+i)=0d0
214  vint(312+i)=0d0
215  ENDIF
216  130 CONTINUE
217 
218 C...Cut on W combines info from two sides.
219  IF(mint(141).NE.0.AND.mint(142).NE.0) THEN
220  w2=-q2(1)-q2(2)+0.5d0*x(1)*pmc(1)*x(2)*pmc(2)/vint(302)-
221  & 2d0*pt(1)*pt(2)*cos(phi(1)-phi(2))+2d0*
222  & sqrt((0.5d0*x(1)*pmc(1)/vint(301))**2+q2(1)-pt(1)**2)*
223  & sqrt((0.5d0*x(2)*pmc(2)/vint(301))**2+q2(2)-pt(2)**2)
224  IF(w2.LT.w2min) goto 120
225  IF(ckin(78).GT.0d0.AND.w2.GT.ckin(78)**2) goto 120
226  pms1=-q2(1)
227  pms2=-q2(2)
228  ELSEIF(mint(141).NE.0) THEN
229  w2=(vint(302)+pms(1))*x(1)+pms(2)*(1d0-x(1))
230  pms1=-q2(1)
231  pms2=pms(2)
232  ELSEIF(mint(142).NE.0) THEN
233  w2=(vint(302)+pms(2))*x(2)+pms(1)*(1d0-x(2))
234  pms1=pms(1)
235  pms2=-q2(2)
236  ENDIF
237 
238 C...Store kinematics info for photon(s) in subsystem cm frame.
239  vint(2)=w2
240  vint(1)=sqrt(w2)
241  vint(291)=0d0
242  vint(292)=0d0
243  vint(293)=0.5d0*sqrt((w2-pms1-pms2)**2-4d0*pms1*pms2)/vint(1)
244  vint(294)=0.5d0*(w2+pms1-pms2)/vint(1)
245  vint(295)=sign(sqrt(abs(pms1)),pms1)
246  vint(296)=0d0
247  vint(297)=0d0
248  vint(298)=-vint(293)
249  vint(299)=0.5d0*(w2+pms2-pms1)/vint(1)
250  vint(300)=sign(sqrt(abs(pms2)),pms2)
251 
252 C...Assign weight for photon flux; different for transverse and
253 C...longitudinal photons. Flag incoming unresolved photon.
254  wtgaga=1d0
255  DO 140 i=1,2
256  IF(mint(140+i).NE.0) THEN
257  wtgaga=wtgaga*2d0*(paru(101)/paru(2))*
258  & log(xmax(i)/xmin(i))*log(q2max(i)/q2min(i))
259  IF(mstp(16).EQ.0) THEN
260  xy=x(i)
261  ELSE
262  wtgaga=wtgaga*x(i)/y(i)
263  xy=y(i)
264  ENDIF
265  IF(isub.EQ.132.OR.isub.EQ.134.OR.isub.EQ.136) THEN
266  IF((mint(11).EQ.22).and.
267  & (mint(12).EQ.2212.or.mint(12).EQ.2112)) THEN
268  wtgaga=wtgaga*(1d0/(1d0+(q2(i)/xy**2/ebeamenucl**2))*
269  & (1d0-xy-(q2(i)/4d0/ebeamenucl**2)))/
270  & q2(i)/xy**2/ebeamenucl*
271  & (ebeamenucl*xy-q2(i)/2d0/massp)*xy*q2(i)
272  ELSE
273  wtgaga=wtgaga*(1d0-xy)
274  ENDIF
275  ELSEIF(i.EQ.1.AND.(isub.EQ.139.OR.isub.EQ.140)) THEN
276  wtgaga=wtgaga*(1d0-xy)
277  ELSEIF(i.EQ.2.AND.(isub.EQ.138.OR.isub.EQ.140)) THEN
278  wtgaga=wtgaga*(1d0-xy)
279  ELSEIF((mint(11).EQ.22).and.
280  & (mint(12).EQ.2212.or.mint(12).EQ.2112)) THEN
281  wtgaga=wtgaga*(0.5d0*((ebeamenucl*xy-q2(i)/2d0/
282  & massp)/q2(i)/xy**2/ebeamenucl*
283  & (xy**2*(1d0-(2d0*masse**2/q2(i)))+
284  & (2d0/(1d0+(q2(i)/xy**2/ebeamenucl**2)))*
285  & (1d0-xy-(q2(i)/4d0/ebeamenucl**2))))*
286  & xy*q2(i))
287  ELSE
288  wtgaga=wtgaga*(0.5d0*(1d0+(1d0-xy)**2)-
289  & pms(i)*xy**2/q2(i))
290  ENDIF
291  IF(mint(106+i).EQ.0) mint(14+i)=22
292  ENDIF
293  140 CONTINUE
294  vint(319)=wtgaga
295  mint(143)=loop
296 
297 C...Update pTmin and cross section information.
298  IF(mstp(82).LE.1) THEN
299  ptmn=parp(81)*(vint(1)/parp(89))**parp(90)
300  ELSE
301  ptmn=parp(82)*(vint(1)/parp(89))**parp(90)
302  ENDIF
303  vint(149)=4d0*ptmn**2/vint(2)
304  vint(154)=ptmn
305  CALL pyxtot
306 
307 C...ECA...Routine modified for rad-corrections
308 C...Generate photons inside leptons and
309 C...calculate photon flux weight factors.
310  ELSEIF(igaga.EQ.5) THEN
311  write(*,*).EQ."In pygaga with IGAGA5"
312  isub=mint(1)
313  mint(15)=0
314  mint(16)=0
315 
316 C...Generate phase space point and check against cuts.
317  loop=0
318  121 loop=loop+1
319  DO 131 i=1,2
320  IF(mint(140+i).NE.0) THEN
321 C...Pick x and Q2
322  mint(199)=0
323  geny=mcset_ymin*(mcset_ymax/mcset_ymin)**pyr(0)
324  genq2=mcset_q2min*(mcset_q2max/mcset_q2min)**pyr(0)
325  genx = genq2/geny/(4.*ebeam*pbeam)
326  genw2 = massp**2.+(genq2*(1./genx-1.))
327  gennu=geny*ebeamenucl
328  geneprim = ebeamenucl - gennu
329  write(*,*) geny, genq2, genx, gennu, genphi, ebeamenucl
330 C....Check to have sensible ranges for variables
331  minq2 = pms(1) * geny**2. / (1.- geny)
332  if (genq2.lt.minq2) then
333  goto 121
334  endif
335 C... check x and Q2 go toghether
336  if (genq2.gt.(2.*gennu*massp)) then
337  goto 121
338  endif
339 C temp = ((genQ2-minq2)/(2D0*ebeamEnucl*geneprim))-1D0
340  temp=(genq2/(2.*ebeamenucl*geneprim))-1.
341  if ((temp.le.1.00).and.(temp.ge.-1.0)) then
342  etheta = acos(temp)
343  genthe = sngl(etheta)
344  genphi=paru(2)*pyr(0)
345  phi(i)=genphi
346  ppt=tan(etheta)
347  ppx=ppt*cos(phi(i))
348  ppy=ppt*sin(phi(i))
349  else
350  goto 121
351  endif
352 
353  if ((genw2.lt.ckin(77)**2).or.
354  & (ckin(78).gt.0.and.genw2.gt.ckin(78)**2)) then
355  goto 121
356  endif
357  write(*,*) geny, genq2, genx, gennu, genphi, ebeamenucl
358 
359  ntries=0
360  122 if (qedrad.eq.1) then
361  write(*,*) geny, genq2, genx, gennu, genphi, ebeamenucl
362  call radgen_event
363  endif
364  if (qedrad.eq.0) then
365  y(i)=geny
366  q2(i)=genq2
367  elseif ((mcradcor_ebrems.eq.mcradcor_ebrems).and.
368  & (mcradcor_thetabrems.eq.mcradcor_thetabrems)) then
369  y(i)=dble(mcradcor_nutrue)/dble(mcset_enebeam)
370  q2(i)=dble(mcradcor_q2true)
371  write(*,*) geny, genq2, genx, gennu, genphi, ebeamenucl
372  write(*,*) mcradcor_nutrue, mcset_enebeam, mcradcor_q2true
373  else
374  write(*,*)"I go to 122 again"
375  write(*,*) mcradcor_thetabrems,mcradcor_ebrems,
376  & mcevent_ievent
377  goto 122
378  endif
379 C write(*,*) 'geny, genq2, genx, genW2, genNu',
380 C & geny, genq2, genx, genW2, genNu
381  x(i)=((pmc(3)*y(i))-q2(i))/pmc(i)
382 C P.L. ...An event with W^2_T<4will be generated new by RADGEN at the
383 C ...same kinematic point, the number of tries needed by RADGEN is
384 C ...counted and saved in the variable rcweight!
385  IF (qedrad.ne.0) then
386 C IF(mcradcor_cType.eq.'qela') then
387 C GOTO 122
388 C ENDIF
389 C IF(mcradcor_cType.eq.'elas') then
390 C GOTO 122
391 C ENDIF
392  IF(dble(mcradcor_w2true).LT.
393  & (ckin(77)**2-1.d-4*abs(ckin(77)**2))) THEN
394  mint(199)=mint(199)+1
395  goto 122
396  ENDIF
397  ENDIF
398  ntries=ntries+1
399  IF(ntries.ge.20) goto 121
400 
401 C ...... New try to implement weights directly into Pythia
402  sigobs=0.0d0
403  sigtrue=0.0d0
404  rccorr=1.0d0
405  if (qedrad.eq.1) then
406  call mkf2(genq2,genx,
407  + mcset_tara,mcset_tarz,py6f2,py6f1)
408  sigobs=pyth_xsec(geny,genx,genq2,py6f1, py6f2)
409  IF(mcradcor_ebrems.eq.0) THEN
410  IF (sig1g.gt.0.d0) then
411  rccorr=(tbor+tine)/sig1g/(dble(mint(199))+1.0d0)
412  ELSE
413  rccorr=0.d0
414  ENDIF
415  ELSEIF(mcradcor_ebrems.gt.0) THEN
416  call mkf2(dble(mcradcor_q2true),dble(mcradcor_xtrue),
417  + mcset_tara,mcset_tarz,py6f2,py6f1)
418  sigtrue=pyth_xsec(dble( mcradcor_ytrue), dble(mcradcor_xtrue),
419  + dble(mcradcor_q2true), py6f1, py6f2)
420  IF ((sig1g.gt.0.d0).and.(sigtrue.gt.0.d0)) then
421  rccorr=(tbor+tine)/sig1g*sigobs/sigtrue/
422  & (dble(mint(199))+1.0d0)
423  ELSE
424  rccorr=0.d0
425  ENDIF
426  ENDIF
427  ENDIF
428  IF(x(i).GT.(xmax(i)+1.d-4*abs(xmax(i)))) THEN
429  goto 121
430  ENDIF
431 C...Cuts on internal consistency in x and Q2.
432  IF(q2(i).LT.x(i)**2*pms(i)/(1d0-x(i))) then
433  goto 121
434  endif
435  IF(q2(i).GT.(1d0-x(i))*(vint(302)-2d0*pms(3-i))-
436  & (2d0-x(i)**2)*pms(i)/(1d0-x(i))) THEN
437  goto 121
438  ENDIF
439 C...Cuts on y and theta.
440  IF(y(i).LT.ckin(71+2*i).OR.y(i).GT.ckin(72+2*i)) THEN
441  goto 121
442  ENDIF
443  rat=((1d0-x(i))*q2(i)-x(i)**2*pms(i))/
444  & ((1d0-x(i))**2*(vint(302)-2d0*pms(3-i)-2d0*pms(i)))
445  theta(i)=2d0*asin(sqrt(max(0d0,min(1d0,rat))))
446  IF(theta(i).LT.ckin(67+2*i)) THEN
447  goto 121
448  ENDIF
449  IF(ckin(68+2*i).GT.0d0.AND.theta(i).GT.ckin(68+2*i))
450  & goto 121
451 
452 C...Phi angle isotropic. Reconstruct pT.
453 C PT(I)=SQRT(((1D0-X(I))*PMC(I))**2/(4D0*VINT(302))-
454 C & PMS(I))*SIN(THETA(I))
455  temp=((1d0-x(i))*pmc(i))**2/(4d0*vint(302))-pms(i)
456  pt(i)=(sqrt(temp))*sin(theta(i))
457 C ... try 'new' phi
458  IF ((qedrad.ne.0).and.(mcradcor_ebrems.gt.0)) then
459  emom=sqrt(geneprim**2-masse**2)
460  phi(i)=atan2((emom*ppy+dplabg(2)),(emom*ppx+dplabg(1)))
461  IF (phi(i).lt.0) THEN
462  phi(i)=phi(i)+paru(2)
463  ENDIF
464  ENDIF
465 C...Store info on variables selected, for documentation purposes.
466  vint(2+i)=-sqrt(q2(i))
467  vint(304+i)=x(i)
468  vint(306+i)=q2(i)
469  vint(308+i)=y(i)
470  vint(310+i)=theta(i)
471  vint(312+i)=phi(i)
472  ELSE
473  vint(304+i)=1d0
474  vint(306+i)=0d0
475  vint(308+i)=1d0
476  vint(310+i)=0d0
477  vint(312+i)=0d0
478  ENDIF
479  131 CONTINUE
480 
481 C...Cut on W combines info from two sides.
482  IF(mint(141).NE.0.AND.mint(142).NE.0) THEN
483  w2=-q2(1)-q2(2)+0.5d0*x(1)*pmc(1)*x(2)*pmc(2)/vint(302)-
484  & 2d0*pt(1)*pt(2)*cos(phi(1)-phi(2))+2d0*
485  & sqrt((0.5d0*x(1)*pmc(1)/vint(301))**2+q2(1)-pt(1)**2)*
486  & sqrt((0.5d0*x(2)*pmc(2)/vint(301))**2+q2(2)-pt(2)**2)
487  IF(w2.LT.w2min) THEN
488  goto 121
489  ENDIF
490  IF(ckin(78).GT.0d0.AND.w2.GT.ckin(78)**2) goto 121
491  pms1=-q2(1)
492  pms2=-q2(2)
493  ELSEIF(mint(141).NE.0) THEN
494  w2=(vint(302)+pms(1))*x(1)+pms(2)*(1d0-x(1))
495  pms1=-q2(1)
496  pms2=pms(2)
497  ELSEIF(mint(142).NE.0) THEN
498  w2=(vint(302)+pms(2))*x(2)+pms(1)*(1d0-x(2))
499  pms1=pms(1)
500  pms2=-q2(2)
501  ENDIF
502 
503 C...Store kinematics info for photon(s) in subsystem cm frame.
504  vint(2)=w2
505  vint(1)=sqrt(w2)
506  vint(291)=0d0
507  vint(292)=0d0
508  vint(293)=0.5d0*sqrt((w2-pms1-pms2)**2-4d0*pms1*pms2)/vint(1)
509  vint(294)=0.5d0*(w2+pms1-pms2)/vint(1)
510  vint(295)=sign(sqrt(abs(pms1)),pms1)
511  vint(296)=0d0
512  vint(297)=0d0
513  vint(298)=-vint(293)
514  vint(299)=0.5d0*(w2+pms2-pms1)/vint(1)
515  vint(300)=sign(sqrt(abs(pms2)),pms2)
516 
517 C...Assign weight for photon flux; different for transverse and
518 C...longitudinal photons. Flag incoming unresolved photon.
519  wtgaga=1d0
520  DO 141 i=1,2
521  IF(mint(140+i).NE.0) THEN
522  wtgaga=wtgaga*2d0*(paru(101)/paru(2))*
523  & (log(real(mcset_ymax))-log(real(mcset_ymin)))*
524  & (log(real(mcset_q2max))-log(real(mcset_q2min)))
525  xy=y(i)
526  IF(isub.EQ.132.OR.isub.EQ.134.OR.isub.EQ.136) THEN
527  IF((mint(11).EQ.22).and.
528  & (mint(12).EQ.2212.or.mint(12).EQ.2112)) THEN
529  xxy=xy*ebeamenucl/ebeamenucl
530  wtgaga=wtgaga*(1d0/(1d0+(q2(i)/xxy**2/ebeamenucl**2))*
531  & (1d0-xxy-(q2(i)/4d0/ebeamenucl**2)))/
532  & q2(i)/xxy**2/ebeamenucl*
533  & (ebeamenucl*xxy-q2(i)/2d0/massp)*xxy*q2(i)
534  ELSE
535  wtgaga=wtgaga*(1d0-xy)
536  ENDIF
537  ELSEIF(i.EQ.1.AND.(isub.EQ.139.OR.isub.EQ.140)) THEN
538  wtgaga=wtgaga*(1d0-xy)
539  ELSEIF(i.EQ.2.AND.(isub.EQ.138.OR.isub.EQ.140)) THEN
540  wtgaga=wtgaga*(1d0-xy)
541  ELSEIF((mint(11).EQ.22).and.
542  & (mint(12).EQ.2212.or.mint(12).EQ.2112)) THEN
543  xxy=xy*ebeamenucl/ebeamenucl
544  wtgaga=wtgaga*(0.5d0*((ebeamenucl*xxy-q2(i)/2d0/
545  & massp)/q2(i)/xxy**2/ebeamenucl*
546  & (xxy**2*(1d0-(2d0*masse**2/q2(i)))+
547  & (2d0/(1d0+(q2(i)/xxy**2/ebeamenucl**2)))*
548  & (1d0-xxy-(q2(i)/4d0/ebeamenucl**2))))*xxy*q2(i))
549  ELSE
550  wtgaga=wtgaga*(0.5d0*(1d0+(1d0-xy)**2)-
551  & pms(i)*xy**2/q2(i))
552  ENDIF
553  IF(mint(106+i).EQ.0) mint(14+i)=22
554  ENDIF
555  141 CONTINUE
556  wtgaga=wtgaga*rccorr
557  vint(319)=wtgaga
558  mint(143)=loop
559 C...Update pTmin and cross section information.
560  IF(mstp(82).LE.1) THEN
561  ptmn=parp(81)*(vint(1)/parp(89))**parp(90)
562  ELSE
563  ptmn=parp(82)*(vint(1)/parp(89))**parp(90)
564  ENDIF
565  vint(149)=4d0*ptmn**2/vint(2)
566  vint(154)=ptmn
567  CALL pyxtot
568 
569 C...Reconstruct kinematics of photons inside leptons.
570  ELSEIF(igaga.EQ.4) THEN
571 
572 C...Make place for incoming particles and scattered leptons.
573  move=3
574  IF(mint(141).NE.0.AND.mint(142).NE.0) move=4
575  mint(4)=mint(4)+move
576  DO 160 i=mint(84)-move,mint(83)+1,-1
577  IF(k(i,1).EQ.21) THEN
578  DO 150 j=1,5
579  k(i+move,j)=k(i,j)
580  p(i+move,j)=p(i,j)
581  v(i+move,j)=v(i,j)
582  150 CONTINUE
583  IF(k(i,3).GT.mint(83).AND.k(i,3).LE.mint(84))
584  & k(i+move,3)=k(i,3)+move
585  IF(k(i,4).GT.mint(83).AND.k(i,4).LE.mint(84))
586  & k(i+move,4)=k(i,4)+move
587  IF(k(i,5).GT.mint(83).AND.k(i,5).LE.mint(84))
588  & k(i+move,5)=k(i,5)+move
589  ENDIF
590  160 CONTINUE
591  DO 170 i=mint(84)+1,n
592  IF(k(i,3).GT.mint(83).AND.k(i,3).LE.mint(84))
593  & k(i,3)=k(i,3)+move
594  170 CONTINUE
595 
596 C...Fill in incoming particles.
597  DO 190 i=mint(83)+1,mint(83)+move
598  DO 180 j=1,5
599  k(i,j)=0
600  p(i,j)=0d0
601  v(i,j)=0d0
602  180 CONTINUE
603  190 CONTINUE
604  DO 200 i=1,2
605  k(mint(83)+i,1)=21
606  IF(mint(140+i).NE.0) THEN
607  k(mint(83)+i,2)=mint(140+i)
608  p(mint(83)+i,5)=vint(302+i)
609  ELSE
610  k(mint(83)+i,2)=mint(10+i)
611  p(mint(83)+i,5)=vint(2+i)
612  ENDIF
613  p(mint(83)+i,3)=0.5d0*sqrt((pmc(3)**2-4d0*pms(1)*pms(2))/
614  & vint(302))*(-1d0)**(i+1)
615  p(mint(83)+i,4)=0.5d0*pmc(i)/vint(301)
616  200 CONTINUE
617 
618 C...New mother-daughter relations in documentation section.
619  IF(mint(141).NE.0.AND.mint(142).NE.0) THEN
620  k(mint(83)+1,4)=mint(83)+3
621  k(mint(83)+1,5)=mint(83)+5
622  k(mint(83)+2,4)=mint(83)+4
623  k(mint(83)+2,5)=mint(83)+6
624  k(mint(83)+3,3)=mint(83)+1
625  k(mint(83)+5,3)=mint(83)+1
626  k(mint(83)+4,3)=mint(83)+2
627  k(mint(83)+6,3)=mint(83)+2
628  ELSEIF(mint(141).NE.0) THEN
629  k(mint(83)+1,4)=mint(83)+3
630  k(mint(83)+1,5)=mint(83)+4
631  k(mint(83)+2,4)=mint(83)+5
632  k(mint(83)+3,3)=mint(83)+1
633  k(mint(83)+4,3)=mint(83)+1
634  k(mint(83)+5,3)=mint(83)+2
635  ELSEIF(mint(142).NE.0) THEN
636  k(mint(83)+1,4)=mint(83)+4
637  k(mint(83)+2,4)=mint(83)+3
638  k(mint(83)+2,5)=mint(83)+5
639  k(mint(83)+3,3)=mint(83)+2
640  k(mint(83)+4,3)=mint(83)+1
641  k(mint(83)+5,3)=mint(83)+2
642  ENDIF
643 
644 C...Fill scattered lepton(s).
645  DO 210 i=1,2
646  IF(mint(140+i).NE.0) THEN
647  lsc=mint(83)+min(i+2,move)
648  k(lsc,1)=21
649  k(lsc,2)=mint(140+i)
650  p(lsc,1)=pt(i)*cos(phi(i))
651  p(lsc,2)=pt(i)*sin(phi(i))
652  p(lsc,4)=(1d0-x(i))*p(mint(83)+i,4)
653  p(lsc,3)=sqrt(p(lsc,4)**2-pms(i))*cos(theta(i))*
654  & (-1d0)**(i-1)
655  p(lsc,5)=vint(302+i)
656  ENDIF
657  210 CONTINUE
658 
659 C...Find incoming four-vectors to subprocess.
660  k(n+1,1)=21
661  IF(mint(141).NE.0) THEN
662  DO 220 j=1,4
663  p(n+1,j)=p(mint(83)+1,j)-p(mint(83)+3,j)
664  220 CONTINUE
665  ELSE
666  DO 230 j=1,4
667  p(n+1,j)=p(mint(83)+1,j)
668  230 CONTINUE
669  ENDIF
670  k(n+2,1)=21
671  IF(mint(142).NE.0) THEN
672  DO 240 j=1,4
673  p(n+2,j)=p(mint(83)+2,j)-p(mint(83)+move,j)
674  240 CONTINUE
675  ELSE
676  DO 250 j=1,4
677  p(n+2,j)=p(mint(83)+2,j)
678  250 CONTINUE
679  ENDIF
680 
681 C...Define boost and rotation between hadronic subsystem and
682 C...collision rest frame; boost hadronic subsystem to this frame.
683  DO 260 j=1,3
684  beta(j)=(p(n+1,j)+p(n+2,j))/(p(n+1,4)+p(n+2,4))
685  260 CONTINUE
686  CALL pyrobo(n+1,n+2,0d0,0d0,-beta(1),-beta(2),-beta(3))
687  bphi=pyangl(p(n+1,1),p(n+1,2))
688  CALL pyrobo(n+1,n+2,0d0,-bphi,0d0,0d0,0d0)
689  btheta=pyangl(p(n+1,3),p(n+1,1))
690  CALL pyrobo(mint(83)+move+1,n,btheta,bphi,beta(1),beta(2),
691  & beta(3))
692 
693 C...Add on scattered leptons to final state.
694  DO 280 i=1,2
695  IF(mint(140+i).NE.0) THEN
696  lsc=mint(83)+min(i+2,move)
697  n=n+1
698  DO 270 j=1,5
699  k(n,j)=k(lsc,j)
700  p(n,j)=p(lsc,j)
701  v(n,j)=v(lsc,j)
702  270 CONTINUE
703  k(n,1)=1
704  k(n,3)=lsc
705  ENDIF
706  280 CONTINUE
707  ENDIF
708 
709  290 CONTINUE
710  RETURN
711  END