* * $Id$ * * $Log$ * Revision 1.1.1.1 1995/10/24 10:21:04 cernlib * Geant * * #include "geant321/pilot.h" *CMZ : 3.21/02 29/03/94 15.41.39 by S.Giani *-- Author : SUBROUTINE GENXPT(IPPP,NFL,AVERN) C C *** GENERATION OF X- AND PT- VALUES FOR ALL PRODUCED PARTICLES *** C *** NVE 02-MAY-1988 CERN GENEVA *** C C ORIGIN : H.FESEFELDT 11-OCT-1987 C C A SIMPLE SINGLE VARIABLE DESCRIPTION E D3S/DP3= F(Q) WITH C Q**2 = (M*X)**2 + PT**2 IS USED. FINAL STATE KINEMATIC IS PRODUCED C BY AN FF-TYPE ITERATIVE CASCADE METHOD C #include "geant321/s_defcom.inc" #include "geant321/s_genio.inc" C REAL MASPAR,LAMB,NUCSUP DIMENSION MASPAR(8),BP(8),PTEX(8),C1PAR(5),G1PAR(5),TAVAI(2), $ SIDE(MXGKCU),IAVAI(2),BINL(20),DNDL(20),TWSUP(8), $ NUCSUP(6),PSUP(6),IPAX(100) DIMENSION RNDM(3) DATA MASPAR/0.75,0.70,0.65,0.60,0.50,0.40,0.75,0.20/ DATA BP/3.50,3.50,3.50,6.00,5.00,4.00,3.50,3.50/ DATA PTEX/1.70,1.70,1.50,1.70,1.40,1.20,1.70,1.20/ DATA C1PAR/0.6,0.6,0.35,0.15,0.10/ DATA G1PAR/2.6,2.6,1.80,1.30,1.20/ DATA BINL/0.,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.11,1.25 $ ,1.43,1.67,2.0,2.5,3.33,5.00,10.00/ DATA TWSUP/1.,1.,0.7,0.5,0.3,0.2,0.1,0.0/ DATA NUCSUP/1.00,0.7,0.5,0.4,0.35,0.3/ DATA PSUP/3.,6.,20.,50.,100.,1000./ C C** CALL HIGSEL(ISEL) IF(ISEL.EQ.1) THEN CALL HIGXPT(IPPP,NFL,AVERN) RETURN ENDIF C** C** FOR ANNIHILATION INTERACTIONS INTRODUCE PROPER KINEMATICS C** CALL CORANH(NIHIL,NFL) C** C** C** CHECK FIRST MASS-INDICES C** EK=ENP(5) EN=ENP(6) P=ENP(7) S=ENP(8) RS=ENP(9) NT=0 DO 1 I=1,100 IF(IPA(I).EQ.0) GOTO 1 NT=NT+1 IPA(NT)=IPA(I) 1 CONTINUE CALL VZERO(IPA(NT+1),MXGKCU-NT) CALL UCOPY(IPA(1),IPAX(1),100) C** C** FOR LOW MULTIPLICITY USE TWO-BODY RESONANCE MODEL OR SINGLE/DOUBLE C** DIFFRACTION MODEL (--> TWOCLU (--> TWOB (--> COSCAT))) C** CFA=0.025*((ATNO2-1.)/120.)*EXP(-(ATNO2-1.)/120.) IF(NIHIL.GT.0) GOTO 200 IF(NT.GE.8) GOTO 200 IF(EK.LT.1.) GOTO 60 CALL GRNDM(RNDM,1) RAN=RNDM(1) IF(IPART.GE.10.AND.IPART.LE.13.AND.RAN.LT.0.5) GOTO 200 CALL GRNDM(RNDM,1) RAN=RNDM(1) WSUP=TWSUP(NT) IF(RAN.GT.WSUP) GOTO 200 60 CALL UCOPY(IPAX,IPA,100) CALL TWOCLU(IPPP,NFL,AVERN) GO TO 9999 C** C** SET EFFECTIVE 4-MOMENTUM OF PRIMARY PARTICLE C** 200 MX =MXGKPV-20 MX1=MX+1 MX2=MX+2 MX3=MX+3 MX4=MX+4 MX5=MX+5 MX6=MX+6 MX7=MX+7 MX8=MX+8 MX9=MX+9 PV( 1,MXGKPV-1)=P*PX PV( 2,MXGKPV-1)=P*PY PV( 3,MXGKPV-1)=P*PZ PV( 4,MXGKPV-1)=EN PV( 5,MXGKPV-1)=AMAS PV( 6,MXGKPV-1)=NCH PV( 7,MXGKPV-1)=TOF PV( 8,MXGKPV-1)=IPART PV( 9,MXGKPV-1)=0. PV(10,MXGKPV-1)=USERW IER(49)=IER(49)+1 C** C** SOME RANDOMISATION OF ORDER OF FINAL STATE PARTICLES C** DO 201 I=3,NT CALL GRNDM(RNDM,1) IPX=IFIX(3.+RNDM(1)*(NT-2.)) IF(IPX.GT.NT) IPX=NT IPA1=IPA(IPX) IPA(IPX)=IPA(I) 201 IPA(I) =IPA1 C** C** DISTRIBUTE IN FORWARD AND BACKWARD HEMISPHERE IN CMS C** SIDE(1)= 1. SIDE(2)=-1. NTB=1 TARG=0. IF(IPART.LT.10.OR.IPART.GT.13) GOTO 53 CALL GRNDM(RNDM,1) IF(RNDM(1).LT.0.7) GOTO 53 IPA1=IPA(1) IPA(1)=IPA(2) IPA(2)=IPA1 53 LEAD=0 IF(IPART.LT.10.OR.IPART.EQ.14.OR.IPART.EQ.16) GOTO 532 IPA1=ABS(IPA(1)) IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 531 LEAD=IPA1 GOTO 532 531 IPA1=ABS(IPA(2)) IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 532 LEAD=IPA1 532 DO 3 I=1,NT IF(I.LE.2) GOTO 54 SIDE(I)= 1. CALL GRNDM(RNDM,1) IF(RNDM(1).LT.0.5) SIDE(I)=-1. IF(SIDE(I).LT.-0.5) NTB=NTB+1 54 CONTINUE 3 CONTINUE TB=2.*NTB CALL GRNDM(RNDM,1) IF(RS.LT.(2.0+RNDM(1))) TB=(2.*NTB+NT)/2. C** C** ADD PARTICLES FROM INTRANUCLEAR CASCADE C** AFC=0.312+0.200*LOG(LOG(S))+S**1.5/6000. IF(AFC.GT.0.75) AFC=0.75 XTARG=AFC*(ATNO2**0.33 -1.0)*TB IF(XTARG.LE.0.) XTARG=0.01 CALL POISSO(XTARG,NTARG) NT2=NT+NTARG IF(NT2.LE.MX) GOTO 2 NT2=MX NTARG=NT2-NT 2 CONTINUE IF (NPRT(4)) WRITE(NEWBCD,3001) NTARG,NT NT1=NT+1 IF(NTARG.EQ.0) GOTO 51 C** C** CHECK NUMBER OF EXTRA NUCLEONS AND PIONS C** DO 881 IPX=1,6 IF(P.LE.PSUP(IPX)) GOTO 882 881 CONTINUE IPX=6 882 DO 4 I=NT1,NT2 CALL GRNDM(RNDM,1) RAN=RNDM(1) IF(RAN.LT.NUCSUP(IPX)) GOTO 52 CALL GRNDM(RNDM,1) IPA(I)=-(7+IFIX(RNDM(1)*3.0)) GOTO 4 52 IPA(I)=-16 PNRAT=1.-ZNO2/ATNO2 CALL GRNDM(RNDM,1) IF(RNDM(1).GT.PNRAT) IPA(I)=-14 TARG=TARG+1. 4 SIDE(I)=-2. NT=NT2 C** C** CHOOSE MASSES AND CHARGES FOR ALL PARTICLES C** 51 DO 5 I=1,NT IPA1=ABS(IPA(I)) PV(5,I)=RMASS(IPA1) PV(6,I)=RCHARG(IPA1) PV(7,I)=1. IF(PV(5,I).LT.0.) PV(7,I)=-1. PV(5,I)=ABS(PV(5,I)) 5 CONTINUE C** C** CHECK AVAILABLE KINETIC ENERGY, IN THIS MODEL CONSERVATION OF C** KINETIC ENERGY IN FORWARD AND BACKWARD HEMISPHERE IS ASSUMED C** 6 IF(NT.LE.1) GOTO 60 TAVAI(1)=RS/2. TAVAI(2)=(TARG+1.)*RS/2. IAVAI(1)=0 IAVAI(2)=0 DO 7 I=1,NT L=1 IF(SIDE(I).LT.0.) L=2 IAVAI(L)=IAVAI(L)+1 TAVAI(L)=TAVAI(L)-ABS(PV(5,I)) 7 CONTINUE NTH=NT IF(NTH.GT.10) NTH=10 IF (NPRT(4)) $ WRITE(NEWBCD,3002) TAVAI,IAVAI,(IPA(I),SIDE(I),I=1,NTH) IF(IAVAI(1).LE.0) GOTO 60 IF(IAVAI(2).LE.0) GOTO 60 IF(TAVAI(1).GT.0.) GOTO 11 CALL GRNDM(RNDM,1) ISKIP=IFIX(RNDM(1)*(IAVAI(1)-1))+1 IS=0 DO 10 I=1,NT II=NT-I+1 IF(SIDE(II).LT.0.) GOTO 10 IS=IS+1 IF(IS.NE.ISKIP) GOTO 10 IF(II.EQ.NT) GOTO 9 NT1=II+1 NT2=NT DO 8 J=NT1,NT2 IPA(J-1)=IPA(J) SIDE(J-1)=SIDE(J) DO 71 K=1,10 71 PV(K,J-1)=PV(K,J) 8 CONTINUE GOTO 9 10 CONTINUE 9 IPA(NT)=0 SIDE(NT)=0. NT=NT-1 GOTO 6 11 IF(TAVAI(2).GT.0.) GOTO 15 CALL GRNDM(RNDM,1) ISKIP=IFIX(RNDM(1)*(IAVAI(2)-1))+1 IS=0 DO 14 I=1,NT II=NT-I+1 IF(SIDE(II).GT.0.) GOTO 14 IS=IS+1 IF(IS.NE.ISKIP) GOTO 14 IF(SIDE(II).LT.-1.5) NTARG=NTARG-1 IF(NTARG.LT.0) NTARG=0 IF(II.EQ.NT) GOTO 13 NT1=II+1 NT2=NT DO 12 J=NT1,NT2 IPA(J-1)=IPA(J) SIDE(J-1)=SIDE(J) DO 74 K=1,10 74 PV(K,J-1)=PV(K,J) 12 CONTINUE GOTO 13 14 CONTINUE 13 IPA(NT)=0 SIDE(NT)=0. NT=NT-1 GOTO 6 15 IF(NT.LE.1) GOTO 60 IF(NT.EQ.MX) GOTO 29 NT1=NT+1 NT2=MX DO 28 I=NT1,NT2 28 IPA(I)=0 29 CONTINUE C** C** NOW THE PREPARATION IS FINISHED. C** DEFINE INITIAL STATE VECTORS FOR LORENTZ TRANSFORMATIONS. C** PV( 1,MX1)=0. PV( 2,MX1)=0. PV( 3,MX1)=P PV( 4,MX1)=SQRT(P*P+AMAS*AMAS) PV( 5,MX1)=ABS(AMAS) PV( 1,MX2)=0. PV( 2,MX2)=0. PV( 3,MX2)=0. PV( 4,MX2)=MP PV( 5,MX2)=MP PV( 1,MX4)=0. PV( 2,MX4)=0. PV( 3,MX4)=0. PV( 4,MX4)=MP*(1.+TARG) PV( 5,MX4)=PV(4,MX4) PV( 1,MX8)=0. PV( 2,MX8)=0. PV( 3,MX8)=0. PV( 1,MX9)=1. PV( 2,MX9)=0. PV( 3,MX9)=0. CALL ADD(MX1,MX2,MX3) CALL ADD(MX4,MX1,MX4) CALL LOR(MX1,MX3,MX1) CALL LOR(MX2,MX3,MX2) C** C** MAIN LOOP FOR 4-MOMENTUM GENERATION , SEE PITHA-REPORT (AACHEN) C** FOR A DETAILED DESCRIPTION OF THE METHOD. C** CALL GRNDM(RNDM,1) PHI=RNDM(1)*TWPI EKIN1=0. EKIN2=0. DO 39 J=1,10 PV(J,MX5)=0. 39 PV(J,MX6)=0. NPG=0 TARG1=0. DO 16 III=1,NT I=NT-III+1 IPA1=ABS(IPA(I)) C** C** COUNT NUMBER OF BACKWARD NUCLEONS C** IF(I.EQ.2) GOTO 301 IF(SIDE(I).LT.-1.5.AND.IPA1.GE.14) GOTO 301 GOTO 38 301 NPG=NPG+1 IF(NPG.GT.18) GOTO 38 SIDE(I)=-3. TARG1=TARG1+1. GOTO 16 38 J=3 IF(IPA1.LT.14) J=2 IF(IPA1.LT.10) J=1 IF(I.LE.2) J=J+3 IF(SIDE(I).LT.-1.5) J=7 IF(J.EQ.7.AND.IPA1.GE.14) J=8 C** C** SET PT - AND PHI VALUES, THEY ARE CHANGED SOMEWHAT IN THE ITERATION C** LOOP, SET MASS PARAMETER FOR LAMBDA FRAGMENTATION MODEL C** CALL GRNDM(RNDM,1) RAN=RNDM(1) BPP=BP(J) BPE=PTEX(J) PT2=-LOG(1.-RAN)/BPP ASPAR=MASPAR(J) PT2=PT2**BPE PT =SQRT(PT2) IF(PT.LT.0.001) PT=0.001 PV(1,I)=PT*COS(PHI) PV(2,I)=PT*SIN(PHI) PV(10,I)=PT BINL(1)=0. RLMAX=1./PV(10,I) DO 73 J=2,20 73 BINL(J)=RLMAX*(J-1)/19. ET=PV(4,MX1) IF(SIDE(I).LT.0.) THEN ET=PV(4,MX2) ENDIF DNDL(1)=0. NTRIAL=0 C** C** START OF BIG ITERATION LOOP C** 30 NTRIAL=NTRIAL+1 IF(NTRIAL.GT. 2) GOTO 169 DO 17 L=2,20 DNDL(L)=0. X=(BINL(L)+BINL(L-1))/2. IF(PV(10,I).LT.0.001) PV(10,I)=0.001 IF(X.GT.1./PV(10,I)) GOTO 17 DX=BINL(L)-BINL(L-1) DNDL(L)=ASPAR/SQRT((1.+(ASPAR*X)**2)**3) DNDL(L)=ET*DNDL(L)/SQRT((X*PV(10,I)*ET)**2+PV(10,I)**2 * +PV(5,I)**2) DNDL(L)=DNDL(L)*DX 17 DNDL(L)=DNDL(L-1)+DNDL(L) NTRI=0 31 CALL GRNDM(RNDM,1) RAN=RNDM(1)*DNDL(20) DO 18 L=2,20 IF(RAN.LT.DNDL(L)) GOTO 19 18 CONTINUE C** C** START OF SMALL ITERATION LOOP C** 19 NTRI=NTRI+1 CALL GRNDM(RNDM,1) RAN=RNDM(1) DX=BINL(L)-BINL(L-1) LAMB=BINL(L-1)+RAN*DX/2. X=PV(10,I)*LAMB IF(X.GT.1.) X=1. X=X*SIDE(I)/ABS(SIDE(I)) PV(3,I)=X*ET PV(4,I)=PV(3,I)**2+PV(10,I)**2+PV(5,I)**2 PV(4,I)=SQRT(PV(4,I)) IF(SIDE(I).LT.0.) GOTO 165 IF(I.GT.2) GOTO 20 EKIN=TAVAI(1)-EKIN1 CALL NORMAL(RAN) IF(EKIN.LT.0.) EKIN=0.04*ABS(RAN) PV(4,I)=ABS(PV(5,I))+EKIN RNVE=ABS(PV(4,I)**2-PV(5,I)**2) PP=SQRT(RNVE) CALL LENGTX(I,PP1) C IF (PP1 .GE. 1.0E-6) GO TO 8000 CALL GRNDM(RNDM,2) RTHNVE=PI*RNDM(1) PHINVE=TWPI*RNDM(2) PV(1,I)=PP*SIN(RTHNVE)*COS(PHINVE) PV(2,I)=PP*SIN(RTHNVE)*SIN(PHINVE) PV(3,I)=PP*COS(RTHNVE) GO TO 8001 8000 CONTINUE PV(1,I)=PV(1,I)*PP/PP1 PV(2,I)=PV(2,I)*PP/PP1 PV(3,I)=PV(3,I)*PP/PP1 8001 CONTINUE C CALL ADD(MX5,I,MX5) GOTO 16 20 EKIN=EKIN1+PV(4,I)-ABS(PV(5,I)) IF(EKIN.LT.0.95*TAVAI(1)) GOTO 161 IF(NTRI.GT. 5) GOTO 167 PV(10,I)=PV(10,I)*0.9 PV( 1,I)=PV( 1,I)*0.9 PV( 2,I)=PV( 2,I)*0.9 DNDL(20)=DNDL(20)*0.9 IF((TAVAI(2)-ABS(PV(5,I))).LT.0.) GOTO 31 SIDE(I)=-1. TAVAI(1)=TAVAI(1)+ABS(PV(5,I)) TAVAI(2)=TAVAI(2)-ABS(PV(5,I)) GOTO 31 161 CALL ADD(MX5,I,MX5) EKIN1=EKIN1+PV(4,I)-ABS(PV(5,I)) GOTO 163 165 EKIN=EKIN2+PV(4,I)-ABS(PV(5,I)) XXX=0.95+0.05*TARG/20. IF(XXX.GT.0.999) X=0.999 IF(EKIN.LT.XXX*TAVAI(2)) GOTO 166 IF(NTRI.GT. 5) GOTO 167 PV(10,I)=PV(10,I)*0.9 PV( 1,I)=PV( 1,I)*0.9 PV( 2,I)=PV( 2,I)*0.9 DNDL(20)=DNDL(20)*0.9 IF((TAVAI(1)-ABS(PV(5,I))).LT.0.) GOTO 31 SIDE(I)=+1. TAVAI(1)=TAVAI(1)-ABS(PV(5,I)) TAVAI(2)=TAVAI(2)+ABS(PV(5,I)) GOTO 31 166 CALL ADD(MX6,I,MX6) EKIN2=EKIN2+PV(4,I)-ABS(PV(5,I)) 163 CALL ADD(MX5,MX6,MX7) PV(3,MX7)=0. CALL ANG(MX7,MX9,COST,PHIS) IF(PV(2,MX7).LT.0.) PHIS=TWPI-PHIS CALL NORMAL(RAN) RAN=RAN*PI/12. PHI=PHIS+PI+RAN IF(PHI.GT.TWPI) PHI=PHI-TWPI IF(PHI.LT.0.) PHI=TWPI-PHI GOTO 16 C** C** PARTICLE MOMENTUM ZERO, REDUCE KINETIC ENERGY OF ALL OTHER C** 167 EKIN1=0. EKIN2=0. DO 162 J=1,10 PV(J,MX5)=0. 162 PV(J,MX6)=0. II=I+1 DO 168 L=II,NT IF(ABS(IPA(L)).GE.14.AND.SIDE(L).LT.0.) GOTO 168 PV(4,L)=PV(4,L)*0.95+0.05*ABS(PV(5,L)) IF(PV(4,L).LT.ABS(PV(5,L))) PV(4,L)=ABS(PV(5,L)) RNVE=ABS(PV(4,L)**2-PV(5,L)**2) PP=SQRT(RNVE) CALL LENGTX(L,PP1) C IF (PP1 .GE. 1.0E-6) GO TO 8002 CALL GRNDM(RNDM,2) RTHNVE=PI*RNDM(1) PHINVE=TWPI*RNDM(2) PV(1,L)=PP*SIN(RTHNVE)*COS(PHINVE) PV(2,L)=PP*SIN(RTHNVE)*SIN(PHINVE) PV(3,L)=PP*COS(RTHNVE) GO TO 8003 8002 CONTINUE PV(1,L)=PV(1,L)*PP/PP1 PV(2,L)=PV(2,L)*PP/PP1 PV(3,L)=PV(3,L)*PP/PP1 8003 CONTINUE C PV(10,L)=SQRT(PV(1,L)**2+PV(2,L)**2) IF(SIDE(L).LT.0.) GOTO 164 EKIN1=EKIN1+PV(4,L)-ABS(PV(5,L)) CALL ADD(MX5,L,MX5) GOTO 168 164 EKIN2=EKIN2+PV(4,L)-ABS(PV(5,L)) CALL ADD(MX6,L,MX6) 168 CONTINUE C *** NEXT STMT. CHANGED TO PREVENT FROM INFINITE LOOPING *** C************* GOTO 38 GO TO 30 C** C** SKIP PARTICLE, IF NOT ENOUGH ENERGY C** 169 IPA(I)=0 DO 170 J=1,10 170 PV(J,I)=0. GOTO 163 16 CONTINUE NTRI=0 II=0 DO 320 I=1,NT IF(IPA(I).EQ.0) GOTO 320 II=II+1 IPA(II)=IPA(I) SIDE(II)=SIDE(I) DO 321 J=1,10 321 PV(J,II)=PV(J,I) 320 CONTINUE NT=II C** C** BACKWARD NUCLEONS PRODUCED WITH A CLUSTER MODEL C** CALL LOR(MX4,MX3,MX7) CALL SUB(MX7,MX5,MX7) CALL SUB(MX7,MX6,MX7) IF(TARG1.GT.1.5) GOTO 310 322 I=2 CALL NORMAL(RAN) EKIN=TAVAI(2)-EKIN2 EKINM=RS/2.-MP IF(EKIN.GT.EKINM) EKIN=EKINM CALL NORMAL(RAN) IF(EKIN.LT.0.04) EKIN=0.04*ABS(RAN) PV(4,I)=ABS(PV(5,I))+EKIN RNVE=ABS(PV(4,I)**2-PV(5,I)**2) PP=SQRT(RNVE) CALL LENGTX(MX7,PP1) C IF (PP1 .GE. 1.0E-6) GO TO 8004 CALL GRNDM(RNDM,2) RTHNVE=PI*RNDM(1) PHINVE=TWPI*RNDM(2) PV(1,I)=PP*SIN(RTHNVE)*COS(PHINVE) PV(2,I)=PP*SIN(RTHNVE)*SIN(PHINVE) PV(3,I)=PP*COS(RTHNVE) GO TO 8005 8004 CONTINUE PV(1,I)=PV(1,MX7)*PP/PP1 PV(2,I)=PV(2,MX7)*PP/PP1 PV(3,I)=PV(3,MX7)*PP/PP1 8005 CONTINUE C CALL ADD(MX6,I,MX6) GOTO 330 310 ITARG1=IFIX(TARG1+0.1) IF(ITARG1.GT.5) ITARG1=5 RMB0=0. NPG=0 DO 311 I=1,NT IF(SIDE(I).GT.-2.5) GOTO 311 NPG=NPG+1 RMB0=RMB0+ABS(PV(5,I)) 311 CONTINUE IF(NPG.LT.2) GOTO 322 CALL GRNDM(RNDM,1) RAN=RNDM(1) RMB=-LOG(1.-RAN) GPAR=G1PAR(ITARG1) CPAR=C1PAR(ITARG1) RMB=RMB0+RMB**CPAR/GPAR PV(5,MX7)=RMB IF(PV(5,MX7).GT.PV(4,MX7)) PV(5,MX7)=PV(4,MX7) RNVE=ABS(PV(4,MX7)**2-PV(5,MX7)**2) PP=SQRT(RNVE) CALL LENGTX(MX7,PP1) C IF (PP1 .GE. 1.0E-6) GO TO 8006 CALL GRNDM(RNDM,2) RTHNVE=PI*RNDM(1) PHINVE=TWPI*RNDM(2) PV(1,MX7)=PP*SIN(RTHNVE)*COS(PHINVE) PV(2,MX7)=PP*SIN(RTHNVE)*SIN(PHINVE) PV(3,MX7)=PP*COS(RTHNVE) GO TO 8007 8006 CONTINUE PV(1,MX7)=PV(1,MX7)*PP/PP1 PV(2,MX7)=PV(2,MX7)*PP/PP1 PV(3,MX7)=PV(3,MX7)*PP/PP1 8007 CONTINUE C I=MX7 IF (NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5) PV(1,MX7)=-PV(1,MX7) PV(2,MX7)=-PV(2,MX7) PV(3,MX7)=-PV(3,MX7) KGENEV=1 TECM=PV(5,MX7) NPG=0 DO 312 I=1,NT IF(SIDE(I).GT.-2.5)GOTO 312 NPG=NPG+1 AMASS(NPG)=ABS(PV(5,I)) 312 CONTINUE CALL PHASP NPG=0 DO 314 I=1,NT IF(SIDE(I).GT.-2.5) GOTO 314 NPG=NPG+1 PV(1,I)=PCM(1,NPG) PV(2,I)=PCM(2,NPG) PV(3,I)=PCM(3,NPG) PV(4,I)=PCM(4,NPG) CALL LOR(I,MX7,I) CALL ADD(MX6,I,MX6) 314 CONTINUE 330 IF (NPRT(4)) $ WRITE(NEWBCD,2002) NTRIAL,EKIN1,EKIN2,TAVAI(1),TAVAI(2) 175 IF (.NOT.NPRT(4)) GOTO 36 CALL ADD(MX5,MX6,MX7) EKIN1=PV(4,MX1)+PV(4,MX2) EKIN2=PV(4,MX5)+PV(4,MX6) WRITE(NEWBCD,2000) EKIN1,EKIN2 I=MX1 WRITE(NEWBCD,2001) I,(PV(J,I),J=1,4) I=MX2 WRITE(NEWBCD,2001) I,(PV(J,I),J=1,4) I=MX5 WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5) I=MX6 WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5) DO 37 I=1,NT 37 WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I) C** C** LORENTZ TRANSFORMATION IN LAB SYSTEM C** 36 IF(NT.LE.2) GOTO 60 TARG=0. DO 601 I=1,NT IF(PV(5,I).GT.0.5) TARG=TARG+1. CALL LOR(I,MX2,I) 601 CONTINUE IF(TARG.LT.0.5) TARG=1. IF(LEAD.EQ.0) GOTO 6085 DO 6081 I=1,NT IF(ABS(IPA(I)).EQ.LEAD) GOTO 6085 6081 CONTINUE I=1 IF(LEAD.GE.14.AND.ABS(IPA(2)).GE.14) I=2 IF(LEAD.LT.14.AND.ABS(IPA(2)).LT.14) I=2 IPA(I)=LEAD EKIN=PV(4,I)-ABS(PV(5,I)) PV(5,I)=RMASS(LEAD) PV(7,I)=1. IF(PV(5,I).LT.0.) PV(7,I)=-1. PV(5,I)=ABS(PV(5,I)) PV(6,I)=RCHARG(LEAD) PV(4,I)=PV(5,I)+EKIN CALL LENGTX(I,PP) RNVE=ABS(PV(4,I)**2-PV(5,I)**2) PP1=SQRT(RNVE) PV(1,I)=PP1*PV(1,I)/PP PV(2,I)=PP1*PV(2,I)/PP PV(3,I)=PP1*PV(3,I)/PP 6085 KGENEV=1 PV(1,MX4)=0. PV(2,MX4)=0. PV(3,MX4)=P PV(4,MX4)=SQRT(P*P+AMAS*AMAS) PV(5,MX4)=ABS(AMAS) EKIN0=PV(4,MX4)-PV(5,MX4) PV(1,MX5)=0. PV(2,MX5)=0. PV(3,MX5)=0. PV(4,MX5)=MP*TARG PV(5,MX5)=PV(4,MX5) EKIN=PV(4,MX4)+PV(4,MX5) I=MX4 IF (NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5) I=MX5 IF (NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5) CALL ADD(MX4,MX5,MX6) CALL LOR(MX4,MX6,MX4) CALL LOR(MX5,MX6,MX5) TECM=PV(4,MX4)+PV(4,MX5) NPG=NT PV(1,MX8)=0. PV(2,MX8)=0. PV(3,MX8)=0. PV(4,MX8)=0. PV(5,MX8)=0. EKIN1=0. DO 598 I=1,NPG IF (NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I) CALL ADD(MX8,I,MX8) EKIN1=EKIN1+PV(4,I)-PV(5,I) EKIN=EKIN-PV(5,I) IF(I.GT.18) GOTO 598 AMASS(I)=PV(5,I) 598 CONTINUE IF(NPG.GT.18) GOTO 597 CALL PHASP EKIN=0. DO 599 I=1,NPG PV(1,MX7)=PCM(1,I) PV(2,MX7)=PCM(2,I) PV(3,MX7)=PCM(3,I) PV(4,MX7)=PCM(4,I) PV(5,MX7)=AMASS(I) CALL LOR(MX7,MX5,MX7) 599 EKIN=EKIN+PV(4,MX7)-PV(5,MX7) CALL ANG(MX8,MX4,COST,TETA) IF (NPRT(4)) WRITE(NEWBCD,2003) TETA,EKIN0,EKIN1,EKIN C** C** MAKE SHURE, THAT KINETIC ENERGIES ARE CORRECT. C** EKIN= KINETIC ENERGY THEORETICALLY C** EKIN1= KINETIC ENERGY SIMULATED C** 597 IF(EKIN1.EQ.0.) GOTO 600 PV(1,MX7)=0. PV(2,MX7)=0. PV(3,MX7)=0. PV(4,MX7)=0. PV(5,MX7)=0. WGT=EKIN/EKIN1 EKIN1=0. DO 602 I=1,NT EKIN=PV(4,I)-PV(5,I) EKIN=EKIN*WGT PV(4,I)=EKIN+PV(5,I) RNVE=ABS(PV(4,I)**2-PV(5,I)**2) PP=SQRT(RNVE) CALL LENGTX(I,PP1) C IF (PP1 .GE. 1.0E-6) GO TO 8008 CALL GRNDM(RNDM,2) RTHNVE=PI*RNDM(1) PHINVE=TWPI*RNDM(2) PV(1,I)=PP*SIN(RTHNVE)*COS(PHINVE) PV(2,I)=PP*SIN(RTHNVE)*SIN(PHINVE) PV(3,I)=PP*COS(RTHNVE) GO TO 8009 8008 CONTINUE PV(1,I)=PV(1,I)*PP/PP1 PV(2,I)=PV(2,I)*PP/PP1 PV(3,I)=PV(3,I)*PP/PP1 8009 CONTINUE C EKIN1=EKIN1+EKIN CALL ADD(MX7,I,MX7) 602 CONTINUE CALL ANG(MX7,MX4,COST,TETA) IF (NPRT(4)) WRITE(NEWBCD,2003) TETA,EKIN0,EKIN1 C** C** ROTATE IN DIRECTION OF Z-AXIS, THIS DOES DISTURB IN SOME WAY OUR C** INCLUSIVE DISTRIBUTIONS, BUT IT IS NESSACARY FOR MOMENTUM CONSER- C** VATION. C** 600 PV(1,MX7)=0. PV(2,MX7)=0. PV(3,MX7)=0. PV(4,MX7)=0. PV(5,MX7)=0. DO 596 I=1,NT 596 CALL ADD(MX7,I,MX7) C** C** SOME SMEARING IN TRANSVERSE DIRECTION FROM FERMI MOTION C** * call rannor(ran1,ran2) CALL GRNDM(RNDM,2) RY=RNDM(1) RZ=RNDM(2) RX=6.283185*RZ A1=SQRT(-2.*LOG(RY)) RAN1=A1*SIN(RX) RAN2=A1*COS(RX) PV(1,MX7)=PV(1,MX7)+RAN1*0.020*TARG PV(2,MX7)=PV(2,MX7)+RAN2*0.020*TARG CALL DEFS(MX4,MX7,MX8) PV(1,MX7)=0. PV(2,MX7)=0. PV(3,MX7)=0. PV(4,MX7)=0. PV(5,MX7)=0. DO 595 I=1,NT CALL TRAC(I,MX8,I) 595 CALL ADD(MX7,I,MX7) CALL ANG(MX7,MX4,COST,TETA) IF (NPRT(4)) WRITE(NEWBCD,2003) TETA C** C** ROTATE IN DIRECTION OF PRIMARY PARTICLE, SUBTRACT BINDING ENERGIES C** AND MAKE SOME FURTHER CORRECTIONS IF REQUIRED (STEEP, STEEQ) C** DEKIN=0. NPIONS=0 EK1=0. DO 21 I=1,NT CALL DEFS1(I,MXGKPV-1,I) IF (NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I) IF(ATNO2.LT.1.5) GOTO 21 CALL LENGTX(I,PP) EKIN=PV(4,I)-ABS(PV(5,I)) CALL NORMAL(RAN) EKIN=EKIN-CFA*(1.+0.5*RAN) IF (EKIN .LT. 1.0E-6) EKIN=1.0E-6 CALL STEEQ(XXH,I) DEKIN=DEKIN+EKIN*(1.-XXH) EKIN=EKIN*XXH IF(ABS(IPA(I)).GE.7.AND.ABS(IPA(I)).LE.9) NPIONS=NPIONS+1 IF(ABS(IPA(I)).GE.7.AND.ABS(IPA(I)).LE.9) EK1=EK1+EKIN PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I)))) PV(4,I)=EKIN+ABS(PV(5,I)) C IF (PP .GE. 1.0E-6) GO TO 8010 CALL GRNDM(RNDM,2) RTHNVE=PI*RNDM(1) PHINVE=TWPI*RNDM(2) PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE) PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE) PV(3,I)=PP1*COS(RTHNVE) GO TO 8011 8010 CONTINUE PV(1,I)=PV(1,I)*PP1/PP PV(2,I)=PV(2,I)*PP1/PP PV(3,I)=PV(3,I)*PP1/PP 8011 CONTINUE C 21 CONTINUE IF(EK1.EQ.0.) GOTO 23 IF(NPIONS.EQ.0) GOTO 23 DEKIN=1.+DEKIN/EK1 DO 22 I=1,NT IF(ABS(IPA(I)).LT.7.OR.ABS(IPA(I)).GT.9) GOTO 22 CALL LENGTX(I,PP) EKIN=PV(4,I)-ABS(PV(5,I)) EKIN=EKIN*DEKIN IF (EKIN .LT. 1.0E-6) EKIN=1.0E-6 PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I)))) PV(4,I)=EKIN+ABS(PV(5,I)) C IF (PP .GE. 1.0E-6) GO TO 8012 CALL GRNDM(RNDM,2) RTHNVE=PI*RNDM(1) PHINVE=TWPI*RNDM(2) PV(1,I)=PP1*SIN(RTHNVE)*COS(PHINVE) PV(2,I)=PP1*SIN(RTHNVE)*SIN(PHINVE) PV(3,I)=PP1*COS(RTHNVE) GO TO 8013 8012 CONTINUE PV(1,I)=PV(1,I)*PP1/PP PV(2,I)=PV(2,I)*PP1/PP PV(3,I)=PV(3,I)*PP1/PP 8013 CONTINUE C 22 CONTINUE C** C** ADD BLACK TRACK PARTICLES, THE TOTAL NUMBER OF PARTICLES PRODUCED C** IS RESTRICTED TO 198, THIS MAY HAVE INFLUENCE ON VERY HIGH ENERGY C** FIRST PROTONS AND NEUTRONS C** 23 IF(ATNO2.LT.1.5) GOTO 40 CALL SELFAB(SPROB) TEX=ENP(1) SPALL=TARG IF(TEX.LT.0.001) GOTO 445 BLACK=(1.5+1.25*TARG)*ENP(1)/(ENP(1)+ENP(3)) CALL POISSO(BLACK,NBL) IF (NPRT(4)) WRITE(NEWBCD,3003) NBL,TEX IF(IFIX(TARG)+NBL.GT.ATNO2) NBL=ATNO2-TARG IF(NT+NBL.GT.MXGKPV-10) NBL=MXGKPV-10-NT IF(NBL.LE.0) GOTO 445 EKIN=TEX/NBL EKIN2=0. CALL STEEP(XX) DO 441 I=1,NBL CALL GRNDM(RNDM,1) IF(RNDM(1).LT.SPROB) GOTO 441 IF(NT.EQ.MXGKPV-10) GOTO 441 IF(EKIN2.GT.TEX) GOTO 443 CALL GRNDM(RNDM,1) RAN1=RNDM(1) CALL NORMAL(RAN2) EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2) IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1) EKIN1=EKIN1*XX EKIN2=EKIN2+EKIN1 IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1) IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6 IPA1=16 PNRAT=1.-ZNO2/ATNO2 CALL GRNDM(RNDM,3) IF(RNDM(1).GT.PNRAT) IPA1=14 NT=NT+1 SPALL=SPALL+1. COST=-1.+RNDM(2)*2. SINT=SQRT(ABS(1.-COST*COST)) PHI=TWPI*RNDM(3) IPA(NT)=-IPA1 SIDE(NT)=-4. PV(5,NT)=ABS(RMASS(IPA1)) PV(6,NT)=RCHARG(IPA1) PV(7,NT)=1. PV(4,NT)=EKIN1+PV(5,NT) RNVE=ABS(PV(4,NT)**2-PV(5,NT)**2) PP=SQRT(RNVE) PV(1,NT)=PP*SINT*SIN(PHI) PV(2,NT)=PP*SINT*COS(PHI) PV(3,NT)=PP*COST 441 CONTINUE 443 IF(ATNO2.LT.10.) GOTO 445 IF(EK.GT.2.0) GOTO 445 II=NT+1 KK=0 EKA=EK IF(EKA.GT.1.) EKA=EKA*EKA IF(EKA.LT.0.1) EKA=0.1 IKA=3.6*EXP((ZNO2**2/ATNO2-35.56)/6.45)/EKA IF(IKA.LE.0) GO TO 445 DO 444 I=1,NT II=II-1 IF(IPA(II).NE.-14) GOTO 444 IPA(II)=-16 IPA1 = 16 PV(5,II)=ABS(RMASS(IPA1)) PV(6,II)=RCHARG(IPA1) KK=KK+1 IF(KK.GT.IKA) GOTO 445 444 CONTINUE C** C** THEN ALSO DEUTERONS, TRITONS AND ALPHAS C** 445 TEX=ENP(3) IF(TEX.LT.0.001) GOTO 40 BLACK=(1.5+1.25*TARG)*ENP(3)/(ENP(1)+ENP(3)) CALL POISSO(BLACK,NBL) IF(NT+NBL.GT.MXGKPV-10) NBL=MXGKPV-10-NT IF(NBL.LE.0) GOTO 40 EKIN=TEX/NBL EKIN2=0. CALL STEEP(XX) IF (NPRT(4)) WRITE(NEWBCD,3004) NBL,TEX DO 442 I=1,NBL CALL GRNDM(RNDM,1) IF(RNDM(1).LT.SPROB) GOTO 442 IF(NT.EQ.MXGKPV-10) GOTO 442 IF(EKIN2.GT.TEX) GOTO 40 CALL GRNDM(RNDM,1) RAN1=RNDM(1) CALL NORMAL(RAN2) EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2) IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1) EKIN1=EKIN1*XX EKIN2=EKIN2+EKIN1 IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1) IF (EKIN1 .LT. 0.0) EKIN1=1.0E-6 CALL GRNDM(RNDM,3) COST=-1.+RNDM(1)*2. SINT=SQRT(ABS(1.-COST*COST)) PHI=TWPI*RNDM(2) RAN=RNDM(3) IPA(NT+1)=-30 IF(RAN.GT.0.60) IPA(NT+1)=-31 IF(RAN.GT.0.90) IPA(NT+1)=-32 SIDE(NT+1)=-4. PV(5,NT+1)=(ABS(IPA(NT+1))-28)*MP SPALL=SPALL+PV(5,NT+1)*1.066 IF(SPALL.GT.ATNO2) GOTO 40 NT=NT+1 PV(6,NT)=1. IF(IPA(NT).EQ.-32) PV(6,NT)=2. PV(7,NT)=1. PV(4,NT)=PV(5,NT)+EKIN1 RNVE=ABS(PV(4,NT)**2-PV(5,NT)**2) PP=SQRT(RNVE) PV(1,NT)=PP*SINT*SIN(PHI) PV(2,NT)=PP*SINT*COS(PHI) PV(3,NT)=PP*COST 442 CONTINUE C** C** STORE ON EVENT COMMON C** 40 CALL GRNDM(RNDM,1) IF(RS.GT.(4.+RNDM(1))) GOTO 42 DO 41 I=1,NT CALL LENGTX(I,ETB) IF(ETB.LT.P) GOTO 41 ETF=P PV(4,I)=SQRT(PV(5,I)**2+ETF**2) ETF=ETF/ETB PV(1,I)=PV(1,I)*ETF PV(2,I)=PV(2,I)*ETF PV(3,I)=PV(3,I)*ETF 41 CONTINUE 42 EKIN=PV(4,MXGKPV)-ABS(PV(5,MXGKPV)) EKIN1=PV(4,MXGKPV-1)-ABS(PV(5,MXGKPV-1)) EKIN2=0. CALL TDELAY(TOF1) CALL GRNDM(RNDM,1) RAN=RNDM(1) TOF=TOF-TOF1*LOG(RAN) DO 44 I=1,NT EKIN2=EKIN2+PV(4,I)-ABS(PV(5,I)) IF(PV(7,I).LT.0.) PV(5,I)=-PV(5,I) PV(7,I)=TOF PV(8,I)=ABS(IPA(I)) PV(9,I)=0. 44 PV(10,I)=0. IF (NPRT(4)) WRITE(NEWBCD,2006) NT,EKIN,ENP(1),ENP(3),EKIN1,EKIN2 INTCT=INTCT+1. CALL SETCUR(NT) NTK=NTK+1 IF(NT.EQ.1) GO TO 9999 DO 50 II=2,NT I=II-1 IF(NTOT.LT.NSIZE/12) GOTO 43 GO TO 9999 43 CALL SETTRK(I) 50 CONTINUE C 2002 FORMAT(' *GENXPT* PRODUCTION OF FINAL STATE KINEMATIC AFTER ',I3, $ ' TRIALS. KINETIC ENERGIES ',2F6.2,' OUT OF ',2F6.2) 2000 FORMAT(' *GENXPT* CMS PARAMETERS OF FINAL STATE PARTICLES,', $ ' ENERGIES IN INITIAL AND FINAL STATE ',2F6.2) 2001 FORMAT(' *GENXPT* TRACK',2X,I3,2X,10F8.3,2X,I3,2X,F4.0) 2003 FORMAT(' *GENXPT* TETA,EKIN0,EKIN1,EKIN ',4F10.4) 2006 FORMAT(' *GENXPT* COMP.',1X,I5,1X,5F7.2) 3001 FORMAT(' *GENXPT* NUCLEAR EXCITATION',I5, $ ' PARTICLES PRODUCED IN ADDITION TO ',I5,' NORMAL PARTICLES') 3002 FORMAT(' *GENXPT* AVAILABLE ENERGIES ',2F10.4, $ ' FOR ',2I3,' PARTICLES IN BEAM/TARGET FRAGM. REGION', $ ' WITH IPA/SIDE ARRAY '/ $ 1H ,5X,10(I3,2X,F3.0,4X)) 3003 FORMAT(' *GENXPT* ',I3,' BLACK TRACK PARTICLES PRODUCED', $ ' WITH TOTAL KINETIC ENERGY OF ',F8.3,' GEV') 3004 FORMAT(' *GENXPT* ',I5,' HEAVY FRAGMENTS PRODUCED', $ ' WITH TOTAL ENERGY OF',F8.4,' GEV') C 9999 CONTINUE C END