* * $Id$ * * $Log$ * Revision 1.1.1.1 1995/10/24 10:20:59 cernlib * Geant * * #include "geant321/pilot.h" *CMZ : 3.21/02 29/03/94 15.41.38 by S.Giani *-- Author : SUBROUTINE TWOB(IPPP,NFL,AVERN) C C *** GENERATION OF MOMENTA FOR ELAST. AND QUASI ELAST. 2 BODY REACT. *** C *** NVE 04-MAY-1988 CERN GENEVA *** C C ORIGIN : H.FESEFELDT 15-SEP-1987 C C THE SIMPLE FORMULA DS/D|T| = S0* EXP(-B*|T|) IS USED. C THE B VALUES ARE PARAMETRIZATIONS FROM EXPERIMENTAL DATA . C NOT AVAILABLE VALUES ARE TAKEN FROM THOSE OF SIMILAR REACTIONS C #include "geant321/s_defcom.inc" DIMENSION RNDM(3) C C DATA CB/3./ DATA CB/0.01/ C C --- STATEMENT FUNCTIONS --- BPP(X)=4.225+1.795*LOG(X) C C** C** FOR DIFFRACTION SCATTERING ON HEAVY NUCLEI USE BETTER ROUTINE C** "COSCAT" C TARMAS=RMASS(14) IF (NFL .EQ. 2) TARMAS=RMASS(16) ENP(8)=RMASS(IPPP)**2+TARMAS**2+2.0*TARMAS*ENP(6) ENP(9)=SQRT(ENP(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.) IF(ATNO2.LT.1.5) GOTO 500 IPA1=ABS(IPA(1)) IPA2=ABS(IPA(2)) RMC=RMASS(IPA1) RMD=RMASS(IPA2) RCHC=RCHARG(IPA1) RCHD=RCHARG(IPA2) IF(ABS(RMC-AMAS).GT.0.001) GOTO 500 RMNVE=RMASS(14) IF (NFL .EQ. 2) RMNVE=RMASS(16) IF(ABS(RMD-RMNVE).GT.0.001) GOTO 500 IF(ABS(RCHC-NCH).GT.0.5) GOTO 500 IF(NFL.EQ.1.AND.RCHD.LT.0.5) GOTO 500 IF(NFL.EQ.2.AND.ABS(RCHD).GT.0.5) GOTO 500 IF(ENP(1).GT.0.0001.OR.ENP(3).GT.0.0001) GOTO 500 CALL COSCAT GO TO 9999 C** C** SET EFFECTIVE 4-MOMENTUM OF INITIAL PARTICLE C** 500 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(47)=IER(47)+1 IF(NPRT(4)) $ WRITE(NEWBCD,4001) (PV(J,MXGKPV-1),J=1,10),IPA(1),IPA(2) DO 2 J=1,6 2 PV(J,1)=PV(J,MXGKPV-1) PV(7,1)=1. IF(PV(5,1).LT.0.) PV(7,1)=-1. PV(5,1)=ABS(PV(5,1)) NT=1 C** C** TWO-BODY SCATTERING POSSIBLE?? IF NOT, CONTINUE WITH ORIGINAL C** PARTICLE, BUT SPEND THE NUCLEAR EVAPORATION ENERGY C** IF(P.LT.0.1) GOTO 200 IF(RS.LT.0.01) GOTO 200 C** C** CALCULATE SLOPE B FOR ELASTIC SCATTERING ON PROTON/NEUTRON C** B=BPP(P) IF(B.LT.CB) B=CB IF(ABS(IPA(2)).GT.13) GOTO 9 IPA(2)=14 CALL GRNDM(RNDM,1) IF(RNDM(1).LT.0.5) IPA(2)=16 C** C** SET MASSES AND MOMENTA FOR FINAL STATE PARTICLES C** 9 RMC=RMASS(ABS(IPA(1))) RMD=RMASS(ABS(IPA(2))) PV(6,1)=RCHARG(ABS(IPA(1))) PV(6,2)=RCHARG(ABS(IPA(2))) PF=(S+RMD*RMD-RMC*RMC)**2 - 4*S*RMD*RMD IF(NPRT(4)) WRITE(NEWBCD,4002) RMC,RMD,PV(6,1),PV(6,2),RS,S,PF IF(PF.LT.0.001) GO TO 9999 PF=SQRT(PF)/(2.*RS) C** C** SET BEAM AND TARGET IN CMS C** PV(1,3)=0. PV(2,3)=0. PV(3,3)=P PV(5,3)=ABS(AMAS) PV(4,3)=SQRT(P*P+AMAS*AMAS) PV(1,4)=0. PV(2,4)=0. PV(3,4)=0. RMNVE=RMASS(14) IF (NFL .EQ. 2) RMNVE=RMASS(16) PV(4,4)=RMNVE PV(5,4)=RMNVE C** C** TRANSFORM INTO CMS. C** CALL ADD(3,4,10) CALL LOR(3,10,3) CALL LOR(4,10,4) C** C** SET FINAL STATE MASSES AND ENERGIES IN CMS C** PV(5,1)=ABS(RMC) PV(5,2)=ABS(RMD) PV(7,1)=1. PV(7,2)=1. IF(RMC.LT.0.) PV(7,1)=-1. IF(RMD.LT.0.) PV(7,2)=-1. PV(4,1)=SQRT(PF*PF+PV(5,1)*PV(5,1)) PV(4,2)=SQRT(PF*PF+PV(5,2)*PV(5,2)) C** C** SET |T| AND |TMIN| C** CALL GRNDM(RNDM,2) CALL LENGTX(3,PIN) BTRANG=B*4.*PIN*PF C** C** SIMPLY A PROTECTION AGAINST EXPONENT OVERFLOW 1.E20 IS BIG ENOUGH C** EXINDT=-1. IF(BTRANG.LT.46) EXINDT=EXINDT+EXP(-BTRANG) TDN=LOG(1.+RNDM(1)*EXINDT)/BTRANG C** C** CACULATE (SIN(TETA/2.)**2 AND COS(TETA), SET AZIMUTH ANGLE PHI C** CTET=1.+2.*TDN IF(ABS(CTET).GT.1.) CTET=SIGN(1.,CTET) STET=SQRT((1.-CTET)*(1.+CTET)) PHI=RNDM(2)*TWPI C** C** CALCULATE FINAL STATE MOMENTA IN CMS C** PV(1,1)=PF*STET*SIN(PHI) PV(2,1)=PF*STET*COS(PHI) PV(3,1)=PF*CTET PV(1,2)=-PV(1,1) PV(2,2)=-PV(2,1) PV(3,2)=-PV(3,1) C** C** TRANSFORM INTO LAB C** DO 11 I=1,2 CALL LOR(I,4,I) CALL DEFS1(I,MXGKPV-1,I) IF(ATNO2.LT.1.5) GOTO 11 CALL LENGTX(I,PP) IF(PP.LT.0.001) GOTO 11 EKIN=PV(4,I)-ABS(PV(5,I)) CALL NORMAL(RAN) EKIN=EKIN-CFA*(1.+0.5*RAN) IF(EKIN.LT.0.0001) EKIN=0.0001 PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I)))) PV(4,I)=EKIN+ABS(PV(5,I)) PV(1,I)=PV(1,I)*PP1/PP PV(2,I)=PV(2,I)*PP1/PP PV(3,I)=PV(3,I)*PP1/PP 11 CONTINUE NT=2 C** C** ADD BLACK TRACK PARTICLES . C** HERE THE PROCEDURE IS SOMEWHAT DIFFERENT AS IN 'TWOCLU' AND 'GENXPT' C** THE REASON IS, THAT WE HAVE TO SIMULATE ALSO THE NUCLEAR REACTIONS C** AT LOW ENERGIES LIKE A(H,P)B, A(H,P P)B, A(H,N)B E.T.C. C** 200 IF(ENP(1).LE.0.0001.AND.ENP(3).LE.0.0001) GOTO 40 SPALL=0. TEX=ENP(1) IF(TEX.LT.0.0001) GOTO 445 BLACK=TEX/0.02 CALL POISSO(BLACK,NBL) IF(NBL.GT.ATNO2) NBL=ATNO2 IF(ENP(1).GT.0.0001.AND.NBL.LE.0) NBL=1 IF (NPRT(4)) WRITE(NEWBCD,3003) NBL,TEX 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 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.) EKIN1=0.0001 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 PV(5,NT)=ABS(RMASS(IPA1)) PV(6,NT)=RCHARG(IPA1) PV(7,NT)=2. PV(4,NT)=EKIN1+PV(5,NT) PP=SQRT(ABS(PV(4,NT)**2-PV(5,NT)**2)) 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.0001) GOTO 40 NBL=IFIX(2.*LOG(ATNO2)) 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 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.010*LOG(RAN1) EKIN1=EKIN1*XX EKIN2=EKIN2+EKIN1 IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1) IF(EKIN1.LT.0.) EKIN1=0.0001 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 INVE=ABS(IPA(NT+1)) PV(5,NT+1)=RMASS(INVE) SPALL=SPALL+PV(5,NT+1)*1.066 IF(SPALL.GT.ATNO2) GOTO 40 NT=NT+1 PV(6,NT)=RCHARG(INVE) PV(7,NT)=2. PV(4,NT)=PV(5,NT)+EKIN1 PP=SQRT(ABS(PV(4,NT)**2-PV(5,NT)**2)) 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 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 1 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. 1 PV(10,I)=0. IF (NPRT(4)) WRITE(NEWBCD,1003) NT,EKIN,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 RETURN 43 CALL SETTRK(I) 50 CONTINUE C 1002 FORMAT(' *TWOB* ',5F10.4,10X,5F10.4/1H ,7X,5F10.4,10X,5F10.4/ $ ' LAB SYSTEM FINAL STATE FOUR VECTORS') 1003 FORMAT(' *TWOB* COMPARISON',2X,I5,2X,3F10.4) 4001 FORMAT(' *TWOB* ',10F10.4,2X,2I3) 4002 FORMAT(' *TWOB* ',7F10.4) 3003 FORMAT(' *TWOB* ',I3,' BLACK TRACK PARTICLES PRODUCED', $ ' WITH TOTAL KINETIC ENERGY OF ',F8.3,' GEV') 3004 FORMAT(' *TWOB* ',I5,' HEAVY FRAGMENTS PRODUCED', $ ' WITH TOTAL ENERGY OF',F8.4,' GEV') C 9999 CONTINUE C END