pythiaeRHIC/pythia/pygaga.f

0001 C*********************************************************************
0002
0003 C...PYGAGA
0004 C...For lepton beams it gives photon-hadron or photon-photon systems
0005 C...to be treated with the ordinary machinery and combines this with a
0006 C...description of the lepton -> lepton + photon branching.
0007
0008       SUBROUTINE PYGAGA(IGAGA,WTGAGA)
0009
0010 C...Double precision and integer declarations.
0011       IMPLICIT DOUBLE PRECISION(A-H, O-Z)
0012       IMPLICIT INTEGER(I-N)
0013       INTEGER PYK,PYCHGE,PYCOMP
0014 C...Commonblocks.
0015       COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5)
0016       COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200)
0017       COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4)
0018       COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200)
0019       COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)
0020       COMMON/PYINT1/MINT(400),VINT(400)
0021       COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3)
0022       DOUBLE PRECISION minq2, rccorr, pyth_xsec, temp, etheta
0023       DOUBLE PRECISION ppt, ppx, ppy
0024       include "mcRadCor.inc"
0025       include "mc_set.inc"
0026       include "radgen.inc"
0027       include "phiout.inc"
0028       include "py6strf.inc"
0029
0030       SAVE /PYJETS/,/PYDAT1/,/PYDAT2/,/PYSUBS/,/PYPARS/,/PYINT1/,
0031      &/PYINT5/
0032 C...Local variables and data statement.
0033       DIMENSION PMS(2),XMIN(2),XMAX(2),Q2MIN(2),Q2MAX(2),PMC(3),
0034      &X(2),Q2(2),Y(2),THETA(2),PHI(2),PT(2),BETA(3)
0035       SAVE PMS,XMIN,XMAX,Q2MIN,Q2MAX,PMC,X,Q2,THETA,PHI,PT,W2MIN,
0036      & YMIN,YMAX
0037       DATA EPS/1D-4/
0038
0039 C...Initialize generation of photons inside leptons.
0040       IF(IGAGA.EQ.1) THEN
0041
0042 C...Save quantities on incoming lepton system.
0043         VINT(301)=VINT(1)
0044         VINT(302)=VINT(2)
0045         PMS(1)=VINT(303)**2
0046         IF(MINT(141).EQ.0) PMS(1)=SIGN(VINT(3)**2,VINT(3))
0047         PMS(2)=VINT(304)**2
0048         IF(MINT(142).EQ.0) PMS(2)=SIGN(VINT(4)**2,VINT(4))
0049         PMC(3)=VINT(302)-PMS(1)-PMS(2)
0050         W2MIN=MAX(CKIN(77),2D0*CKIN(3),2D0*CKIN(5))**2
0051
0052 C...Calculate range of x and Q2 values allowed in generation.
0053         DO 100 I=1,2
0054           PMC(I)=VINT(302)+PMS(I)-PMS(3-I)
0055           IF(MINT(140+I).NE.0) THEN
0056             XMIN(I)=MAX(CKIN(59+2*I),EPS)
0057             XMAX(I)=MIN(CKIN(60+2*I),1D0-2D0*VINT(301)*SQRT(PMS(I))/
0058      &      PMC(I),1D0-EPS)
0059             YMIN=MAX(CKIN(71+2*I),EPS)
0060             YMAX=MIN(CKIN(72+2*I),1D0-EPS)
0061             IF(CKIN(64+2*I).GT.0D0) XMIN(I)=MAX(XMIN(I),
0062      &      (YMIN*PMC(3)-CKIN(64+2*I))/PMC(I))
0063             XMAX(I)=MIN(XMAX(I),(YMAX*PMC(3)-CKIN(63+2*I))/PMC(I))
0064             THEMIN=MAX(CKIN(67+2*I),0D0)
0065             THEMAX=MIN(CKIN(68+2*I),PARU(1))
0066             IF(CKIN(68+2*I).LT.0D0) THEMAX=PARU(1)
0067             Q2MIN(I)=MAX(CKIN(63+2*I),XMIN(I)**2*PMS(I)/(1D0-XMIN(I))+
0068      &      ((1D0-XMAX(I))*(VINT(302)-2D0*PMS(3-I))-
0069      &      2D0*PMS(I)/(1D0-XMAX(I)))*SIN(THEMIN/2D0)**2,0D0)
0070             Q2MAX(I)=XMAX(I)**2*PMS(I)/(1D0-XMAX(I))+
0071      &      ((1D0-XMIN(I))*(VINT(302)-2D0*PMS(3-I))-
0072      &      2D0*PMS(I)/(1D0-XMIN(I)))*SIN(THEMAX/2D0)**2
0073             IF(CKIN(64+2*I).GT.0D0) Q2MAX(I)=MIN(CKIN(64+2*I),Q2MAX(I))
0074 C...W limits when lepton on one side only.
0075             IF(MINT(143-I).EQ.0) THEN
0076               XMIN(I)=MAX(XMIN(I),(W2MIN-PMS(3-I))/PMC(I))
0077               IF(CKIN(78).GT.0D0) XMAX(I)=MIN(XMAX(I),
0078      &        (CKIN(78)**2-PMS(3-I))/PMC(I))
0079             ENDIF
0080           ENDIF
0081   100   CONTINUE
0082
0083 C...W limits when lepton on both sides.
0084         IF(MINT(141).NE.0.AND.MINT(142).NE.0) THEN
0085           IF(CKIN(78).GT.0D0) XMAX(1)=MIN(XMAX(1),
0086      &    (CKIN(78)**2+PMC(3)-PMC(2)*XMIN(2))/PMC(1))
0087           IF(CKIN(78).GT.0D0) XMAX(2)=MIN(XMAX(2),
0088      &    (CKIN(78)**2+PMC(3)-PMC(1)*XMIN(1))/PMC(2))
0089           IF(IABS(MINT(141)).NE.IABS(MINT(142))) THEN
0090             XMIN(1)=MAX(XMIN(1),(PMS(1)-PMS(2)+VINT(302)*(W2MIN-
0091      &      PMS(1)-PMS(2))/(PMC(2)*XMAX(2)+PMS(1)-PMS(2)))/PMC(1))
0092             XMIN(2)=MAX(XMIN(2),(PMS(2)-PMS(1)+VINT(302)*(W2MIN-
0093      &      PMS(1)-PMS(2))/(PMC(1)*XMAX(1)+PMS(2)-PMS(1)))/PMC(2))
0094           ELSE
0095             XMIN(1)=MAX(XMIN(1),W2MIN/(VINT(302)*XMAX(2)))
0096             XMIN(2)=MAX(XMIN(2),W2MIN/(VINT(302)*XMAX(1)))
0097           ENDIF
0098         ENDIF
0099
0100 C...Q2 and W values and photon flux weight factors for initialization.
0101       ELSEIF(IGAGA.EQ.2) THEN
0102         ISUB=MINT(1)
0103         MINT(15)=0
0104         MINT(16)=0
0105
0106 C...W value for photon on one or both sides, and for processes
0107 C...with gamma-gamma cross section peaked at small shat.
0108         IF(MINT(141).NE.0.AND.MINT(142).EQ.0) THEN
0109           VINT(2)=VINT(302)+PMS(1)-PMC(1)*(1D0-XMAX(1))
0110         ELSEIF(MINT(141).EQ.0.AND.MINT(142).NE.0) THEN
0111           VINT(2)=VINT(302)+PMS(2)-PMC(2)*(1D0-XMAX(2))
0112         ELSEIF(ISUB.GE.137.AND.ISUB.LE.140) THEN
0113           VINT(2)=MAX(CKIN(77)**2,12D0*MAX(CKIN(3),CKIN(5))**2)
0114           IF(CKIN(78).GT.0D0) VINT(2)=MIN(VINT(2),CKIN(78)**2)
0115         ELSE
0116           VINT(2)=XMAX(1)*XMAX(2)*VINT(302)
0117           IF(CKIN(78).GT.0D0) VINT(2)=MIN(VINT(2),CKIN(78)**2)
0118         ENDIF
0119         VINT(1)=SQRT(MAX(0D0,VINT(2)))
0120
0121 C...Upper estimate of photon flux weight factor.
0122 C...Initialization Q2 scale. Flag incoming unresolved photon.
0123         WTGAGA=1D0
0124         DO 110 I=1,2
0125           IF(MINT(140+I).NE.0) THEN
0126            IF(MSTP(199).EQ.1) then
0127             WTGAGA=WTGAGA*2D0*(PARU(101)/PARU(2))*
0128      &      (LOG(real(mcSet_YMax/mcSet_YMin)))*
0129      &      (LOG(real(mcSet_Q2Max/mcSet_Q2Min)))
0130            ELSE
0131             WTGAGA=WTGAGA*2D0*(PARU(101)/PARU(2))*
0132      &      LOG(XMAX(I)/XMIN(I))*LOG(Q2MAX(I)/Q2MIN(I))
0133            ENDIF
0134             IF(ISUB.EQ.99.AND.MINT(106+I).EQ.4.AND.MINT(109-I).EQ.3)
0135      &      THEN
0136               Q2INIT=5D0+Q2MIN(3-I)
0137             ELSEIF(ISUB.EQ.99.AND.MINT(106+I).EQ.4) THEN
0138               Q2INIT=PMAS(PYCOMP(113),1)**2+Q2MIN(3-I)
0139             ELSEIF(ISUB.EQ.132.OR.ISUB.EQ.134.OR.ISUB.EQ.136) THEN
0140               Q2INIT=MAX(CKIN(1),2D0*CKIN(3),2D0*CKIN(5))**2/3D0
0141             ELSEIF((ISUB.EQ.138.AND.I.EQ.2).OR.
0142      &      (ISUB.EQ.139.AND.I.EQ.1)) THEN
0143               Q2INIT=VINT(2)/3D0
0144             ELSEIF(ISUB.EQ.140) THEN
0145               Q2INIT=VINT(2)/2D0
0146             ELSE
0147               Q2INIT=Q2MIN(I)
0148             ENDIF
0149             VINT(2+I)=-SQRT(MAX(Q2MIN(I),MIN(Q2MAX(I),Q2INIT)))
0150             IF(MSTP(14).EQ.0.OR.(ISUB.GE.131.AND.ISUB.LE.140))
0151      &      MINT(14+I)=22
0152             VINT(306+I)=VINT(2+I)**2
0153           ENDIF
0154   110   CONTINUE
0155         VINT(320)=WTGAGA
0156
0157 C...Update pTmin and cross section information.
0158         IF(MSTP(82).LE.1) THEN
0159           PTMN=PARP(81)*(VINT(1)/PARP(89))**PARP(90)
0160         ELSE
0161           PTMN=PARP(82)*(VINT(1)/PARP(89))**PARP(90)
0162         ENDIF
0163         VINT(149)=4D0*PTMN**2/VINT(2)
0164         VINT(154)=PTMN
0165         CALL PYXTOT
0166         VINT(318)=VINT(317)
0167
0168 C...Generate photons inside leptons and
0169 C...calculate photon flux weight factors.
0170       ELSEIF(IGAGA.EQ.3) THEN
0171         ISUB=MINT(1)
0172         MINT(15)=0
0173         MINT(16)=0
0174
0175 C...Generate phase space point and check against cuts.
0176         LOOP=0
0177   120   LOOP=LOOP+1
0178         DO 130 I=1,2
0179           IF(MINT(140+I).NE.0) THEN
0180 C...Pick x and Q2
0181             X(I)=XMIN(I)*(XMAX(I)/XMIN(I))**PYR(0)
0182             Q2(I)=Q2MIN(I)*(Q2MAX(I)/Q2MIN(I))**PYR(0)
0183 C...Cuts on internal consistency in x and Q2.
0184             IF(Q2(I).LT.X(I)**2*PMS(I)/(1D0-X(I))) GOTO 120
0185             IF(Q2(I).GT.(1D0-X(I))*(VINT(302)-2D0*PMS(3-I))-
0186      &      (2D0-X(I)**2)*PMS(I)/(1D0-X(I))) GOTO 120
0187 C...Cuts on y and theta.
0188             Y(I)=(PMC(I)*X(I)+Q2(I))/PMC(3)
0189             IF(Y(I).LT.CKIN(71+2*I).OR.Y(I).GT.CKIN(72+2*I)) GOTO 120
0190             RAT=((1D0-X(I))*Q2(I)-X(I)**2*PMS(I))/
0191      &      ((1D0-X(I))**2*(VINT(302)-2D0*PMS(3-I)-2D0*PMS(I)))
0192             THETA(I)=2D0*ASIN(SQRT(MAX(0D0,MIN(1D0,RAT))))
0193             IF(THETA(I).LT.CKIN(67+2*I)) GOTO 120
0194             IF(CKIN(68+2*I).GT.0D0.AND.THETA(I).GT.CKIN(68+2*I))
0195      &      GOTO 120
0196
0197 C...Phi angle isotropic. Reconstruct pT.
0198             PHI(I)=PARU(2)*PYR(0)
0199             PT(I)=SQRT(((1D0-X(I))*PMC(I))**2/(4D0*VINT(302))-
0200      &      PMS(I))*SIN(THETA(I))
0201
0202 C...Store info on variables selected, for documentation purposes.
0203             VINT(2+I)=-SQRT(Q2(I))
0204             VINT(304+I)=X(I)
0205             VINT(306+I)=Q2(I)
0206             VINT(308+I)=Y(I)
0207             VINT(310+I)=THETA(I)
0208             VINT(312+I)=PHI(I)
0209           ELSE
0210             VINT(304+I)=1D0
0211             VINT(306+I)=0D0
0212             VINT(308+I)=1D0
0213             VINT(310+I)=0D0
0214             VINT(312+I)=0D0
0215           ENDIF
0216   130   CONTINUE
0217
0218 C...Cut on W combines info from two sides.
0219         IF(MINT(141).NE.0.AND.MINT(142).NE.0) THEN
0220           W2=-Q2(1)-Q2(2)+0.5D0*X(1)*PMC(1)*X(2)*PMC(2)/VINT(302)-
0221      &    2D0*PT(1)*PT(2)*COS(PHI(1)-PHI(2))+2D0*
0222      &    SQRT((0.5D0*X(1)*PMC(1)/VINT(301))**2+Q2(1)-PT(1)**2)*
0223      &    SQRT((0.5D0*X(2)*PMC(2)/VINT(301))**2+Q2(2)-PT(2)**2)
0224           IF(W2.LT.W2MIN) GOTO 120
0225           IF(CKIN(78).GT.0D0.AND.W2.GT.CKIN(78)**2) GOTO 120
0226           PMS1=-Q2(1)
0227           PMS2=-Q2(2)
0228         ELSEIF(MINT(141).NE.0) THEN
0229           W2=(VINT(302)+PMS(1))*X(1)+PMS(2)*(1D0-X(1))
0230           PMS1=-Q2(1)
0231           PMS2=PMS(2)
0232         ELSEIF(MINT(142).NE.0) THEN
0233           W2=(VINT(302)+PMS(2))*X(2)+PMS(1)*(1D0-X(2))
0234           PMS1=PMS(1)
0235           PMS2=-Q2(2)
0236         ENDIF
0237
0238 C...Store kinematics info for photon(s) in subsystem cm frame.
0239         VINT(2)=W2
0240         VINT(1)=SQRT(W2)
0241         VINT(291)=0D0
0242         VINT(292)=0D0
0243         VINT(293)=0.5D0*SQRT((W2-PMS1-PMS2)**2-4D0*PMS1*PMS2)/VINT(1)
0244         VINT(294)=0.5D0*(W2+PMS1-PMS2)/VINT(1)
0245         VINT(295)=SIGN(SQRT(ABS(PMS1)),PMS1)
0246         VINT(296)=0D0
0247         VINT(297)=0D0
0248         VINT(298)=-VINT(293)
0249         VINT(299)=0.5D0*(W2+PMS2-PMS1)/VINT(1)
0250         VINT(300)=SIGN(SQRT(ABS(PMS2)),PMS2)
0251
0252 C...Assign weight for photon flux; different for transverse and
0253 C...longitudinal photons. Flag incoming unresolved photon.
0254         WTGAGA=1D0
0255         DO 140 I=1,2
0256           IF(MINT(140+I).NE.0) THEN
0257             WTGAGA=WTGAGA*2D0*(PARU(101)/PARU(2))*
0258      &      LOG(XMAX(I)/XMIN(I))*LOG(Q2MAX(I)/Q2MIN(I))
0259             IF(MSTP(16).EQ.0) THEN
0260               XY=X(I)
0261             ELSE
0262               WTGAGA=WTGAGA*X(I)/Y(I)
0263               XY=Y(I)
0264             ENDIF
0265             IF(ISUB.EQ.132.OR.ISUB.EQ.134.OR.ISUB.EQ.136) THEN
0266               IF((MINT(11).EQ.22).and.
0267      &           (MINT(12).EQ.2212.or.MINT(12).EQ.2112)) THEN
0268                WTGAGA=WTGAGA*(1D0/(1D0+(Q2(I)/XY**2/ebeamEnucl**2))*
0269      &               (1D0-XY-(Q2(I)/4D0/ebeamEnucl**2)))/
0270      &               Q2(I)/XY**2/ebeamEnucl*
0271      &               (ebeamEnucl*XY-Q2(I)/2D0/massp)*XY*Q2(I)
0272               ELSE
0273                 WTGAGA=WTGAGA*(1D0-XY)
0274               ENDIF
0275             ELSEIF(I.EQ.1.AND.(ISUB.EQ.139.OR.ISUB.EQ.140)) THEN
0276               WTGAGA=WTGAGA*(1D0-XY)
0277             ELSEIF(I.EQ.2.AND.(ISUB.EQ.138.OR.ISUB.EQ.140)) THEN
0278               WTGAGA=WTGAGA*(1D0-XY)
0279             ELSEIF((MINT(11).EQ.22).and.
0280      &             (MINT(12).EQ.2212.or.MINT(12).EQ.2112)) THEN
0281               WTGAGA=WTGAGA*(0.5D0*((ebeamEnucl*XY-Q2(I)/2D0/
0282      &               massp)/Q2(I)/XY**2/ebeamEnucl*
0283      &               (XY**2*(1D0-(2D0*masse**2/Q2(I)))+
0284      &               (2D0/(1D0+(Q2(I)/XY**2/ebeamEnucl**2)))*
0285      &               (1D0-XY-(Q2(I)/4D0/ebeamEnucl**2))))*
0286      &                XY*Q2(I))
0287             ELSE
0288               WTGAGA=WTGAGA*(0.5D0*(1D0+(1D0-XY)**2)-
0289      &        PMS(I)*XY**2/Q2(I))
0290             ENDIF
0291             IF(MINT(106+I).EQ.0) MINT(14+I)=22
0292           ENDIF
0293   140   CONTINUE
0294         VINT(319)=WTGAGA
0295         MINT(143)=LOOP
0296
0297 C...Update pTmin and cross section information.
0298         IF(MSTP(82).LE.1) THEN
0299           PTMN=PARP(81)*(VINT(1)/PARP(89))**PARP(90)
0300         ELSE
0301           PTMN=PARP(82)*(VINT(1)/PARP(89))**PARP(90)
0302         ENDIF
0303         VINT(149)=4D0*PTMN**2/VINT(2)
0304         VINT(154)=PTMN
0305         CALL PYXTOT
0306
0307 C...ECA...Routine modified for rad-corrections
0308 C...Generate photons inside leptons and
0309 C...calculate photon flux weight factors.
0310       ELSEIF(IGAGA.EQ.5) THEN
0311         write(*,*)"In pygaga with IGAGA.EQ.5"
0312         ISUB=MINT(1)
0313         MINT(15)=0
0314         MINT(16)=0
0315
0316 C...Generate phase space point and check against cuts.
0317         LOOP=0
0318   121   LOOP=LOOP+1
0319         DO 131 I=1,2
0320           IF(MINT(140+I).NE.0) THEN
0321 C...Pick x and Q2
0322             MINT(199)=0
0323             geny=mcSet_YMin*(mcSet_YMax/mcSet_YMin)**PYR(0)
0324             genQ2=mcSet_Q2Min*(mcSet_Q2Max/mcSet_Q2Min)**PYR(0)
0325             genx = genQ2/geny/(4.*ebeam*pbeam)
0326             genW2 = massp**2.+(genQ2*(1./genx-1.))
0327             gennu=geny*ebeamEnucl
0328             geneprim = ebeamEnucl - gennu
0329             write(*,*) geny, genq2, genx, gennu, genphi, ebeamEnucl
0330 C....Check to have sensible ranges for variables
0331             minq2 = PMS(1) * geny**2. / (1.- geny)
0332             if (genQ2.lt.minq2) then
0333                GOTO 121
0334             endif
0335 C... check x and Q2 go toghether
0336             if (genQ2.gt.(2.*gennu*massp)) then
0337                GOTO 121
0338             endif
0339 C            temp = ((genQ2-minq2)/(2D0*ebeamEnucl*geneprim))-1D0
0340             temp=(genQ2/(2.*ebeamEnucl*geneprim))-1.
0341             if ((temp.le.1.00).and.(temp.ge.-1.0)) then
0342                 etheta = acos(temp)
0343                 genthe = sngl(etheta)
0344                 genphi=PARU(2)*PYR(0)
0345                 PHI(I)=genphi
0346                 ppt=tan(etheta)
0347                 ppx=ppt*cos(PHI(I))
0348                 ppy=ppt*sin(PHI(I))
0349             else
0350                 GOTO 121
0351             endif
0352
0353             if ((genW2.lt.CKIN(77)**2).or.
0354      &          (CKIN(78).gt.0.and.genW2.gt.CKIN(78)**2)) then
0355                GOTO 121
0356             endif
0357             write(*,*) geny, genq2, genx, gennu, genphi, ebeamEnucl
0358
0359             ntries=0
0360   122       if (qedrad.eq.1) then
0361             write(*,*) geny, genq2, genx, gennu, genphi, ebeamEnucl
0362               call radgen_event
0363             endif
0364             if (qedrad.eq.0) then
0365              Y(I)=geny
0366              Q2(I)=genq2
0367             elseif ((mcRadCor_EBrems.eq.mcRadCor_EBrems).and.
0368      &          (mcRadCor_ThetaBrems.eq.mcRadCor_ThetaBrems)) then
0369              Y(I)=dble(mcRadCor_NuTrue)/dble(mcSet_EneBeam)
0370              Q2(I)=dble(mcRadCor_Q2True)
0371             write(*,*) geny, genq2, genx, gennu, genphi, ebeamEnucl
0372             write(*,*) mcRadCor_NuTrue, mcSet_EneBeam, mcRadCor_Q2True
0373             else
0374              write(*,*)"I go to 122 again"
0375              write(*,*) mcRadCor_ThetaBrems,mcRadCor_EBrems,
0376      &                  mcEvent_iEvent
0377              GOTO 122
0378             endif
0379 C            write(*,*) 'geny, genq2, genx, genW2, genNu',
0380 C     &                  geny, genq2, genx, genW2, genNu
0381             X(I)=((PMC(3)*Y(I))-Q2(I))/PMC(I)
0382 C P.L. ...An event with W^2_T<4will be generated new by RADGEN at the
0383 C      ...same kinematic point, the number of tries needed by RADGEN is
0384 C      ...counted and saved in the variable rcweight!
0385              IF (qedrad.ne.0) then
0386 C              IF(mcradcor_cType.eq.'qela') then
0387 C                 GOTO 122
0388 C              ENDIF
0389 C               IF(mcradcor_cType.eq.'elas') then
0390 C                 GOTO 122
0391 C              ENDIF
0392               IF(dble(mcRadCor_W2True).LT.
0393      &               (CKIN(77)**2-1.D-4*abs(CKIN(77)**2))) THEN
0394                 MINT(199)=MINT(199)+1
0395                 GOTO 122
0396               ENDIF
0397              ENDIF
0398              ntries=ntries+1
0399              IF(ntries.ge.20) GOTO 121
0400
0401 C ...... New try to implement weights directly into Pythia
0402           sigobs=0.0D0
0403           sigtrue=0.0D0
0404           rccorr=1.0D0
0405         if (qedrad.eq.1) then
0406           call MKF2(genq2,genx,
0407      +              mcSet_TarA,mcSet_TarZ,py6f2,py6f1)
0408           sigobs=pyth_xsec(geny,genx,genq2,py6f1, py6f2)
0409           IF(mcRadCor_EBrems.eq.0) THEN
0410            IF (sig1g.gt.0.D0) then
0411             rccorr=(tbor+tine)/sig1g/(DBLE(MINT(199))+1.0D0)
0412            ELSE
0413             rccorr=0.D0
0414            ENDIF
0415           ELSEIF(mcRadCor_EBrems.gt.0) THEN
0416           call MKF2(dble(mcRadCor_Q2True),dble(mcRadCor_XTrue),
0417      +           mcSet_TarA,mcSet_TarZ,py6f2,py6f1)
0418           sigtrue=pyth_xsec(dble( mcRadCor_YTrue), dble(mcRadCor_XTrue),
0419      +            dble(mcRadCor_Q2True), py6f1, py6f2)
0420            IF ((sig1g.gt.0.D0).and.(sigtrue.gt.0.D0)) then
0421             rccorr=(tbor+tine)/sig1g*sigobs/sigtrue/
0422      &             (DBLE(MINT(199))+1.0D0)
0423            ELSE
0424             rccorr=0.D0
0425            ENDIF
0426           ENDIF
0427         ENDIF
0428             IF(X(I).GT.(XMAX(I)+1.D-4*abs(XMAX(I)))) THEN
0429                GOTO 121
0430             ENDIF
0431 C...Cuts on internal consistency in x and Q2.
0432             IF(Q2(I).LT.X(I)**2*PMS(I)/(1D0-X(I))) then
0433                GOTO 121
0434             endif
0435             IF(Q2(I).GT.(1D0-X(I))*(VINT(302)-2D0*PMS(3-I))-
0436      &      (2D0-X(I)**2)*PMS(I)/(1D0-X(I))) THEN
0437               GOTO 121
0438             ENDIF
0439 C...Cuts on y and theta.
0440             IF(Y(I).LT.CKIN(71+2*I).OR.Y(I).GT.CKIN(72+2*I)) THEN
0441                GOTO 121
0442             ENDIF
0443             RAT=((1D0-X(I))*Q2(I)-X(I)**2*PMS(I))/
0444      &      ((1D0-X(I))**2*(VINT(302)-2D0*PMS(3-I)-2D0*PMS(I)))
0445             THETA(I)=2D0*ASIN(SQRT(MAX(0D0,MIN(1D0,RAT))))
0446             IF(THETA(I).LT.CKIN(67+2*I)) THEN
0447               GOTO 121
0448             ENDIF
0449             IF(CKIN(68+2*I).GT.0D0.AND.THETA(I).GT.CKIN(68+2*I))
0450      &      GOTO 121
0451
0452 C...Phi angle isotropic. Reconstruct pT.
0453 C            PT(I)=SQRT(((1D0-X(I))*PMC(I))**2/(4D0*VINT(302))-
0454 C     &      PMS(I))*SIN(THETA(I))
0455             temp=((1D0-X(I))*PMC(I))**2/(4D0*VINT(302))-PMS(I)
0456             PT(I)=(SQRT(temp))*SIN(THETA(I))
0457 C ... try 'new' phi
0458             IF ((qedrad.ne.0).and.(mcRadCor_EBrems.gt.0)) then
0459              emom=sqrt(geneprim**2-masse**2)
0460              PHI(I)=atan2((emom*ppy+dplabg(2)),(emom*ppx+dplabg(1)))
0461              IF (PHI(I).lt.0) THEN
0462                PHI(I)=PHI(I)+PARU(2)
0463               ENDIF
0464             ENDIF
0465 C...Store info on variables selected, for documentation purposes.
0466             VINT(2+I)=-SQRT(Q2(I))
0467             VINT(304+I)=X(I)
0468             VINT(306+I)=Q2(I)
0469             VINT(308+I)=Y(I)
0470             VINT(310+I)=THETA(I)
0471             VINT(312+I)=PHI(I)
0472           ELSE
0473             VINT(304+I)=1D0
0474             VINT(306+I)=0D0
0475             VINT(308+I)=1D0
0476             VINT(310+I)=0D0
0477             VINT(312+I)=0D0
0478           ENDIF
0479   131   CONTINUE
0480
0481 C...Cut on W combines info from two sides.
0482         IF(MINT(141).NE.0.AND.MINT(142).NE.0) THEN
0483           W2=-Q2(1)-Q2(2)+0.5D0*X(1)*PMC(1)*X(2)*PMC(2)/VINT(302)-
0484      &    2D0*PT(1)*PT(2)*COS(PHI(1)-PHI(2))+2D0*
0485      &    SQRT((0.5D0*X(1)*PMC(1)/VINT(301))**2+Q2(1)-PT(1)**2)*
0486      &    SQRT((0.5D0*X(2)*PMC(2)/VINT(301))**2+Q2(2)-PT(2)**2)
0487           IF(W2.LT.W2MIN) THEN
0488             GOTO 121
0489           ENDIF
0490           IF(CKIN(78).GT.0D0.AND.W2.GT.CKIN(78)**2) GOTO 121
0491           PMS1=-Q2(1)
0492           PMS2=-Q2(2)
0493         ELSEIF(MINT(141).NE.0) THEN
0494           W2=(VINT(302)+PMS(1))*X(1)+PMS(2)*(1D0-X(1))
0495           PMS1=-Q2(1)
0496           PMS2=PMS(2)
0497         ELSEIF(MINT(142).NE.0) THEN
0498           W2=(VINT(302)+PMS(2))*X(2)+PMS(1)*(1D0-X(2))
0499           PMS1=PMS(1)
0500           PMS2=-Q2(2)
0501         ENDIF
0502
0503 C...Store kinematics info for photon(s) in subsystem cm frame.
0504         VINT(2)=W2
0505         VINT(1)=SQRT(W2)
0506         VINT(291)=0D0
0507         VINT(292)=0D0
0508         VINT(293)=0.5D0*SQRT((W2-PMS1-PMS2)**2-4D0*PMS1*PMS2)/VINT(1)
0509         VINT(294)=0.5D0*(W2+PMS1-PMS2)/VINT(1)
0510         VINT(295)=SIGN(SQRT(ABS(PMS1)),PMS1)
0511         VINT(296)=0D0
0512         VINT(297)=0D0
0513         VINT(298)=-VINT(293)
0514         VINT(299)=0.5D0*(W2+PMS2-PMS1)/VINT(1)
0515         VINT(300)=SIGN(SQRT(ABS(PMS2)),PMS2)
0516
0517 C...Assign weight for photon flux; different for transverse and
0518 C...longitudinal photons. Flag incoming unresolved photon.
0519         WTGAGA=1D0
0520         DO 141 I=1,2
0521           IF(MINT(140+I).NE.0) THEN
0522             WTGAGA=WTGAGA*2D0*(PARU(101)/PARU(2))*
0523      &      (LOG(real(mcSet_YMax))-LOG(real(mcSet_YMin)))*
0524      &      (LOG(real(mcSet_Q2Max))-LOG(real(mcSet_Q2Min)))
0525           XY=Y(I)
0526             IF(ISUB.EQ.132.OR.ISUB.EQ.134.OR.ISUB.EQ.136) THEN
0527               IF((MINT(11).EQ.22).and.
0528      &             (MINT(12).EQ.2212.or.MINT(12).EQ.2112)) THEN
0529                XXY=XY*ebeamEnucl/ebeamEnucl
0530                WTGAGA=WTGAGA*(1D0/(1D0+(Q2(I)/XXY**2/ebeamEnucl**2))*
0531      &               (1D0-XXY-(Q2(I)/4D0/ebeamEnucl**2)))/
0532      &               Q2(I)/XXY**2/ebeamEnucl*
0533      &               (ebeamEnucl*XXY-Q2(I)/2D0/massp)*XXY*Q2(I)
0534               ELSE
0535                 WTGAGA=WTGAGA*(1D0-XY)
0536               ENDIF
0537             ELSEIF(I.EQ.1.AND.(ISUB.EQ.139.OR.ISUB.EQ.140)) THEN
0538               WTGAGA=WTGAGA*(1D0-XY)
0539             ELSEIF(I.EQ.2.AND.(ISUB.EQ.138.OR.ISUB.EQ.140)) THEN
0540               WTGAGA=WTGAGA*(1D0-XY)
0541             ELSEIF((MINT(11).EQ.22).and.
0542      &             (MINT(12).EQ.2212.or.MINT(12).EQ.2112)) THEN
0543               XXY=XY*ebeamEnucl/ebeamEnucl
0544               WTGAGA=WTGAGA*(0.5D0*((ebeamEnucl*XXY-Q2(I)/2D0/
0545      &               massp)/Q2(I)/XXY**2/ebeamEnucl*
0546      &               (XXY**2*(1D0-(2D0*masse**2/Q2(I)))+
0547      &               (2D0/(1D0+(Q2(I)/XXY**2/ebeamEnucl**2)))*
0548      &               (1D0-XXY-(Q2(I)/4D0/ebeamEnucl**2))))*XXY*Q2(I))
0549             ELSE
0550               WTGAGA=WTGAGA*(0.5D0*(1D0+(1D0-XY)**2)-
0551      &        PMS(I)*XY**2/Q2(I))
0552             ENDIF
0553             IF(MINT(106+I).EQ.0) MINT(14+I)=22
0554           ENDIF
0555  141   CONTINUE
0556         WTGAGA=WTGAGA*rccorr
0557         VINT(319)=WTGAGA
0558         MINT(143)=LOOP
0559 C...Update pTmin and cross section information.
0560         IF(MSTP(82).LE.1) THEN
0561           PTMN=PARP(81)*(VINT(1)/PARP(89))**PARP(90)
0562         ELSE
0563           PTMN=PARP(82)*(VINT(1)/PARP(89))**PARP(90)
0564         ENDIF
0565         VINT(149)=4D0*PTMN**2/VINT(2)
0566         VINT(154)=PTMN
0567         CALL PYXTOT
0568
0569 C...Reconstruct kinematics of photons inside leptons.
0570       ELSEIF(IGAGA.EQ.4) THEN
0571
0572 C...Make place for incoming particles and scattered leptons.
0573         MOVE=3
0574         IF(MINT(141).NE.0.AND.MINT(142).NE.0) MOVE=4
0575         MINT(4)=MINT(4)+MOVE
0576         DO 160 I=MINT(84)-MOVE,MINT(83)+1,-1
0577           IF(K(I,1).EQ.21) THEN
0578             DO 150 J=1,5
0579               K(I+MOVE,J)=K(I,J)
0580               P(I+MOVE,J)=P(I,J)
0581               V(I+MOVE,J)=V(I,J)
0582   150       CONTINUE
0583             IF(K(I,3).GT.MINT(83).AND.K(I,3).LE.MINT(84))
0584      &      K(I+MOVE,3)=K(I,3)+MOVE
0585             IF(K(I,4).GT.MINT(83).AND.K(I,4).LE.MINT(84))
0586      &      K(I+MOVE,4)=K(I,4)+MOVE
0587             IF(K(I,5).GT.MINT(83).AND.K(I,5).LE.MINT(84))
0588      &      K(I+MOVE,5)=K(I,5)+MOVE
0589           ENDIF
0590   160   CONTINUE
0591         DO 170 I=MINT(84)+1,N
0592           IF(K(I,3).GT.MINT(83).AND.K(I,3).LE.MINT(84))
0593      &    K(I,3)=K(I,3)+MOVE
0594   170   CONTINUE
0595
0596 C...Fill in incoming particles.
0597         DO 190 I=MINT(83)+1,MINT(83)+MOVE
0598           DO 180 J=1,5
0599             K(I,J)=0
0600             P(I,J)=0D0
0601             V(I,J)=0D0
0602   180     CONTINUE
0603   190   CONTINUE
0604         DO 200 I=1,2
0605           K(MINT(83)+I,1)=21
0606           IF(MINT(140+I).NE.0) THEN
0607             K(MINT(83)+I,2)=MINT(140+I)
0608             P(MINT(83)+I,5)=VINT(302+I)
0609           ELSE
0610             K(MINT(83)+I,2)=MINT(10+I)
0611             P(MINT(83)+I,5)=VINT(2+I)
0612           ENDIF
0613           P(MINT(83)+I,3)=0.5D0*SQRT((PMC(3)**2-4D0*PMS(1)*PMS(2))/
0614      &    VINT(302))*(-1D0)**(I+1)
0615           P(MINT(83)+I,4)=0.5D0*PMC(I)/VINT(301)
0616   200   CONTINUE
0617
0618 C...New mother-daughter relations in documentation section.
0619         IF(MINT(141).NE.0.AND.MINT(142).NE.0) THEN
0620           K(MINT(83)+1,4)=MINT(83)+3
0621           K(MINT(83)+1,5)=MINT(83)+5
0622           K(MINT(83)+2,4)=MINT(83)+4
0623           K(MINT(83)+2,5)=MINT(83)+6
0624           K(MINT(83)+3,3)=MINT(83)+1
0625           K(MINT(83)+5,3)=MINT(83)+1
0626           K(MINT(83)+4,3)=MINT(83)+2
0627           K(MINT(83)+6,3)=MINT(83)+2
0628         ELSEIF(MINT(141).NE.0) THEN
0629           K(MINT(83)+1,4)=MINT(83)+3
0630           K(MINT(83)+1,5)=MINT(83)+4
0631           K(MINT(83)+2,4)=MINT(83)+5
0632           K(MINT(83)+3,3)=MINT(83)+1
0633           K(MINT(83)+4,3)=MINT(83)+1
0634           K(MINT(83)+5,3)=MINT(83)+2
0635         ELSEIF(MINT(142).NE.0) THEN
0636           K(MINT(83)+1,4)=MINT(83)+4
0637           K(MINT(83)+2,4)=MINT(83)+3
0638           K(MINT(83)+2,5)=MINT(83)+5
0639           K(MINT(83)+3,3)=MINT(83)+2
0640           K(MINT(83)+4,3)=MINT(83)+1
0641           K(MINT(83)+5,3)=MINT(83)+2
0642         ENDIF
0643
0644 C...Fill scattered lepton(s).
0645         DO 210 I=1,2
0646           IF(MINT(140+I).NE.0) THEN
0647             LSC=MINT(83)+MIN(I+2,MOVE)
0648             K(LSC,1)=21
0649             K(LSC,2)=MINT(140+I)
0650             P(LSC,1)=PT(I)*COS(PHI(I))
0651             P(LSC,2)=PT(I)*SIN(PHI(I))
0652             P(LSC,4)=(1D0-X(I))*P(MINT(83)+I,4)
0653             P(LSC,3)=SQRT(P(LSC,4)**2-PMS(I))*COS(THETA(I))*
0654      &      (-1D0)**(I-1)
0655             P(LSC,5)=VINT(302+I)
0656           ENDIF
0657   210   CONTINUE
0658
0659 C...Find incoming four-vectors to subprocess.
0660         K(N+1,1)=21
0661         IF(MINT(141).NE.0) THEN
0662           DO 220 J=1,4
0663             P(N+1,J)=P(MINT(83)+1,J)-P(MINT(83)+3,J)
0664   220     CONTINUE
0665         ELSE
0666           DO 230 J=1,4
0667             P(N+1,J)=P(MINT(83)+1,J)
0668   230     CONTINUE
0669         ENDIF
0670         K(N+2,1)=21
0671         IF(MINT(142).NE.0) THEN
0672           DO 240 J=1,4
0673             P(N+2,J)=P(MINT(83)+2,J)-P(MINT(83)+MOVE,J)
0674   240     CONTINUE
0675         ELSE
0676           DO 250 J=1,4
0677             P(N+2,J)=P(MINT(83)+2,J)
0678   250     CONTINUE
0679         ENDIF
0680
0681 C...Define boost and rotation between hadronic subsystem and
0682 C...collision rest frame; boost hadronic subsystem to this frame.
0683         DO 260 J=1,3
0684           BETA(J)=(P(N+1,J)+P(N+2,J))/(P(N+1,4)+P(N+2,4))
0685   260   CONTINUE
0686         CALL PYROBO(N+1,N+2,0D0,0D0,-BETA(1),-BETA(2),-BETA(3))
0687         BPHI=PYANGL(P(N+1,1),P(N+1,2))
0688         CALL PYROBO(N+1,N+2,0D0,-BPHI,0D0,0D0,0D0)
0689         BTHETA=PYANGL(P(N+1,3),P(N+1,1))
0690         CALL PYROBO(MINT(83)+MOVE+1,N,BTHETA,BPHI,BETA(1),BETA(2),
0691      &  BETA(3))
0692
0693 C...Add on scattered leptons to final state.
0694         DO 280 I=1,2
0695           IF(MINT(140+I).NE.0) THEN
0696             LSC=MINT(83)+MIN(I+2,MOVE)
0697             N=N+1
0698             DO 270 J=1,5
0699               K(N,J)=K(LSC,J)
0700               P(N,J)=P(LSC,J)
0701               V(N,J)=V(LSC,J)
0702   270       CONTINUE
0703             K(N,1)=1
0704             K(N,3)=LSC
0705           ENDIF
0706   280   CONTINUE
0707       ENDIF
0708
0709   290 CONTINUE
0710       RETURN
0711       END