File indexing completed on 2024-04-06 12:14:04
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022
0023
0024
0025 SUBROUTINE phoGPHERAepo(imod)
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
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
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
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
0086
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
0095 call idmass(12,ams)
0096 ELEM = dble(ams)
0097 ELEM2 = ELEM**2
0098
0099 Q2MN = dble(qdmin)
0100 Q2MX = dble(qdmax)
0101
0102 YMIN2=dble(ydmin)
0103 YMAX2=dble(ydmax)
0104 XIMAX = LOG(YMAX2)
0105 XIMIN = LOG(YMIN2)
0106 XIDEL = XIMAX-XIMIN
0107
0108 THMIN2=dble(themin*pi/180.)
0109 THMAX2=dble(themax*pi/180.)
0110
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
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
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
0167
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
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
0201 if(ish.ge.2)WRITE(ifch,*) 'phoGPHERAepo: event with Q2,Y:',Q2,YY
0202
0203
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
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
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
0231 PFTHE = ACOS(PFIN(3)/PFIN(4))
0232
0233 IF(PFIN(4).GT.EEMIN2) THEN
0234 IF((PFTHE.LT.THMIN2).OR.(PFTHE.GT.THMAX2)) GOTO 175
0235 ENDIF
0236
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
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
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
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
0257 engy=sngl(GGECM)
0258 xbjevt=GQ2/GWD
0259 qsqevt=GQ2
0260
0261 rgampr(1) = P2(1)
0262 rgampr(2) = P2(2)
0263 rgampr(3) = P2(3)
0264 rgampr(4) = P2(4)
0265
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
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
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
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
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
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
0340
0341 endif
0342
0343
0344 END
0345
0346
0347
0348 DOUBLE PRECISION FUNCTION ALLM97(Q2,W)
0349
0350
0351
0352
0353
0354
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
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
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
0395 M02 = 0.31985D0
0396 LAM2 = 0.065270D0
0397 Q02 = 0.46017D0 +LAM2
0398
0399
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
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
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
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
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
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
0469 subroutine lepexp(rxbj,rqsq)
0470
0471
0472
0473
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
0534 subroutine fremny(wp1,wm1,pnx,pny,sm,ic1,ic2,ic3,ic4,coord,ey0)
0535
0536
0537
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
0647 subroutine psadis(iret)
0648
0649
0650
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
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
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
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
1150
1151
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
1274
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
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
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
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
1605
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
1617
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
1625
1626
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
1636
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
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
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
1947
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
1956
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
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
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
2156 subroutine psaevc
2157
2158 include 'epos.inc'
2159
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
2297 function psdbom(s,t,u,qq,long)
2298
2299
2300
2301
2302
2303
2304
2305
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
2318 function psdbin(q1,qq,s,m1,long)
2319
2320
2321
2322
2323
2324
2325
2326
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
2355 function psdbor(q1,qq,s,long)
2356
2357
2358
2359
2360
2361
2362
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
2400 subroutine psdint(s,qq,sds,sdn,sdb,sdt,sdr,m1,long)
2401
2402
2403
2404
2405
2406
2407
2408
2409
2410
2411
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
2522 function psdnsi(q1,qq,s,long)
2523
2524
2525
2526
2527
2528
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
2542 function psdrga(qq,s,s2min,j)
2543
2544
2545
2546
2547
2548
2549
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
2582 function psdres(qq,s,s2min,j)
2583
2584
2585
2586
2587
2588
2589
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
2602
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
2639 function psds(q1,qq,s,j,long)
2640
2641
2642
2643
2644
2645
2646
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
2743 function psdsj(q1,xx,s,qt,t,u,qq,j,long)
2744
2745
2746
2747
2748
2749
2750
2751
2752
2753
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
2763 function psdsin(q1,qq,s,m1,long)
2764
2765
2766
2767
2768
2769
2770
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
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
2872
2873
2874
2875
2876
2877
2878
2879
2880
2881
2882
2883
2884
2885
2886
2887
2888
2889
2890
2891
2892
2893
2894
2895
2896
2897
2898
2899
2900
2901
2902
2903
2904
2905
2906
2907
2908
2909
2910
2911
2912
2913
2914
2915
2916
2917
2918
2919
2920
2921
2922
2923
2924
2925
2926
2927
2928
2929
2930
2931
2932
2933
2934
2935
2936
2937
2938
2939
2940
2941
2942
2943
2944
2945
2946
2947
2948
2949
2950
2951
2952
2953
2954
2955
2956
2957
2958
2959
2960
2961
2962
2963
2964
2965
2966
2967
2968
2969
2970
2971
2972
2973
2974
2975
2976
2977
2978
2979
2980
2981
2982
2983
2984
2985
2986
2987
2988
2989
2990
2991
2992
2993
2994
2995
2996
2997
2998
2999
3000
3001
3002
3003
3004
3005
3006
3007
3008
3009
3010
3011
3012
3013
3014
3015
3016
3017
3018
3019
3020
3021
3022
3023
3024
3025
3026
3027
3028
3029
3030
3031
3032
3033
3034
3035
3036
3037
3038
3039
3040
3041
3042
3043
3044
3045
3046
3047
3048
3049
3050
3051
3052
3053
3054
3055
3056
3057
3058
3059
3060
3061
3062
3063
3064
3065
3066
3067
3068
3069
3070
3071
3072
3073
3074
3075
3076
3077
3078
3079
3080
3081
3082
3083
3084
3085
3086
3087
3088
3089
3090
3091
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101
3102
3103
3104
3105
3106
3107
3108
3109
3110
3111
3112
3113
3114
3115
3116
3117
3118
3119
3120
3121
3122
3123
3124
3125
3126
3127
3128
3129
3130
3131
3132
3133
3134
3135
3136
3137
3138
3139
3140
3141
3142
3143
3144
3145
3146
3147
3148
3149
3150
3151
3152
3153
3154
3155
3156
3157
3158
3159
3160
3161
3162
3163
3164
3165
3166
3167
3168
3169
3170
3171
3172
