Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:07:29

0001 c  reshuffled from sem, sto, sha
0002 
0003 c    contains DIS, and unused 3P stuff
0004 c             ---      ---------------
0005 
0006 
0007 
0008 
0009 c###########################################################################
0010 c###########################################################################
0011 c###########################################################################
0012 c###########################################################################
0013 c
0014 c                              DIS
0015 c
0016 c###########################################################################
0017 c###########################################################################
0018 c###########################################################################
0019 c###########################################################################
0020 
0021 
0022 
0023 
0024 CDECK  ID>, PHO_GPHERA from phojet by Ralph Engel
0025       SUBROUTINE phoGPHERAepo(imod)
0026 C**********************************************************************
0027 C
0028 C     interface to call PHOJET (variable energy run) with
0029 C     HERA kinematics, photon as particle 2
0030 C
0031 C     equivalent photon approximation to get photon flux
0032 C
0033 C     input:     imod=0        Main initialization
0034 C                     1        Event initialization
0035 C             from /photrans/ and /lept1/
0036 C                EE1=ebeam     proton energy (LAB system)
0037 C                EE2=elepti    electron energy (LAB system)
0038 C             from /psar12/:
0039 C                YMIN2=ydmin    lower limit of Y
0040 C                        (energy fraction taken by photon from electron)
0041 C                YMAX2=ydmax    upper limit of Y
0042 C                Q2MN2=qdmin   lower limit of photon virtuality
0043 C                Q2MX2=qdmax   upper limit of photon virtuality
0044 C
0045 C**********************************************************************
0046       include 'epos.inc'
0047       include 'epos.incsem'
0048       common/photrans/phoele(4),ebeam
0049       double precision pgampr,rgampr
0050       common/cgampr/pgampr(5),rgampr(4)
0051 
0052 c PHOJET common
0053       PARAMETER ( phoPI   = 3.14159265359D0 )
0054 
0055       DOUBLE PRECISION EE1,EE2,PROM2,YMIN2,YMAX2,THMIN2,THMAX2
0056      &                ,ELEM,ELEM2,Q2MN,Q2MX,XIMAX,XIMIN,XIDEL,DELLY
0057      &                ,FLUXT,FLUXL,Q2LOW,Y,FFT,FFL,AY,AY2,YY,WGMAX
0058      &                ,ECMIN2,ECMAX2,Q22MIN,Q22AVE,Q22AV2,Q22MAX,AN2MIN
0059      &                ,AN2MAX,YY2MIN,YY2MAX,EEMIN2,gamNsig0
0060      &                ,Q21MIN,Q21MAX,AN1MIN,AN1MAX,YY1MIN,YY1MAX
0061       SAVE EE1,EE2,PROM2,YMIN2,YMAX2,THMIN2,THMAX2
0062      &                ,ELEM,ELEM2,Q2MN,Q2MX,XIMAX,XIMIN,XIDEL,DELLY
0063      &                ,FLUXT,FLUXL,Q2LOW,Y,FFT,FFL,AY,AY2,YY,WGMAX
0064      &                ,ECMIN2,ECMAX2,Q22MIN,Q22AVE,Q22AV2,Q22MAX,AN2MIN
0065      &                ,AN2MAX,YY2MIN,YY2MAX,EEMIN2,gamNsig0
0066      &                ,Q21MIN,Q21MAX,AN1MIN,AN1MAX,YY1MIN,YY1MAX
0067       INTEGER ITRY,ITRW
0068       SAVE ITRY,ITRW
0069 
0070 C LOCAL VARIABLES
0071       DOUBLE PRECISION PINI(5),PFIN(5),GGECM,PFTHE,YQ2,YEFF,Q2LOG,WGH,Q2
0072      &                ,WEIGHT,Q2E,E1Y,PHI,COF,SIF,WGY,DAY,P1(5),P2(4)
0073      &                ,drangen,ALLM97,gamNsig
0074 
0075 
0076 
0077       if(imod.eq.0)then        !initialization
0078 
0079         EE1=dble(ebeam)         !proton energy (LAB system)
0080         EE2=dble(elepti)        !electron energy (LAB system)
0081 
0082 
0083         if(ish.ge.2)WRITE(ifch,*) 'phoGPHERAepo: energy to process'
0084      *                          ,EE1,EE2
0085 C  assign particle momenta according to HERA kinematics
0086 C  proton data
0087 
0088         if(idtarg.ne.0)then
0089           call idmass(idtarg,ams)
0090         else
0091           call idmass(1120,ams)
0092         endif
0093         PROM2 = dble(ams)**2
0094 C electron data
0095         call idmass(12,ams)
0096         ELEM = dble(ams)
0097         ELEM2 = ELEM**2
0098 C
0099         Q2MN = dble(qdmin)
0100         Q2MX = dble(qdmax)
0101 C
0102         YMIN2=dble(ydmin)
0103         YMAX2=dble(ydmax)
0104         XIMAX = LOG(YMAX2)
0105         XIMIN = LOG(YMIN2)
0106         XIDEL = XIMAX-XIMIN
0107 C
0108         THMIN2=dble(themin*pi/180.)
0109         THMAX2=dble(themax*pi/180.)
0110 C
0111         IF(Q2MN.GT.ELEM2*YMIN2**2/(1.D0-YMIN2))
0112      &  WRITE(*,'(/1X,A,1P2E11.4)')
0113      &  'phoGPHERAepo: lower Q2 cutoff larger than kin. limit:',
0114      &  Q2MN,ELEM2*YMIN2**2/(1.D0-YMIN2)
0115 C
0116         IF(ish.GE.6)THEN
0117           Max_tab = 50
0118           DELLY = LOG(YMAX2/YMIN2)/DBLE(Max_tab-1)
0119           FLUXT = 0.D0
0120           FLUXL = 0.D0
0121 
0122           WRITE(ifch,'(1X,A,I5)')
0123      &  'phoGPHERAepo: table of photon flux (trans/long)',Max_tab
0124           DO 100 I=1,Max_tab
0125             Y = EXP(XIMIN+DELLY*DBLE(I-1))
0126             Q2LOW = MAX(Q2MN,ELEM2*Y**2/(1.D0-Y))
0127             FFT = ((1.D0+(1.D0-Y)**2)/Y*LOG(Q2MX/Q2LOW)
0128      &        -2.D0*ELEM2*Y*(1.D0/Q2LOW-1.D0/Q2MX))/(2.D0*phoPI*137.D0)
0129             FFL = 2.D0*(1.D0-Y)/Y*LOG(Q2MX/Q2LOW)/(2.D0*phoPI*137.D0)
0130             FLUXT = FLUXT + Y*FFT
0131             FLUXL = FLUXL + Y*FFL
0132             WRITE(ifch,'(5X,1P3E14.4)') Y,FFT,FFL
0133  100      CONTINUE
0134           FLUXT = FLUXT*DELLY
0135           FLUXL = FLUXL*DELLY
0136           WRITE(ifch,'(1X,A,1P2E12.4)')
0137      &  'PHOGPHERA: integrated flux (trans./long.):',FLUXT,FLUXL
0138         ENDIF
0139 C
0140 
0141         YY = YMIN2
0142         Q2LOW = MAX(Q2MN,ELEM2*YY**2/(1.D0-YY))
0143         WGMAX = (1.D0+(1.D0-YY)**2)*LOG(Q2MX/Q2LOW)
0144      &       -2.D0*ELEM2*YY*(1.D0/Q2LOW-1.D0/Q2MX)*YY
0145         WGMAX = WGMAX+2.D0*(1.D0-YY)*LOG(Q2MX/Q2LOW)
0146 
0147         ECMIN2 = dble(egymin)**2
0148         ECMAX2 = dble(egymax)**2
0149         EEMIN2 = dble(elomin)
0150         AY = 0.D0
0151         AY2 = 0.D0
0152         Q22MIN = 1.D30
0153         Q22AVE = 0.D0
0154         Q22AV2 = 0.D0
0155         Q22MAX = 0.D0
0156         AN2MIN = 1.D30
0157         AN2MAX = 0.D0
0158         YY2MIN = 1.D30
0159         YY2MAX = 0.D0
0160         ITRY = 0
0161         ITRW = 0
0162         gamNsig0 = 5d0 * ALLM97(Q2LOW,WGMAX)
0163 
0164       elseif(imod.eq.1)then     !event
0165 
0166 C
0167 C  sample y
0168         ITRY = ITRY+1
0169         ntry=0
0170  175    CONTINUE
0171         ntry=ntry+1
0172         IF(ntry.ge.1000) THEN
0173             WRITE(*,'(1X,A,2E12.5,2(1X,1A,1X,3E13.5))')
0174      &        'phoGPHERAepo: problem with cuts:',PFIN(4),EEMIN2,'|'
0175      &       ,THMIN2,PFTHE,THMAX2,'|',ECMIN2,GGECM,ECMAX2
0176             call utstop("Problem with cuts in phoGPHERAepo !&")
0177         ENDIF
0178            ITRW = ITRW+1
0179           YY = EXP(XIDEL*drangen(AY)+XIMIN)
0180           YEFF = 1.D0+(1.D0-YY)**2+2.D0*(1.D0-YY)
0181           Q2LOW = MAX(Q2MN,ELEM2*YY**2/(1.D0-YY))
0182           Q2LOG = LOG(Q2MX/Q2LOW)
0183           WGH = YEFF*Q2LOG-2.D0*ELEM2*YY**2*(1.D0/Q2LOW-1.D0/Q2MX)
0184           IF(WGMAX.LE.WGH) THEN
0185             WRITE(*,'(1X,A,3E12.5)')
0186      &        'phoGPHERAepo: inconsistent weight:',YY,WGMAX,WGH
0187             call utstop("Problem with YY in phoGPHERAepo !$")
0188           ENDIF
0189         IF(drangen(AY2)*WGMAX.GT.WGH) GOTO 175
0190 C  sample Q2
0191  185    CONTINUE
0192           Q2 = Q2LOW*EXP(Q2LOG*drangen(YY))
0193           WEIGHT = (YEFF-2.D0*ELEM2*YY**2/Q2)/YEFF
0194           IF(WEIGHT.GE.1d0) THEN
0195             WRITE(*,'(1X,A,3E12.5)')
0196      &        'phoGPHERAepo: inconsistent weight:',YY,Q2,YEFF,WEIGHT
0197             call utstop("Problem with Q2 in phoGPHERAepo !$")
0198           ENDIF
0199         IF(WEIGHT.LT.drangen(Q2)) GOTO 185
0200 C
0201         if(ish.ge.2)WRITE(ifch,*) 'phoGPHERAepo: event with Q2,Y:',Q2,YY
0202 
0203 C  incoming electron
0204         PINI(1) = 0.D0
0205         PINI(2) = 0.D0
0206         PINI(3) = sqrt((EE2+ELEM)*(EE2-ELEM))
0207         PINI(4) = EE2
0208         PINI(5) = ELEM
0209 C  outgoing electron
0210         YQ2 = SQRT((1.D0-YY)*Q2)
0211         Q2E = Q2/(4.D0*EE2)
0212         E1Y = EE2*(1.D0-YY)
0213         phi=2d0*phoPI*drangen(E1Y)
0214         COF=cos(phi)
0215         SIF=sin(phi)
0216         PFIN(1) = YQ2*COF
0217         PFIN(2) = YQ2*SIF
0218         PFIN(4) = E1Y+Q2E
0219         PFIN(5) = ELEM
0220 c        PFIN(3) = E1Y+Q2E
0221         PFIN(3)=(PFIN(4)+sqrt(YQ2*YQ2+ELEM2))
0222      *         *(PFIN(4)-sqrt(YQ2*YQ2+ELEM2))
0223         if(PFIN(3).ge.0d0)then
0224           PFIN(3) = sqrt(PFIN(3))
0225         else
0226           PFIN(3) = E1Y+Q2E
0227         endif
0228         GQ2 = sngl(Q2)
0229         GWD = 4.*ebeam*elepti*sngl(YY)
0230 C  polar angle
0231         PFTHE = ACOS(PFIN(3)/PFIN(4))
0232 C  electron tagger
0233         IF(PFIN(4).GT.EEMIN2) THEN
0234           IF((PFTHE.LT.THMIN2).OR.(PFTHE.GT.THMAX2)) GOTO 175
0235         ENDIF
0236 C  photon momentum
0237         P2(1) = -PFIN(1)
0238         P2(2) = -PFIN(2)
0239         P2(3) = PINI(3)-PFIN(3)
0240         P2(4) = PINI(4)-PFIN(4)
0241 C  proton momentum
0242         P1(1) = 0.D0
0243         P1(2) = 0.D0
0244         P1(3) = -SQRT(EE1**2-PROM2)
0245         P1(4) = EE1
0246         P1(5) = sqrt(prom2)
0247 C  ECMS cut
0248         GGECM = (P1(4)+P2(4))**2-(P1(1)+P2(1))**2
0249      &         -(P1(2)+P2(2))**2-(P1(3)+P2(3))**2
0250         IF((GGECM.LT.ECMIN2).OR.(GGECM.GT.ECMAX2)) GOTO 175
0251         GGECM = SQRT(GGECM)
0252 C accept A2 and W according to gamma-p cross section (function of F2)
0253         gamNsig=ALLM97(Q2,GGECM)/gamNsig0
0254         if(gamNsig.ge.1d0)print *,'R>1 in DIS',gamNsig
0255         if(drangen(gamNsig).gt.gamNsig)goto 175 !no interaction
0256 C output
0257         engy=sngl(GGECM)
0258         xbjevt=GQ2/GWD
0259         qsqevt=GQ2
0260 c gamma
0261         rgampr(1) = P2(1)
0262         rgampr(2) = P2(2)
0263         rgampr(3) = P2(3)
0264         rgampr(4) = P2(4)
0265 c boost gamma in proton rest frame to get the rotation vector
0266         call utlob2(1,P1(1),P1(2),P1(3),P1(4),P1(5)
0267      *               ,rgampr(1),rgampr(2),rgampr(3),rgampr(4),99)
0268 c array to define boost needed to put proton from lab to rest frame
0269         pgampr(1) = P1(1)
0270         pgampr(2) = P1(2)
0271         pgampr(3) = P1(3)
0272         pgampr(4) = P1(4)
0273         pgampr(5) = P1(5)
0274 c electron
0275         elepto=sngl(PFIN(4))
0276         phoele(1) = sngl(PFIN(1))
0277         phoele(2) = sngl(PFIN(2))
0278         phoele(3) = sngl(PFIN(3))
0279         phoele(4) = sngl(PFIN(4))
0280 
0281         if(ish.ge.2)then        !statistic
0282 C  statistics
0283           AY = AY+YY
0284           AY2 = AY2+YY*YY
0285           YY1MIN = YY2MIN
0286           YY1MAX = YY2MAX
0287           YY2MIN = MIN(YY2MIN,YY)
0288           YY2MAX = MAX(YY2MAX,YY)
0289           Q21MIN = Q22MIN
0290           Q21MAX = Q22MAX
0291           Q22MIN = MIN(Q22MIN,Q2)
0292           Q22MAX = MAX(Q22MAX,Q2)
0293           Q22AVE = Q22AVE+Q2
0294           Q22AV2 = Q22AV2+Q2*Q2
0295           AN1MIN = AN2MIN
0296           AN1MAX = AN2MAX
0297           AN2MIN = MIN(AN2MIN,PFTHE)
0298           AN2MAX = MAX(AN2MAX,PFTHE)
0299         endif
0300 C
0301       elseif(ish.ge.2)then     !statistic
0302 
0303         NITER=nrevt
0304 
0305         WGY = WGMAX*DBLE(ITRY)/DBLE(ITRW)/(137.D0*2.D0*phoPI)
0306         WGY = WGY*LOG(YMAX2/YMIN2)
0307         AY  = AY/DBLE(NITER)
0308         AY2 = AY2/DBLE(NITER)
0309         Q22AVE = Q22AVE/DBLE(NITER)
0310         Q22AV2 = Q22AV2/DBLE(NITER)
0311         if(NITER.gt.1)then
0312           DAY = SQRT((AY2-AY**2)/DBLE(NITER))
0313           Q22AV2 = SQRT((Q22AV2-Q22AVE**2)/DBLE(NITER))
0314         else
0315           DAY = 0d0
0316           Q22AV2 = 0d0
0317         endif
0318         SIGMAX = 1d0
0319         WEIGHT = WGY*SIGMAX*DBLE(NITER)/DBLE(ITRY)
0320 C  output of histograms
0321         WRITE(ifch,'(//1X,A,/1X,A,1PE12.3,A,/1X,A)')
0322      &'=========================================================',
0323      &' *****   simulated cross section: ',WEIGHT,' mb  *****',
0324      &'========================================================='
0325         WRITE(ifch,'(//1X,A,3I10)')
0326      &  'PHOGPHERA:SUMMARY:NITER,ITRY,ITRW',NITER,ITRY,ITRW
0327         WRITE(ifch,'(1X,A,1P2E12.4)') 'EFFECTIVE WEIGHT (FLUX,TOTAL)',
0328      &  WGY,WEIGHT
0329         WRITE(ifch,'(1X,A,1P2E12.4)') 'AVERAGE Y,DY                 ',AY
0330      &                                                              ,DAY
0331         WRITE(ifch,'(1X,A,1P2E12.4)') 'SAMPLED Y RANGE PHOTON       ',
0332      &  YY2MIN,YY2MAX
0333         WRITE(ifch,'(1X,A,1P2E12.4)') 'AVERAGE Q2,DQ2               ',
0334      &  Q22AVE,Q22AV2
0335         WRITE(ifch,'(1X,A,1P2E12.4)') 'SAMPLED Q2 RANGE PHOTON      ',
0336      &  Q22MIN,Q22MAX
0337         WRITE(ifch,'(1X,A,1P4E12.4)') 'SAMPLED THETA RANGE ELECTRON ',
0338      &  AN2MIN,AN2MAX,phoPI-AN2MAX,phoPI-AN2MIN
0339 C
0340 
0341       endif
0342 
0343 
0344       END
0345 
0346 
0347 CDECK  ID>, from PHO_ALLM97 in Phojet
0348       DOUBLE PRECISION FUNCTION ALLM97(Q2,W)
0349 C**********************************************************************
0350 C
0351 C     ALLM97 parametrization for gamma*-p cross section
0352 C     (for F2 see comments, code adapted from V. Shekelyan, H1)
0353 C
0354 C**********************************************************************
0355 
0356       IMPLICIT NONE
0357 
0358       SAVE
0359 
0360       DOUBLE PRECISION Q2,W
0361       DOUBLE PRECISION M02,M12,LAM2,M22
0362       DOUBLE PRECISION S11,S12,S13,A11,A12,A13,B11,B12,B13
0363       DOUBLE PRECISION S21,S22,S23,A21,A22,A23,B21,B22,B23
0364       DOUBLE PRECISION ALFA,XMP2,W2,Q02,S,T,T0,Z,CIN,
0365      &                 AP,BP,AR,BR,XP,XR,SR,SP,F2P,F2R
0366       DATA ALFA,XMP2 /112.2D0 , .8802D0 /
0367 
0368       W2=W*W
0369       ALLM97 = 0.D0
0370 
0371 C  pomeron
0372       S11   =   0.28067D0
0373       S12   =   0.22291D0
0374       S13   =   2.1979D0
0375       A11   =  -0.0808D0
0376       A12   =  -0.44812D0
0377       A13   =   1.1709D0
0378       B11   =   0.60243D0
0379       B12   =   1.3754D0
0380       B13   =   1.8439D0
0381       M12   =  49.457D0
0382 
0383 C  reggeon
0384       S21   =   0.80107D0
0385       S22   =   0.97307D0
0386       S23   =   3.4942D0
0387       A21   =   0.58400D0
0388       A22   =   0.37888D0
0389       A23   =   2.6063D0
0390       B21   =   0.10711D0
0391       B22   =   1.9386D0
0392       B23   =   0.49338D0
0393       M22   =   0.15052D0
0394 C
0395       M02   =   0.31985D0
0396       LAM2  =   0.065270D0
0397       Q02   =   0.46017D0 +LAM2
0398 
0399 C
0400       S=0.
0401       T=LOG((Q2+Q02)/LAM2)
0402       T0=LOG(Q02/LAM2)
0403       IF(Q2.GT.0.D0) S=LOG(T/T0)
0404       Z=1.D0
0405 
0406       IF(Q2.GT.0.D0) Z=(W2-XMP2)/(Q2+W2-XMP2)
0407 
0408       IF(S.LT.0.01D0) THEN
0409 
0410 C   pomeron part
0411 
0412         XP=1.D0 /(1.D0 +(W2-XMP2)/(Q2+M12))
0413 
0414         AP=A11
0415         BP=B11**2
0416 
0417         SP=S11
0418         F2P=SP*XP**AP*Z**BP
0419 
0420 C   reggeon part
0421 
0422         XR=1.D0 /(1.D0 +(W2-XMP2)/(Q2+M22))
0423 
0424         AR=A21
0425         BR=B21**2
0426 
0427         SR=S21
0428         F2R=SR*XR**AR*Z**BR
0429 
0430       ELSE
0431 
0432 C   pomeron part
0433 
0434         XP=1.D0 /(1.D0 +(W2-XMP2)/(Q2+M12))
0435 
0436         AP=A11+(A11-A12)*(1.D0 /(1.D0 +S**A13)-1.D0 )
0437 
0438         BP=B11**2+B12**2*S**B13
0439 
0440         SP=S11+(S11-S12)*(1.D0 /(1.D0 +S**S13)-1.D0 )
0441 
0442         F2P=SP*XP**AP*Z**BP
0443 
0444 C   reggeon part
0445 
0446         XR=1.D0 /(1.D0 +(W2-XMP2)/(Q2+M22))
0447 
0448         AR=A21+A22*S**A23
0449         BR=B21**2+B22**2*S**B23
0450 
0451         SR=S21+S22*S**S23
0452         F2R=SR*XR**AR*Z**BR
0453 
0454       ENDIF
0455 
0456 *     F2 = (F2P+F2R)*Q2/(Q2+M02)
0457 
0458       CIN=ALFA/(Q2+M02)*(1.D0 +4.D0*XMP2*Q2/(Q2+W2-XMP2)**2)/Z
0459       ALLM97 = CIN*(F2P+F2R)
0460 
0461       END
0462 
0463 
0464 
0465 
0466 
0467 
0468 c-----------------------------------------------------------------------
0469       subroutine lepexp(rxbj,rqsq)
0470 c-----------------------------------------------------------------------
0471 c     generates x_bjorken and q**2 according to an experimental
0472 c     distribution ( given in array xq(nxbj,nqsq) ).
0473 c-----------------------------------------------------------------------
0474       parameter (nxbj=10,nqsq=10)
0475       parameter (xbjmin=0.,qsqmin=4.)
0476       parameter (xbjwid=0.025, qsqwid=4.)
0477       dimension xq(nxbj,nqsq),vxq(nxbj*nqsq)
0478       equivalence (xq(1,1),vxq(1))
0479 
0480       data (vxq(i),i=1,50)/
0481      &         1304.02,   366.40,    19.84,    10.79,     6.42,
0482      &            4.54,     4.15,     3.38,     2.03,     1.56,
0483      &          241.63,  1637.26,   427.36,   164.51,    73.72,
0484      &           43.07,    20.73,    12.78,     9.34,     5.83,
0485      &            0.01,   724.66,   563.79,   275.08,   176.13,
0486      &          106.44,    85.82,    54.52,    37.12,    28.65,
0487      &            0.01,   202.40,   491.10,   245.13,   157.07,
0488      &          104.43,    61.05,    49.42,    37.84,    26.79,
0489      &            0.01,     3.77,   316.38,   226.92,   133.45,
0490      &           90.30,    63.67,    48.42,    35.73,    28.04/
0491       data (vxq(i),i=51,100)/
0492      &            0.01,     0.01,   153.74,   213.09,   114.14,
0493      &           76.26,    60.02,    43.15,    43.47,    25.60,
0494      &            0.01,     0.01,    39.31,   185.74,   108.56,
0495      &           88.40,    47.29,    39.35,    31.80,    22.91,
0496      &            0.01,     0.01,     0.01,   104.61,   107.01,
0497      &           66.24,    45.34,    37.45,    33.44,    23.78,
0498      &            0.01,     0.01,     0.01,    56.58,    99.39,
0499      &           67.78,    43.28,    35.98,    34.63,    18.31,
0500      &            0.01,     0.01,     0.01,    13.56,    76.25,
0501      &           64.30,    42.80,    28.56,    21.19,    20.75 /
0502 
0503       data init/0/
0504       init=init+1
0505       if(init.eq.1) then
0506       n=nxbj*nqsq
0507       sum=0.
0508       do 1 i=1,n
0509       sum=sum+vxq(i)
0510 1     continue
0511       do 2 i=2,n
0512 2     vxq(i)=vxq(i)+vxq(i-1)
0513       do 3 i=1,n
0514 3     vxq(i)=vxq(i)/sum
0515       endif
0516 
0517       n=nxbj*nqsq
0518       r=rangen()
0519       call utloc(vxq,n,r,iloc)
0520       if(iloc.ge.n) iloc=iloc-1
0521       i=mod(iloc,nxbj)+1
0522       if(i.eq.0) i=nxbj
0523       j=iloc/nxbj + 1
0524       dxint=vxq(1)
0525       if(iloc.gt.0) dxint=vxq(iloc+1)-vxq(iloc)
0526       dxbj=xbjwid*abs(r-vxq(iloc+1))/dxint
0527       dy=qsqwid*rangen()
0528       rxbj=xbjmin+xbjwid*float(i-1)+dxbj
0529       rqsq=qsqmin+qsqwid*float(j-1)+dy
0530       return
0531       end
0532 
0533 c-----------------------------------------------------------------------
0534       subroutine fremny(wp1,wm1,pnx,pny,sm,ic1,ic2,ic3,ic4,coord,ey0)
0535 c-----------------------------------------------------------------------
0536 c  treats remnant from deep inelastic process;
0537 c-----------------------------------------------------------------------
0538       include 'epos.inc'
0539       include 'epos.incsem'
0540       dimension coord(6),ic(2),ep(4),ey(3),ey0(3),ep3(4)
0541       double precision  ept(4),ept1(4)
0542 
0543       call utpri('fremny',ish,ishini,5)
0544       if(ish.ge.5)write(ifch,*)'writing remnant'
0545       if(ish.ge.5)write(ifch,*)
0546      *'wp1,wm1,pnx,pny,sm,ic1,ic2,ic3,ic4,coord,ey0:'
0547       if(ish.ge.5)write(ifch,*)
0548      *wp1,wm1,pnx,pny,sm,ic1,ic2,ic3,ic4,coord,ey0
0549 
0550         if(ic3.eq.0.and.ic4.eq.0)then
0551 
0552       ic(1)=ic1
0553       ic(2)=ic2
0554       nptl=nptl+1
0555       ep3(3)=pnx
0556       ep3(4)=pny
0557       ep3(2)=(wp1-wm1)/2
0558       ep3(1)=(wp1+wm1)/2
0559       call pstrans(ep3,ey0,-1)
0560       pptl(1,nptl)=ep3(3)
0561       pptl(2,nptl)=ep3(4)
0562       pptl(3,nptl)=ep3(2)
0563       pptl(4,nptl)=ep3(1)
0564       pptl(5,nptl)=sqrt(sm)
0565       idptl(nptl)=idtra(ic,0,0,3)
0566       iorptl(nptl)=1
0567       istptl(nptl)=0
0568       jorptl(nptl)=0
0569       do i=1,4
0570       xorptl(i,nptl)=coord(i)
0571       enddo
0572       tivptl(1,nptl)=coord(5)
0573       tivptl(2,nptl)=coord(6)
0574       ityptl(nptl)=40
0575 
0576        if(ish.ge.7)then
0577       write(ifch,*)'proj: nptl, mass**2,id',nptl,sm,idptl(nptl)
0578       write(ifch,*)'ept',
0579      *pptl(1,nptl),pptl(2,nptl),pptl(3,nptl),pptl(4,nptl)
0580        endif
0581 
0582         else
0583 
0584       ic(1)=ic1
0585       ic(2)=ic2
0586       nptl=nptl+1
0587       idptl(nptl)=idtra(ic,0,0,3)
0588       istptl(nptl)=20
0589       iorptl(nptl)=1
0590       jorptl(nptl)=0
0591       do i=1,4
0592       xorptl(i,nptl)=coord(i)
0593       enddo
0594       tivptl(1,nptl)=coord(5)
0595       tivptl(2,nptl)=coord(6)
0596       ityptl(nptl)=40
0597 
0598       ic(1)=ic3
0599       ic(2)=ic4
0600       idptl(nptl+1)=idtra(ic,0,0,3)
0601       istptl(nptl+1)=20
0602       iorptl(nptl+1)=1
0603       jorptl(nptl+1)=0
0604       do i=1,4
0605       xorptl(i,nptl+1)=coord(i)
0606       enddo
0607       tivptl(1,nptl+1)=coord(5)
0608       tivptl(2,nptl+1)=coord(6)
0609       ityptl(nptl+1)=40
0610 
0611       ep3(3)=pnx
0612       ep3(4)=pny
0613       ep3(2)=(wp1-wm1)/2
0614       ep3(1)=(wp1+wm1)/2
0615       call pstrans(ep3,ey0,-1)   !boost to hadronic c.m.s.
0616       ept(1)=ep3(3)
0617       ept(2)=ep3(4)
0618       ept(3)=ep3(2)
0619       ept(4)=ep3(1)
0620       do i=1,4
0621         ept1(i)=ep3(i)
0622       enddo
0623 
0624       sww=sqrt(sm)
0625       call psdeftr(sm,ept1,ey)
0626       ep(1)=.5*sww
0627       ep(2)=.5*sww
0628       ep(3)=0.
0629       ep(4)=0.
0630       call pstrans(ep,ey,1)
0631       pptl(1,nptl)=ep(3)
0632       pptl(2,nptl)=ep(4)
0633       pptl(3,nptl)=ep(2)
0634       pptl(4,nptl)=ep(1)
0635       do i=1,4
0636         pptl(i,nptl+1)=ept(i)-pptl(i,nptl)
0637       enddo
0638       nptl=nptl+1
0639         endif
0640 
0641       if(ish.ge.5)write(ifch,*)'fremny: final nptl',nptl
0642       call utprix('fremny',ish,ishini,5)
0643       return
0644       end
0645 
0646 c-----------------------------------------------------------------------
0647       subroutine psadis(iret)
0648 c-----------------------------------------------------------------------
0649 c psadis - DIS interaction
0650 c-----------------------------------------------------------------------
0651       double precision ept(4),ept1(4),xx,wpt(2),eprt,pl,plprt,psutz
0652      *,psuds
0653       dimension ep3(4),ey(3),ey0(3),bx(6),
0654      *qmin(2),iqc(2),nqc(2),ncc(2,2),gdv(2),gds(2),dfp(4)
0655       parameter (mjstr=20000)
0656       common /psar29/ eqj(4,mjstr),iqj(mjstr),ncj(2,mjstr),ioj(mjstr),nj
0657       common /psar30/ iorj(mjstr),ityj(mjstr),bxj(6,mjstr),q2j(mjstr)
0658       double precision pgampr,rgampr
0659       common/cgampr/pgampr(5),rgampr(4)
0660       parameter (ntim=1000)
0661       common/cprt/pprt(5,ntim),q2prt(ntim),idaprt(2,ntim),idprt(ntim)
0662      &,iorprt(ntim),jorprt(ntim),nprtj
0663       common/ciptl/iptl
0664       include 'epos.inc'
0665       include 'epos.incsem'
0666 
0667       call utpri('psadis',ish,ishini,3)
0668       if(ish.ge.3)write (ifch,*)'engy,elepti,iolept:'
0669       if(ish.ge.3)write (ifch,*)engy,elepti,iolept
0670       nptl=nptl+1
0671       idptl(nptl)=1220
0672       istptl(nptl)=1
0673       nptlh=nptl
0674       iptl=nptl
0675       s00=1.
0676 
0677       pptl(1,nptl)=0.
0678       pptl(2,nptl)=0.
0679       pptl(3,nptl)=-engy/2
0680       pptl(4,nptl)=engy/2
0681       pptl(5,nptl)=0
0682 
0683 1     continue
0684       if(iolept.eq.1)then
0685         wtot=engy**2
0686         engypr=wtot/4./elepti
0687         gdv01=psdh(ydmax*wtot,qdmin,iclpro,0)
0688         gdv02=psdh(ydmax*wtot,qdmin,iclpro,1)
0689         gds01=psdsh(ydmax*wtot,qdmin,iclpro,dqsh,0)
0690         gds02=psdsh(ydmax*wtot,qdmin,iclpro,dqsh1,1)
0691         gb0=(1.+(1.-ydmax)**2)*(gdv01+gds01)
0692      *  +2.*(1.-ydmax)*(gdv02+gds02)
0693 
0694 2       continue
0695         qq=qdmin*(qdmax/qdmin)**rangen()
0696         yd=ydmin*(ydmax/ydmin)**rangen()
0697         wd=yd*wtot
0698         if(ish.ge.4)write (ifch,*)'qq,wd,yd,ydmin,ydmax:'
0699         if(ish.ge.4)write (ifch,*)qq,wd,yd,ydmin,ydmax
0700         if(wd.lt.qq)goto 2
0701 
0702         gdv(1)=psdh(wd,qq,iclpro,0)
0703         gdv(2)=psdh(wd,qq,iclpro,1)
0704         gds(1)=psdsh(wd,qq,iclpro,dqsh,0)
0705         gds(2)=psdsh(wd,qq,iclpro,dqsh1,1)
0706         gbtr=(1.+(1.-yd)**2)*(gdv(1)+gds(1))
0707         gblong=2.*(1.-yd)*(gdv(2)+gds(2))
0708 c        gblong=0.   !???????
0709         gb=(gbtr+gblong)/gb0*.7
0710         if(ish.ge.4)then
0711           if(gb.gt.1.)write(ifmt,*)'gb,qq,yd,wd',gb,qq,yd,wd
0712           write (ifch,*)'gb,gdv,gds,gdv0,gds0,yd:'
0713           write (ifch,*)gb,gdv,gds,gdv01,gds01,
0714      *    gdv02,gds02,yd
0715         endif
0716         if(rangen().gt.gb)goto 2
0717 
0718         long=int(rangen()+gblong/(gbtr+gblong))
0719         elepto=qq/elepti/4.+elepti*(1.-yd)
0720         costhet=1.-qq/elepti/elepto/2.
0721         theta=acos(costhet)
0722         if(theta/pi*180..lt.themin)goto 2
0723         if(theta/pi*180..gt.themax)goto 2
0724         if(elepto.lt.elomin)goto 2
0725         if(ish.ge.3)write (ifch,*)'theta,elepto,elepti,iclpro:'
0726         if(ish.ge.3)write (ifch,*)theta/pi*180.,elepto,elepti,iclpro
0727         xbjevt=qq/wd
0728         qsqevt=qq
0729 
0730         call pscs(bcos,bsin)
0731         rgampr(1)=-elepto*sin(theta)*bcos
0732         rgampr(2)=-elepto*sin(theta)*bsin
0733         rgampr(3)=elepti-elepto*costhet
0734         rgampr(4)=elepti-elepto
0735 
0736         pgampr(1)=rgampr(1)
0737         pgampr(2)=rgampr(2)
0738         pgampr(3)=rgampr(3)-engypr
0739         pgampr(4)=rgampr(4)+engypr
0740         sm2=pgampr(4)*pgampr(4)
0741      *  -pgampr(1)*pgampr(1)-pgampr(2)*pgampr(2)-pgampr(3)*pgampr(3)
0742         pgampr(5)=sqrt(sm2)
0743         call utlob2(1,pgampr(1),pgampr(2),pgampr(3),pgampr(4),pgampr(5)
0744      *  ,rgampr(1),rgampr(2),rgampr(3),rgampr(4),40)
0745         if(ish.ge.4)write (ifch,*)'rgampr:',rgampr
0746 
0747       elseif(iolept.lt.0)then
0748 21      call lepexp(xbjevt,qsq)
0749         qq=qsq
0750         wd=qq/xbjevt
0751         if(qq.lt.qdmin.or.qq.gt.qdmax)goto21
0752 
0753         gdv(1)=psdh(wd,qq,iclpro,0)
0754         gdv(2)=psdh(wd,qq,iclpro,1)
0755         gds(1)=psdsh(wd,qq,iclpro,dqsh,0)
0756         gds(2)=psdsh(wd,qq,iclpro,dqsh1,1)
0757         yd=wd/engy**2
0758         gbtr=(1.+(1.-yd)**2)*(gdv(1)+gds(1))
0759         gblong=2.*(1.-yd)*(gdv(2)+gds(2))
0760         gblong=0. !????????????
0761         long=int(rangen()+gblong/(gbtr+gblong))
0762       else
0763         stop'wrong iolept'
0764       endif
0765       if(ish.ge.3)write (ifch,*)'qq,xbj,wd,gdv,gds,dqsh:'
0766       if(ish.ge.3)write (ifch,*)qq,xbjevt,wd,gdv,gds,dqsh
0767 
0768       egyevt=sqrt(wd-qq)
0769       pmxevt=.5*egyevt
0770 
0771       wp0=sqrt(qq)       !breit frame
0772       wm0=(wd-qq)/wp0
0773       ey0(1)=egyevt/wp0  !boost to the hadronic c.m.s.
0774       ey0(2)=1.
0775       ey0(3)=1.
0776       do i=1,6
0777         bx(i)=0.
0778       enddo
0779 
0780       if(long.eq.0)then
0781         sdmin=qq/(1.-sqrt(q2ini/qq))
0782         sqmin=sdmin
0783       else
0784         sdmin=4.*max(q2min,qcmass**2)+qq  !minimal mass for born
0785         xmm=(5.*sdmin-qq)/4.
0786         sqmin=1.1*(xmm+sqrt(xmm**2-qq*(sdmin-qq-4.*q2ini)))
0787      *  /2./(1.-4.*q2ini/(sdmin-qq))
0788       endif
0789       if(long.eq.1.and.wd.lt.1.001*sdmin)goto 1
0790 
0791       proja=210000.
0792       projb=0.
0793       call fremnu(ammin,proja,projb,proja,projb,
0794      *icp1,icp2,icp3,icp4)
0795 
0796       nj=0
0797 
0798       if((rangen().lt.gdv(long+1)/(gdv(long+1)+gds(long+1)).or.
0799      *egyevt.lt.1.0001*(ammin+sqrt(sdmin-qq))).and.
0800      *(long.eq.0.or.wd.gt.sqmin))then
0801         if(long.eq.0)then
0802           xd=qq/wd
0803           tu=psdfh4(xd,q2min,0.,iclpro,1)/2.25
0804           td=psdfh4(xd,q2min,0.,iclpro,2)/9.
0805           gdv0=(tu+td)*4.*pi**2*alfe/qq
0806      *    *sngl(psuds(qq,1)/psuds(q2min,1))
0807           if(ish.ge.4)write (ifch,*)'gdv0:',gdv0,sdmin
0808 
0809           if(rangen().lt.gdv0/gdv(1).or.wd.le.1.0001*sdmin)then    !?????
0810             if(ish.ge.3)write (ifch,*)'no cascade,gdv0,gdv',gdv0,gdv
0811             if(rangen().lt.tu/(tu+td))then
0812               iq1=1
0813               izh=3
0814             else
0815               iq1=2
0816               izh=6
0817             endif
0818             jq=1
0819             if(ish.ge.8)write (ifch,*)'before call timsh2: ',
0820      *      'qq,egyevt,iq1',qq,egyevt,iq1
0821             call timsh2(qq,0.,egyevt,iq1,-iq1,iq1,-iq1)
0822 
0823             nj=nj+1
0824             iqj(nj)=izh
0825             nqc(1)=nj
0826             nqc(2)=0
0827 
0828             ep3(1)=pprt(4,2)
0829             ep3(2)=pprt(3,2)
0830             ep3(3)=0.
0831             ep3(4)=0.
0832             call pstrans(ep3,ey0,1)
0833             do i=1,4
0834               eqj(i,nj)=ep3(i)
0835             enddo
0836             s0h=0.
0837             c0h=1.
0838             s0xh=0.
0839             c0xh=1.
0840             call psreti(nqc,jq,1,ey0,s0xh,c0xh,s0h,c0h)
0841             goto 17
0842           endif
0843         endif
0844 
0845         call psdint(wd,qq,sds0,sdn0,sdb0,sdt0,sdr0,1,long)
0846         if(ish.ge.3)write (ifch,*)'wd,qq,sds0,sdn0,sdb0,sdt0,sdr0:'
0847         if(ish.ge.3)write (ifch,*)wd,qq,sds0,sdn0,sdb0,sdt0
0848         gb10=(sdn0+sdt0)*(1.-qq/wd)
0849         xdmin=sqmin/wd
0850 
0851 3       continue
0852         xd=(xdmin-qq/wd)/((xdmin-qq/wd)/(1.-qq/wd))
0853      *  **rangen()+qq/wd
0854         call psdint(xd*wd,qq,sds,sdn,sdb,sdt,sdr,1,long)
0855         if(ish.ge.3)write (ifch,*)'wdhard,qq,sds,sdn,sdb,sdt:'
0856         if(ish.ge.3)write (ifch,*)xd*wd,qq,sds,sdn,sdb,sdt
0857         tu=psdfh4(xd,q2min,0.,iclpro,1)
0858         td=psdfh4(xd,q2min,0.,iclpro,2)
0859         gb1=(sdn*(tu/2.25+td/9.)+sdt*(tu+td)/4.5)
0860      *  *(1.-qq/wd/xd)/gb10
0861         if(gb1.gt.1..and.ish.ge.1)write(ifmt,*)'gb1,xd,wd,qq,sdt0,sdt',
0862      *   gb1,xd,wd,qq,sdt0,sdt
0863         if(ish.ge.6)write (ifch,*)'gb1,xd,wd,qq,sdt0,sdt:'
0864         if(ish.ge.6)write (ifch,*)gb1,xd,wd,qq,sdt0,sdt
0865         if(rangen().gt.gb1)goto 3
0866 
0867         gdres=(sdt-sds)/4.5
0868         gdrga=sdr/4.5
0869         gdsin=sds/4.5
0870         dtu=tu*(sdn/2.25+sdt/4.5)
0871         dtd=td*(sdn/9.+sdt/4.5)
0872         if(rangen().lt.dtu/(dtu+dtd))then
0873           iq1=1
0874           izh=3
0875           gdbor=sdb/2.25
0876           gdnon=sdn/2.25
0877         else
0878           iq1=2
0879           izh=6
0880           gdbor=sdb/9.
0881           gdnon=sdn/9.
0882         endif
0883 
0884         wpi=wp0
0885         wmi=(xd*wd-qq)/wpi
0886         iqc(2)=iq1
0887         nj=nj+1
0888         iqj(nj)=izh
0889         eqj(1,nj)=.5*(wm0-wmi)
0890         eqj(2,nj)=-eqj(1,nj)
0891         eqj(3,nj)=0.
0892         eqj(4,nj)=0.
0893         ncc(1,2)=nj
0894         ncc(2,2)=0
0895         if(ish.ge.3)write (ifch,*)'wp0,wm0,wpi,wmi,iqc(2),eqj'
0896         if(ish.ge.3)write (ifch,*)wp0,wm0,wpi,wmi,iqc(2),eqj(2,nj)
0897 
0898       else
0899         xdmin=sdmin/wd
0900         xpmax=((egyevt-ammin)**2+qq)/wd
0901         iq1=int(3.*rangen()+1.)*(2.*int(.5+rangen())-1.)
0902 
0903         aks=rangen()
0904         if(long.eq.0.and.aks.lt.dqsh/gds(1).and.
0905      *  egyevt.gt.ammin+sqrt(s00))then
0906           if(ish.ge.3)write (ifch,*)'no cascade for q_s',
0907      *    aks,dqsh/gds(1)
0908           xd=qq/wd
0909           xpmin=xd+s00/wd
0910           jcasc=0
0911           if(iq1.gt.0)then
0912             jq=1
0913           else
0914             jq=2
0915           endif
0916         else
0917           jcasc=1
0918           call psdint(xpmax*wd,qq,sds0,sdn0,sdb0,sdt0,sdr0,0,long)
0919           call psdint(xpmax*wd,qq,sdsq0,sdnq0,sdbq0,sdtq0,sdrq0,1,long)
0920           if(ish.ge.3)write (ifch,*)
0921      *    'xpmax*wd,qq,sds0,sdn0,sdb0,sdt0,sdr0:'
0922           if(ish.ge.3)write (ifch,*)
0923      *    xpmax*wd,qq,sds0,sdn0,sdb0,sdt0,sdr0
0924         gb10=sdt0*fzeroGluZZ(0.,iclpro)+(sdnq0+sdtq0)
0925      *        *fzeroSeaZZ(0.,iclpro)
0926           gb10=gb10*15.
0927 
0928 4         xd=xdmin*(xpmax/xdmin)**rangen()
0929           xpmin=xd
0930           call psdint(xd*wd,qq,sds,sdn,sdb,sdt,sdr,0,long)
0931           call psdint(xd*wd,qq,sdsq,sdnq,sdbq,sdtq,sdrq,1,long)
0932           if(ish.ge.3)write (ifch,*)'xd*wd,qq,sds,sdn,sdb,sdt,sdr:'
0933           if(ish.ge.3)write (ifch,*)xd*wd,qq,sds,sdn,sdb,sdt,sdr
0934           wwg=sdt*fzeroGluZZ(xd,iclpro)
0935           wwq=(sdnq+sdtq)*fzeroSeaZZ(xd,iclpro)
0936           gb12=(wwq+wwg)/gb10*(xpmax/xd)**dels
0937           if(gb12.gt.1..and.ish.ge.1)write(ifmt,*)
0938      *    'gb12,xpmax*wd,xd*wd,sdt0,sdnq0+sdtq0,sdt,sdnq+sdtq',
0939      *    gb12,xpmax*wd,xd*wd,sdt0,sdnq0+sdtq0,sdt,sdnq+sdtq,
0940      *    wwq,wwg,(xpmax/xd)**dels,gb10
0941           if(ish.ge.5)write (ifch,*)'gb12,xd,xpmax,wwq,wwg:'
0942           if(ish.ge.5)write (ifch,*)gb12,xd,xpmax,wwq,wwg
0943           if(rangen().gt.gb12)goto 4
0944         endif
0945 
0946         if(jcasc.ne.0)then
0947           gb20=(1.-xd/xpmax)**betpom*sdt*(1.-glusea)+
0948      *    EsoftQZero(xd/xpmax)*(sdnq+sdtq)*glusea
0949         else
0950           gb20=EsoftQZero(xd/xpmax)
0951         endif
0952         if(1.+2.*(-alpqua)+dels.ge.0.)then
0953           xpminl=(1.-xpmax)**(alplea(iclpro)+1.)
0954           xpmaxl=(1.-xpmin)**(alplea(iclpro)+1.)
0955 
0956 5         xp=1.-(xpminl+(xpmaxl-xpminl)*rangen())**
0957      *    (1./(alplea(iclpro)+1.))
0958           if(jcasc.ne.0)then
0959             gb2=((1.-xd/xp)**betpom*sdt*(1.-glusea)+
0960      *      EsoftQZero(xd/xp)*(sdnq+sdtq)*glusea)*(xp/xpmax)**
0961      *      (1.+2.*(-alpqua)+dels)/gb20
0962           else
0963            gb2=EsoftQZero(xd/xp)*(xp/xpmax)**(1.+2.*(-alpqua)+dels)/gb20
0964           endif
0965           if(gb2.gt.1..and.ish.ge.1)then
0966             write(ifmt,*)'gb2,xp:',gb2,xp
0967 c            read (*,*)
0968           endif
0969           if(rangen().gt.gb2)goto 5
0970         else
0971           xpmaxl=xpmax**(2.+2.*(-alpqua)+dels)
0972           xpminl=xpmin**(2.+2.*(-alpqua)+dels)
0973 
0974 6         xp=(xpminl+(xpmaxl-xpminl)*rangen())**
0975      *    (1./(2.+2.*(-alpqua)+dels))
0976           if(jcasc.ne.0)then
0977             gb21=((1.-xd/xp)**betpom*sdt*(1.-glusea)+
0978      *      EsoftQZero(xd/xp)*(sdnq+sdtq)*glusea)*
0979      *      ((1.-xp)/(1.-xd))**alplea(iclpro)/gb20
0980           else
0981           gb21=EsoftQZero(xd/xp)*((1.-xp)/(1.-xd))**alplea(iclpro)/gb20
0982           endif
0983           if(gb21.gt.1..and.ish.ge.1)then
0984             write(ifmt,*)'gb21,xp:',gb21,xp
0985 c            read (*,*)
0986           endif
0987           if(rangen().gt.gb21)goto 6
0988         endif
0989 
0990         wwh=xd*wd-qq
0991         wwsh=xp*wd-qq
0992         ammax=(egyevt-sqrt(wwsh))**2
0993 22      call fremnx(ammax,ammin,sm,icp3,icp4,iret)
0994         if(iret.ne.0.and.ish.ge.1)write(ifmt,*)'iret.ne.0!'
0995      *                                         ,ammax,ammin**2
0996         wmn=(1.-xp)*wd/wp0
0997         wpn=sm/wmn
0998         pnx=0.
0999         pny=0.
1000         wpp=wp0-wpn
1001         wmp=wm0-wmn
1002         if(ish.ge.5)write(ifch,*)'wp0,wm0,wpn,wmn,wpp,wmp:'
1003         if(ish.ge.5)write(ifch,*)wp0,wm0,wpn,wmn,wpp,wmp
1004 
1005         if(jcasc.eq.0.or.rangen().lt.wwq/(wwg+wwq).
1006      *  and.xd*wd.gt.sqmin.and.wwsh.gt.
1007      *  (sqrt(wwh)+sqrt(s00))**2)then
1008           zgmin=xd/xp
1009           zgmax=1./(1.+wp0/xd/wd/(wpp-wwh/wmp))
1010           if(zgmin.gt.zgmax)goto 22
1011 23        zg=zgmin-rangen()*(zgmin-zgmax)
1012           if(rangen().gt.zg**dels*((1.-xd/xp/zg)/ (1.-xd/xp))**betpom)
1013      *    goto 23
1014           xg=xd/zg             !w- share for the struck quark
1015           wmq=wd/wp0*(xg-xd)   !w- for its counterpart
1016           wpq=s00/wmq          !1. gev^2 / wmq
1017           wmq=0.
1018           wpp=wpp-wpq
1019           wmp=wmp-wmq
1020           sxx=wpp*wmp
1021           if(ish.ge.5)write (ifch,*)'wpq,wmq,wpp,wmp,sxx:'
1022           if(ish.ge.5)write (ifch,*)wpq,wmq,wpp,wmp,sxx
1023 
1024           if(jcasc.eq.0)then
1025             if(ish.ge.6)write (ifch,*)'before call timsh2: qq,sxx,iq1',
1026      *      qq,sxx,iq1
1027             call timsh2(qq,0.,sqrt(sxx),iq1,-iq1,iq1,-iq1)
1028             ept(1)=.5*(wpp+wmp)
1029             ept(2)=.5*(wpp-wmp)
1030             ept(3)=0.
1031             ept(4)=0.
1032             call psdeftr(sxx,ept,ey)
1033             ep3(1)=pprt(4,2)
1034             ep3(2)=pprt(3,2)
1035             ep3(3)=0.
1036             ep3(4)=0.
1037 
1038             call pstrans(ep3,ey,1)
1039             wmp=ep3(1)-ep3(2)
1040             goto 24
1041           endif
1042         else
1043           iq1=0
1044           sxx=wpp*wmp
1045         endif
1046 
1047         if(ish.ge.3)write (ifch,*)'wwh,wwsh,sxx,wpp,wmp:',
1048      *  wwh,wwsh,sxx,wpp,wmp
1049 
1050         wpi=wpp
1051         wmi=wwh/wpp
1052         wmp=wmp-wmi
1053 24      call fremny(wpn,wmn,pnx,pny,sm,icp1,icp2,icp3,icp4,bx,ey0)
1054 
1055         if((-alpqua).eq.-1.)stop'dis does not work for 1/x'
1056 25      aks=rangen()
1057         z=.5*aks**(1./(1.+(-alpqua)))
1058         if(z.lt.1.e-5.or.rangen().gt.(2.*(1.-z))**(-alpqua))goto 25
1059         if(rangen().gt..5)z=1.-z
1060         wm2=wmp*z
1061         wm1=wmp-wm2
1062 
1063         iqc(2)=iq1
1064         nj=nj+1
1065         iqj(nj)=-int(2.*rangen()+1.)
1066         iqj(nj+1)=-iqj(nj)
1067         eqj(1,nj)=.5*wm1
1068         eqj(2,nj)=-eqj(1,nj)
1069         eqj(3,nj)=0.
1070         eqj(4,nj)=0.
1071         eqj(1,nj+1)=.5*wm2
1072         eqj(2,nj+1)=-eqj(1,nj+1)
1073         eqj(3,nj+1)=0.
1074         eqj(4,nj+1)=0.
1075         nj=nj+1
1076 
1077         if(iq1.eq.0)then
1078           ncc(1,2)=nj-1
1079           ncc(2,2)=nj
1080           gdres=sdt-sds
1081           gdrga=sdr
1082           gdsin=sds
1083           gdbor=sdb
1084           gdnon=sdn
1085         else
1086           nj=nj+1
1087           if(iabs(iq1).eq.3)then
1088             iqj(nj)=-iq1*4/3
1089           else
1090             iqj(nj)=-iq1
1091           endif
1092           eqj(1,nj)=.5*(wpq+wmq)
1093           eqj(2,nj)=.5*(wpq-wmq)
1094           eqj(3,nj)=0.
1095           eqj(4,nj)=0.
1096           if(iq1.gt.0)then
1097             ncj(1,nj)=nj-1
1098             ncj(1,nj-1)=nj
1099             ncj(2,nj)=0
1100             ncj(2,nj-1)=0
1101           else
1102             ncj(1,nj)=nj-2
1103             ncj(1,nj-2)=nj
1104             ncj(2,nj)=0
1105             ncj(2,nj-2)=0
1106           endif
1107 
1108           if(jcasc.eq.0)then
1109             if(iq1.gt.0)then
1110               nqc(1)=nj-2
1111               nqc(2)=0
1112             else
1113               nqc(1)=nj-1
1114               nqc(2)=0
1115             endif
1116             s0h=0.
1117             c0h=1.
1118             s0xh=0.
1119             c0xh=1.
1120             call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
1121             goto 17
1122           else
1123             gdres=(sdtq-sdsq)/4.5
1124             gdrga=sdrq/4.5
1125             gdsin=sdsq/4.5
1126             gdbor=sdbq/4.5
1127             gdnon=sdnq/4.5
1128             if(iq1.gt.0)then
1129               ncc(1,2)=nj-2
1130               ncc(2,2)=0
1131             else
1132               ncc(1,2)=nj-1
1133               ncc(2,2)=0
1134             endif
1135           endif
1136         endif
1137 
1138         if(ish.ge.3)write (ifch,*)'wpn,wmn,wpi,wmi,wm1,wm2,nj'
1139         if(ish.ge.3)write (ifch,*)wpn,wmn,wpi,wmi,wm1,wm2,nj
1140       endif
1141 
1142       si=wpi*wmi+qq
1143       qmin(2)=q2min                 !effective momentum cutoff below
1144       s2min=max(4.*qq,16.*q2min)    !mass cutoff for born scattering
1145 
1146       if(rangen().gt.gdres/(gdres+gdsin+gdnon).or.
1147      *si.lt.(s2min+qq))goto 12
1148 
1149 c---------------------------------------
1150 c hard pomeron (resolved photon)
1151 c---------------------------------------
1152       if(ish.ge.3)write(ifmt,*)'resolved,gdrga,gdres',gdrga,gdres
1153 
1154       jj=1
1155       if(rangen().gt.gdrga/gdres.and.si.gt.1.1*s2min+qq)then
1156         if(ish.ge.3)write(ifmt,*)'dir-res,si,qq',si,qq
1157         pt=0.
1158         pt2=0.
1159         iqc(1)=0
1160         ept(1)=.5*(wpi+wmi)
1161         ept(2)=.5*(wpi-wmi)
1162         ept(3)=0.
1163         ept(4)=0.
1164         wpt(1)=wpi           !lc+ for the current jet emission
1165         wpt(2)=si/wpi        !lc- for the current jet emission
1166 
1167         qqmin=max(q2min,s2min/(si/qq-1.))
1168         qqmax=min(si/2.,si-s2min)
1169         qmin(1)=qqmin
1170         xmax=1.
1171         xmin=(s2min+qq)/si
1172         if(qqmin.ge.qqmax.or.xmin.ge.xmax)stop'min>max'
1173         gb0=psjti(qmin(1),qq,si-qq,7,iqc(2),1)*psfap(1.d0,0,1)
1174 
1175         ncc(1,1)=0
1176         ncc(2,1)=0
1177         jgamma=1
1178         ntry=0
1179         xmin1=0.
1180         xmin2=0.
1181         xmax1=0.
1182         djl=0.
1183         goto 9
1184       else
1185         if(ish.ge.3)write(ifmt,*)'res,si,qq',si,qq
1186         qmin(1)=q2min                 !effective momentum cutoff above
1187         si=si-qq
1188         zmin=s2min/si
1189         dft0=psdfh4(zmin,q2min,qq,0,0)*psjti(q2min,qq,si,0,iqc(2),1)
1190      *  +(psdfh4(zmin,q2min,qq,0,1)+psdfh4(zmin,q2min,qq,0,2)+
1191      *  psdfh4(zmin,q2min,qq,0,3))*psjti(q2min,qq,si,7,iqc(2),1)
1192 
1193 7       continue
1194         z=zmin**rangen()
1195         do i=1,4
1196           dfp(i)=psdfh4(z,q2min,qq,0,i-1)
1197         enddo
1198         dfp(1)=dfp(1)*psjti(q2min,qq,z*si,0,iqc(2),1)
1199         dfptot=dfp(1)
1200         if(iqc(2).eq.0)then
1201           sjq=psjti(q2min,qq,z*si,1,0,1)
1202           do i=2,4
1203             dfp(i)=dfp(i)*sjq
1204             dfptot=dfptot+dfp(i)
1205           enddo
1206         else
1207           sjqqp=psjti(q2min,qq,z*si,1,2,1)
1208           do i=2,4
1209             if(iabs(iqc(2)).eq.i-1)then
1210               dfp(i)=dfp(i)*(psjti(q2min,qq,z*si,1,1,1)+
1211      *        psjti(q2min,qq,z*si,1,-1,1))/2.
1212             else
1213               dfp(i)=dfp(i)*sjqqp
1214             endif
1215             dfptot=dfptot+dfp(i)
1216           enddo
1217         endif
1218 
1219         if(rangen().gt.dfptot/dft0)goto 7
1220 
1221         wpq=wpi*(1.-z)
1222         wpi=wpi*z
1223         aks=dfptot*rangen()
1224         if(aks.lt.dfp(1))then
1225           iqc(1)=0
1226           nj=nj+1
1227           ncc(1,1)=nj
1228           ncc(2,1)=nj+1
1229 
1230           iqj(nj)=-int(2.*rangen()+1.)
1231           iqj(nj+1)=-iqj(nj)
1232           wpq1=wpq*rangen()
1233           eqj(1,nj)=.5*wpq1
1234           eqj(2,nj)=eqj(1,nj)
1235           eqj(3,nj)=0.
1236           eqj(4,nj)=0.
1237           eqj(1,nj+1)=.5*(wpq-wpq1)
1238           eqj(2,nj+1)=eqj(1,nj+1)
1239           eqj(3,nj+1)=0.
1240           eqj(4,nj+1)=0.
1241           nj=nj+1
1242 
1243         else
1244           if(aks.lt.dfp(1)+dfp(2))then
1245             iqc(1)=1
1246           elseif(aks.lt.dfp(1)+dfp(2)+dfp(3))then
1247             iqc(1)=2
1248           else
1249             iqc(1)=3
1250           endif
1251           iqc(1)=iqc(1)*(2*int(2.*rangen())-1)
1252           nj=nj+1
1253           ncc(1,1)=nj
1254           ncc(2,1)=0
1255 
1256           iqj(nj)=-iqc(1)
1257           eqj(1,nj)=.5*wpq
1258           eqj(2,nj)=eqj(1,nj)
1259           eqj(3,nj)=0.
1260           eqj(4,nj)=0.
1261         endif
1262 
1263         ept(1)=.5*(wpi+wmi)
1264         ept(2)=.5*(wpi-wmi)
1265         ept(3)=0.
1266         ept(4)=0.
1267         jgamma=0
1268         ntry=0
1269       endif
1270 
1271 8     continue
1272 
1273 c ladder rung
1274 c---------------------------------------
1275       pt2=ept(3)**2+ept(4)**2
1276       pt=sqrt(pt2)
1277 
1278       wpt(1)=ept(1)+ept(2)              !lc+ for the current jet emissi
1279       wpt(2)=ept(1)-ept(2)              !lc- for the current jet emissi
1280 
1281       s2min=max(qmin(1),16.*qmin(2))    !mass cutoff for born
1282       s2min=max(s2min,4.*qq)
1283 
1284       if(jj.eq.1)then
1285         wwmin=2.*s2min-2.*pt*sqrt(q2ini)
1286         wwmin=(wwmin+sqrt(wwmin**2+4.*pt2*(s2min-q2ini)))
1287      *  /(1.-q2ini/s2min)/2.
1288         sj=psjti(qmin(1),qq,si,iqc(1),iqc(2),1)   !total jet
1289         sj2=psjti1(q2min,qmin(1),qq,si,iqc(2),iqc(1),1)
1290         if(ish.ge.3)write(ifch,*)'resolved - si,wwmin,s2min,sj,sj2:'
1291         if(ish.ge.3)write(ifch,*)si,wwmin,s2min,sj,sj2
1292         if(sj.eq.0.)stop'sj=0'
1293         if(rangen().gt.sj2/sj.and.si.gt.1.1*wwmin)goto 26
1294         jj=2
1295       endif
1296       sj=psjti1(qmin(2),qmin(1),qq,si,iqc(2),iqc(1),1)
1297       sjb=psbint(qmin(1),qmin(2),qq,si,iqc(1),iqc(2),1) !born parton-parton
1298       wwmin=17./16*s2min-2.*pt*sqrt(q2ini)
1299       wwmin=(wwmin+sqrt(wwmin**2+pt2*(s2min/4.-4.*q2ini)))
1300      */(1.-16.*q2ini/s2min)/2.
1301       if(rangen().lt.sjb/sj.or.si.lt.1.1*wwmin)goto 10
1302 
1303 26    continue
1304       wpt(jj)=wpt(jj)-pt2/wpt(3-jj)
1305 
1306       if(jj.eq.1)then
1307         discr=(si+2.*pt*sqrt(q2ini))**2-4.*q2ini*(2.*si+pt2)
1308         if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr,si,pt,wwmin',
1309      *  discr,si,pt,wwmin
1310         discr=sqrt(discr)
1311         qqmax=(si+2.*pt*sqrt(q2ini)+discr)/2./(2.+pt2/si)
1312       else
1313         discr=(si+2.*pt*sqrt(q2ini))**2-4.*q2ini*(17.*si+pt2)
1314         if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr,si,pt,wwmin',
1315      *  discr,si,pt,wwmin
1316         discr=sqrt(discr)
1317         qqmax=(si+2.*pt*sqrt(q2ini)+discr)/2./(17.+pt2/si)
1318       endif
1319       qqmin=2.*q2ini*si/(si+2.*pt*sqrt(q2ini)+discr)
1320       if(jj.eq.1.and.s2min.gt.qqmin.or.
1321      *jj.eq.2.and.s2min.gt.16.*qqmin)then
1322         xmm=.5*(si-s2min+2.*pt*sqrt(q2ini))
1323         discr=xmm**2-q2ini*(si+pt2)
1324         if(discr.lt.0..and.ish.ge.1)write(ifmt,*)'discr1,si,pt,wwmin',
1325      *  discr,si,pt,wwmin
1326         qqmin=q2ini*si/(xmm+sqrt(discr))
1327       endif
1328 
1329       xmin=1.-q2ini/qqmin
1330       xmax=1.-q2ini/qqmax
1331       if(ish.ge.6)write(ifch,*)'qqmin,qqmax,xmin,xmax',
1332      *qqmin,qqmax,xmin,xmax
1333       if(qqmin.lt.qmin(jj))then
1334         qqmin=qmin(jj)
1335         xmi=max(1.-((pt*sqrt(qqmin)+sqrt(pt2*qqmin+
1336      *  si*(si-s2min-qqmin*(1.+pt2/si))))/si)**2,
1337      *  (s2min+qqmin*(1.+pt2/si)-2.*pt*sqrt(qqmin))/si)
1338         xmin=max(xmin,xmi)
1339         if(xmin.le.0.)xmin=(s2min+qqmin*(1.+pt2/si))/si
1340         if(ish.ge.6)write(ifch,*)'qqmin,qmin(jj),xmin,s2min',
1341      *  qqmin,qmin(jj),xmin,s2min
1342       endif
1343 
1344       qm0=qmin(jj)
1345       xm0=1.-q2ini/qm0
1346       if(xm0.gt.xmax.or.xm0.lt.xmin)then
1347         xm0=.5*(xmax+xmin)
1348       endif
1349 c      s2max=xm0*si
1350       s2max=xm0*si-qm0*(1.+pt2/si)+2.*pt*sqrt(q2ini)  !new ladder mass squared
1351       xx=xm0
1352 
1353       if(jj.eq.1)then
1354         sj0=psjti(qm0,qq,s2max,0,iqc(2),1)*psfap(xx,iqc(1),0)+
1355      *  psjti(qm0,qq,s2max,7,iqc(2),1)*psfap(xx,iqc(1),1)
1356         gb0=sj0/log(q2ini/qcdlam)*sngl(psuds(qm0,iqc(1)))*qm0*2.
1357       else
1358         sj0=psjti1(qm0,qmin(1),qq,s2max,0,iqc(1),1)*psfap(xx,iqc(2),0)
1359      *  +psjti1(qm0,qmin(1),qq,s2max,7,iqc(1),1)*psfap(xx,iqc(2),1)
1360         gb0=sj0/log(q2ini/qcdlam)*sngl(psuds(qm0,iqc(2)))*qm0*2.
1361       endif
1362       if(gb0.le.0.)then
1363          write(ifmt,*)'gb0.le.0.  si,qq,pt2:',si,qq,pt2
1364          iret=1
1365          goto 9999
1366       endif
1367       if(xm0.le..5)then
1368         gb0=gb0*xm0**(1.-delh)
1369       else
1370         gb0=gb0*(1.-xm0)*2.**delh
1371       endif
1372 
1373       xmin2=max(.5,xmin)
1374       xmin1=xmin**delh                 !xmin, xmax are put into powe
1375       xmax1=min(xmax,.5)**delh       !to simulate x value below
1376       if(xmin.ge..5)then
1377         djl=1.
1378       elseif(xmax.lt..5)then
1379         djl=0.
1380       else
1381         djl=1./(1.+((2.*xmin)**delh-1.)/delh/
1382      *  log(2.*(1.-xmax)))
1383       endif
1384 
1385       ntry=0
1386 9     continue
1387       ntry=ntry+1
1388       if(ntry.ge.10000)then
1389         print *,"ntry.ge.10000"
1390         iret=1
1391         goto 9999
1392       endif
1393       if(jgamma.ne.1)then
1394         if(rangen().gt.djl)then        !lc momentum share in the cur
1395           x=(xmin1+rangen()*(xmax1-xmin1))**(1./delh)
1396         else
1397           x=1.-(1.-xmin2)*((1.-xmax)/(1.-xmin2))**rangen()
1398         endif
1399         q2=qqmin/(1.+rangen()*(qqmin/qqmax-1.))
1400         qt2=q2*(1.-x)
1401         if(ish.ge.6)write(ifch,*)'jj,q2,x,qt2',jj,q2,x,qt2
1402         if(qt2.lt.q2ini)goto 9
1403       else
1404         x=xmin+rangen()*(xmax-xmin)
1405         q2=qqmin*(qqmax/qqmin)**rangen()
1406         qt2=(q2-x*qq)*(1.-x)
1407         if(ish.ge.6)write(ifch,*)'jj,q2,x,qt2',jj,q2,x,qt2
1408         if(qt2.lt.0.)goto 9
1409       endif
1410 
1411       qt=sqrt(qt2)
1412       call pscs(bcos,bsin)
1413 c ep3 is now 4-vector for s-channel gluon produced in current ladder run
1414       ep3(3)=qt*bcos
1415       ep3(4)=qt*bsin
1416       ptnew=(ept(3)-ep3(3))**2+(ept(4)-ep3(4))**2
1417       if(jj.eq.1)then
1418         s2min2=max(q2,s2min)
1419       else
1420         s2min2=max(s2min,16.*q2)
1421       endif
1422 
1423       if(jgamma.ne.1)then
1424         s2=x*si-q2*(1.+pt2/si)-ptnew+pt2+qt2  !new ladder mass squared
1425         if(s2.lt.s2min2)goto 9      !rejection in case of too low mass
1426         xx=x
1427 
1428         if(jj.eq.1)then
1429           sj1=psjti(q2,qq,s2,0,iqc(2),1)
1430           if(iqc(1).ne.0)then
1431             sj2=psjti(q2,qq,s2,iqc(1),iqc(2),1)
1432           elseif(iqc(2).eq.0)then
1433             sj2=psjti(q2,qq,s2,1,0,1)
1434           else
1435             sj2=psjti(q2,qq,s2,1,1,1)/6.+
1436      *      psjti(q2,qq,s2,-1,1,1)/6.+
1437      *      psjti(q2,qq,s2,2,1,1)/1.5
1438           endif
1439         else
1440           sj1=psjti1(q2,qmin(1),qq,s2,0,iqc(1),1)
1441           if(iqc(2).ne.0)then
1442             sj2=psjti1(q2,qmin(1),qq,s2,iqc(2),iqc(1),1)
1443           elseif(iqc(1).eq.0)then
1444             sj2=psjti1(q2,qmin(1),qq,s2,1,0,1)
1445           else
1446             sj2=psjti1(q2,qmin(1),qq,s2,1,1,1)/6.+
1447      *      psjti1(q2,qmin(1),qq,s2,-1,1,1)/6.+
1448      *      psjti1(q2,qmin(1),qq,s2,2,1,1)/1.5
1449           endif
1450         endif
1451 c gb7 is the rejection function for x and q**2 simulation
1452         gb7=(sj1*psfap(xx,iqc(jj),0)+sj2*psfap(xx,iqc(jj),1))
1453      *  /log(qt2/qcdlam)*sngl(psuds(q2,iqc(jj)))*q2/gb0
1454 
1455         if(x.le..5)then
1456           gb7=gb7*x**(1.-delh)
1457         else
1458           gb7=gb7*(1.-x)*2.**delh
1459         endif
1460       else
1461         s2=x*si-q2               !new ladder mass squared
1462         if(s2.lt.s2min2)goto 9   !rejection in case of too low mass
1463 
1464         sj1=0.
1465         xx=x
1466         if(iqc(2).eq.0)then
1467           sj2=psjti(q2,qq,s2,1,0,1)
1468         else
1469           sj2=psjti(q2,qq,s2,1,1,1)/naflav/2.+
1470      *    psjti(q2,qq,s2,-1,1,1)/naflav/2.+
1471      *    psjti(q2,qq,s2,2,1,1)*(1.-1./naflav)
1472         endif
1473         gb7=sj2*psfap(xx,0,1)/gb0  !????*(1.-x*qq/q2)
1474       endif
1475       if(gb7.gt.1..or.gb7.lt.0..and.ish.ge.1)write(ifmt,*)'gb7,q2,x,gb0'
1476      *,gb7,q2,x,gb0
1477       if(rangen().gt.gb7)goto 9
1478 
1479       if(ish.ge.6)write(ifch,*)'res: jj,iqc,ncc:',
1480      *jj,iqc(jj),ncc(1,jj),ncc(2,jj)
1481 
1482       nqc(2)=0
1483       iqnew=iqc(jj)
1484       if(jgamma.ne.1)then
1485         if(rangen().lt.sj1/(sj1+sj2))then
1486           if(iqc(jj).eq.0)then
1487             jt=1
1488             jq=int(1.5+rangen())
1489             nqc(1)=ncc(jq,jj)
1490           else
1491             jt=2
1492             if(iqc(jj).gt.0)then
1493               jq=1
1494             else
1495               jq=2
1496             endif
1497             nqc(1)=0
1498             iqnew=0
1499           endif
1500           iq1=iqc(jj)
1501         else
1502           if(iqc(jj).ne.0)then
1503             iq1=0
1504             jt=3
1505             if(iqc(jj).gt.0)then
1506               jq=1
1507             else
1508               jq=2
1509             endif
1510             nqc(1)=ncc(1,jj)
1511           else
1512             jt=4
1513             jq=int(1.5+rangen())
1514             iq1=int(naflav*rangen()+1.)*(3-2*jq)
1515             nqc(1)=ncc(jq,jj)
1516             iqnew=-iq1
1517           endif
1518         endif
1519       else
1520         jt=5
1521         jq=int(1.5+rangen())
1522         iq1=int(naflav*rangen()+1.)*(3-2*jq)
1523         iqnew=-iq1
1524         nqc(1)=0
1525       endif
1526       eprt=max(1.d0*qt,
1527      *.5d0*((1.d0-x)*wpt(jj)+qt2/(1.d0-x)/wpt(jj)))
1528       pl=((1.d0-x)*wpt(jj)-eprt)*(3-2*jj)
1529       zeta=sqrt(qt2/si)/sqrt(x*(1.-x))
1530       if(iq1.eq.0)then
1531         iq2ini=9
1532         jo=iq2ini
1533         if(zeta.gt.zetacut)jo=-jo
1534       else
1535         iq2ini=iq1
1536         jo=iq2ini
1537       endif
1538 27    call timsh1(q2,sngl(eprt),iq2ini,jo)
1539       amprt=pprt(5,1)**2
1540       plprt=eprt**2-amprt-qt2
1541       if(plprt.lt.-1d-6)goto 27
1542       ep3(1)=eprt
1543       ep3(2)=dsqrt(max(0.d0,plprt))
1544       if(pl.lt.0.d0)ep3(2)=-ep3(2)
1545       ey(1)=1.
1546       ey(2)=1.
1547       ey(3)=1.
1548       do i=1,4
1549         ept1(i)=ept(i)-ep3(i)
1550       enddo
1551       call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
1552       if(ish.ge.6)then
1553       write(ifch,*)'q2,amprt,qt2',q2,amprt,qt2
1554       write(ifch,*)'eprt,plprt',eprt,plprt
1555       write(ifch,*)'ep3',ep3
1556       write(ifch,*)'ept',ept
1557       write(ifch,*)'ept1',ept1
1558       endif
1559       s2new=psnorm(ept1)
1560 
1561       if(s2new.gt.s2min2)then
1562         if(jj.eq.1)then
1563           gb=psjti(q2,qq,s2new,iqnew,iqc(2),1)
1564         else
1565           gb=psjti1(q2,qmin(1),qq,s2new,iqnew,iqc(1),1)
1566         endif
1567         if(iqnew.eq.0)then
1568           gb=gb/sj1
1569         else
1570           gb=gb/sj2
1571         endif
1572         if(ish.ge.1)then
1573           if(gb.gt.1.)write (ifch,*)'gb,s2new,s2,q2,iqnew',
1574      *    gb,s2new,s2,q2,iqnew
1575         endif
1576         if(rangen().gt.gb)goto 9
1577       else
1578         goto 9
1579       endif
1580       jgamma=0
1581 
1582       call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
1583 
1584       if(jt.eq.1)then
1585         ncc(jq,jj)=nqc(2)
1586       elseif(jt.eq.2)then
1587         ncc(jq,jj)=ncc(1,jj)
1588         ncc(3-jq,jj)=nqc(1)
1589       elseif(jt.eq.3)then
1590         ncc(1,jj)=nqc(2)
1591       elseif(jt.eq.4)then
1592         ncc(1,jj)=ncc(3-jq,jj)
1593       elseif(jt.eq.5)then
1594         ncc(1,jj)=nqc(1)
1595         ncc(2,jj)=0
1596       endif
1597       iqc(jj)=iqnew
1598       if(ish.ge.6)write(ifch,*)'qt2,amprt,ncc:',
1599      *qt2,amprt,ncc(1,jj),ncc(2,jj)
1600 
1601       do i=1,4
1602         ept(i)=ept1(i)
1603       enddo
1604 c c.m. energy squared, minimal  4-momentum transfer square and gluon 4-v
1605 c for the next ladder run
1606       qmin(jj)=q2
1607       si=s2new
1608       if(ish.ge.3)write (ifch,*)'res: new jet - iqj,ncj,ep3,ept',
1609      *iqj(nj),ncj(1,nj),ncj(2,nj),ep3,ept
1610 
1611       goto 8            !next simulation step will be considered
1612 
1613 10    continue
1614       if(ish.ge.3)write(ifch,*)'res: iqc,si,ept:',iqc,si,ept
1615 
1616 c highest virtuality subprocess in the ladder
1617 c---------------------------------------
1618       qqs=max(qmin(1)/4.,4.*qmin(2))
1619       qqs=max(qqs,qq)
1620       call psabor(si,qqs,iqc,ncc,ept,1,nptlh,bx)
1621       goto 17
1622 
1623 12    continue
1624 c---------------------------------------
1625 c hard pomeron (direct photon)
1626 c---------------------------------------
1627       ept(1)=.5*(wpi+wmi)
1628       ept(2)=.5*(wpi-wmi)
1629       ept(3)=0.
1630       ept(4)=0.
1631       if(ish.ge.3)write (ifch,*)'direct photon - ept,si,qq:',ept,si,qq
1632 
1633 13    continue
1634 
1635 c ladder rung
1636 c---------------------------------------
1637       pt2=ept(3)**2+ept(4)**2
1638       pt=sqrt(pt2)
1639       wpt(1)=ept(1)+ept(2)
1640       wpt(2)=si/wpt(1)
1641 
1642       gdbor=psdbin(qmin(2),qq,si,iqc(2),long)
1643       gdtot=psdsin(qmin(2),qq,si,iqc(2),long)
1644       if(iqc(2).ne.0)then
1645         if(ish.ge.8)write (ifch,*)'qmin(2),qq,si',qmin(2),qq,si
1646         gdnon=psdnsi(qmin(2),qq,si,long)
1647         if(iabs(iqc(2)).eq.1.or.iabs(iqc(2)).eq.4)then
1648           gdbor=gdbor/2.25
1649           gdtot=gdnon/2.25+gdtot/4.5
1650         else
1651           gdbor=gdbor/9.
1652           gdtot=gdnon/9.+gdtot/4.5
1653         endif
1654       else
1655         gdnon=0.
1656       endif
1657 
1658       if(long.ne.0.or.qmin(2).ge.qq)then
1659         s2min=qq+4.*max(qmin(2),qcmass**2)
1660         wwmin=(5.*s2min-qq)/4.-2.*pt*sqrt(q2ini)
1661         wwmin=(wwmin+sqrt(wwmin**2-(qq-pt2)*(s2min-qq-4.*q2ini)))
1662      *  /2./(1.-4.*q2ini/(s2min-qq))
1663       else
1664         s2min=qq/(1.-sqrt(q2ini/qq))
1665         wwmin=s2min+qq-2.*pt*sqrt(q2ini)
1666         wwmin=(wwmin+sqrt(wwmin**2-4.*(qq-pt2)*(qq-q2ini)))
1667      *  /2./(1.-q2ini/qq)
1668       endif
1669 
1670       if(ish.ge.3)write(ifch,*)'si,s2min,wwmin,qmin(2),gdtot,gdbor:'
1671       if(ish.ge.3)write(ifch,*)si,s2min,wwmin,qmin(2),gdtot,gdbor
1672 
1673       if((rangen().lt.gdbor/gdtot.or.si.lt.1.1*wwmin).and.
1674      *(long.eq.0.and.qmin(2).lt.qq.or.iqc(2).eq.0))goto 15
1675       if(si.lt.1.1*wwmin)stop'si<1.1*wwmin'
1676 
1677       qqmax=0.
1678       qqmin=0.
1679 
1680       xmm=si+2.*sqrt(q2ini)*pt-qq
1681       discr=xmm**2-4.*q2ini*(5.*si-qq+pt2)
1682       if(discr.lt.0.)goto 29
1683       discr=sqrt(discr)
1684       qqmax=(xmm+discr)/2./(5.-(qq-pt2)/si)
1685       qqmin=2.*q2ini*si/(xmm+discr)
1686 
1687 29    continue
1688       if(4.*qqmin.lt.s2min-qq.or.long.eq.0.and.
1689      *qmin(2).lt.qq)then
1690         xmm=si-s2min+2.*sqrt(q2ini)*pt
1691         qqmin=2.*q2ini*si/(xmm+sqrt(xmm**2-4.*q2ini*(si-qq+pt2)))
1692       endif
1693       xmin=1.-q2ini/qqmin
1694 
1695       if(qqmin.lt.qmin(2))then
1696         qqmin=qmin(2)
1697         xmi=max(1.-((pt*sqrt(qqmin)+sqrt(pt2*qqmin+
1698      *  si*(si-s2min-qqmin*(1.-(qq-pt2)/si))))/si)**2,
1699      *  (s2min+qqmin*(1.-(qq-pt2)/si)-2.*pt*sqrt(qqmin))/si)
1700         xmin=max(xmin,xmi)
1701       endif
1702       if(xmin.le.qq/si)xmin=1.001*qq/si
1703 
1704       if(long.eq.0.and.qmin(2).lt.qq)qqmax=max(qqmax,qq)
1705       xmax=1.-q2ini/qqmax
1706 
1707       if(ish.ge.6)write(ifch,*)'qqmax,qqmin,xmax,xmin:',
1708      *qqmax,qqmin,xmax,xmin
1709       if(qqmax.lt.qqmin)stop'qqmax<qqmin'
1710 
1711       qm0=qqmin
1712       xm0=1.-q2ini/qm0
1713       s2max=si*xm0-qm0*(1.-qq/si)
1714 
1715       sds=psdsin(qm0,qq,s2max,0,long)/4.5
1716       sdv=psdsin(qm0,qq,s2max,1,long)/4.5
1717 
1718       sdn=psdnsi(qm0,qq,s2max,long)
1719       if(iqc(2).eq.0)then
1720         sdn=sdn/4.5
1721       elseif(iabs(iqc(2)).eq.1.or.iabs(iqc(2)).eq.4)then
1722         sdn=sdn/2.25
1723       else
1724         sdn=sdn/9.
1725       endif
1726       sdv=sdv+sdn
1727       xx=xm0
1728 
1729       sj0=sds*psfap(xx,iqc(2),0)+sdv*psfap(xx,iqc(2),1)
1730       gb0=sj0/log(q2ini/qcdlam)*sngl(psuds(qm0,iqc(2)))*qm0*5.
1731       if(gb0.le.0.)then
1732          write(ifmt,*)'gb0.le.0.  si,qq,pt2:',si,qq,pt2
1733          iret=1
1734          goto 9999
1735       endif
1736 
1737       if(xm0.le..5)then
1738         gb0=gb0*(xm0-qq/si)/(1.-2.*qq/si)
1739       else
1740         gb0=gb0*(1.-xm0)
1741       endif
1742 
1743       xmin2=max(.5,xmin)
1744       xmax1=min(xmax,.5)
1745       if(xmin.ge..5)then
1746         djl=1.
1747       elseif(xmax.lt..5)then
1748         djl=0.
1749       else
1750         djl=1./(1.-(1.-2.*qq/si)*log((.5-qq/si)/(xmin-qq/si))/
1751      *  log(2.*(1.-xmax)))
1752       endif
1753 
1754 14    continue
1755       if(rangen().gt.djl)then        !lc momentum share in the cur
1756         x=(xmin-qq/si)*((xmax1-qq/si)/(xmin-qq/si))**rangen()+qq/si
1757       else
1758         x=1.-(1.-xmin2)*((1.-xmax)/(1.-xmin2))**rangen()
1759       endif
1760 
1761       q2=qqmin/(1.+rangen()*(qqmin/qqmax-1.))
1762 
1763       qt2=q2*(1.-x)
1764       if(ish.ge.9)write(ifch,*)'q2,x,qt2,qq,qqmin,qqmax:',
1765      *q2,x,qt2,qq,qqmin,qqmax
1766       if(qt2.lt.q2ini)goto 14   !p_t check
1767 
1768       if(long.ne.0.or.q2.ge.qq)then
1769         s2min2=max(4.*q2+qq,s2min)
1770       else
1771         s2min2=s2min
1772       endif
1773       qt=sqrt(qt2)
1774       call pscs(bcos,bsin)
1775 c ep3 is now 4-vector for s-channel gluon produced in current ladder run
1776       ep3(3)=qt*bcos
1777       ep3(4)=qt*bsin
1778       ptnew=(ept(3)-ep3(3))**2+(ept(4)-ep3(4))**2
1779 
1780       s2=x*si-ptnew+pt2-q2*(x-(qq-pt2)/si)
1781       if(s2.lt.s2min2)goto 14   !check of the kinematics
1782       sds=psdsin(q2,qq,s2,0,long)/4.5
1783       sdv0=psdsin(q2,qq,s2,1,long)/4.5
1784       if(ish.ge.8)write (ifch,*)'q2,qq,s2',q2,qq,s2
1785       sdn0=psdnsi(q2,qq,s2,long)
1786 
1787       if(iqc(2).eq.0)then
1788         sdn=sdn0/4.5
1789       else
1790         if(iabs(iqc(2)).eq.1.or.iabs(iqc(2)).eq.4)then
1791           sdn=sdn0/2.25
1792         else
1793           sdn=sdn0/9.
1794         endif
1795       endif
1796       sdv=sdv0+sdn
1797 
1798       xx=x
1799       sj1=sds*psfap(xx,iqc(2),0)
1800       sj2=sdv*psfap(xx,iqc(2),1)
1801 
1802 c gb7 is the rejection function for x and q**2 simulation.
1803       gb7=(sj1+sj2)/log(qt2/qcdlam)*sngl(psuds(q2,iqc(2)))/gb0*q2
1804       if(x.le..5)then
1805         gb7=gb7*(x-qq/si)/(1.-2.*qq/si)
1806       else
1807         gb7=gb7*(1.-x)
1808       endif
1809 
1810       if(gb7.gt.1..and.ish.ge.1)write(ifmt,*)'gb7,q2,x,qt2,iqc(2),'
1811      * ,'gb0,sj1,sj2',gb7,q2,x,qt2,iqc(2),gb0,sj1,sj2
1812       if(ish.ge.3)write (ifch,*)'gb7,q2,x,qt2,iqc(2),gb0,sj1,sj2,long',
1813      * gb7,q2,x,qt2,iqc(2),gb0,sj1,sj2,long
1814       if(rangen().gt.gb7)goto 14
1815 
1816 
1817       if(ish.ge.6)write(ifch,*)'iqc,ncc:',iqc(2),ncc(1,2),ncc(2,2)
1818       iqcnew=iqc(2)
1819       nqc(2)=0         !emitted parton color connections
1820       if(rangen().lt.sj1/(sj1+sj2).or.(long.ne.0.or.q2.ge.qq).and.
1821      *s2.lt.1.5*s2min2)then
1822         if(iqc(2).eq.0)then
1823           jt=1
1824           jq=int(1.5+rangen())
1825           nqc(1)=ncc(jq,2)
1826         else
1827           jt=2
1828           if(iqc(2).gt.0)then
1829             jq=1
1830           else
1831             jq=2
1832           endif
1833           nqc(1)=0
1834         endif
1835         iq1=iqc(2)
1836         iqcnew=0
1837 
1838       else
1839         if(iqc(2).ne.0)then
1840           jt=3
1841           iq1=0
1842           if(iqc(2).gt.0)then
1843             jq=1
1844           else
1845             jq=2
1846           endif
1847           nqc(1)=ncc(1,2)
1848 
1849         else
1850           tu=sdn0/2.25+sdv0
1851           if(naflav.eq.4)tu=tu*2.
1852           td=sdn0/9.+sdv0
1853           if(rangen().lt.tu/(tu+2.*td))then
1854             if(naflav.eq.3)then
1855               iq1=1
1856             else
1857               iq1=1+3*int(.5+rangen())
1858             endif
1859           else
1860             iq1=int(2.5+rangen())
1861           endif
1862           jq=int(1.5+rangen())
1863           iq1=iq1*(3-2*jq)
1864           iqcnew=-iq1
1865           jt=4
1866           nqc(1)=ncc(jq,2)
1867         endif
1868       endif
1869 
1870       eprt=max(1.d0*qt,
1871      *.5d0*((1.d0-x)*wpt(2)+qt2/(1.d0-x)/wpt(2)))
1872       pl=eprt-(1.d0-x)*wpt(2)
1873       zeta=sqrt(qq/si)/sqrt(x*(1.-x))
1874       if(iq1.eq.0)then
1875         iq2ini=9
1876         jo=iq2ini
1877         if(zeta.gt.zetacut)jo=-jo
1878       else
1879         iq2ini=iq1
1880         jo=iq2ini
1881       endif
1882 28    call timsh1(q2,sngl(eprt),iq2ini,jo)
1883       amprt=pprt(5,1)**2
1884       plprt=eprt**2-amprt-qt2
1885       if(plprt.lt.-1d-6)goto 28
1886       ep3(1)=eprt
1887       ep3(2)=dsqrt(max(0.d0,plprt))
1888       if(pl.lt.0.d0)ep3(2)=-ep3(2)
1889       ey(1)=1.
1890       ey(2)=1.
1891       ey(3)=1.
1892       do i=1,4
1893         ept1(i)=ept(i)-ep3(i)
1894       enddo
1895       call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
1896       call psrotat(ep3,s0xh,c0xh,s0h,c0h)
1897       s2new=psnorm(ept1)+qq
1898 
1899       if((long.ne.0.or.q2.ge.qq).and.iqcnew.ne.0)then
1900         xmm=(5.*s2min2-qq)/4.-2.*sqrt(ptnew*q2ini)
1901         s2min2=1.1*(xmm+sqrt(xmm**2-(qq-ptnew)*
1902      *  (s2min2-qq-4.*q2ini)))/2./(1.-4.*q2ini/(s2min2-qq))
1903       endif
1904       if(s2new.gt.s2min2)then
1905         sds1=psdsin(q2,qq,s2new,iqcnew,long)/4.5
1906         if(iqcnew.eq.0)then
1907           gb=sds1/sds
1908         else
1909           if(ish.ge.8)write (ifch,*)'q2,qq,s2new',q2,qq,s2new
1910           sdn1=psdnsi(q2,qq,s2new,long)
1911           if(iabs(iqcnew).eq.1.or.iabs(iqcnew).eq.4)then
1912             sdn1=sdn1/2.25
1913             sdv=sdv0+sdn0/2.25
1914           else
1915             sdn1=sdn1/9.
1916             sdv=sdv0+sdn0/9.
1917           endif
1918           gb=.9999*(sds1+sdn1)/sdv
1919         endif
1920         if(ish.ge.3.and.gb.gt.1..and.ish.ge.1)write(ifmt,*)'gbs2',gb
1921         if(rangen().gt.gb)goto 14
1922       else
1923         goto 14
1924       endif
1925 
1926       call psreti(nqc,jq,1,ey,s0xh,c0xh,s0h,c0h)
1927 
1928       iqc(2)=iqcnew
1929       if(jt.eq.1)then      !current parton color connections
1930         ncc(jq,2)=nqc(2)
1931       elseif(jt.eq.2)then
1932         ncc(jq,2)=ncc(1,2)
1933         ncc(3-jq,2)=nqc(1)
1934       elseif(jt.eq.3)then
1935         ncc(1,2)=nqc(2)
1936       elseif(jt.eq.4)then
1937         ncc(1,2)=ncc(3-jq,2)
1938         ncc(2,2)=0
1939       endif
1940 
1941       do i=1,4
1942         ept(i)=ept1(i)
1943       enddo
1944       if(ish.ge.3)write (ifch,*)'new jet - iqj,ncj,ep3,ept',
1945      *iqj(nj),ncj(1,nj),ncj(2,nj),ep3,ept
1946 c c.m. energy squared, minimal  4-momentum transfer square and gluon 4-v
1947 c for the next ladder run
1948       qmin(2)=q2
1949       si=s2new
1950       goto 13            !next simulation step will be considered
1951 
1952 15    continue
1953       if(ish.ge.3)write(ifch,*)'iqc,si,qmin(2),nj:',
1954      *iqc(2),si,qmin(2),nj
1955 c highest virtuality subprocess in the ladder
1956 c---------------------------------------
1957       gb01=0.
1958       tmax=0.
1959       tmin=si
1960       if(iqc(2).eq.0.and.si.gt.qq+4.*max(qcmass**2,qmin(2)))then
1961         qminn=max(qcmass**2,qmin(2))
1962         tmin1=2.*qminn/(1.-qq/si)/(1.+sqrt(1.-4.*qminn/(si-qq)))
1963         tmin=tmin1
1964         tmax=si/2.
1965         fb01=psdbom(si,si/2.,si/2.,qq,long)
1966         if(long.eq.0)fb01=fb01*si/2.
1967         gb01=fb01/log(qminn/qcdlam)*sngl(psuds(qminn,iqc(2)))/si**2
1968         gb0=gb01
1969       else
1970         tmin1=0.
1971       endif
1972 
1973       if(long.eq.0.and.qmin(2).lt.qq)then
1974         tmax=max(tmax,qq)
1975         tmin=max(qmin(2),
1976      *  2.*q2ini/(1.-qq/si)/(1.+sqrt(1.-4.*q2ini/(si-qq))))
1977         ze=qq/si+tmin/si*(1.-qq/si)
1978         xx=ze
1979         qt2=tmin*(1.-ze)
1980         if(qt2.lt..999*q2ini.and.ish.ge.1)write(ifmt,*)'bor-dir:qt20'
1981      *                                                 ,qt2
1982         gb0=gb01+psfap(xx,iqc(2),1)/log(qt2/qcdlam)
1983      *  *sngl(psuds(tmin,iqc(2))/psuds(tmin,1)*psuds(qq,1))
1984      *  /si*(1.-tmin*qq/si**2/ze)
1985       endif
1986       gb0=gb0*2.
1987 
1988       call psdeftr(si-qq,ept,ey)
1989 
1990       if(ish.ge.6)write(ifch,*)'tmin,tmax,qq,si-qq,gb0:'
1991       if(ish.ge.6)write(ifch,*)tmin,tmax,qq,si-qq,psnorm(ept),gb0
1992 
1993 c------------------------------------------------
1994 16    continue
1995       if(long.eq.0)then
1996         t=tmin*(tmax/tmin)**rangen()
1997       else
1998         t=tmin+(tmax-tmin)*rangen()
1999       endif
2000 
2001       u=si-t
2002       ze=qq/si+t/si*(1.-qq/si)
2003       qt2=t*(1.-ze)
2004       if(t.le.qq.and.long.eq.0)then
2005         xx=ze
2006         gb=psfap(xx,iqc(2),1)/log(qt2/qcdlam)*sngl(psuds(t,iqc(2))
2007      *  /psuds(t,1)*psuds(qq,1))/si*(1.-t*qq/si**2/ze)/gb0
2008       else
2009         gb=0.
2010       endif
2011 
2012       gb1=0.
2013       if(iqc(2).eq.0..and.si.gt.qq+4.*max(qcmass**2,qmin(2)).
2014      *and.qt2.gt.qcmass**2.and.t.le.si/2..and.t.ge.tmin1)then
2015         fb1=psdbom(si,t,u,qq,long)
2016         if(long.eq.0)fb1=fb1*t
2017         gb1=fb1/log(qt2/qcdlam)*sngl(psuds(qt2,iqc(2)))/si**2/gb0
2018 c        gb1=0.  !???????????????????????
2019         gb=gb+gb1
2020       endif
2021 
2022 
2023       if(ish.ge.6)write(ifch,*)'gb,t,iqc(2),si,qq,qmin(2),long:',
2024      *gb,t,iqc(2),si,qq,qmin(2),long
2025       if (ish.ge.1) then
2026         if(gb.gt.1.)write(*,*)'gb,gb1,gb0,gb01',
2027      *  ',t,iqc(2),si,qq,qmin(2),long:',
2028      *  gb,gb1,gb0,gb01,fb1,fb01,t,iqc(2),si,qq,qmin(2),long
2029       endif
2030 
2031       if(rangen().gt.gb)goto 16
2032       if(ish.ge.3)write(ifch,*)'born:t,qt2:',t,qt2
2033 
2034       nqc(2)=0
2035       if(iqc(2).eq.0)then
2036         jq=int(1.5+rangen())
2037         jq2=3-jq
2038         if(rangen().gt.gb1/gb)then
2039           iq1=(1+int(3.*rangen()))*(3-2*jq)
2040         else
2041           iq1=4*(3-2*jq)
2042         endif
2043         iq2=-iq1                             !quark flavors
2044         nqc(1)=ncc(jq,2)
2045       else
2046         if(iqc(2).gt.0)then
2047           jq=1
2048         else
2049           jq=2
2050         endif
2051         jq2=jq
2052         iq1=0
2053         iq2=iqc(2)
2054         nqc(1)=ncc(1,2)
2055       endif
2056 
2057       call pscs(bcos,bsin)
2058       z=sngl(psutz(dble(si-qq),dble(qt2),dble(qt2)))
2059       if(t.lt..5*si)z=1.-z
2060       wp3=z*sqrt(si-qq)
2061       wm3=qt2/wp3
2062       if(iabs(iq1).eq.4)qt2=qt2-qcmass**2
2063       qt=sqrt(qt2)
2064       ep3(1)=.5*(wp3+wm3)
2065       ep3(2)=.5*(wp3-wm3)
2066       ep3(3)=qt*bcos
2067       ep3(4)=qt*bsin
2068       call psdefrot(ep3,s0xh,c0xh,s0h,c0h)
2069       zeta=2.
2070       if(iq1.eq.0)then
2071         iq2ini1=9
2072         jo1=iq2ini1
2073         if(zeta.gt.zetacut)jo1=-jo1
2074       else
2075         iq2ini1=iq1
2076         jo1=iq2ini1
2077       endif
2078       if(iq2.eq.0)then
2079         iq2ini2=9
2080         jo2=iq2ini2
2081         if(zeta.gt.zetacut)jo2=-jo2
2082       else
2083         iq2ini2=iq2
2084         jo2=iq2ini2
2085       endif
2086       if(ish.ge.5)write (ifch,*)'jq,jt,iq2ini1,iq2ini2',
2087      *jq,jt,iq2ini1,iq2ini2
2088 
2089       if(t.lt.qq.and.iabs(iq1).ne.4)then
2090         qq1=t*(1.-ze)
2091         qq2=qq
2092       else
2093         qq1=qt2
2094         qq2=qt2
2095       endif
2096       call timsh2(qq1,qq2,sqrt(si-qq),iq2ini1,iq2ini2,jo1,jo2)
2097       nfprt=1
2098       call psreti(nqc,jq,nfprt,ey,s0xh,c0xh,s0h,c0h)
2099 
2100       if(iqc(2).eq.0)then
2101         nqc(1)=ncc(3-jq,2)
2102         nqc(2)=0
2103       else
2104         nqc(1)=nqc(2)
2105         nqc(2)=0
2106       endif
2107 
2108       nfprt=2
2109       call psreti(nqc,jq2,nfprt,ey,s0xh,c0xh,s0h,c0h)
2110 
2111 17    continue
2112       if(ish.ge.3)write (ifch,*)'nj',nj
2113       if(nj.gt.0)then
2114 
2115           ityj(i)=30
2116           iorj(i)=nptlh
2117         do n=1,nj
2118           do i=1,4
2119             ep3(i)=eqj(i,n)
2120           enddo
2121           call pstrans(ep3,ey0,-1)         !boost to the c.m. system
2122           do i=1,4
2123             eqj(i,n)=ep3(i)
2124           enddo
2125           do l=1,6
2126             bxj(l,n)=bx(l)
2127           enddo
2128           ityj(n)=0
2129           iorj(n)=1
2130         enddo
2131       endif
2132       call psjarr(jfl)       !kinky strings formation
2133       if(ish.ge.3)write (ifch,*)'jfl',jfl
2134       if(jfl.eq.0)then
2135         iret=1
2136       else
2137         iret=0
2138         ep3(4)=egyevt
2139         ep3(2)=0.
2140         ep3(3)=0.
2141         ep3(1)=0.
2142         do i=2,nptl
2143           do l=1,4
2144             ep3(l)=ep3(l)-pptl(l,i)
2145           enddo
2146         enddo
2147         if(ish.ge.3)write(ifch,*)'energy-momentum balance:'
2148         if(ish.ge.3)write(ifch,*)ep3
2149         if(abs(ep3(4)).gt.3.e-2)write(*,*)'energy-momentum balance:',ep3
2150       endif
2151  9999 call utprix('psadis',ish,ishini,3)
2152       return
2153       end
2154 
2155 c-----------------------------------------------------------------------
2156       subroutine psaevc
2157 c-----------------------------------------------------------------------
2158       include 'epos.inc'
2159 c structure functions calculation
2160       logical lcalc
2161       double precision xx,xmax
2162       dimension evs(21,21,135)
2163       common /psar2/  edmax,epmax
2164       common /psar31/ evk0(21,21,54)
2165       common /psar32/ evk(21,21,135)
2166       common/producetab/ producetables              !used to link with CRMC
2167       logical producetables
2168       include 'epos.incsem'
2169 
2170       inquire(file=fnie(1:nfnie),exist=lcalc)
2171       if(lcalc)then
2172        if(inicnt.eq.1)then
2173         write(ifmt,'(3a)')'read from ',fnie(1:nfnie),' ...'
2174         open(1,file=fnie(1:nfnie),status='old')
2175         read (1,*)qcdlam0,q2min0,q2ini0,naflav0,epmax0
2176         if(qcdlam0.ne.qcdlam)write(ifmt,'(a)')'iniev: wrong qcdlam'
2177         if(q2min0 .ne.q2min )write(ifmt,'(a)')'iniev: wrong q2min'
2178         if(q2ini0 .ne.q2ini )write(ifmt,'(a)')'iniev: wrong q2ini'
2179         if(naflav0.ne.naflav)write(ifmt,'(a)')'iniev: wrong naflav'
2180         if(epmax0 .ne.epmax )write(ifmt,'(a)')'iniev: wrong epmax'
2181         if(qcdlam0.ne.qcdlam.or.q2min0.ne.q2min.or.q2ini0.ne.q2ini
2182      *  .or.naflav0.ne.naflav.or.epmax0.ne.epmax)then
2183            write(6,'(//a//)')'   iniev has to be reinitialized!!!'
2184            stop
2185         endif
2186         read (1,*)evk0,evk
2187         close(1)
2188        endif
2189        goto 101
2190 
2191       elseif(.not.producetables)then
2192         write(ifmt,*) "Missing epos.iniev file !"        
2193         write(ifmt,*) "Please correct the defined path ",
2194      &"or force production ..."
2195         stop
2196 
2197       endif
2198 
2199       write(ifmt,'(a)')'iniev does not exist -> calculate tables  ...'
2200       xmax=1.d0-2.d0*q2ini/epmax
2201       do l=1,27
2202         if(l.le.12)then
2203           xx=.1d0*exp(l-13.d0)
2204         elseif(l.le.21)then
2205           xx=.1d0*(l-12.d0)
2206         else
2207           xx=1.d0-.1d0*(10.d0*(1.d0-xmax))**((l-21)/6.)
2208         endif
2209 
2210         qmin=max(1.d0*q2min,q2ini/(1.-xx))
2211       do i=1,21
2212         qq=qmin*(.5*epmax/qmin)**((i-1)/20.)
2213       do j=1,21
2214         qj=qmin*(qq/qmin)**((j-1)/20.)
2215         if(l.eq.27.or.i.eq.1.or.j.eq.21)then
2216           evk0(i,j,l)=0.
2217           evk0(i,j,l+27)=0.
2218           do k=1,5
2219             evk(i,j,l+27*(k-1))=0.
2220           enddo
2221         else
2222           do k=1,2
2223             evk0(i,j,l+27*(k-1))=log(psev0(qj,qq,xx,k))
2224           enddo
2225         endif
2226       enddo
2227       enddo
2228       enddo
2229 
2230       n=1
2231 
2232 1     n=n+1
2233       write(ifmt,2)n
2234 2     format(5x,i2,'-th order contribution')
2235 
2236       do l=1,26
2237         write(ifmt,*)'l',l
2238         if(l.le.12)then
2239           xx=.1d0*exp(l-13.d0)
2240         elseif(l.le.21)then
2241           xx=.1d0*(l-12.d0)
2242         else
2243           xx=1.d0-.1d0*(10.d0*(1.d0-xmax))**((l-21)/6.)
2244         endif
2245 
2246         qmin=max(1.d0*q2min,q2ini/(1.d0-xx))
2247       do i=2,21
2248         qq=qmin*(.5*epmax/qmin)**((i-1)/20.)
2249       do j=1,20
2250         qj=qmin*(qq/qmin)**((j-1)/20.)
2251         do m=1,3
2252         do k=1,2
2253           if(m.ne.3)then
2254             ev=psev(qj,qq,xx,m,k,n)
2255             ev0=psevi0(qj,qq,xx,m,k)
2256             evs(i,j,l+27*(m-1)+54*(k-1))=log((ev+ev0)/psfap(xx,m-1,k-1)
2257      *      /log(log(qq*(1.d0-xx)/qcdlam)/log(qj*(1.d0-xx)/qcdlam))*4.5)
2258           elseif(k.ne.1)then
2259             evs(i,j,l+108)=log((psev(qj,qq,xx,m,k,n)+
2260      *      psevi0(qj,qq,xx,2,2))/psfap(xx,2,2)
2261      *      /log(log(qq*(1.d0-xx)/qcdlam)/log(qj*(1.d0-xx)/qcdlam))*4.5)
2262           endif
2263         enddo
2264         enddo
2265       enddo
2266       enddo
2267       enddo
2268 
2269       jec=0
2270       do i=2,21
2271       do j=1,20
2272       do l=1,26
2273       do k=1,5
2274         if(n.eq.2.or.evs(i,j,l+27*(k-1)).ne.0..and.
2275      *  abs(1.-evk(i,j,l+27*(k-1))/evs(i,j,l+27*(k-1))).gt.1.e-2)then
2276           jec=1
2277           evk(i,j,l+27*(k-1))=evs(i,j,l+27*(k-1))
2278         endif
2279       enddo
2280       enddo
2281       enddo
2282       enddo
2283 
2284       if(jec.ne.0)goto 1
2285 
2286       write(ifmt,'(a)')'write to iniev ...'
2287       open(1,file=fnie(1:nfnie),status='unknown')
2288       write (1,*)qcdlam,q2min,q2ini,naflav,epmax
2289       write (1,*)evk0,evk
2290       close(1)
2291 
2292 101   continue
2293       return
2294       end
2295 
2296 c------------------------------------------------------------------------
2297       function psdbom(s,t,u,qq,long)
2298 c-----------------------------------------------------------------------
2299 c psdbom - integrand for DIS c-quark cross-sections (matrix element squared)
2300 c s  - total c.m. energy squared for the scattering (for n=2: s+qq),
2301 c t  - invariant variable for the scattering |(p1-p3)**2|
2302 c u  - invariant variable for the scattering |(p1-p4)**2|
2303 c qq - photon virtuality
2304 c long: 0 - contr. to (F2-F_L), 1 - contr. to F_L
2305 c-----------------------------------------------------------------------
2306       include 'epos.incsem'
2307       if(long.eq.0)then       !F2-F_L
2308         psdbom=(2.*(t/u+u/t)*(qq**2+(s-qq)**2)/s**2+
2309      *  4.*(qcmass*s/t/u)**2*(qq-2.*qcmass**2)+
2310      *  8.*qcmass**2/t/u*(s-2.*qq))   *2.    !=4.5/2.25
2311       else                    !F_L_C
2312         psdbom=16.*qq*((s-qq)/s**2-qcmass**2/t/u) *2.  !=4.5/2.25
2313       endif
2314       return
2315       end
2316 
2317 c------------------------------------------------------------------------
2318       function psdbin(q1,qq,s,m1,long)
2319 c-----------------------------------------------------------------------
2320 c psdbin - DIS born cross-section
2321 c q1      - virtuality cutoff for current end of the ladder
2322 c qq      - photon virtuality
2323 c s=2(pq) - s_true + qq
2324 c s2min   - mass cutoff for born scattering
2325 c m1       - incoming parton type (0 - g, 1,2 - q)
2326 c-----------------------------------------------------------------------
2327       double precision xx
2328       include 'epos.incsem'
2329       include 'epos.inc'
2330 
2331       psdbin=0.
2332       q2mass=qcmass**2
2333       s2min=4.*max(q1,q2mass)+qq
2334       if(m1.eq.0.and.s.gt.s2min.and.(idisco.eq.0.or.idisco.eq.2))then
2335         tmax=s/2.
2336         qtq=4.*max(q2mass,q1)/(s-qq)
2337         if(qtq.lt.1.)then
2338           tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
2339         else
2340           tmin=.5*s
2341         endif
2342         psdbin=psdbin+psdbor(q1,qq,s,long)*(1./tmin-1./tmax)
2343       endif
2344 
2345       if(long.eq.0.and.q1.lt.qq.and.s.gt.qq/(1.-q2ini/qq)
2346      *.and.(idisco.eq.0.or.idisco.eq.1))then
2347         m=min(1,iabs(m1))+1
2348         xx=qq/s
2349         psdbin=psdbin+psevi0(q1,qq,xx,m,2)*4.*pi**2*alfe/s
2350       endif
2351       return
2352       end
2353 
2354 c------------------------------------------------------------------------
2355       function psdbor(q1,qq,s,long)
2356 c-----------------------------------------------------------------------
2357 c psdbor - DIS born cross-section
2358 c q1      - virtuality cutoff for current end of the ladder
2359 c qq      - photon virtuality
2360 c s=2(pq) - s_true + qq
2361 c s2min   - mass cutoff for born scattering
2362 c-----------------------------------------------------------------------
2363       common /ar3/    x1(7),a1(7)
2364       include 'epos.inc'
2365       include 'epos.incsem'
2366       double precision psuds
2367 
2368       psdbor=0.
2369       q2mass=qcmass**2
2370       qtq=4.*max(q2mass,q1)/(s-qq)
2371       j=0   !Gluon
2372 
2373       tmax=s/2.
2374       if(qtq.lt.1.)then
2375         tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
2376       else
2377         tmin=.5*s
2378       endif
2379       if(tmax.lt.tmin.and.ish.ge.1)write(ifmt,*)'s,q1,qq,tmin,tmax',
2380      *s,q1,qq,tmin,tmax
2381 
2382       ft=0.
2383       do i=1,7
2384       do m=1,2
2385         t=2.*tmin/(1.+tmin/tmax+(2*m-3)*x1(i)*(1.-tmin/tmax))
2386         u=s-t
2387 
2388         qt=t*u/s*(1.-qq/s)
2389         if(qt.lt..999*max(q2mass,q1).and.ish.ge.1)
2390      &  write(ifmt,*)'psdbor:qt,q1',qt,q1
2391         fb=psdbom(s,t,u,qq,long)*t**2
2392         ft=ft+a1(i)*fb*pssalf(qt/qcdlam)*sngl(psuds(qt,j))
2393       enddo
2394       enddo
2395       psdbor=ft/s**2*pi**2*alfe/sngl(psuds(q1,j))
2396       return
2397       end
2398 
2399 c------------------------------------------------------------------------
2400       subroutine psdint(s,qq,sds,sdn,sdb,sdt,sdr,m1,long)
2401 c-----------------------------------------------------------------------
2402 c psdint - dis cross-sections interpolation - for minimal
2403 c effective momentum cutoff in the ladder
2404 c s   - total c.m. energy squared for the ladder,
2405 c qq  - photon virtuality,
2406 c sds - dis singlet cross-section,
2407 c sdn - dis nonsinglet cross-section,
2408 c sdb - dis born cross-section,
2409 c sdt - dis singlet+resolved cross-section,
2410 c m1  - parton type at current end of the ladder (0 - g, 1,2 - q)
2411 c-----------------------------------------------------------------------
2412       double precision xx
2413       dimension wk(3),wj(3)
2414       common /psar2/  edmax,epmax
2415       common /psar27/ csds(21,26,4),csdt(21,26,2),csdr(21,26,2)
2416       include 'epos.incsem'
2417       include 'epos.inc'
2418 
2419       sds=0.
2420       sdn=0.
2421       sdt=0.
2422       sdr=0.
2423       sdb=psdbin(q2min,qq,s,m1,long)
2424 
2425       m=min(1,iabs(m1))+1
2426       qlj=log(qq/q2min)*2.+1.
2427       j=int(qlj)
2428       if(j.lt.1)j=1
2429       if(j.gt.19)j=19
2430       wj(2)=qlj-j
2431       wj(3)=wj(2)*(wj(2)-1.)*.5
2432       wj(1)=1.-wj(2)+wj(3)
2433       wj(2)=wj(2)-2.*wj(3)
2434 
2435       s2min=4.*max(q2min,qcmass**2)+qq
2436       if(m1.ne.0)s2min=s2min/(1.-4.*q2ini/(s2min-qq))
2437       if(s.le.s2min.or.idisco.ne.0.and.idisco.ne.2)goto 1
2438 
2439       qtq=4.*max(q2min,qcmass**2)/(s-qq)
2440       if(qtq.lt.1.)then
2441         tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
2442       else
2443         tmin=.5*s
2444       endif
2445       tmax=s/2.
2446 
2447       sl=log(s/s2min)/log(edmax/s2min)*25.+1.
2448       k=int(sl)
2449       if(k.lt.1)k=1
2450       if(k.gt.24)k=24
2451       wk(2)=sl-k
2452       wk(3)=wk(2)*(wk(2)-1.)*.5
2453       wk(1)=1.-wk(2)+wk(3)
2454       wk(2)=wk(2)-2.*wk(3)
2455 
2456       do k1=1,3
2457         k2=k+k1-1
2458       do j1=1,3
2459         sds=sds+csds(j+j1-1,k2,m+2*long)*wj(j1)*wk(k1)
2460       enddo
2461       enddo
2462       if(m.eq.1)then
2463         sds=exp(sds)*(1./tmin-1./tmax)
2464       else
2465         sds=max(sds,0.)
2466       endif
2467 
2468 1     continue
2469       s2min=max(4.*qq,16.*q2min)+qq
2470       if(s.le.s2min.or.long.ne.0.or.idisco.ne.0.and.idisco.ne.3)then
2471         sdt=sds
2472         goto 2
2473       endif
2474 
2475       sl=log(s/s2min)/log(edmax/s2min)*25.+1.
2476       k=int(sl)
2477       if(k.lt.1)k=1
2478       if(k.gt.24)k=24
2479       wk(2)=sl-k
2480       wk(3)=wk(2)*(wk(2)-1.)*.5
2481       wk(1)=1.-wk(2)+wk(3)
2482       wk(2)=wk(2)-2.*wk(3)
2483 
2484       do k1=1,3
2485         k2=k+k1-1
2486       do j1=1,3
2487         sdr=sdr+csdr(j+j1-1,k2,m)*wj(j1)*wk(k1)
2488         sdt=sdt+csdt(j+j1-1,k2,m)*wj(j1)*wk(k1)
2489       enddo
2490       enddo
2491 
2492       sdr=max(sdr,0.)
2493       sdt=max(sds,sds+sdt)
2494       sdt=sdt+sdr
2495 
2496 2     continue
2497       if(long.eq.0.and.q2min.lt.qq.and.s.gt.qq/(1.-q2ini/qq)
2498      *.and.(idisco.eq.0.or.idisco.eq.1))then
2499         xx=qq/s
2500         dsi=psevi(q2min,qq,xx,m,2)*4.*pi**2*alfe/s
2501         if(m1.eq.0)then
2502           sds=sds+dsi
2503           sdt=sdt+dsi
2504         else
2505           dnsi=psevi(q2min,qq,xx,3,2)*4.*pi**2*alfe/s
2506           sdn=sdn+dnsi
2507           sds=sds+max(dsi-dnsi,0.)
2508           sdt=sdt+max(dsi-dnsi,0.)
2509         endif
2510       endif
2511 
2512       if(m1.eq.0)then
2513         sds=max(sds,sdb)
2514         sdt=max(sdt,sdb)
2515       else
2516         sdn=max(sdn,sdb)
2517       endif
2518       return
2519       end
2520 
2521 c-----------------------------------------------------------------------
2522       function psdnsi(q1,qq,s,long)
2523 c-----------------------------------------------------------------------
2524 c psdnsi - DIS nonsinglet cross-section interpolation
2525 c q1 - effective momentum cutoff for current end of the ladder,
2526 c qq - photon virtuality,
2527 c s - total c.m. energy squared for the ladder,
2528 c-----------------------------------------------------------------------
2529       double precision xx
2530       include 'epos.incsem'
2531       include 'epos.inc'
2532 
2533       psdnsi=0.
2534       if(long.eq.0.and.q1.lt.qq.and.s.gt.qq/(1.-q2ini/qq))then
2535         xx=qq/s
2536         psdnsi=psdnsi+max(0.,psevi(q1,qq,xx,3,2)*4.*pi**2*alfe/s)
2537       endif
2538       return
2539       end
2540 
2541 c-----------------------------------------------------------------------
2542       function psdrga(qq,s,s2min,j)
2543 c-----------------------------------------------------------------------
2544 c psdrga - DIS resolved cross-section (photon sf)
2545 c qq    - photon virtuality
2546 c s     - total c.m. energy squared for the process
2547 c s2min - mass cutoff for born scattering
2548 c j     - parton type at current end of the ladder (0 - g, 1,2 etc. - q)
2549 c-----------------------------------------------------------------------
2550       common /ar3/   x1(7),a1(7)
2551       include 'epos.incsem'
2552 
2553       psdrga=0.
2554       if(s.le.s2min)return
2555 
2556       xmin=s2min/s
2557       do i=1,7
2558       do m=1,2
2559         z=xmin**(.5+(m-1.5)*x1(i))
2560         tu=psdfh4(z,q2min,qq,0,1)
2561         td=psdfh4(z,q2min,qq,0,2)
2562         ts=psdfh4(z,q2min,qq,0,3)
2563         tg=psdfh4(z,q2min,qq,0,0)
2564         if(j.eq.0)then
2565           sj=tg*psjti(q2min,qq,z*s,0,j,1)+
2566      *    (tu+td+ts)*psjti(q2min,qq,z*s,1,j,1)
2567         else
2568           sj=tg*psjti(q2min,qq,z*s,0,j,1)+
2569      *    (tu+td)*(psjti(q2min,qq,z*s,1,1,1)/4.+
2570      *    psjti(q2min,qq,z*s,-1,1,1)/4.+
2571      *    psjti(q2min,qq,z*s,2,1,1)/2.)+
2572      *    ts*psjti(q2min,qq,z*s,2,1,1)
2573         endif
2574         psdrga=psdrga+a1(i)*sj
2575       enddo
2576       enddo
2577       psdrga=-psdrga*log(xmin)*alfe/2.  *4.5 !mean e^2 is taken out
2578       return
2579       end
2580 
2581 c-----------------------------------------------------------------------
2582       function psdres(qq,s,s2min,j)
2583 c-----------------------------------------------------------------------
2584 c psdres - DIS resolved photon cross-section
2585 c qq    - photon virtuality
2586 c s     - total w squared for the ladder (s+qq)
2587 c s2min - mass cutoff for born scattering
2588 c j     - parton type at current end of the ladder (0 - g, 1,2 etc. - q)
2589 c-----------------------------------------------------------------------
2590       double precision xx
2591       common /ar3/   x1(7),a1(7)
2592       include 'epos.inc'
2593       include 'epos.incsem'
2594 
2595       psdres=0.
2596       if(s.le.s2min+qq)return
2597 
2598       qmin=max(q2min,s2min/(s/qq-1.))
2599       qmax=min(s-s2min,s/2.)
2600 
2601 c numerical integration over transverse momentum squared;
2602 c gaussian integration is used
2603       do i=1,7
2604       do m=1,2
2605         qi=2.*qmin/(1.+qmin/qmax+(2*m-3)*x1(i)*(1.-qmin/qmax))
2606 
2607         zmax=min(1.,qi/qq)
2608         zmin=(max(qi,s2min)+qi)/s
2609 
2610         fsj=0.
2611         if(zmax.gt.zmin)then
2612           do i1=1,7
2613           do m1=1,2
2614             z=.5*(zmax+zmin+(2*m1-3)*x1(i1)*(zmax-zmin))
2615             s2=z*s-qi
2616             xx=z
2617             if(j.eq.0)then
2618               sj=psfap(xx,0,1)*psjti(qi,qq,s2,1,j,1)
2619             else
2620               sj=psfap(xx,0,1)*(psjti(qi,qq,s2,1,1,1)/6.+
2621      *        psjti(qi,qq,s2,-1,1,1)/6.+
2622      *        psjti(qi,qq,s2,2,1,1)/1.5)
2623             endif
2624             fsj=fsj+a1(i1)*sj*qi    !????????(qi-z*qq)
2625           enddo
2626           enddo
2627           fsj=fsj*(zmax-zmin)
2628         elseif(ish.ge.1)then
2629           write(ifmt,*)'psdres:zmax,zmin',zmax,zmin
2630         endif
2631         psdres=psdres+a1(i)*fsj
2632       enddo
2633       enddo
2634       psdres=psdres*(1./qmin-1./qmax)*alfe*.75/pi  !alpha_s -> 6 alpha_e
2635       return
2636       end
2637 
2638 c------------------------------------------------------------------------
2639       function psds(q1,qq,s,j,long)
2640 c-----------------------------------------------------------------------
2641 c psds - DIS singlet cross-section
2642 c q1      - virtuality cutoff for current end of the ladder
2643 c qq      - photon virtuality
2644 c s=2(pq) - s_true + qq
2645 c s2min   - mass cutoff for born scattering
2646 c-----------------------------------------------------------------------
2647       double precision xxe,xmax,xmin,xmax1,xmin1
2648       common /ar3/    x1(7),a1(7)
2649       include 'epos.inc'
2650       include 'epos.incsem'
2651 
2652       psds=0.
2653       q2mass=qcmass**2
2654       s2min=4.*max(q1,q2mass)
2655       smin=(s2min+qq)/(1.-4.*q2ini/s2min)
2656       if(s.le.1.001*smin)return
2657 
2658       xmax=.5d0*(1.d0+qq/s+dsqrt((1.d0-qq/s)**2-16.d0*q2ini/s))
2659       xmin=max(1.d0+qq/s-xmax,1.d0*(s2min+qq)/s)
2660       if(xmin.gt.xmax.and.ish.ge.1)write(ifmt,*)'xmin,xmax,q1,qq,s,smin'
2661      *,xmin,xmax,q1,qq,s,smin
2662 
2663       fx1=0.
2664       fx2=0.
2665       if(xmax.gt..9d0)then
2666         xmin1=max(xmin,.9d0)
2667         do i=1,7
2668         do m=1,2
2669           xxe=1.d0-(1.d0-xmax)*((1.d0-xmin1)/(1.d0-xmax))**
2670      *    (.5d0-x1(i)*(m-1.5))
2671           xx=xxe
2672 
2673           sh=xx*s
2674           qtmin=max(1.d0*max(q2mass,q1),q2ini/(1.d0-xxe))
2675           qtq=4.*qtmin/(sh-qq)
2676 
2677           tmin=.5*sh*qtq/(1.+sqrt(1.-qtq))
2678           tmax=.5*sh
2679           if(tmin.gt.tmax.and.ish.ge.1)write(ifmt,*)'psds:tmin,tmax'
2680      &                                              ,tmin,tmax
2681 
2682           ft=0.
2683           do i1=1,7
2684           do m1=1,2
2685             t=.5*(tmin+tmax+(2*m1-3)*x1(i1)*(tmin-tmax))
2686             u=sh-t
2687             qt=t*u/sh*(1.-qq/sh)
2688             if(qt.lt.qtmin.and.ish.ge.1)write(ifmt,*)'psds:qt,qtmin'
2689      &                                               ,qt,qtmin
2690 
2691             fb=psdsj(q1,xxe,sh,qt,t,u,qq,j,long)
2692             ft=ft+a1(i1)*fb*pssalf(qt/qcdlam)
2693           enddo
2694           enddo
2695           ft=ft*(tmax-tmin)
2696           fx1=fx1+a1(i)*ft*(1.-xx)/sh**2
2697         enddo
2698         enddo
2699         fx1=fx1*log((1.d0-xmin1)/(1.d0-xmax))
2700       endif
2701 
2702       if(xmin.lt..9d0)then
2703         xmax1=min(xmax,.9d0)
2704         do i=1,7
2705         do m=1,2
2706           xxe=xmin*(xmax1/xmin)**(.5-x1(i)*(m-1.5))
2707           xx=xxe
2708 
2709           sh=xx*s
2710           qtmin=max(1.d0*max(q2mass,q1),q2ini/(1.d0-xxe))
2711           qtq=4.*qtmin/(sh-qq)
2712 
2713           tmin=.5*sh*qtq/(1.+sqrt(1.-qtq))
2714           tmax=.5*sh
2715           if(tmin.gt.tmax.and.ish.ge.1)write(ifmt,*)'psds:tmin,tmax'
2716      *                                              ,tmin,tmax
2717 
2718           ft=0.
2719           do i1=1,7
2720           do m1=1,2
2721             t=(.5*(tmin+tmax+(2*m1-3)*x1(i1)*
2722      *      (tmin-tmax)))
2723             u=sh-t
2724             qt=t*u/sh*(1.-qq/sh)
2725             if(qt.lt.qtmin.and.ish.ge.1)write(ifmt,*)'psds:qt,qtmin'
2726      *                                               ,qt,qtmin
2727 
2728             fb=psdsj(q1,xxe,sh,qt,t,u,qq,j,long)
2729             ft=ft+a1(i1)*fb*pssalf(qt/qcdlam)
2730           enddo
2731           enddo
2732           ft=ft*(tmax-tmin)
2733           fx2=fx2+a1(i)*ft*xx/sh**2
2734         enddo
2735         enddo
2736         fx2=fx2*log(xmax1/xmin)
2737       endif
2738       psds=(fx1+fx2)*pi**2*alfe
2739       return
2740       end
2741 
2742 c-----------------------------------------------------------------------
2743       function psdsj(q1,xx,s,qt,t,u,qq,j,long)
2744 c-----------------------------------------------------------------------
2745 c psdsj - integrand for dis singlet cross-section
2746 c q1 - virtuality cutoff for current end of the ladder
2747 c xx - lc momentum ratio between initial (j) and final (l) partons
2748 c s  - c.m. energy squared for the born scattering,
2749 c t  - invariant variable for the born scattering |(p1-p3)**2|
2750 c u  - invariant variable for the born scattering |(p1-p4)**2|
2751 c qq - photon virtuality
2752 c j  - initial parton at the end of the ladder (0 - g, 1,2 - q)
2753 c-----------------------------------------------------------------------
2754       double precision xx
2755       include 'epos.incsem'
2756 
2757       fb=psdbom(s,t,u,qq,long)
2758       psdsj=psevi(q1,qt,xx,min(1,iabs(j))+1,1)*fb
2759       return
2760       end
2761 
2762 c------------------------------------------------------------------------
2763       function psdsin(q1,qq,s,m1,long)
2764 c-----------------------------------------------------------------------
2765 c psdsin - DIS singlet cross-section interpolation
2766 c q1 - effective momentum cutoff for current end of the ladder,
2767 c qq - photon virtuality,
2768 c s -  total c.m. energy squared for the ladder,
2769 c m1 - parton type at current end of the ladder (0 - g, 1,2 - q)
2770 c-----------------------------------------------------------------------
2771       double precision xx
2772       dimension wi(3),wj(3),wk(3)
2773       common /psar2/  edmax,epmax
2774       common /psar25/ csdsi(21,21,104)
2775       include 'epos.incsem'
2776       include 'epos.inc'
2777 
2778       psdsin=0.
2779       m=min(1,iabs(m1))+1
2780 
2781       q2mass=qcmass**2
2782       s2min=4.*max(q2min,q2mass)+qq
2783       sdmin=4.*max(q1,q2mass)+qq
2784       if(m1.ne.0)then
2785         s2min=s2min/(1.-4.*q2ini/(s2min-qq))
2786         sdmin=sdmin/(1.-4.*q2ini/(sdmin-qq))
2787       endif
2788 c      if(s.le.1.e8*sdmin)goto 2  !????????????????
2789       if(s.le.sdmin)goto 2
2790 
2791       qmin=q2min
2792       qmax=(s-qq)/4.
2793       if(m1.ne.0)qmax=(s-qq+sqrt((s-qq)**2-16.*s*q2ini))/8.
2794       qtq=4.*max(q2mass,q1)/(s-qq)
2795       if(qtq.lt.1.)then
2796         tmin=.5*s*qtq/(1.+sqrt(1.-qtq))
2797       else
2798         tmin=.5*s
2799       endif
2800       tmax=s/2.
2801 
2802       qlj=log(qq/q2min)*2.+1.
2803       j=int(qlj)
2804       if(j.lt.1)j=1
2805       if(j.gt.19)j=19
2806       wj(2)=qlj-j
2807       wj(3)=wj(2)*(wj(2)-1.)*.5
2808       wj(1)=1.-wj(2)+wj(3)
2809       wj(2)=wj(2)-2.*wj(3)
2810 
2811       qli=log(q1/qmin)/log(qmax/qmin)*20.+1.
2812       i=int(qli)
2813       if(i.lt.1)i=1
2814       if(i.gt.19)i=19
2815       wi(2)=qli-i
2816       wi(3)=wi(2)*(wi(2)-1.)*.5
2817       wi(1)=1.-wi(2)+wi(3)
2818       wi(2)=wi(2)-2.*wi(3)
2819 
2820       sl=log(s/s2min)/log(edmax/s2min)*25.+1.
2821       k=int(sl)
2822       if(k.lt.1)k=1
2823       if(k.gt.24)k=24
2824       wk(2)=sl-k
2825       wk(3)=wk(2)*(wk(2)-1.)*.5
2826       wk(1)=1.-wk(2)+wk(3)
2827       wk(2)=wk(2)-2.*wk(3)
2828 
2829       dsin1=0.
2830       do k1=1,3
2831         k2=k+k1-1+26*(m-1)+52*long
2832       do i1=1,3
2833       do j1=1,3
2834         dsin1=dsin1+csdsi(i+i1-1,j+j1-1,k2)*wi(i1)*wj(j1)*wk(k1)
2835       enddo
2836       enddo
2837       enddo
2838       if(m1.eq.0)then
2839         psdsin=psdsin+exp(dsin1)*(1./tmin-1./tmax)
2840       else
2841         psdsin=psdsin+max(0.,dsin1)
2842       endif
2843 
2844 2     continue
2845       if(long.eq.0.and.q1.lt.qq.and.s.gt.qq/(1.-q2ini/qq))then
2846         xx=qq/s
2847         dsi=psevi(q1,qq,xx,m,2)*4.*pi**2*alfe/s
2848         if(m1.eq.0)then
2849           psdsin=psdsin+max(dsi,0.)
2850         else
2851           dnsi=psevi(q1,qq,xx,3,2)*4.*pi**2*alfe/s
2852           psdsin=psdsin+max(dsi-dnsi,0.)
2853         endif
2854       endif
2855       return
2856       end
2857 
2858 
2859 
2860 
2861 
2862 
2863 
2864 
2865 
2866 
2867 
2868 
2869 
2870 
2871 c###########################################################################
2872 c###########################################################################
2873 c###########################################################################
2874 c###########################################################################
2875 c
2876 c                        unused 3P
2877 c
2878 c###########################################################################
2879 c###########################################################################
2880 c###########################################################################
2881 c###########################################################################
2882 
2883 
2884 cc------------------------------------------------------------------------
2885 c      function psvy(xpp0,xpr0,xpm0,xmr0,b,iqq)
2886 cc-----------------------------------------------------------------------
2887 cc psvy - 3p-contributions to the interaction eikonal
2888 cc xpp - lc+ for the pomeron,
2889 cc xpr - lc+ for the remnant,
2890 cc xpm - lc- for the pomeron,
2891 cc xpr - lc- for the remnant,
2892 cc b   - impact parameter,
2893 cc iqq=1  - Y-proj-uncut
2894 cc iqq=2  - Y-proj-1-cut
2895 cc iqq=3  - Y-proj-2-cut
2896 cc iqq=4  - Y-proj-soft-cut
2897 cc iqq=5  - Y-proj-gss-cut
2898 cc iqq=6  - Y-proj-qss-cut
2899 cc iqq=7  - Y-proj-ssg-cut
2900 cc iqq=8  - Y-proj-ssq-cut
2901 cc iqq=9  - Y-proj-difr
2902 cc iqq=-1 - Y-targ-uncut
2903 cc iqq=-2 - Y-targ-1-cut
2904 cc iqq=-3 - Y-targ-2-cut
2905 cc iqq=-4 - Y-targ-soft-cut
2906 cc iqq=-5 - Y-targ-gss-cut
2907 cc iqq=-6 - Y-targ-qss-cut
2908 cc iqq=-7 - Y-targ-ssg-cut
2909 cc iqq=-8 - Y-targ-ssq-cut
2910 cc iqq=-9 - Y-targ-difr
2911 cc------------------------------------------------------------------------
2912 c
2913 c      psvy=0.
2914 c      return
2915 c      end
2916 c
2917 cc------------------------------------------------------------------------
2918 c      function psvx(xpp,xpr,xpm,xmr,b,iqq)
2919 cc-----------------------------------------------------------------------
2920 cc psvx - 4p-contributions to the interaction eikonal
2921 cc xpp - lc+ for the pomeron,
2922 cc xpr - lc+ for the remnant,
2923 cc xpm - lc- for the pomeron,
2924 cc xpr - lc- for the remnant,
2925 cc b   - impact parameter,
2926 cc iqq=0  - X-uncut
2927 cc iqq=1  - X-1-cut
2928 cc iqq=2  - X-Y+cut
2929 cc iqq=-2 - X-Y-cut
2930 cc iqq=3  - X-1-cut-soft
2931 cc iqq=4  - X-1-cut-gss
2932 cc iqq=-4 - X-1-cut-gss
2933 cc iqq=5  - X-1-cut-qss
2934 cc iqq=-5 - X-1-cut-qss
2935 cc iqq=6  - X-difr+
2936 cc iqq=-6 - X-difr-
2937 cc------------------------------------------------------------------------
2938 c
2939 c      psvx=0.
2940 c      return
2941 c      end
2942 c
2943 c
2944 c
2945 c
2946 c
2947 c
2948 c
2949 c
2950 c
2951 c
2952 cc------------------------------------------------------------------------
2953 c      function psftig(x,xpomr,jj)
2954 cc-----------------------------------------------------------------------
2955 c
2956 c        psftig=0.
2957 c      end
2958 c
2959 c
2960 cc------------------------------------------------------------------------
2961 c      function psftih(zz,ddd)
2962 cc-----------------------------------------------------------------------
2963 c
2964 c      psftih=0.
2965 c      return
2966 c      end
2967 c
2968 cc------------------------------------------------------------------------
2969 c      function psftij(zz,ddd)
2970 cc------------------------------------------------------------------------
2971 c
2972 c      psftij=0.
2973 c      return
2974 c      end
2975 c
2976 cc------------------------------------------------------------------------
2977 c      function psftik(xp,del1,rh1,rh2,rh3,rp,z)
2978 cc-----------------------------------------------------------------------
2979 c
2980 c      psftik=0.
2981 c      return
2982 c      end
2983 c
2984 cc------------------------------------------------------------------------
2985 c      function psftim(xp1,xp,del1,rh1,rh2,rh3,rp,z)
2986 cc-----------------------------------------------------------------------
2987 c
2988 c      psftim=0.
2989 c      return
2990 c      end
2991 c
2992 cc------------------------------------------------------------------------
2993 c      function psftil(xp,del1,rh1,rh2,rh3,rp,z)
2994 cc------------------------------------------------------------------------
2995 c
2996 c      psftil=0.
2997 c      return
2998 c      end
2999 c
3000 cc------------------------------------------------------------------------
3001 c      function psftigt(x)
3002 cc-----------------------------------------------------------------------
3003 c
3004 c      psftigt=0.
3005 c       return
3006 c       end
3007 c
3008 cc------------------------------------------------------------------------
3009 c      function psftist(x)
3010 cc-----------------------------------------------------------------------
3011 c      psftist=0.
3012 c       return
3013 c       end
3014 c
3015 cc------------------------------------------------------------------------
3016 c      function psftig1(x,xpomr,jj)
3017 cc-----------------------------------------------------------------------
3018 c      psftig1=0.
3019 c      return
3020 c      end
3021 c
3022 cc------------------------------------------------------------------------
3023 c      function psftis1(x,xpomr,jj)
3024 cc-----------------------------------------------------------------------
3025 c      psftis1=0.
3026 c      return
3027 c      end
3028 c
3029 cc------------------------------------------------------------------------
3030 c      function psd3p1(xd,xpomr,qq,jj)
3031 cc-----------------------------------------------------------------------
3032 cc psd3p1 - df2difr/dx_pomr
3033 cc xd - bjorken x,
3034 cc xpomr - pomeron x,
3035 cc qq - photon virtuality
3036 cc jj=1 - 1st order
3037 cc-----------------------------------------------------------------------
3038 c
3039 c      psd3p1=0.
3040 c
3041 c      return
3042 c      end
3043 c
3044 cc------------------------------------------------------------------------
3045 c      function psv3p(sy,xpp,xpm,zb)
3046 cc-----------------------------------------------------------------------
3047 cc psv3p - 3p-contributions to the interaction eikonal
3048 cc sy  - energy squared for the hard interaction,
3049 cc xpp - lc+ for the sh pomeron,
3050 cc xpm - lc- for the sh pomeron,
3051 cc z   - impact parameter factor, z=exp(-b**2/4*rp),
3052 cc------------------------------------------------------------------------
3053 c
3054 c      psv3p=0.
3055 c      return
3056 c      end
3057 c
3058 cc------------------------------------------------------------------------
3059 c      function psfro(xpomr,zb,iclp,icdpro)
3060 cc-----------------------------------------------------------------------
3061 cc psfro - generalized froissaron between proj. and 3p-vertex
3062 cc xpp   - lc+ for the proj. side,
3063 cc xpomr - lc+ for the vertex,
3064 cc zb    - impact parameter factor, z=exp(-b**2/4*rp),
3065 cc------------------------------------------------------------------------
3066 c      psfro=0.
3067 c      return
3068 c      end
3069 c
3070 cc------------------------------------------------------------------------
3071 c      function psv2(xpomr,zb,iclp,iqq)
3072 cc-----------------------------------------------------------------------
3073 cc psv2 - 2-pom contribution to the froissaron
3074 cc xpomr - lc+ for the vertex,
3075 cc zb    - impact parameter factor, z=exp(-b**2/4*rp),
3076 cc------------------------------------------------------------------------
3077 c
3078 c      psv2=0.
3079 c      return
3080 c      end
3081 c
3082 cc------------------------------------------------------------------------
3083 c      function psvfro(xpp,xpr,xpomr,zb,iclp,icdpro,iqq)
3084 cc-----------------------------------------------------------------------
3085 cc psvfro - effective froissaron contributions
3086 cc xpomr - lc+ for the vertex,
3087 cc zb    - impact parameter factor, z=exp(-b**2/4*rp),
3088 cc iqq=0 - total uncut
3089 cc iqq=1 - total 1-cut
3090 cc iqq=2 - total 2-cut
3091 cc iqq=3 - soft 1-cut
3092 cc iqq=4 - gg 1-cut
3093 cc iqq=5 - qg 1-cut
3094 cc iqq=6 - difr
3095 cc------------------------------------------------------------------------
3096 c
3097 c      psvfro=0.
3098 c      return
3099 c      end
3100 c
3101 cc------------------------------------------------------------------------
3102 c      function psfroin(xpomr,z,iclp,icdpro)
3103 cc-----------------------------------------------------------------------
3104 cc psfroin - interpolation of effective froissaron corrections
3105 cc xpomr - lc+ for the 3p-vertex,
3106 cc z   - impact parameter factor, z=exp(-b**2/4*rp),
3107 cc-----------------------------------------------------------------------
3108 c
3109 c      psfroin=0.
3110 c      return
3111 c      end
3112 c
3113 cc------------------------------------------------------------------------
3114 c      function psvnorm(b)
3115 cc-----------------------------------------------------------------------
3116 cc psvnorm - X-contribution normalization
3117 cc b   - impact parameter
3118 cc-----------------------------------------------------------------------
3119 c
3120 c      psvnorm=0.
3121 c      return
3122 c      end
3123 c
3124 cc------------------------------------------------------------------------
3125 c      function psvxb(dxpp,dxpr,dxpm,dxmr,dxpomr,bb1,bb2
3126 c     *,iclp,iclt,icdpro,icdtar,iqq)
3127 cc-----------------------------------------------------------------------
3128 cc psvxb   - integrand for X-contributions
3129 cc dxpomr  - lc+ for the vertex,
3130 cc bb1,bb2 - impact parameters to proj(targ),
3131 cc iqq=0   - uncut
3132 cc iqq=1   - 1-cut
3133 cc iqq=2   - Y-cut
3134 cc iqq=3   - soft 1-cut
3135 cc iqq=4   - gg 1-cut
3136 cc iqq=5   - qg 1-cut
3137 cc iqq=6   - difr
3138 cc------------------------------------------------------------------------
3139 c      double precision dxpp,dxpr,dxpm,dxmr,dxpomr
3140 c
3141 c      psvxb=0.
3142 c      return
3143 c      end
3144 c
3145 cc------------------------------------------------------------------------
3146 c      function pscoef(zz,alp1,alp2,iclp,iqq)
3147 cc-----------------------------------------------------------------------
3148 cc pscoef - integrated vertexes
3149 cc zz=xpomr/xpr(xpp)
3150 cc iqq=0 - 2-cut
3151 cc iqq=1 - 1-uncut
3152 cc iqq=2 - 2-uncut
3153 cc------------------------------------------------------------------------
3154 c
3155 c      pscoef=0.
3156 c      return
3157 c      end
3158 c
3159 cc------------------------------------------------------------------------
3160 c      function pscoefi(z,i1,i2,iclp,iqq)
3161 cc-----------------------------------------------------------------------
3162 cc pscoefi - interpolation of integrated vertexes
3163 cc z=xpomr/xpr(xpp)
3164 cc iqq=0 - 2-cut
3165 cc iqq=1 - 1-uncut
3166 cc iqq=2 - 2-uncut
3167 cc------------------------------------------------------------------------
3168 c      pscoefi=0.
3169 c      return
3170 c      end
3171 c
3172 c