C======================================================================C C C C QQQ GGG SSSS JJJJJJJ EEEEEEE TTTTTTT C C Q Q G G S S J E T C C Q Q G S J E T C C Q Q G GGG SSSS J EEEEE T C C Q Q Q G G S J E T C C Q Q G G S S J J E T C C QQQ QQ GGG SSSS JJJ EEEEEEE T C C C C C C----------------------------------------------------------------------C C C C QUARK - GLUON - STRING - MODEL C C C C HIGH ENERGY HADRON INTERACTION PROGRAM C C C C BY C C C C N. N. KALMYKOV AND S. S. OSTAPCHENKO C C C C MOSCOW STATE UNIVERSITY, MOSCOW, RUSSIA C C e-mail: serg@eas.npi.msu.su C C----------------------------------------------------------------------C C SUBROUTINE VERSION TO BE LINKED WITH C C C O R S I K A C C KARLSRUHE AIR SHOWER SIMULATION PROGRAM C C WITH MODIFICATIONS C C BY C C D. HECK IK3 FZK KARLSRUHE C C----------------------------------------------------------------------C C last modification: feb 21, 1997 C C----------------------------------------------------------------------C C======================================================================= SUBROUTINE PSAINI c Common initialization procedure c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG CHARACTER *7 TY LOGICAL LCALC,LSECT ******************************************** DIMENSION EQ(17),MIJ(17,17,4),NIJ(17,17,4),CSJET(17,17,68), *CS1(17,17,68),GZ0(2),GZ1(3) COMMON /XSECT/ GSECT(10,5,4) COMMON /AREA1/ IA(2),ICZ,ICP COMMON /AREA5/ RD(2),CR1(2),CR2(2),CR3(2) ******************************************** COMMON /AREA6/ PI,BM,AM COMMON /AREA7/ RP1 COMMON /AREA10/ STMASS,AM0,AMN,AMK,AMC,AMLAMC,AMLAM,AMETA COMMON /AREA15/ FP(5),RQ(5),CD(5) COMMON /AREA16/ CC(5) COMMON /AREA17/ DEL,RS,RS0,FS,ALFP,RR,SH,DELH COMMON /AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0 COMMON /AREA19/ AHL(5) ******************************************** COMMON /AREA22/ SJV0,FJS0(5,3) ******************************************** COMMON /AREA23/ RJV(50) COMMON /AREA24/ RJS(50,5,10) COMMON /AREA27/ FP0(5) COMMON /AREA28/ ARR(4) COMMON /AREA29/ CSTOT(17,17,68) COMMON /AREA30/ CS0(17,17,68) COMMON /AREA31/ CSBORN(17,68) COMMON /AREA32/ CSQ(17,2,2),CSBQ(17,2,2) COMMON /AREA33/ FSUD(10,2) COMMON /AREA34/ QRT(10,101,2) COMMON /AREA35/ SJV(10,5),FJS(10,5,15) COMMON /AREA39/ JCALC COMMON /AREA41/ TY(5) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG ******************************************** COMMON /AREA44/ GZ(10,5,4) c Auxiliary common blocks to calculate hadron-nucleus cross-sections COMMON /AR1/ ANORM COMMON /AR2/ RRR,RRRM ******************************************** c------------------------------------------------- WRITE(MONIOU,100) 100 FORMAT(' ', * '====================================================', * /,' ','| |', * /,' ','| QUARK GLUON STRING JET MODEL |', * /,' ','| |', * /,' ','| HADRONIC INTERACTION MONTE CARLO |', * /,' ','| BY |', * /,' ','| N.N. KALMYKOV AND S.S. OSTAPCHENKO |', * /,' ','| |', * /,' ','| e-mail: serg@eas.npi.msu.su |', * /,' ','| |', * /,' ','| last modification: feb. 21, 1997 by D. Heck |', * /,' ','====================================================', * /) IF(DEBUG.GE.1)WRITE (MONIOU,210) 210 FORMAT(2X,'PSAINI - MAIN INITIALIZATION PROCEDURE') c AHL(i) - parameter for the energy sharing procedure (govern leading hadronic state c inelasticity for primary pion, nucleon, kaon, D-meson, Lambda_C correspondingly) AHL(1)=1.D0-2.D0*ARR(1) AHL(2)=1.D0-ARR(1)-ARR(2) AHL(3)=1.D0-ARR(1)-ARR(3) AHL(4)=1.D0-ARR(1)-ARR(4) AHL(5)=AHL(2)+ARR(1)-ARR(4) c------------------------------------------------- c 1/CC(i) = C_i - shower enhancement coefficients for one vertex c (C_ab=C_a*C_b) (i - ICZ) CC(2)=1.D0/DSQRT(CD(2)) CC(1)=1.D0/CC(2)/CD(1) CC(3)=1.D0/CC(2)/CD(3) CC(4)=1.D0/CC(2)/CD(4) CC(5)=1.D0/CC(2)/CD(5) c FP0(i) - vertex constant (FP_ij=FP0_i*FP0_j) for pomeron-hadron interaction (i - ICZ) FP0(2)=DSQRT(FP(2)) FP0(1)=FP(1)/FP0(2) FP0(3)=FP(3)/FP0(2) FP0(4)=FP(4)/FP0(2) FP0(5)=FP(5)/FP0(2) c SH - hard interaction effective squared (SH=pi*R_h>2, R_h>2=4/Q0>2) SH=4.D0/QT0*PI c Auxiliary constants for the hard interaction AQT0=DLOG(4.D0*QT0) QLOG=DLOG(QT0/ALM) QLL=DLOG(QLOG) ******************************************** INQUIRE(FILE='QGSDATA4',EXIST=LCALC) IF(LCALC)then IF(DEBUG.GE.1)WRITE (MONIOU,211) 211 FORMAT(2X,'PSAINI: HARD CROSS SECTION RATIOS READOUT FROM THE' * ' FILE QGSDATA4') OPEN(1,FILE='QGSDATA4',STATUS='OLD') READ (1,*)CSBORN,CS0,CSTOT,CSQ,CSBQ, * FSUD,QRT,SJV,FJS,RJV,RJS,GZ,GSECT CLOSE(1) ELSE ******************************************** IF(DEBUG.GE.1)WRITE (MONIOU,201) 201 FORMAT(2X,'PSAINI: HARD CROSS SECTIONS CALCULATION') c-------------------------------------------------- c Hard pomeron inclusive cross sections calculation c-------------------------------------------------- c EQ(I) - energy squared tabulation (Q0>2, 4*Q0>2, ...) DO 1 I=1,17 1 EQ(I)=QT0*4.D0**FLOAT(I-1) DO 2 I=1,17 c QI - effective momentum (Qt**2/(1-z)**2) cutoff for the Born process QI=EQ(I) c M, L define parton types (1-g, 2-q) DO 2 M=1,2 DO 2 L=1,2 c K defines c.m. energy squared for the process (for current energy tabulation) DO 2 K=1,17 K1=K+17*(M-1)+34*(L-1) IF(K.LE.I.OR.K.EQ.2)THEN CSBORN(I,K1)=0.D0 ELSE c SK - c.m. energy squared for the hard interaction SK=EQ(K) c CSBORN(I,K1) - Born cross-section (2->2 process) - procedure PSBORN CSBORN(I,K1)=PSBORN(QI,SK,M-1,L-1) ENDIF 2 CONTINUE c Cross-sections initialization DO 3 I=1,17 DO 3 J=1,17 N=MAX(I,J) DO 3 M=1,2 DO 3 L=1,2 ML=M+2*L-2 DO 3 K=1,17 K1=K+17*(M-1)+34*(L-1) CSJET(I,J,K1)=0.D0 IF(K.LE.N.OR.K.EQ.2)THEN CSTOT(I,J,K1)=-80.D0 CS0(I,J,K1)=-80.D0 MIJ(I,J,ML)=K+1 NIJ(I,J,ML)=K+1 ELSE CSTOT(I,J,K1)=DLOG(CSBORN(N,K1)) CS0(I,J,K1)=CSTOT(I,J,K1) ENDIF 3 CONTINUE c N-maximal number of ladder runs taken into account N=2 4 CONTINUE IF(DEBUG.GE.2)WRITE (MONIOU,202)N,EQ(MIJ(1,1,1)),EQ(NIJ(1,1,1)) 202 FORMAT(2X,'PSAINI: NUMBER OF LADDER RUNS TO BE CONSIDERED:',I2/ * 4X,'MINIMAL MASSES SQUARED FOR THE UNORDERED AND STRICTLY', * ' ORDERED LADDERS:'/4X,E10.3,3X,E10.3) DO 6 I=1,17 c QI - effective momentum cutoff for upper end of the ladder QI=EQ(I) DO 6 J=1,17 c QJ - effective momentum cutoff for lower end of the ladder QJ=EQ(J) c QQ - maximal effective momentum cutoff QQ=MAX(QI,QJ) c S2MIN - minimal energy squared for 2->2 subprocess S2MIN=MAX(QQ,4.D0*QT0) SM=DSQRT(QT0/S2MIN) c SMIN - minimal energy squared for 2->3 subprocess SMIN=S2MIN*(1.D0+SM)/(1.D0-SM) c M, L define parton types (1-g, 2-q) DO 6 M=1,2 DO 6 L=1,2 ML=M+2*L-2 c KMIN corresponds to minimal energy at which more runs are to be considered - c stored in array NIJ(I,J,ML) - for strictly ordered ladder KMIN=NIJ(I,J,ML) IF(KMIN.LE.17)THEN DO 5 K=KMIN,17 SK=EQ(K) IF(SK.LE.SMIN)THEN NIJ(I,J,ML)=NIJ(I,J,ML)+1 ELSE K1=K+17*(M-1)+34*(L-1) c CS1(I,J,K1) - cross-section for strictly ordered ladder (highest virtuality run c is the lowest one) - procedure PSJET1 CS1(I,J,K1)=PSJET1(QI,QJ,SK,S2MIN,M-1,L) ENDIF 5 CONTINUE ENDIF 6 CONTINUE DO 8 I=1,17 DO 8 J=1,17 DO 8 M=1,2 DO 8 L=1,2 ML=M+2*L-2 KMIN=NIJ(I,J,ML) IF(KMIN.LE.17)THEN DO 7 K=KMIN,17 K1=K+17*(M-1)+34*(L-1) c CSJ - cross-section for strictly ordered ladder (highest virtuality run is the c lowest one) - Born contribution is added CSJ=CS1(I,J,K1)+CSBORN(MAX(I,J),K1) IF(DEBUG.GE.2)WRITE (MONIOU,204)CSJ,EXP(CS0(I,J,K1)) 204 FORMAT(2X,'PSAINI: NEW AND OLD VALUES OF THE CONTRIBUTION', * ' OF THE STRICTLY ORDERED LADDER:'/4X,E10.3,3X,E10.3) IF(CSJ.EQ.0.D0.OR.ABS(1.D0-EXP(CS0(I,J,K1))/CSJ).LT.1.D-2)THEN NIJ(I,J,ML)=NIJ(I,J,ML)+1 ELSE c CS0(I,J,K1) - cross-section logarithm for strictly ordered ladder CS0(I,J,K1)=DLOG(CSJ) ENDIF 7 CONTINUE ENDIF 8 CONTINUE DO 10 I=1,17 QI=EQ(I) DO 10 J=1,17 QJ=EQ(J) QQ=MAX(QI,QJ) S2MIN=MAX(QQ,4.D0*QT0) SM=DSQRT(QT0/S2MIN) c SMIN - minimal energy squared for 2->3 subprocess SMIN=S2MIN*(1.D0+SM)/(1.D0-SM) DO 10 M=1,2 DO 10 L=1,2 ML=M+2*L-2 c KMIN corresponds to minimal energy at which more runs are to be considered c stored in array MIJ(I,J,ML) - for any ordering in the ladder KMIN=MIJ(I,J,ML) IF(KMIN.LE.17)THEN DO 9 K=KMIN,17 SK=EQ(K) IF(SK.LE.SMIN)THEN MIJ(I,J,ML)=MIJ(I,J,ML)+1 ELSE K1=K+17*(M-1)+34*(L-1) c CS1(I,J,K1) - cross-section for any ordering in the ladder (highest virtuality c run is somewhere in the middle; runs above and below it are strictly ordered c towards highest effective momentum run) - procedure PSJET CS1(I,J,K1)=PSJET(QI,QJ,SK,S2MIN,M-1,L) ENDIF 9 CONTINUE ENDIF 10 CONTINUE DO 12 I=1,17 DO 12 J=1,17 DO 12 M=1,2 DO 12 L=1,2 ML=M+2*L-2 c KMIN corresponds to minimal energy at which more runs are to be considered KMIN=MIJ(I,J,ML) IF(KMIN.LE.17)THEN DO 11 K=KMIN,17 K1=K+17*(M-1)+34*(L-1) K2=K+17*(L-1)+34*(M-1) CSJ=CS1(I,J,K1)+EXP(CS0(J,I,K2)) IF(CSJ.EQ.0.D0.OR.ABS(1.D0-EXP(CSTOT(I,J,K1))/CSJ).LT.1.D-2) * MIJ(I,J,ML)=MIJ(I,J,ML)+1 IF(DEBUG.GE.2)WRITE (MONIOU,203)CSJ,EXP(CSTOT(I,J,K1)) 203 FORMAT(2X,'PSAINI: NEW AND OLD VALUES OF THE UNORDERED LADDER', * ' CROSS SECTION:'/4X,E10.3,3X,E10.3) 11 CSTOT(I,J,K1)=DLOG(CSJ) ENDIF 12 CONTINUE c One more run N=N+1 DO 13 L=1,4 13 IF(MIJ(1,1,L).LE.17.OR.NIJ(1,1,L).LE.17)GOTO 4 c Logarithms of the Born cross-section are calculated - to be interpolated in the c PSBINT procedure DO 14 I=1,17 DO 14 K=1,17 DO 14 M=1,2 DO 14 L=1,2 K1=K+17*(M-1)+34*(L-1) IF(K.LE.I.OR.K.EQ.2)THEN CSBORN(I,K1)=-80.D0 ELSE CSBORN(I,K1)=DLOG(CSBORN(I,K1)) ENDIF 14 CONTINUE c Total and Born hard cross-sections logarithms for minimal cutoff (QT0) - to be c interpolated in the PSJINT0 procedure DO 15 M=1,2 DO 15 L=1,2 DO 15 K=1,17 IF(K.LE.2)THEN CSQ(K,M,L)=-80.D0 CSBQ(K,M,L)=-80.D0 ELSE K1=K+17*(M-1)+34*(L-1) CSBQ(K,M,L)=CSBORN(1,K1) CSQ(K,M,L)=CSTOT(1,1,K1) ENDIF 15 CONTINUE c------------------------------------------------- c FSUD(K,M)=-ln(SUD) - timelike Sudakov formfactor logarithm - procedure c PSUDT(QMAX,M-1), M=1 - g, M=2 - q DO 17 M=1,2 FSUD(1,M)=0.D0 DO 17 K=2,10 c QMAX is the maximal effective momentum ( Qt**2/z**2/(1-z)**2 in case of the timelike c evolution ) QMAX=QTF*4.D0**(1.D0+K) 17 FSUD(K,M)=PSUDT(QMAX,M-1) c QRT(K,L,M) - effective momentum logarithm for timelike branching ( ln QQ/16/QTF ) c for given QMAX (defined by K, QLMAX = ln QMAX/16/QTF ) and a number c of random number values (defined by L) - to be interpolated by the PSQINT c procedure; M=1 - g, M=2 - q DO 18 M=1,2 DO 18 K=1,10 QLMAX=1.38629D0*(K-1) QRT(K,1,M)=0.D0 QRT(K,101,M)=QLMAX DO 18 I=1,99 IF(K.EQ.1)THEN QRT(K,I+1,M)=0.D0 ELSE QRT(K,I+1,M)=PSROOT(QLMAX,.01D0*I,M) ENDIF 18 CONTINUE c------------------------------------------------- IF(DEBUG.GE.2)WRITE (MONIOU,205) 205 FORMAT(2X,'PSAINI: PRETABULATION OF THE INTERACTION EIKONALS') c------------------------------------------------- ************************************************************************ c------------------------------------------------- c Interaction cross sections c Factors for interaction eikonals calculation c (convolution of the hard cross-sections with partons structure functions) c - to be used in the PSPSFAZ procedure c------------------------------------------------- IA(1)=1 c------------------------------------------------- DO 21 IE=1,10 c Energy of the interaction (per nucleon) E0N=10.D0**IE c------------------------------------------------- c Energy dependent factors: c WP0, WM0 - initial light cone momenta for the interaction (E+-p) S=2.D0*E0N*AMN c Y0 - total rapidity range for the interaction Y0=DLOG(S) c Type of the incident hadron (icz = 1: pion, 2: nucleon, 3: kaon, etc DO 21 ICZ=1,5 c RS - soft pomeron elastic scattering slope (lambda_ab) RS=RQ(ICZ)+ALFP*Y0 c RS0 - initial slope (sum of the pomeron-hadron vertices slopes squared - R_ab) RS0=RQ(ICZ) c FS - factor for pomeron eikonal calculation c (gamma_ab * s**del /lambda_ab * C_ab FS=FP(ICZ)*EXP(Y0*DEL)/RS*CD(ICZ) c RP1 - factor for the impact parameter dependence of the eikonal ( in fm>2 ) RP1=RS*4.D0*.0391D0/AM**2 c Factor for cross-sections calculation ( in mb ) G0=PI*RP1/CD(ICZ)*AM**2*10.D0 c SJV - valence-valence cross-section (divided by 8*pi*lambda_ab) SJV(IE,ICZ)=PSHARD(S,ICZ) SJV0=SJV(IE,ICZ) DO 19 I=1,5 DO 19 M=1,3 Z=.2D0*I c Eikonals for gluon-gluon and valence-gluon semihard interactions c (m=1 - gg, 2 - qg, 3 - gq); c Z - impact parameter factor ( exp(-b**2/R_p) ) M1=M+3*(ICZ-1) FJS(IE,I,M1)=DLOG(PSFSH(S,Z,ICZ,M-1)/Z) FJS0(I,M)=FJS(IE,I,M1) 19 CONTINUE DO 20 IIA=1,4 c Target mass number IA(2) IA(2)=4**(IIA-1) IF(DEBUG.GE.1)WRITE (MONIOU,206)E0N,TY(ICZ),IA(2) 206 FORMAT(2X,'PSAINI: INITIAL PARTICLE ENERGY:',E10.3,2X, *'ITS TYPE:',A7,2X,'TARGET MASS NUMBER:',I2) c------------------------------------------------- c Nuclear radii IF(IA(2).GT.10)THEN c RD - Wood-Saxon density radius (fit to the data of Murthy et al.) RD(2)=0.7D0*FLOAT(IA(2))**.446/AM ELSE c RD - gaussian density radius (for light nucleus) RD(2)=.9D0*FLOAT(IA(2))**.3333/AM ENDIF IF(IA(2).EQ.1)THEN c Hadron-proton interaction c BM - impact parameter cutoff value BM=2.D0*DSQRT(RP1) c XXFZ - impact parameter integration for the hadron-nucleon interaction eikonal; c GZ0 - total and absorptive cross-sections (up to a factor); first parameter is c used only in case of hadron-nucleus interaction (to make convolution with target c nucleus profile function) CALL XXFZ(0.D0,GZ0) write (*,*)gz0 c GTOT - total cross-section GTOT=G0*GZ0(1) c GABS - cut pomerons cross-section GABS=G0*GZ0(2)*.5D0 c GD0 - cross-section for the cut between pomerons GD0=GTOT-GABS c GDP - projectile diffraction cross section GDP=(1.D0-CC(ICZ))*CC(2)*GD0 c GDT - target diffraction cross section GDT=(1.D0-CC(2))*CC(ICZ)*GD0 c GDD - double diffractive cross section GDD=(1.D0-CC(ICZ))*(1.D0-CC(2))*GD0 c GIN - inelastic cross section GIN=GABS+GDP+GDT+GDD GEL=GD0*CC(ICZ)*CC(2) c IF(DEBUG.GE.1)WRITE (MONIOU,225)GTOT,GIN,GEL,GDP,GDT,GDD c 225 FORMAT(2X,'PSAINI: HADRON-PROTON CROSS SECTIONS:'/ * 4X,'GTOT=',E10.3,2X,'GIN=',E10.3,2X,'GEL=',E10.3/4X, * 'GDIFR_PROJ=',E10.3,2X,'GDIFR_TARG=',E10.3,2X, * 'G_DOUBLE_DIFR',E10.3) c GZ - probability to have target diffraction GZ(IE,ICZ,IIA)=GDT/GIN C?????? GSECT(IE,ICZ,IIA)=LOG(GIN) C?????? ELSE c Hadron-nucleus interaction c BM - impact parameter cutoff value BM=RD(2)+DLOG(29.D0) c RRR - Wood-Saxon radius for the target nucleus RRR=RD(2) c RRRM - auxiliary parameter for numerical integration RRRM=RRR+DLOG(9.D0) c ANORM - nuclear density normalization factor multiplied by RP1 ANORM=1.5D0/PI/RRR**3/(1.D0+(PI/RRR)**2)*RP1 c GAU(GZ) - cross sections calculation ( integration over impact parameters less than c BM ) CALL XXGAU(GZ1) c GAU1(GZ) - cross sections calculation ( integration over impact c parameters greater than BM ) CALL XXGAU1(GZ1) c GIN - total inelastic cross section GIN=GZ1(1)+GZ1(2)+GZ1(3) c IF(DEBUG.GE.1)WRITE (MONIOU,224) * GIN*10.D0,GZ1(1)*10.D0,GZ1(2)*10.D0 c 224 FORMAT(2X,'PSAINI: HADRON-NUCLEUS CROSS SECTIONS:'/ * 4X,'GIN=',E10.3,2X,'GDIFR_TARG=',E10.3,2X, * 'GDIFR_PROJ=',E10.3) c GZ - probability to have target diffraction GZ(IE,ICZ,IIA)=GZ1(1)/GIN C?????? GIN=GIN*10. GSECT(IE,ICZ,IIA)=LOG(GIN) C?????? ENDIF 20 CONTINUE 21 CONTINUE c Rejection functions calculation - to be interpolated in the RJINT procedure DO 23 I=1,50 c Rapidity range tabulation for the hard interaction YJ=AQT0+.5D0*I c Rejection function for valence quark energy distribution RJV(I)=PSREJV(EXP(YJ)) DO 22 J=1,5 DO 22 M=1,2 Z=.2D0*J DO 22 ICZ=1,5 c RS0 - initial slope (sum of the pomeron-hadron vertices slopes squared - R_ab) RS0=RQ(ICZ) M1=M+2*(ICZ-1) c Rejection function for semihard block energy distribution (m=1 - gg, c 2 - qg) RJS(I,J,M1)=PSREJS(EXP(YJ),Z,M-1) 22 CONTINUE 23 CONTINUE IF(DEBUG.GE.1)WRITE (MONIOU,212) 212 FORMAT(2X,'PSAINI: HARD CROSS SECTIONS ARE WRITTEN TO THE FILE' * ,' QGSDATA4') OPEN(1,FILE='QGSDATA4',STATUS='unknown') WRITE (1,*)CSBORN,CS0,CSTOT,CSQ,CSBQ, * FSUD,QRT,SJV,FJS,RJV,RJS,GZ,GSECT CLOSE(1) ENDIF ************************************************************************ IF(DEBUG.GE.3)WRITE (MONIOU,218) 218 FORMAT(2X,'PSAINI - END') RETURN END C======================================================================= FUNCTION PSAPINT(X,J,L) c PSAPINT - integrated Altarelli-Parisi function c X - light cone momentum share value, c J - type of initial parton (0 - g, 1 - q) c L - type of final parton (0 - g, 1 - q) C----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG SAVE IF(DEBUG.GE.2)WRITE (MONIOU,201)X,J,L 201 FORMAT(2X,'PSAPINT: X=',E10.3,2X,'J= ',I1,2X,'L= ',I1) IF(J.EQ.0)THEN IF(L.EQ.0)THEN PSAPINT=6.D0*(DLOG(X/(1.D0-X))-X**3/3.D0+X**2/2.D0-2.D0*X) ELSE PSAPINT=3.D0*(X+X**3/1.5D0-X*X) ENDIF ELSE IF(L.EQ.0)THEN PSAPINT=(DLOG(X)-X+.25D0*X*X)/.375D0 ELSE Z=1.D0-X PSAPINT=-(DLOG(Z)-Z+.25D0*Z*Z)/.375D0 ENDIF ENDIF IF(DEBUG.GE.2)WRITE (MONIOU,202)PSAPINT 202 FORMAT(2X,'PSAPINT=',E10.3) RETURN END C======================================================================= SUBROUTINE PSASET c Common model parameters setting c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG CHARACTER*7 TY COMMON /AREA15/ FP(5),RQ(5),CD(5) COMMON /AREA17/ DEL,RS,RS0,FS,ALFP,RR,SH,DELH COMMON /AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0 COMMON /AREA25/ AHV(5) COMMON /AREA26/ FACTORK COMMON /AREA41/ TY(5) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG IF(DEBUG.GE.1)WRITE (MONIOU,210) 210 FORMAT(2X,'PSASET - COMMON MODEL PARAMETERS SETTING') c Soft pomeron parameters: c DEL - overcriticity, c ALFP - trajectory slope; c FP(i) - vertices for pomeron-hadrons interaction (gamma(i)*gamma(proton)), c RQ(i) - vertices slopes (R(i)**2+R(proton)**2), c CD(i) - shower enhancement coefficients c (i=1,...5 - pion,proton,kaon,D-meson,Lambda_C ), c (Kaidalov et al., Sov.J.Nucl.Phys.,1984 - proton and pion parameters) DEL=.07D0 ALFP=.21D0 FP(1)=2.43D0 RQ(1)=2.4D0 CD(1)=1.6D0 FP(2)=3.64D0 RQ(2)=3.56D0 CD(2)=1.5D0 FP(3)=1.75D0 RQ(3)=2.D0 CD(3)=1.7D0 FP(4)=1.21D0 RQ(4)=1.78D0 CD(4)=2.0D0 FP(5)=2.43D0 RQ(5)=2.4D0 CD(5)=2.0D0 c------------------------------------------------- c Hard interaction parameters: c ALM - Lambda_QCD squared, c QT0 - Q**2 cutoff, c RR - vertex constant square for soft pomeron interaction with the hard block (r**2),; c BET - gluon structure function parameter for the soft pomeron ((1-x)**BET), c AMJ0 - jet mass, c QTF - Q**2 cutoff for the timelike evolution, c FACTORK - K-factor value; c DELH is not a parameter of the model; it is used only for energy sharing c procedure - initially energy is shared according to s**DELH dependence c for the hard interaction cross-section and then rejection is used according c to real Sigma_hard(s) dependence. ALM=.04D0 RR=.35D0 QT0=4.D0 BET=1.D0 DELH=0.25D0 AMJ0=0.D0 QTF=.5D0 FACTORK=2.D0 c------------------------------------------------- c Valence quark structure functions for the hard scattering c (~1/sqrt(x)*(1-x)**AHV(i), i=1,...5 corresponds to pion, nucleon etc.) AHV(1)=1.5D0 AHV(2)=2.5D0 AHV(3)=2.D0 AHV(4)=4.D0 AHV(5)=5.D0 c Initial particle types TY(1)='pion ' TY(2)='nucleon' TY(3)='kaon ' TY(4)='D-meson' TY(5)='LambdaC' RETURN END C======================================================================= FUNCTION PSBINT(QQ,S,M,L) C PSBINT - Born cross-section interpolation c QQ - effective momentum cutoff for the scattering, c S - total c.m. energy squared for the scattering, c M - parton type at current end of the ladder (1 - g, 2 - q) c L - parton type at opposite end of the ladder (1 - g, 2 - q) C----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION WI(3),WK(3) COMMON /AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0 COMMON /AREA31/ CSJ(17,68) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG SAVE IF(DEBUG.GE.2)WRITE (MONIOU,201)QQ,S,M,L 201 FORMAT(2X,'PSBINT: QQ=',E10.3,2X,'S= ',E10.3,2X,'M= ',I1,2X, * 'L= ',I1) PSBINT=0.D0 IF(S.LE.MAX(4.D0*QT0,QQ))THEN IF(DEBUG.GE.3)WRITE (MONIOU,202)PSBINT 202 FORMAT(2X,'PSBINT=',E10.3) RETURN ENDIF ML=17*(M-1)+34*(L-1) QLI=DLOG(QQ/QT0)/1.38629d0 SL=DLOG(S/QT0)/1.38629d0 SQL=SL-QLI I=INT(QLI) K=INT(SL) IF(I.GT.13)I=13 IF(SQL.GT.10.D0)THEN IF(K.GT.14)K=14 WI(2)=QLI-I WI(3)=WI(2)*(WI(2)-1.D0)*.5D0 WI(1)=1.D0-WI(2)+WI(3) WI(2)=WI(2)-2.D0*WI(3) WK(2)=SL-K WK(3)=WK(2)*(WK(2)-1.D0)*.5D0 WK(1)=1.D0-WK(2)+WK(3) WK(2)=WK(2)-2.D0*WK(3) DO 1 I1=1,3 DO 1 K1=1,3 1 PSBINT=PSBINT+CSJ(I+I1,K+K1+ML)*WI(I1)*WK(K1) PSBINT=EXP(PSBINT) ELSEIF(SQL.LT.1.D0.AND.I.NE.0)THEN SQ=(S/QQ-1.D0)/3.D0 WI(2)=QLI-I WI(3)=WI(2)*(WI(2)-1.D0)*.5D0 WI(1)=1.D0-WI(2)+WI(3) WI(2)=WI(2)-2.D0*WI(3) DO 2 I1=1,3 I2=I+I1 K2=I2+1+ML 2 PSBINT=PSBINT+CSJ(I2,K2)*WI(I1) PSBINT=EXP(PSBINT)*SQ ELSEIF(K.EQ.1)THEN SQ=(S/QT0/4.D0-1.D0)/3.D0 WI(2)=QLI WI(1)=1.D0-QLI DO 3 I1=1,2 3 PSBINT=PSBINT+CSJ(I1,3+ML)*WI(I1) PSBINT=EXP(PSBINT)*SQ ELSEIF(K.LT.15)THEN KL=INT(SQL) IF(I+KL.GT.12)I=12-KL IF(I+KL.EQ.1)KL=2 WI(2)=QLI-I WI(3)=WI(2)*(WI(2)-1.D0)*.5D0 WI(1)=1.D0-WI(2)+WI(3) WI(2)=WI(2)-2.D0*WI(3) WK(2)=SQL-KL WK(3)=WK(2)*(WK(2)-1.D0)*.5D0 WK(1)=1.D0-WK(2)+WK(3) WK(2)=WK(2)-2.D0*WK(3) DO 4 I1=1,3 I2=I+I1 DO 4 K1=1,3 K2=I2+KL+K1-1+ML 4 PSBINT=PSBINT+CSJ(I2,K2)*WI(I1)*WK(K1) PSBINT=EXP(PSBINT) ELSE K=15 IF(I.GT.K-3)I=K-3 WI(2)=QLI-I WI(3)=WI(2)*(WI(2)-1.D0)*.5D0 WI(1)=1.D0-WI(2)+WI(3) WI(2)=WI(2)-2.D0*WI(3) WK(2)=SL-K WK(1)=1.D0-WK(2) DO 5 I1=1,3 DO 5 K1=1,2 5 PSBINT=PSBINT+CSJ(I+I1,K+K1+ML)*WI(I1)*WK(K1) PSBINT=EXP(PSBINT) ENDIF IF(DEBUG.GE.3)WRITE (MONIOU,202)PSBINT RETURN END C======================================================================= FUNCTION PSBORN(QQ,S,IQ1,IQ2) c PSFBORN -hard 2->2 parton scattering Born cross-section c S is the c.m. energy square for the scattering process, c IQ1 - parton type at current end of the ladder (0 - g, 1,2 - q) c IQ2 - parton type at opposite end of the ladder (0 - g, 1,2 - q) c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG COMMON /AREA6/ PI,BM,AM COMMON /AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0 COMMON /AREA26/ FACTORK COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG COMMON /AR3/ X1(7),A1(7) SAVE IF(DEBUG.GE.2)WRITE (MONIOU,201)QQ,S,IQ1,IQ2 201 FORMAT(2X,'PSBORN: QQ=',E10.3,2X,'S= ',E10.3,2X,'IQ1= ',I1,2X, * 'IQ2= ',I1) TMIN=S*(.5D0-DSQRT(.25D0-QT0/S)) TMIN=MAX(TMIN,S*QQ/(S+QQ)) IF(IQ1*IQ2.EQ.0)THEN IQ=IQ2 ELSE IQ=2 ENDIF PSBORN=0.D0 DO 1 I=1,7 DO 1 M=1,2 T=2.D0*TMIN/(1.D0+2.D0*TMIN/S-X1(I)*(2*M-3)*(1.D0-2.D0*TMIN/S)) QT=T*(1.D0-T/S) FB=PSFBORN(S,T,IQ1,IQ)+PSFBORN(S,S-T,IQ1,IQ) 1 PSBORN=PSBORN+A1(I)*FB/DLOG(QT/ALM)**2*T**2 PSBORN=PSBORN*(.5D0/TMIN-1.D0/S)*FACTORK*PI**3/2.25D0**2/S**2 IF(IQ1.EQ.0.AND.IQ2.EQ.0)PSBORN=PSBORN*.5D0 IF(DEBUG.GE.3)WRITE (MONIOU,202)PSBORN 202 FORMAT(2X,'PSBORN=',E10.3) RETURN END C======================================================================= SUBROUTINE PSCAJET(QQ,IQ1,QV,ZV,QM,IQV,LDAU,LPAR,JQ) c Final state emission process (all branchings as well as parton masses c are determined) C----------------------------------------------------------------------- c QQ - maximal effective momentum transfer for the first branching c IQ1, IQ2 - initial jet flavours in forward and backward direction c (0 - for gluon) c QV(i,j) - effective momentum for the branching of the parton in i-th row c on j-th level (0 - in case of no branching) - to be determined c ZV(i,j) - Z-value for the branching of the parton in i-th row c on j-th level - to be determined c QM(i,j) - mass squared for the parton in i-th row c on j-th level - to be determined c IQV(i,j) - flavour for the parton in i-th row on j-th level c - to be determined c LDAU(i,j) - first daughter row for the branching of the parton in i-th row c on j-th level - to be determined c LPAR(i,j) - the parent row for the parton in i-th row c on j-th level - to be determined IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION QMAX(30,50),IQM(2),LNV(50), * QV(30,50),ZV(30,50),QM(30,50),IQV(30,50), * LDAU(30,49),LPAR(30,50) COMMON /AREA11/ B10 COMMON /AREA18/ ALM,QT0,QLOG,QLL,AQT0,QTF,BET,AMJ0 COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG IF(DEBUG.GE.2)WRITE (MONIOU,201)QQ,IQ1,JQ 201 FORMAT(2X,'PSCAJET: QQ=',E10.3,2X,'IQ1= ',I1,2X,'JQ=',I1) DO 1 I=2,20 1 LNV(I)=0 LNV(1)=1 QMAX(1,1)=QQ IQV(1,1)=IQ1 NLEV=1 NROW=1 2 QLMAX=DLOG(QMAX(NROW,NLEV)/QTF/16.D0) IQ=MIN(1,IABS(IQV(NROW,NLEV)))+1 IF(PSRAN(B10).GT.PSUDINT(QLMAX,IQ))THEN Q=PSQINT(QLMAX,PSRAN(B10),IQ) Z=PSZSIM(Q,IQ) LL=LNV(NLEV+1)+1 LDAU(NROW,NLEV)=LL LPAR(LL,NLEV+1)=NROW LPAR(LL+1,NLEV+1)=NROW LNV(NLEV+1)=LL+1 IF(IQ.NE.1)THEN IF((3-2*JQ)*IQV(NROW,NLEV).GT.0)THEN IQM(1)=0 IQM(2)=IQV(NROW,NLEV) ELSE IQM(2)=0 IQM(1)=IQV(NROW,NLEV) Z=1.D0-Z ENDIF ELSE ********************************************************* WG=PSFAP(Z,0,0) ********************************************************* WG=WG/(WG+PSFAP(Z,0,1)) IF(PSRAN(B10).LT.WG)THEN IQM(1)=0 IQM(2)=0 ELSE IQM(1)=INT(3.D0*PSRAN(B10)+1.D0)*(3-2*JQ) IQM(2)=-IQM(1) ENDIF IF(PSRAN(B10).LT..5D0)Z=1.D0-Z ENDIF QV(NROW,NLEV)=Q ZV(NROW,NLEV)=Z NROW=LL NLEV=NLEV+1 QMAX(NROW,NLEV)=Q*Z**2 QMAX(NROW+1,NLEV)=Q*(1.D0-Z)**2 IQV(NROW,NLEV)=IQM(1) IQV(NROW+1,NLEV)=IQM(2) IF(DEBUG.GE.3)WRITE (MONIOU,203)NLEV,NROW,Q,Z 203 FORMAT(2X,'PSCAJET: NEW BRANCHING AT LEVEL NLEV=',I2, * ' NROW=',I2/4X,' EFFECTIVE MOMENTUM Q=',E10.3,2X,' Z=',E10.3) GOTO 2 ELSE QV(NROW,NLEV)=0.D0 ZV(NROW,NLEV)=0.D0 QM(NROW,NLEV)=AMJ0 IF(DEBUG.GE.3)WRITE (MONIOU,204)NLEV,NROW 204 FORMAT(2X,'PSCAJET: NEW FINAL JET AT LEVEL NLEV=',I2, * ' NROW=',I2) ENDIF 4 CONTINUE IF(NLEV.EQ.1)THEN IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'PSCAJET - END') RETURN ENDIF LPROW=LPAR(NROW,NLEV) IF(LDAU(LPROW,NLEV-1).EQ.NROW)THEN NROW=NROW+1 GOTO 2 ELSE Z=ZV(LPROW,NLEV-1) QM(LPROW,NLEV-1)=Z*(1.D0-Z)*QV(LPROW,NLEV-1) * +QM(NROW-1,NLEV)/Z+QM(NROW,NLEV)/(1.D0-Z) NROW=LPROW NLEV=NLEV-1 IF(DEBUG.GE.3)WRITE (MONIOU,205)NLEV,NROW,QM(LPROW,NLEV) 205 FORMAT(2X,'PSCAJET: JET MASS AT LEVEL NLEV=',I2, * ' NROW=',I2,' - QM=',E10.3) GOTO 4 ENDIF END C======================================================================= SUBROUTINE PSCONF c Simulation of the interaction configuration: impact parameter, nucleons positions, c numbers of cut soft pomerons and semihard blocks, their connections. c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG c XA(56,3),XB(56,3) - arrays for projectile and target nucleons positions recording, c FHARD(i) give the factors to the scattering amplitude due to c valence quark-gluon (i=1), gluon-valence quark (i=2) and c valence quark-valence quark (i=3) interactions DIMENSION XA(56,3),XB(56,3),FHARD(3) COMMON /AREA1/ IA(2),ICZ,ICP COMMON /AREA2/ S,Y0,WP0,WM0 COMMON /AREA6/ PI,BM,AM c Arrays for interaction configuration recording: c LQA(i) (LQB(j)) - numbers of cut soft pomerons, connected to i-th projectile c (j-th target) nucleon (hadron); c LHA(i) (LHB(j)) - the same for hard pomerons numbers; c IAS(k) (IBS(k)) - number (position in array) of the projectile (target) nucleon, c connected to k-th block of soft pomerons; c NQS(k) - number of soft pomerons in k-th block; c IAH(k) (IBH(k)) - number (position in array) of the projectile (target) nucleon, c connected to k-th hard pomeron; c ZH(k) - impact parameter between the two nucleons connected to k-th hard pomeron c (more exactly exp(-b**2/RP1)); c LVA(i)=1 if valence quark from i-th nucleon (i=1 for hadron) is involved into c the hard interaction and LVA(i)=0 otherwise, LVB(j) - similar. COMMON /AREA9/ LQA(56),LQB(56),NQS(1000),IAS(1000),IBS(1000), * LHA(56),LHB(56),ZH(1000),IAH(1000),IBH(1000), * IQH(1000),LVA(56),LVB(56) COMMON /AREA11/ B10 c NSP - number of secondary particles COMMON /AREA12/ NSP COMMON /AREA16/ CC(5) COMMON /AREA40/ JDIFR COMMON /AREA43/ MONIOU ************************************************** COMMON /AREA45/ GDT ************************************************** COMMON /AREA99/ NWT COMMON /DEBUG/ DEBUG SAVE DIMENSION IWT(56) IF(DEBUG.GE.1)WRITE (MONIOU,201) 201 FORMAT(2X,'PSCONF - CONFIGURATION OF THE INTERACTION') NSP=0 IF(IA(1).EQ.1)THEN ************************************************** IF(JDIFR.EQ.1.AND.PSRAN(B10).LT.GDT)THEN c Target diffraction IF(IA(2).NE.1)THEN c ICT - partner target nucleon type (proton - 2 or neutron - 3) ICT=INT(2.5+PSRAN(B10)) ELSE c Target proton ICT=2 ENDIF WPI=WP0 WMI=WM0 c write (*,*)'difr' CALL XXDTG(WPI,WMI,ICP,ICT,0) RETURN ENDIF ************************************************** c For hadron projectile we have given position in transverse plane; c initially primary hadron is positioned at (X,Y)=(0,0) DO 1 I=1,3 1 XA(1,I)=0.D0 ENDIF c------------------------------------------------- c Inelastic interaction at B (Kaidalov proposed 0.5 value for ALMPT-1, c Sov.J.Nucl.Phys.,1987)) ALMPT=1.7D0 c------------------------------------------------- c Parameters for nuclear spectator part fragmentation: c RMIN - coupling radius squared (fm>2), c EMAX - relative critical energy ( divided per mean excitation energy (~12.5 Mev)), c EEV - relative evaporation energy ( divided per mean excitation energy (~12.5 Mev)) RMIN=3.35D0 EMAX=.11D0 EEV=.25D0 c------------------------------------------------- c DMMIN(i) - minimal diffractive mass for low-mass diffraction for pion, nucleon, c kaon, D-meson, Lambda_C corresp. DMMIN(1)=.76D0 DMMIN(2)=1.24D0 DMMIN(3)=.89D0 DMMIN(4)=2.01D0 DMMIN(5)=2.45D0 c Proton, kaon, pion, D-meson, Lambda, Lambda_C, eta masses AMN=.939D0 AMK=.496D0 AM0=.14D0 AMC=1.868D0 AMLAM=1.116D0 AMLAMC=2.27D0 AMETA=.548D0 c------------------------------------------------- c B10 - initial value of the pseudorandom number, c PI - pi-number c AM - diffusive radius for the Saxon-Wood nuclear density parametrization B10=.43876194D0 PI=3.1416D0 AM=.523D0 C STMASS - minimal string mass to produce secondary particles STMASS=4.D0*AM0**2 c Here and below all radii, distances and so on are divided by AM. RMIN=RMIN/AM**2 TYQ(1)='DD' TYQ(2)='UU' TYQ(3)='C ' TYQ(4)='S ' TYQ(5)='UD ' TYQ(6)='D ' TYQ(7)='U ' TYQ(8)='G ' TYQ(9)='u ' TYQ(10)='d ' TYQ(11)='ud' TYQ(12)='s ' TYQ(13)='c ' TYQ(14)='uu' TYQ(15)='dd' IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXASET - END') RETURN END C======================================================================= SUBROUTINE XXDDFR(WP0,WM0,ICP,ICT) c Double diffractive dissociation c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION EP3(4),EP1(4),EP2(4),EY(3) COMMON /AREA1/ IA(2),ICZ,ICP0 COMMON /AREA2/ S,Y0,WP00,WM00 COMMON /AREA8/ WWM,BE(4),DC(5),DETA,ALMPT COMMON /AREA10/ STMASS,AM(7) COMMON /AREA11/ B10 COMMON /AREA21/ DMMIN(5) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG SAVE IF(DEBUG.GE.2)WRITE (MONIOU,201)ICP,ICT,WP0,WM0 201 FORMAT(2X,'XXDDFR - LEADING CLUSTERS HADRONIZATION:' * /4X,'CLUSTER TYPES ICP=',I2,2X, * 'ICT=',I2/4X,'AVAILABLE LIGHT CONE MOMENTA: WP0=',E10.3, * ' WM0=',E10.3) DO 100 I=1,3 100 EY(I)=1.D0 SD0=WP0*WM0 IF(SD0.LT.0.D0)SD0=0.D0 DDMIN1=DMMIN(ICZ) DDMIN2=DMMIN(2) DDMAX1=MIN(5.D0,DSQRT(SD0)-DDMIN2) IF(DDMAX1.LT.DDMIN1)THEN c Registration of too slow "leading" hadron if its energy is insufficient for c diffractive exhitation IF(DSQRT(SD0).LT.AM(ICZ)+AM(2))THEN IF(WP0.GT.0.D0.AND.(AM(ICZ)+AM(2))**2/WP0.LT..5D0*WM00)THEN SD0=(AM(ICZ)+AM(2))**2 WM0=SD0/WP0 ELSE IF(DEBUG.GE.3)WRITE (MONIOU,202) RETURN ENDIF ENDIF EP3(3)=0.D0 EP3(4)=0.D0 XW=XXTWDEC(SD0,AM(ICZ)**2,AM(2)**2) WP1=XW*WP0 WM1=AM(ICZ)**2/WP1 EP3(1)=.5D0*(WP1+WM1) EP3(2)=.5D0*(WP1-WM1) CALL XXREG(EP3,ICP) WM2=WM0-WM1 WP2=AM(2)**2/WM2 EP3(1)=.5D0*(WP2+WM2) EP3(2)=.5D0*(WP2-WM2) CALL XXREG(EP3,ICT) WP0=0.D0 WM0=0.D0 IF(DEBUG.GE.3)WRITE (MONIOU,202) RETURN ENDIF DMASS1=(DDMIN1/(1.D0-PSRAN(B10)*(1.D0-DDMIN1/DDMAX1)))**2 DDMAX2=MIN(5.D0,DSQRT(SD0)-DSQRT(DMASS1)) DMASS2=(DDMIN2/(1.D0-PSRAN(B10)*(1.D0-DDMIN2/DDMAX2)))**2 WPD1=WP0*XXTWDEC(SD0,DMASS1,DMASS2) WMD1=DMASS1/WPD1 WMD2=WM0-WMD1 WPD2=DMASS2/WMD2 IF(ICP.NE.0)IS=IABS(ICP)/ICP IF(ICZ.EQ.5)THEN ICH1=ICP ICH2=0 AMH1=AM(5)**2 AMH2=AM(1)**2 PTMAX=PSLAM(DMASS1,AMH1,AMH2) IF(PTMAX.LT.0.)PTMAX=0. IF(PTMAX.LT.BE(4)**2)THEN 1 PTI=PTMAX*PSRAN(B10) IF(PSRAN(B10).GT.EXP(-DSQRT(PTI)/BE(4)))GOTO 1 ELSE 2 PTI=(BE(4)*DLOG(PSRAN(B10)*PSRAN(B10)))**2 IF(PTI.GT.PTMAX)GOTO 2 ENDIF AMT1=AMH1+PTI AMT2=AMH2+PTI Z=XXTWDEC(DMASS1,AMT1,AMT2) WP1=WPD1*Z WM1=AMT1/WP1 EP3(1)=.5D0*(WP1+WM1) EP3(2)=.5D0*(WP1-WM1) PT=DSQRT(PTI) CALL PSCS(C,S) EP3(3)=PT*C EP3(4)=PT*S CALL XXREG(EP3,ICH1) WP1=WPD1*(1.D0-Z) WM1=AMT2/WP1 EP3(1)=.5D0*(WP1+WM1) EP3(2)=.5D0*(WP1-WM1) EP3(3)=-PT*C EP3(4)=-PT*S CALL XXREG(EP3,ICH2) GOTO 3 ENDIF IF(ICZ.EQ.1)THEN IF(ICP.NE.0)THEN IC1=ICP*(1-3*INT(.5D0+PSRAN(B10))) IC2=-ICP-IC1 ELSE IC1=INT(1.5D0+PSRAN(B10))*(2*INT(.5D0+PSRAN(B10))-1) IC2=-IC1 ENDIF ELSEIF(ICZ.EQ.2)THEN IF(PSRAN(B10).GT..33333D0)THEN IC1=3*IS IC2=ICP-IS ELSE IC1=ICP+4*IS IC2=4*IS-ICP ENDIF ELSEIF(ICZ.EQ.3)THEN IC1=-4*IS IC2=ICP-3*IS ELSEIF(ICZ.EQ.4)THEN IC1=5*IS IC2=ICP-9*IS ENDIF CALL XXGENER(WPD1,WMD1,EY,0.D0,1.D0,0.D0,1.D0,IC1,IC2) 3 CONTINUE IS=IABS(ICT)/ICT IF(PSRAN(B10).GT..33333D0)THEN IC1=3*IS IC2=ICT-IS ELSE IC1=ICT+4*IS IC2=4*IS-ICT ENDIF CALL XXGENER(WPD2,WMD2,EY,0.D0,1.D0,0.D0,1.D0,IC2,IC1) IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXDDFR - END') RETURN END C======================================================================= SUBROUTINE XXDEC2(EP,EP1,EP2,WW,A,B) c Two particle decay c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG dimension ep(4),ep1(4),ep2(4),EY(3) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG COMMON /AREA11/ B10 IF(DEBUG.GE.2)WRITE (MONIOU,201) 201 FORMAT(2X,'XXDEC2 - TWO PARTICLE DECAY') PL=PSLAM(WW,A,B) EP1(1)=DSQRT(PL+A) EP2(1)=DSQRT(PL+B) PL=DSQRT(PL) COSZ=2.D0*PSRAN(B10)-1.D0 PT=PL*DSQRT(1.D0-COSZ**2) EP1(2)=PL*COSZ CALL PSCS(C,S) EP1(3)=PT*C EP1(4)=PT*S do 1 I=2,4 1 EP2(I)=-EP1(I) CALL PSDEFTR(WW,EP,EY) CALL PSTRANS(EP1,EY) CALL PSTRANS(EP2,EY) IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXDEC2 - END') RETURN END C======================================================================= SUBROUTINE XXDEC3(EP,EP1,EP2,EP3,SWW,AM1,AM2,AM3) c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION EP(4),EP1(4),EP2(4),EP3(4),EPT(4),EY(3) COMMON/AREA11/B10 COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG IF(DEBUG.GE.2)WRITE (MONIOU,201) 201 FORMAT(2X,'XXDEC3 - THREE PARTICLE DECAY') AM12=AM1**2 AM23=(AM2+AM3)**2 AM32=(AM2-AM3)**2 S23MAX=(SWW-AM1)**2 EMAX=.25D0*(SWW+(AM12-AM23)/SWW)**2 GB0=DSQRT((EMAX-AM12)/EMAX*(1.D0-AM23/S23MAX) * *(1.D0-AM32/S23MAX)) 1 P1=PSRAN(B10)*(EMAX-AM12) E1=DSQRT(P1+AM12) S23=SWW**2+AM12-2.D0*E1*SWW GB=DSQRT(P1*(1.D0-AM23/S23)*(1.D0-AM32/S23))/E1/GB0 IF(PSRAN(B10).GT.GB)GOTO 1 P1=DSQRT(P1) EP1(1)=E1 COSZ=2.D0*PSRAN(B10)-1.D0 PT=P1*DSQRT(1.D0-COSZ**2) EP1(2)=P1*COSZ CALL PSCS(C,S) EP1(3)=PT*C EP1(4)=PT*S do 2 I=2,4 2 EPT(I)=-EP1(I) EPT(1)=SWW-EP1(1) CALL PSDEFTR(SWW**2,EP,EY) CALL PSTRANS(EP1,EY) CALL PSTRANS(EPT,EY) CALL XXDEC2(EPT,EP2,EP3,S23,AM2**2,AM3**2) IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXDEC3 - END') RETURN END C======================================================================= SUBROUTINE XXDPR(WP0,WM0,ICP,ICT,LQ2) c Projectile hadron dissociation c Leading hadronic state hadronization c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION EP3(4),EP1(4),EP2(4),EY(3) COMMON /AREA1/ IA(2),ICZ,ICP0 COMMON /AREA2/ S,Y0,WP00,WM00 COMMON /AREA8/ WWM,BE(4),DC(5),DETA,ALMPT COMMON /AREA10/ STMASS,AM(7) COMMON /AREA11/ B10 COMMON /AREA17/ DEL,RS,RS0,FS,ALFP,RR,SH,DELH COMMON /AREA21/ DMMIN(5) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG SAVE IF(DEBUG.GE.2)WRITE (MONIOU,201)ICP,ICT,WP0,WM0 201 FORMAT(2X,'XXDPR - LEADING (PROJECTILE) CLUSTER HADRONIZATION:' * /4X,'CLUSTER TYPE ICP=',I2,2X,'TARGET TYPE ', * 'ICT=',I2/4X,'AVAILABLE LIGHT CONE MOMENTA: WP0=',E10.3, * ' WM0=',E10.3) DO 100 I=1,3 100 EY(I)=1.D0 SD0=WP0*WM0 IF(SD0.LT.0.D0)SD0=0.D0 DDMAX=MIN(5.D0,DSQRT(SD0)-AM(2)) DDMIN=DMMIN(ICZ) IF(DDMAX.LT.DDMIN)THEN c Registration of too slow "leading" hadron if its energy is insufficient for c diffractive exhitation EP3(3)=0.D0 EP3(4)=0.D0 IF(LQ2.NE.0)THEN WPI=WP0 IF(AM(ICZ)**2.GT.WPI*WM0)THEN IF(WPI.GT.0.D0.AND.AM(ICZ)**2/WPI.LT..5D0*WM00)THEN WMI=AM(ICZ)**2/WPI WM0=WMI ELSE RETURN ENDIF ENDIF WM0=WM0-WMI WP0=0.D0 EP3(1)=.5D0*(WPI+WMI) EP3(2)=.5D0*(WPI-WMI) CALL XXREG(EP3,ICP) IF(DEBUG.GE.3)WRITE (MONIOU,202) RETURN ELSE IF(DSQRT(SD0).LT.AM(ICZ)+AM(2))THEN IF(WP0.GT.0.D0.AND.(AM(ICZ)+AM(2))**2/WP0.LT..5D0*WM00) * THEN SD0=(AM(ICZ)+AM(2))**2 WM0=SD0/WP0 ELSE IF(DEBUG.GE.3)WRITE (MONIOU,202) RETURN ENDIF ENDIF XW=XXTWDEC(SD0,AM(ICZ)**2,AM(2)**2) WP1=XW*WP0 WM1=AM(ICZ)**2/WP1 EP3(1)=.5D0*(WP1+WM1) EP3(2)=.5D0*(WP1-WM1) CALL XXREG(EP3,ICP) WM2=WM0-WM1 WP2=AM(2)**2/WM2 EP3(1)=.5D0*(WP2+WM2) EP3(2)=.5D0*(WP2-WM2) CALL XXREG(EP3,ICT) WP0=0.D0 WM0=0.D0 ENDIF IF(DEBUG.GE.3)WRITE (MONIOU,202) RETURN ENDIF IF(ICP.NE.0)IS=IABS(ICP)/ICP DMASS=DDMIN**2/(1.D0-PSRAN(B10)*(1.D0-(DDMIN/DDMAX)))**2 IF(LQ2.NE.0)THEN WPD=WP0 WMD=DMASS/WPD WM0=WM0-WMD WP0=0.D0 ELSE IF(ICZ.EQ.5)THEN WPD=WP0*XXTWDEC(SD0,DMASS,AM(2)**2) WMD=DMASS/WPD WM2=WM0-WMD WP2=AM(2)**2/WM2 EP3(1)=.5D0*(WP2+WM2) EP3(2)=.5D0*(WP2-WM2) EP3(3)=0.D0 EP3(4)=0.D0 CALL XXREG(EP3,ICT) ELSE PTMAX=PSLAM(SD0,DMASS,AM(2)**2) IF(PTMAX.LT.0.)PTMAX=0. PTI=-1.D0/RS*DLOG(1.D0-PSRAN(B10)*(1.D0-EXP(-RS*PTMAX))) AMT1=DMASS+PTI AMT2=AM(2)**2+PTI WPD=WP0*XXTWDEC(SD0,AMT1,AMT2) WMD=AMT1/WPD WM2=WM0-WMD WP2=AMT2/WM2 PT=DSQRT(PTI) CALL PSCS(CCOS,SSIN) EP3(3)=PT*CCOS EP3(4)=PT*SSIN EP3(1)=.5D0*(WP2+WM2) EP3(2)=.5D0*(WP2-WM2) CALL XXREG(EP3,ICT) EP3(3)=-EP3(3) EP3(4)=-EP3(4) EP3(1)=.5D0*(WPD+WMD) EP3(2)=.5D0*(WPD-WMD) CALL PSDEFTR(DMASS,EP3,EY) WPD=DSQRT(DMASS) WMD=WPD ENDIF WP0=0.D0 WM0=0.D0 ENDIF IF(ICZ.EQ.5)THEN ICH1=ICP ICH2=0 AMH1=AM(5)**2 AMH2=AM(1)**2 PTMAX=PSLAM(DMASS,AMH1,AMH2) IF(PTMAX.LT.0.)PTMAX=0. IF(PTMAX.LT.BE(4)**2)THEN 1 PTI=PTMAX*PSRAN(B10) IF(PSRAN(B10).GT.EXP(-DSQRT(PTI)/BE(4)))GOTO 1 ELSE 2 PTI=(BE(4)*DLOG(PSRAN(B10)*PSRAN(B10)))**2 IF(PTI.GT.PTMAX)GOTO 2 ENDIF AMT1=AMH1+PTI AMT2=AMH2+PTI Z=XXTWDEC(DMASS,AMT1,AMT2) WP1=WPD*Z WM1=AMT1/WP1 EP3(1)=.5D0*(WP1+WM1) EP3(2)=.5D0*(WP1-WM1) PT=DSQRT(PTI) CALL PSCS(C,S) EP3(3)=PT*C EP3(4)=PT*S CALL XXREG(EP3,ICH1) WP1=WPD*(1.D0-Z) WM1=AMT2/WP1 EP3(1)=.5D0*(WP1+WM1) EP3(2)=.5D0*(WP1-WM1) EP3(3)=-PT*C EP3(4)=-PT*S CALL XXREG(EP3,ICH2) IF(DEBUG.GE.3)WRITE (MONIOU,202) RETURN ENDIF IF(ICZ.EQ.1)THEN IF(ICP.NE.0)THEN IC1=ICP*(1-3*INT(.5D0+PSRAN(B10))) IC2=-ICP-IC1 ELSE IC1=INT(1.5D0+PSRAN(B10))*(2*INT(.5D0+PSRAN(B10))-1) IC2=-IC1 ENDIF ELSEIF(ICZ.EQ.2)THEN IF(PSRAN(B10).GT..33333D0)THEN IC1=3*IS IC2=ICP-IS ELSE IC1=ICP+4*IS IC2=4*IS-ICP ENDIF ELSEIF(ICZ.EQ.3)THEN IC1=-4*IS IC2=ICP-3*IS ELSEIF(ICZ.EQ.4)THEN IC1=5*IS IC2=ICP-9*IS ENDIF CALL XXGENER(WPD,WMD,EY,0.D0,1.D0,0.D0,1.D0, * IC1,IC2) IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXDPR - END') RETURN END C======================================================================= SUBROUTINE XXDTG(WP0,WM0,ICP,ICT,LQ1) c Target nucleon dissociation c Leading hadronic state hadronization c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION EP3(4),EY(3) COMMON /AREA1/ IA(2),ICZ,ICP0 COMMON /AREA2/ S,Y0,WP00,WM00 COMMON /AREA10/ STMASS,AM(7) COMMON /AREA11/ B10 COMMON /AREA17/ DEL,RS,RS0,FS,ALFP,RR,SH,DELH COMMON /AREA21/ DMMIN(5) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG SAVE IF(DEBUG.GE.2)WRITE (MONIOU,201)ICT,ICP,WP0,WM0 201 FORMAT(2X,'XXDTG - LEADING (TARGET) CLUSTER HADRONIZATION:' * /4X,'CLUSTER TYPE ICT=',I2,2X,'PROJECTILE TYPE ', * 'ICP=',I2/4X,'AVAILABLE LIGHT CONE MOMENTA: WP0=',E10.3, * ' WM0=',E10.3) DO 100 I=1,3 100 EY(I)=1.D0 SD0=WP0*WM0 IF(SD0.LT.0.D0)SD0=0.D0 DDMIN=DMMIN(2) DDMAX=MIN(5.D0,DSQRT(SD0)-AM(ICZ)) IF(DDMAX.LT.DDMIN)THEN c Registration of too slow "leading" hadron if its energy is insufficient for c diffractive exhitation EP3(3)=0.D0 EP3(4)=0.D0 IF(LQ1.NE.0)THEN WMI=WM0 IF( WP0.LE.0.D0.OR.AM(2)**2.GT.WMI*WP0)RETURN WPI=AM(2)**2/WMI WP0=WP0-WPI WM0=0.D0 EP3(1)=.5D0*(WPI+WMI) EP3(2)=.5D0*(WPI-WMI) CALL XXREG(EP3,ICT) IF(DEBUG.GE.3)WRITE (MONIOU,202) RETURN ELSE IF(DSQRT(SD0).LT.AM(ICZ)+AM(2))THEN IF(WP0.GT.0.D0.AND.(AM(ICZ)+AM(2))**2/WP0.LT..5D0*WM00) * THEN SD0=(AM(ICZ)+AM(2))**2 WM0=SD0/WP0 ELSE IF(DEBUG.GE.3)WRITE (MONIOU,202) RETURN ENDIF ENDIF XW=XXTWDEC(SD0,AM(ICZ)**2,AM(2)**2) WP1=XW*WP0 WM1=AM(ICZ)**2/WP1 EP3(1)=.5D0*(WP1+WM1) EP3(2)=.5D0*(WP1-WM1) CALL XXREG(EP3,ICP) WM2=WM0-WM1 WP2=AM(2)**2/WM2 EP3(1)=.5D0*(WP2+WM2) EP3(2)=.5D0*(WP2-WM2) CALL XXREG(EP3,ICT) WP0=0.D0 WM0=0.D0 ENDIF IF(DEBUG.GE.3)WRITE (MONIOU,202) RETURN ENDIF DMASS=(DDMIN/(1.D0-PSRAN(B10)*(1.D0-DDMIN/DDMAX)))**2 IF(LQ1.NE.0)THEN WMD=WM0 WPD=DMASS/WMD WP0=WP0-WPD WM0=0.D0 ELSE PTMAX=PSLAM(SD0,DMASS,AM(ICZ)**2) IF(PTMAX.LT.0.)PTMAX=0. PTI=-1.D0/RS*DLOG(1.D0-PSRAN(B10)*(1.D0-EXP(-RS*PTMAX))) AMT1=DMASS+PTI AMT2=AM(ICZ)**2+PTI WMD=WM0*XXTWDEC(SD0,AMT1,AMT2) WPD=AMT1/WMD WP2=WP0-WPD WM2=AMT2/WP2 PT=DSQRT(PTI) CALL PSCS(CCOS,SSIN) EP3(3)=PT*CCOS EP3(4)=PT*SSIN EP3(1)=.5D0*(WP2+WM2) EP3(2)=.5D0*(WP2-WM2) CALL XXREG(EP3,ICP) EP3(3)=-EP3(3) EP3(4)=-EP3(4) EP3(1)=.5D0*(WPD+WMD) EP3(2)=.5D0*(WPD-WMD) CALL PSDEFTR(DMASS,EP3,EY) WPD=DSQRT(DMASS) WMD=WPD WP0=0.D0 WM0=0.D0 ENDIF IS=IABS(ICT)/ICT IF(PSRAN(B10).GT..33333D0)THEN IC1=3*IS IC2=ICT-IS ELSE IC1=ICT+4*IS IC2=4*IS-ICT ENDIF CALL XXGENER(WPD,WMD,EY, * 0.D0,1.D0,0.D0,1.D0,IC2,IC1) IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXDTG - END') RETURN END C======================================================================= SUBROUTINE XXFAU(B,GZ) c Integrands for hadron-hadron and hadron-nucleus cross-sections calculation c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION GZ(3),GZ0(2) COMMON /AREA1/ IA(2),ICZ,ICP COMMON /AREA16/ CC(5) COMMON /AR1/ ANORM COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG IF(DEBUG.GE.2)WRITE (MONIOU,201) 201 FORMAT(2X,'XXFAU - INTEGRANDS FOR HADRON-HADRON AND ' * 'HADRON-NUCLEUS CROSS-SECTIONS CALCULATION') CALL XXFZ(B,GZ0) DO 1 L=1,2 1 GZ0(L)=GZ0(L)*CC(2)*ANORM*.5D0 AB=FLOAT(IA(2)) GZ1=(1.D0-GZ0(1))**AB GZ2=(1.D0-GZ0(2))**AB GZ3=(1.D0-CC(2)*GZ0(2)-2.D0*(1.D0-CC(2))*GZ0(1))**AB GZ(1)=CC(ICZ)**2*(GZ2-GZ3) GZ(2)=CC(ICZ)*(1.D0-CC(ICZ))*(1.D0+GZ2-2.D0*GZ1) GZ(3)=CC(ICZ)*(1.D0-GZ2) IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXFAU - END') RETURN END C======================================================================= SUBROUTINE XXFRAG(SA,NA,RC) c Connected nucleon clasters extraction - used for the nuclear spectator part c multifragmentation: c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION SA(56,3) COMMON /AREA13/ NSF,IAF(56) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG SAVE IF(DEBUG.GE.2)WRITE (MONIOU,201)NA 201 FORMAT(2X,'XXFRAG-MULTIFRAGMENTATION: NUCLEUS MASS NUMBER: NA=' * ,I2) IF(DEBUG.GE.3)THEN WRITE (MONIOU,203) 203 FORMAT(2X,'NUCLEONS COORDINATES:') 204 FORMAT(2X,3E10.3) DO 205 I=1,NA 205 WRITE (MONIOU,204)(SA(I,L),L=1,3) ENDIF NI=1 NG=1 J=0 1 J=J+1 J1=NI+1 DO 4 I=J1,NA RI=0.D0 DO 2 M=1,3 2 RI=RI+(SA(J,M)-SA(I,M))**2 IF(RI.GT.RC)GOTO 4 NI=NI+1 NG=NG+1 IF(I.EQ.NI)GOTO 4 DO 3 M=1,3 S0=SA(NI,M) SA(NI,M)=SA(I,M) 3 SA(I,M)=S0 4 CONTINUE IF(J.LT.NI.AND.NA-NI.GT.0)GOTO 1 NSF=NSF+1 IAF(NSF)=NG IF(DEBUG.GE.3)WRITE (MONIOU,206)NSF,IAF(NSF) 206 FORMAT(2X,'XXFRAG: FRAGMENT N',I2,2X,'FRAGMENT MASS - ',I2) NG=1 J=NI NI=NI+1 IF(NA-NI)6,5,1 5 NSF=NSF+1 IAF(NSF)=1 IF(DEBUG.GE.3)WRITE (MONIOU,206)NSF,IAF(NSF) 6 CONTINUE IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXFRAG - END') RETURN END C======================================================================= SUBROUTINE XXFRAGM(NS,XA) c Fragmentation of the spectator part of the nucleus c XA(56,3) - arrays for spectator nucleons positions c NS - total number of spectators c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION XA(56,3) INTEGER DEBUG COMMON /AREA1/ IA(2),ICZ,ICP COMMON /AREA3/ RMIN,EMAX,EEV COMMON /AREA11/ B10 c NSF - number of secondary fragments; c IAF(i) - mass of the i-th fragment COMMON /AREA13/ NSF,IAF(56) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG SAVE IF(DEBUG.GE.2)WRITE (MONIOU,201)NS 201 FORMAT(2X,'XXFRAGM: NUMBER OF SPECTATORS: NS=',I2) NSF=0 IF(NS-1)6,1,2 c Single spectator nucleon is recorded 1 NSF=NSF+1 IAF(NSF)=1 IF(DEBUG.GE.3)WRITE (MONIOU,205) 205 FORMAT(2X,'XXFRAGM - SINGLE SPECTATOR') GOTO 6 2 EEX=0.D0 c EEX - spectator part excitation energy; calculated as the sum of excitations c from all wounded nucleons ( including diffractively excited ) DO 3 I=1,IA(1)-NS c Partial excitation is simulated according to distribution f(E) ~ 1/sqrt(E) c * exp(-E/(2*)), for sqrt(E) we have then normal distribution 3 EEX=EEX+(PSRAN(B10)+PSRAN(B10)+PSRAN(B10)+ * PSRAN(B10)+PSRAN(B10)-2.5D0)**2*2.4D0 IF(DEBUG.GE.3)WRITE (MONIOU,203)EEX 203 FORMAT(2X,'XXFRAGM: EXCITATION ENERGY: EEX=',E10.3) c If the excitation energy per spectator is larger than EMAX c multifragmentation takes place ( percolation algorithm is used for it ) IF(EEX/NS.GT.EMAX)THEN c Multifragmentation CALL XXFRAG(XA,NS,RMIN) ELSE c Otherwise average number of eveporated nucleons equals EEX/EEV, where c EEV - mean excitation energy carried out by one nucleon NF=IXXSON(NS,EEX/EEV,PSRAN(B10)) NSF=NSF+1 c Recording of the fragment produced IAF(NSF)=NS-NF IF(DEBUG.GE.3)WRITE (MONIOU,206)IAF(NSF) 206 FORMAT(2X,'XXFRAGM - EVAPORATION: MASS NUMBER OF THE FRAGMENT:' * ,I2) c Some part of excitation energy is carried out by alphas; we determine the c number of alphas simply as NF/4 NAL=NF/4 IF(NAL.NE.0)THEN c Recording of the evaporated alphas DO 4 I=1,NAL NSF=NSF+1 4 IAF(NSF)=4 ENDIF NF=NF-4*NAL IF(NF.NE.0)THEN c Recording of the evaporated nucleons DO 5 I=1,NF NSF=NSF+1 5 IAF(NSF)=1 ENDIF IF(DEBUG.GE.3)WRITE (MONIOU,204)NF,NAL 204 FORMAT(2X,'XXFRAGM - EVAPORATION: NUMBER OF NUCLEONS NF=',I2, * 'NUMBER OF ALPHAS NAL=',I2) ENDIF 6 CONTINUE IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXFRAGM - END') RETURN END C======================================================================= SUBROUTINE XXFZ(B,GZ) c Hadron-hadron and hadron-nucleus cross sections calculation c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION GZ(2),FHARD(3) COMMON /AREA1/ IA(2),ICZ,ICP COMMON /AREA2/ S,Y0,WP0,WM0 COMMON /AREA7/ RP1 COMMON /AR3/ X1(7),A1(7) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG IF(DEBUG.GE.2)WRITE (MONIOU,201) 201 FORMAT(2X,'XXFZ - HADRONIC CROSS-SECTIONS CALCULATION') DO 1 L=1,2 1 GZ(L)=0.D0 E1=EXP(-1.D0) DO 2 I1=1,7 DO 2 M=1,2 Z=.5D0+X1(I1)*(M-1.5D0) S1=DSQRT(RP1*Z) ZV1=EXP(-Z) S2=DSQRT(RP1*(1.D0-DLOG(Z))) ZV2=E1*Z C?????????? C VV1=EXP(-PSFAZ(ZV1,FSOFT,FHARD,FSHARD))*(1.D0-FHARD(1) C * -FHARD(2)-FHARD(3)) C VV2=EXP(-PSFAZ(ZV2,FSOFT,FHARD,FSHARD))*(1.D0-FHARD(1) C * -FHARD(2)-FHARD(3)) VV1=EXP(-PSFAZ(ZV1,FSOFT,FHARD,FSHARD)-FHARD(1) * -FHARD(2)-FHARD(3)) VV2=EXP(-PSFAZ(ZV2,FSOFT,FHARD,FSHARD)-FHARD(1) * -FHARD(2)-FHARD(3)) c??????????? IF(IA(2).EQ.1)THEN CG1=1.D0 CG2=1.D0 ELSE CG1=XXROT(B,S1) CG2=XXROT(B,S2) ENDIF DO 2 L=1,2 2 GZ(L)=GZ(L)+ A1(I1)*(CG1*(1.D0-VV1**L)+CG2*(1.D0-VV2**L)/Z) IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXFZ - END') RETURN END C======================================================================= SUBROUTINE XXGAU(GZ) c Impact parameter integration for impact parameters BM - c for hadron-hadron and hadron-nucleus cross-sections calculation c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION GZ(3),GZ0(3) COMMON /AREA6/ PI,BM,AM COMMON /AR5/ X5(2),A5(2) COMMON /AR2/ R,RM COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG IF(DEBUG.GE.2)WRITE (MONIOU,201) 201 FORMAT(2X,'XXGAU1 - NUCLEAR CROSS-SECTIONS CALCULATION') DO 1 I=1,2 B=BM+X5(I) CALL XXFAU(B,GZ0) DO 1 L=1,3 1 GZ(L)=GZ(L)+GZ0(L)*A5(I)*EXP(X5(I))*B*2.D0*PI*AM*AM IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXGAU1 - END') RETURN END C======================================================================= SUBROUTINE XXGENER(WP0,WM0,EY0,S0X,C0X,S0,C0,IC1,IC2) c To simulate the fragmentation of the string into secondary hadrons c The algorithm conserves energy-momentum; c WP0, WM0 are initial longitudinal momenta ( E+p, E-p ) of the quarks c at the ends of the string; IC1, IC2 - their types c The following partons types are used: 1 - u, -1 - U, 2 - d, -2 - D, c 3 - ud, -3 - UD, 4 - s, -4 - S, 5 - c, -5 - C, c 6 - uu, 7 - dd, -6 - UU, -7 - DD c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG CHARACTER *2 TYQ DIMENSION WP(2),IC(2),EPT(4),EP(4),EY(3),EY0(3) c WP(1), WP(2) - current longitudinal momenta of the partons at the string c ends, IC(1), IC(2) - their types COMMON /AREA8/ WWM,BEP,BEN,BEK,BEC,DC(5),DETA,ALMPT COMMON /AREA10/ STMASS,AM0,AMN,AMK,AMC,AMLAMC,AMLAM,AMETA COMMON /AREA11/ B10 COMMON /AREA19/ AHL(5) ******************************************************** COMMON /AREA21/ DMMIN(5) ******************************************************** COMMON /AREA28/ ARR(4) COMMON /AREA42/ TYQ(15) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG SAVE IF(DEBUG.GE.2)WRITE (MONIOU,201)TYQ(8+IC1),TYQ(8+IC2), * WP0,WM0,EY0,S0X,C0X,S0,C0 201 FORMAT(2X,'XXGENER: PARTON FLAVORS AT THE ENDS OF THE STRING:', * 2X,A2,2X,A2/4X,'LIGHT CONE MOMENTA OF THE STRING: ',E10.3, * 2X,E10.3/4X,'EY0=',3E10.3/4X, * 'S0X=',E10.3,2X,'C0X=',E10.3,2X,'S0=',E10.3,2X,'C0=',E10.3) WW=WP0*WM0 EPT(1)=.5D0*(WP0+WM0) EPT(2)=.5D0*(WP0-WM0) EPT(3)=0.D0 EPT(4)=0.D0 IC(1)=IC1 IC(2)=IC2 1 SWW=DSQRT(WW) CALL PSDEFTR(WW,EPT,EY) J=INT(2.D0*PSRAN(B10))+1 IF(DEBUG.GE.3)THEN IQT=8+IC(J) WRITE (MONIOU,203)J,TYQ(IQT),WW 203 FORMAT(2X,'XXGENER: CURRENT PARTON FLAVOR AT THE END ',I1, * ' OF THE STRING: ',A2/4X,' STRING MASS: ',E10.3) ENDIF IAB=IABS(IC(J)) IS=IC(J)/IAB IF(IAB.GT.5)IAB=3 IAJ=IABS(IC(3-J)) IF(IAJ.GT.5)IAJ=3 IF(IAJ.EQ.3)THEN RESTM=AMN ELSEIF(IAJ.EQ.4)THEN RESTM=AMK ELSEIF(IAJ.EQ.5)THEN RESTM=AMC ELSE RESTM=AM0 ENDIF IF(IAB.LE.2.AND.SWW.GT.RESTM+2.D0*AM0+WWM.OR. *IAB.EQ.3.AND.SWW.GT.RESTM+AM0+AMN+WWM.OR. *IAB.EQ.4.AND.SWW.GT.RESTM+AM0+AMK+WWM.OR. *IAB.EQ.5.AND.SWW.GT.RESTM+AM0+AMC+WWM)THEN IF(IAB.LE.2)THEN IF(SWW.GT.RESTM+2.D0*AMC.AND.PSRAN(B10).LT.DC(3))THEN c D-meson generation RESTM=(RESTM+AMC)**2 BET=BEC AMI=AMC**2 ALF=ALMPT-ARR(4) BLF=AHL(4) IC0=IC(J)-9*IS IC(J)=5*IS ELSEIF(SWW.GT.RESTM+2.D0*AMN.AND.PSRAN(B10).LT.DC(1))THEN c Nucleon generation RESTM=(RESTM+AMN)**2 BET=BEN AMI=AMN**2 ALF=ALMPT-ARR(2) BLF=AHL(2) IC0=IC(J)+IS IC(J)=-3*IS ELSEIF(SWW.GT.RESTM+2.D0*AMK.AND.PSRAN(B10).LT.DC(2))THEN c Kaon generation RESTM=(RESTM+AMK)**2 BET=BEK AMI=AMK**2 ALF=ALMPT-ARR(3) BLF=AHL(3) IC0=IC(J)+3*IS IC(J)=4*IS ELSEIF(SWW.GT.RESTM+AMETA+AM0.AND.PSRAN(B10).LT.DETA)THEN c Eta generation RESTM=(RESTM+AM0)**2 BET=BEK AMI=AMETA**2 ALF=ALMPT-ARR(1) BLF=AHL(1) IC0=10 ELSE c Pion generation RESTM=(RESTM+AM0)**2 BET=BEP AMI=AM0**2 ALF=ALMPT-ARR(1) BLF=AHL(1) IF(PSRAN(B10).LT..3333D0)THEN IC0=0 ELSE IC0=3*IS-2*IC(J) IC(J)=3*IS-IC(J) ENDIF ENDIF ELSEIF(IAB.EQ.3)THEN IF(SWW.GT.RESTM+AMC+AMLAMC.AND.PSRAN(B10).LT.DC(5).AND. * IABS(IC(J)).EQ.3)THEN c Lambda_C generation RESTM=(RESTM+AMC)**2 BET=BEC AMI=AMLAMC**2 ALF=ALMPT-ARR(4) BLF=AHL(5) IC0=9*IS IC(J)=-5*IS ELSEIF(SWW.GT.RESTM+AMK+AMLAM.AND.PSRAN(B10).LT.DC(4).AND. * IABS(IC(J)).EQ.3)THEN c Lambda generation RESTM=(RESTM+AMK)**2 BET=BEK AMI=AMLAM**2 ALF=ALMPT-ARR(3) BLF=AHL(2)+ARR(1)-ARR(3) IC0=6*IS IC(J)=-4*IS ELSE c Nucleon generation RESTM=(RESTM+AM0)**2 BET=BEN AMI=AMN**2 ALF=ALMPT-ARR(1) BLF=AHL(2) IF(IABS(IC(J)).EQ.3)THEN IC0=IS*INT(2.5D0+PSRAN(B10)) IC(J)=IS-IC0 ELSE IC0=IC(J)-4*IS IC(J)=IC0-4*IS ENDIF ENDIF ELSEIF(IAB.EQ.4)THEN IF(SWW.GT.RESTM+AMN+AMLAM.AND.PSRAN(B10).LT.DC(1))THEN c Lambda generation RESTM=(RESTM+AMN)**2 BET=BEN AMI=AMLAM**2 ALF=ALMPT-ARR(2) BLF=AHL(2)+ARR(1)-ARR(3) IC0=6*IS IC(J)=-3*IS ELSE c Kaon generation RESTM=(RESTM+AM0)**2 BET=BEP AMI=AMK**2 ALF=ALMPT-ARR(1) BLF=AHL(3) IC(J)=IS*INT(1.5D0+PSRAN(B10)) IC0=-3*IS-IC(J) ENDIF ELSEIF(IAB.EQ.5)THEN IF(SWW.GT.RESTM+AMN+AMLAMC.AND.PSRAN(B10).LT.DC(1))THEN c Lambda_C generation RESTM=(RESTM+AMN)**2 BET=BEN AMI=AMLAMC**2 ALF=ALMPT-ARR(2) BLF=AHL(5) IC0=9*IS IC(J)=-3*IS ELSE c D-meson generation RESTM=(RESTM+AM0)**2 BET=BEP AMI=AMC**2 ALF=ALMPT-ARR(1) BLF=AHL(4) IC(J)=IS*INT(1.5D0+PSRAN(B10)) IC0=9*IS-IC(J) ENDIF ENDIF ******************************************************** PTMAX=PSLAM(WW,RESTM,AMI) IF(PTMAX.LT.0.)PTMAX=0. IF(PTMAX.LT.BET**2)THEN 2 PTI=PTMAX*PSRAN(B10) IF(PSRAN(B10).GT.EXP(-DSQRT(PTI)/BET))GOTO 2 ELSE 3 PTI=(BET*DLOG(PSRAN(B10)*PSRAN(B10)))**2 IF(PTI.GT.PTMAX)GOTO 3 ENDIF AMT=AMI+PTI RESTM1=RESTM+PTI ******************************************************** c ALF=ALF+2.*PTI ZMIN=DSQRT(AMT/WW) ZMAX=XXTWDEC(WW,AMT,RESTM1) Z1=(1.-ZMAX)**ALF Z2=(1.-ZMIN)**ALF 4 Z=1.-(Z1+(Z2-Z1)*PSRAN(B10))**(1./ALF) IF(PSRAN(B10).GT.(Z/ZMAX)**BLF)GOTO 4 WP(J)=Z*SWW WP(3-J)=AMT/WP(J) EP(1)=.5D0*(WP(1)+WP(2)) EP(2)=.5D0*(WP(1)-WP(2)) PTI=DSQRT(PTI) CALL PSCS(C,S) EP(3)=PTI*C EP(4)=PTI*S EPT(1)=SWW-EP(1) DO 5 I=2,4 5 EPT(I)=-EP(I) WW=PSNORM(EPT) IF(WW.LT.RESTM)GOTO 4 CALL PSTRANS(EP,EY) CALL PSTRANS(EPT,EY) IF(S0X.NE.0.D0.OR.S0.NE.0.D0)THEN CALL PSROTAT(EP,S0X,C0X,S0,C0) ENDIF IF(EY0(1)*EY0(2)*EY0(3).NE.1.D0)THEN CALL PSTRANS(EP,EY0) ENDIF CALL XXREG(EP,IC0) ELSE AMI2=RESTM**2 BET=BEP IF(IAB.LE.2.AND.IAJ.LE.2)THEN AMI=AM0**2 IC0=-IC(1)-IC(2) IF(IC0.NE.0)THEN IC(J)=IC0*INT(.5D0+PSRAN(B10)) IC(3-J)=IC0-IC(J) ELSE IF(PSRAN(B10).LT..2D0)THEN IC(J)=0 IC(3-J)=0 ELSE IC(J)=3*IS-2*IC(J) IC(3-J)=-IC(J) ENDIF ENDIF ELSEIF(IAB.EQ.3.OR.IAJ.EQ.3)THEN IF(IAB.EQ.3)THEN AMI=AMN**2 IF(IABS(IC(J)).EQ.3)THEN IF(IAJ.EQ.3)THEN IF(IABS(IC(3-J)).EQ.3)THEN IC(J)=IS*INT(2.5D0+PSRAN(B10)) IC(3-J)=-IC(J) ELSE IC(3-J)=IC(3-J)+4*IS IC(J)=5*IS+IC(3-J) ENDIF ELSEIF(IAJ.LT.3)THEN IF(PSRAN(B10).LT..3333D0)THEN IC(J)=IC(3-J)+IS IC(3-J)=0 ELSE IC(J)=IS*(4-IAJ) IC(3-J)=IS*(3-2*IAJ) ENDIF ELSEIF(IAJ.EQ.4)THEN IC(J)=IS*INT(2.5D0+PSRAN(B10)) IC(3-J)=-IC(J)-2*IS ELSEIF(IAJ.EQ.5)THEN IC(J)=IS*INT(2.5D0+PSRAN(B10)) IC(3-J)=-IC(J)+10*IS ENDIF ELSE IC(J)=IC(J)-4*IS IC0=IC(J)-4*IS IF(IAJ.EQ.3)THEN IC(3-J)=IC0-IS ELSEIF(IAJ.LT.3)THEN IC(3-J)=-IC(3-J)-IC0 ELSEIF(IAJ.EQ.4)THEN IC(3-J)=IC0-3*IS ELSEIF(IAJ.EQ.5)THEN IC(3-J)=IC0+9*IS ENDIF ENDIF ELSE IF(IABS(IC(3-J)).EQ.3)THEN IF(IAB.LT.3)THEN AMI=AM0**2 IF(PSRAN(B10).LT..3333D0)THEN IC(3-J)=IC(J)+IS IC(J)=0 ELSE IC(3-J)=IS*(4-IAB) IC(J)=IS*(3-2*IAB) ENDIF ELSEIF(IAB.EQ.4)THEN AMI=AMK**2 IC(3-J)=IS*INT(2.5D0+PSRAN(B10)) IC(J)=-IC(3-J)-2*IS ELSEIF(IAB.EQ.5)THEN AMI=AMC**2 IC(3-J)=IS*INT(2.5D0+PSRAN(B10)) IC(J)=-IC(3-J)+10*IS ENDIF ELSE IC(3-J)=IC(3-J)-4*IS IC0=IC(3-J)-4*IS IF(IAB.LT.3)THEN AMI=AM0**2 IC(J)=-IC0-IC(J) ELSEIF(IAB.EQ.4)THEN AMI=AMK**2 IC(J)=IC0-3*IS ELSEIF(IAB.EQ.5)THEN AMI=AMC**2 IC(J)=IC0+9*IS ENDIF ENDIF ENDIF ELSEIF(IAB.EQ.4.OR.IAJ.EQ.4)THEN IF(IAB.EQ.4)THEN AMI=AMK**2 IF(IAJ.EQ.4)THEN IC(J)=-IS*INT(4.5D0+PSRAN(B10)) IC(3-J)=-IC(J) ELSEIF(IAJ.EQ.5)THEN IC(J)=-IS*INT(4.5D0+PSRAN(B10)) IC(3-J)=-IC(J)-12*IS ELSE IC0=IC(3-J)+INT(.6667D0+PSRAN(B10))*(-3*IS-2*IC(3-J)) IC(J)=IC0-3*IS IC(3-J)=IC0-IC(3-J) ENDIF ELSE IF(IAB.LE.2)THEN AMI=AM0**2 IC0=IC(J)+INT(.6667D0+PSRAN(B10))*(3*IS-2*IC(J)) IC(J)=IC0-IC(J) IC(3-J)=IC0+3*IS ELSEIF(IAB.EQ.5)THEN AMI=AMC**2 IC(3-J)=IS*INT(4.5D0+PSRAN(B10)) IC(J)=-IC(3-J)+12*IS ENDIF ENDIF ELSEIF(IAB.EQ.5.OR.IAJ.EQ.5)THEN IF(IAB.EQ.5)THEN AMI=AMC**2 IF(IAJ.EQ.5)THEN IC(J)=IS*INT(7.5D0+PSRAN(B10)) IC(3-J)=-IC(J) ELSE IC0=IC(3-J)+INT(.6667D0+PSRAN(B10))*(-3*IS-2*IC(3-J)) IC(J)=IC0+9*IS IC(3-J)=IC0-IC(3-J) ENDIF ELSE AMI=AM0**2 IC0=IC(J)+INT(.6667D0+PSRAN(B10))*(3*IS-2*IC(J)) IC(J)=IC0-IC(J) IC(3-J)=IC0-9*IS ENDIF ENDIF PTMAX=PSLAM(WW,AMI2,AMI) IF(PTMAX.LT.0.)PTMAX=0. IF(PTMAX.LT.BET**2)THEN 6 PTI=PTMAX*PSRAN(B10) IF(PSRAN(B10).GT.EXP(-DSQRT(PTI)/BET))GOTO 6 ELSE 7 PTI=(BET*DLOG(PSRAN(B10)*PSRAN(B10)))**2 IF(PTI.GT.PTMAX)GOTO 7 ENDIF AMT1=AMI+PTI AMT2=AMI2+PTI Z=XXTWDEC(WW,AMT1,AMT2) WP(J)=Z*SWW WP(3-J)=AMT1/WP(J) EP(1)=.5D0*(WP(1)+WP(2)) EP(2)=.5D0*(WP(1)-WP(2)) PTI=DSQRT(PTI) CALL PSCS(C,S) EP(3)=PTI*C EP(4)=PTI*S EPT(1)=SWW-EP(1) DO 8 I=2,4 8 EPT(I)=-EP(I) CALL PSTRANS(EP,EY) CALL PSTRANS(EPT,EY) IF(S0X.NE.0.D0.OR.S0.NE.0.D0)THEN CALL PSROTAT(EP,S0X,C0X,S0,C0) CALL PSROTAT(EPT,S0X,C0X,S0,C0) ENDIF IF(EY0(1)*EY0(2)*EY0(3).NE.1.D0)THEN CALL PSTRANS(EP,EY0) CALL PSTRANS(EPT,EY0) ENDIF CALL XXREG(EP,IC(J)) CALL XXREG(EPT,IC(3-J)) IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXGENER - END') RETURN ENDIF GOTO 1 END C======================================================================= SUBROUTINE XXJETSIM c Procedure for jet hadronization - each gluon is c considered to be splitted into quark-antiquark pair and usual soft c strings are assumed to be formed between quark and antiquark c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION EP(4),EP1(4),ey(3) COMMON /AREA10/ STMASS,AM(7) COMMON /AREA11/ B10 COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG COMMON /AREA46/ EPJET(4,2,1000),IPJET(2,1000) COMMON /AREA47/ NJTOT IF(DEBUG.GE.2)WRITE (MONIOU,201)NJTOT 201 FORMAT(2X,'XXJETSIM: TOTAL NUMBER OF JETS NJTOT=',I4) IF(NJTOT.EQ.0)RETURN DO 2 NJ=1,NJTOT DO 1 I=1,4 EP1(I)=EPJET(I,1,NJ) 1 EP(I)=EP1(I)+EPJET(I,2,NJ) PT3=DSQRT(EP1(3)**2+EP1(4)**2) PT4=DSQRT(EPJET(3,2,NJ)**2+EPJET(4,2,NJ)**2) c Invariant mass square for the jet WW=PSNORM(EP) SWW=DSQRT(WW) CALL PSDEFTR(WW,EP,EY) CALL PSTRANS1(EP1,EY) CALL PSDEFROT(EP1,S0X,C0X,S0,C0) 2 CALL XXGENER(SWW,SWW,EY,S0X,C0X,S0,C0,IPJET(1,NJ),IPJET(2,NJ)) IF(DEBUG.GE.3)WRITE (MONIOU,202) 202 FORMAT(2X,'XXJETSIM - END') RETURN END C======================================================================= SUBROUTINE XXREG(EP0,IC) c Registration of the produced hadron; c EP - 4-momentum, c IC - hadron type c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION EP(4),EP0(4) COMMON /AREA4/ EY0(3) COMMON /AREA10/ STMASS,AM0,AMN,AMK,AMC,AMLAMC,AMLAM,AMETA COMMON /AREA11/ B10 COMMON /AREA12/ NSH COMMON /AREA14/ ESP(4,15000),ICH(15000) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG IF(DEBUG.GE.2)WRITE (MONIOU,201)IC,EP0 201 FORMAT(2X,'XXREG: IC=',I2,2X,'C.M. 4-MOMENTUM:',2X,4(E10.3,1X)) pt=dsqrt(ep0(3)**2+ep0(4)**2) c if(pt.gt.11.d0)write (MONIOU,*)'pt,ic,ep',pt,ic,ep0 c if(pt.gt.11.d0)write (*,*)'pt,ic,ep',pt,ic,ep0 NSH=NSH+1 IF (NSH .GT. 15000) THEN WRITE(MONIOU,*)'XXREG: TOO MUCH SECONDARY PARTICLES' WRITE(MONIOU,*)'XXREG: NSH = ',NSH STOP ENDIF DO 4 I=1,4 4 EP(I)=EP0(I) CALL PSTRANS(EP,EY0) IF(DEBUG.GE.3)WRITE (MONIOU,202)EP 202 FORMAT(2X,'XXREG: LAB. 4-MOMENTUM:',2X,4(E10.3,1X)) ICH(NSH)=IC DO 3 I=1,4 3 ESP(I,NSH)=EP(I) IF(DEBUG.GE.3)WRITE (MONIOU,203) 203 FORMAT(2X,'XXREG - END') RETURN END C======================================================================= FUNCTION XXROT(S,B) c Convolution of nuclear profile functions (axial angle integration) c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG COMMON /AR8/ X2(4),A2 COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG IF(DEBUG.GE.2)WRITE (MONIOU,201)B 201 FORMAT(2X,'XXROT - AXIAL ANGLE INTEGRATION OF THE ', * 'NUCLEAR PROFILE FUNCTION'/4X, * 'IMPACT PARAMETER B=',E10.3,2X,'NUCLEON COORDINATE S=',E10.3) XXROT=0. DO 1 I=1,4 SB1=B**2+S**2-2.*B*S*(2.*X2(I)-1.) SB2=B**2+S**2-2.*B*S*(1.-2.*X2(I)) 1 XXROT=XXROT+(XXT(SB1)+XXT(SB2)) XXROT=XXROT*A2 IF(DEBUG.GE.3)WRITE (MONIOU,202)XXROT 202 FORMAT(2X,'XXROT=',E10.3) RETURN END C======================================================================= SUBROUTINE XXSTR(WPI0,WMI0,WP0,WM0,IC10,IC120,IC210,IC20) ************************************************** c Fragmentation process for the pomeron ( quarks and antiquarks types at the c ends of the two strings are determined, energy-momentum is shared c between them and strings fragmentation is simulated ) c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG DIMENSION EY(3) COMMON /AREA6/ PI,BM,AMMM COMMON /AREA10/ STMASS,AM(7) COMMON /AREA11/ B10 COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG SAVE IF(DEBUG.GE.2)WRITE (MONIOU,201)WPI0,WMI0,WP0,WM0 201 FORMAT(2X,'XXSTR: WPI0=',E10.3,2X,'WMI0=',E10.3,2X, * 'WP0=',E10.3,2X,'WM0=',E10.3) DO 1 I=1,3 1 EY(I)=1.D0 WPI=WPI0 WMI=WMI0 c Quark-antiquark types (1 - u, 2 - d, -1 - u~, -2 - d~); s- and d- quarks are c taken into consideration at the fragmentation step ************************************************** IF(IC10.EQ.0)THEN IC1=INT(1.5+PSRAN(B10)) IC12=-IC1 ELSEIF(IC10.GT.0)THEN IC1=IC10 IC12=IC120 ELSE IC1=IC120 IC12=IC10 ENDIF IF(IC20.EQ.0)THEN IC2=INT(1.5+PSRAN(B10)) IC21=-IC2 ELSEIF(IC20.gt.0)THEN IC2=IC20 IC21=IC210 ELSE IC2=IC210 IC21=IC20 ENDIF ************************************************** c Longitudinal momenta for the strings WP1=WPI*COS(PI*PSRAN(B10))**2 WM1=WMI*COS(PI*PSRAN(B10))**2 WPI=WPI-WP1 WMI=WMI-WM1 c String masses SM1=WP1*WM1 SM2=WPI*WMI c Too short strings are neglected (energy is given to partner string or to the hadron c (nucleon) to which the pomeron is connected) IF(SM1.GT.STMASS.AND.SM2.GT.STMASS)THEN c Strings fragmentation is simulated - GENER CALL XXGENER(WP1,WM1,EY,0.D0,1.D0,0.D0,1.D0,IC1,IC21) CALL XXGENER(WPI,WMI,EY,0.D0,1.D0,0.D0,1.D0,IC12,IC2) ELSEIF(SM1.GT.STMASS)THEN CALL XXGENER(WP1+WPI,WM1+WMI,EY,0.D0,1.D0,0.D0,1.D0,IC1,IC21) ELSEIF(SM2.GT.STMASS)THEN CALL XXGENER(WPI+WP1,WMI+WM1,EY,0.D0,1.D0,0.D0,1.D0,IC12,IC2) ELSE WP0=WP0+WP1+WPI WM0=WM0+WM1+WMI ENDIF IF(DEBUG.GE.3)WRITE (MONIOU,202)WP0,WM0 202 FORMAT(2X,'XXSTR - RETURNED LIGHT CONE MOMENTA:', * 2X,'WP0=',E10.3,2X,'WM0=',E10.3) RETURN END C======================================================================= FUNCTION XXT(B) c Nuclear profile function value at impact parameter squared B c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG COMMON /AREA6/ PI,BM,AM COMMON /AR2/ R,RM COMMON /AR5/ X5(2),A5(2) COMMON /AR9/ X9(3),A9(3) COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG IF(DEBUG.GE.2)WRITE (MONIOU,201)B 201 FORMAT(2X,'XXT - NUCLEAR PROFILE FUNCTION VALUE AT IMPACT', * ' PARAMETER SQUARED B=',E10.3) XXT=0. ZM=RM**2-B IF(ZM.GT.4.*B)THEN ZM=DSQRT(ZM) ELSE ZM=2.*DSQRT(B) ENDIF DO 1 I=1,3 Z1=ZM*(1.+X9(I))*0.5 Z2=ZM*(1.-X9(I))*0.5 QUQ=DSQRT(B+Z1**2)-R IF (QUQ.LT.85.)XXT=XXT+A9(I)/(1.+EXP(QUQ)) QUQ=DSQRT(B+Z2**2)-R IF (QUQ.LT.85.)XXT=XXT+A9(I)/(1.+EXP(QUQ)) 1 CONTINUE XXT=XXT*ZM*0.5 DT=0. DO 2 I=1,2 Z1=X5(I)+ZM QUQ=DSQRT(B+Z1**2)-R-X5(I) IF (QUQ.LT.85.)DT=DT+A5(I)/(EXP(-X5(I))+EXP(QUQ)) 2 CONTINUE XXT=XXT+DT IF(DEBUG.GE.3)WRITE (MONIOU,202)XXROT 202 FORMAT(2X,'XXT=',E10.3) RETURN END C======================================================================= FUNCTION XXTWDEC(S,A,B) c Kinematical function for two particle decay - C light cone momentum share for c the particle of mass squared A, C B - partner's mass squared, C S - two particle invariant mass c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) INTEGER DEBUG COMMON /AREA43/ MONIOU COMMON /DEBUG/ DEBUG IF(DEBUG.GE.2)WRITE (MONIOU,201)S,A,B 201 FORMAT(2X,'XXTWDEC: S=',E10.3,2X,'A=',E10.3,2X,'B=',E10.3) X=.5D0*(1.D0+(A-B)/S) DX=(X*X-A/S) IF(DX.GT.0.D0)THEN X=X+DSQRT(DX) ELSE X=DSQRT(A/S) ENDIF XXTWDEC=X IF(DEBUG.GE.3)WRITE (MONIOU,202)XXTWDEC 202 FORMAT(2X,'XXTWDEC=',E10.3) RETURN END C======================================================================= DOUBLE PRECISION FUNCTION GAMFUN(Y) C Gamma function : See Abramowitz, page 257, form. 6.4.40 c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION(A-H,O-Z) DOUBLE PRECISION + Y,R,S,T,AFSPL,X, + COEF(10),PI,ZEROD,HALFD,ONED,TWOD,TEND C DATA COEF/8.3333333333333334D-02,-2.7777777777777778D-03, . 7.9365079365079365D-04,-5.9523809523809524D-04, . 8.4175084175084175D-04,-1.9175269175269175D-03, . 6.4102564102564103D-03,-2.9550653594771242D-02, . 0.1796443723688306 ,-0.6962161084529506 / DATA PI/ 3.141592653589793D0/ DATA ZEROD/0.D0/,HALFD/0.5D0/,ONED/1.D0/,TWOD/2.D0/,TEND/10.D0/ C X=Y AFSPL=ONED N=INT(TEND-Y) DO 10 I=0,N AFSPL=AFSPL*X X=X+ONED 10 CONTINUE R=(X-HALFD)* LOG(X)-X+HALFD* LOG(TWOD*PI) S=X T=ZEROD DO 20 I=1,10 T=T+COEF(I)/S S=S*X**2 20 CONTINUE GAMFUN = EXP(R+T)/AFSPL END C======================================================================= BLOCK DATA PSDATA c Constants for numerical integration (Gaussian weights) c----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) COMMON /AR3/ X1(7),A1(7) COMMON /AR5/ X5(2),A5(2) COMMON /AR8/ X2(4),A2 COMMON /AR9/ X9(3),A9(3) DATA X1/.9862838D0,.9284349D0,.8272013D0,.6872929D0,.5152486D0, * .3191124D0,.1080549D0/ DATA A1/.03511946D0,.08015809D0,.1215186D0,.1572032D0, * .1855384D0,.2051985D0,.2152639D0/ DATA X2/.00960736D0,.0842652D0,.222215D0,.402455D0/ DATA A2/.392699D0/ DATA X5/.585786D0,3.41421D0/ DATA A5/.853553D0,.146447D0/ DATA X9/.93247D0,.661209D0,.238619D0/ DATA A9/.171324D0,.360762D0,.467914D0/ END