File indexing completed on 2024-04-06 12:14:10
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
0026
0027
0028
0029
0030
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041
0042
0043
0044
0045
0046 SUBROUTINE QGSPSAINI
0047
0048
0049 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0050 INTEGER DEBUG
0051 CHARACTER *7 TY
0052 LOGICAL LCALC!,LSECT
0053
0054 DIMENSION EQ(17),MIJ(17,17,4),NIJ(17,17,4),CSJET(17,17,68),
0055 *CS1(17,17,68),GZ0(2),GZ1(3)
0056 COMMON /Q_XSECT/ GSECT(10,5,4)
0057 COMMON /Q_AREA1/ IA(2),ICZ,ICP
0058 COMMON /Q_AREA5/ RD(2),CR1(2),CR2(2),CR3(2)
0059
0060 COMMON /Q_AREA6/ PI,BM,AM
0061 COMMON /Q_AREA7/ RP1
0062 COMMON /Q_AREA10/ STMASS,AM0,AMN,AMK,AMC,AMLAMC,AMLAM,AMETA
0063 COMMON /Q_AREA15/ FP(5),RQ(5),CD(5)
0064 COMMON /Q_AREA16/ CC(5)
0065 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALFP,RR,SH,DELH
0066 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
0067 COMMON /Q_AREA19/ AHL(5)
0068
0069 COMMON /Q_AREA22/ SJV0,FJS0(5,3)
0070
0071 COMMON /Q_AREA23/ RJV(50)
0072 COMMON /Q_AREA24/ RJS(50,5,10)
0073 COMMON /Q_AREA27/ FP0(5)
0074 COMMON /Q_AREA28/ ARR(4)
0075 COMMON /Q_AREA29/ CSTOT(17,17,68)
0076 COMMON /Q_AREA30/ CS0(17,17,68)
0077 COMMON /Q_AREA31/ CSBORN(17,68)
0078 COMMON /Q_AREA32/ CSQ(17,2,2),CSBQ(17,2,2)
0079 COMMON /Q_AREA33/ FSUD(10,2)
0080 COMMON /Q_AREA34/ QRT(10,101,2)
0081 COMMON /Q_AREA35/ SJV(10,5),FJS(10,5,15)
0082 COMMON /Q_AREA39/ JCALC
0083 COMMON /Q_AREA40/ JDIFR
0084 COMMON /Q_AREA41/ TY(5)
0085 COMMON /Q_AREA43/ MONIOU
0086 COMMON /Q_DEBUG/ DEBUG
0087
0088 COMMON /Q_AREA44/ GZ(10,5,4),GZP(10,5,4)
0089
0090 COMMON /Q_AR1/ ANORM
0091 COMMON /Q_AR2/ RRR,RRRM
0092
0093
0094 COMMON /Q_AREA48/ QGSASECT(10,6,4)
0095
0096 COMMON /Q_VERSION/VERSION
0097 SAVE
0098
0099
0100 character*500 fndat,fnncs
0101 common/qgsfname/ fndat, fnncs, ifdat, ifncs
0102 common/qgsnfname/ nfndat, nfnncs
0103
0104 WRITE(MONIOU,100)
0105 100 FORMAT(' ',
0106 * '====================================================',
0107 * /,' ','| |',
0108 * /,' ','| QUARK GLUON STRING JET MODEL |',
0109 * /,' ','| |',
0110 * /,' ','| HADRONIC INTERACTION MONTE CARLO |',
0111 * /,' ','| BY |',
0112 * /,' ','| N.N. KALMYKOV AND S.S. OSTAPCHENKO |',
0113 * /,' ','| |',
0114 * /,' ','| e-mail: serg@eas.npi.msu.su |',
0115 * /,' ','| |',
0116 * /,' ','| Publication to be cited when using this program: |',
0117 * /,' ','| N.N. Kalmykov & S.S. Ostapchenko, A.I. Pavlov |',
0118 * /,' ','| Nucl. Phys. B (Proc. Suppl.) 52B (1997) 17 |',
0119 * /,' ','| |',
0120 * /,' ','====================================================',
0121 * /)
0122 IF(DEBUG.GE.1)WRITE (MONIOU,210)
0123 210 FORMAT(2X,'QGSPSAINI - MAIN INITIALIZATION PROCEDURE')
0124
0125 VERSION = 2003.
0126
0127
0128
0129 AHL(1)=1.D0-2.D0*ARR(1)
0130 AHL(2)=1.D0-ARR(1)-ARR(2)
0131 AHL(3)=1.D0-ARR(1)-ARR(3)
0132 AHL(4)=1.D0-ARR(1)-ARR(4)
0133 AHL(5)=AHL(2)+ARR(1)-ARR(4)
0134
0135
0136
0137
0138 CC(2)=1.D0/DSQRT(CD(2))
0139 CC(1)=1.D0/CC(2)/CD(1)
0140 CC(3)=1.D0/CC(2)/CD(3)
0141 CC(4)=1.D0/CC(2)/CD(4)
0142 CC(5)=1.D0/CC(2)/CD(5)
0143
0144
0145 FP0(2)=DSQRT(FP(2))
0146 FP0(1)=FP(1)/FP0(2)
0147 FP0(3)=FP(3)/FP0(2)
0148 FP0(4)=FP(4)/FP0(2)
0149 FP0(5)=FP(5)/FP0(2)
0150
0151
0152 SH=4.D0/QT0*PI
0153
0154 AQT0=DLOG(4.D0*QT0)
0155 QLOG=DLOG(QT0/ALM)
0156 QLL=DLOG(QLOG)
0157
0158
0159 IF(IFDAT.NE.1)THEN
0160 INQUIRE(FILE='QGSDAT01',EXIST=LCALC)
0161 ELSE !ctp
0162 INQUIRE(FILE=FNDAT(1:NFNDAT),EXIST=LCALC)
0163 ENDIF
0164 IF(LCALC)then
0165 IF(DEBUG.GE.1)WRITE (MONIOU,211)
0166 211 FORMAT(2X,'PSAINI: HARD CROSS SECTION RATIOS READOUT FROM THE'
0167 * ,' FILE QGSDAT01')
0168 IF(IFDAT.NE.1)THEN
0169 OPEN(1,FILE='QGSDAT01',STATUS='OLD')
0170 ELSE !ctp
0171 OPEN(IFDAT,FILE=FNDAT(1:NFNDAT),STATUS='OLD')
0172 ENDIF
0173 READ (1,*)CSBORN,CS0,CSTOT,CSQ,CSBQ,
0174 * FSUD,QRT,SJV,FJS,RJV,RJS,GZ,GZP,GSECT
0175 CLOSE(1)
0176 ELSE
0177
0178
0179
0180 WRITE (MONIOU,201)
0181 201 FORMAT(2X,'QGSPSAINI: HARD CROSS SECTIONS CALCULATION')
0182
0183
0184
0185
0186 DO 1 I=1,17
0187 1 EQ(I)=QT0*4.D0**FLOAT(I-1)
0188
0189 DO 2 I=1,17
0190
0191 QI=EQ(I)
0192
0193 DO 2 M=1,2
0194 DO 2 L=1,2
0195
0196 DO 2 K=1,17
0197 K1=K+17*(M-1)+34*(L-1)
0198 IF(K.LE.I.OR.K.EQ.2)THEN
0199 CSBORN(I,K1)=0.D0
0200 ELSE
0201
0202 SK=EQ(K)
0203
0204 CSBORN(I,K1)=QGSPSBORN(QI,SK,M-1,L-1)
0205 ENDIF
0206 2 CONTINUE
0207
0208
0209 DO 3 I=1,17
0210 DO 3 J=1,17
0211 N=MAX(I,J)
0212 DO 3 M=1,2
0213 DO 3 L=1,2
0214 ML=M+2*L-2
0215 DO 3 K=1,17
0216 K1=K+17*(M-1)+34*(L-1)
0217 CSJET(I,J,K1)=0.D0
0218 IF(K.LE.N.OR.K.EQ.2)THEN
0219 CSTOT(I,J,K1)=-80.D0
0220 CS0(I,J,K1)=-80.D0
0221 MIJ(I,J,ML)=K+1
0222 NIJ(I,J,ML)=K+1
0223 ELSE
0224 CSTOT(I,J,K1)=DLOG(CSBORN(N,K1))
0225 CS0(I,J,K1)=CSTOT(I,J,K1)
0226 ENDIF
0227 3 CONTINUE
0228
0229
0230 N=2
0231 4 CONTINUE
0232 IF(DEBUG.GE.2)WRITE (MONIOU,202)N,EQ(MIJ(1,1,1)),EQ(NIJ(1,1,1))
0233 202 FORMAT(2X,'PSAINI: NUMBER OF LADDER RUNS TO BE CONSIDERED:',I2/
0234 * 4X,'MINIMAL MASSES SQUARED FOR THE UNORDERED AND STRICTLY',
0235 * ' ORDERED LADDERS:'/4X,E10.3,3X,E10.3)
0236 DO 6 I=1,17
0237
0238 QI=EQ(I)
0239 DO 6 J=1,17
0240
0241 QJ=EQ(J)
0242
0243 QQ=MAX(QI,QJ)
0244
0245 S2MIN=MAX(QQ,4.D0*QT0)
0246 SM=DSQRT(QT0/S2MIN)
0247
0248 SMIN=S2MIN*(1.D0+SM)/(1.D0-SM)
0249
0250
0251 DO 6 M=1,2
0252 DO 6 L=1,2
0253 ML=M+2*L-2
0254
0255
0256 KMIN=NIJ(I,J,ML)
0257 IF(KMIN.LE.17)THEN
0258 DO 5 K=KMIN,17
0259 SK=EQ(K)
0260 IF(SK.LE.SMIN)THEN
0261 NIJ(I,J,ML)=NIJ(I,J,ML)+1
0262 ELSE
0263 K1=K+17*(M-1)+34*(L-1)
0264
0265
0266 CS1(I,J,K1)=QGSPSJET1(QI,QJ,SK,S2MIN,M-1,L)
0267 ENDIF
0268 5 CONTINUE
0269 ENDIF
0270 6 CONTINUE
0271
0272 DO 8 I=1,17
0273 DO 8 J=1,17
0274 DO 8 M=1,2
0275 DO 8 L=1,2
0276 ML=M+2*L-2
0277 KMIN=NIJ(I,J,ML)
0278 IF(KMIN.LE.17)THEN
0279 DO 7 K=KMIN,17
0280 K1=K+17*(M-1)+34*(L-1)
0281
0282
0283 CSJ=CS1(I,J,K1)+CSBORN(MAX(I,J),K1)
0284 IF(DEBUG.GE.2)WRITE (MONIOU,204)CSJ,EXP(CS0(I,J,K1))
0285 204 FORMAT(2X,'PSAINI: NEW AND OLD VALUES OF THE CONTRIBUTION',
0286 * ' OF THE STRICTLY ORDERED LADDER:'/4X,E10.3,3X,E10.3)
0287 IF(CSJ.EQ.0.D0.OR.ABS(1.D0-EXP(CS0(I,J,K1))/CSJ).LT.1.D-2)THEN
0288 NIJ(I,J,ML)=NIJ(I,J,ML)+1
0289 ELSE
0290
0291 CS0(I,J,K1)=DLOG(CSJ)
0292 ENDIF
0293 7 CONTINUE
0294 ENDIF
0295 8 CONTINUE
0296
0297 DO 10 I=1,17
0298 QI=EQ(I)
0299 DO 10 J=1,17
0300 QJ=EQ(J)
0301 QQ=MAX(QI,QJ)
0302 S2MIN=MAX(QQ,4.D0*QT0)
0303 SM=DSQRT(QT0/S2MIN)
0304
0305 SMIN=S2MIN*(1.D0+SM)/(1.D0-SM)
0306
0307 DO 10 M=1,2
0308 DO 10 L=1,2
0309 ML=M+2*L-2
0310
0311
0312 KMIN=MIJ(I,J,ML)
0313 IF(KMIN.LE.17)THEN
0314 DO 9 K=KMIN,17
0315 SK=EQ(K)
0316 IF(SK.LE.SMIN)THEN
0317 MIJ(I,J,ML)=MIJ(I,J,ML)+1
0318 ELSE
0319 K1=K+17*(M-1)+34*(L-1)
0320
0321
0322
0323 CS1(I,J,K1)=QGSPSJET(QI,QJ,SK,S2MIN,M-1,L)
0324 ENDIF
0325 9 CONTINUE
0326 ENDIF
0327 10 CONTINUE
0328
0329 DO 12 I=1,17
0330 DO 12 J=1,17
0331 DO 12 M=1,2
0332 DO 12 L=1,2
0333 ML=M+2*L-2
0334
0335 KMIN=MIJ(I,J,ML)
0336 IF(KMIN.LE.17)THEN
0337 DO 11 K=KMIN,17
0338 K1=K+17*(M-1)+34*(L-1)
0339 K2=K+17*(L-1)+34*(M-1)
0340 CSJ=CS1(I,J,K1)+EXP(CS0(J,I,K2))
0341 IF(CSJ.EQ.0.D0.OR.ABS(1.D0-EXP(CSTOT(I,J,K1))/CSJ).LT.1.D-2)
0342 * MIJ(I,J,ML)=MIJ(I,J,ML)+1
0343 IF(DEBUG.GE.2)WRITE (MONIOU,203)CSJ,EXP(CSTOT(I,J,K1))
0344 203 FORMAT(2X,'PSAINI: NEW AND OLD VALUES OF THE UNORDERED LADDER',
0345 * ' CROSS SECTION:'/4X,E10.3,3X,E10.3)
0346 11 CSTOT(I,J,K1)=DLOG(CSJ)
0347 ENDIF
0348 12 CONTINUE
0349
0350
0351 N=N+1
0352 DO 13 L=1,4
0353 13 IF(MIJ(1,1,L).LE.17.OR.NIJ(1,1,L).LE.17)GOTO 4
0354
0355
0356
0357 DO 14 I=1,17
0358 DO 14 K=1,17
0359 DO 14 M=1,2
0360 DO 14 L=1,2
0361 K1=K+17*(M-1)+34*(L-1)
0362 IF(K.LE.I.OR.K.EQ.2)THEN
0363 CSBORN(I,K1)=-80.D0
0364 ELSE
0365 CSBORN(I,K1)=DLOG(CSBORN(I,K1))
0366 ENDIF
0367 14 CONTINUE
0368
0369
0370
0371 DO 15 M=1,2
0372 DO 15 L=1,2
0373 DO 15 K=1,17
0374 IF(K.LE.2)THEN
0375 CSQ(K,M,L)=-80.D0
0376 CSBQ(K,M,L)=-80.D0
0377 ELSE
0378 K1=K+17*(M-1)+34*(L-1)
0379 CSBQ(K,M,L)=CSBORN(1,K1)
0380 CSQ(K,M,L)=CSTOT(1,1,K1)
0381 ENDIF
0382 15 CONTINUE
0383
0384
0385
0386
0387 DO 17 M=1,2
0388 FSUD(1,M)=0.D0
0389 DO 17 K=2,10
0390
0391
0392 QMAX=QTF*4.D0**(1.D0+K)
0393 17 FSUD(K,M)=PSUDT(QMAX,M-1)
0394
0395
0396
0397
0398
0399 DO 18 M=1,2
0400 DO 18 K=1,10
0401 QLMAX=1.38629D0*(K-1)
0402 QRT(K,1,M)=0.D0
0403 QRT(K,101,M)=QLMAX
0404 DO 18 I=1,99
0405 IF(K.EQ.1)THEN
0406 QRT(K,I+1,M)=0.D0
0407 ELSE
0408 QRT(K,I+1,M)=PSROOT(QLMAX,.01D0*I,M)
0409 ENDIF
0410 18 CONTINUE
0411
0412
0413 IF(DEBUG.GE.2)WRITE (MONIOU,205)
0414 205 FORMAT(2X,'QGSPSAINI: PRETABULATION OF THE INTERACTION EIKONALS')
0415
0416
0417
0418
0419
0420
0421
0422
0423 IA(1)=1
0424
0425 DO 21 IE=1,10
0426
0427 E0N=10.D0**IE
0428
0429
0430
0431 S=2.D0*E0N*AMN
0432
0433 Y0=DLOG(S)
0434
0435
0436 DO 21 ICZ=1,5
0437
0438 RS=RQ(ICZ)+ALFP*Y0
0439
0440 RS0=RQ(ICZ)
0441
0442
0443 FS=FP(ICZ)*EXP(Y0*DEL)/RS*CD(ICZ)
0444
0445 RP1=RS*4.D0*.0391D0/AM**2
0446
0447 G0=PI*RP1/CD(ICZ)*AM**2*10.D0
0448
0449 SJV(IE,ICZ)=QGSPSHARD(S,ICZ)
0450 SJV0=SJV(IE,ICZ)
0451
0452 DO 19 I=1,5
0453 DO 19 M=1,3
0454 Z=.2D0*I
0455
0456
0457
0458 M1=M+3*(ICZ-1)
0459 FJS(IE,I,M1)=DLOG(QGSPSFSH(S,Z,ICZ,M-1)/Z)
0460 FJS0(I,M)=FJS(IE,I,M1)
0461 19 CONTINUE
0462
0463 DO 20 IIA=1,4
0464
0465 IA(2)=4**(IIA-1)
0466 IF(DEBUG.GE.1)WRITE (MONIOU,206)E0N,TY(ICZ),IA(2)
0467 206 FORMAT(2X,'QGSPSAINI: INITIAL PARTICLE ENERGY:',E10.3,2X,
0468 *'ITS TYPE:',A7,2X,'TARGET MASS NUMBER:',I2)
0469
0470
0471 IF(IA(2).GT.10)THEN
0472
0473 RD(2)=0.7D0*FLOAT(IA(2))**.446/AM
0474 ELSE
0475
0476 RD(2)=.9D0*FLOAT(IA(2))**.3333/AM
0477 ENDIF
0478
0479 IF(IA(2).EQ.1)THEN
0480
0481
0482 BM=2.D0*DSQRT(RP1)
0483
0484
0485
0486
0487 CALL XXFZ(0.D0,GZ0)
0488 if (debug .ge.1) write (moniou,*)gz0
0489
0490 GTOT=G0*GZ0(1)
0491
0492 GABS=G0*GZ0(2)*.5D0
0493
0494 GD0=GTOT-GABS
0495
0496 GDP=(1.D0-CC(ICZ))*CC(2)*GD0
0497
0498 GDT=(1.D0-CC(2))*CC(ICZ)*GD0
0499
0500 GDD=(1.D0-CC(ICZ))*(1.D0-CC(2))*GD0
0501
0502 GIN=GABS+GDP+GDT+GDD
0503 GEL=GD0*CC(ICZ)*CC(2)
0504
0505 IF(DEBUG.GE.1)WRITE (MONIOU,225)GTOT,GIN,GEL,GDP,GDT,GDD
0506
0507 225 FORMAT(2X,'QGSPSAINI: HADRON-PROTON CROSS SECTIONS:'/
0508 * 4X,'GTOT=',E10.3,2X,'GIN=',E10.3,2X,'GEL=',E10.3/4X,
0509 * 'GDIFR_PROJ=',E10.3,2X,'GDIFR_TARG=',E10.3,2X,
0510 * 'G_DOUBLE_DIFR',E10.3)
0511
0512 GZ(IE,ICZ,IIA)=GDT/GIN
0513 GZP(IE,ICZ,IIA)=(GDP+GDD)/GIN ! so00
0514
0515 GSECT(IE,ICZ,IIA)=LOG(GIN)
0516
0517 ELSE
0518
0519
0520
0521 BM=RD(2)+DLOG(29.D0)
0522
0523 RRR=RD(2)
0524
0525 RRRM=RRR+DLOG(9.D0)
0526
0527 ANORM=1.5D0/PI/RRR**3/(1.D0+(PI/RRR)**2)*RP1
0528
0529
0530
0531 CALL XXGAU(GZ1)
0532
0533
0534 CALL XXGAU1(GZ1)
0535
0536 GIN=GZ1(1)+GZ1(2)+GZ1(3)
0537
0538 IF(DEBUG.GE.1)WRITE (MONIOU,224)
0539 * GIN*10.D0,GZ1(1)*10.D0,GZ1(2)*10.D0
0540
0541 224 FORMAT(2X,'QGSPSAINI: HADRON-NUCLEUS CROSS SECTIONS:'/
0542 * 4X,'GIN=',E10.3,2X,'GDIFR_TARG=',E10.3,2X,
0543 * 'GDIFR_PROJ=',E10.3)
0544
0545 GZ(IE,ICZ,IIA)=GZ1(1)/GIN
0546 GZP(IE,ICZ,IIA)=GZ1(2)/GIN ! so00
0547
0548 GIN=GIN*10.
0549 GSECT(IE,ICZ,IIA)=LOG(GIN)
0550
0551 ENDIF
0552 20 CONTINUE
0553 21 CONTINUE
0554
0555
0556 DO 23 I=1,50
0557
0558 YJ=AQT0+.5D0*I
0559
0560 RJV(I)=PSREJV(EXP(YJ))
0561
0562 DO 22 J=1,5
0563 DO 22 M=1,2
0564 Z=.2D0*J
0565 DO 22 ICZ=1,5
0566
0567 RS0=RQ(ICZ)
0568 M1=M+2*(ICZ-1)
0569
0570
0571 RJS(I,J,M1)=PSREJS(EXP(YJ),Z,M-1)
0572 22 CONTINUE
0573 23 CONTINUE
0574
0575
0576 WRITE (MONIOU,212)
0577 212 FORMAT(2X,'PSAINI: HARD CROSS SECTIONS ARE WRITTEN TO THE FILE'
0578 * ,' QGSDAT01')
0579 IF(IFDAT.NE.1)THEN
0580 OPEN(1,FILE='QGSDAT01',STATUS='unknown')
0581 ELSE !ctp
0582 OPEN(IFDAT,FILE=FNDAT(1:NFNDAT),STATUS='unknown')
0583 ENDIF
0584 WRITE (1,*)CSBORN,CS0,CSTOT,CSQ,CSBQ,
0585 * FSUD,QRT,SJV,FJS,RJV,RJS,GZ,GZP,GSECT
0586 CLOSE(1)
0587 ENDIF
0588
0589
0590
0591
0592 IF(IFNCS.NE.2)THEN
0593 INQUIRE(FILE='SECTNU',EXIST=LCALC)
0594 ELSE !ctp
0595 INQUIRE(FILE=FNNCS(1:NFNNCS),EXIST=LCALC)
0596 ENDIF
0597 IF(LCALC)then
0598 IF(DEBUG.GE.1)WRITE (MONIOU,208)
0599 208 FORMAT(2X,'PSAINI: NUCLEAR CROSS SECTIONS READOUT FROM THE FILE'
0600 * ,' SECTNU')
0601 IF(IFNCS.NE.2)THEN
0602 OPEN(2,FILE='SECTNU',STATUS='OLD')
0603 ELSE !ctp
0604 OPEN(IFNCS,FILE=FNNCS(1:NFNNCS),STATUS='OLD')
0605 ENDIF
0606 READ (2,*)QGSASECT
0607 CLOSE(2)
0608 ELSE
0609
0610 NITER=5000 !NUMBER OF ITERATIONS
0611 DO IE=1,10
0612 E0N=10.D0**IE
0613 DO IIA1=1,6
0614 IAP=2**IIA1
0615 DO IIA2=1,4
0616 IAT=4**(IIA2-1)
0617
0618 WRITE (MONIOU,207)E0N,IAP,IAT
0619 207 FORMAT(2X,'QGSPSAINI: INITIAL NUCLEUS ENERGY:',E10.3,2X,
0620 * 'PROJECTILE MASS:',I2,2X,'TARGET MASS:',I2)
0621 CALL XXAINI(E0N,2,IAP,IAT)
0622 CALL CROSSC(NITER,GTOT,GPROD,GABS,GDD,GQEL,GCOH)
0623 IF(DEBUG.GE.1)WRITE (MONIOU,209)
0624 * GTOT,GPROD,GABS,GDD,GQEL,GCOH
0625
0626 209 FORMAT(2X,'GTOT',D10.3,' GPROD',D10.3,' GABS',D10.3/2X,
0627 * 'GDD',D10.3,' GQEL',D10.3,' GCOH',D10.3)
0628 QGSASECT(IE,IIA1,IIA2)=LOG(GPROD)
0629 ENDDO
0630 ENDDO
0631 ENDDO
0632 IF(IFNCS.NE.2)THEN
0633 OPEN(2,FILE='SECTNU',STATUS='UNKNOWN')
0634 ELSE !ctp
0635 OPEN(IFNCS,FILE=FNNCS(1:NFNNCS),STATUS='UNKNOWN')
0636 ENDIF
0637 WRITE (2,*)QGSASECT
0638 CLOSE(2)
0639 ENDIF
0640
0641 IF(DEBUG.GE.3)WRITE (MONIOU,218)
0642 218 FORMAT(2X,'QGSPSAINI - END')
0643 RETURN
0644 END
0645
0646
0647 FUNCTION PSAPINT(X,J,L)
0648
0649
0650
0651
0652
0653 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0654 INTEGER DEBUG
0655 COMMON /Q_AREA43/ MONIOU
0656 COMMON /Q_DEBUG/ DEBUG
0657 SAVE
0658
0659 IF(DEBUG.GE.2)WRITE (MONIOU,201)X,J,L
0660 201 FORMAT(2X,'PSAPINT: X=',E10.3,2X,'J= ',I1,2X,'L= ',I1)
0661 IF(J.EQ.0)THEN
0662 IF(L.EQ.0)THEN
0663 PSAPINT=6.D0*(DLOG(X/(1.D0-X))-X**3/3.D0+X**2/2.D0-2.D0*X)
0664 ELSE
0665 PSAPINT=3.D0*(X+X**3/1.5D0-X*X)
0666 ENDIF
0667 ELSE
0668 IF(L.EQ.0)THEN
0669 PSAPINT=(DLOG(X)-X+.25D0*X*X)/.375D0
0670 ELSE
0671 Z=1.D0-X
0672 PSAPINT=-(DLOG(Z)-Z+.25D0*Z*Z)/.375D0
0673 ENDIF
0674 ENDIF
0675 IF(DEBUG.GE.2)WRITE (MONIOU,202)PSAPINT
0676 202 FORMAT(2X,'PSAPINT=',E10.3)
0677 RETURN
0678 END
0679
0680
0681 SUBROUTINE PSASET
0682
0683
0684 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0685 INTEGER DEBUG
0686 CHARACTER*7 TY
0687 COMMON /Q_AREA15/ FP(5),RQ(5),CD(5)
0688 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALFP,RR,SH,DELH
0689 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
0690 COMMON /Q_AREA25/ AHV(5)
0691 COMMON /Q_AREA26/ FACTORK
0692 COMMON /Q_AREA41/ TY(5)
0693 COMMON /Q_AREA43/ MONIOU
0694 COMMON /Q_DEBUG/ DEBUG
0695 COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX !ctp
0696 DIMENSION XA(210,3),XB(210,3)
0697 SAVE
0698
0699 IF(DEBUG.GE.1)WRITE (MONIOU,210)
0700 210 FORMAT(2X,'PSASET - COMMON MODEL PARAMETERS SETTING')
0701
0702 BQGS=0.d0 !ctp used to link with nexus
0703 BMAXQGS=0.d0 !ctp used to link with nexus
0704 BMAXNEX=-1.d0 !ctp used to link with nexus
0705 BMINNEX=0.d0 !ctp used to link with nexus
0706
0707
0708
0709
0710
0711
0712
0713
0714
0715 DEL=.07D0
0716 ALFP=.21D0
0717
0718 FP(1)=2.43D0
0719 RQ(1)=2.4D0
0720 CD(1)=1.6D0
0721
0722 FP(2)=3.64D0
0723 RQ(2)=3.56D0
0724 CD(2)=1.5D0
0725
0726 FP(3)=1.75D0
0727 RQ(3)=2.D0
0728 CD(3)=1.7D0
0729
0730 FP(4)=1.21D0
0731 RQ(4)=1.78D0
0732 CD(4)=2.0D0
0733
0734 FP(5)=2.43D0
0735 RQ(5)=2.4D0
0736 CD(5)=2.0D0
0737
0738
0739
0740
0741
0742
0743
0744
0745
0746
0747
0748
0749
0750
0751 ALM=.04D0
0752 RR=.35D0 ! produces 76 mbarn for p-pbar at Tevatron energies
0753
0754 QT0=4.D0
0755 BET=1.D0
0756 DELH=0.25D0
0757 AMJ0=0.D0
0758 QTF=.5D0
0759 FACTORK=2.D0
0760
0761
0762
0763
0764 AHV(1)=1.5D0
0765 AHV(2)=2.5D0
0766 AHV(3)=2.D0
0767 AHV(4)=4.D0
0768 AHV(5)=5.D0
0769
0770 TY(1)='pion '
0771 TY(2)='nucleon'
0772 TY(3)='kaon '
0773 TY(4)='D-meson'
0774 TY(5)='LambdaC'
0775 RETURN
0776 END
0777
0778
0779 FUNCTION QGSPSBINT(QQ,S,M,L)
0780
0781
0782
0783
0784
0785
0786 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0787 INTEGER DEBUG
0788 DIMENSION WI(3),WK(3)
0789 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
0790 COMMON /Q_AREA31/ CSJ(17,68)
0791 COMMON /Q_AREA43/ MONIOU
0792 COMMON /Q_DEBUG/ DEBUG
0793 SAVE
0794
0795 IF(DEBUG.GE.2)WRITE (MONIOU,201)QQ,S,M,L
0796 201 FORMAT(2X,'QGSPSBINT: QQ=',E10.3,2X,'S= ',E10.3,2X,'M= ',I1,2X,
0797 * 'L= ',I1)
0798 QGSPSBINT=0.D0
0799 IF(S.LE.MAX(4.D0*QT0,QQ))THEN
0800 IF(DEBUG.GE.3)WRITE (MONIOU,202)QGSPSBINT
0801 202 FORMAT(2X,'QGSPSBINT=',E10.3)
0802 RETURN
0803 ENDIF
0804
0805 ML=17*(M-1)+34*(L-1)
0806 QLI=DLOG(QQ/QT0)/1.38629d0
0807 SL=DLOG(S/QT0)/1.38629d0
0808 SQL=SL-QLI
0809 I=INT(QLI)
0810 K=INT(SL)
0811 IF(I.GT.13)I=13
0812
0813 IF(SQL.GT.10.D0)THEN
0814 IF(K.GT.14)K=14
0815 WI(2)=QLI-I
0816 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
0817 WI(1)=1.D0-WI(2)+WI(3)
0818 WI(2)=WI(2)-2.D0*WI(3)
0819 WK(2)=SL-K
0820 WK(3)=WK(2)*(WK(2)-1.D0)*.5D0
0821 WK(1)=1.D0-WK(2)+WK(3)
0822 WK(2)=WK(2)-2.D0*WK(3)
0823
0824 DO 1 I1=1,3
0825 DO 1 K1=1,3
0826 1 QGSPSBINT=QGSPSBINT+CSJ(I+I1,K+K1+ML)*WI(I1)*WK(K1)
0827 QGSPSBINT=EXP(QGSPSBINT)
0828 ELSEIF(SQL.LT.1.D0.AND.I.NE.0)THEN
0829 SQ=(S/QQ-1.D0)/3.D0
0830 WI(2)=QLI-I
0831 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
0832 WI(1)=1.D0-WI(2)+WI(3)
0833 WI(2)=WI(2)-2.D0*WI(3)
0834
0835 DO 2 I1=1,3
0836 I2=I+I1
0837 K2=I2+1+ML
0838 2 QGSPSBINT=QGSPSBINT+CSJ(I2,K2)*WI(I1)
0839 QGSPSBINT=EXP(QGSPSBINT)*SQ
0840 ELSEIF(K.EQ.1)THEN
0841 SQ=(S/QT0/4.D0-1.D0)/3.D0
0842 WI(2)=QLI
0843 WI(1)=1.D0-QLI
0844
0845 DO 3 I1=1,2
0846 3 QGSPSBINT=QGSPSBINT+CSJ(I1,3+ML)*WI(I1)
0847 QGSPSBINT=EXP(QGSPSBINT)*SQ
0848 ELSEIF(K.LT.15)THEN
0849 KL=INT(SQL)
0850 IF(I+KL.GT.12)I=12-KL
0851 IF(I+KL.EQ.1)KL=2
0852 WI(2)=QLI-I
0853 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
0854 WI(1)=1.D0-WI(2)+WI(3)
0855 WI(2)=WI(2)-2.D0*WI(3)
0856 WK(2)=SQL-KL
0857 WK(3)=WK(2)*(WK(2)-1.D0)*.5D0
0858 WK(1)=1.D0-WK(2)+WK(3)
0859 WK(2)=WK(2)-2.D0*WK(3)
0860
0861 DO 4 I1=1,3
0862 I2=I+I1
0863 DO 4 K1=1,3
0864 K2=I2+KL+K1-1+ML
0865 4 QGSPSBINT=QGSPSBINT+CSJ(I2,K2)*WI(I1)*WK(K1)
0866 QGSPSBINT=EXP(QGSPSBINT)
0867
0868 ELSE
0869 K=15
0870 IF(I.GT.K-3)I=K-3
0871 WI(2)=QLI-I
0872 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
0873 WI(1)=1.D0-WI(2)+WI(3)
0874 WI(2)=WI(2)-2.D0*WI(3)
0875 WK(2)=SL-K
0876 WK(1)=1.D0-WK(2)
0877
0878 DO 5 I1=1,3
0879 DO 5 K1=1,2
0880 5 QGSPSBINT=QGSPSBINT+CSJ(I+I1,K+K1+ML)*WI(I1)*WK(K1)
0881 QGSPSBINT=EXP(QGSPSBINT)
0882 ENDIF
0883 IF(DEBUG.GE.3)WRITE (MONIOU,202)QGSPSBINT
0884 RETURN
0885 END
0886
0887
0888 FUNCTION QGSPSBORN(QQ,S,IQ1,IQ2)
0889
0890
0891
0892
0893
0894 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0895 INTEGER DEBUG
0896 COMMON /Q_AREA6/ PI,BM,AM
0897 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
0898 COMMON /Q_AREA26/ FACTORK
0899 COMMON /Q_AREA43/ MONIOU
0900 COMMON /Q_DEBUG/ DEBUG
0901 COMMON /Q_AR13/ X1(7),A1(7)
0902 SAVE
0903
0904 IF(DEBUG.GE.2)WRITE (MONIOU,201)QQ,S,IQ1,IQ2
0905 201 FORMAT(2X,'QGSPSBORN: QQ=',E10.3,2X,'S= ',E10.3,2X,'IQ1= ',I1
0906 * ,2X,'IQ2= ',I1)
0907 TMIN=S*(.5D0-DSQRT(.25D0-QT0/S))
0908 TMIN=MAX(TMIN,S*QQ/(S+QQ))
0909
0910 IF(IQ1*IQ2.EQ.0)THEN
0911 IQ=IQ2
0912 ELSE
0913 IQ=2
0914 ENDIF
0915
0916 QGSPSBORN=0.D0
0917 DO 1 I=1,7
0918 DO 1 M=1,2
0919 T=2.D0*TMIN/(1.D0+2.D0*TMIN/S-X1(I)*(2*M-3)*(1.D0-2.D0*TMIN/S))
0920 QT=T*(1.D0-T/S)
0921 FB=PSFBORN(S,T,IQ1,IQ)+PSFBORN(S,S-T,IQ1,IQ)
0922 1 QGSPSBORN=QGSPSBORN+A1(I)*FB/DLOG(QT/ALM)**2*T**2
0923 QGSPSBORN=QGSPSBORN*(.5D0/TMIN-1.D0/S)*FACTORK*PI**3/2.25D0**2
0924 & /S**2
0925 IF(IQ1.EQ.0.AND.IQ2.EQ.0)QGSPSBORN=QGSPSBORN*.5D0
0926 IF(DEBUG.GE.3)WRITE (MONIOU,202)QGSPSBORN
0927 202 FORMAT(2X,'QGSPSBORN=',E10.3)
0928 RETURN
0929 END
0930
0931
0932 SUBROUTINE PSCAJET(QQ,IQ1,QV,ZV,QM,IQV,LDAU,LPAR,JQ)
0933
0934
0935
0936
0937
0938
0939
0940
0941
0942
0943
0944
0945
0946
0947
0948
0949
0950
0951 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
0952 INTEGER DEBUG
0953 DIMENSION QMAX(30,50),IQM(2),LNV(50),
0954 * QV(30,50),ZV(30,50),QM(30,50),IQV(30,50),
0955 * LDAU(30,49),LPAR(30,50)
0956
0957 COMMON /Q_AREA11/ B10
0958 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
0959 COMMON /Q_AREA43/ MONIOU
0960 COMMON /Q_DEBUG/ DEBUG
0961
0962 SAVE
0963 EXTERNAL PSRAN
0964
0965 IF(DEBUG.GE.2)WRITE (MONIOU,201)QQ,IQ1,JQ
0966 201 FORMAT(2X,'PSCAJET: QQ=',E10.3,2X,'IQ1= ',I1,2X,'JQ=',I1)
0967
0968 DO 1 I=2,20
0969 1 LNV(I)=0
0970 LNV(1)=1
0971 QMAX(1,1)=QQ
0972 IQV(1,1)=IQ1
0973 NLEV=1
0974 NROW=1
0975
0976 2 QLMAX=DLOG(QMAX(NROW,NLEV)/QTF/16.D0)
0977 IQ=MIN(1,IABS(IQV(NROW,NLEV)))+1
0978
0979 IF(PSRAN(B10).GT.PSUDINT(QLMAX,IQ))THEN
0980 Q=PSQINT(QLMAX,PSRAN(B10),IQ)
0981 Z=PSZSIM(Q,IQ)
0982
0983 LL=LNV(NLEV+1)+1
0984 LDAU(NROW,NLEV)=LL
0985 LPAR(LL,NLEV+1)=NROW
0986 LPAR(LL+1,NLEV+1)=NROW
0987 LNV(NLEV+1)=LL+1
0988
0989 IF(IQ.NE.1)THEN
0990 IF((3-2*JQ)*IQV(NROW,NLEV).GT.0)THEN
0991 IQM(1)=0
0992 IQM(2)=IQV(NROW,NLEV)
0993 ELSE
0994 IQM(2)=0
0995 IQM(1)=IQV(NROW,NLEV)
0996 Z=1.D0-Z
0997 ENDIF
0998 ELSE
0999
1000 WG=QGSPSFAP(Z,0,0)
1001
1002 WG=WG/(WG+QGSPSFAP(Z,0,1))
1003 IF(PSRAN(B10).LT.WG)THEN
1004 IQM(1)=0
1005 IQM(2)=0
1006 ELSE
1007 IQM(1)=INT(3.D0*PSRAN(B10)+1.D0)*(3-2*JQ)
1008 IQM(2)=-IQM(1)
1009 ENDIF
1010 IF(PSRAN(B10).LT..5D0)Z=1.D0-Z
1011 ENDIF
1012
1013 QV(NROW,NLEV)=Q
1014 ZV(NROW,NLEV)=Z
1015
1016 NROW=LL
1017 NLEV=NLEV+1
1018 QMAX(NROW,NLEV)=Q*Z**2
1019 QMAX(NROW+1,NLEV)=Q*(1.D0-Z)**2
1020 IQV(NROW,NLEV)=IQM(1)
1021 IQV(NROW+1,NLEV)=IQM(2)
1022 IF(DEBUG.GE.3)WRITE (MONIOU,203)NLEV,NROW,Q,Z
1023 203 FORMAT(2X,'PSCAJET: NEW BRANCHING AT LEVEL NLEV=',I2,
1024 * ' NROW=',I2/4X,' EFFECTIVE MOMENTUM Q=',E10.3,2X,' Z=',E10.3)
1025 GOTO 2
1026 ELSE
1027
1028 QV(NROW,NLEV)=0.D0
1029 ZV(NROW,NLEV)=0.D0
1030 QM(NROW,NLEV)=AMJ0
1031 IF(DEBUG.GE.3)WRITE (MONIOU,204)NLEV,NROW
1032 204 FORMAT(2X,'PSCAJET: NEW FINAL JET AT LEVEL NLEV=',I2,
1033 * ' NROW=',I2)
1034 ENDIF
1035
1036 4 CONTINUE
1037 IF(NLEV.EQ.1)THEN
1038 IF(DEBUG.GE.3)WRITE (MONIOU,202)
1039 202 FORMAT(2X,'PSCAJET - END')
1040 RETURN
1041 ENDIF
1042 LPROW=LPAR(NROW,NLEV)
1043
1044 IF(LDAU(LPROW,NLEV-1).EQ.NROW)THEN
1045 NROW=NROW+1
1046 GOTO 2
1047 ELSE
1048 Z=ZV(LPROW,NLEV-1)
1049 QM(LPROW,NLEV-1)=Z*(1.D0-Z)*QV(LPROW,NLEV-1)
1050 * +QM(NROW-1,NLEV)/Z+QM(NROW,NLEV)/(1.D0-Z)
1051 NROW=LPROW
1052 NLEV=NLEV-1
1053 IF(DEBUG.GE.3)WRITE (MONIOU,205)NLEV,NROW,QM(LPROW,NLEV)
1054 205 FORMAT(2X,'PSCAJET: JET MASS AT LEVEL NLEV=',I2,
1055 * ' NROW=',I2,' - QM=',E10.3)
1056 GOTO 4
1057 ENDIF
1058 END
1059
1060
1061 SUBROUTINE PSCONF
1062
1063
1064
1065 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1066 INTEGER DEBUG
1067
1068
1069
1070
1071
1072 DIMENSION XA(210,3),XB(210,3),FHARD(3)
1073 COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX !ctp
1074 COMMON /Q_AREA1/ IA(2),ICZ,ICP
1075 COMMON /Q_AREA2/ S,Y0,WP0,WM0
1076 COMMON /Q_AREA6/ PI,BM,AM
1077
1078
1079
1080
1081
1082
1083
1084
1085
1086
1087
1088
1089
1090 COMMON /Q_AREA9/ LQA(56),LQB(56),NQS(1000),IAS(1000),IBS(1000),
1091 * LHA(56),LHB(56),ZH(4000),IAH(4000),IBH(4000),
1092 * IQH(4000),LVA(56),LVB(56)
1093 COMMON /Q_AREA10/ STMASS,AM0,AMN,AMK,AMC,AMLAMC,AMLAM,AMETA
1094 COMMON /Q_AREA11/ B10
1095
1096 COMMON /Q_AREA12/ NSP
1097 COMMON /Q_AREA16/ CC(5)
1098 COMMON /Q_AREA40/ JDIFR
1099 COMMON /Q_AREA43/ MONIOU
1100
1101 COMMON /Q_AREA45/ GDT,GDP !so00
1102
1103 COMMON /Q_AREA99/ NWT
1104 COMMON /Q_DEBUG/ DEBUG
1105
1106
1107 integer ng1evt,ng2evt,ikoevt
1108 real rglevt,sglevt,eglevt,fglevt,typevt
1109 common/c2evt/ng1evt,ng2evt,rglevt,sglevt,eglevt,fglevt,ikoevt
1110 *,typevt !in epos.inc
1111
1112
1113 DIMENSION IWT(56)
1114 SAVE
1115 EXTERNAL PSRAN
1116
1117 IF(DEBUG.GE.1)WRITE (MONIOU,201)
1118 201 FORMAT(2X,'PSCONF - CONFIGURATION OF THE INTERACTION')
1119
1120 100 NSP=0
1121 typevt=1
1122 IF(IA(1).EQ.1)THEN
1123
1124 IF(JDIFR.EQ.1.AND.PSRAN(B10).LT.GDT)THEN
1125
1126 IF(IA(2).NE.1)THEN
1127
1128 ICT=INT(2.5+PSRAN(B10))
1129 ELSE
1130
1131 ICT=2
1132 ENDIF
1133 WPI=WP0
1134 WMI=WM0
1135
1136 CALL XXDTG(WPI,WMI,ICP,ICT,0)
1137 typevt=-4
1138 goto 21 !so00
1139 ELSEIF(ABS(JDIFR).EQ.1.AND.PSRAN(B10).LT.GDP)THEN !so00
1140 IF(IA(2).NE.1)THEN !so00
1141
1142 ICT=INT(2.5+PSRAN(B10)) !so00
1143 ELSE !so00
1144
1145 ICT=2 !so00
1146 ENDIF !so00
1147 IF(DEBUG.GE.2)WRITE (MONIOU,206) !so00
1148 206 FORMAT(2X,'PROJECTILE HADRON DIFFRACTION') !so00
1149 ICP0=ICP !so00
1150 WPI=WP0 !so00
1151 WMI=WM0 !so00
1152 LQ=0 !so00
1153 CALL XXDPR(WPI,WMI,ICP0,ICT,LQ) !so00
1154 typevt=4
1155 goto 21 !so00
1156 ENDIF
1157
1158
1159
1160 DO 1 I=1,3
1161 1 XA(1,I)=0.D0
1162 ENDIF
1163
1164
1165
1166
1167
1168
1169
1170
1171
1172
1173
1174
1175
1176
1177
1178
1179
1180
1181
1182
1183
1184
1185
1186 DO 3 I=1,IA(1)
1187 LHA(I)=0
1188 LVA(I)=0
1189 3 LQA(I)=0
1190 DO 4 I=1,IA(2)
1191 LHB(I)=0
1192 LVB(I)=0
1193 4 LQB(I)=0
1194
1195
1196
1197 5 CONTINUE
1198
1199 IF(IA(2).NE.1)THEN !changed!!!!!!!!! dh 8/10/98
1200
1201
1202
1203 NT=0
1204 DO 6 I=1,IA(2)
1205 6 NT=NT+INT(CC(2)+PSRAN(B10))
1206
1207 IF(NT.EQ.0)GOTO 5
1208 IF(DEBUG.GE.3)WRITE (MONIOU,203)NT
1209 203 FORMAT(2X,'PSCONF: NUMBER OF ACTIVE TARGET NUCLEONS NT=',
1210 * I2)
1211
1212
1213 CALL PSGEA(IA(2),XB,2) !changed!!!!!!!!! 25.03.99
1214
1215
1216
1217 ELSE !changed!!!!!!!!! dh 8/10/98
1218 NT=1 !changed!!!!!!!!! dh 8/10/98
1219 XB(1,1)=0.D0 !changed!!!!!!!!! dh 8/10/98
1220 XB(1,2)=0.D0 !changed!!!!!!!!! dh 8/10/98
1221 ENDIF !changed!!!!!!!!! dh 8/10/98
1222
1223
1224
1225
1226 B=BM*DSQRT(PSRAN(B10))
1227 IF(BMAXNEX.GE.0.D0)THEN
1228 B1=BMINNEX/AM
1229 B2=MIN(BM*AM,BMAXNEX)/AM
1230 if(B1.gt.B2)stop'bmin > bmax in QGSJet'
1231 B=DSQRT(B1*B1+(B2*B2-B1*B1)*PSRAN(B10))
1232 BQGS=B*AM !ctp
1233 ENDIF
1234 IF(DEBUG.GE.2)WRITE (MONIOU,204)B*AM
1235 204 FORMAT(2X,'PSCONF: IMPACT PARAMETER FOR THE INTERACTION:',
1236 * E10.3,' FM')
1237
1238
1239
1240
1241 IF(IA(1).GT.1)CALL PSGEA(IA(1),XA,1)
1242
1243 NW=0
1244 LS=0
1245 NS=0
1246 NHP=0
1247 DO 101 IT = 1,NT
1248 IWT(IT) = 0
1249 101 CONTINUE
1250
1251
1252
1253 DO 14 IN=1,IA(1)
1254 IF(DEBUG.GE.2.AND.ICZ.EQ.2)WRITE (MONIOU,205)IN
1255 205 FORMAT(2X,'PSCONF: ',I2,'-TH PROJECTILE NUCLEON')
1256
1257
1258 IF(IA(1).NE.1.AND.PSRAN(B10).GT.CC(2))GOTO 12
1259
1260 X=XA(IN,1)+B
1261 Y=XA(IN,2)
1262
1263 IQS=0
1264 NW=NW+1
1265
1266
1267 DO 11 M=1,NT
1268
1269 Z=PSDR(X-XB(M,1),Y-XB(M,2))
1270
1271
1272 VV=2.D0*PSFAZ(Z,FSOFT,FHARD,FSHARD)
1273 EV=EXP(-VV)
1274
1275 EH=FHARD(1)+FHARD(2)+FHARD(3)
1276
1277 AKS=PSRAN(B10)
1278
1279
1280 IF(AKS.GT.1.D0-EV*(1.D0-2.D0*EH))GOTO 11
1281 IF(DEBUG.GE.2)WRITE (MONIOU,208)M
1282 208 FORMAT(2X,'PSCONF: INTERACTION WITH',I2,'-TH TARGET NUCLEON')
1283
1284 IWT(M) = 1
1285
1286
1287
1288 IQV=0
1289
1290
1291
1292 SUM=2.D0*EH*EV
1293
1294
1295 IF(AKS.LT.SUM)THEN
1296 AKS1=EH*PSRAN(B10)
1297 IF(AKS1.LT.FHARD(1))THEN
1298
1299 IF(LVA(NW).NE.0)GOTO 11
1300
1301 LVA(NW)=1
1302 IQV=1
1303 ELSEIF(AKS1.LT.FHARD(1)+FHARD(2))THEN
1304
1305 IF(LVB(M).NE.0)GOTO 11
1306
1307 LVB(M)=1
1308 IQV=2
1309 ELSE
1310
1311 IF(LVA(NW)+LVB(M).NE.0)GOTO 11
1312
1313 LVA(NW)=1
1314 LVB(M)=1
1315 IQV=3
1316 ENDIF
1317 N=1
1318
1319 LNH=1
1320 GOTO 22
1321 ENDIF
1322
1323
1324
1325 LNH=0
1326
1327 WH=2.D0*FSHARD/VV
1328
1329
1330
1331 DO 7 N=1,45
1332 EV=EV*VV/N
1333 SUM=SUM+EV
1334 7 IF(AKS.LT.SUM)GOTO 8
1335
1336
1337
1338 8 DO 9 I=1,N
1339 9 LNH=LNH+INT(WH+PSRAN(B10))
1340
1341
1342 AKS1=.5D0*PSRAN(B10)
1343
1344
1345 IF(AKS1.LT.EH)THEN
1346 IF(AKS1.LT.FHARD(1))THEN
1347 IF(LVA(NW).NE.0)GOTO 22
1348
1349 LVA(NW)=1
1350 IQV=1
1351 ELSEIF(AKS1.LT.FHARD(1)+FHARD(2))THEN
1352 IF(LVB(M).NE.0)GOTO 22
1353
1354 LVB(M)=1
1355 IQV=2
1356 ELSE
1357 IF(LVA(NW)+LVB(M).NE.0)GOTO 22
1358
1359 LVA(NW)=1
1360 LVB(M)=1
1361 IQV=3
1362 ENDIF
1363 N=N+1
1364 LNH=LNH+1
1365 ENDIF
1366
1367 22 IQS=1
1368 IF(LNH.NE.0)THEN
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380 N=N-LNH
1381 LHA(NW)=LHA(NW)+LNH
1382 LHB(M)=LHB(M)+LNH
1383 DO 10 I=1,LNH
1384 I1=NHP+I
1385 If (I1 .ge. 4000) then
1386 write(moniou,*)'psconf: I1 > 4000, index out of bounds'
1387 stop
1388 endif
1389 IF(I.EQ.1.AND.IQV.NE.0)THEN
1390 IQH(I1)=IQV
1391 ELSE
1392 IQH(I1)=0
1393 ENDIF
1394 IF(DEBUG.GE.2)WRITE (MONIOU,209)I1,NW,M,IQH(I1)
1395 209 FORMAT(2X,'PSCONF: ',I4,'-TH HARD BLOCK IS CONNECTED TO',1X,
1396 * I2,'-TH PROJECTILE NUCLEON (HADRON) AND'/4X,I2,
1397 * '-TH TARGET NUCLEON; TYPE OF THE SEMIHARD INTERACTION:',I1)
1398 ZH(I1)=Z
1399 IAH(I1)=NW
1400 10 IBH(I1)=M
1401
1402
1403 NHP=NHP+LNH
1404 ENDIF
1405
1406
1407 IF(N.GT.0)THEN
1408
1409 LS=LS+1
1410 IAS(LS)=NW
1411 IBS(LS)=M
1412 LQA(NW)=LQA(NW)+N
1413 LQB(M)=LQB(M)+N
1414 NQS(LS)=N
1415 IF(DEBUG.GE.2)WRITE (MONIOU,210)LS,NW,M,N
1416 210 FORMAT(2X,'PSCONF: ',I4,'-TH SOFT BLOCK IS CONNECTED TO',1X,
1417 * I2,'-TH PROJECTILE NUCLEON (HADRON) AND'/4X,I2,
1418 * '-TH TARGET NUCLEON; NUMBER OF POMERONS IN THE BLOCK NP=',
1419 * I2)
1420 ENDIF
1421 11 CONTINUE
1422
1423
1424 IF(IQS.NE.0)GOTO 14
1425
1426
1427
1428
1429
1430
1431
1432
1433 IF(JDIFR.EQ.1 .AND. IA(1).NE. 1
1434 * .AND.PSRAN(B10).LT.(1.D0-CC(ICZ))*PSV(X,Y,XB,NT))THEN
1435
1436 IF(IA(2).NE.1)THEN
1437
1438 ICT=INT(2.5+PSRAN(B10))
1439 ELSE
1440
1441 ICT=2
1442 ENDIF
1443
1444 IF(DEBUG.GE.2)WRITE(MONIOU,207)IN
1445 207 FORMAT(2X,I2,'-TH PROJECTILE NUCLEON DIFFRACTION')
1446 ICP0=INT(2.5+PSRAN(B10))
1447 WPI=WP0
1448 WMI=WM0
1449 IF(IA(2).EQ.1)THEN
1450 LQ=0
1451 ELSE
1452 LQ=1
1453 ENDIF
1454 CALL XXDPR(WPI,WMI,ICP0,ICT,LQ)
1455 GOTO 14
1456 ENDIF
1457
1458
1459
1460 NW=NW-1
1461 12 CONTINUE
1462
1463
1464
1465 NS=NS+1
1466 IF(NS.NE.IN)THEN
1467 DO 13 L=1,3
1468 13 XA(NS,L)=XA(IN,L)
1469 ENDIF
1470 14 CONTINUE
1471
1472
1473
1474
1475 IF(NS.EQ.IA(1))THEN
1476 IF(DEBUG.GE.3)WRITE (MONIOU,211)
1477 211 FORMAT(2X,'PSCONF: NO ONE NUCLEON (HADRON) INTERACTS - ',
1478 * 'REJECTION')
1479 GOTO 5
1480 ENDIF
1481
1482
1483 if(nhp.gt.1500)then
1484 WRITE (MONIOU,213)NHP
1485 213 FORMAT(2X,'PSCONF: TOO GREAT NUMBER OF HARD POMERONS: NHP=',
1486 * I5,' - REJECTION')
1487 GOTO 100
1488 endif
1489
1490 NWT = 0
1491
1492 DO 102 IT = 1,NT
1493 NWT = NWT + IWT(IT)
1494 102 CONTINUE
1495
1496
1497
1498 CALL XXFRAGM(NS,XA)
1499
1500
1501
1502 IF(NW.NE.0)CALL PSSHAR(LS,NHP,NW,NT)
1503 21 continue !so00
1504 IF(DEBUG.GE.3)WRITE (MONIOU,212)
1505 212 FORMAT(2X,'PSCONF - END')
1506 RETURN
1507 END
1508
1509
1510 SUBROUTINE QGSPSCS(C,S)
1511
1512
1513 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1514 INTEGER DEBUG
1515 COMMON /Q_AREA11/ B10
1516 COMMON /Q_AREA43/ MONIOU
1517 COMMON /Q_DEBUG/ DEBUG
1518 SAVE
1519 EXTERNAL PSRAN
1520
1521 IF(DEBUG.GE.2)WRITE (MONIOU,201)
1522 201 FORMAT(2X,'QGSPSCS - COS(FI) AND SIN(FI) ARE GENERATED',
1523 * ' (0<FI<2*PI)')
1524 1 S1=2.D0*PSRAN(B10)-1.D0
1525 S2=2.D0*PSRAN(B10)-1.D0
1526 S3=S1*S1+S2*S2
1527 IF(S3.GT.1.D0)GOTO 1
1528 S3=DSQRT(S3)
1529 C=S1/S3
1530 S=S2/S3
1531 IF(DEBUG.GE.3)WRITE (MONIOU,202)C,S
1532 202 FORMAT(2X,'QGSPSCS: C=',E10.3,2X,'S=',E10.3)
1533 RETURN
1534 END
1535
1536
1537 SUBROUTINE QGSPSDEFTR(S,EP,EY)
1538
1539
1540
1541 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1542 INTEGER DEBUG
1543 DIMENSION EY(3),EP(4)
1544 COMMON /Q_AREA43/ MONIOU
1545 COMMON /Q_DEBUG/ DEBUG
1546 SAVE
1547
1548 IF(DEBUG.GE.2)WRITE (MONIOU,201)EP,S
1549 201 FORMAT(2X,'QGSPSDEFTR - LORENTZ BOOST PARAMETERS:'/
1550 * 4X,'4-VECTOR EP=',4E10.3/4X,'4-VECTOR SQUARED S=',E10.3)
1551 DO 2 I=1,3
1552 IF(EP(I+1).EQ.0.D0)THEN
1553 EY(I)=1.D0
1554 ELSE
1555 WP=EP(1)+EP(I+1)
1556 WM=EP(1)-EP(I+1)
1557 IF(WM/WP.LT.1.D-8)THEN
1558 WW=S
1559 DO 1 L=1,3
1560 1 IF(L.NE.I)WW=WW+EP(L+1)**2
1561 WM=WW/WP
1562 ENDIF
1563 EY(I)=DSQRT(WM/WP)
1564 EP(1)=WP*EY(I)
1565 EP(I+1)=0.D0
1566 ENDIF
1567 2 CONTINUE
1568 IF(DEBUG.GE.3)WRITE (MONIOU,202)EY
1569 202 FORMAT(2X,'QGSPSDEFTR: LORENTZ BOOST PARAMETERS EY(I)=',2X
1570 * ,3E10.3)
1571 RETURN
1572 END
1573
1574
1575 SUBROUTINE QGSPSDEFROT(EP,S0X,C0X,S0,C0)
1576
1577
1578
1579 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1580 INTEGER DEBUG
1581 DIMENSION EP(4)
1582 COMMON /Q_AREA43/ MONIOU
1583 COMMON /Q_DEBUG/ DEBUG
1584 SAVE
1585
1586 IF(DEBUG.GE.2)WRITE (MONIOU,201)EP
1587 201 FORMAT(2X,'QGSPSDEFROT - SPACIAL ROTATION PARAMETERS'/4X,
1588 * '4-VECTOR EP=',2X,4(E10.3,1X))
1589
1590 PT2=EP(3)**2+EP(4)**2
1591 IF(PT2.NE.0.D0)THEN
1592 PT=DSQRT(PT2)
1593
1594
1595 C0X=EP(3)/PT
1596 S0X=EP(4)/PT
1597
1598 PL=DSQRT(PT2+EP(2)**2)
1599 S0=PT/PL
1600 C0=EP(2)/PL
1601 ELSE
1602 C0X=1.D0
1603 S0X=0.D0
1604 PL=ABS(EP(2))
1605 S0=0.D0
1606 C0=EP(2)/PL
1607 ENDIF
1608
1609 EP(2)=PL
1610 EP(3)=0.D0
1611 EP(4)=0.D0
1612 IF(DEBUG.GE.3)WRITE (MONIOU,202)S0X,C0X,S0,C0,EP
1613 202 FORMAT(2X,'QGSPSDEFROT: SPACIAL ROTATION PARAMETERS'/
1614 * 4X,'S0X=',E10.3,2X,'C0X=',E10.3,2X,'S0=',E10.3,2X,'C0=',E10.3/
1615 * 4X,'ROTATED 4-VECTOR EP=',4(E10.3,1X))
1616 RETURN
1617 END
1618
1619
1620 FUNCTION PSDR(X,Y)
1621
1622
1623 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1624 INTEGER DEBUG
1625 COMMON /Q_AREA7/ RP
1626 COMMON /Q_AREA43/ MONIOU
1627 COMMON /Q_DEBUG/ DEBUG
1628 SAVE
1629
1630 IF(DEBUG.GE.2)WRITE (MONIOU,201)X,Y
1631 201 FORMAT(2X,'PSDR: NUCLEON COORDINATES - X=',E10.3,2X,'Y=',E10.3)
1632 PSDR=EXP(-(X*X+Y*Y)/RP)
1633 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSDR
1634 202 FORMAT(2X,'PSDR=',E10.3)
1635 RETURN
1636 END
1637
1638
1639 FUNCTION QGSPSFAP(X,J,L)
1640
1641
1642
1643
1644
1645 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1646 INTEGER DEBUG
1647 COMMON /Q_AREA43/ MONIOU
1648 COMMON /Q_DEBUG/ DEBUG
1649 SAVE
1650
1651 IF(DEBUG.GE.2)WRITE (MONIOU,201)X,J,L
1652 201 FORMAT(2X,'QGSPSFAP - ALTARELLI-PARISI FUNCTION:',2X,
1653 * 'X=',E10.3,2X,'J=',I1,2X,'L=',I1)
1654 IF(J.EQ.0)THEN
1655 IF(L.EQ.0)THEN
1656 QGSPSFAP=((1.D0-X)/X+X/(1.D0-X)+X*(1.D0-X))*6.d0
1657 ELSE
1658 QGSPSFAP=(X**2+(1.D0-X)**2)*3.d0
1659 ENDIF
1660 ELSE
1661 IF(l.EQ.0)THEN
1662 QGSPSFAP=(1.D0+(1.D0-X)**2)/X/.75D0
1663 ELSE
1664 QGSPSFAP=(X**2+1.D0)/(1.D0-X)/.75D0
1665 ENDIF
1666 ENDIF
1667 IF(DEBUG.GE.3)WRITE (MONIOU,202)QGSPSFAP
1668 202 FORMAT(2X,'QGSPSFAP=',E10.3)
1669 RETURN
1670 END
1671
1672
1673 FUNCTION PSFAZ(Z,FSOFT,FHARD,FSHARD)
1674
1675
1676
1677
1678
1679
1680
1681 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1682 INTEGER DEBUG
1683 DIMENSION FHARD(3)
1684 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALF,RR,SH,DELH
1685 COMMON /Q_AREA22/ SJV,FJS(5,3)
1686 COMMON /Q_AREA43/ MONIOU
1687 COMMON /Q_DEBUG/ DEBUG
1688 SAVE
1689
1690 IF(DEBUG.GE.2)WRITE (MONIOU,201)Z
1691 201 FORMAT(2X,'PSFAZ - HADRON-NUCLEON (NUCLEON-NUCLEON)',
1692 * ' INTERACTION EIKONAL; Z=',E10.3)
1693 FSOFT=FS*Z
1694 FHARD(3)=SJV*Z**(RS/RS0)
1695
1696 JZ=INT(5.D0*Z)
1697 IF(JZ.GT.3)JZ=3
1698 WZ=5.D0*Z-JZ
1699
1700 DO 1 I=1,3
1701 IF(JZ.EQ.0)THEN
1702 FSR=(EXP(FJS(1,I))*WZ+(EXP(FJS(2,I))-2.D0*
1703 * EXP(FJS(1,I)))*WZ*(WZ-1.D0)*.5D0)*Z
1704 ELSE
1705 FSR=EXP(FJS(JZ,I)+(FJS(JZ+1,I)-FJS(JZ,I))*WZ
1706 * +(FJS(JZ+2,I)+FJS(JZ,I)-2.D0*FJS(JZ+1,I))
1707 * *WZ*(WZ-1.D0)*.5D0)*Z
1708 ENDIF
1709 IF(I.NE.1)THEN
1710 FHARD(I-1)=FSR
1711 ELSE
1712 FSHARD=FSR
1713 ENDIF
1714 1 CONTINUE
1715
1716 PSFAZ=FSOFT+FSHARD
1717 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSFAZ,FSOFT,FSHARD,FHARD
1718 202 FORMAT(2X,'PSFAZ=',E10.3,2X,'FSOFT=',E10.3,2X,'FSHARD=',E10.3/4X
1719 * ,'FHARD=',3E10.3)
1720 RETURN
1721 END
1722
1723
1724 FUNCTION PSFBORN(S,T,IQ1,IQ2)
1725
1726
1727
1728
1729
1730
1731 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1732 INTEGER DEBUG
1733 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
1734 COMMON /Q_AREA43/ MONIOU
1735 COMMON /Q_DEBUG/ DEBUG
1736 SAVE
1737
1738 IF(DEBUG.GE.2)WRITE (MONIOU,201)S,T,IQ1,IQ2
1739 201 FORMAT(2X,'PSFBORN - HARD SCATTERING MATRIX ELEMENT SQUARED:'/
1740 * 4X,'S=',E10.3,2X,'|T|=',E10.3,2X,'IQ1=',I2,2X,'IQ2=',I2)
1741 U=S-T
1742 IF(IQ1.EQ.0.AND.IQ2.EQ.0)THEN
1743
1744 PSFBORN=(3.D0-T*U/S**2+S*U/T**2+S*T/U**2)*4.5D0
1745 ELSEIF(IQ1*IQ2.EQ.0)THEN
1746
1747 PSFBORN=(S**2+U**2)/T**2+(S/U+U/S)/2.25D0
1748 ELSEIF(IQ1.EQ.IQ2)THEN
1749
1750 PSFBORN=((S**2+U**2)/T**2+(S**2+T**2)/U**2)/2.25D0
1751 * -S**2/T/U/3.375D0
1752 ELSEIF(IQ1+IQ2.EQ.0)THEN
1753
1754 PSFBORN=((S**2+U**2)/T**2+(U**2+T**2)/S**2)/2.25D0
1755 * -U**2/T/S/3.375D0
1756 ELSE
1757
1758 PSFBORN=(S**2+U**2)/T**2/2.25D0
1759 ENDIF
1760 IF(DEBUG.GE.2)WRITE (MONIOU,202)PSFBORN
1761 202 FORMAT(2X,'PSFBORN=',E10.3)
1762 RETURN
1763 END
1764
1765
1766 FUNCTION QGSPSFSH(S,Z,ICZ,IQQ)
1767
1768
1769
1770
1771
1772
1773 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1774 INTEGER DEBUG
1775 COMMON /Q_AREA6/ PI,BM,AM
1776 COMMON /Q_AREA15/ FP(5),RQ(5),CD(5)
1777 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALF,RR,SH,DELH
1778 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
1779 COMMON /Q_AREA19/ AHL(5)
1780 COMMON /Q_AREA25/ AHV(5)
1781 COMMON /Q_AREA27/ FP0(5)
1782 COMMON /Q_AR13/ X1(7),A1(7)
1783 COMMON /Q_AREA43/ MONIOU
1784 COMMON /Q_DEBUG/ DEBUG
1785 SAVE
1786
1787 IF(DEBUG.GE.2)WRITE (MONIOU,201)S,Z,IQQ,ICZ
1788 201 FORMAT(2X,'QGSPSFSH - SEMIHARD INTERACTION EIKONAL:'/
1789 * 4X,'S=',E10.3,2X,'Z=',E10.3,2X,'IQQ=',I1,2X,'ICZ=',I1)
1790 XMIN=4.D0*QT0/S
1791 XMIN=XMIN**(DELH-DEL)
1792 QGSPSFSH=0.D0
1793 IF(IQQ.EQ.1)THEN
1794 ICV=ICZ
1795 ICQ=2
1796 ELSEIF(IQQ.EQ.2)THEN
1797 ICV=2
1798 ICQ=ICZ
1799 ENDIF
1800 IQ=(IQQ+1)/2
1801
1802
1803 DO 3 I=1,7
1804 DO 3 M=1,2
1805 Z1=(.5D0*(1.D0+XMIN-(2*M-3)*X1(I)*(1.D0-XMIN)))**(1.D0/
1806 * (DELH-DEL))
1807
1808
1809
1810
1811 CALL PSJINT0(Z1*S,SJ,SJB,IQ,0)
1812
1813
1814
1815 GY=2.D0*SH*PSGINT((SJ-SJB)/SH*.5D0)+SJB
1816 IF(DEBUG.GE.3)WRITE (MONIOU,203)Z1*S,GY
1817 203 FORMAT(2X,'QGSPSFSH:',2X,'S_HARD=',E10.3,2X,'SIGMA_HARD=',E10.3)
1818
1819 IF(IQQ.EQ.0)THEN
1820 ST2=0.D0
1821 DO 1 J=1,7
1822 DO 1 K=1,2
1823 XX=.5D0*(1.D0+X1(J)*(2*K-3))
1824 1 ST2=ST2+A1(J)*QGSPSFTILD(Z1**XX,ICZ)*
1825 * QGSPSFTILD(Z1**(1.D0-XX),2)
1826
1827 RH=RS0-ALF*DLOG(Z1)
1828 QGSPSFSH=QGSPSFSH-A1(I)*DLOG(Z1)*GY/Z1**DELH*Z**(RS/RH)/RH*ST2
1829 ELSE
1830
1831 ST2=0.D0
1832 DO 2 J=1,7
1833 DO 2 K=1,2
1834 XX=.5D0*(1.D0+X1(J)*(2*K-3))
1835 XAM=Z1**(DEL+.5D0)
1836 XA=(XAM+(1.D0-XAM)*XX)**(1.D0/(DEL+.5D0))
1837 RH=RS0+ALF*DLOG(XA/Z1)
1838 2 ST2=ST2+A1(J)*(1.D0-XA)**AHV(ICV)*Z**(RS/RH)/RH*
1839 * QGSPSFTILD(Z1/XA,ICQ)
1840 ST2=ST2*(1.D0-XAM)
1841
1842 QGSPSFSH=QGSPSFSH+A1(I)*GY/Z1**DELH*ST2
1843 ENDIF
1844 3 CONTINUE
1845
1846 IF(IQQ.EQ.0)THEN
1847 QGSPSFSH=QGSPSFSH*.125D0*RR*(1.D0-XMIN)/(DELH-DEL)*FP0(ICZ)
1848 * *FP0(2)
1849 * *CD(ICZ)
1850 ELSE
1851 QGSPSFSH=QGSPSFSH*DSQRT(RR)/16.D0*FP0(ICQ)*(1.D0-XMIN)
1852 * /(DELH-DEL)/(DEL+.5D0)*GAMFUN(AHV(ICV)+1.5D0)
1853 * /GAMFUN(AHV(ICV)+1.D0)/PI*CD(ICZ)
1854 IF(ICZ.EQ.2.OR.IQQ.EQ.2)THEN
1855 QGSPSFSH=QGSPSFSH*3.D0
1856 ELSEIF((ICZ-1)*(ICZ-3)*(ICZ-5).EQ.0)THEN
1857 QGSPSFSH=QGSPSFSH*2.D0
1858 ENDIF
1859 ENDIF
1860 IF(DEBUG.GE.3)WRITE (MONIOU,202)QGSPSFSH
1861 202 FORMAT(2X,'QGSPSFSH=',E10.3)
1862 RETURN
1863 END
1864
1865
1866 FUNCTION QGSPSFTILD(Z,ICZ)
1867
1868
1869
1870
1871
1872 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1873 INTEGER DEBUG
1874 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALFP,RR,SH,DELH
1875 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
1876 COMMON /Q_AREA19/ AHL(5)
1877 COMMON /Q_AR13/ X1(7),A1(7)
1878 COMMON /Q_AREA43/ MONIOU
1879 COMMON /Q_DEBUG/ DEBUG
1880 SAVE
1881
1882 IF(DEBUG.GE.2)WRITE (MONIOU,201)Z,ICZ
1883 201 FORMAT(2X,'QGSPSFTILD:',2X,'Z=',E10.3,2X,'ICZ=',I1)
1884 QGSPSFTILD=0.
1885 DO 1 I=1,7
1886 DO 1 M=1,2
1887 XB=1.D0-(1.D0-Z)*(.5D0*(1.D0+(2*M-3)*X1(I)))**(1.D0/
1888 * (AHL(ICZ)+1.D0))
1889 1 QGSPSFTILD=QGSPSFTILD+A1(I)*XB**DEL*(1.D0-Z/XB)**BET
1890 QGSPSFTILD=QGSPSFTILD*.5D0*(1.D0-Z)**(AHL(ICZ)+1.D0)
1891 * /(AHL(ICZ)+1.D0)
1892 IF(DEBUG.GE.3)WRITE (MONIOU,202)QGSPSFTILD
1893 202 FORMAT(2X,'QGSPSFTILD=',E10.3)
1894 RETURN
1895 END
1896
1897
1898 SUBROUTINE PSGEA(IA,XA,JJ)
1899
1900
1901
1902 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1903 INTEGER DEBUG
1904
1905 DIMENSION XA(210,3)
1906 COMMON /Q_AREA5/ RD(2),CA1(2),CA2(2),CA3(2)
1907 COMMON /Q_AREA11/ B10
1908 COMMON /Q_AREA43/ MONIOU
1909 COMMON /Q_DEBUG/ DEBUG
1910 SAVE
1911 EXTERNAL PSRAN
1912
1913 IF(DEBUG.GE.2)WRITE (MONIOU,201)JJ,IA
1914 201 FORMAT(2X,'PSGEA - CONFIGURATION OF THE NUCLEUS ',I1,';',2X,
1915 * 'COORDINATES FOR ',I2,' NUCLEONS')
1916
1917 IF(IA.GE.10)THEN !this line had been changed!!!!!!! dh 8/10/98
1918
1919 DO 7 I=1,IA
1920 1 ZUK=PSRAN(B10)*CA1(JJ)-1.D0
1921 IF(ZUK)2,2,3
1922 2 TT=RD(JJ)*(PSRAN(B10)**.3333D0-1.D0)
1923 GOTO 6
1924 3 IF(ZUK.GT.CA2(JJ))GOTO 4
1925 TT=-DLOG(PSRAN(B10))
1926 GOTO 6
1927 4 IF(ZUK.GT.CA3(JJ))GOTO 5
1928 TT=-DLOG(PSRAN(B10))-DLOG(PSRAN(B10))
1929 GOTO 6
1930 5 TT=-DLOG(PSRAN(B10))-DLOG(PSRAN(B10))-DLOG(PSRAN(B10))
1931 6 IF(PSRAN(B10).GT.1.D0/(1.D0+EXP(-ABS(TT))))GOTO 1
1932 RIM=TT+RD(JJ)
1933 Z=RIM*(2.D0*PSRAN(B10)-1.D0)
1934 RIM=DSQRT(RIM*RIM-Z*Z)
1935 XA(I,3)=Z
1936 CALL QGSPSCS(C,S)
1937 XA(I,1)=RIM*C
1938 7 XA(I,2)=RIM*S
1939 ELSE
1940
1941 DO 9 L=1,3
1942 SUMM=0.D0
1943 DO 8 I=1,IA-1
1944 J=IA-I
1945 AKS=RD(JJ)*(PSRAN(B10)+PSRAN(B10)+PSRAN(B10)-1.5D0)
1946 K=J+1
1947 XA(K,L)=SUMM-AKS*SQRT(FLOAT(J)/K)
1948 8 SUMM=SUMM+AKS/SQRT(FLOAT(J*K))
1949 9 XA(1,L)=SUMM
1950 ENDIF
1951 IF(DEBUG.GE.3)THEN
1952 WRITE (MONIOU,203)
1953 DO 206 I=1,IA
1954 206 WRITE (MONIOU,204)I,(XA(I,L),L=1,3)
1955 WRITE (MONIOU,202)
1956 ENDIF
1957 202 FORMAT(2X,'PSGEA - END')
1958 203 FORMAT(2X,'PSGEA: POSITIONS OF THE NUCLEONS')
1959 204 FORMAT(2X,'PSGEA: ',I2,' - ',3(E10.3,1X))
1960 RETURN
1961 END
1962
1963
1964 FUNCTION PSGINT(Z)
1965
1966
1967
1968 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1969 INTEGER DEBUG
1970 COMMON /Q_AR13/ X1(7),A1(7)
1971 COMMON /Q_AREA43/ MONIOU
1972 COMMON /Q_DEBUG/ DEBUG
1973 SAVE
1974
1975 F(Z,X)=(1.-EXP(-.5*Z*(1.+X)))/(1.+X)
1976
1977 IF(DEBUG.GE.2)WRITE (MONIOU,201)Z
1978 201 FORMAT(2X,'PSGINT:',2X,'Z=',E10.3)
1979 PSGINT=0.
1980 DO 5 I=1,7
1981 5 PSGINT=PSGINT+A1(I)*(F(Z,X1(I))+F(Z,-X1(I)))
1982 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSGINT
1983 202 FORMAT(2X,'PSGINT=',E10.3)
1984 RETURN
1985 END
1986
1987
1988 FUNCTION QGSPSHARD(S,ICZ)
1989
1990
1991
1992
1993 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
1994 INTEGER DEBUG
1995 COMMON /Q_AR13/ X1(7),A1(7)
1996 COMMON /Q_AREA6/ PI,BM,AM
1997 COMMON /Q_AREA15/ FP(5),RQ(5),CD(5)
1998 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALF,RR,SH,DELH
1999 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
2000 COMMON /Q_AREA19/ AHL(5)
2001 COMMON /Q_AREA25/ AHV(5)
2002 COMMON /Q_AREA43/ MONIOU
2003 COMMON /Q_DEBUG/ DEBUG
2004 SAVE
2005
2006 IF(DEBUG.GE.2)WRITE (MONIOU,201)S,ICZ
2007 201 FORMAT(2X,'QGSPSHARD - HARD QUARK-QUARK INTERACTION CROSS',
2008 * ' SECTION:',
2009 * 2X,'S=',E10.3,2X,'ICZ=',I1)
2010 XMIN=4.D0*QT0/S
2011 XMIN=XMIN**(DELH+.5D0)
2012 QGSPSHARD=0.D0
2013
2014
2015 DO 2 I=1,7
2016 DO 2 M=1,2
2017 Z1=(.5D0*(1.D0+XMIN-(2*M-3)*X1(I)*(1.D0-XMIN)))**(1.D0/
2018 * (DELH+.5D0))
2019
2020 ST2=0.D0
2021 DO 1 J=1,7
2022 DO 1 K=1,2
2023 XX=.5D0*(1.D0+X1(J)*(2*K-3))
2024 ST2=ST2+A1(J)*(1.D0-Z1**XX)**AHV(ICZ)*
2025 * (1.D0-Z1**(1.D0-XX))**AHV(2)
2026 1 CONTINUE
2027
2028
2029
2030
2031
2032 CALL PSJINT0(Z1*S,SJ,SJB,1,1)
2033
2034
2035
2036 GY=2.D0*SH*PSGINT((SJ-SJB)/SH*.5D0)+SJB
2037
2038 IF(DEBUG.GE.3)WRITE (MONIOU,203)Z1*S,GY
2039 203 FORMAT(2X,'QGSPSHARD:',2X,'S_HARD=',E10.3,2X,'SIGMA_HARD=',E10.3)
2040 QGSPSHARD=QGSPSHARD-A1(I)*DLOG(Z1)*GY/Z1**DELH*ST2
2041 2 CONTINUE
2042
2043 QGSPSHARD=QGSPSHARD*(1.D0-XMIN)/(.5D0+DELH)*.25D0
2044 QGSPSHARD=QGSPSHARD/(GAMFUN(AHV(ICZ)+1.D0)*GAMFUN(AHV(2)+1.D0)
2045 * *PI)*GAMFUN(AHV(ICZ)+1.5D0)*GAMFUN(AHV(2)+1.5D0)
2046
2047 IF(ICZ.EQ.2)THEN
2048 QGSPSHARD=QGSPSHARD*9.D0
2049 ELSEIF((ICZ-1)*(ICZ-3)*(ICZ-5).EQ.0)THEN
2050 QGSPSHARD=QGSPSHARD*6.D0
2051 ELSE
2052 QGSPSHARD=QGSPSHARD*3.D0
2053 ENDIF
2054
2055
2056
2057
2058 QGSPSHARD=QGSPSHARD/(8.D0*PI*RS0)*CD(ICZ)
2059 IF(DEBUG.GE.2)WRITE (MONIOU,202)QGSPSHARD
2060 202 FORMAT(2X,'QGSPSHARD=',E10.3)
2061 RETURN
2062 END
2063
2064
2065 SUBROUTINE PSHOT(WP0,WM0,Z,IPC,EPC,IZP,IZT,ICZ,IQQ)
2066
2067
2068
2069
2070
2071
2072
2073 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
2074 INTEGER DEBUG
2075 CHARACTER*2 TYQ
2076 DIMENSION EPT(4),EP3(4),EPJ(4),EPJ1(4),EY(3),
2077 * QMIN(2),WP(2),IQC(2),IQP(2),
2078 * IPC(2,2),EPC(8,2),IQJ(2),EQJ(4,2),IPQ(2,2),EPQ(8,2),
2079 * ebal(4),
2080 * QV1(30,50),ZV1(30,50),QM1(30,50),IQV1(30,50),
2081 * LDAU1(30,49),LPAR1(30,50),
2082 * QV2(30,50),ZV2(30,50),QM2(30,50),IQV2(30,50),
2083 * LDAU2(30,49),LPAR2(30,50)!,EP(4,2),EPT0(4)
2084 COMMON /Q_AREA6/ PI,BM,AMMM
2085 COMMON /Q_AREA8/ WWM,BE(4),DC(5),DETA,ALMPT
2086 COMMON /Q_AREA10/ STMASS,AM(7)
2087 COMMON /Q_AREA11/ B10
2088 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALF,RR,SH,DELH
2089 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
2090 COMMON /Q_AREA42/ TYQ(15)
2091 COMMON /Q_AREA43/ MONIOU
2092 COMMON /Q_AREA46/ EPJET(4,2,15000),IPJET(2,15000)
2093 COMMON /Q_AREA47/ NJTOT
2094 COMMON /Q_DEBUG/ DEBUG
2095 SAVE
2096 EXTERNAL PSRAN
2097
2098 IF(DEBUG.GE.1)WRITE (MONIOU,201)IQQ,WP0,WM0
2099 201 FORMAT(2X,'PSHOT - SEMIHARD INTERACTION SIMULATION:'/
2100 * 4X,'TYPE OF THE INTERACTION:',I2/
2101 * 4X,'INITIAL LIGHT CONE MOMENTA:',2E10.3)
2102
2103 NJTOT0=NJTOT
2104 IZP0=IZP
2105 IZT0=IZT
2106
2107 301 S=WP0*WM0
2108 NJTOT=NJTOT0
2109 IZP=IZP0
2110 IZT=IZT0
2111
2112 IF(IQQ.EQ.3)THEN
2113
2114 WPI=WP0
2115 WMI=WM0
2116
2117
2118
2119
2120
2121 CALL PSJINT0(S,SJ,SJB,1,1)
2122
2123
2124
2125 GY=2.D0*SH*PSGINT((SJ-SJB)/SH*.5D0)+SJB
2126
2127 ELSE
2128
2129
2130
2131
2132
2133 XMIN=4.D0*(QT0+AMJ0)/S
2134 XMIN1=XMIN**(DELH-DEL)
2135
2136
2137 IQ=(IQQ+1)/2
2138
2139
2140 IF(IQQ.EQ.0)THEN
2141 GB0=-DLOG(XMIN)*(1.D0-DSQRT(XMIN))**(2.D0*BET)
2142 ELSE
2143 GB0=(1.D0-XMIN)**BET
2144 ENDIF
2145
2146
2147
2148
2149
2150 CALL PSJINT0(S,SJ,SJB,IQ,0)
2151
2152
2153
2154 GY0=2.D0*SH*PSGINT((SJ-SJB)/SH*.5D0)+SJB
2155 GB0=GB0*GY0/S**DELH/RS0*Z
2156
2157
2158
2159
2160
2161
2162
2163
2164
2165
2166
2167
2168
2169 1 ZPM=(XMIN1+PSRAN(B10)*(1.D0-XMIN1))**(1.D0/(DELH-DEL))
2170
2171
2172
2173
2174 CALL PSJINT0(ZPM*S,SJ,SJB,IQ,0)
2175 YJ=DLOG(ZPM*S)
2176
2177 RH=RS0-ALF*DLOG(ZPM)
2178
2179 IF(IQQ.EQ.0)THEN
2180
2181 XP=ZPM**PSRAN(B10)
2182 XM=ZPM/XP
2183
2184
2185 GBYJ=-DLOG(ZPM)*((1.-XP)*(1.-XM))**BET
2186
2187 WPI=WP0*XP
2188 WMI=WM0*XM
2189 ELSE
2190 IF(IQQ.EQ.1)THEN
2191 WPI=WP0
2192 WMI=WM0*ZPM
2193 ELSE
2194 WPI=WP0*ZPM
2195 WMI=WM0
2196 ENDIF
2197 GBYJ=(1.D0-ZPM)**BET
2198 ENDIF
2199
2200
2201
2202
2203 GY=2.D0*SH*PSGINT((SJ-SJB)/SH*.5D0)+SJB
2204
2205
2206
2207
2208 GBYJ=GBYJ*GY*EXP(-DELH*YJ)/GB0*Z**(RS/RH)/RH
2209 IF(PSRAN(B10).GT.GBYJ)GOTO 1
2210 ENDIF
2211
2212 S=WPI*WMI
2213
2214 IF(DEBUG.GE.2)WRITE (MONIOU,203)S
2215 203 FORMAT(2X,'PSHOT: MASS SQUARED FOR THE HARD PARTON-PARTON',
2216 * ' INTERACTION:',E10.3)
2217
2218
2219
2220
2221
2222 DO 302 I=1,8
2223 DO 302 M=1,2
2224 302 EPC(I,M)=0.D0
2225
2226 IF((IQQ-1)*(IQQ-3).EQ.0)THEN
2227 CALL PSVDEF(IZP,IC1,ICZ)
2228 IQC(1)=IC1
2229 IPC(1,1)=0
2230 IPC(2,1)=0
2231 ELSE
2232 IQC(1)=0
2233 IPC(1,1)=-INT(2.D0*PSRAN(B10)+1.D0)
2234 IPC(2,1)=-IPC(1,1)
2235 WP1=WP0-WPI
2236 WP2=WP1*PSRAN(B10)
2237 WP1=WP1-WP2
2238 EPC(1,1)=.5D0*WP1
2239 EPC(2,1)=EPC(1,1)
2240 EPC(5,1)=.5D0*WP2
2241 EPC(6,1)=EPC(5,1)
2242 ENDIF
2243
2244 IF((IQQ-2)*(IQQ-3).EQ.0)THEN
2245 CALL PSVDEF(IZT,IC1,2)
2246 IQC(2)=IC1
2247 IPC(1,2)=0
2248 IPC(2,2)=0
2249 ELSE
2250 IQC(2)=0
2251 IPC(1,2)=-INT(2.D0*PSRAN(B10)+1.D0)
2252 IPC(2,2)=-IPC(1,2)
2253 WM1=WM0-WMI
2254 WM2=WM1*PSRAN(B10)
2255 WM1=WM1-WM2
2256 EPC(1,2)=.5D0*WM1
2257 EPC(2,2)=-EPC(1,2)
2258 EPC(5,2)=.5D0*WM2
2259 EPC(6,2)=-EPC(5,2)
2260 ENDIF
2261
2262 EPT(1)=.5D0*(WPI+WMI)
2263 EPT(2)=.5D0*(WPI-WMI)
2264 EPT(3)=0.D0
2265 EPT(4)=0.D0
2266
2267 QMIN(1)=QT0
2268 QMIN(2)=QT0
2269 DO 303 L=1,2
2270 DO 303 M=1,2
2271 IPQ(L,M)=IPC(L,M)
2272 DO 303 I=1,4
2273 303 EPQ(I+4*(L-1),M)=EPC(I+4*(L-1),M)
2274
2275 QMINN=MAX(QMIN(1),QMIN(2))
2276 SI=QGSPSNORM(EPT)
2277
2278 5 CONTINUE
2279
2280
2281 IF(DEBUG.GE.2)WRITE (MONIOU,208)ILAD, SI,IQC,EPT
2282 208 FORMAT(2X,'PSHOT: ',I2,'-TH HARD LADDER;',
2283 * ' MASS SQUARED FOR THE LADDDER:',E10.3/
2284 * 4X,'LADDER END FLAVORS:',2I3/4X,
2285 * 'LADDER 4-MOMENTUM: ',4E10.3)
2286
2287 ebal(1)=.5*(wp0+wm0)-ept(1)
2288 ebal(2)=.5*(wp0-wm0)-ept(2)
2289 ebal(3)=0.d0-ept(3)
2290 ebal(4)=0.d0-ept(4)
2291 do 503 l=1,4
2292 do 501 m=1,2
2293 ebal(l)=ebal(l)-epq(l,m)
2294 501 if(iqc(m).eq.0) ebal(l)=ebal(l)-epq(l+4,m)
2295 if(njtot.ne.0)then
2296 do 502 i=1,njtot
2297 do 502 m=1,2
2298 502 ebal(l)=ebal(l)-epjet(l,m,i)
2299 endif
2300 503 continue
2301
2302
2303 PT2=EPT(3)**2+EPT(4)**2
2304 PT=DSQRT(PT2)
2305 WW=SI+PT2
2306 SWW=DSQRT(WW)
2307
2308 IQP(1)=MIN(1,IABS(IQC(1)))
2309 IQP(2)=MIN(1,IABS(IQC(2)))
2310
2311
2312 WP(1)=EPT(1)+EPT(2)
2313 WP(2)=EPT(1)-EPT(2)
2314
2315 S2MIN=MAX(QMINN,4.D0*(QT0+AMJ0))
2316
2317
2318 WWMIN=(S2MIN+(PT-DSQRT(QT0))**2+(QT0+AMJ0)*(DSQRT(S2MIN/QT0)-
2319 * 1.D0))/(1.D0-DSQRT(QT0/S2MIN))
2320
2321
2322
2323 SJ=PSJINT(QMIN(1),QMIN(2),SI,IQP(1)+1,IQP(2)+1)
2324 SJB=QGSPSBINT(QMINN,SI,IQP(1)+1,IQP(2)+1)
2325
2326 IF(DEBUG.GE.2)WRITE (MONIOU,251)S2MIN,WWMIN,SJ,SJB
2327 251 FORMAT(2X,'PSHOT: KINEMATICAL BOUNDS S2MIN=',E10.3,
2328 * 2X,'WWMIN=',E10.3/4X,'JET CROSS SETION SJ=',E10.3,
2329 * 2X,'BORN CROSS SECTION SJB=',E10.3)
2330
2331 IF(PSRAN(B10).LT.SJB/SJ.
2332 * OR.WW.LT.1.2D0*WWMIN)GOTO 12
2333
2334 IF((SJ-SJB)/SJ.GT..1D0)THEN
2335 SJ1=PSJINT1(QMIN(1),QMIN(2),SI,IQP(1)+1,IQP(2)+1)
2336 SJ2=PSJINT1(QMIN(2),QMIN(1),SI,IQP(2)+1,IQP(1)+1)
2337 DSJ=(SJ2-SJ1)/(SJ-SJB)*.5D0
2338 ELSE
2339 DSJ=0.D0
2340 ENDIF
2341
2342
2343 JJ=INT(1.5D0+DSJ+PSRAN(B10))
2344
2345 AQ=-(SI+AMJ0+2.D0*PT*DSQRT(QT0))/WW
2346 BQ=(QT0+AMJ0)/WW
2347 CQ=QT0/WW
2348 PQ=-AQ**2/3.D0+BQ
2349 QQ=AQ**3/13.5D0-AQ*BQ/3.D0+CQ
2350 PQ=DSQRT(-PQ/3.D0)
2351 COSQ=-.5D0*QQ/PQ**3
2352 FQ=ATAN(1.D0/COSQ**2-1.D0)
2353 IF(COSQ.LT.0.D0)FQ=PI-FQ
2354 FQ=FQ/3.D0
2355
2356
2357
2358 XMIN=1.D0+AQ/3.D0-2.D0*PQ*COS(FQ)
2359 XMAX=1.D0+AQ/3.D0-PQ*(DSQRT(3.D0)*SIN(FQ)-COS(FQ))
2360
2361
2362
2363 QQMAX=QT0/(1.D0-XMAX)**2
2364 QQMIN=QT0/(1.D0-XMIN)**2
2365
2366 IF(QQMIN.LT.S2MIN)THEN
2367 XMM=(SI-S2MIN+AMJ0+2.D0*PT*DSQRT(QT0))/WW*.5D0
2368 XMIN=1.D0-XMM-DSQRT(XMM*XMM-(QT0+AMJ0)/WW)
2369 QQMIN=QT0/(1.D0-XMIN)**2
2370
2371 IF(QQMIN.LT.QMIN(JJ))THEN
2372 QQMIN=QMIN(JJ)
2373 XMM1=WW-2.D0*PT*DSQRT(QQMIN)+QQMIN
2374 XMM=(SI-S2MIN+AMJ0)/XMM1*.5D0
2375 XMIN=1.D0-XMM-DSQRT(XMM*XMM-AMJ0/XMM1)
2376 ENDIF
2377 ENDIF
2378
2379
2380 XM0=MAX(.5D0,1.D0-DSQRT(QT0/QMIN(JJ)))
2381 IF(XM0.GT..95D0*XMAX.OR.XM0.LT.1.05D0*XMIN)
2382 * XM0=.5D0*(XMAX+XMIN)
2383 QM0=QT0/(1.D0-XM0)**2
2384 S2MAX=XM0*WW
2385
2386 SJ0=PSJINT(QM0,QMIN(3-JJ),S2MAX,1,IQP(3-JJ)+1)*
2387 * QGSPSFAP(XM0,IQP(JJ),0)+
2388 * PSJINT(QM0,QMIN(3-JJ),S2MAX,2,IQP(3-JJ)+1)
2389 * *QGSPSFAP(XM0,IQP(JJ),1)
2390
2391 GB0=SJ0*QM0/QLOG*QGSPSUDS(QM0,IQP(JJ))*1.5D0
2392 IF(XM0.LE..5D0)THEN
2393 GB0=GB0*XM0**(1.D0-DELH)
2394 ELSE
2395 GB0=GB0*(1.D0-XM0)*2.D0**DELH
2396 ENDIF
2397
2398 XMIN2=MAX(.5D0,XMIN)
2399 XMIN1=XMIN**DELH
2400 XMAX1=MIN(XMAX,.5D0)**DELH
2401 IF(XMIN.GE..5D0)THEN
2402 DJL=1.D0
2403 ELSEIF(XMAX.LT..5D0)THEN
2404 DJL=0.D0
2405 ELSE
2406 DJL=1.D0/(1.D0+((2.D0*XMIN)**DELH-1.D0)/DELH/
2407 * DLOG(2.D0*(1.D0-XMAX)))
2408 ENDIF
2409
2410 7 CONTINUE
2411
2412
2413 IF(PSRAN(B10).GT.DJL)THEN
2414 X=(XMIN1+PSRAN(B10)*(XMAX1-XMIN1))**(1.D0/DELH)
2415 ELSE
2416 X=1.D0-(1.D0-XMIN2)*((1.D0-XMAX)/(1.D0-XMIN2))**PSRAN(B10)
2417 ENDIF
2418
2419
2420
2421
2422 QQ=QQMIN/(1.D0+PSRAN(B10)*(QQMIN/QQMAX-1.D0))
2423
2424 IF(DEBUG.GE.2)WRITE (MONIOU,253)QQ,QQMIN,QQMAX
2425 253 FORMAT(2X,'PSHOT: QQ=',E10.3,2X,'QQMIN=',E10.3,2X,
2426 * 'QQMAX=',E10.3)
2427
2428 QT2=QQ*(1.D0-X)**2
2429 IF(QT2.LT.QT0)GOTO 7
2430
2431 IF(QQ.GT.QMINN)THEN
2432 QMIN2=QQ
2433 ELSE
2434 QMIN2=QMINN
2435 ENDIF
2436
2437 QT=DSQRT(QT2)
2438 CALL QGSPSCS(CCOS,SSIN)
2439
2440 EP3(3)=QT*CCOS
2441 EP3(4)=QT*SSIN
2442 PT2=(EPT(3)-EP3(3))**2+(EPT(4)-EP3(4))**2
2443 S2MIN2=MAX(S2MIN,QMIN2)
2444
2445 ZMIN=(QT2+AMJ0)/WW/(1.D0-X)
2446
2447
2448 S2=X*(1.D0-ZMIN)*WW-PT2
2449
2450
2451 IF(S2.LT.S2MIN2)GOTO 7
2452
2453 SJ1=PSJINT(QQ,QMIN(3-JJ),S2,1,IQP(3-jj)+1)
2454 * *QGSPSFAP(X,IQP(JJ),0)
2455 SJ2=PSJINT(QQ,QMIN(3-JJ),S2,2,IQP(3-jj)+1)
2456 * *QGSPSFAP(X,IQP(JJ),1)
2457
2458
2459
2460
2461
2462 GB7=(SJ1+SJ2)/DLOG(QT2/ALM)*QQ*QGSPSUDS(QQ,IQP(JJ))/GB0
2463
2464
2465 IF(X.LE..5D0)THEN
2466 GB7=GB7*X**(1.D0-DELH)
2467 ELSE
2468 GB7=GB7*(1.D0-X)*2.D0**DELH
2469 ENDIF
2470
2471 IF(PSRAN(B10).GT.GB7)GOTO 7
2472
2473 IF(PSRAN(B10).LT.SJ1/(SJ1+SJ2))THEN
2474 IF(IQC(JJ).EQ.0)THEN
2475 JT=1
2476 JQ=INT(1.5D0+PSRAN(B10))
2477 IQJ(1)=IPQ(JQ,JJ)
2478 IQJ(2)=0
2479 DO 31 I=1,4
2480 EQJ(I,1)=EPQ(I+4*(JQ-1),JJ)
2481 31 EQJ(I,2)=0.D0
2482 ELSE
2483 JT=2
2484 IF(IQC(JJ).GT.0)THEN
2485 JQ=1
2486 ELSE
2487 JQ=2
2488 ENDIF
2489 IQJ(1)=0
2490 DO 32 I=1,4
2491 32 EQJ(I,1)=0.D0
2492
2493 IPQ(JQ,JJ)=IPQ(1,JJ)
2494 DO 135 I=1,4
2495 135 EPQ(I+4*(JQ-1),JJ)=EPQ(I,JJ)
2496 ENDIF
2497 IQ1=IQC(JJ)
2498 IQC(JJ)=0
2499
2500 ELSE
2501 IF(IQP(JJ).NE.0)THEN
2502 IQ1=0
2503 JT=3
2504 IF(IQC(JJ).GT.0)THEN
2505 JQ=1
2506 ELSE
2507 JQ=2
2508 ENDIF
2509 IQJ(1)=IPQ(1,JJ)
2510 IQJ(2)=0
2511 DO 33 I=1,4
2512 EQJ(I,1)=EPQ(I,JJ)
2513 33 EQJ(I,2)=0.D0
2514
2515 ELSE
2516 IQ1=INT(3.D0*PSRAN(B10)+1.D0)*(2*INT(.5D0+PSRAN(B10))-1)
2517 IQC(JJ)=-IQ1
2518 JT=4
2519 IF(IQ1.GT.0)THEN
2520 JQ=1
2521 ELSE
2522 JQ=2
2523 ENDIF
2524 IQJ(1)=IPQ(JQ,JJ)
2525 DO 34 I=1,4
2526 34 EQJ(I,1)=EPQ(I+4*(JQ-1),JJ)
2527 ENDIF
2528 ENDIF
2529 IF(DEBUG.GE.3)WRITE (MONIOU,240)JT
2530
2531 CALL PSCAJET(QT2,IQ1,QV1,ZV1,QM1,IQV1,
2532 * LDAU1,LPAR1,JQ)
2533 Z=(QT2+QM1(1,1))/WW/(1.D0-X)
2534 SI=X*(1.D0-Z)*WW-PT2
2535
2536 IF(SI.GT.S2MIN2)THEN
2537 IQ=MIN(1,IABS(IQC(JJ)))+1
2538 GB=PSJINT(QQ,QMIN(3-JJ),SI,IQ,IQP(3-JJ)+1)/
2539 * PSJINT(QQ,QMIN(3-JJ),S2,IQ,IQP(3-JJ)+1)
2540 IF(PSRAN(B10).GT.GB)GOTO 301
2541 ELSE
2542 GOTO 301
2543 ENDIF
2544
2545 WP3=WP(JJ)*(1.D0-X)
2546 WM3=(QT2+QM1(1,1))/WP3
2547 EP3(1)=.5D0*(WP3+WM3)
2548 EP3(2)=.5D0*(WP3-WM3)*(3-2*JJ)
2549
2550 PT3=DSQRT(EP3(3)**2+EP3(4)**2)
2551
2552 CALL PSREC(EP3,QV1,ZV1,QM1,IQV1,LDAU1,LPAR1,IQJ,EQJ,JFL,JQ)
2553 IF(JFL.EQ.0)GOTO 301
2554
2555 IF(JT.EQ.1)THEN
2556 IPQ(JQ,JJ)=IQJ(2)
2557 DO 35 I=1,4
2558 35 EPQ(I+4*(JQ-1),JJ)=EQJ(I,2)
2559
2560 IF(IPC(JQ,JJ).EQ.0)THEN
2561 IPC(JQ,JJ)=IQJ(1)
2562 DO 36 I=1,4
2563 36 EPC(I+4*(JQ-1),JJ)=EQJ(I,1)
2564 ENDIF
2565
2566 ELSEIF(JT.EQ.2)THEN
2567 IPQ(3-JQ,JJ)=IQJ(1)
2568 DO 37 I=1,4
2569 37 EPQ(I+4*(2-JQ),JJ)=EQJ(I,1)
2570
2571 ELSEIF(JT.EQ.3)THEN
2572 IPQ(1,JJ)=IQJ(2)
2573 DO 38 I=1,4
2574 38 EPQ(I,JJ)=EQJ(I,2)
2575
2576 IF(IPC(JQ,JJ).EQ.0)THEN
2577 IPC(JQ,JJ)=IQJ(1)
2578 DO 39 I=1,4
2579 39 EPC(I+4*(JQ-1),JJ)=EQJ(I,1)
2580 ENDIF
2581
2582 ELSEIF(JT.EQ.4)THEN
2583 IF(IPC(JQ,JJ).EQ.0)THEN
2584 IPC(JQ,JJ)=IQJ(1)
2585 DO 40 I=1,4
2586 40 EPC(I+4*(JQ-1),JJ)=EQJ(I,1)
2587 ENDIF
2588 IF(JQ.EQ.1)THEN
2589 IPQ(1,JJ)=IPQ(2,JJ)
2590 DO 30 I=1,4
2591 30 EPQ(I,JJ)=EPQ(I+4,JJ)
2592 ENDIF
2593 ENDIF
2594
2595 IF(IABS(IQ1).EQ.3)THEN
2596 IQQQ=8+IQ1/3*4
2597 ELSE
2598 IQQQ=8+IQ1
2599 ENDIF
2600 IF(DEBUG.GE.2)WRITE (MONIOU,209)TYQ(IQQQ),QT2,EP3
2601 209 FORMAT(2X,'PSHOT: NEW JET FLAVOR:',A2,
2602 * ' PT SQUARED FOR THE JET:',E10.3/
2603 * 4X,'JET 4-MOMENTUM:',4E10.3)
2604 DO 8 I=1,4
2605 8 EPT(I)=EPT(I)-EP3(I)
2606
2607
2608 QMIN(JJ)=QQ
2609 QMINN=QMIN2
2610
2611
2612 GOTO 5
2613
2614
2615
2616
2617
2618 12 CONTINUE
2619 IF(DEBUG.GE.2)WRITE (MONIOU,211)SI
2620 211 FORMAT(2X,'PSHOT: HIGHEST VIRTUALITY SUBPROCESS IN THE LADDER'/
2621 * 4X,'MASS SQUARED FOR THE PROCESS:',E10.3)
2622
2623 XMIN=QMINN/(QMINN+SI)
2624 XMIN1=.5D0-DSQRT(.25D0-(QT0+AMJ0)/SI)
2625 XMIN=MAX(XMIN,XMIN1)
2626 TMIN=SI*XMIN
2627
2628 IF(IQC(1).NE.0.OR.IQC(2).NE.0)THEN
2629 GB0=TMIN**2/DLOG(TMIN*(1.D0-XMIN)/ALM)**2*
2630 * PSFBORN(SI,TMIN,IQC(1),IQC(2))
2631 ELSE
2632 GB0=.25D0*SI**2/DLOG(TMIN*(1.D0-XMIN)/ALM)**2*
2633 * PSFBORN(SI,.5D0*SI,IQC(1),IQC(2))
2634 ENDIF
2635
2636
2637
2638
2639 13 Q2=TMIN/(1.D0-PSRAN(B10)*(1.D0-2.D0*TMIN/SI))
2640 Z=Q2/SI
2641 QT2=Q2*(1.D0-Z)
2642 IF(PSRAN(B10).LT..5D0)THEN
2643 JM=2
2644 TQ=SI-Q2
2645 ELSE
2646 JM=1
2647 TQ=Q2
2648 ENDIF
2649
2650 GB=Q2**2/DLOG(QT2/ALM)**2/GB0*
2651 * PSFBORN(SI,TQ,IQC(1),IQC(2))
2652 IF(DEBUG.GE.3)WRITE (MONIOU,241)Q2,GB
2653 241 FORMAT(2X,'PSHOT: Q2=',E10.3,' GB=',E10.3)
2654
2655 IF(PSRAN(B10).GT.GB)GOTO 13
2656
2657 IF(IQC(1).EQ.0.AND.IQC(2).EQ.0)THEN
2658 JQ=INT(1.5D0+PSRAN(B10))
2659 IQJ(1)=IPQ(JQ,JM)
2660 DO 51 I=1,4
2661 51 EQJ(I,1)=EPQ(I+4*(JQ-1),JM)
2662
2663 IF(PSRAN(B10).LT..5D0)THEN
2664 JT=1
2665 IF(IPQ(3-JQ,JM)*IPQ(JQ,3-JM).NE.0)THEN
2666 IPJ=IPQ(3-JQ,JM)
2667 IPJ1=IPQ(JQ,3-JM)
2668 IF(IABS(IPJ).EQ.3)IPJ=IPJ*4/3
2669 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
2670 DO 52 I=1,4
2671 EPJ(I)=EPQ(I+4*(2-JQ),JM)
2672 52 EPJ1(I)=EPQ(I+4*(JQ-1),3-JM)
2673 CALL PSJDEF(IPJ,IPJ1,EPJ,EPJ1,JFL)
2674 IF(JFL.EQ.0)GOTO 301
2675 ELSEIF(IPQ(3-JQ,JM).NE.0)THEN
2676 IPC(JQ,3-JM)=IPQ(3-JQ,JM)
2677 DO 53 I=1,4
2678 53 EPC(I+4*(JQ-1),3-JM)=EPQ(I+4*(2-JQ),JM)
2679 ELSEIF(IPQ(JQ,3-JM).NE.0)THEN
2680 IPC(3-JQ,JM)=IPQ(JQ,3-JM)
2681 DO 54 I=1,4
2682 54 EPC(I+4*(2-JQ),JM)=EPQ(I+4*(JQ-1),3-JM)
2683 ENDIF
2684
2685 IQJ(2)=0
2686 DO 55 I=1,4
2687 55 EQJ(I,2)=0.D0
2688
2689 ELSE
2690 JT=2
2691 IQJ(2)=IPQ(3-JQ,3-JM)
2692 DO 56 I=1,4
2693 56 EQJ(I,2)=EPQ(I+4*(2-JQ),3-JM)
2694 ENDIF
2695
2696 ELSEIF(IQC(1)*IQC(2).EQ.0)THEN
2697 IF(IQC(1)+IQC(2).GT.0)THEN
2698 JQ=1
2699 ELSE
2700 JQ=2
2701 ENDIF
2702
2703 IF(PSRAN(B10).LT..5D0)THEN
2704 IF(IQC(JM).EQ.0)THEN
2705 JT=3
2706 IQJ(1)=IPQ(JQ,JM)
2707 IQJ(2)=0
2708 DO 57 I=1,4
2709 EQJ(I,1)=EPQ(I+4*(JQ-1),JM)
2710 57 EQJ(I,2)=0.D0
2711
2712 IF(IPQ(3-JQ,JM)*IPQ(1,3-JM).NE.0)THEN
2713 IPJ=IPQ(3-JQ,JM)
2714 IPJ1=IPQ(1,3-JM)
2715 IF(IABS(IPJ).EQ.3)IPJ=IPJ*4/3
2716 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
2717 DO 58 I=1,4
2718 EPJ(I)=EPQ(I+4*(2-JQ),JM)
2719 58 EPJ1(I)=EPQ(I,3-JM)
2720 CALL PSJDEF(IPJ,IPJ1,EPJ,EPJ1,JFL)
2721 IF(JFL.EQ.0)GOTO 301
2722 ELSEIF(IPQ(3-JQ,JM).NE.0)THEN
2723 IPC(JQ,3-JM)=IPQ(3-JQ,JM)
2724 DO 59 I=1,4
2725 59 EPC(I+4*(JQ-1),3-JM)=EPQ(I+4*(2-JQ),JM)
2726 ELSEIF(IPQ(1,3-JM).NE.0)THEN
2727 IPC(3-JQ,JM)=IPQ(1,3-JM)
2728 DO 60 I=1,4
2729 60 EPC(I+4*(2-JQ),JM)=EPQ(I,3-JM)
2730 ENDIF
2731
2732 ELSE
2733 JT=4
2734 IQJ(1)=0
2735 DO 61 I=1,4
2736 61 EQJ(I,1)=0.D0
2737
2738 IF(IPQ(1,JM)*IPQ(3-JQ,3-JM).NE.0)THEN
2739 IPJ=IPQ(1,JM)
2740 IPJ1=IPQ(3-JQ,3-JM)
2741 IF(IABS(IPJ).EQ.3)IPJ=IPJ*4/3
2742 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
2743 DO 62 I=1,4
2744 EPJ(I)=EPQ(I,JM)
2745 62 EPJ1(I)=EPQ(I+4*(2-JQ),3-JM)
2746 CALL PSJDEF(IPJ,IPJ1,EPJ,EPJ1,JFL)
2747 IF(JFL.EQ.0)GOTO 301
2748 ELSEIF(IPQ(3-JQ,3-JM).NE.0)THEN
2749 IPC(JQ,JM)=IPQ(3-JQ,3-JM)
2750 DO 63 I=1,4
2751 63 EPC(I+4*(JQ-1),JM)=EPQ(I+4*(2-JQ),3-JM)
2752 ELSEIF(IPQ(1,JM).NE.0)THEN
2753 IPC(3-JQ,3-JM)=IPQ(1,JM)
2754 DO 64 I=1,4
2755 64 EPC(I+4*(2-JQ),3-JM)=EPQ(I,JM)
2756 ENDIF
2757 ENDIF
2758
2759 ELSE
2760 IF(IQC(JM).EQ.0)THEN
2761 JT=5
2762 IQJ(2)=IPQ(3-JQ,JM)
2763 IQJ(1)=IPQ(1,3-JM)
2764 DO 65 I=1,4
2765 EQJ(I,2)=EPQ(I+4*(2-JQ),JM)
2766 65 EQJ(I,1)=EPQ(I,3-JM)
2767 ELSE
2768 JT=6
2769 IQJ(1)=IPQ(JQ,3-JM)
2770 DO 66 I=1,4
2771 66 EQJ(I,1)=EPQ(I+4*(JQ-1),3-JM)
2772 ENDIF
2773 ENDIF
2774
2775 ELSEIF(IQC(1)*IQC(2).GT.0)THEN
2776 JT=7
2777 IF(IQC(1).GT.0)THEN
2778 JQ=1
2779 ELSE
2780 JQ=2
2781 ENDIF
2782 IQJ(1)=IPQ(1,3-JM)
2783 DO 67 I=1,4
2784 67 EQJ(I,1)=EPQ(I,3-JM)
2785
2786 ELSE
2787 JT=8
2788 IF(IQC(JM).GT.0)THEN
2789 JQ=1
2790 ELSE
2791 JQ=2
2792 ENDIF
2793 IQJ(1)=0
2794 DO 68 I=1,4
2795 68 EQJ(I,1)=0.D0
2796
2797 IF(IPQ(1,JM)*IPQ(1,3-JM).NE.0)THEN
2798 IPJ=IPQ(1,JM)
2799 IPJ1=IPQ(1,3-JM)
2800 IF(IABS(IPJ).EQ.3)IPJ=IPJ*4/3
2801 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
2802 DO 69 I=1,4
2803 EPJ(I)=EPQ(I,JM)
2804 69 EPJ1(I)=EPQ(I,3-JM)
2805 CALL PSJDEF(IPJ,IPJ1,EPJ,EPJ1,JFL)
2806 IF(JFL.EQ.0)GOTO 301
2807 ELSEIF(IPQ(1,3-JM).NE.0)THEN
2808 IPC(JQ,JM)=IPQ(1,3-JM)
2809 DO 70 I=1,4
2810 70 EPC(I+4*(JQ-1),JM)=EPQ(I,3-JM)
2811 ELSEIF(IPQ(1,JM).NE.0)THEN
2812 IPC(3-JQ,3-JM)=IPQ(1,JM)
2813 DO 71 I=1,4
2814 71 EPC(I+4*(2-JQ),3-JM)=EPQ(I,JM)
2815 ENDIF
2816 ENDIF
2817 IF(JT.NE.8)THEN
2818 JQ2=JQ
2819 ELSE
2820 JQ2=3-JQ
2821 ENDIF
2822 IF(DEBUG.GE.3)WRITE (MONIOU,240)JT
2823 240 FORMAT(2X,'PSHOT: COLOUR CONNECTION JT=:',I1)
2824
2825 CALL PSCAJET(QT2,IQC(JM),QV1,ZV1,QM1,IQV1,
2826 * LDAU1,LPAR1,JQ)
2827 CALL PSCAJET(QT2,IQC(3-JM),QV2,ZV2,QM2,IQV2,
2828 * LDAU2,LPAR2,JQ2)
2829
2830 AMT1=QT2+QM1(1,1)
2831 AMT2=QT2+QM2(1,1)
2832
2833 IF(DSQRT(SI).GT.DSQRT(AMT1)+DSQRT(AMT2))THEN
2834 Z=XXTWDEC(SI,AMT1,AMT2)
2835 ELSE
2836 GOTO 301
2837 ENDIF
2838
2839 CALL QGSPSDEFTR(SI,EPT,EY)
2840
2841 WP3=Z*DSQRT(SI)
2842 WM3=(QT2+QM1(1,1))/WP3
2843 EP3(1)=.5D0*(WP3+WM3)
2844 EP3(2)=.5D0*(WP3-WM3)
2845 QT=DSQRT(QT2)
2846 CALL QGSPSCS(CCOS,SSIN)
2847
2848 EP3(3)=QT*CCOS
2849 EP3(4)=QT*SSIN
2850
2851 CALL QGSPSTRANS(EP3,EY)
2852 PT3=DSQRT(EP3(3)**2+EP3(4)**2)
2853
2854 CALL PSREC(EP3,QV1,ZV1,QM1,IQV1,LDAU1,LPAR1,IQJ,EQJ,JFL,JQ)
2855 IF(JFL.EQ.0)GOTO 301
2856
2857 if(iabs(IQC(JM)).eq.3)then
2858 iqqq=8+IQC(JM)/3*4
2859 else
2860 iqqq=8+IQC(JM)
2861 endif
2862 IF(DEBUG.GE.2)WRITE (MONIOU,209)TYQ(IQQQ),QT2
2863
2864 WP3=(1.D0-Z)*DSQRT(SI)
2865 WM3=(QT2+QM2(1,1))/WP3
2866 EP3(1)=.5D0*(WP3+WM3)
2867 EP3(2)=.5D0*(WP3-WM3)
2868 EP3(3)=-QT*CCOS
2869 EP3(4)=-QT*SSIN
2870 CALL QGSPSTRANS(EP3,EY)
2871 PT3=DSQRT(EP3(3)**2+EP3(4)**2)
2872
2873 IF(JT.EQ.1)THEN
2874 IF(IPC(JQ,JM).EQ.0)THEN
2875 IPC(JQ,JM)=IQJ(1)
2876 DO 72 I=1,4
2877 72 EPC(I+4*(JQ-1),JM)=EQJ(I,1)
2878 ENDIF
2879
2880 IQJ(1)=IQJ(2)
2881 IQJ(2)=IPQ(3-JQ,3-JM)
2882 DO 73 I=1,4
2883 EQJ(I,1)=EQJ(I,2)
2884 73 EQJ(I,2)=EPQ(I+4*(2-JQ),3-JM)
2885
2886 ELSEIF(JT.EQ.2)THEN
2887 IF(IPC(JQ,JM).EQ.0)THEN
2888 IPC(JQ,JM)=IQJ(1)
2889 DO 74 I=1,4
2890 74 EPC(I+4*(JQ-1),JM)=EQJ(I,1)
2891 ENDIF
2892 IF(IPC(3-JQ,3-JM).EQ.0)THEN
2893 IPC(3-JQ,3-JM)=IQJ(2)
2894 DO 75 I=1,4
2895 75 EPC(I+4*(2-JQ),3-JM)=EQJ(I,2)
2896 ENDIF
2897
2898 IQJ(2)=IPQ(3-JQ,JM)
2899 IQJ(1)=IPQ(JQ,3-JM)
2900 DO 76 I=1,4
2901 EQJ(I,2)=EPQ(I+4*(2-JQ),JM)
2902 76 EQJ(I,1)=EPQ(I+4*(JQ-1),3-JM)
2903
2904 ELSEIF(JT.EQ.3)THEN
2905 IF(IPC(JQ,JM).EQ.0)THEN
2906 IPC(JQ,JM)=IQJ(1)
2907 DO 77 I=1,4
2908 77 EPC(I+4*(JQ-1),JM)=EQJ(I,1)
2909 ENDIF
2910 IQJ(1)=IQJ(2)
2911 DO 78 I=1,4
2912 78 EQJ(I,1)= EQJ(I,2)
2913
2914 ELSEIF(JT.EQ.4)THEN
2915 IQJ(2)=IQJ(1)
2916 IQJ(1)=IPQ(JQ,3-JM)
2917 DO 79 I=1,4
2918 EQJ(I,2)=EQJ(I,1)
2919 79 EQJ(I,1)=EPQ(I+4*(JQ-1),3-JM)
2920
2921 ELSEIF(JT.EQ.5)THEN
2922 IF(IPC(3-JQ,JM).EQ.0)THEN
2923 IPC(3-JQ,JM)=IQJ(2)
2924 DO 80 I=1,4
2925 80 EPC(I+4*(2-JQ),JM)=EQJ(I,2)
2926 ENDIF
2927 IF(IPC(JQ,3-JM).EQ.0)THEN
2928 IPC(JQ,3-JM)=IQJ(1)
2929 DO 81 I=1,4
2930 81 EPC(I+4*(JQ-1),3-JM)=EQJ(I,1)
2931 ENDIF
2932
2933 IQJ(1)=IPQ(JQ,JM)
2934 DO 82 I=1,4
2935 82 EQJ(I,1)=EPQ(I+4*(JQ-1),JM)
2936
2937 ELSEIF(JT.EQ.6)THEN
2938 IF(IPC(JQ,3-JM).EQ.0)THEN
2939 IPC(JQ,3-JM)=IQJ(1)
2940 DO 83 I=1,4
2941 83 EPC(I+4*(JQ-1),3-JM)=EQJ(I,1)
2942 ENDIF
2943
2944 IQJ(2)=IPQ(3-JQ,3-JM)
2945 IQJ(1)=IPQ(1,JM)
2946 DO 84 I=1,4
2947 EQJ(I,2)=EPQ(I+4*(2-JQ),3-JM)
2948 84 EQJ(I,1)=EPQ(I,JM)
2949
2950 ELSEIF(JT.EQ.7)THEN
2951 IF(IPC(JQ,3-JM).EQ.0)THEN
2952 IPC(JQ,3-JM)=IQJ(1)
2953 DO 85 I=1,4
2954 85 EPC(I+4*(JQ-1),3-JM)=EQJ(I,1)
2955 ENDIF
2956 IQJ(1)=IPQ(1,JM)
2957 DO 86 I=1,4
2958 86 EQJ(I,1)= EPQ(I,JM)
2959 ENDIF
2960
2961 CALL PSREC(EP3,QV2,ZV2,QM2,IQV2,LDAU2,LPAR2,IQJ,EQJ,JFL,JQ2)
2962 IF(JFL.EQ.0)GOTO 301
2963
2964 if(iabs(IQC(3-JM)).eq.3)then
2965 iqqq=8+IQC(3-JM)/3*4
2966 else
2967 iqqq=8+IQC(3-JM)
2968 endif
2969 IF(DEBUG.GE.2)WRITE (MONIOU,209)TYQ(IQQQ),QT2
2970 IF(DEBUG.GE.2)WRITE (MONIOU,212)NJTOT
2971 212 FORMAT(2X,'PSHOT: TOTAL NUMBER OF JETS:',I2)
2972
2973 IF(JT.EQ.1)THEN
2974 IF(IPC(3-JQ,3-JM).EQ.0)THEN
2975 IPC(3-JQ,3-JM)=IQJ(2)
2976 DO 87 I=1,4
2977 87 EPC(I+4*(2-JQ),3-JM)=EQJ(I,2)
2978 ENDIF
2979
2980 ELSEIF(JT.EQ.2)THEN
2981 IF(IPC(3-JQ,JM).EQ.0)THEN
2982 IPC(3-JQ,JM)=IQJ(2)
2983 DO 88 I=1,4
2984 88 EPC(I+4*(2-JQ),JM)=EQJ(I,2)
2985 ENDIF
2986 IF(IPC(JQ,3-JM).EQ.0)THEN
2987 IPC(JQ,3-JM)=IQJ(1)
2988 DO 89 I=1,4
2989 89 EPC(I+4*(JQ-1),3-JM)=EQJ(I,1)
2990 ENDIF
2991
2992 ELSEIF(JT.EQ.4)THEN
2993 IF(IPC(JQ,3-JM).EQ.0)THEN
2994 IPC(JQ,3-JM)=IQJ(1)
2995 DO 90 I=1,4
2996 90 EPC(I+4*(JQ-1),3-JM)=EQJ(I,1)
2997 ENDIF
2998
2999 ELSEIF(JT.EQ.5)THEN
3000 IF(IPC(JQ,JM).EQ.0)THEN
3001 IPC(JQ,JM)=IQJ(1)
3002 DO 91 I=1,4
3003 91 EPC(I+4*(JQ-1),JM)=EQJ(I,1)
3004 ENDIF
3005
3006 ELSEIF(JT.EQ.6)THEN
3007 IF(IPC(3-JQ,3-JM).EQ.0)THEN
3008 IPC(3-JQ,3-JM)=IQJ(2)
3009 DO 92 I=1,4
3010 92 EPC(I+4*(2-JQ),3-JM)=EQJ(I,2)
3011 ENDIF
3012 IF(IPC(JQ,JM).EQ.0)THEN
3013 IPC(JQ,JM)=IQJ(1)
3014 DO 93 I=1,4
3015 93 EPC(I+4*(JQ-1),JM)=EQJ(I,1)
3016 ENDIF
3017
3018 ELSEIF(JT.EQ.7)THEN
3019 IF(IPC(JQ,JM).EQ.0)THEN
3020 IPC(JQ,JM)=IQJ(1)
3021 DO 94 I=1,4
3022 94 EPC(I+4*(JQ-1),JM)=EQJ(I,1)
3023 ENDIF
3024 ENDIF
3025
3026
3027 IF(DEBUG.GE.3)WRITE (MONIOU,217)
3028 217 FORMAT(2X,'PSHOT - END')
3029 ebal(1)=.5*(wp0+wm0)
3030 ebal(2)=.5*(wp0-wm0)
3031 ebal(3)=0.d0
3032 ebal(4)=0.d0
3033 do 500 i=1,njtot
3034 do 500 m=1,2
3035 do 500 l=1,4
3036 500 ebal(l)=ebal(l)-epjet(l,m,i)
3037
3038 RETURN
3039 END
3040
3041
3042 SUBROUTINE PSJDEF(IPJ,IPJ1,EPJ,EPJ1,JFL)
3043
3044
3045
3046
3047 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3048 INTEGER DEBUG
3049 DIMENSION EPJ(4),EPJ1(4),EPT(4)
3050 COMMON /Q_AREA10/ STMASS,AM(7)
3051 COMMON /Q_AREA43/ MONIOU
3052 COMMON /Q_DEBUG/ DEBUG
3053 COMMON /Q_AREA46/ EPJET(4,2,15000),IPJET(2,15000)
3054 COMMON /Q_AREA47/ NJTOT
3055 SAVE
3056
3057
3058
3059
3060
3061
3062
3063
3064
3065 IF(DEBUG.GE.2)WRITE (MONIOU,201)IPJ,IPJ1,EPJ,EPJ1
3066 201 FORMAT(2X,'PSJDEF: PARTON FLAVORS',
3067 * ': IPJ=',I2,2X,'IPJ1=',I2/
3068 * 4X,'PARTON 4-MOMENTA:',2X,4(E10.3,1X))
3069 DO 1 I=1,4
3070 1 EPT(I)=EPJ(I)+EPJ1(I)
3071
3072
3073 WW=QGSPSNORM(EPt)
3074
3075 IF(IABS(IPJ).LE.2)THEN
3076 AM1=AM(1)
3077 ELSEIF(IABS(IPJ).EQ.4)THEN
3078 AM1=AM(3)
3079 ELSE
3080 AM1=AM(2)
3081 ENDIF
3082 IF(IABS(IPJ1).LE.2)THEN
3083 AM2=AM(1)
3084 ELSEIF(IABS(IPJ1).EQ.4)THEN
3085 AM2=AM(3)
3086 ELSE
3087 AM2=AM(2)
3088 ENDIF
3089 AMJ=(AM1+AM2)**2
3090
3091 IF(AMJ.GT.WW)THEN
3092 JFL=0
3093 RETURN
3094 ELSE
3095 JFL=1
3096 ENDIF
3097
3098 NJTOT=NJTOT+1
3099 IF( NJTOT . GT. 15000 ) THEN
3100 WRITE(MONIOU,*)'PSJDEF: TOO MANY JETS'
3101 WRITE(MONIOU,*)'PSJDEF: NJTOT = ',NJTOT
3102 STOP
3103 ENDIF
3104 IPJET(1,NJTOT)=IPJ
3105 IPJET(2,NJTOT)=IPJ1
3106 DO 2 I=1,4
3107 EPJET(I,1,NJTOT)=EPJ(I)
3108 2 EPJET(I,2,NJTOT)=EPJ1(I)
3109
3110 IF(DEBUG.GE.3)WRITE (MONIOU,202)
3111 202 FORMAT(2X,'PSJDEF - END')
3112 RETURN
3113 END
3114
3115
3116 FUNCTION QGSPSJET(Q1,Q2,S,S2MIN,J,L)
3117
3118
3119
3120
3121
3122
3123
3124
3125
3126 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3127 INTEGER DEBUG
3128 COMMON /Q_AREA6/ PI,BM,AM
3129 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALF,RR,SH,DELH
3130 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
3131 COMMON/AR13/X1(7),A1(7)
3132 COMMON /Q_AREA43/ MONIOU
3133 COMMON /Q_DEBUG/ DEBUG
3134 SAVE
3135
3136 IF(DEBUG.GE.2)WRITE (MONIOU,201)S,Q1,Q2,S2MIN,J,L
3137 201 FORMAT(2X,'QGSPSJET - UNORDERED LADDER CROSS SECTION:'/
3138 * 4X,'S=',E10.3,2X,'Q1=',E10.3,2X,'Q2=',E10.3,2X,'S2MIN=',
3139 * E10.3,2X,'J=',I1,2X,'L=',I1)
3140 QGSPSJET=0.D0
3141
3142 P=DSQRT(1.D0-3.D0*QT0/S)
3143 COSF=(1.D0-18.D0*QT0/S)/P**3
3144 FI=ATAN(1.D0/COSF**2-1.D0)
3145 IF(COSF.LT.0.D0)FI=PI-FI
3146 FI=FI/3.D0
3147 ZMAX=(2.D0-P*(DSQRT(3.D0)*SIN(FI)-COS(FI)))/3.D0
3148 ZMIN=(1.D0-P*COS(FI))/1.5D0
3149
3150 IF(QT0/(1.D0-ZMIN)**2.LT.S2MIN)
3151 * ZMIN=.5D0*(1.D0+S2MIN/S-DSQRT((1.D0-S2MIN/S)**2-4.D0*QT0/S))
3152
3153
3154 IF(1.D0-ZMIN.LT.DSQRT(QT0/Q1))THEN
3155 QMIN=QT0/(1.D0-ZMIN)**2
3156 ELSE
3157 QMIN=Q1
3158 ENDIF
3159
3160 QMAX=QT0/(1.D0-ZMAX)**2
3161 SUD0=QGSPSUDS(QMIN,J)
3162
3163
3164 IF(DEBUG.GE.3)WRITE (MONIOU,203)QMIN,QMAX
3165 203 FORMAT(2X,'QGSPSJET:',2X,'QMIN=',E10.3,2X,'QMAX=',E10.3)
3166 IF(QMAX.GT.QMIN)THEN
3167
3168
3169
3170 DO 3 I=1,7
3171 DO 3 M=1,2
3172 QI=2.D0*QMIN/(1.D0+QMIN/QMAX+(2*M-3)*X1(I)*(1.D0-QMIN/QMAX))
3173
3174 ZMAX=(1.D0-DSQRT(QT0/QI))**DELH
3175 ZMIN=((QI+MAX(QI,S2MIN))/(QI+S))**DELH
3176
3177 FSJ=0.D0
3178
3179 IF(DEBUG.GE.3)WRITE (MONIOU,204)QI,ZMIN,ZMAX
3180 204 FORMAT(2X,'QGSPSJET:',2X,'QI=',E10.3,2X,'ZMIN=',E10.3,2X,
3181 * 'ZMAX=',E10.3)
3182 IF(ZMAX.GT.ZMIN)THEN
3183 DO 2 I1=1,7
3184 DO 2 M1=1,2
3185 Z=(.5D0*(ZMAX+ZMIN+(2*M1-3)*X1(I1)*(ZMAX-ZMIN)))**
3186 * (1.D0/DELH)
3187 QT=QI*(1.D0-Z)**2
3188 S2=Z*S-QI*(1.D0-Z)
3189
3190 SJ=0.D0
3191 DO 1 K=1,2
3192 1 SJ=SJ+PSJINT(QI,Q2,S2,K,L)*QGSPSFAP(Z,J,K-1)*Z
3193 2 FSJ=FSJ+A1(I1)*SJ/DLOG(QT/ALM)/Z**DELH
3194 FSJ=FSJ*(ZMAX-ZMIN)
3195 ENDIF
3196
3197 3 QGSPSJET=QGSPSJET+A1(I)*FSJ*QI*QGSPSUDS(QI,J)
3198 QGSPSJET=QGSPSJET*(1.D0/QMIN-1.D0/QMAX)/SUD0/DELH/18.D0
3199 ENDIF
3200 IF(DEBUG.GE.3)WRITE (MONIOU,202)QGSPSJET
3201 202 FORMAT(2X,'QGSPSJET=',E10.3)
3202 RETURN
3203 END
3204
3205
3206 FUNCTION QGSPSJET1(Q1,Q2,S,S2MIN,J,L)
3207
3208
3209
3210
3211
3212
3213
3214
3215
3216 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3217 INTEGER DEBUG
3218 COMMON /Q_AREA6/ PI,BM,AM
3219 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALF,RR,SH,DELH
3220 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
3221 COMMON/AR13/X1(7),A1(7)
3222 COMMON /Q_AREA43/ MONIOU
3223 COMMON /Q_DEBUG/ DEBUG
3224 SAVE
3225
3226 IF(DEBUG.GE.2)WRITE (MONIOU,201)S,Q1,Q2,S2MIN,J,L
3227 201 FORMAT(2X,'QGSPSJET1 - STRICTLY ORDERED LADDER CROSS SECTION:'/
3228 * 4X,'S=',E10.3,2X,'Q1=',E10.3,2X,'Q2=',E10.3,2X,'S2MIN=',
3229 * E10.3,2X,'J=',I1,2X,'L=',I1)
3230 QGSPSJET1=0.D0
3231
3232 P=DSQRT(1.D0-3.D0*QT0/S)
3233 COSF=(1.D0-18.D0*QT0/S)/P**3
3234 FI=ATAN(1.D0/COSF**2-1.D0)
3235 IF(COSF.LT.0.D0)FI=PI-FI
3236 FI=FI/3.D0
3237 ZMAX=(2.D0-P*(DSQRT(3.D0)*SIN(FI)-COS(FI)))/3.D0
3238 ZMIN=(1.D0-P*COS(FI))/1.5D0
3239
3240 IF(QT0/(1.D0-ZMIN)**2.LT.S2MIN)
3241 * ZMIN=.5D0*(1.D0+S2MIN/S-DSQRT((1.D0-S2MIN/S)**2-4.D0*QT0/S))
3242
3243
3244 IF(1.D0-ZMIN.LT.DSQRT(QT0/Q1))THEN
3245 QMIN=QT0/(1.D0-ZMIN)**2
3246 ELSE
3247 QMIN=Q1
3248 ENDIF
3249
3250 QMAX=QT0/(1.D0-ZMAX)**2
3251 SUD0=QGSPSUDS(QMIN,J)
3252
3253
3254 IF(DEBUG.GE.3)WRITE (MONIOU,203)QMIN,QMAX
3255 203 FORMAT(2X,'QGSPSJET1:',2X,'QMIN=',E10.3,2X,'QMAX=',E10.3)
3256 IF(QMAX.GT.QMIN)THEN
3257
3258
3259
3260 DO 3 I=1,7
3261 DO 3 M=1,2
3262 QI=2.D0*QMIN/(1.D0+QMIN/QMAX+(2*M-3)*X1(I)*(1.D0-QMIN/QMAX))
3263
3264 ZMAX=(1.D0-DSQRT(QT0/QI))**DELH
3265 ZMIN=((QI+MAX(QI,S2MIN))/(QI+S))**DELH
3266
3267 FSJ=0.D0
3268
3269 IF(DEBUG.GE.3)WRITE (MONIOU,204)QI,ZMIN,ZMAX
3270 204 FORMAT(2X,'QGSPSJET1:',2X,'QI=',E10.3,2X,'ZMIN=',E10.3,2X,
3271 * 'ZMAX=',E10.3)
3272 IF(ZMAX.GT.ZMIN)THEN
3273 DO 2 I1=1,7
3274 DO 2 M1=1,2
3275 Z=(.5D0*(ZMAX+ZMIN+(2*M1-3)*X1(I1)*(ZMAX-ZMIN)))**
3276 * (1.D0/DELH)
3277 QT=QI*(1.D0-Z)**2
3278 S2=Z*S-QI*(1.D0-Z)
3279
3280 SJ=0.D0
3281 DO 1 K=1,2
3282 1 SJ=SJ+PSJINT1(QI,Q2,S2,K,L)*QGSPSFAP(Z,J,K-1)*Z
3283 2 FSJ=FSJ+A1(I1)*SJ/DLOG(QT/ALM)/Z**DELH
3284 FSJ=FSJ*(ZMAX-ZMIN)
3285 ENDIF
3286
3287 3 QGSPSJET1=QGSPSJET1+A1(I)*FSJ*QI*QGSPSUDS(QI,J)
3288 QGSPSJET1=QGSPSJET1*(1.D0/QMIN-1.D0/QMAX)/SUD0/DELH/18.D0
3289 ENDIF
3290 IF(DEBUG.GE.3)WRITE (MONIOU,202)QGSPSJET1
3291 202 FORMAT(2X,'QGSPSJET1=',E10.3)
3292 RETURN
3293 END
3294
3295
3296 FUNCTION PSJINT(Q1,Q2,S,M,L)
3297
3298
3299
3300
3301
3302
3303
3304
3305 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3306 INTEGER DEBUG
3307 DIMENSION WI(3),WJ(3),WK(3)
3308 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
3309 COMMON /Q_AREA29/ CSJ(17,17,68)
3310 COMMON /Q_AREA43/ MONIOU
3311 COMMON /Q_DEBUG/ DEBUG
3312 SAVE
3313
3314 IF(DEBUG.GE.2)WRITE (MONIOU,201)S,Q1,Q2,M,L
3315 201 FORMAT(2X,'PSJINT - UNORDERED LADDER CROSS SECTION INTERPOL.:'/
3316 * 4X,'S=',E10.3,2X,'Q1=',E10.3,2X,'Q2=',E10.3,2X,
3317 * 2X,'M=',I1,2X,'L=',I1)
3318 PSJINT=0.D0
3319 QQ=MAX(Q1,Q2)
3320 IF(S.LE.MAX(4.D0*QT0,QQ))THEN
3321 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSJINT
3322 202 FORMAT(2X,'PSJINT=',E10.3)
3323 RETURN
3324 ENDIF
3325
3326 ML=17*(M-1)+34*(L-1)
3327 QLI=DLOG(Q1/QT0)/1.38629D0
3328 QLJ=DLOG(Q2/QT0)/1.38629D0
3329 SL=DLOG(S/QT0)/1.38629D0
3330 SQL=SL-MAX(QLI,QLJ)
3331 I=INT(QLI)
3332 J=INT(QLJ)
3333 K=INT(SL)
3334 IF(I.GT.13)I=13
3335 IF(J.GT.13)J=13
3336
3337 IF(SQL.GT.10.D0)THEN
3338 IF(K.GT.14)K=14
3339 IF(I.GT.K-3)I=K-3
3340 IF(J.GT.K-3)J=K-3
3341 WI(2)=QLI-I
3342 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
3343 WI(1)=1.D0-WI(2)+WI(3)
3344 WI(2)=WI(2)-2.D0*WI(3)
3345 WJ(2)=QLJ-J
3346 WJ(3)=WJ(2)*(WJ(2)-1.D0)*.5D0
3347 WJ(1)=1.D0-WJ(2)+WJ(3)
3348 WJ(2)=WJ(2)-2.D0*WJ(3)
3349 WK(2)=SL-K
3350 WK(3)=WK(2)*(WK(2)-1.D0)*.5D0
3351 WK(1)=1.D0-WK(2)+WK(3)
3352 WK(2)=WK(2)-2.D0*WK(3)
3353
3354 DO 1 I1=1,3
3355 DO 1 J1=1,3
3356 DO 1 K1=1,3
3357 1 PSJINT=PSJINT+CSJ(I+I1,J+J1,K+K1+ML)*WI(I1)*WJ(J1)*WK(K1)
3358 PSJINT=EXP(PSJINT)
3359 ELSEIF(SQL.LT.1.D0.AND.I+J.NE.0)THEN
3360 SQ=(S/MAX(Q1,Q2)-1.D0)/3.D0
3361 WI(2)=QLI-I
3362 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
3363 WI(1)=1.D0-WI(2)+WI(3)
3364 WI(2)=WI(2)-2.D0*WI(3)
3365 WJ(2)=QLJ-J
3366 WJ(3)=WJ(2)*(WJ(2)-1.D0)*.5D0
3367 WJ(1)=1.D0-WJ(2)+WJ(3)
3368 WJ(2)=WJ(2)-2.D0*WJ(3)
3369
3370 DO 2 I1=1,3
3371 I2=I+I1
3372 DO 2 J1=1,3
3373 J2=J+J1
3374 K2=MAX(I2,J2)+1+ML
3375 2 PSJINT=PSJINT+CSJ(I2,J2,K2)*WI(I1)*WJ(J1)
3376 PSJINT=EXP(PSJINT)*SQ
3377 ELSEIF(K.EQ.1)THEN
3378 SQ=(S/QT0/4.D0-1.D0)/3.D0
3379 WI(2)=QLI
3380 WI(1)=1.D0-QLI
3381 WJ(2)=QLJ
3382 WJ(1)=1.D0-QLJ
3383
3384 DO 3 I1=1,2
3385 DO 3 J1=1,2
3386 3 PSJINT=PSJINT+CSJ(I1,J1,3+ML)*WI(I1)*WJ(J1)
3387 PSJINT=EXP(PSJINT)*SQ
3388 ELSEIF(K.LT.15)THEN
3389 KL=INT(SQL)
3390 IF(I+KL.GT.12)I=12-KL
3391 IF(J+KL.GT.12)J=12-KL
3392 IF(I+J+KL.EQ.1)KL=2
3393 WI(2)=QLI-I
3394 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
3395 WI(1)=1.D0-WI(2)+WI(3)
3396 WI(2)=WI(2)-2.D0*WI(3)
3397 WJ(2)=QLJ-J
3398 WJ(3)=WJ(2)*(WJ(2)-1.D0)*.5D0
3399 WJ(1)=1.D0-WJ(2)+WJ(3)
3400 WJ(2)=WJ(2)-2.D0*WJ(3)
3401 WK(2)=SQL-KL
3402 WK(3)=WK(2)*(WK(2)-1.D0)*.5D0
3403 WK(1)=1.D0-WK(2)+WK(3)
3404 WK(2)=WK(2)-2.D0*WK(3)
3405
3406 DO 4 I1=1,3
3407 I2=I+I1
3408 DO 4 J1=1,3
3409 J2=J+J1
3410 DO 4 K1=1,3
3411 K2=MAX(I2,J2)+KL+K1-1+ML
3412 4 PSJINT=PSJINT+CSJ(I2,J2,K2)*WI(I1)*WJ(J1)*WK(K1)
3413 PSJINT=EXP(PSJINT)
3414 ELSE
3415 K=15
3416 IF(I.GT.K-3)I=K-3
3417 IF(J.GT.K-3)J=K-3
3418 WI(2)=QLI-I
3419 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
3420 WI(1)=1.D0-WI(2)+WI(3)
3421 WI(2)=WI(2)-2.D0*WI(3)
3422 WJ(2)=QLJ-J
3423 WJ(3)=WJ(2)*(WJ(2)-1.D0)*.5D0
3424 WJ(1)=1.D0-WJ(2)+WJ(3)
3425 WJ(2)=WJ(2)-2.D0*WJ(3)
3426 WK(2)=SL-K
3427 WK(1)=1.D0-WK(2)
3428
3429 DO 5 I1=1,3
3430 DO 5 J1=1,3
3431 DO 5 K1=1,2
3432 5 PSJINT=PSJINT+CSJ(I+I1,J+J1,K+K1+ML)*WI(I1)*WJ(J1)*WK(K1)
3433 PSJINT=EXP(PSJINT)
3434 ENDIF
3435 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSJINT
3436 RETURN
3437 END
3438
3439
3440 SUBROUTINE PSJINT0(S,SJ,SJB,M,L)
3441
3442
3443
3444
3445
3446
3447
3448
3449 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3450 INTEGER DEBUG
3451 DIMENSION WK(3)
3452 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
3453 COMMON /Q_AREA32/ CSJ(17,2,2),CSB(17,2,2)
3454 COMMON /Q_AREA43/ MONIOU
3455 COMMON /Q_DEBUG/ DEBUG
3456 SAVE
3457
3458 IF(DEBUG.GE.2)WRITE (MONIOU,201)S,M,L
3459 201 FORMAT(2X,'PSJINT0 - HARD CROSS SECTION INTERPOLATION:'/
3460 * 4X,'S=',E10.3,2X,'M=',I1,2X,'L=',I1)
3461 SJ=0.D0
3462 SJB=0.D0
3463 IF(S.LE.4.D0*QT0)THEN
3464 IF(DEBUG.GE.3)WRITE (MONIOU,202)SJ,SJB
3465 202 FORMAT(2X,'PSJINT0: SJ=',E10.3,2X,'SJB=',E10.3)
3466 RETURN
3467 ENDIF
3468
3469 SL=DLOG(S/QT0)/1.38629d0
3470 K=INT(SL)
3471 IF(K.EQ.1)THEN
3472 SQ=(S/QT0/4.D0-1.D0)/3.D0
3473 SJB=EXP(CSB(3,M+1,L+1))*SQ
3474 SJ=EXP(CSJ(3,M+1,L+1))*SQ
3475 ELSE
3476 IF(K.GT.14)K=14
3477 WK(2)=SL-K
3478 WK(3)=WK(2)*(WK(2)-1.D0)*.5D0
3479 WK(1)=1.D0-WK(2)+WK(3)
3480 WK(2)=WK(2)-2.D0*WK(3)
3481
3482 DO 1 K1=1,3
3483 SJ=SJ+CSJ(K+K1,M+1,L+1)*WK(K1)
3484 1 SJB=SJB+CSB(K+K1,M+1,L+1)*WK(K1)
3485 SJB=EXP(SJB)
3486 SJ=EXP(SJ)
3487 ENDIF
3488 IF(DEBUG.GE.3)WRITE (MONIOU,202)SJ,SJB
3489 RETURN
3490 END
3491
3492
3493 FUNCTION PSJINT1(Q1,Q2,S,M,L)
3494
3495
3496
3497
3498
3499
3500
3501
3502 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3503 INTEGER DEBUG
3504 DIMENSION WI(3),WJ(3),WK(3)
3505 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
3506 COMMON /Q_AREA30/ CSJ(17,17,68)
3507 COMMON /Q_AREA43/ MONIOU
3508 COMMON /Q_DEBUG/ DEBUG
3509 SAVE
3510
3511 IF(DEBUG.GE.2)WRITE (MONIOU,201)S,Q1,Q2,M,L
3512 201 FORMAT(2X,'PSJINT1 - STRICTLY ORDERED LADDER CROSS SECTION',
3513 * ' INTERPOLATION:'/
3514 * 4X,'S=',E10.3,2X,'Q1=',E10.3,2X,'Q2=',E10.3,2X,
3515 * 4X,'M=',I1,2X,'L=',I1)
3516 PSJINT1=0.D0
3517 QQ=MAX(Q1,Q2)
3518 IF(S.LE.MAX(4.D0*QT0,QQ))THEN
3519 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSJINT1
3520 202 FORMAT(2X,'PSJINT1=',E10.3)
3521 RETURN
3522 ENDIF
3523
3524 ML=17*(M-1)+34*(L-1)
3525 QLI=DLOG(Q1/QT0)/1.38629d0
3526 QLJ=DLOG(Q2/QT0)/1.38629d0
3527 SL=DLOG(S/QT0)/1.38629d0
3528 SQL=SL-MAX(QLI,QLJ)
3529 I=INT(QLI)
3530 J=INT(QLJ)
3531 K=INT(SL)
3532 IF(I.GT.13)I=13
3533 IF(J.GT.13)J=13
3534
3535 IF(SQL.GT.10.D0)THEN
3536 IF(K.GT.14)K=14
3537 IF(I.GT.K-3)I=K-3
3538 IF(J.GT.K-3)J=K-3
3539 WI(2)=QLI-I
3540 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
3541 WI(1)=1.D0-WI(2)+WI(3)
3542 WI(2)=WI(2)-2.D0*WI(3)
3543 WJ(2)=QLJ-J
3544 WJ(3)=WJ(2)*(WJ(2)-1.D0)*.5D0
3545 WJ(1)=1.D0-WJ(2)+WJ(3)
3546 WJ(2)=WJ(2)-2.D0*WJ(3)
3547 WK(2)=SL-K
3548 WK(3)=WK(2)*(WK(2)-1.D0)*.5D0
3549 WK(1)=1.D0-WK(2)+WK(3)
3550 WK(2)=WK(2)-2.D0*WK(3)
3551
3552 DO 1 I1=1,3
3553 DO 1 J1=1,3
3554 DO 1 K1=1,3
3555 1 PSJINT1=PSJINT1+CSJ(I+I1,J+J1,K+K1+ML)*WI(I1)*WJ(J1)*WK(K1)
3556 PSJINT1=EXP(PSJINT1)
3557 ELSEIF(SQL.LT.1.D0.AND.I+J.NE.0)THEN
3558 SQ=(S/MAX(Q1,Q2)-1.D0)/3.D0
3559 WI(2)=QLI-I
3560 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
3561 WI(1)=1.D0-WI(2)+WI(3)
3562 WI(2)=WI(2)-2.D0*WI(3)
3563 WJ(2)=QLJ-J
3564 WJ(3)=WJ(2)*(WJ(2)-1.D0)*.5D0
3565 WJ(1)=1.D0-WJ(2)+WJ(3)
3566 WJ(2)=WJ(2)-2.D0*WJ(3)
3567
3568 DO 2 I1=1,3
3569 I2=I+I1
3570 DO 2 J1=1,3
3571 J2=J+J1
3572 K2=MAX(I2,J2)+1+ML
3573 2 PSJINT1=PSJINT1+CSJ(I2,J2,K2)*WI(I1)*WJ(J1)
3574 PSJINT1=EXP(PSJINT1)*SQ
3575 ELSEIF(K.EQ.1)THEN
3576 SQ=(S/QT0/4.D0-1.D0)/3.D0
3577 WI(2)=QLI
3578 WI(1)=1.D0-QLI
3579 WJ(2)=QLJ
3580 WJ(1)=1.D0-QLJ
3581
3582 DO 3 I1=1,2
3583 DO 3 J1=1,2
3584 3 PSJINT1=PSJINT1+CSJ(I1,J1,3+ML)*WI(I1)*WJ(J1)
3585 PSJINT1=EXP(PSJINT1)*SQ
3586 ELSEIF(K.LT.15)THEN
3587 KL=INT(SQL)
3588 IF(I+KL.GT.12)I=12-KL
3589 IF(J+KL.GT.12)J=12-KL
3590 IF(I+J+KL.EQ.1)KL=2
3591
3592 WI(2)=QLI-I
3593 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
3594 WI(1)=1.D0-WI(2)+WI(3)
3595 WI(2)=WI(2)-2.D0*WI(3)
3596 WJ(2)=QLJ-J
3597 WJ(3)=WJ(2)*(WJ(2)-1.D0)*.5D0
3598 WJ(1)=1.D0-WJ(2)+WJ(3)
3599 WJ(2)=WJ(2)-2.D0*WJ(3)
3600 WK(2)=SQL-KL
3601 WK(3)=WK(2)*(WK(2)-1.D0)*.5D0
3602 WK(1)=1.D0-WK(2)+WK(3)
3603 WK(2)=WK(2)-2.D0*WK(3)
3604
3605 DO 4 I1=1,3
3606 I2=I+I1
3607 DO 4 J1=1,3
3608 J2=J+J1
3609 DO 4 K1=1,3
3610 K2=MAX(I2,J2)+KL+K1-1+ML
3611 4 PSJINT1=PSJINT1+CSJ(I2,J2,K2)*WI(I1)*WJ(J1)*WK(K1)
3612 PSJINT1=EXP(PSJINT1)
3613 ELSE
3614 K=15
3615 IF(I.GT.K-3)I=K-3
3616 IF(J.GT.K-3)J=K-3
3617 WI(2)=QLI-I
3618 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
3619 WI(1)=1.D0-WI(2)+WI(3)
3620 WI(2)=WI(2)-2.D0*WI(3)
3621 WJ(2)=QLJ-J
3622 WJ(3)=WJ(2)*(WJ(2)-1.D0)*.5D0
3623 WJ(1)=1.D0-WJ(2)+WJ(3)
3624 WJ(2)=WJ(2)-2.D0*WJ(3)
3625 WK(2)=SL-K
3626 WK(1)=1.D0-WK(2)
3627
3628 DO 5 I1=1,3
3629 DO 5 J1=1,3
3630 DO 5 K1=1,2
3631 5 PSJINT1=PSJINT1+CSJ(I+I1,J+J1,K+K1+ML)*WI(I1)*WJ(J1)*WK(K1)
3632 PSJINT1=EXP(PSJINT1)
3633 ENDIF
3634 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSJINT1
3635 RETURN
3636 END
3637
3638
3639 FUNCTION QGSPSLAM(S,A,B)
3640
3641
3642
3643
3644
3645 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3646 INTEGER DEBUG
3647 COMMON /Q_AREA43/ MONIOU
3648 COMMON /Q_DEBUG/ DEBUG
3649 SAVE
3650
3651 IF(DEBUG.GE.2)WRITE (MONIOU,201)S,A,B
3652 201 FORMAT(2X,'QGSPSLAM - KINEMATICAL FUNCTION S=',E10.3,2X,'A=',
3653 * E10.3,2X,'B=',E10.3)
3654 QGSPSLAM=.25D0/S*(S+A-B)**2-A
3655 IF(DEBUG.GE.3)WRITE (MONIOU,202)QGSPSLAM
3656 202 FORMAT(2X,'QGSPSLAM=',E10.3)
3657 RETURN
3658 END
3659
3660
3661 FUNCTION QGSPSNORM(EP)
3662
3663
3664 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3665 INTEGER DEBUG
3666 DIMENSION EP(4)
3667 COMMON /Q_AREA43/ MONIOU
3668 COMMON /Q_DEBUG/ DEBUG
3669 SAVE
3670
3671 IF(DEBUG.GE.2)WRITE (MONIOU,201)EP
3672 201 FORMAT(2X,'QGSPSNORM - 4-VECTOR SQUARED FOR ',
3673 * 'EP=',4(E10.3,1X))
3674 QGSPSNORM=EP(1)**2
3675 DO 1 I=1,3
3676 1 QGSPSNORM=QGSPSNORM-EP(I+1)**2
3677 IF(DEBUG.GE.3)WRITE (MONIOU,202)QGSPSNORM
3678 202 FORMAT(2X,'QGSPSNORM=',E10.3)
3679 RETURN
3680 END
3681
3682
3683 SUBROUTINE PSREC(EP,QV,ZV,QM,IQV,LDAU,LPAR,IQJ,EQJ,JFL,JQ)
3684
3685
3686
3687
3688
3689
3690
3691
3692
3693
3694
3695
3696
3697
3698 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3699 INTEGER DEBUG
3700 DIMENSION EP(4),EP3(4),EPV(4,30,50),
3701 * QV(30,50),ZV(30,50),QM(30,50),IQV(30,50),
3702 * LDAU(30,49),LPAR(30,50),
3703 * IQJ(2),EQJ(4,2),IPQ(2,30,50),EPQ(8,30,50),
3704 * EPJ(4),EPJ1(4)
3705 COMMON /Q_AREA43/ MONIOU
3706 COMMON /Q_DEBUG/ DEBUG
3707 SAVE
3708
3709 IF(DEBUG.GE.2)WRITE (MONIOU,201)JQ,EP,IQJ
3710 201 FORMAT(2X,'PSREC - JET RECONSTRUCTURING: JQ=',I1/
3711 * 4X,'JET 4-MOMENTUM EP=',4(E10.3,1X)/4X,'IQJ=',2I2)
3712 JFL = 1
3713 DO 1 I=1,4
3714 EPV(I,1,1)=EP(I)
3715 1 EPQ(I,1,1)=EQJ(I,1)
3716 IPQ(1,1,1)=IQJ(1)
3717
3718 IF(IQV(1,1).EQ.0)THEN
3719 DO 2 I=1,4
3720 2 EPQ(I+4,1,1)=EQJ(I,2)
3721 IPQ(2,1,1)=IQJ(2)
3722 ENDIF
3723
3724 NLEV=1
3725 NROW=1
3726
3727 3 CONTINUE
3728
3729 IF(QV(NROW,NLEV).EQ.0.D0)THEN
3730 IPJ=IQV(NROW,NLEV)
3731 IF(IPJ.NE.0)THEN
3732 IF(EPQ(1,NROW,NLEV).NE.0.D0)THEN
3733 IF(IABS(IPJ).EQ.3)IPJ=IPJ*4/3
3734 DO 4 I=1,4
3735 EPJ(I)=EPV(I,NROW,NLEV)
3736 4 EPJ1(I)=EPQ(I,NROW,NLEV)
3737 IPJ1=IPQ(1,NROW,NLEV)
3738 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
3739 CALL PSJDEF(IPJ,IPJ1,EPJ,EPJ1,JFL)
3740 IF(DEBUG.GE.3)WRITE (MONIOU,211)IPJ,IPJ1,JFL
3741 211 FORMAT(2X,'PSREC - NEW STRING FLAVOURS: ',2I3,' JFL=',I1)
3742 IF(JFL.EQ.0)RETURN
3743 ELSE
3744 IPQ(1,NROW,NLEV)=IPJ
3745 DO 5 I=1,4
3746 5 EPQ(I,NROW,NLEV)=EPV(I,NROW,NLEV)
3747 IF(DEBUG.GE.3)WRITE (MONIOU,212)IPJ,
3748 * (EPV(I,NROW,NLEV),I=1,4),JFL
3749 212 FORMAT(2X,'PSREC: NEW FINAL JET FLAVOR: ',I3,2X,
3750 * 'JET 4-MOMENTUM:', 4(E10.3,1X),' JFL=',I1)
3751 ENDIF
3752
3753 ELSE
3754 IPJ=INT(2.D0*PSRAN(B10)+1.D0)*(3-2*JQ)
3755 DO 6 I=1,4
3756 6 EPJ(I)=.5D0*EPV(I,NROW,NLEV)
3757
3758 DO 9 M=1,2
3759 IF(EPQ(1+4*(M-1),NROW,NLEV).NE.0.D0)THEN
3760 DO 7 I=1,4
3761 7 EPJ1(I)=EPQ(4*(M-1)+I,NROW,NLEV)
3762 IPJ1=IPQ(M,NROW,NLEV)
3763 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
3764 CALL PSJDEF(IPJ,IPJ1,EPJ,EPJ1,JFL)
3765 IF(JFL.EQ.0)RETURN
3766 ELSE
3767 IPQ(M,NROW,NLEV)=IPJ
3768 DO 8 I=1,4
3769 8 EPQ(4*(M-1)+I,NROW,NLEV)=EPJ(I)
3770 ENDIF
3771 9 IPJ=-IPJ
3772 ENDIF
3773
3774 IF(DEBUG.GE.3)WRITE (MONIOU,204)NLEV,NROW,IQV(NROW,NLEV),
3775 * (EPV(I,NROW,NLEV),I=1,4)
3776 204 FORMAT(2X,'PSREC: FINAL JET AT LEVEL NLEV=',I2,
3777 * ' NROW=',I2/4X,'JET FLAVOR: ',I3,2X,'JET 4-MOMENTUM:',
3778 * 4(E10.3,1X))
3779 ELSE
3780
3781 DO 10 I=1,4
3782 10 EP3(I)=EPV(I,NROW,NLEV)
3783 CALL QGSPSDEFROT(EP3,S0X,C0X,S0,C0)
3784 Z=ZV(NROW,NLEV)
3785 QT2=(Z*(1.D0-Z))**2*QV(NROW,NLEV)
3786 LDROW=LDAU(NROW,NLEV)
3787
3788 WP0=EP3(1)+EP3(2)
3789 WPI=Z*WP0
3790 WMI=(QT2+QM(LDROW,NLEV+1))/WPI
3791 EP3(1)=.5D0*(WPI+WMI)
3792 EP3(2)=.5D0*(WPI-WMI)
3793 QT=DSQRT(QT2)
3794 CALL QGSPSCS(C,S)
3795 EP3(3)=QT*C
3796 EP3(4)=QT*S
3797 CALL QGSPSROTAT(EP3,S0X,C0X,S0,C0)
3798
3799 DO 11 I=1,4
3800 11 EPV(I,LDROW,NLEV+1)=EP3(I)
3801 IF(DEBUG.GE.3)WRITE (MONIOU,206)NLEV+1,LDROW,EP3
3802 206 FORMAT(2X,'PSREC: JET AT LEVEL NLEV=',I2,
3803 * ' NROW=',I2/4X,'JET 4-MOMENTUM:',4(E10.3,1X))
3804
3805 WPI=(1.D0-Z)*WP0
3806 WMI=(QT2+QM(LDROW+1,NLEV+1))/WPI
3807 EP3(1)=.5D0*(WPI+WMI)
3808 EP3(2)=.5D0*(WPI-WMI)
3809 EP3(3)=-QT*C
3810 EP3(4)=-QT*S
3811 CALL QGSPSROTAT(EP3,S0X,C0X,S0,C0)
3812 IF(DEBUG.GE.3)WRITE (MONIOU,206)NLEV+1,LDROW+1,EP3
3813
3814 DO 12 I=1,4
3815 12 EPV(I,LDROW+1,NLEV+1)=EP3(I)
3816
3817 IF(IQV(NROW,NLEV).EQ.0)THEN
3818 IF(IQV(LDROW,NLEV+1).NE.0)THEN
3819 IPQ(1,LDROW,NLEV+1)=IPQ(1,NROW,NLEV)
3820 IPQ(1,LDROW+1,NLEV+1)=IPQ(2,NROW,NLEV)
3821 DO 13 I=1,4
3822 EPQ(I,LDROW,NLEV+1)=EPQ(I,NROW,NLEV)
3823 13 EPQ(I,LDROW+1,NLEV+1)=EPQ(I+4,NROW,NLEV)
3824 ELSE
3825 IPQ(1,LDROW,NLEV+1)=IPQ(1,NROW,NLEV)
3826 IPQ(2,LDROW,NLEV+1)=0
3827 IPQ(1,LDROW+1,NLEV+1)=0
3828 IPQ(2,LDROW+1,NLEV+1)=IPQ(2,NROW,NLEV)
3829 DO 14 I=1,4
3830 EPQ(I,LDROW,NLEV+1)=EPQ(I,NROW,NLEV)
3831 EPQ(I+4,LDROW,NLEV+1)=0.D0
3832 EPQ(I,LDROW+1,NLEV+1)=0.D0
3833 14 EPQ(I+4,LDROW+1,NLEV+1)=EPQ(I+4,NROW,NLEV)
3834 ENDIF
3835 ELSE
3836 IF(IQV(LDROW,NLEV+1).EQ.0)THEN
3837 IPQ(1,LDROW,NLEV+1)=IPQ(1,NROW,NLEV)
3838 IPQ(2,LDROW,NLEV+1)=0
3839 IPQ(1,LDROW+1,NLEV+1)=0
3840 DO 15 I=1,4
3841 EPQ(I,LDROW,NLEV+1)=EPQ(I,NROW,NLEV)
3842 EPQ(I+4,LDROW,NLEV+1)=0.D0
3843 15 EPQ(I,LDROW+1,NLEV+1)=0.D0
3844 ELSE
3845 IPQ(1,LDROW,NLEV+1)=0
3846 IPQ(1,LDROW+1,NLEV+1)=0
3847 IPQ(2,LDROW+1,NLEV+1)=IPQ(1,NROW,NLEV)
3848 DO 16 I=1,4
3849 EPQ(I,LDROW,NLEV+1)=0.D0
3850 EPQ(I,LDROW+1,NLEV+1)=0.D0
3851 16 EPQ(I+4,LDROW+1,NLEV+1)=EPQ(I,NROW,NLEV)
3852 ENDIF
3853 ENDIF
3854
3855 NROW=LDROW
3856 NLEV=NLEV+1
3857 GOTO 3
3858 ENDIF
3859
3860 17 CONTINUE
3861 IF(NLEV.EQ.1)THEN
3862 IQJ(1)=IPQ(1,1,1)
3863 DO 18 I=1,4
3864 18 EQJ(I,1)=EPQ(I,1,1)
3865 IF(IQV(1,1).EQ.0)THEN
3866 IQJ(2)=IPQ(2,1,1)
3867 DO 19 I=1,4
3868 19 EQJ(I,2)=EPQ(I+4,1,1)
3869 ENDIF
3870 IF(DEBUG.GE.3)WRITE (MONIOU,202)iqj
3871 202 FORMAT(2X,'PSREC - END',2x,'iqj=',2i2)
3872 RETURN
3873 ENDIF
3874
3875 LPROW=LPAR(NROW,NLEV)
3876
3877 IF(LDAU(LPROW,NLEV-1).EQ.NROW)THEN
3878 IF(IQV(NROW,NLEV).EQ.0)THEN
3879 IF(EPQ(1,LPROW,NLEV-1).EQ.0.D0)THEN
3880 IPQ(1,LPROW,NLEV-1)=IPQ(1,NROW,NLEV)
3881 DO 20 I=1,4
3882 20 EPQ(I,LPROW,NLEV-1)=EPQ(I,NROW,NLEV)
3883 ENDIF
3884 IPQ(1,NROW+1,NLEV)=IPQ(2,NROW,NLEV)
3885 DO 21 I=1,4
3886 21 EPQ(I,NROW+1,NLEV)=EPQ(I+4,NROW,NLEV)
3887 ELSE
3888 IF(IQV(LPROW,NLEV-1).EQ.0)THEN
3889 IF(EPQ(1,LPROW,NLEV-1).EQ.0.D0)THEN
3890 IPQ(1,LPROW,NLEV-1)=IPQ(1,NROW,NLEV)
3891 DO 22 I=1,4
3892 22 EPQ(I,LPROW,NLEV-1)=EPQ(I,NROW,NLEV)
3893 ENDIF
3894 ELSE
3895 IPQ(1,NROW+1,NLEV)=IPQ(1,NROW,NLEV)
3896 DO 23 I=1,4
3897 23 EPQ(I,NROW+1,NLEV)=EPQ(I,NROW,NLEV)
3898 ENDIF
3899 ENDIF
3900 NROW=NROW+1
3901 GOTO 3
3902
3903 ELSE
3904 IF(IQV(NROW,NLEV).EQ.0)THEN
3905 IF(IQV(LPROW,NLEV-1).EQ.0)THEN
3906 IF(EPQ(5,LPROW,NLEV-1).EQ.0.D0)THEN
3907 IPQ(2,LPROW,NLEV-1)=IPQ(2,NROW,NLEV)
3908 DO 24 I=1,4
3909 24 EPQ(I+4,LPROW,NLEV-1)=EPQ(I+4,NROW,NLEV)
3910 ENDIF
3911 ELSE
3912 IF(EPQ(1,LPROW,NLEV-1).EQ.0.D0)THEN
3913 IPQ(1,LPROW,NLEV-1)=IPQ(2,NROW,NLEV)
3914 DO 25 I=1,4
3915 25 EPQ(I,LPROW,NLEV-1)=EPQ(I+4,NROW,NLEV)
3916 ENDIF
3917 ENDIF
3918 ELSE
3919 IF(IQV(LPROW,NLEV-1).EQ.0.AND.
3920 * EPQ(5,LPROW,NLEV-1).EQ.0.D0)THEN
3921 IPQ(2,LPROW,NLEV-1)=IPQ(1,NROW,NLEV)
3922 DO 26 I=1,4
3923 26 EPQ(I+4,LPROW,NLEV-1)=EPQ(I,NROW,NLEV)
3924 ENDIF
3925 ENDIF
3926
3927 NROW=LPROW
3928 NLEV=NLEV-1
3929 GOTO 17
3930 ENDIF
3931 END
3932
3933
3934 FUNCTION PSREJS(S,Z,IQQ)
3935
3936
3937
3938
3939
3940
3941 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
3942 INTEGER DEBUG
3943 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALF,RR,SH,DELH
3944 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
3945 COMMON /Q_AR13/ X1(7),A1(7)
3946 COMMON /Q_AREA43/ MONIOU
3947 COMMON /Q_DEBUG/ DEBUG
3948 SAVE
3949
3950 IF(DEBUG.GE.2)WRITE (MONIOU,201)S,Z,IQQ
3951 201 FORMAT(2X,'PSREJS - REJECTION FUNCTION TABULATION: '/
3952 * 4X,'S=',E10.3,2X,'Z=',E10.3,2X,'IQQ=',I1)
3953 XMIN=4.D0*(QT0+AMJ0)/S
3954 XMIN=XMIN**(DELH-DEL)
3955 PSREJS=0.D0
3956
3957
3958 DO 2 I=1,7
3959 DO 2 M=1,2
3960 Z1=(.5D0*(1.D0+XMIN-(2*M-3)*X1(I)*(1.D0-XMIN)))**(1.D0/
3961 *(DELH-DEL))
3962
3963
3964
3965
3966
3967 YJ=DLOG(Z1*S)
3968 CALL PSJINT0(Z1*S,SJ,SJB,IQQ,0)
3969
3970
3971
3972 GY=2.D0*SH*PSGINT((SJ-SJB)/SH*.5D0)+SJB
3973 RH=RS0-ALF*DLOG(Z1)
3974
3975 IF(IQQ.NE.0)THEN
3976 PSREJS=PSREJS+A1(I)*GY/(Z1*S)**DELH*Z**(RS0/RH)/RH*
3977 * (1.D0-Z1)*BET
3978 ELSE
3979 ST2=0.D0
3980 DO 1 J=1,7
3981 1 ST2=ST2+A1(J)*((1.D0-Z1**(.5D0*(1.D0+X1(J))))*
3982 * (1.D0-Z1**(.5D0*(1.D0-X1(J)))))**BET
3983
3984 PSREJS=PSREJS-A1(I)*DLOG(Z1)*GY/(Z1*S)**DELH*Z**(RS0/RH)/RH*ST2
3985 ENDIF
3986 2 CONTINUE
3987 PSREJS=DLOG(PSREJS*(1.D0-XMIN)/Z)
3988 IF(DEBUG.GE.2)WRITE (MONIOU,202)PSREJS
3989 202 FORMAT(2X,'PSREJS=',E10.3)
3990 RETURN
3991 END
3992
3993
3994 FUNCTION PSREJV(S)
3995
3996
3997
3998
3999 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
4000 INTEGER DEBUG
4001 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALF,RR,SH,DELH
4002 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
4003 COMMON /Q_AREA43/ MONIOU
4004 COMMON /Q_DEBUG/ DEBUG
4005 SAVE
4006
4007 IF(DEBUG.GE.2)WRITE (MONIOU,201)S
4008 201 FORMAT(2X,'PSREJV - REJECTION FUNCTION TABULATION: ',
4009 * 'S=',E10.3)
4010
4011
4012
4013
4014 CALL PSJINT0(S,SJ,SJB,1,1)
4015
4016
4017
4018
4019 GY=2.D0*SH*PSGINT((SJ-SJB)/SH*.5D0)+SJB
4020 PSREJV=DLOG(GY/S**DELH)
4021 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSREJV
4022 202 FORMAT(2X,'PSREJV=',E10.3)
4023 RETURN
4024 END
4025
4026
4027 FUNCTION PSRJINT(YJ,Z0,IQQ)
4028
4029
4030
4031
4032
4033 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
4034 INTEGER DEBUG
4035 DIMENSION A(3)
4036 COMMON /Q_AREA1/ IA(2),ICZ,ICP
4037 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALF,RR,SH,DELH
4038 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
4039 COMMON /Q_AREA23/ RJV(50)
4040 COMMON /Q_AREA24/ RJS(50,5,10)
4041 COMMON /Q_AREA43/ MONIOU
4042 COMMON /Q_DEBUG/ DEBUG
4043 SAVE
4044
4045 IF(DEBUG.GE.2)WRITE (MONIOU,201)YJ,Z0,IQQ
4046 201 FORMAT(2X,'PSRJINT - REJECTION FUNCTION INTERPOLATION:'/
4047 * 4X,'YJ=',E10.3,2X,'Z0=',E10.3,2X,'IQQ=',I1)
4048 YY=(YJ-AQT0)*2.D0
4049
4050 JY=MIN(48,INT(YY)) ! modified 15.oct.03 D.H.
4051
4052 IF(IQQ.EQ.3)THEN
4053 IF(JY.EQ.0)THEN
4054 PSRJINT=EXP(RJV(1))*YY+(EXP(RJV(2))-2.D0*
4055 * EXP(RJV(1)))*YY*(YY-1.D0)*.5D0
4056 ELSE
4057 PSRJINT=EXP(RJV(JY)+(RJV(JY+1)-RJV(JY))*(YY-JY)
4058 * +(RJV(JY+2)+RJV(JY)-2.D0*RJV(JY+1))*(YY-JY)*
4059 * (YY-JY-1.D0)*.5D0)
4060 ENDIF
4061 ELSE
4062 Z=Z0**(RS/RS0)
4063 IQ=(IQQ+1)/2+1+2*(ICZ-1)
4064 JZ=INT(5.D0*Z)
4065 IF(JZ.GT.3)JZ=3
4066 WZ=5.D0*Z-JZ
4067
4068 IF(JZ.EQ.0)THEN
4069 I1=2
4070 ELSE
4071 I1=1
4072 ENDIF
4073
4074 DO 1 I=I1,3
4075 J1=JZ+I-1
4076 IF(JY.EQ.0)THEN
4077 A(I)=EXP(RJS(1,J1,IQ))*YY+(EXP(RJS(2,J1,IQ))-2.D0*
4078 * EXP(RJS(1,J1,IQ)))*YY*(YY-1.D0)*.5D0
4079 IF(A(I).GT.0.D0)THEN
4080 A(I)=DLOG(A(I))
4081 ELSE
4082 A(I)=-80.D0
4083 ENDIF
4084 ELSE
4085 A(I)=RJS(JY,J1,IQ)+(RJS(JY+1,J1,IQ)-
4086 * RJS(JY,J1,IQ))*(YY-JY)
4087 * +(RJS(JY+2,J1,IQ)+RJS(JY,J1,IQ)-2.D0*
4088 * RJS(JY+1,J1,IQ))*(YY-JY)*(YY-JY-1.D0)*.5D0
4089 ENDIF
4090 1 CONTINUE
4091
4092 IF(JZ.NE.0)THEN
4093 PSRJINT=EXP(A(1)+(A(2)-A(1))*WZ+(A(3)+A(1)-2.D0*A(2))*WZ*
4094 * (WZ-1.D0)*.5D0)*Z
4095 ELSE
4096 PSRJINT=(EXP(A(2))*WZ+(EXP(A(3))-2.D0*EXP(A(2)))*WZ*
4097 * (WZ-1.D0)*.5D0)*Z
4098 IF(PSRJINT.LE.0.D0)PSRJINT=1.D-10
4099 ENDIF
4100 ENDIF
4101 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSRJINT
4102 202 FORMAT(2X,'PSRJINT=',E10.3)
4103 RETURN
4104 END
4105
4106
4107 FUNCTION PSROOT(QLMAX,G,J)
4108
4109
4110
4111
4112
4113
4114
4115 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
4116 INTEGER DEBUG
4117 COMMON /Q_AREA43/ MONIOU
4118 COMMON /Q_DEBUG/ DEBUG
4119 SAVE
4120
4121 IF(DEBUG.GE.2)WRITE (MONIOU,201)QLMAX,G,J
4122 201 FORMAT(2X,'PSQINT - BRANCHING MOMENTUM TABULATION:'/
4123 * 4X,'QLMAX=',E10.3,2X,'G=',E10.3,2X,'J=',I1)
4124 QL0=0.D0
4125 QL1=QLMAX
4126 F0=-G
4127 F1=1.D0-G
4128 SUD0=-DLOG(PSUDINT(QLMAX,J))
4129
4130 1 QL2=QL1-(QL1-QL0)*F1/(F1-F0)
4131 IF(QL2.LT.0.D0)THEN
4132 QL2=0.D0
4133 F2=-G
4134 ELSEIF(QL2.GT.QLMAX)THEN
4135 QL2=QLMAX
4136 F2=1.D0-G
4137 ELSE
4138 F2=-DLOG(PSUDINT(QL2,J))/SUD0-G
4139 ENDIF
4140
4141 IF(ABS(F2).GT.1.D-3)THEN
4142 QL0=QL1
4143 QL1=QL2
4144 F0=F1
4145 F1=F2
4146 GOTO 1
4147 ELSE
4148 PSROOT=QL2
4149 ENDIF
4150 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSROOT
4151 202 FORMAT(2X,'PSROOT=',E10.3)
4152 RETURN
4153 END
4154
4155
4156 SUBROUTINE QGSPSROTAT(EP,S0X,C0X,S0,C0)
4157
4158
4159 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
4160 INTEGER DEBUG
4161 DIMENSION EP(4),EP1(3)
4162 COMMON /Q_AREA43/ MONIOU
4163 COMMON /Q_DEBUG/ DEBUG
4164 SAVE
4165
4166 IF(DEBUG.GE.2)WRITE (MONIOU,201)EP,S0X,C0X,S0,C0
4167 201 FORMAT(2X,'QGSPSROTAT - SPACIAL ROTATION:'/4X,
4168 * '4-VECTOR EP=',4(E10.3,1X)/4X,'S0X=',E10.3,'C0X=',E10.3,
4169 * 2X,'S0=',E10.3,'C0=',E10.3)
4170 EP1(3)=EP(4)
4171 EP1(2)=EP(2)*S0+EP(3)*C0
4172 EP1(1)=EP(2)*C0-EP(3)*S0
4173
4174 EP(2)=EP1(1)
4175 EP(4)=EP1(2)*S0X+EP1(3)*C0X
4176 EP(3)=EP1(2)*C0X-EP1(3)*S0X
4177 IF(DEBUG.GE.3)WRITE (MONIOU,202)EP
4178 202 FORMAT(2X,'QGSPSROTAT: ROTATED 4-VECTOR EP=',
4179 * 2X,4E10.3)
4180 RETURN
4181 END
4182
4183
4184 FUNCTION PSQINT(QLMAX,G,J)
4185
4186
4187
4188
4189
4190
4191 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
4192 INTEGER DEBUG
4193 DIMENSION WI(3),WK(3)
4194 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
4195 COMMON /Q_AREA34/ QRT(10,101,2)
4196 COMMON /Q_AREA43/ MONIOU
4197 COMMON /Q_DEBUG/ DEBUG
4198 SAVE
4199
4200 IF(DEBUG.GE.2)WRITE (MONIOU,201)QLMAX,G,J
4201 201 FORMAT(2X,'PSQINT - BRANCHING MOMENTUM INTERPOLATION:'/
4202 * 4X,'QLMAX=',E10.3,2X,'G=',E10.3,2X,'J=',I1)
4203 QLI=QLMAX/1.38629d0
4204 SUD0=1.D0/PSUDINT(QLMAX,J)
4205 SL=100.D0*DLOG(1.D0-G*(1.D0-SUD0))/DLOG(SUD0)
4206 I=INT(QLI)
4207 K=INT(SL)
4208 IF(K.GT.98)K=98
4209 WK(2)=SL-K
4210 WK(3)=WK(2)*(WK(2)-1.D0)*.5D0
4211 WK(1)=1.D0-WK(2)+WK(3)
4212 WK(2)=WK(2)-2.D0*WK(3)
4213 PSQINT=0.D0
4214
4215 IF(I.GT.7)I=7
4216 WI(2)=QLI-I
4217 WI(3)=WI(2)*(WI(2)-1.D0)*.5D0
4218 WI(1)=1.D0-WI(2)+WI(3)
4219 WI(2)=WI(2)-2.D0*WI(3)
4220
4221 DO 1 K1=1,3
4222 DO 1 I1=1,3
4223 1 PSQINT=PSQINT+QRT(I+I1,K+K1,J)*WI(I1)*WK(K1)
4224 IF(PSQINT.LE.0.D0)PSQINT=0.D0
4225 PSQINT=16.D0*QTF*EXP(PSQINT)
4226 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSQINT
4227 202 FORMAT(2X,'PSQINT=',E10.3)
4228 RETURN
4229 END
4230
4231
4232 SUBROUTINE PSSHAR(LS,NHP,NW,NT)
4233
4234
4235
4236
4237
4238
4239 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
4240 INTEGER DEBUG
4241
4242 DIMENSION WP(56),WM(56),WHA(4000),WHB(4000),LHA0(56),
4243 * LHB0(56),IZP(56),IZT(56),WP0H(56),WM0H(56),
4244 * WPP(2),WMM(2),EP3(4),LQA0(56),LQB0(56),IPC(2,2),EPC(8,2),
4245 * ILA(56),ILB(56),ELA(4,56),ELB(4,56),EP(4),EP1(4)
4246 COMMON /Q_AREA1/ IA(2),ICZ,ICP
4247 COMMON /Q_AREA2/ S,Y0,WP0,WM0
4248 COMMON /Q_AREA9/ LQA(56),LQB(56),NQS(1000),IAS(1000),
4249 * IBS(1000),LHA(56),LHB(56),ZH(4000),IAH(4000),IBH(4000),
4250 * IQH(4000),LVA(56),LVB(56)
4251 COMMON /Q_AREA10/ STMASS,AM(7)
4252 COMMON /Q_AREA11/ B10
4253 COMMON /Q_AREA12/ NSH
4254 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALFP,RR,SH,DELH
4255 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
4256 COMMON /Q_AREA19/ AHL(5)
4257 COMMON /Q_AREA20/ WPPP
4258 COMMON /Q_AREA25/ AHV(5)
4259 COMMON /Q_AREA43/ MONIOU
4260 COMMON /Q_DEBUG/ DEBUG
4261 COMMON /Q_AREA47/ NJTOT
4262
4263
4264 integer ng1evt,ng2evt,ikoevt
4265 real rglevt,sglevt,eglevt,fglevt,typevt
4266 common/c2evt/ng1evt,ng2evt,rglevt,sglevt,eglevt,fglevt,ikoevt
4267 *,typevt !in epos.inc
4268
4269 SAVE
4270 EXTERNAL PSRAN
4271 DATA xdiv /100.D0/
4272
4273 IF(DEBUG.GE.1)WRITE (MONIOU,201)NW,NT,NHP,LS
4274 201 FORMAT(2X,'PSSHARE - ENERGY SHARING PROCEDURE'/
4275 * 4X,'NUMBER OF WOUNDED PROJECTILE NUCLEONS(HADRONS) NW=',I2/
4276 * 4X,'NUMBER OF TARGET NUCLEONS IN THE ACTIVE STATE NT=',I2/
4277 * 4X,'NUMBER OF SEMIHARD BLOCKS NHP=',I3/
4278 * 4X,'NUMBER OF SOFT POMERON BLOCKS LS=',I3)
4279 NSH1=NSH
4280 DO 101 I=1,NW
4281 101 LQA0(I)=LQA(I)
4282 DO 102 I=1,NT
4283 102 LQB0(I)=LQB(I)
4284
4285 100 NSH=NSH1
4286 NJTOT=0
4287 DO 103 I=1,NW
4288 103 LQA(I)=LQA0(I)
4289 DO 104 I=1,NT
4290 104 LQB(I)=LQB0(I)
4291
4292
4293 IF(IA(1).NE.1)THEN
4294
4295 DO 1 I=1,NW
4296 1 IZP(I)=INT(2.5+PSRAN(B10))
4297 ELSE
4298
4299 IZP(1)=ICP
4300 ENDIF
4301 IF(IA(2).NE.1)THEN
4302
4303 DO 2 I=1,NT
4304 2 IZT(I)=INT(2.5+PSRAN(B10))
4305 ELSE
4306
4307 IZT(1)=2
4308 ENDIF
4309
4310
4311
4312 WREJ=.0001D0
4313
4314 3 CONTINUE
4315
4316 IF(NHP.NE.0)THEN
4317 IF(DEBUG.GE.3)WRITE (MONIOU,211)NHP
4318 211 FORMAT(2X,'PSSHARE: NUMBER OF HARD POMERONS NHP=',I3)
4319
4320
4321
4322
4323
4324
4325
4326
4327
4328
4329 GBH0=.6D0
4330
4331 NREJ=0
4332 NHP1=NHP
4333
4334 DO 5 IH=1,NHP1
4335 IF(DEBUG.GE.3)WRITE (MONIOU,212)IH
4336 212 FORMAT(2X,'PSSHARE: GBH-INI; CONTRIBUTION FROM ',I3,
4337 * '-TH HARD POMERON')
4338
4339
4340
4341
4342
4343
4344
4345 IQQ=IQH(IH)
4346 Z=ZH(IH)
4347 I=IAH(IH)
4348 J=IBH(IH)
4349
4350
4351 ZA=1.D0/LHA(I)
4352 ZB=1.D0/LHB(J)
4353
4354 SI=ZA*ZB*S
4355
4356 IF(SI.LT.4.D0*(QT0+AMJ0))THEN
4357
4358
4359
4360
4361 NHP=NHP-1
4362 LHA(I)=LHA(I)-1
4363 LHB(J)=LHB(J)-1
4364
4365 IF(IQQ.EQ.1)THEN
4366 LVA(I)=0
4367 ELSEIF(IQQ.EQ.2)THEN
4368 LVB(J)=0
4369 ELSEIF(IQQ.EQ.3)THEN
4370 LVA(I)=0
4371 LVB(J)=0
4372 ENDIF
4373
4374 IF(NHP.GE.IH)THEN
4375 DO 4 IH1=IH,NHP
4376 IQH(IH1)=IQH(IH1+1)
4377 ZH(IH1)=ZH(IH1+1)
4378 IAH(IH1)=IAH(IH1+1)
4379 4 IBH(IH1)=IBH(IH1+1)
4380 ENDIF
4381
4382
4383 GOTO 3
4384 ENDIF
4385
4386
4387 YI=DLOG(SI)
4388 IF(YI.GT.17.D0)YI=17.D0
4389
4390 GBH0=GBH0/PSRJINT(YI,Z,IQQ)
4391 GBH0 = GBH0/xdiv
4392 5 CONTINUE
4393 IF(DEBUG.GE.3)WRITE (MONIOU,213)
4394 213 FORMAT(2X,'PSSHARE: GBH-INI - END')
4395
4396
4397
4398
4399
4400
4401
4402
4403
4404
4405
4406
4407 6 DO 7 I=1,NW
4408 LHA0(I)=LHA(I)-LVA(I)
4409 7 WP(I)=WP0
4410
4411 DO 8 I=1,NT
4412 LHB0(I)=LHB(I)-LVB(I)
4413 8 WM(I)=WM0
4414
4415
4416
4417
4418 DO 10 I=1,NW
4419 IF(LVA(I).NE.0)THEN
4420 9 XW=PSRAN(B10)**(1.D0/(.5D0+DELH))
4421 IF(PSRAN(B10).GT.(1.D0-XW)**AHV(ICZ))GOTO 9
4422 IF(DEBUG.GE.3)WRITE (MONIOU,214)I,XW
4423 214 FORMAT(2X,'PSSHARE: ',I2,'-TH PROJ. NUCLEON (HADRON); LIGHT',
4424 * ' CONE MOMENTUM SHARE XW=',E10.3)
4425
4426 WP0H(I)=XW*WP(I)
4427
4428 WP(I)=WP(I)*(1.D0-XW)
4429 ENDIF
4430 10 CONTINUE
4431
4432
4433
4434 DO 12 I=1,NT
4435 IF(LVB(I).NE.0)THEN
4436 11 XW=PSRAN(B10)**(1.D0/(.5D0+DELH))
4437 IF(PSRAN(B10).GT.(1.D0-XW)**AHV(2))GOTO 11
4438 IF(DEBUG.GE.3)WRITE (MONIOU,215)I,XW
4439 215 FORMAT(2X,'PSSHARE: ',I2,'-TH TARGET NUCLEON (HADRON); LIGHT',
4440 * ' CONE MOMENTUM SHARE XW=',E10.3)
4441
4442 WM0H(I)=XW*WM(I)
4443
4444 WM(I)=WM(I)*(1.D0-XW)
4445 ENDIF
4446 12 CONTINUE
4447
4448
4449 GBH=GBH0
4450
4451
4452
4453 DO 18 IH=1,NHP1
4454
4455
4456
4457
4458
4459 IQQ=IQH(IH)
4460 Z=ZH(IH)
4461 I=IAH(IH)
4462 J=IBH(IH)
4463
4464 IF((IQQ-3)*(IQQ-1).EQ.0)THEN
4465
4466
4467 WHA(IH)=WP0H(I)
4468 ELSE
4469
4470 LHA0(I)=LHA0(I)-1
4471
4472
4473
4474
4475 BPI=1.D0/(1.D0+AHL(ICZ)+
4476 * (1.D0+DELH)*LHA0(I))
4477
4478
4479 15 XW=1.-PSRAN(B10)**BPI
4480
4481 IF(PSRAN(B10).GT.XW**DELH)GOTO 15
4482
4483 WHA(IH)=WP(I)*XW
4484
4485 WP(I)=WP(I)*(1.D0-XW)
4486 ENDIF
4487
4488 IF((IQQ-3)*(IQQ-2).EQ.0)THEN
4489
4490
4491 WHB(IH)=WM0H(J)
4492 ELSE
4493
4494 LHB0(J)=LHB0(J)-1
4495 BPI=1.D0/(1.D0+AHL(2)+(1.D0+DELH)
4496 * *LHB0(J))
4497
4498
4499 16 XW=1.-PSRAN(B10)**BPI
4500 IF(PSRAN(B10).GT.XW**DELH)GOTO 16
4501
4502 WHB(IH)=WM(J)*XW
4503
4504 WM(J)=WM(J)*(1.D0-XW)
4505 ENDIF
4506
4507
4508 SW=WHA(IH)*WHB(IH)
4509 IF(SW.LT.4.D0*(QT0+AMJ0))THEN
4510
4511 NREJ=NREJ+1
4512
4513 IF(NREJ.GT.30)THEN
4514
4515
4516
4517
4518 NHP=NHP-1
4519 LHA(I)=LHA(I)-1
4520 LHB(J)=LHB(J)-1
4521
4522 IF(IQQ.EQ.1)THEN
4523 LVA(I)=0
4524 ELSEIF(IQQ.EQ.2)THEN
4525 LVB(J)=0
4526 ELSEIF(IQQ.EQ.3)THEN
4527 LVA(I)=0
4528 LVB(J)=0
4529 ENDIF
4530
4531 IF(NHP.GE.IH)THEN
4532 DO 17 IH1=IH,NHP
4533 IQH(IH1)=IQH(IH1+1)
4534 ZH(IH1)=ZH(IH1+1)
4535 IAH(IH1)=IAH(IH1+1)
4536 17 IBH(IH1)=IBH(IH1+1)
4537 ENDIF
4538 GOTO 3
4539
4540
4541
4542
4543 ELSE
4544 GOTO 6
4545 ENDIF
4546 ENDIF
4547 IF(DEBUG.GE.3)WRITE (MONIOU,216)IH,WHA(IH),WHB(IH),WP(I),WM(J)
4548 216 FORMAT(2X,'PSSHARE: ',I3,'-TH SEMIHARD BLOCK; LIGHT',
4549 * ' CONE MOMENTA SHARES:',2E10.3/
4550 * 4X,'REMAINED LIGHT CONE MOMENTA:',2E10.3)
4551
4552 YH=DLOG(SW)
4553
4554
4555
4556
4557
4558 GBH=GBH*PSRJINT(YH,Z,IQQ)
4559 GBH = GBH * xdiv
4560 18 CONTINUE
4561
4562
4563
4564
4565
4566
4567 IF(DEBUG.GE.2)WRITE (MONIOU,217)1.D0-GBH,NHP
4568 217 FORMAT(2X,'PSSHARE: REJECTION PROBABILITY:',E10.3,
4569 * 2X,'NUMBER OF SEMIHARD BLOCKS:',I3)
4570 IF(PSRAN(B10).GT.GBH)THEN
4571 NREJ=NREJ+1
4572
4573 IF(NREJ.GT.30)THEN
4574 IF(DEBUG.GE.2)WRITE (MONIOU,218)
4575 218 FORMAT(2X,'PSSHARE: MORE THAN 30 REJECTIONS - HARD POMERON',
4576 * ' NUMBER IS PUT DOWN')
4577
4578
4579
4580
4581 LNH=1+NHP/20
4582 DO 19 IHP=NHP-LNH+1,NHP
4583 IIH=IAH(IHP)
4584 JIH=IBH(IHP)
4585 IQQ=IQH(IHP)
4586
4587 IF(IQQ.EQ.1)THEN
4588 LVA(IIH)=0
4589 ELSEIF(IQQ.EQ.2)THEN
4590 LVB(JIH)=0
4591 ELSEIF(IQQ.EQ.3)THEN
4592 LVA(IIH)=0
4593 LVB(JIH)=0
4594 ENDIF
4595
4596 LHA(IIH)=LHA(IIH)-1
4597 19 LHB(JIH)=LHB(JIH)-1
4598
4599 NHP=NHP-LNH
4600 GOTO 3
4601
4602
4603
4604 ELSE
4605 GOTO 6
4606 ENDIF
4607 ENDIF
4608
4609
4610 DO 31 I=1,NW
4611 31 LHA0(I)=LHA(I)
4612 DO 32 I=1,NT
4613 32 LHB0(I)=LHB(I)
4614
4615
4616
4617
4618
4619 DO 20 IH=1,NHP
4620 IQQ=IQH(IH)
4621 Z=ZH(IH)
4622 I=IAH(IH)
4623 J=IBH(IH)
4624
4625 LHA0(I)=LHA0(I)-1
4626 LHB0(J)=LHB0(J)-1
4627
4628
4629 WPI=WHA(IH)
4630 WMI=WHB(IH)
4631 IF(DEBUG.GE.2)WRITE (MONIOU,219)IH,IQQ,WPI,WMI,WP(I),WM(J)
4632 219 FORMAT(2X,'PSSHARE: ',I3,
4633 * '-TH HARD BLOCK; TYPE OF THE INTERACTION:',I1/
4634 * 4X,'INITIAL LIGHT CONE MOMENTA:',2E10.3/
4635 * 4X,'REMAINED LIGHT CONE MOMENTA:',2E10.3)
4636
4637
4638
4639 CALL PSHOT(WPI,WMI,Z,IPC,EPC,IZP(I),IZT(J),ICZ,IQQ)
4640 IF(IQQ.EQ.1.OR.IQQ.EQ.3)THEN
4641 IF((IABS(IZP(I)).GT.5.OR.IABS(IZP(I)).EQ.3).AND.
4642 * IZP(I).GT.0.OR.IABS(IZP(I)).NE.3.AND.
4643 * IABS(IZP(I)).LE.5.AND.IZP(I).LT.0)THEN
4644 JQ=1
4645 ELSE
4646 JQ=2
4647 ENDIF
4648 ILA(I)=IPC(JQ,1)
4649 DO 330 L=1,4
4650 330 ELA(L,I)=EPC(L+4*(JQ-1),1)
4651 ENDIF
4652 IF(IQQ.EQ.2.OR.IQQ.EQ.3)THEN
4653 IF((IABS(IZT(J)).GT.5.OR.IABS(IZT(J)).EQ.3).AND.
4654 * IZT(J).GT.0.OR.IABS(IZT(J)).NE.3.AND.
4655 * IABS(IZT(J)).LE.5.AND.IZT(J).LT.0)THEN
4656 JQ=1
4657 ELSE
4658 JQ=2
4659 ENDIF
4660 ILB(J)=IPC(JQ,2)
4661 DO 331 L=1,4
4662 331 ELB(L,J)=EPC(L+4*(JQ-1),2)
4663 ENDIF
4664 IF(IQQ.EQ.3.AND.ILA(I)+ILB(J).EQ.0)NIAS=J
4665
4666
4667
4668
4669
4670
4671
4672
4673
4674
4675
4676
4677 IF(LQA(I)+LHA0(I).EQ.0.AND.LQB(J)+LHB0(J).EQ.0)THEN
4678 IF(LVA(I).EQ.0.AND.LVB(J).EQ.0)THEN
4679 CALL XXDDFR(WP(I),WM(J),IZP(I),IZT(J))
4680 ELSEIF(LVA(I).EQ.0)THEN
4681 CALL XXDPR(WP(I),WM(J),IZP(I),IZT(J),1)
4682 IF(ILB(J).NE.0)THEN
4683 DO 341 L=1,4
4684 341 EP1(L)=ELB(L,J)
4685 EP(1)=.5D0*WM(J)
4686 EP(2)=-EP(1)
4687 EP(3)=0.D0
4688 EP(4)=0.D0
4689 IPJ1=ILB(J)
4690 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
4691 CALL PSJDEF(IZT(J),IPJ1,EP,EP1,JFL)
4692 IF(JFL.EQ.0)GOTO 100
4693 ENDIF
4694 ELSEIF(LVB(J).EQ.0)THEN
4695 CALL XXDTG(WP(I),WM(J),IZP(I),IZT(J),1)
4696 IF(ILA(I).NE.0)THEN
4697 DO 342 L=1,4
4698 342 EP1(L)=ELA(L,I)
4699 EP(1)=.5D0*WP(I)
4700 EP(2)=EP(1)
4701 EP(3)=0.D0
4702 EP(4)=0.D0
4703 IPJ1=ILA(I)
4704 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
4705 CALL PSJDEF(IZP(I),IPJ1,EP,EP1,JFL)
4706 IF(JFL.EQ.0)GOTO 100
4707 ENDIF
4708 ELSE
4709 IF(ILA(I).NE.0)THEN
4710 DO 343 L=1,4
4711 343 EP1(L)=ELA(L,I)
4712 EP(1)=.5D0*WP(I)
4713 EP(2)=EP(1)
4714 EP(3)=0.D0
4715 EP(4)=0.D0
4716 IPJ1=ILA(I)
4717 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
4718 CALL PSJDEF(IZP(I),IPJ1,EP,EP1,JFL)
4719 IF(JFL.EQ.0)GOTO 100
4720 ENDIF
4721 IF(ILB(J).NE.0)THEN
4722 DO 351 L=1,4
4723 351 EP1(L)=ELB(L,J)
4724 EP(1)=.5D0*WM(J)
4725 EP(2)=-EP(1)
4726 EP(3)=0.D0
4727 EP(4)=0.D0
4728 IPJ1=ILB(J)
4729 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
4730 CALL PSJDEF(IZT(J),IPJ1,EP,EP1,JFL)
4731 IF(JFL.EQ.0)GOTO 100
4732 ENDIF
4733 ENDIF
4734 ELSEIF(LQA(I)+LHA0(I).EQ.0)THEN
4735 IF(LVA(I).EQ.0)THEN
4736 CALL XXDPR(WP(I),WM(J),IZP(I),IZT(J),LQB(J)+LHB0(J))
4737 ELSE
4738 IF(ILA(I).NE.0)THEN
4739 DO 344 L=1,4
4740 344 EP1(L)=ELA(L,I)
4741 EP(1)=.5D0*WP(I)
4742 EP(2)=EP(1)
4743 EP(3)=0.D0
4744 EP(4)=0.D0
4745 IPJ1=ILA(I)
4746 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
4747 CALL PSJDEF(IZP(I),IPJ1,EP,EP1,JFL)
4748 IF(JFL.EQ.0)GOTO 100
4749 ENDIF
4750 ENDIF
4751 ELSEIF(LQB(J)+LHB0(J).EQ.0)THEN
4752 IF(LVB(J).EQ.0)THEN
4753 CALL XXDTG(WP(I),WM(J),IZP(I),IZT(J),LQA(I)+LHA0(I))
4754 ELSE
4755 IF(ILB(J).NE.0)THEN
4756 DO 345 L=1,4
4757 345 EP1(L)=ELB(L,J)
4758 EP(1)=.5D0*WM(J)
4759 EP(2)=-EP(1)
4760 EP(3)=0.D0
4761 EP(4)=0.D0
4762 IPJ1=ILB(J)
4763 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
4764 CALL PSJDEF(IZT(J),IPJ1,EP,EP1,JFL)
4765 IF(JFL.EQ.0)GOTO 100
4766 ENDIF
4767 ENDIF
4768 ENDIF
4769
4770 20 CONTINUE
4771
4772
4773
4774
4775 ELSE
4776
4777
4778 DO 21 I=1,NW
4779 21 WP(I)=WP0
4780 DO 22 I=1,NT
4781 22 WM(I)=WM0
4782 ENDIF
4783
4784 IF(LS.NE.0)THEN
4785
4786
4787
4788 DO 28 IS=1,LS
4789
4790
4791
4792
4793
4794
4795
4796
4797
4798
4799 I=IAS(IS)
4800 J=IBS(IS)
4801 LQ1=LQA(I)
4802 LQ2=LQB(J)
4803 WPN=WP(I)
4804 WMN=WM(J)
4805 NP=NQS(IS)
4806 IF(DEBUG.GE.3)WRITE (MONIOU,222)IS,I,J,NP
4807 222 FORMAT(2X,'PSSHARE: ',I3,'-TH SOFT POMERON BLOCK IS',
4808 * ' CONNECTED TO ',I2,
4809 * '-TH PROJECTILE NUCLEON'/4x,'(HADRON) AND ',I2,
4810 * '-TH TARGET NUCLEON'/
4811 * 4X,'NUMBER OF CUT SOFT POMERONS IN THE BLOCK:',I2)
4812
4813
4814 DO 27 IP=1,NP
4815
4816
4817
4818 14 JPP=0
4819
4820 IF(LQ1.EQ.1.AND.WPN.EQ.WP0.AND.PSRAN(B10).LT.WPPP
4821 * .AND.LVB(J).EQ.0)THEN !!!!!!!!!!!!!!!!!!so-07.03.99
4822
4823
4824
4825
4826 YW=1.D0+PSRAN(B10)*(Y0-2.D0)
4827 IF(DEBUG.GE.3)WRITE (MONIOU,223)YW
4828 223 FORMAT(2X,'PSSHARE: TRIPLE POMERON CONTRIBUTION YW=',E10.3)
4829
4830
4831 XPW=EXP(-YW)
4832 JPP=1
4833
4834
4835 ELSE
4836 LQ1=LQ1-1
4837
4838
4839
4840 BPI=1.D0/(1.D0+AHL(ICZ)+(1.D0+DEL)*LQ1)
4841 23 XPW=1.-PSRAN(B10)**BPI
4842
4843 IF(PSRAN(B10).GT.XPW**DEL)GOTO 23
4844 ENDIF
4845
4846 LQ2=LQ2-1
4847
4848
4849 BPI=1.D0/(1.D0+AHL(2)+(1.D0+DEL)*LQ2)
4850 24 XMW=1.-PSRAN(B10)**BPI
4851
4852 IF(PSRAN(B10).GT.XMW**DEL)GOTO 24
4853
4854
4855
4856
4857 IF(JPP.EQ.1.AND.XPW*XMW*WPN*WMN.LT.2.72D0)THEN
4858 LQ2=LQ2+1
4859 GOTO 14
4860 ENDIF
4861
4862
4863
4864
4865 WPI=WPN*XPW
4866 WPN=WPN-WPI
4867 WMI=WMN*XMW
4868 WMN=WMN-WMI
4869
4870
4871
4872 IF(LQ1.EQ.0.AND.LVA(I).EQ.0)THEN
4873 CALL IXXDEF(IZP(I),IC11,IC12,ICZ)
4874 ELSE
4875 IC11=0
4876 IC12=0
4877 ENDIF
4878 IF(LQ2.EQ.0.AND.LVB(J).EQ.0)THEN
4879 CALL IXXDEF(IZT(J),IC21,IC22,2)
4880 ELSE
4881 IC21=0
4882 IC22=0
4883 ENDIF
4884
4885
4886
4887
4888
4889 IF(DEBUG.GE.3)WRITE (MONIOU,224)IP,WPI,WMI
4890 224 FORMAT(2X,'PSSHARE: ',I2,'-TH SOFT POMERON IN THE BLOCK'/
4891 * 4X,'LIGHT CONE MOMENTA FOR THE POMERON:',2E10.3)
4892 CALL XXSTR(WPI,WMI,WPN,WMN,IC11,IC12,IC22,IC21)
4893
4894
4895
4896
4897 IF(JPP.EQ.1)THEN
4898 IF(PSRAN(B10).LT..5D0)THEN
4899 SW=WPN*WMN
4900 IF(WPN.LT.0.D0.OR.WMN.LT.0.D0.OR.
4901 * SW.LT.(AM(ICZ)+AM(2))**2)THEN
4902
4903 if (debug.ge.1)
4904 * write (*,*)'difr,i,j,WPn,WMn,sw,lq1,lq2',
4905 * i,j,WPn,WMn,sw,lq1,lq2
4906 NREJ=NREJ+1
4907 GOTO 100
4908 ENDIF
4909 typevt=3 !high mass diffraction
4910
4911 IF(LQ2.EQ.0)THEN
4912 CALL XXDTG(WPN,WMN,IZP(I),IZT(J),0)
4913 ELSE
4914 WP1=WPN
4915 WM1=AM(ICZ)**2/WP1
4916 EP3(1)=.5D0*(WP1+WM1)
4917 EP3(2)=.5D0*(WP1-WM1)
4918 EP3(3)=0.D0
4919 EP3(4)=0.D0
4920 CALL XXREG(EP3,IZP(I))
4921 WMN=WMN-WM1
4922 WPN=0.D0
4923 ENDIF
4924 GOTO 30
4925 ELSE
4926
4927
4928 IF(DEBUG.GE.3)WRITE (MONIOU,225)
4929 225 FORMAT(2X,'PSSHARE: TRIPLE POMERON CONRITRIBUTION WITH 3 CUT',
4930 *' POMERONS')
4931 WMM(1)=1.D0/WPI
4932 WMN=WMN-WMM(1)
4933
4934 WMM(2)=WMM(1)*PSRAN(B10)
4935 WMM(1)=WMM(1)-WMM(2)
4936 LQ1=2
4937 DO 26 L=1,2
4938 LQ1=LQ1-1
4939
4940 BPI=(1.D0+DEL)*LQ1+1.D0+AHL(ICZ)
4941 BPI=1.D0/BPI
4942 25 XPW=1.-PSRAN(B10)**BPI
4943 IF(PSRAN(B10).GT.XPW**DEL)GOTO 25
4944 WPP(L)=WPN*XPW
4945 WPN=WPN*(1.D0-XPW)
4946
4947 26 CALL XXSTR(WPP(L),WMM(L),WPN,WMN,0,0,0,0)
4948 SW=WPN*WMN
4949 IF(WPN.LT.0.D0.OR.WMN.LT.0.D0.OR.
4950 * SW.LT.(AM(ICZ)+AM(2))**2)THEN
4951 NREJ=NREJ+1
4952 GOTO 100
4953 ENDIF
4954 ENDIF
4955 ENDIF
4956
4957 27 CONTINUE
4958
4959
4960
4961
4962
4963
4964
4965
4966
4967
4968
4969
4970 IF(LQ1.EQ.0.AND.LQ2.EQ.0)THEN
4971 IF(LVA(I).EQ.0.AND.LVB(J).EQ.0)THEN
4972 CALL XXDDFR(WPN,WMN,IZP(I),IZT(J))
4973 ELSEIF(LVA(I).EQ.0)THEN
4974 CALL XXDPR(WPN,WMN,IZP(I),IZT(J),1)
4975 IF(ILB(J).NE.0)THEN
4976 DO 346 L=1,4
4977 346 EP1(L)=ELB(L,J)
4978 EP(1)=.5D0*WMN
4979 EP(2)=-EP(1)
4980 EP(3)=0.D0
4981 EP(4)=0.D0
4982 IPJ1=ILB(J)
4983 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
4984 CALL PSJDEF(IZT(J),IPJ1,EP,EP1,JFL)
4985 IF(JFL.EQ.0)GOTO 100
4986 ENDIF
4987 ELSEIF(LVB(J).EQ.0)THEN
4988 CALL XXDTG(WPN,WMN,IZP(I),IZT(J),1)
4989 IF(ILA(I).NE.0)THEN
4990 DO 347 L=1,4
4991 347 EP1(L)=ELA(L,I)
4992 EP(1)=.5D0*WPN
4993 EP(2)=EP(1)
4994 EP(3)=0.D0
4995 EP(4)=0.D0
4996 IPJ1=ILA(I)
4997 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
4998 CALL PSJDEF(IZP(I),IPJ1,EP,EP1,JFL)
4999 IF(JFL.EQ.0)GOTO 100
5000 ENDIF
5001 ELSE
5002 IF(ILA(I).NE.0)THEN
5003 DO 348 L=1,4
5004 348 EP1(L)=ELA(L,I)
5005 EP(1)=.5D0*WPN
5006 EP(2)=EP(1)
5007 EP(3)=0.D0
5008 EP(4)=0.D0
5009 IPJ1=ILA(I)
5010 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
5011 CALL PSJDEF(IZP(I),IPJ1,EP,EP1,JFL)
5012 IF(JFL.EQ.0)GOTO 100
5013 ENDIF
5014 IF(ILB(J).NE.0)THEN
5015 DO 349 L=1,4
5016 349 EP1(L)=ELB(L,J)
5017 EP(1)=.5D0*WMN
5018 EP(2)=-EP(1)
5019 EP(3)=0.D0
5020 EP(4)=0.D0
5021 IPJ1=ILB(J)
5022 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
5023 CALL PSJDEF(IZT(J),IPJ1,EP,EP1,JFL)
5024 IF(JFL.EQ.0)GOTO 100
5025 ENDIF
5026 ENDIF
5027
5028 ELSEIF(LQ1.EQ.0)THEN
5029 IF(LVA(I).EQ.0)THEN
5030 CALL XXDPR(WPN,WMN,IZP(I),IZT(J),LQ2)
5031 ELSE
5032 IF(ILA(I).NE.0)THEN
5033 DO 350 L=1,4
5034 350 EP1(L)=ELA(L,I)
5035 EP(1)=.5D0*WPN
5036 EP(2)=EP(1)
5037 EP(3)=0.D0
5038 EP(4)=0.D0
5039 IPJ1=ILA(I)
5040 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
5041 CALL PSJDEF(IZP(I),IPJ1,EP,EP1,JFL)
5042 IF(JFL.EQ.0)GOTO 100
5043 ENDIF
5044 ENDIF
5045
5046 ELSEIF(LQ2.EQ.0)THEN
5047 IF(LVB(J).EQ.0)THEN
5048 CALL XXDTG(WPN,WMN,IZP(I),IZT(J),LQ1)
5049 ELSE
5050 IF(ILB(J).NE.0)THEN
5051 DO 352 L=1,4
5052 352 EP1(L)=ELB(L,J)
5053 EP(1)=.5D0*WMN
5054 EP(2)=-EP(1)
5055 EP(3)=0.D0
5056 EP(4)=0.D0
5057 IPJ1=ILB(J)
5058 IF(IABS(IPJ1).EQ.3)IPJ1=IPJ1*4/3
5059 CALL PSJDEF(IZT(J),IPJ1,EP,EP1,JFL)
5060 IF(JFL.EQ.0)GOTO 100
5061 ENDIF
5062 ENDIF
5063 ENDIF
5064
5065
5066
5067
5068
5069 30 LQA(I)=LQ1
5070 LQB(J)=LQ2
5071 WP(I)=WPN
5072 28 WM(J)=WMN
5073 ENDIF
5074
5075
5076
5077 IF(IA(1).EQ.1.AND.LVA(1).NE.0.AND.ILA(1).EQ.0)THEN
5078 EP(1)=.5D0*WP(1)
5079 EP(2)=EP(1)
5080 EP(3)=0.D0
5081 EP(4)=0.D0
5082 EP1(1)=.5D0*WM(NIAS)
5083 EP1(2)=-EP1(1)
5084 EP1(3)=0.D0
5085 EP1(4)=0.D0
5086 CALL PSJDEF(IZP(1),IZT(NIAS),EP,EP1,JFL)
5087 IF(JFL.EQ.0)GOTO 100
5088 ENDIF
5089
5090 CALL XXJETSIM
5091
5092 IF(DEBUG.GE.3)WRITE (MONIOU,227)
5093 227 FORMAT(2X,'PSSHARE - END')
5094 RETURN
5095 END
5096
5097
5098 SUBROUTINE QGSPSTRANS(EP,EY)
5099
5100
5101
5102 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5103 INTEGER DEBUG
5104 DIMENSION EY(3),EP(4)
5105 COMMON /Q_AREA43/ MONIOU
5106 COMMON /Q_DEBUG/ DEBUG
5107 SAVE
5108
5109 IF(DEBUG.GE.2)WRITE (MONIOU,201)EP,EY
5110 201 FORMAT(2X,'QGSPSTRANS - LORENTZ BOOST FOR 4-VECTOR'/4X,'EP=',
5111 * 2X,4(E10.3,1X)/4X,'BOOST PARAMETERS EY=',3E10.3)
5112
5113 DO 1 I=1,3
5114 IF(EY(4-I).NE.1.D0)THEN
5115 WP=(EP(1)+EP(5-I))/EY(4-I)
5116 WM=(EP(1)-EP(5-I))*EY(4-I)
5117 EP(1)=.5D0*(WP+WM)
5118 EP(5-I)=.5D0*(WP-WM)
5119 ENDIF
5120 1 CONTINUE
5121 IF(DEBUG.GE.3)WRITE (MONIOU,202)EP
5122 202 FORMAT(2X,'QGSPSTRANS: TRANSFORMED 4-VECTOR EP=',
5123 * 2X,4(E10.3,1X))
5124 RETURN
5125 END
5126
5127
5128 SUBROUTINE QGSPSTRANS1(EP,EY)
5129
5130
5131
5132 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5133 INTEGER DEBUG
5134 DIMENSION EY(3),EP(4)
5135 COMMON /Q_AREA43/ MONIOU
5136 COMMON /Q_DEBUG/ DEBUG
5137 SAVE
5138
5139 IF(DEBUG.GE.2)WRITE (MONIOU,201)EP,EY
5140 201 FORMAT(2X,'QGSPSTRANS1 - LORENTZ BOOST FOR 4-VECTOR'/4X,'EP=',
5141 * 2X,4(E10.3,1X)/4X,'BOOST PARAMETERS EY=',3E10.3)
5142
5143 DO 2 I=1,3
5144 IF(EY(I).NE.1.D0)THEN
5145 WP=(EP(1)+EP(I+1))*EY(I)
5146 WM=(EP(1)-EP(I+1))/EY(I)
5147 EP(1)=.5D0*(WP+WM)
5148 EP(I+1)=.5D0*(WP-WM)
5149 ENDIF
5150 2 CONTINUE
5151 IF(DEBUG.GE.3)WRITE (MONIOU,202)EP
5152 202 FORMAT(2X,'QGSPSTRANS1: TRANSFORMED 4-VECTOR EP=',
5153 * 2X,4(E10.3,1X))
5154 RETURN
5155 END
5156
5157
5158 FUNCTION PSUDINT(QLMAX,J)
5159
5160
5161
5162
5163 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5164 INTEGER DEBUG
5165 DIMENSION WK(3)
5166 COMMON /Q_AREA33/ FSUD(10,2)
5167 COMMON /Q_AREA43/ MONIOU
5168 COMMON /Q_DEBUG/ DEBUG
5169 SAVE
5170
5171 IF(DEBUG.GE.2)WRITE (MONIOU,201)J,QLMAX
5172 201 FORMAT(2X,'PSUDINT - SPACELIKE FORM FACTOR INTERPOLATION:'/
5173 * 4X,'PARTON TYPE J=',
5174 * I1,2X,'MOMENTUM LOGARITHM QLMAX=',E10.3)
5175 QL=QLMAX/1.38629d0
5176
5177 IF(QL.LE.0.D0)THEN
5178 PSUDINT=1.D0
5179 ELSE
5180 K=INT(QL)
5181 IF(K.GT.7)K=7
5182 WK(2)=QL-K
5183 WK(3)=WK(2)*(WK(2)-1.D0)*.5D0
5184 WK(1)=1.D0-WK(2)+WK(3)
5185 WK(2)=WK(2)-2.D0*WK(3)
5186
5187 PSUDINT=0.D0
5188 DO 1 K1=1,3
5189 1 PSUDINT=PSUDINT+FSUD(K+K1,J)*WK(K1)
5190 IF(PSUDINT.LE.0.D0)PSUDINT=0.D0
5191 PSUDINT=EXP(-PSUDINT)
5192 ENDIF
5193 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSUDINT
5194 202 FORMAT(2X,'PSUDINT=',E10.3)
5195 RETURN
5196 END
5197
5198
5199 FUNCTION QGSPSUDS(Q,J)
5200
5201
5202
5203
5204 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5205 INTEGER DEBUG
5206 COMMON /Q_AREA6/ PI,BM,AM
5207 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
5208 COMMON /Q_AREA43/ MONIOU
5209 COMMON /Q_DEBUG/ DEBUG
5210 SAVE
5211
5212 IF(DEBUG.GE.2)WRITE (MONIOU,201)J,Q
5213 201 FORMAT(2X,'QGSPSUDS - SPACELIKE FORM FACTOR: PARTON TYPE J=',
5214 * I1,2X,'MOMENTUM Q=',E10.3)
5215 IF(Q.GT.QT0)THEN
5216 QLM=DLOG(Q/ALM)
5217 QGSPSUDS=(QLM*DLOG(QLM/QLOG)-DLOG(Q/QT0))/9.D0
5218
5219 IF(J.EQ.0)THEN
5220 QGSPSUDS=QGSPSUDS*6.D0
5221 ELSE
5222 QGSPSUDS=QGSPSUDS/.375D0
5223 ENDIF
5224 QGSPSUDS=EXP(-QGSPSUDS)
5225
5226 ELSE
5227 QGSPSUDS=1.D0
5228 ENDIF
5229 IF(DEBUG.GE.3)WRITE (MONIOU,202)QGSPSUDS
5230 202 FORMAT(2X,'QGSPSUDS=',E10.3)
5231 RETURN
5232 END
5233
5234
5235 FUNCTION PSUDT(QMAX,J)
5236
5237
5238
5239
5240 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5241 INTEGER DEBUG
5242 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
5243 COMMON/AR13/X1(7),A1(7)
5244 COMMON /Q_AREA43/ MONIOU
5245 COMMON /Q_DEBUG/ DEBUG
5246 SAVE
5247
5248 IF(DEBUG.GE.2)WRITE (MONIOU,201)J,QMAX
5249 201 FORMAT(2X,'PSUDT - TIMELIKE FORM FACTOR: PARTON TYPE J=',
5250 * I1,2X,'MOMENTUM QMAX=',E10.3)
5251 PSUDT=0.D0
5252 QLMAX=DLOG(DLOG(QMAX/16.D0/ALM))
5253 QFL=DLOG(DLOG(QTF/ALM))
5254
5255
5256
5257 DO 1 I=1,7
5258 DO 1 M=1,2
5259 QTL=.5D0*(QLMAX+QFL+(2*M-3)*X1(I)*(QLMAX-QFL))
5260 QT=ALM*EXP(EXP(QTL))
5261 IF(QT.GE.QMAX/16.D0)QT=QMAX/16.0001D0
5262 ZMIN=.5D0-DSQRT((.25D0-DSQRT(QT/QMAX)))
5263 ZMAX=1.D0-ZMIN
5264 IF(J.EQ.0)THEN
5265
5266 AP=(PSAPINT(ZMAX,0,0)-PSAPINT(ZMIN,0,0)+
5267 * PSAPINT(ZMAX,0,1)-PSAPINT(ZMIN,0,1))*.5D0
5268
5269 ELSE
5270 AP=PSAPINT(ZMAX,1,0)-PSAPINT(ZMIN,1,0)
5271 ENDIF
5272 1 PSUDT=PSUDT+A1(I)*AP
5273 PSUDT=PSUDT*(QLMAX-QFL)/9.D0
5274 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSUDT
5275 202 FORMAT(2X,'PSUDT=',E10.3)
5276 RETURN
5277 END
5278
5279
5280 FUNCTION PSV(X,Y,XB,IB)
5281
5282
5283
5284 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5285 INTEGER DEBUG
5286
5287 DIMENSION XB(210,3),FHARD(3)
5288 COMMON /Q_AREA43/ MONIOU
5289 COMMON /Q_DEBUG/ DEBUG
5290 SAVE
5291
5292 IF(DEBUG.GE.2)WRITE (MONIOU,201)X,Y,IB
5293 201 FORMAT(2X,'PSV - EIKONAL FACTOR: NUCLEON COORDINATES X=',
5294 * E10.3,2X,'Y=',E10.3/4X,'NUMBER OF ACTIVE TARGET NUCLEONS IB='
5295 * ,I2)
5296 DV=0.D0
5297
5298 DO 1 M=1,IB
5299 Z=PSDR(X-XB(M,1),Y-XB(M,2))
5300 DV=DV+PSFAZ(Z,FSOFT,FHARD,FSHARD)+FHARD(1)+FHARD(2)+FHARD(3)
5301 1 CONTINUE
5302 PSV=(1.D0-EXP(-DV))**2
5303
5304
5305
5306
5307
5308
5309
5310 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSV
5311 202 FORMAT(2X,'PSV=',E10.3)
5312 RETURN
5313 END
5314
5315
5316 SUBROUTINE PSVDEF(ICH,IC1,ICZ)
5317
5318
5319
5320 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5321 INTEGER DEBUG
5322 COMMON /Q_AREA11/ B10
5323 COMMON /Q_AREA43/ MONIOU
5324 COMMON /Q_DEBUG/ DEBUG
5325 SAVE
5326 EXTERNAL PSRAN
5327
5328 IF(DEBUG.GE.2)WRITE (MONIOU,201)ICH,ICZ
5329 201 FORMAT(2X,'PSVDEF: HADRON TYPE ICH=',I2,' AUXILLIARY TYPE ICZ='
5330 * ,I1)
5331
5332 IS=IABS(ICH)/ICH
5333 IF(ICZ.EQ.1)THEN
5334 IC1=ICH*(1-3*INT(.5+PSRAN(B10)))
5335 ICH=-IC1-ICH
5336 ELSEIF(ICZ.EQ.2)THEN
5337 IF(PSRAN(B10).GT..33333D0.OR.ICH.LT.0)THEN
5338 IC1=ICH-IS
5339 ICH=3*IS
5340 ELSE
5341 IC1=4*IS-ICH
5342 ICH=ICH+4*IS
5343 ENDIF
5344 ELSEIF(ICZ.EQ.3)THEN
5345 IC1=ICH-3*IS
5346 ICH=-4*IS
5347 ELSEIF(ICZ.EQ.4)THEN
5348 IC1=ICH-9*IS
5349 ICH=5*IS
5350 ENDIF
5351 IF(DEBUG.GE.3)WRITE (MONIOU,202)IC1,ICH
5352 202 FORMAT(2X,'PSVDEF-END: QUARK FLAVOR IC1=',I2,
5353 * 'DIQUARK TYPE ICH=',I2)
5354 RETURN
5355 END
5356
5357
5358 FUNCTION PSZSIM(QQ,J)
5359
5360
5361
5362
5363
5364 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5365 INTEGER DEBUG
5366 COMMON /Q_AREA11/ B10
5367 COMMON /Q_AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0
5368 COMMON /Q_AREA43/ MONIOU
5369 COMMON /Q_DEBUG/ DEBUG
5370 SAVE
5371 EXTERNAL PSRAN
5372
5373
5374 IF(DEBUG.GE.2)WRITE (MONIOU,201)QQ,J
5375 201 FORMAT(2X,'PSZSIM - Z-SHARE SIMULATION: QQ=',E10.3,2X,'J=',I1)
5376 ZMIN=.5D0-DSQRT(.25D0-DSQRT(QTF/QQ))
5377 QLF=DLOG(QTF/ALM)
5378
5379 1 CONTINUE
5380 IF(J.EQ.1)THEN
5381 PSZSIM=.5D0*(2.D0*ZMIN)**PSRAN(B10)
5382
5383 GB=PSZSIM*(QGSPSFAP(PSZSIM,0,0)+QGSPSFAP(PSZSIM,0,1))/7.5D0
5384
5385 ELSE
5386 PSZSIM=ZMIN*((1.D0-ZMIN)/ZMIN)**PSRAN(B10)
5387 GB=PSZSIM*QGSPSFAP(PSZSIM,1,0)*.375D0
5388 ENDIF
5389 QT=QQ*(PSZSIM*(1.D0-PSZSIM))**2
5390 GB=GB/DLOG(QT/ALM)*QLF
5391 IF(DEBUG.GE.3)WRITE (MONIOU,203)QT,GB
5392 203 FORMAT(2X,'PSZSIM: QT=',E10.3,2X,'GB=',E10.3)
5393 IF(PSRAN(B10).GT.GB)GOTO 1
5394 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSZSIM
5395 202 FORMAT(2X,'PSZSIM=',E10.3)
5396 RETURN
5397 END
5398
5399
5400 SUBROUTINE IXXDEF(ICH,IC1,IC2,ICZ)
5401
5402
5403
5404 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5405 INTEGER DEBUG
5406 COMMON /Q_AREA11/ B10
5407 COMMON /Q_AREA43/ MONIOU
5408 COMMON /Q_DEBUG/ DEBUG
5409 SAVE
5410 EXTERNAL PSRAN
5411
5412 IF(DEBUG.GE.2)WRITE (MONIOU,201)ICH,ICZ
5413 201 FORMAT(2X,'IXXDEF: HADRON TYPE ICH=',I2,' AUXILLIARY TYPE ICZ='
5414 * ,I1)
5415 IS=IABS(ICH)/ICH
5416 IF(ICZ.EQ.1)THEN
5417 IC1=ICH*(1-3*INT(.5+PSRAN(B10)))
5418 ICH1=ICH*INT(.5D0+PSRAN(B10))
5419 IC2=-IC1*IABS(ICH1)-(ICH+IC1)*IABS(ICH-ICH1)
5420
5421 ELSEIF(ICZ.EQ.2)THEN
5422
5423 IC1=INT(1.3333+PSRAN(B10))
5424
5425 ICH1=(2-IC1)*INT(PSRAN(B10)+.5)+2
5426
5427
5428 IC2=(3-ICH1)*(2-IC1)-2
5429
5430 IF(IABS(ICH).EQ.3)THEN
5431 IC1=3-IC1
5432 IC2=-3-IC2
5433 ICH1=5-ICH1
5434 ENDIF
5435 IF(ICH.LT.0)THEN
5436 IC1=-IC1
5437 IC2=-IC2
5438 ICH1=-ICH1
5439 ENDIF
5440
5441 ELSEIF(ICZ.EQ.3)THEN
5442 IC1=ICH-3*IS
5443 IC2=-IS*INT(1.5+PSRAN(B10))
5444 ICH1=3*IS-IC2
5445 ELSEIF(ICZ.EQ.4)THEN
5446 IC1=ICH-9*IS
5447 IC2=IS*INT(1.5+PSRAN(B10))
5448 ICH1=9*IS-IC2
5449 ELSEIF(ICZ.EQ.5)THEN
5450 IC1=IS*INT(1.5+PSRAN(B10))
5451 IC2=-IC1
5452 ICH1=ICH
5453 ENDIF
5454
5455 ICH=ICH1
5456 IF(DEBUG.GE.3)WRITE (MONIOU,202)IC1,IC2,ICH
5457 202 FORMAT(2X,'IXXDEF-END: PARTON FLAVORS IC1=',I2,' IC2=',I2,
5458 * 'NEW HADRON TYPE ICH=',I2)
5459 RETURN
5460 END
5461
5462
5463 FUNCTION IXXSON(NS,AW,G)
5464
5465
5466
5467
5468
5469 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5470 INTEGER DEBUG
5471 COMMON /Q_AREA43/ MONIOU
5472 COMMON /Q_DEBUG/ DEBUG
5473 SAVE
5474
5475 IF(DEBUG.GE.2)WRITE (MONIOU,201)AW,NS-1,G
5476 201 FORMAT(2X,'IXXSON - POISSON DITR.: AVERAGE AW=',E10.3,
5477 * ' MAXIMAL VALUE NS=',I2,' RANDOM NUMBER G=',E10.3)
5478 W=EXP(-AW)
5479 SUMM=W
5480 DO 1 I=1,NS
5481 J = I
5482 IF(G.LT.SUMM)GOTO 2
5483 W=W*AW/I
5484 1 SUMM=SUMM+W
5485 2 IXXSON=J-1
5486 IF(DEBUG.GE.3)WRITE (MONIOU,202)IXXSON
5487 202 FORMAT(2X,'IXXSON=',I2)
5488 RETURN
5489 END
5490
5491
5492 SUBROUTINE XXAINI(E0N,ICP0,IAP,IAT)
5493
5494
5495 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5496 INTEGER DEBUG
5497
5498 DIMENSION WK(3),WA(3)
5499
5500 COMMON /Q_AREA1/ IA(2),ICZ,ICP
5501 COMMON /Q_AREA2/ S,Y0,WP0,WM0
5502 COMMON /Q_AREA4/ EY0(3)
5503 COMMON /Q_AREA5/ RD(2),CR1(2),CR2(2),CR3(2)
5504 COMMON /Q_AREA6/ PI,BM,AM
5505 COMMON /Q_AREA7/ RP1
5506 COMMON /Q_AREA10/ STMASS,AM0,AMN,AMK,AMC,AMLAMC,AMLAM,AMETA
5507 COMMON /Q_AREA15/ FP(5),RQ(5),CD(5)
5508 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALFP,RR,SH,DELH
5509 COMMON /Q_AREA22/ SJV,FJS(5,3)
5510 COMMON /Q_AREA35/ SJV0(10,5),FJS0(10,5,15)
5511 COMMON /Q_AREA43/ MONIOU
5512
5513 COMMON /Q_AREA44/ GZ(10,5,4),GZP(10,5,4)
5514 COMMON /Q_AREA45/ GDT,GDP !so00
5515
5516 COMMON /Q_DEBUG/ DEBUG
5517 COMMON /Q_QGSNEX1/ XA,XB,BQGS,BMAXQGS,BMAXNEX,BMINNEX !ctp
5518 DIMENSION XA(210,3),XB(210,3)
5519 SAVE
5520
5521 IF(DEBUG.GE.1)WRITE (MONIOU,201)ICP0,IAP,IAT,E0N
5522 201 FORMAT(2X,'XXAINI - MINIINITIALIZATION: PARTICLE TYPE ICP0=',
5523 * I1,2X,'PROJECTILE MASS NUMBER IAP=',I2/4X,
5524 * 'TARGET MASS NUMBER IAT=',I2,' INTERACTION ENERGY E0N=',E10.3)
5525 ICP=ICP0
5526 IA(1)=IAP
5527 IA(2)=IAT
5528
5529
5530 IF(IABS(ICP).LT.6)THEN
5531 ICZ=IABS(ICP)/2+1
5532 ELSE
5533 ICZ=(IABS(ICP)+1)/2
5534 ENDIF
5535
5536
5537
5538 S=2.D0*E0N*AMN
5539 WP0=DSQRT(S)
5540 WM0=WP0
5541
5542 Y0=DLOG(S)
5543
5544 RS=RQ(ICZ)+ALFP*Y0
5545
5546 RS0=RQ(ICZ)
5547
5548 FS=FP(ICZ)*EXP(Y0*DEL)/RS*CD(ICZ)
5549
5550 RP1=RS*4.D0*.0391D0/AM**2
5551
5552 EY0(2)=1.D0
5553 EY0(3)=1.D0
5554 EY0(1)=DSQRT(AMN/E0N/2.D0)
5555
5556
5557
5558 DO 1 I=1,2
5559
5560 RD(I)=0.7D0*FLOAT(IA(I))**.446/AM
5561 CR1(I)=1.D0+3.D0/RD(I)+6.D0/RD(I)**2+6.D0/RD(I)**3
5562 CR2(I)=3.D0/RD(I)
5563 CR3(I)=3.D0/RD(I)+6.D0/RD(I)**2
5564 IF(IA(I).LT.10.AND.IA(I).NE.1)THEN
5565
5566 RD(I)=.9D0*FLOAT(IA(I))**.3333/AM
5567 IF(IA(I).EQ.2)RD(I)=3.16D0
5568
5569 RD(I)=RD(I)*DSQRT(2.D0*IA(I)/(IA(I)-1.))
5570 ENDIF
5571 1 CONTINUE
5572
5573 GDT=0.D0
5574
5575
5576
5577 IF(IA(1).NE.1)THEN
5578
5579
5580
5581 BM=RD(1)+RD(2)+5.D0
5582 ELSE
5583
5584
5585 BM=RD(2)+5.D0
5586 ENDIF
5587
5588 BMAXQGS=BM*AM !ctp
5589
5590 YE=DLOG10(E0N)
5591 IF(YE.LT.1.D0)YE=1.D0
5592 JE=INT(YE)
5593 IF(JE.GT.8)JE=8
5594
5595
5596 WK(2)=YE-JE
5597 WK(3)=WK(2)*(WK(2)-1.D0)*.5D0
5598 WK(1)=1.D0-WK(2)+WK(3)
5599 WK(2)=WK(2)-2.D0*WK(3)
5600
5601 SJV=SJV0(JE,ICZ)*WK(1)+SJV0(JE+1,ICZ)*WK(2)+SJV0(JE+2,ICZ)*WK(3)
5602
5603 DO 2 I=1,5
5604 DO 2 M=1,3
5605 M1=M+3*(ICZ-1)
5606 2 FJS(I,M)=FJS0(JE,I,M1)*WK(1)+FJS0(JE+1,I,M1)*WK(2)+
5607 *FJS0(JE+2,I,M1)*WK(3)
5608
5609 GDT=0.D0
5610 GDP=0.D0 !so00
5611 IF(IA(1).EQ.1)THEN
5612 YA=IA(2)
5613 YA=DLOG(YA)/1.38629D0+1.D0
5614 JA=MIN(INT(YA),2)
5615 WA(2)=YA-JA
5616 WA(3)=WA(2)*(WA(2)-1.D0)*.5D0
5617 WA(1)=1.D0-WA(2)+WA(3)
5618 WA(2)=WA(2)-2.D0*WA(3)
5619 DO 3 I=1,3
5620 DO 3 M=1,3
5621 GDP=GDP+GZP(JE+I-1,ICZ,JA+M-1)*WK(I)*WA(M) !so00
5622 3 GDT=GDT+GZ(JE+I-1,ICZ,JA+M-1)*WK(I)*WA(M)
5623 ENDIF
5624
5625
5626
5627 IF(DEBUG.GE.3)WRITE (MONIOU,202)
5628 202 FORMAT(2X,'XXAINI - END')
5629 RETURN
5630 END
5631
5632
5633 SUBROUTINE XXASET
5634
5635
5636 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5637 INTEGER DEBUG
5638 CHARACTER *2 TYQ
5639 COMMON /Q_AREA3/ RMIN,EMAX,EEV
5640 COMMON /Q_AREA6/ PI,BM,AM
5641 COMMON /Q_AREA8/ WWM,BE(4),DC(5),DETA,ALMPT
5642 COMMON /Q_AREA10/ STMASS,AM0,AMN,AMK,AMC,AMLAMC,AMLAM,AMETA
5643 COMMON /Q_AREA11/ B10
5644 COMMON /Q_AREA20/ WPPP
5645 COMMON /Q_AREA21/ DMMIN(5)
5646 COMMON /Q_AREA28/ ARR(4)
5647 COMMON /Q_AREA40/ JDIFR
5648 COMMON /Q_AREA42/ TYQ(15)
5649 COMMON /Q_AREA43/ MONIOU
5650 COMMON /Q_DEBUG/ DEBUG
5651 SAVE
5652
5653 IF(DEBUG.GE.1)WRITE (MONIOU,201)
5654 201 FORMAT(2X,'XXASET - HADRONIZATION PARAMETERS SETTING')
5655
5656 ARR(1)=0.5D0
5657 ARR(2)=-.5D0
5658 ARR(3)=0.D0
5659 ARR(4)=-2.D0
5660
5661
5662 WPPP=0.4d0
5663
5664 JDIFR=1
5665
5666
5667
5668
5669
5670
5671 DC(1)=.06D0
5672 DC(2)=.10D0
5673
5674 DC(3)=.000D0
5675 DC(4)=.36D0
5676
5677 DC(5)=.0D0
5678
5679 DETA=.11111D0
5680
5681
5682 WWM=.53D0
5683
5684
5685 BE(1)=.22D0
5686 BE(2)=.35D0
5687 BE(3)=.29D0
5688 BE(4)=.40D0
5689
5690
5691
5692 ALMPT=1.7D0
5693
5694
5695
5696
5697
5698
5699 RMIN=3.35D0
5700 EMAX=.11D0
5701 EEV=.25D0
5702
5703
5704
5705
5706 DMMIN(1)=.76D0
5707 DMMIN(2)=1.24D0
5708 DMMIN(3)=.89D0
5709 DMMIN(4)=2.01D0
5710 DMMIN(5)=2.45D0
5711
5712 AMN=.939D0
5713 AMK=.496D0
5714 AM0=.14D0
5715 AMC=1.868D0
5716 AMLAM=1.116D0
5717 AMLAMC=2.27D0
5718 AMETA=.548D0
5719
5720
5721
5722
5723
5724 B10=.43876194D0
5725 PI=3.1416D0
5726 AM=.523D0
5727
5728
5729 STMASS=4.D0*AM0**2
5730
5731 RMIN=RMIN/AM**2
5732
5733 TYQ(1)='DD'
5734 TYQ(2)='UU'
5735 TYQ(3)='C '
5736 TYQ(4)='S '
5737 TYQ(5)='UD'
5738 TYQ(6)='D '
5739 TYQ(7)='U '
5740 TYQ(8)='G '
5741 TYQ(9)='u '
5742 TYQ(10)='d '
5743 TYQ(11)='ud'
5744 TYQ(12)='s '
5745 TYQ(13)='c '
5746 TYQ(14)='uu'
5747 TYQ(15)='dd'
5748 IF(DEBUG.GE.3)WRITE (MONIOU,202)
5749 202 FORMAT(2X,'XXASET - END')
5750 RETURN
5751 END
5752
5753
5754 SUBROUTINE XXDDFR(WP0,WM0,ICP,ICT)
5755
5756
5757 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5758 INTEGER DEBUG
5759 DIMENSION EP3(4),EY(3)!,EP1(4),EP2(4)
5760 COMMON /Q_AREA1/ IA(2),ICZ,ICP0
5761 COMMON /Q_AREA2/ S,Y0,WP00,WM00
5762 COMMON /Q_AREA8/ WWM,BE(4),DC(5),DETA,ALMPT
5763 COMMON /Q_AREA10/ STMASS,AM(7)
5764 COMMON /Q_AREA11/ B10
5765 COMMON /Q_AREA21/ DMMIN(5)
5766 COMMON /Q_AREA43/ MONIOU
5767 COMMON /Q_DEBUG/ DEBUG
5768 SAVE
5769 EXTERNAL PSRAN
5770
5771 IF(DEBUG.GE.2)WRITE (MONIOU,201)ICP,ICT,WP0,WM0
5772 201 FORMAT(2X,'XXDDFR - LEADING CLUSTERS HADRONIZATION:'
5773 * /4X,'CLUSTER TYPES ICP=',I2,2X,
5774 * 'ICT=',I2/4X,'AVAILABLE LIGHT CONE MOMENTA: WP0=',E10.3,
5775 * ' WM0=',E10.3)
5776 DO 100 I=1,3
5777 100 EY(I)=1.D0
5778
5779 SD0=WP0*WM0
5780 IF(SD0.LT.0.D0)SD0=0.D0
5781 DDMIN1=DMMIN(ICZ)
5782 DDMIN2=DMMIN(2)
5783 DDMAX1=MIN(5.D0,DSQRT(SD0)-DDMIN2)
5784
5785 IF(DDMAX1.LT.DDMIN1)THEN
5786
5787
5788 IF(DSQRT(SD0).LT.AM(ICZ)+AM(2))THEN
5789 IF(WP0.GT.0.D0.AND.(AM(ICZ)+AM(2))**2/WP0.LT..5D0*WM00)THEN
5790 SD0=(AM(ICZ)+AM(2))**2
5791 WM0=SD0/WP0
5792 ELSE
5793 IF(DEBUG.GE.3)WRITE (MONIOU,202)
5794 RETURN
5795 ENDIF
5796 ENDIF
5797
5798 EP3(3)=0.D0
5799 EP3(4)=0.D0
5800 XW=XXTWDEC(SD0,AM(ICZ)**2,AM(2)**2)
5801 WP1=XW*WP0
5802 WM1=AM(ICZ)**2/WP1
5803 EP3(1)=.5D0*(WP1+WM1)
5804 EP3(2)=.5D0*(WP1-WM1)
5805 CALL XXREG(EP3,ICP)
5806 WM2=WM0-WM1
5807 WP2=AM(2)**2/WM2
5808 EP3(1)=.5D0*(WP2+WM2)
5809 EP3(2)=.5D0*(WP2-WM2)
5810 CALL XXREG(EP3,ICT)
5811 WP0=0.D0
5812 WM0=0.D0
5813 IF(DEBUG.GE.3)WRITE (MONIOU,202)
5814 RETURN
5815 ENDIF
5816
5817 DMASS1=(DDMIN1/(1.D0-PSRAN(B10)*(1.D0-DDMIN1/DDMAX1)))**2
5818 DDMAX2=MIN(5.D0,DSQRT(SD0)-DSQRT(DMASS1))
5819 DMASS2=(DDMIN2/(1.D0-PSRAN(B10)*(1.D0-DDMIN2/DDMAX2)))**2
5820
5821 WPD1=WP0*XXTWDEC(SD0,DMASS1,DMASS2)
5822 WMD1=DMASS1/WPD1
5823 WMD2=WM0-WMD1
5824 WPD2=DMASS2/WMD2
5825
5826 IF(ICP.NE.0)IS=IABS(ICP)/ICP
5827 IF(ICZ.EQ.5)THEN
5828 ICH1=ICP
5829 ICH2=0
5830 AMH1=AM(5)**2
5831 AMH2=AM(1)**2
5832
5833 PTMAX=QGSPSLAM(DMASS1,AMH1,AMH2)
5834 IF(PTMAX.LT.0.)PTMAX=0.
5835 IF(PTMAX.LT.BE(4)**2)THEN
5836 1 PTI=PTMAX*PSRAN(B10)
5837 IF(PSRAN(B10).GT.EXP(-DSQRT(PTI)/BE(4)))GOTO 1
5838 ELSE
5839 2 PTI=(BE(4)*DLOG(PSRAN(B10)*PSRAN(B10)))**2
5840 IF(PTI.GT.PTMAX)GOTO 2
5841 ENDIF
5842 AMT1=AMH1+PTI
5843 AMT2=AMH2+PTI
5844 Z=XXTWDEC(DMASS1,AMT1,AMT2)
5845 WP1=WPD1*Z
5846 WM1=AMT1/WP1
5847 EP3(1)=.5D0*(WP1+WM1)
5848 EP3(2)=.5D0*(WP1-WM1)
5849 PT=DSQRT(PTI)
5850 CALL QGSPSCS(C,S)
5851 EP3(3)=PT*C
5852 EP3(4)=PT*S
5853 CALL XXREG(EP3,ICH1)
5854
5855 WP1=WPD1*(1.D0-Z)
5856 WM1=AMT2/WP1
5857 EP3(1)=.5D0*(WP1+WM1)
5858 EP3(2)=.5D0*(WP1-WM1)
5859 EP3(3)=-PT*C
5860 EP3(4)=-PT*S
5861 CALL XXREG(EP3,ICH2)
5862 GOTO 3
5863 ENDIF
5864
5865 IF(ICZ.EQ.1)THEN
5866 IF(ICP.NE.0)THEN
5867 IC1=ICP*(1-3*INT(.5D0+PSRAN(B10)))
5868 IC2=-ICP-IC1
5869 ELSE
5870 IC1=INT(1.5D0+PSRAN(B10))*(2*INT(.5D0+PSRAN(B10))-1)
5871 IC2=-IC1
5872 ENDIF
5873 ELSEIF(ICZ.EQ.2)THEN
5874 IF(PSRAN(B10).GT..33333D0)THEN
5875 IC1=3*IS
5876 IC2=ICP-IS
5877 ELSE
5878 IC1=ICP+4*IS
5879 IC2=4*IS-ICP
5880 ENDIF
5881 ELSEIF(ICZ.EQ.3)THEN
5882 IC1=-4*IS
5883 IC2=ICP-3*IS
5884 ELSEIF(ICZ.EQ.4)THEN
5885 IC1=5*IS
5886 IC2=ICP-9*IS
5887 ENDIF
5888 CALL XXGENER(WPD1,WMD1,EY,0.D0,1.D0,0.D0,1.D0,IC1,IC2)
5889
5890 3 CONTINUE
5891 IS=IABS(ICT)/ICT
5892 IF(PSRAN(B10).GT..33333D0)THEN
5893 IC1=3*IS
5894 IC2=ICT-IS
5895 ELSE
5896 IC1=ICT+4*IS
5897 IC2=4*IS-ICT
5898 ENDIF
5899 CALL XXGENER(WPD2,WMD2,EY,0.D0,1.D0,0.D0,1.D0,IC2,IC1)
5900 IF(DEBUG.GE.3)WRITE (MONIOU,202)
5901 202 FORMAT(2X,'XXDDFR - END')
5902 RETURN
5903 END
5904
5905
5906 SUBROUTINE XXDEC2(EP,EP1,EP2,WW,A,B)
5907
5908
5909 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5910 INTEGER DEBUG
5911 dimension ep(4),ep1(4),ep2(4),EY(3)
5912 COMMON /Q_AREA43/ MONIOU
5913 COMMON /Q_DEBUG/ DEBUG
5914 COMMON /Q_AREA11/ B10
5915 SAVE
5916 EXTERNAL PSRAN
5917
5918 IF(DEBUG.GE.2)WRITE (MONIOU,201)
5919 201 FORMAT(2X,'XXDEC2 - TWO PARTICLE DECAY')
5920
5921 PL=QGSPSLAM(WW,A,B)
5922 EP1(1)=DSQRT(PL+A)
5923 EP2(1)=DSQRT(PL+B)
5924 PL=DSQRT(PL)
5925 COSZ=2.D0*PSRAN(B10)-1.D0
5926 PT=PL*DSQRT(1.D0-COSZ**2)
5927 EP1(2)=PL*COSZ
5928 CALL QGSPSCS(C,S)
5929 EP1(3)=PT*C
5930 EP1(4)=PT*S
5931 do 1 I=2,4
5932 1 EP2(I)=-EP1(I)
5933 CALL QGSPSDEFTR(WW,EP,EY)
5934 CALL QGSPSTRANS(EP1,EY)
5935 CALL QGSPSTRANS(EP2,EY)
5936 IF(DEBUG.GE.3)WRITE (MONIOU,202)
5937 202 FORMAT(2X,'XXDEC2 - END')
5938 RETURN
5939 END
5940
5941
5942 SUBROUTINE XXDEC3(EP,EP1,EP2,EP3,SWW,AM1,AM2,AM3)
5943
5944
5945 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5946 INTEGER DEBUG
5947 DIMENSION EP(4),EP1(4),EP2(4),EP3(4),EPT(4),EY(3)
5948 COMMON/AREA11/B10
5949 COMMON /Q_AREA43/ MONIOU
5950 COMMON /Q_DEBUG/ DEBUG
5951 SAVE
5952 EXTERNAL PSRAN
5953
5954 IF(DEBUG.GE.2)WRITE (MONIOU,201)
5955 201 FORMAT(2X,'XXDEC3 - THREE PARTICLE DECAY')
5956 AM12=AM1**2
5957 AM23=(AM2+AM3)**2
5958 AM32=(AM2-AM3)**2
5959 S23MAX=(SWW-AM1)**2
5960 EMAX=.25D0*(SWW+(AM12-AM23)/SWW)**2
5961 GB0=DSQRT((EMAX-AM12)/EMAX*(1.D0-AM23/S23MAX)
5962 * *(1.D0-AM32/S23MAX))
5963 1 P1=PSRAN(B10)*(EMAX-AM12)
5964 E1=DSQRT(P1+AM12)
5965 S23=SWW**2+AM12-2.D0*E1*SWW
5966 GB=DSQRT(P1*(1.D0-AM23/S23)*(1.D0-AM32/S23))/E1/GB0
5967 IF(PSRAN(B10).GT.GB)GOTO 1
5968
5969 P1=DSQRT(P1)
5970 EP1(1)=E1
5971 COSZ=2.D0*PSRAN(B10)-1.D0
5972 PT=P1*DSQRT(1.D0-COSZ**2)
5973 EP1(2)=P1*COSZ
5974 CALL QGSPSCS(C,S)
5975 EP1(3)=PT*C
5976 EP1(4)=PT*S
5977 do 2 I=2,4
5978 2 EPT(I)=-EP1(I)
5979 EPT(1)=SWW-EP1(1)
5980 CALL QGSPSDEFTR(SWW**2,EP,EY)
5981 CALL QGSPSTRANS(EP1,EY)
5982 CALL QGSPSTRANS(EPT,EY)
5983
5984 CALL XXDEC2(EPT,EP2,EP3,S23,AM2**2,AM3**2)
5985 IF(DEBUG.GE.3)WRITE (MONIOU,202)
5986 202 FORMAT(2X,'XXDEC3 - END')
5987 RETURN
5988 END
5989
5990
5991 SUBROUTINE XXDPR(WP0,WM0,ICP,ICT,LQ2)
5992
5993
5994
5995 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
5996 INTEGER DEBUG
5997 DIMENSION EP3(4),EY(3)!,EP1(4),EP2(4)
5998 COMMON /Q_AREA1/ IA(2),ICZ,ICP0
5999 COMMON /Q_AREA2/ S,Y0,WP00,WM00
6000 COMMON /Q_AREA8/ WWM,BE(4),DC(5),DETA,ALMPT
6001 COMMON /Q_AREA10/ STMASS,AM(7)
6002 COMMON /Q_AREA11/ B10
6003 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALFP,RR,SH,DELH
6004 COMMON /Q_AREA21/ DMMIN(5)
6005 COMMON /Q_AREA43/ MONIOU
6006 COMMON /Q_DEBUG/ DEBUG
6007 SAVE
6008 EXTERNAL PSRAN
6009
6010
6011 IF(DEBUG.GE.2)WRITE (MONIOU,201)ICP,ICT,WP0,WM0
6012 201 FORMAT(2X,'XXDPR - LEADING (PROJECTILE) CLUSTER HADRONIZATION:'
6013 * /4X,'CLUSTER TYPE ICP=',I2,2X,'TARGET TYPE ',
6014 * 'ICT=',I2/4X,'AVAILABLE LIGHT CONE MOMENTA: WP0=',E10.3,
6015 * ' WM0=',E10.3)
6016 DO 100 I=1,3
6017 100 EY(I)=1.D0
6018
6019 SD0=WP0*WM0
6020 IF(SD0.LT.0.D0)SD0=0.D0
6021 DDMAX=MIN(5.D0,DSQRT(SD0)-AM(2))
6022 DDMIN=DMMIN(ICZ)
6023
6024 IF(DDMAX.LT.DDMIN)THEN
6025
6026
6027 EP3(3)=0.D0
6028 EP3(4)=0.D0
6029
6030 IF(LQ2.NE.0)THEN
6031 WPI=WP0
6032 IF(AM(ICZ)**2.GT.WPI*WM0)THEN
6033 IF(WPI.GT.0.D0.AND.AM(ICZ)**2/WPI.LT..5D0*WM00)THEN
6034 WMI=AM(ICZ)**2/WPI
6035 WM0=WMI
6036 ELSE
6037 RETURN
6038 ENDIF
6039
6040 ELSE
6041 WMI=AM(ICZ)**2/WPI
6042
6043 ENDIF
6044 WM0=WM0-WMI
6045 WP0=0.D0
6046 EP3(1)=.5D0*(WPI+WMI)
6047 EP3(2)=.5D0*(WPI-WMI)
6048 CALL XXREG(EP3,ICP)
6049 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6050 RETURN
6051 ELSE
6052
6053 IF(DSQRT(SD0).LT.AM(ICZ)+AM(2))THEN
6054 IF(WP0.GT.0.D0.AND.(AM(ICZ)+AM(2))**2/WP0.LT..5D0*WM00)
6055 * THEN
6056 SD0=(AM(ICZ)+AM(2))**2
6057 WM0=SD0/WP0
6058 ELSE
6059 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6060 RETURN
6061 ENDIF
6062 ENDIF
6063 XW=XXTWDEC(SD0,AM(ICZ)**2,AM(2)**2)
6064 WP1=XW*WP0
6065 WM1=AM(ICZ)**2/WP1
6066 EP3(1)=.5D0*(WP1+WM1)
6067 EP3(2)=.5D0*(WP1-WM1)
6068 CALL XXREG(EP3,ICP)
6069 WM2=WM0-WM1
6070 WP2=AM(2)**2/WM2
6071 EP3(1)=.5D0*(WP2+WM2)
6072 EP3(2)=.5D0*(WP2-WM2)
6073 CALL XXREG(EP3,ICT)
6074 WP0=0.D0
6075 WM0=0.D0
6076 ENDIF
6077 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6078 RETURN
6079 ENDIF
6080
6081 IF(ICP.NE.0)IS=IABS(ICP)/ICP
6082
6083 DMASS=DDMIN**2/(1.D0-PSRAN(B10)*(1.D0-(DDMIN/DDMAX)))**2
6084
6085 IF(LQ2.NE.0)THEN
6086 WPD=WP0
6087 WMD=DMASS/WPD
6088 WM0=WM0-WMD
6089 WP0=0.D0
6090 ELSE
6091 IF(ICZ.EQ.5)THEN
6092 WPD=WP0*XXTWDEC(SD0,DMASS,AM(2)**2)
6093 WMD=DMASS/WPD
6094 WM2=WM0-WMD
6095 WP2=AM(2)**2/WM2
6096 EP3(1)=.5D0*(WP2+WM2)
6097 EP3(2)=.5D0*(WP2-WM2)
6098 EP3(3)=0.D0
6099 EP3(4)=0.D0
6100 CALL XXREG(EP3,ICT)
6101 ELSE
6102 PTMAX=QGSPSLAM(SD0,DMASS,AM(2)**2)
6103 IF(PTMAX.LT.0.)PTMAX=0.
6104 PTI=-1.D0/RS*DLOG(1.D0-PSRAN(B10)*(1.D0-EXP(-RS*PTMAX)))
6105
6106 AMT1=DMASS+PTI
6107 AMT2=AM(2)**2+PTI
6108 WPD=WP0*XXTWDEC(SD0,AMT1,AMT2)
6109 WMD=AMT1/WPD
6110 WM2=WM0-WMD
6111 WP2=AMT2/WM2
6112 PT=DSQRT(PTI)
6113 CALL QGSPSCS(CCOS,SSIN)
6114 EP3(3)=PT*CCOS
6115 EP3(4)=PT*SSIN
6116 EP3(1)=.5D0*(WP2+WM2)
6117 EP3(2)=.5D0*(WP2-WM2)
6118 CALL XXREG(EP3,ICT)
6119 EP3(3)=-EP3(3)
6120 EP3(4)=-EP3(4)
6121 EP3(1)=.5D0*(WPD+WMD)
6122 EP3(2)=.5D0*(WPD-WMD)
6123 CALL QGSPSDEFTR(DMASS,EP3,EY)
6124 WPD=DSQRT(DMASS)
6125 WMD=WPD
6126 ENDIF
6127 WP0=0.D0
6128 WM0=0.D0
6129 ENDIF
6130
6131 IF(ICZ.EQ.5)THEN
6132 ICH1=ICP
6133 ICH2=0
6134 AMH1=AM(5)**2
6135 AMH2=AM(1)**2
6136
6137 PTMAX=QGSPSLAM(DMASS,AMH1,AMH2)
6138 IF(PTMAX.LT.0.)PTMAX=0.
6139 IF(PTMAX.LT.BE(4)**2)THEN
6140 1 PTI=PTMAX*PSRAN(B10)
6141 IF(PSRAN(B10).GT.EXP(-DSQRT(PTI)/BE(4)))GOTO 1
6142 ELSE
6143 2 PTI=(BE(4)*DLOG(PSRAN(B10)*PSRAN(B10)))**2
6144 IF(PTI.GT.PTMAX)GOTO 2
6145 ENDIF
6146 AMT1=AMH1+PTI
6147 AMT2=AMH2+PTI
6148 Z=XXTWDEC(DMASS,AMT1,AMT2)
6149 WP1=WPD*Z
6150 WM1=AMT1/WP1
6151 EP3(1)=.5D0*(WP1+WM1)
6152 EP3(2)=.5D0*(WP1-WM1)
6153 PT=DSQRT(PTI)
6154 CALL QGSPSCS(C,S)
6155 EP3(3)=PT*C
6156 EP3(4)=PT*S
6157 CALL XXREG(EP3,ICH1)
6158
6159 WP1=WPD*(1.D0-Z)
6160 WM1=AMT2/WP1
6161 EP3(1)=.5D0*(WP1+WM1)
6162 EP3(2)=.5D0*(WP1-WM1)
6163 EP3(3)=-PT*C
6164 EP3(4)=-PT*S
6165 CALL XXREG(EP3,ICH2)
6166 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6167 RETURN
6168 ENDIF
6169
6170 IF(ICZ.EQ.1)THEN
6171 IF(ICP.NE.0)THEN
6172 IC1=ICP*(1-3*INT(.5D0+PSRAN(B10)))
6173 IC2=-ICP-IC1
6174 ELSE
6175 IC1=INT(1.5D0+PSRAN(B10))*(2*INT(.5D0+PSRAN(B10))-1)
6176 IC2=-IC1
6177 ENDIF
6178 ELSEIF(ICZ.EQ.2)THEN
6179 IF(PSRAN(B10).GT..33333D0)THEN
6180 IC1=3*IS
6181 IC2=ICP-IS
6182 ELSE
6183 IC1=ICP+4*IS
6184 IC2=4*IS-ICP
6185 ENDIF
6186 ELSEIF(ICZ.EQ.3)THEN
6187 IC1=-4*IS
6188 IC2=ICP-3*IS
6189 ELSEIF(ICZ.EQ.4)THEN
6190 IC1=5*IS
6191 IC2=ICP-9*IS
6192 ENDIF
6193 CALL XXGENER(WPD,WMD,EY,0.D0,1.D0,0.D0,1.D0,
6194 * IC1,IC2)
6195 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6196 202 FORMAT(2X,'XXDPR - END')
6197 RETURN
6198 END
6199
6200
6201 SUBROUTINE XXDTG(WP0,WM0,ICP,ICT,LQ1)
6202
6203
6204
6205 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6206 INTEGER DEBUG
6207 DIMENSION EP3(4),EY(3)
6208 COMMON /Q_AREA1/ IA(2),ICZ,ICP0
6209 COMMON /Q_AREA2/ S,Y0,WP00,WM00
6210 COMMON /Q_AREA10/ STMASS,AM(7)
6211 COMMON /Q_AREA11/ B10
6212 COMMON /Q_AREA17/ DEL,RS,RS0,FS,ALFP,RR,SH,DELH
6213 COMMON /Q_AREA21/ DMMIN(5)
6214 COMMON /Q_AREA43/ MONIOU
6215 COMMON /Q_DEBUG/ DEBUG
6216 SAVE
6217 EXTERNAL PSRAN
6218
6219
6220 IF(DEBUG.GE.2)WRITE (MONIOU,201)ICT,ICP,WP0,WM0
6221 201 FORMAT(2X,'XXDTG - LEADING (TARGET) CLUSTER HADRONIZATION:'
6222 * /4X,'CLUSTER TYPE ICT=',I2,2X,'PROJECTILE TYPE ',
6223 * 'ICP=',I2/4X,'AVAILABLE LIGHT CONE MOMENTA: WP0=',E10.3,
6224 * ' WM0=',E10.3)
6225 DO 100 I=1,3
6226 100 EY(I)=1.D0
6227
6228 SD0=WP0*WM0
6229 IF(SD0.LT.0.D0)SD0=0.D0
6230 DDMIN=DMMIN(2)
6231 DDMAX=MIN(5.D0,DSQRT(SD0)-AM(ICZ))
6232
6233 IF(DDMAX.LT.DDMIN)THEN
6234
6235
6236 EP3(3)=0.D0
6237 EP3(4)=0.D0
6238
6239 IF(LQ1.NE.0)THEN
6240 WMI=WM0
6241 IF( WP0.LE.0.D0.OR.AM(2)**2.GT.WMI*WP0)RETURN
6242 WPI=AM(2)**2/WMI
6243 WP0=WP0-WPI
6244 WM0=0.D0
6245 EP3(1)=.5D0*(WPI+WMI)
6246 EP3(2)=.5D0*(WPI-WMI)
6247 CALL XXREG(EP3,ICT)
6248 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6249 RETURN
6250 ELSE
6251
6252 IF(DSQRT(SD0).LT.AM(ICZ)+AM(2))THEN
6253 IF(WP0.GT.0.D0.AND.(AM(ICZ)+AM(2))**2/WP0.LT..5D0*WM00)
6254 * THEN
6255 SD0=(AM(ICZ)+AM(2))**2
6256 WM0=SD0/WP0
6257 ELSE
6258 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6259 RETURN
6260 ENDIF
6261 ENDIF
6262 XW=XXTWDEC(SD0,AM(ICZ)**2,AM(2)**2)
6263 WP1=XW*WP0
6264 WM1=AM(ICZ)**2/WP1
6265 EP3(1)=.5D0*(WP1+WM1)
6266 EP3(2)=.5D0*(WP1-WM1)
6267 CALL XXREG(EP3,ICP)
6268 WM2=WM0-WM1
6269 WP2=AM(2)**2/WM2
6270 EP3(1)=.5D0*(WP2+WM2)
6271 EP3(2)=.5D0*(WP2-WM2)
6272 CALL XXREG(EP3,ICT)
6273 WP0=0.D0
6274 WM0=0.D0
6275 ENDIF
6276 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6277 RETURN
6278 ENDIF
6279
6280 DMASS=(DDMIN/(1.D0-PSRAN(B10)*(1.D0-DDMIN/DDMAX)))**2
6281 IF(LQ1.NE.0)THEN
6282 WMD=WM0
6283 WPD=DMASS/WMD
6284 WP0=WP0-WPD
6285 WM0=0.D0
6286 ELSE
6287 PTMAX=QGSPSLAM(SD0,DMASS,AM(ICZ)**2)
6288 IF(PTMAX.LT.0.)PTMAX=0.
6289 PTI=-1.D0/RS*DLOG(1.D0-PSRAN(B10)*(1.D0-EXP(-RS*PTMAX)))
6290
6291 AMT1=DMASS+PTI
6292 AMT2=AM(ICZ)**2+PTI
6293 WMD=WM0*XXTWDEC(SD0,AMT1,AMT2)
6294 WPD=AMT1/WMD
6295 WP2=WP0-WPD
6296 WM2=AMT2/WP2
6297 PT=DSQRT(PTI)
6298 CALL QGSPSCS(CCOS,SSIN)
6299 EP3(3)=PT*CCOS
6300 EP3(4)=PT*SSIN
6301 EP3(1)=.5D0*(WP2+WM2)
6302 EP3(2)=.5D0*(WP2-WM2)
6303 CALL XXREG(EP3,ICP)
6304 EP3(3)=-EP3(3)
6305 EP3(4)=-EP3(4)
6306 EP3(1)=.5D0*(WPD+WMD)
6307 EP3(2)=.5D0*(WPD-WMD)
6308 CALL QGSPSDEFTR(DMASS,EP3,EY)
6309 WPD=DSQRT(DMASS)
6310 WMD=WPD
6311 WP0=0.D0
6312 WM0=0.D0
6313 ENDIF
6314
6315 IS=IABS(ICT)/ICT
6316 IF(PSRAN(B10).GT..33333D0)THEN
6317 IC1=3*IS
6318 IC2=ICT-IS
6319 ELSE
6320 IC1=ICT+4*IS
6321 IC2=4*IS-ICT
6322 ENDIF
6323 CALL XXGENER(WPD,WMD,EY,
6324 * 0.D0,1.D0,0.D0,1.D0,IC2,IC1)
6325 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6326 202 FORMAT(2X,'XXDTG - END')
6327 RETURN
6328 END
6329
6330
6331 SUBROUTINE XXFAU(B,GZ)
6332
6333
6334 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6335 INTEGER DEBUG
6336 DIMENSION GZ(3),GZ0(2)
6337 COMMON /Q_AREA1/ IA(2),ICZ,ICP
6338 COMMON /Q_AREA16/ CC(5)
6339 COMMON /Q_AR1/ ANORM
6340 COMMON /Q_AREA43/ MONIOU
6341 COMMON /Q_DEBUG/ DEBUG
6342 SAVE
6343
6344 IF(DEBUG.GE.2)WRITE (MONIOU,201)
6345 201 FORMAT(2X,'XXFAU - INTEGRANDS FOR HADRON-HADRON AND '
6346 * ,'HADRON-NUCLEUS CROSS-SECTIONS CALCULATION')
6347
6348 CALL XXFZ(B,GZ0)
6349 DO 1 L=1,2
6350 1 GZ0(L)=GZ0(L)*CC(2)*ANORM*.5D0
6351
6352 AB=FLOAT(IA(2))
6353
6354 GZ1=(1.D0-GZ0(1))**AB
6355 GZ2=(1.D0-GZ0(2))**AB
6356 GZ3=(1.D0-CC(2)*GZ0(2)-2.D0*(1.D0-CC(2))*GZ0(1))**AB
6357
6358
6359 GZ(1)=CC(ICZ)**2*(GZ2-GZ3)
6360 GZ(2)=CC(ICZ)*(1.D0-CC(ICZ))*(1.D0+GZ2-2.D0*GZ1)
6361 GZ(3)=CC(ICZ)*(1.D0-GZ2)
6362 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6363 202 FORMAT(2X,'XXFAU - END')
6364 RETURN
6365 END
6366
6367
6368 SUBROUTINE XXFRAG(SA,NA,RC)
6369
6370
6371
6372 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6373 INTEGER DEBUG
6374
6375 DIMENSION SA(210,3)
6376 COMMON /Q_AREA13/ NSF,IAF(56)
6377 COMMON /Q_AREA43/ MONIOU
6378 COMMON /Q_DEBUG/ DEBUG
6379 SAVE
6380
6381 IF(DEBUG.GE.2)WRITE (MONIOU,201)NA
6382 201 FORMAT(2X,'XXFRAG-MULTIFRAGMENTATION: NUCLEUS MASS NUMBER: NA='
6383 * ,I2)
6384 IF(DEBUG.GE.3)THEN
6385 WRITE (MONIOU,203)
6386 203 FORMAT(2X,'NUCLEONS COORDINATES:')
6387 204 FORMAT(2X,3E10.3)
6388 DO 205 I=1,NA
6389 205 WRITE (MONIOU,204)(SA(I,L),L=1,3)
6390 ENDIF
6391
6392 NI=1
6393 NG=1
6394 J=0
6395 1 J=J+1
6396 J1=NI+1
6397 DO 4 I=J1,NA
6398 RI=0.D0
6399 DO 2 M=1,3
6400 2 RI=RI+(SA(J,M)-SA(I,M))**2
6401 IF(RI.GT.RC)GOTO 4
6402 NI=NI+1
6403 NG=NG+1
6404 IF(I.EQ.NI)GOTO 4
6405 DO 3 M=1,3
6406 S0=SA(NI,M)
6407 SA(NI,M)=SA(I,M)
6408 3 SA(I,M)=S0
6409 4 CONTINUE
6410 IF(J.LT.NI.AND.NA-NI.GT.0)GOTO 1
6411 NSF=NSF+1
6412 IAF(NSF)=NG
6413 IF(DEBUG.GE.3)WRITE (MONIOU,206)NSF,IAF(NSF)
6414 206 FORMAT(2X,'XXFRAG: FRAGMENT N',I2,2X,'FRAGMENT MASS - ',I2)
6415 NG=1
6416 J=NI
6417 NI=NI+1
6418 IF(NA-NI)6,5,1
6419 5 NSF=NSF+1
6420 IAF(NSF)=1
6421 IF(DEBUG.GE.3)WRITE (MONIOU,206)NSF,IAF(NSF)
6422 6 CONTINUE
6423 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6424 202 FORMAT(2X,'XXFRAG - END')
6425 RETURN
6426 END
6427
6428
6429 SUBROUTINE XXFRAGM(NS,XA)
6430
6431
6432
6433
6434 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6435
6436 DIMENSION XA(210,3)
6437 INTEGER DEBUG
6438 COMMON /Q_AREA1/ IA(2),ICZ,ICP
6439 COMMON /Q_AREA3/ RMIN,EMAX,EEV
6440 COMMON /Q_AREA11/ B10
6441
6442
6443 COMMON /Q_AREA13/ NSF,IAF(56)
6444 COMMON /Q_AREA43/ MONIOU
6445 COMMON /Q_DEBUG/ DEBUG
6446 SAVE
6447 EXTERNAL PSRAN
6448
6449 IF(DEBUG.GE.2)WRITE (MONIOU,201)NS
6450 201 FORMAT(2X,'XXFRAGM: NUMBER OF SPECTATORS: NS=',I2)
6451
6452 NSF=0
6453
6454 IF(NS-1)6,1,2
6455
6456 1 NSF=NSF+1
6457 IAF(NSF)=1
6458 IF(DEBUG.GE.3)WRITE (MONIOU,205)
6459 205 FORMAT(2X,'XXFRAGM - SINGLE SPECTATOR')
6460 GOTO 6
6461 2 EEX=0.D0
6462
6463
6464 DO 3 I=1,IA(1)-NS
6465
6466
6467 3 EEX=EEX+(PSRAN(B10)+PSRAN(B10)+PSRAN(B10)+
6468 * PSRAN(B10)+PSRAN(B10)-2.5D0)**2*2.4D0
6469 IF(DEBUG.GE.3)WRITE (MONIOU,203)EEX
6470 203 FORMAT(2X,'XXFRAGM: EXCITATION ENERGY: EEX=',E10.3)
6471
6472
6473
6474 IF(EEX/NS.GT.EMAX)THEN
6475
6476 CALL XXFRAG(XA,NS,RMIN)
6477 ELSE
6478
6479
6480
6481 NF=IXXSON(NS,EEX/EEV,PSRAN(B10))
6482 NSF=NSF+1
6483
6484 IAF(NSF)=NS-NF
6485 IF(DEBUG.GE.3)WRITE (MONIOU,206)IAF(NSF)
6486 206 FORMAT(2X,'XXFRAGM - EVAPORATION: MASS NUMBER OF THE FRAGMENT:'
6487 * ,I2)
6488
6489
6490
6491 NAL=NF/4
6492 IF(NAL.NE.0)THEN
6493
6494 DO 4 I=1,NAL
6495 NSF=NSF+1
6496 4 IAF(NSF)=4
6497 ENDIF
6498
6499 NF=NF-4*NAL
6500 IF(NF.NE.0)THEN
6501
6502 DO 5 I=1,NF
6503 NSF=NSF+1
6504 5 IAF(NSF)=1
6505 ENDIF
6506 IF(DEBUG.GE.3)WRITE (MONIOU,204)NF,NAL
6507 204 FORMAT(2X,'XXFRAGM - EVAPORATION: NUMBER OF NUCLEONS NF=',I2,
6508 * 'NUMBER OF ALPHAS NAL=',I2)
6509 ENDIF
6510 6 CONTINUE
6511 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6512 202 FORMAT(2X,'XXFRAGM - END')
6513 RETURN
6514 END
6515
6516
6517 SUBROUTINE XXFZ(B,GZ)
6518
6519
6520 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6521 INTEGER DEBUG
6522 DIMENSION GZ(2),FHARD(3)
6523 COMMON /Q_AREA1/ IA(2),ICZ,ICP
6524 COMMON /Q_AREA2/ S,Y0,WP0,WM0
6525 COMMON /Q_AREA7/ RP1
6526 COMMON /Q_AR13/ X1(7),A1(7)
6527 COMMON /Q_AREA43/ MONIOU
6528 COMMON /Q_DEBUG/ DEBUG
6529 SAVE
6530
6531 IF(DEBUG.GE.2)WRITE (MONIOU,201)
6532 201 FORMAT(2X,'XXFZ - HADRONIC CROSS-SECTIONS CALCULATION')
6533
6534 DO 1 L=1,2
6535 1 GZ(L)=0.D0
6536 E1=EXP(-1.D0)
6537
6538 DO 2 I1=1,7
6539 DO 2 M=1,2
6540 Z=.5D0+X1(I1)*(M-1.5D0)
6541 S1=DSQRT(RP1*Z)
6542 ZV1=EXP(-Z)
6543 S2=DSQRT(RP1*(1.D0-DLOG(Z)))
6544 ZV2=E1*Z
6545
6546
6547
6548
6549
6550
6551 VV1=EXP(-PSFAZ(ZV1,FSOFT,FHARD,FSHARD)-FHARD(1)
6552 * -FHARD(2)-FHARD(3))
6553 VV2=EXP(-PSFAZ(ZV2,FSOFT,FHARD,FSHARD)-FHARD(1)
6554 * -FHARD(2)-FHARD(3))
6555
6556
6557 IF(IA(2).EQ.1)THEN
6558 CG1=1.D0
6559 CG2=1.D0
6560 ELSE
6561 CG1=XXROT(B,S1)
6562 CG2=XXROT(B,S2)
6563 ENDIF
6564
6565 DO 2 L=1,2
6566 2 GZ(L)=GZ(L)+ A1(I1)*(CG1*(1.D0-VV1**L)+CG2*(1.D0-VV2**L)/Z)
6567 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6568 202 FORMAT(2X,'XXFZ - END')
6569 RETURN
6570 END
6571
6572
6573 SUBROUTINE XXGAU(GZ)
6574
6575
6576
6577 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6578 INTEGER DEBUG
6579 DIMENSION GZ(3),GZ0(3)
6580 COMMON /Q_AREA6/ PI,BM,AM
6581 COMMON /Q_AR13/ X1(7),A1(7)
6582 COMMON /Q_AR2/ R,RM
6583 COMMON /Q_AREA43/ MONIOU
6584 COMMON /Q_DEBUG/ DEBUG
6585 SAVE
6586
6587 IF(DEBUG.GE.2)WRITE (MONIOU,201)
6588 201 FORMAT(2X,'XXGAU - NUCLEAR CROSS-SECTIONS CALCULATION')
6589
6590 DO 1 I=1,3
6591 1 GZ(I)=0.D0
6592
6593 DO 2 I=1,7
6594 DO 2 M=1,2
6595 B=BM*DSQRT(.5D0+X1(I)*(M-1.5D0))
6596 CALL XXFAU(B,GZ0)
6597 DO 2 L=1,3
6598 2 GZ(L)=GZ(L)+GZ0(L)*A1(I)
6599 DO 3 L=1,3
6600 3 GZ(L)=GZ(L)*(BM*AM)**2*PI*.5D0
6601 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6602 202 FORMAT(2X,'XXGAU - END')
6603 RETURN
6604 END
6605
6606
6607 SUBROUTINE XXGAU1(GZ)
6608
6609
6610
6611 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6612 INTEGER DEBUG
6613 DIMENSION GZ(3),GZ0(3)
6614 COMMON /Q_AREA6/ PI,BM,AM
6615 COMMON /Q_AR15/ X5(2),A5(2)
6616 COMMON /Q_AR2/ R,RM
6617 COMMON /Q_AREA43/ MONIOU
6618 COMMON /Q_DEBUG/ DEBUG
6619 SAVE
6620
6621 IF(DEBUG.GE.2)WRITE (MONIOU,201)
6622 201 FORMAT(2X,'XXGAU1 - NUCLEAR CROSS-SECTIONS CALCULATION')
6623
6624 DO 1 I=1,2
6625 B=BM+X5(I)
6626 CALL XXFAU(B,GZ0)
6627 DO 1 L=1,3
6628 1 GZ(L)=GZ(L)+GZ0(L)*A5(I)*EXP(X5(I))*B*2.D0*PI*AM*AM
6629 IF(DEBUG.GE.3)WRITE (MONIOU,202)
6630 202 FORMAT(2X,'XXGAU1 - END')
6631 RETURN
6632 END
6633
6634
6635 SUBROUTINE XXGENER(WP0,WM0,EY0,S0X,C0X,S0,C0,IC1,IC2)
6636
6637
6638
6639
6640
6641
6642
6643
6644 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
6645 INTEGER DEBUG
6646 CHARACTER *2 TYQ
6647 DIMENSION WP(2),IC(2),EPT(4),EP(4),EY(3),EY0(3)
6648
6649
6650 COMMON /Q_AREA8/ WWM,BEP,BEN,BEK,BEC,DC(5),DETA,ALMPT
6651 COMMON /Q_AREA10/ STMASS,AM0,AMN,AMK,AMC,AMLAMC,AMLAM,AMETA
6652 COMMON /Q_AREA11/ B10
6653 COMMON /Q_AREA19/ AHL(5)
6654
6655 COMMON /Q_AREA21/ DMMIN(5)
6656
6657 COMMON /Q_AREA28/ ARR(4)
6658 COMMON /Q_AREA42/ TYQ(15)
6659 COMMON /Q_AREA43/ MONIOU
6660 COMMON /Q_DEBUG/ DEBUG
6661 SAVE
6662 EXTERNAL PSRAN
6663
6664 IF(DEBUG.GE.2)WRITE (MONIOU,201)TYQ(8+IC1),TYQ(8+IC2),
6665 * WP0,WM0,EY0,S0X,C0X,S0,C0
6666 201 FORMAT(2X,'XXGENER: PARTON FLAVORS AT THE ENDS OF THE STRING:',
6667 * 2X,A2,2X,A2/4X,'LIGHT CONE MOMENTA OF THE STRING: ',E10.3,
6668 * 2X,E10.3/4X,'EY0=',3E10.3/4X,
6669 * 'S0X=',E10.3,2X,'C0X=',E10.3,2X,'S0=',E10.3,2X,'C0=',E10.3)
6670
6671 WW=WP0*WM0
6672 EPT(1)=.5D0*(