* * $Id$ * * $Log$ * Revision 1.2 1996/04/18 15:30:15 ravndal * PCM index overflow protection * * Revision 1.1.1.1 1995/10/24 10:21:13 cernlib * Geant * * #include "geant321/pilot.h" *CMZ : 3.21/02 29/03/94 15.41.40 by S.Giani *FCA : 05/01/99 12:38:13 by Federico Carminati * Added protection for uninitialised variables *-- Author : SUBROUTINE HIGCLU(IPPP,NFL,AVERN) C C *** GENERATION OF X- AND PT- VALUES FOR ALL PRODUCED PARTICLES *** C *** NVE 01-AUG-1988 CERN GENEVA *** C C ORIGIN : H.FESEFELDT (11-OCT-1987) C C A SIMPLE TWO CLUSTER MODEL IS USED C THIS SHOULD BE SUFFICIENT FOR LOW ENERGY INTERACTIONS C C #include "geant321/s_defcom.inc" #include "geant321/s_genio.inc" C REAL NUCSUP DIMENSION SIDE(MXGKCU),C1PAR(5),G1PAR(5),NUCSUP(6) DIMENSION RNDM(3) DATA C1PAR/0.6,0.6,0.35,0.15,0.10/ DATA G1PAR/2.6,2.6,1.8,1.30,1.20/ DATA NUCSUP/1.0,0.7,0.5,0.4,0.35,0.3/ DIMENSION PSUP(6) DATA PSUP/3.,6.,20.,50.,100.,1000./ C DATA CB/3.0/ DATA CB/0.01/ BPP(X)=4.000+1.600*LOG(X) C 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 EK=ENP(5) EN=ENP(6) P=ENP(7) S=ENP(8) RS=ENP(9) CFA=0.025*((ATNO2-1.)/120.)*EXP(-(ATNO2-1.)/120.) SPALL=0 IF(P.LT.0.001) GOTO 60 NT=0 C** C** CHECK MASS-INDICES FOR ALL PARTICLES C** 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) C** C** SET THE EFFECTICE 4-MOMENTUM-VECTOR FOR INTERACTION C** 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(48)=IER(48)+1 C** C** DISTRIBUTE PARTICLES IN FORWARD AND BACKWARD HEMISPHERE OF CMS C** OF THE HADRON NUCLEON INTERACTION C** SIDE(1)= 1. SIDE(2)=-1. TARG=0. IFOR=1 IBACK=1 DO 3 I=1,NT IF (I .LE. 2) GO TO 78 SIDE(I)=1. CALL GRNDM(RNDM,1) IF (RNDM(1) .LT. 0.5) SIDE(I)=-1. IF (SIDE(I) .LT. 0.) GO TO 76 C C --- PARTICLE IN FORWARD HEMISPHERE --- 77 CONTINUE IFOR=IFOR+1 IF (IFOR .LE. 18) GO TO 78 C C --- CHANGE IT TO BACKWARD --- SIDE(I)=-1. IFOR=IFOR-1 IBACK=IBACK+1 GO TO 78 C C --- PARTICLE IN BACKWARD HEMISPHERE --- 76 CONTINUE IBACK=IBACK+1 IF (IBACK .LE. 18) GO TO 78 C C --- CHANGE IT TO FORWARD --- SIDE(I)=1. IBACK=IBACK-1 IFOR=IFOR+1 C** C** SUPPRESSION OF CHARGED PIONS FOR VARIOUS REASONS C** 78 IF(IPART.EQ.15.OR.IPART.GE.17) GOTO 3 IF(ABS(IPA(I)).GE.10) GOTO 3 IF(ABS(IPA(I)).EQ. 8) GOTO 3 CALL GRNDM(RNDM,1) IF(RNDM(1).GT.(10.-P)/6.) GOTO 3 CALL GRNDM(RNDM,1) IF(RNDM(1).GT.ATNO2/300.) GOTO 3 IPA(I)=14 CALL GRNDM(RNDM,1) IF(RNDM(1).GT.ZNO2/ATNO2) IPA(I)=16 TARG=TARG+1. 3 CONTINUE TB=2.*IBACK CALL GRNDM(RNDM,1) IF(RS.LT.(2.0+RNDM(1))) TB=(2.*IBACK+NT)/2. C** C** NUCLEONS + SOME PIONS FROM INTRANUCLEAR CASCADE C** AFC=0.312+0.200*LOG(LOG(S))+S**1.5/6000. IF(AFC.GT.0.50) AFC= 0.50 XTARG=AFC*(ATNO2**0.33-1.0)*TB IF(XTARG.LE.0.) XTARG=0.01 DO 881 IPX=1,6 IF(P.LE.PSUP(IPX)) GOTO 882 881 CONTINUE IPX=6 882 XPNHMF = XTARG*NUCSUP(IPX) XSHHMF = XTARG - XPNHMF IF(XSHHMF.LT.0.01) XSHHMF=0.01 IF(XPNHMF.LT.0.01) XPNHMF=0.01 SSHHMF=0.5*XSHHMF SPNHMF=0.9*XPNHMF RSHHMF=SSHHMF**2/XSHHMF RPNHMF=SPNHMF**2/XPNHMF IF(RSHHMF.LT.1.1) THEN CALL POISSO(XSHHMF,NSHHMF) GOTO 541 ELSE RSHHMF=XSHHMF/(RSHHMF-1.) IF(RSHHMF.LE.20.) THEN CALL SVGAM7(RSHHMF,XHMF) ELSE KRSHMF=IFIX(RSHHMF+0.5) CALL SVERL2(KRSHMF,XHMF) ENDIF XSHHMF=XHMF*XSHHMF/RSHHMF CALL POISSO(XSHHMF,NSHHMF) ENDIF 541 IF(RPNHMF.LE.1.1) THEN CALL POISSO(XPNHMF,NPNHMF) GOTO 542 ELSE RPNHMF=XPNHMF/(RPNHMF-1.) IF(RPNHMF.LE.20.) THEN CALL SVGAM7(RPNHMF,XHMF) ELSE KRPHMF=IFIX(RPNHMF+0.5) CALL SVERL2(KRPHMF,XHMF) ENDIF XPNHMF=XHMF*XPNHMF/RPNHMF CALL POISSO(XPNHMF,NPNHMF) ENDIF 542 NTARG=NSHHMF+NPNHMF NT2=NT+NTARG IF(NT2.LE.MXGKPV-30) GOTO 2 NT2=MXGKPV-30 NTARG=NT2-NT 2 CONTINUE IF(NPRT(4)) *WRITE(NEWBCD,3001) NTARG,NT NT1=NT+1 IF(NTARG.EQ.0) GOTO 51 DO 4 I=NT1,NT2 IF(NPNHMF.GT.0) GOTO 52 CALL GRNDM(RNDM,1) IPA(I)=-(7+IFIX(RNDM(1)*3.0)) SIDE(I)=-1. 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. SIDE(I)=-2. NPNHMF=NPNHMF-1 4 CONTINUE 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** MARK LEADING STRANGE PARTICLES C** LEAD=0 IF(IPART.LT.10.OR.IPART.EQ.14.OR.IPART.EQ.16) GOTO 6 IPA1=ABS(IPA(1)) IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 531 LEAD=IPA1 GOTO 6 531 IPA1=ABS(IPA(2)) IF(IPA1.LT.10.OR.IPA1.EQ.14.OR.IPA1.EQ.16) GOTO 6 LEAD=IPA1 C** C** CHECK AVAILABLE KINETIC ENERGY , CHANGE HEMISPHERE FOR PARTICLES C** UNTIL IT FITS C** 6 IF(NT.LE.1) GOTO 60 TAVAI=0. DO 7 I=1,NT IF(SIDE(I).LT.-1.5) GOTO 7 TAVAI=TAVAI+ABS(PV(5,I)) 7 CONTINUE IF(TAVAI.LT.RS) GOTO 12 IF(NPRT(4)) *WRITE(NEWBCD,3002) (IPA(I),I=1,20),(SIDE(I),I=1,20),TAVAI,RS 3002 FORMAT(' *HIGCLU* CHECK AVAILABLE ENERGIES'/ * 1H ,20I5/1H ,20F5.0/1H ,'TAVAI,RS ',2F10.3) DO 10 I=1,NT II=NT-I+1 IF(SIDE(II).LT.-1.5) GOTO 10 IF(II.EQ.NT) GOTO 11 NT1=II+1 NT2=NT DO 8 J=NT1,NT2 IPA(J-1)=IPA(J) SIDE(J-1)=SIDE(J) DO 8 K=1,10 8 PV(K,J-1)=PV(K,J) GOTO 11 10 CONTINUE 11 SIDE(NT)=0. IPA(NT)=0 NT=NT-1 GOTO 6 12 IF(NT.LE.1) GOTO 60 B=BPP(P) IF(B.LT.CB) B=CB C** C** CHOOSE MASSES FOR THE 3 CLUSTER: 1. FORWARD CLUSTER C** 2. BACKWARD MESON CLUSTER 3. BACKWARD NUCLEON CLUSTER C** RMC0=0. RMD0=0. RME0=0. NTC=0 NTD=0 NTE=0 DO 31 I=1,NT IF(SIDE(I).GT.0.) RMC0=RMC0+ABS(PV(5,I)) IF(SIDE(I).GT.0.) NTC =NTC +1 IF(SIDE(I).LT.0..AND.SIDE(I).GT.-1.5) RMD0=RMD0+ABS(PV(5,I)) IF( SIDE(I).LT.-1.5) RME0=RME0+ABS(PV(5,I)) IF(SIDE(I).LT.0..AND.SIDE(I).GT.-1.5) NTD =NTD +1 IF( SIDE(I).LT.-1.5) NTE =NTE +1 31 CONTINUE 32 CALL GRNDM(RNDM,1) RAN=RNDM(1) RMC=RMC0 IF(NTC.LE.1) GOTO 33 NTC1=NTC IF(NTC1.GT.5) NTC1=5 RMC=-LOG(1.-RAN) GPAR=G1PAR(NTC1) CPAR=C1PAR(NTC1) DUMNVE=GPAR IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10 RMC=RMC0+RMC**CPAR/DUMNVE 33 RMD=RMD0 IF(NTD.LE.1) GOTO 34 NTD1=NTD IF(NTD1.GT.5) NTD1=5 CALL GRNDM(RNDM,1) RAN=RNDM(1) RMD=-LOG(1.-RAN) GPAR=G1PAR(NTD1) CPAR=C1PAR(NTD1) DUMNVE=GPAR IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10 RMD=RMD0+RMD**CPAR/DUMNVE 34 IF(RMC+RMD.LE.RS) GOTO 35 IF (RMC.LE.RMC0.AND.RMD.LE.RMD0) THEN HNRMDC = 0.999*RS/(RMC+RMD) RMD = RMD*HNRMDC RMC = RMC*HNRMDC ELSE RMC=0.1*RMC0+0.9*RMC RMD=0.1*RMD0+0.9*RMD END IF GOTO 34 35 IF(NTE.LE.0) GOTO 38 RME=RME0 IF(NTE.EQ.1) GOTO 38 NTE1=NTE IF(NTE1.GT.5) NTE1=5 CALL GRNDM(RNDM,1) RAN=RNDM(1) RME=-LOG(1.-RAN) GPAR=G1PAR(NTE1) CPAR=C1PAR(NTE1) DUMNVE=GPAR IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10 RME=RME0+RME**CPAR/DUMNVE C** C** SET BEAM , TARGET OF FIRST INTERACTION IN CMS C** 38 PV( 1,MX1) =0. PV( 2,MX1) =0. PV( 3,MX1) =P PV( 5,MX1) =ABS(AMAS) PV( 4,MX1) =SQRT(P*P+AMAS*AMAS) PV( 1,MX2) =0. PV( 2,MX2) =0. PV( 3,MX2) =0. PV( 4,MX2) =MP PV( 5,MX2) =MP C** TRANSFORM INTO CMS. CALL ADD(MX1,MX2,MX) CALL LOR(MX1,MX,MX1) CALL LOR(MX2,MX,MX2) PF=(S+RMD*RMD-RMC*RMC)**2 - 4*S*RMD*RMD IF(PF.LT.0.0001) PF=0.0001 DUMNVE=2.0*RS IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10 PF=SQRT(PF)/DUMNVE IF(NPRT(4)) WRITE(6,2002) PF,RMC,RMD,RS C** C** SET FINAL STATE MASSES AND ENERGIES IN CMS C** PV(5,MX3) =RMC PV(5,MX4) =RMD PV(4,MX3) =SQRT(PF*PF+RMC*RMC) PV(4,MX4) =SQRT(PF*PF+RMD*RMD) C** C** SET |T| AND |TMIN| C** T=-1.0E10 CALL GRNDM(RNDM,1) IF (B .NE. 0.0) T=LOG(1.-RNDM(1))/B CALL LENGTX(MX1,PIN) TACMIN=(PV(4,MX1) -PV(4,MX3))**2 -(PIN-PF)**2 C** C** CACULATE (SIN(TETA/2.)**2 AND COS(TETA), SET AZIMUTH ANGLE PHI C** DUMNVE=4.0*PIN*PF IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10 CTET=-(T-TACMIN)/DUMNVE CTET=1.0-2.0*CTET IF (CTET .GT. 1.0) CTET=1.0 IF (CTET .LT. -1.0) CTET=-1.0 DUMNVE=1.0-CTET*CTET IF (DUMNVE .LT. 0.0) DUMNVE=0.0 STET=SQRT(DUMNVE) CALL GRNDM(RNDM,1) PHI=RNDM(1)*TWPI C** C** CALCULATE FINAL STATE MOMENTA IN CMS C** PV(1,MX3) =PF*STET*SIN(PHI) PV(2,MX3) =PF*STET*COS(PHI) PV(3,MX3) =PF*CTET PV(1,MX4) =-PV(1,MX3) PV(2,MX4) =-PV(2,MX3) PV(3,MX4) =-PV(3,MX3) C** C** SIMULATE BACKWARD NUCLEON CLUSTER IN LAB. SYSTEM AND TRANSFORM IN C** CMS. C** IF(NTE.EQ.0) GOTO 28 GA=1.2 EKIT1=0.04 EKIT2=0.6 IF(EK.GT.5.) GOTO 666 EKIT1=EKIT1*EK**2/25. EKIT2=EKIT2*EK**2/25. 666 A=(1.-GA)/(EKIT2**(1.-GA)-EKIT1**(1.-GA)) DO 29 I=1,NT IF(SIDE(I).GT.-1.5) GOTO 29 CALL GRNDM(RNDM,3) RAN=RNDM(1) EKIT=(RAN*(1.-GA)/A+EKIT1**(1.-GA))**(1./(1.-GA)) PV(4,I)=EKIT+PV(5,I) DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2) PP=SQRT(DUMNVE) RAN=RNDM(2) COST=LOG(2.23*RAN+0.383)/0.96 IF (COST .LT. -1.0) COST=-1.0 IF (COST .GT. 1.0) COST=1.0 DUMNVE=1.0-COST*COST IF (DUMNVE .LT. 0.0) DUMNVE=0.0 SINT=SQRT(DUMNVE) PHI=TWPI*RNDM(3) PV(1,I)=PP*SINT*SIN(PHI) PV(2,I)=PP*SINT*COS(PHI) PV(3,I)=PP*COST CALL LOR(I,MX,I) 29 CONTINUE C** C** FRAGMENTATION OF FORWARD CLUSTER AND BACKWARD MESON CLUSTER C** 28 PV(1,1)=PV(1,MX3) PV(2,1)=PV(2,MX3) PV(3,1)=PV(3,MX3) PV(4,1)=PV(4,MX3) PV(1,2)=PV(1,MX4) PV(2,2)=PV(2,MX4) PV(3,2)=PV(3,MX4) PV(4,2)=PV(4,MX4) DO 17 I=MX5,MX6 DO 16 J=1,3 16 PV(J,I)=-PV(J,I-2) DO 17 J=4,5 17 PV(J,I)= PV(J,I-2) KGENEV=1 IF(NTC.LE.1) GOTO 26 TECM=PV(5,MX3) NPG=0 DO 18 I=1,NT IF(SIDE(I).LT.0.) GOTO 18 IF(NPG.EQ.18) THEN SIDE(I)=-SIDE(I) GOTO 18 ENDIF NPG=NPG+1 AMASS(NPG)=ABS(PV(5,I)) 18 CONTINUE IF(NPRT(4)) WRITE(NEWBCD,2004) TECM,NPG,(AMASS(I),I=1,NPG) CALL PHASP NPG=0 DO 19 I=1,NT IF(SIDE(I).LT.0.OR.NPG.GE.18) GOTO 19 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) IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5) CALL LOR(I,MX5,I) IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I) 19 CONTINUE 26 IF(NTD.LE.1) GOTO 27 TECM=PV(5,MX4) NPG=0 DO 20 I=1,NT IF(SIDE(I).GT.0..OR.SIDE(I).LT.-1.5) GOTO 20 IF(NPG.EQ.18) THEN SIDE(I)=-2. PV(4,I)=ABS(PV(5,I)) DO 48 J=1,3 PV(J,I)=0. 48 CONTINUE GOTO 20 ENDIF NPG=NPG+1 AMASS(NPG)=ABS(PV(5,I)) 20 CONTINUE IF(NPRT(4)) WRITE(NEWBCD,2004) TECM,NPG,(AMASS(I),I=1,NPG) CALL PHASP NPG=0 DO 21 I=1,NT IF(SIDE(I).GT.0..OR.SIDE(I).LT.-1.5) GOTO 21 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) IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,5) CALL LOR(I,MX6,I) IF(NPRT(4)) WRITE(NEWBCD,2001) I,(PV(J,I),J=1,10),IPA(I),SIDE(I) 21 CONTINUE C** C** LORENTZ TRANSFORMATION IN LAB SYSTEM C** 27 TARG=0. DO 36 I=1,NT IF(PV(5,I).GT.0.5) TARG=TARG+1. CALL LOR(I,MX2,I) 36 CONTINUE IF(TARG.LT.0.5) TARG=1. C** C** SOMETIMES THE LEADING STRANGE PARTICLES ARE LOST , SET THEM BACK C** 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) DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2) PP1=SQRT(DUMNVE) C IF (PP .GE. 1.0E-6) GO TO 8000 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 8001 8000 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 8001 CONTINUE C C** FOR VARIOUS REASONS, THE ENERGY BALANCE IS NOT SUFFICIENT, C** CHECK THAT, ENERGY BALANCE, ANGLE OF FINAL SYSTEM E.T.C. 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** THE 3. CLUSTER IS NOT PRODUCED WITHIN PROPER KINEMATICS!!! 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) DUMNVE=ABS(PV(4,I)**2-PV(5,I)**2) PP=SQRT(DUMNVE) CALL LENGTX(I,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,I)=PP*SIN(RTHNVE)*COS(PHINVE) PV(2,I)=PP*SIN(RTHNVE)*SIN(PHINVE) PV(3,I)=PP*COS(RTHNVE) GO TO 8003 8002 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 8003 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, SEE COMMENTS IN 'GENXPT' 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) * 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 C** DEKIN=0. NPIONS=0 EK1=0. DO 25 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 25 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 8004 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 8005 8004 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 8005 CONTINUE C 25 CONTINUE IF(EK1.EQ.0.) GOTO 23 IF(NPIONS.LE.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 8006 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 8007 8006 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 8007 CONTINUE C 22 CONTINUE 23 IF(ATNO2.LT.1.5) GOTO 40 C** C** ADD BLACK TRACK PARTICLES C** CALL HIGHAB(SPROB) CALL GRNDM(RNDM,1) IF(RNDM(1).LT.SPROB) GOTO 40 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-2) NBL=MXGKPV-2-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-2) 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.0+RNDM(2)*2.0 DUMNVE=1.0-COST*COST IF (DUMNVE .LT. 0.0) DUMNVE=0.0 SINT=SQRT(DUMNVE) 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) DUMNVE=ABS(PV(4,NT)**2-PV(5,NT)**2) PP=SQRT(DUMNVE) 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 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-2) NBL=MXGKPV-2-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-2) 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.005*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.0+RNDM(1)*2.0 DUMNVE=1.0-COST*COST IF (DUMNVE .LT. 0.0) DUMNVE=0.0 SINT=SQRT(DUMNVE) 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 DUMNVE=ABS(PV(4,NT)**2-PV(5,NT)**2) PP=SQRT(DUMNVE) 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)*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) DUMNVE=ETB IF (DUMNVE .EQ. 0.0) DUMNVE=1.0E-10 ETF=ETF/DUMNVE 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 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. CALL GHETUN(NT) DO 45 I=1,NT EKIN2=EKIN2+PV(4,I)-ABS(PV(5,I)) 45 CONTINUE EKIN2=(EKIN2-EKIN)/EKIN IF(NPRT(4)) $ WRITE(NEWBCD,2006) NT,EKIN,ENP(1),ENP(3),EKIN1,EKIN2 IF(EKIN2.GT.0.2) GOTO 60 INTCT=INTCT+1. NMODE=3 IF(SPALL.LT.0.5.AND.ATNO2.GT.1.5) NMODE=14 CALL SETCUR(NT) NTK=NTK+1 IF(NT.EQ.1) GOTO 300 DO 50 II=2,NT I=II-1 IF(NTOT.LT.NSIZE/12) GOTO 43 GO TO 9999 43 CALL SETTRK(I) 50 CONTINUE 300 CONTINUE GO TO 9999 C** C** IT IS NOT POSSIBLE TO PRODUCE A PROPER TWO CLUSTER FINAL STATE. C** CONTINUE WITH QUASI ELASTIC SCATTERING C** 60 IF(NPRT(4)) WRITE(NEWBCD,2005) DO 61 I=3,MXGKCU 61 IPA(I)=0 IPA(1)=IPART IPA(2)=14 IF(NFL.EQ.2) IPA(2)=16 CALL TWOB(IPPP,NFL,AVERN) GO TO 9999 C 2000 FORMAT(' *HIGCLU* CMS PARAMETERS OF FINAL STATE PARTICLES', $ ' AFTER ',I3,' TRIALS') 2001 FORMAT(' *HIGCLU* TRACK',2X,I3,2X,10F8.2,2X,I3,2X,F3.0) 2002 FORMAT(' *HIGCLU* MOMENTUM ',F8.3,' MASSES ',2F8.4,' RS ',F8.4) 2003 FORMAT(' *HIGCLU* TETA,EKIN0,EKIN1,EKIN ',4F10.4) 2004 FORMAT(' *HIGCLU* TECM,NPB,MASSES: ',F10.4,1X,I3,1X,8F10.4/ $ 1H ,26X,15X,8F10.4) 2005 FORMAT(' *HIGCLU* NUMBER OF FINAL STATE PARTICLES', $ ' LESS THAN 2 ==> CONTINUE WITH 2-BODY SCATTERING') 2006 FORMAT(' *HIGCLU* COMP.',1X,I5,1X,5F7.2) 3001 FORMAT(' *HIGCLU* NUCLEAR EXCITATION ',I5,' PARTICLES PRODUCED', $ ' IN ADDITION TO',I5,' NORMAL PARTICLES') 3003 FORMAT(' *HIGCLU* ',I3,' BLACK TRACK PARTICLES PRODUCED', $ ' WITH TOTAL KINETIC ENERGY OF ',F8.3,' GEV') 3004 FORMAT(' *HIGCLU* ',I5,' HEAVY FRAGMENTS WITH TOTAL ENERGY OF ', $ F8.4,' GEV') C 9999 CONTINUE END