SUBROUTINE NIHILA C----------------------------------------------------------------------- C (AN)NIHILA(TION) C C TREATES ANNIHILATION OF ANTINUCLEONS WITH FREE NUCLEONS C MOMENTA CONSERVED IN ALL DIRECTIONS C ENERGY CONSERVED BY MULTPLICATION OF ALL MOMENTA WITH A CORRECTION C FACTOR, CONSERVING MOMENTUM BALANCE C THIS SUBROUTINE IS CALLED FROM BOX60, BOX61, BOX62, AND BOX63 C C DESIGN : D. HECK IK3 FZK KARLSRUHE C----------------------------------------------------------------------- IMPLICIT DOUBLE PRECISION (A-H,O-Z) *KEEP,ANNI. COMMON /ANNI/ CAN,CANN DOUBLE PRECISION CAN(50),CANN(50) *KEEP,CONST. COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER *KEEP,MULT. COMMON /MULT/ EKINL,MSMM,MULTMA,MULTOT DOUBLE PRECISION EKINL INTEGER MSMM,MULTMA(37,13),MULTOT(37,13) *KEEP,PAM. COMMON /PAM/ PAMA,SIGNUM DOUBLE PRECISION PAMA(6000),SIGNUM(6000) *KEEP,PARPAR. COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C, * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14), * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH INTEGER ITYPE,LEVL *KEEP,PARPAE. DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE), * (CURPAR(4), PHI ), (CURPAR(5), H ), * (CURPAR(6), T ), (CURPAR(7), X ), * (CURPAR(8), Y ), (CURPAR(9), CHI ), * (CURPAR(10),BETA), (CURPAR(11),GCM ), * (CURPAR(12),ECM ) *KEEP,RANDPA. COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR DOUBLE PRECISION FAC,U1,U2 REAL RD(3000) INTEGER ISEED(103,10),NSEQ LOGICAL KNOR *KEEP,RUNPAR. COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB, * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN, * MONIOU,MDEBUG,NUCNUC, * CETAPE, * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL, * N1STTR,MDBASE, * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR, * FIX1I,FMUADD,FNKG,FPRINT,FDBASE * ,GHEISH,GHESIG COMMON /RUNPAC/ DSN,HOST,USER DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB REAL STEPFC INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC, * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL, * N1STTR,MDBASE INTEGER CETAPE CHARACTER*79 DSN CHARACTER*20 HOST,USER LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR, * FIX1I,FMUADD,FNKG,FPRINT,FDBASE * ,GHEISH,GHESIG *KEEP,VKIN. COMMON /VKIN/ BETACM DOUBLE PRECISION BETACM *KEND. DOUBLE PRECISION E(10),PHIPAR(10),PL(10),PTR(10),PTSQ(10) DOUBLE PRECISION PX(10),PY(10) DIMENSION ISEQ(10),NTYP(10) C----------------------------------------------------------------------- IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9) 444 FORMAT(' NIHILA: CURPAR=',1P,9E10.3) IREP = 0 C RANDOM DECISION FOR ANNIHILATION WITH PROTON OR NEUTRON 20 CONTINUE CALL RMMAR( RD,2,1 ) IREP = IREP + 1 C AFTER THE 5TH TRY, QUIT THE ANNIHILATION WITHOUT ANY PION GENERATED IF ( IREP .GT. 5 ) GOTO 999 IF ( RD(1) .LE. 0.5 ) THEN C----------------------------------------------------------------------- C ANTIPROTON - PROTON AND ANTINEUTRON - NEUTRON ANNIHILATION IF ( RD(2) .LE. CANN(1) ) THEN C ANNIHILATION INTO PI+, PI- NPIPOS = 1 NPINEG = 1 NPIZ = 0 ELSEIF ( RD(2) .LE. CANN(2) ) THEN C ANNIHILATION INTO PI+, PI-, PI0 NPIPOS = 1 NPINEG = 1 NPIZ = 1 ELSEIF ( RD(2) .LE. CANN(3) ) THEN C ANNIHILATION INTO PI+, PI-, 2PI0 NPIPOS = 1 NPINEG = 1 NPIZ = 2 ELSEIF ( RD(2) .LE. CANN(4) ) THEN C ANNIHILATION INTO PI+, PI-, 3PI0 NPIPOS = 1 NPINEG = 1 NPIZ = 3 ELSEIF ( RD(2) .LE. CANN(5) ) THEN C ANNIHILATION INTO PI+, PI-, 4PI0 NPIPOS = 1 NPINEG = 1 NPIZ = 4 ELSEIF ( RD(2) .LE. CANN(6) ) THEN C ANNIHILATION INTO 2PI+, 2PI- NPIPOS = 2 NPINEG = 2 NPIZ = 0 ELSEIF ( RD(2) .LE. CANN(7) ) THEN C ANNIHILATION INTO 2PI+, 2PI-, PI0 NPIPOS = 2 NPINEG = 2 NPIZ = 1 ELSEIF ( RD(2) .LE. CANN(8) ) THEN C ANNIHILATION INTO 2PI+, 2PI-, 2PI0 NPIPOS = 2 NPINEG = 2 NPIZ = 2 ELSEIF ( RD(2) .LE. CANN(9) ) THEN C ANNIHILATION INTO 2PI+, 2PI-, 3PI0 NPIPOS = 2 NPINEG = 2 NPIZ = 3 ELSEIF ( RD(2) .LE. CANN(10) ) THEN C ANNIHILATION INTO 3PI+, 3PI- NPIPOS = 3 NPINEG = 3 NPIZ = 0 ELSEIF ( RD(2) .LE. CANN(11) ) THEN C ANNIHILATION INTO 3PI+, 3PI-, PI0 NPIPOS = 3 NPINEG = 3 NPIZ = 1 ELSEIF ( RD(2) .LE. CANN(12) ) THEN C ANNIHILATION INTO 3PI+, 3PI-, 2PI0 NPIPOS = 3 NPINEG = 3 NPIZ = 2 ELSE C ANNIHILATION INTO 4PI0 NPIPOS = 0 NPINEG = 0 NPIZ = 4 ENDIF ELSE C----------------------------------------------------------------------- C ANTIPROTON - NEUTRON (OR ANTINEUTRON - PROTON) ANNIHILATION IF ( RD(2) .LE. CANN(13) ) THEN C ANNIHILATION INTO PI-, PI0 NPIPOS = 0 NPINEG = 1 NPIZ = 1 ELSEIF ( RD(2) .LE. CANN(14) ) THEN C ANNIHILATION INTO PI-, 2PI0 NPIPOS = 0 NPINEG = 1 NPIZ = 2 ELSEIF ( RD(2) .LE. CANN(15) ) THEN C ANNIHILATION INTO PI-, 3PI0 NPIPOS = 0 NPINEG = 1 NPIZ = 3 ELSEIF ( RD(2) .LE. CANN(16) ) THEN C ANNIHILATION INTO PI-, 4PI0 NPIPOS = 0 NPINEG = 1 NPIZ = 4 ELSEIF ( RD(2) .LE. CANN(17) ) THEN C ANNIHILATION INTO PI-, 5PI0 NPIPOS = 0 NPINEG = 1 NPIZ = 5 ELSEIF ( RD(2) .LE. CANN(18) ) THEN C ANNIHILATION INTO PI+, 2PI- NPIPOS = 1 NPINEG = 2 NPIZ = 0 ELSEIF ( RD(2) .LE. CANN(19) ) THEN C ANNIHILATION INTO PI+, 2PI-, PI0 NPIPOS = 1 NPINEG = 2 NPIZ = 1 ELSEIF ( RD(2) .LE. CANN(20) ) THEN C ANNIHILATION INTO PI+, 2PI-, 2PI0 NPIPOS = 1 NPINEG = 2 NPIZ = 2 ELSEIF ( RD(2) .LE. CANN(21) ) THEN C ANNIHILATION INTO PI+, 2PI-, 3PI0 NPIPOS = 1 NPINEG = 2 NPIZ = 3 ELSEIF ( RD(2) .LE. CANN(22) ) THEN C ANNIHILATION INTO PI+, 2PI-, 4PI0 NPIPOS = 1 NPINEG = 2 NPIZ = 4 ELSEIF ( RD(2) .LE. CANN(23) ) THEN C ANNIHILATION INTO 2PI+, 3PI- NPIPOS = 2 NPINEG = 3 NPIZ = 0 ELSEIF ( RD(2) .LE. CANN(24) ) THEN C ANNIHILATION INTO 2PI+, 3PI-, PI0 NPIPOS = 2 NPINEG = 3 NPIZ = 1 ELSEIF ( RD(2) .LE. CANN(25) ) THEN C ANNIHILATION INTO 2PI+, 3PI-, 2PI0 NPIPOS = 2 NPINEG = 3 NPIZ = 2 ELSEIF ( RD(2) .LE. CANN(26) ) THEN C ANNIHILATION INTO 2PI+, 3PI-, 3PI0 NPIPOS = 2 NPINEG = 3 NPIZ = 3 ELSE C ANNIHILATION INTO 3PI+, 4PI- NPIPOS = 3 NPINEG = 4 NPIZ = 0 ENDIF C CHARGE INVERSION IF ANTINEUTRON ANNIHILATES WITH PROTON IF ( ITYPE .EQ. 25 ) THEN NPINEG = NPINEG - 1 NPIPOS = NPIPOS + 1 ENDIF ENDIF NPI = NPIPOS + NPINEG + NPIZ FNPI = 1.D0 / NPI GCMI = 1.D0 / GCM C----------------------------------------------------------------------- C CHARGE ASSIGNMENT DO 26 I = 1,NPI IF ( I .LE. NPIZ ) THEN C NEUTRAL PIONS NTYP(I) = 7 ELSEIF ( I .LE. NPIZ+NPIPOS ) THEN C POSITIVE PIONS NTYP(I) = 8 ELSE C NEGATIVE PIONS NTYP(I) = 9 ENDIF 26 CONTINUE C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C KINEMATIC CALCULATIONS ISCALE = 0 27 CONTINUE ISCALE = ISCALE + 1 C AFTER THE 5TH TRY, TAKE A NEW SET OF PIONS IF ( ISCALE .GT. 5 ) GOTO 20 C DISTRIBUTUION OF TRANSVERSE MOMENTA PTR CORECT = 5.3333333333D0 * FNPI**1.5D0 * ECM DO 28 I = 1,NPI PTR(I) = PTRANS(DUMMY) * (1.33333333D0 + CORECT) 28 CONTINUE SUMPX = 0.D0 SUMPY = 0.D0 CALL RMMAR( RD,NPI,1 ) DO 29 I = 1,NPI C SELECT EMISSION ANGLE BY REDUCED RESIDUAL DIRECTION IF ( SUMPX .NE. 0.D0 .OR. SUMPY .NE. 0.D0 ) THEN PHISUM = ATAN2( SUMPY, SUMPX ) ELSE PHISUM = 0.D0 ENDIF PHIPAR(I) = PHISUM + PI + PI * (NPI+1-I) * (2.*RD(I)-1.) * FNPI PX(I) = COS( PHIPAR(I) ) * PTR(I) PY(I) = SIN( PHIPAR(I) ) * PTR(I) SUMPX = SUMPX + PX(I) SUMPY = SUMPY + PY(I) 29 CONTINUE C CORRECTION OF TRANSVERSE MOMENTA TO KEEP TRANSVERSE MOMENTUM BALANCE SUMPT2 = ECM DPX = SUMPX * FNPI DPY = SUMPY * FNPI DO 30 I = 1,NPI PX(I) = PX(I) - DPX PY(I) = PY(I) - DPY IF ( PX(I) .NE. 0.D0 .OR. PY(I) .NE. 0.D0 ) THEN PHIPAR(I) = ATAN2( PY(I), PX(I) ) ELSE PHIPAR(I) = 0.D0 ENDIF PTSQ(I) = PX(I)**2 + PY(I)**2 SUMPT2 = SUMPT2 - SQRT( PAMA(NTYP(I))**2 + PTSQ(I) ) 30 CONTINUE C CHECK, IF C.M. ENERGY IS EXHAUSTED BY TRANSVERSE MOMENTA C IF SO, TRY ANOTHER SET OF TRANSVERSE MOMENTA IF ( SUMPT2 .LE. 0.D0 ) GOTO 27 C DISTRIBUTION OF LONGITUDINAL MOMENTA PL C SUM1PL IS SUM OF ABS. VALUES OF LONGITUDINAL MOMENTA F = SUMPT2 * FNPI SUM1PL = 0.D0 DO 31 I = 1,NPI FWHM = F + 0.5D0 * SQRT( PTSQ(I) ) PL(I) = ABS( RANNOR(0.D0,FWHM) ) SUM1PL = SUM1PL + PL(I) C SET SEQUENCE COUNTER ISEQ(I) = I 31 CONTINUE C SORT ISEQ IN DECREASING SIZE OF THE LONGITUDINAL MOMENTUM DO 33 I = 1,NPI DO 32 K = I+1,NPI IF ( PL(ISEQ(I)) .LT. PL(ISEQ(K)) ) THEN IHELP = ISEQ(I) ISEQ(I) = ISEQ(K) ISEQ(K) = IHELP ENDIF 32 CONTINUE 33 CONTINUE C TRY TO BALANCE LONGITUDINAL MOMENTA (TO MINIMIZE CORRECTIONS) C START WITH LONG. MOMENTA IN FORWARD/BACKWARD DIRECTION BY RANDOM CALL RMMAR( RD,1,1 ) IF ( RD(1) .LE. 0.5 ) PL(ISEQ(1)) = - PL(ISEQ(1)) SUMPL = PL(ISEQ(1)) DO 34 I = 2,NPI SUM1PL = SUM1PL - PL(ISEQ(I)) C IF THERE IS NOT ENOUGH MOMENTUM LEFT, SELECT FORWARD/BACKWARD TO C BALANCE MOMENTUM, ELSE CHOOSE DIRECTION BY RANDOM IF ( PL(ISEQ(I))+ABS(SUMPL) .GT. SUM1PL ) THEN IF ( PL(ISEQ(I))*SUMPL .GT. 0.D0 ) PL(ISEQ(I)) = -PL(ISEQ(I)) ELSE CALL RMMAR( RD,1,1 ) IF ( RD(1) .LE. 0.5 ) PL(ISEQ(I)) = - PL(ISEQ(I)) ENDIF SUMPL = SUMPL + PL(ISEQ(I)) 34 CONTINUE C CORRECTION OF LONGITUDINAL MOMENTA TO KEEP MOMENTUM BALANCE DPL = SUMPL * FNPI DO 35 I = 1,NPI PL(I) = PL(I) - DPL 35 CONTINUE C ITERATIVE CORRECTION OF ALL MOMENTA TO KEEP ENERGY BALANCE IREPET = 0 36 CONTINUE IREPET = IREPET + 1 IF ( IREPET .GT. 10 ) GOTO 27 ETOT = 0.D0 C CHECK ENERGY CONSERVATION DO 37 I = 1,NPI PTSQ(I) = PX(I)**2 + PY(I)**2 E(I) = SQRT( PAMA(NTYP(I))**2 + PTSQ(I) + PL(I)**2 ) ETOT = ETOT + E(I) 37 CONTINUE ECORR = ECM / ETOT - 1.D0 C LOOK WHETHER ENERGY IS CONSERVED WITHIN 1 % IF ( ABS(ECORR) .GT. .01D0 ) THEN C FACTOR IS MODIFIED WITH EMPIRICAL TERM 1/GCM FOR FASTER CONVERGENCE FACT = (0.5D0+GCMI) * ECORR * 0.02D0 * NPI DO 38 I = 1,NPI PX(I) = PX(I) * ( FACT + 1.D0 ) PY(I) = PY(I) * ( FACT + 1.D0 ) PL(I) = PL(I) * ( FACT * 20.D0 + 1.D0 ) 38 CONTINUE GOTO 36 ENDIF C LORENTZ TRANSFORMATION FROM C.M. TO LAB. FRAME DO 40 K = 5,8 SECPAR(K) = CURPAR(K) 40 CONTINUE DO 41 I = 1,NPI PLLAB = GCM * ( PL(I) + BETACM * E(I) ) IF ( PLLAB .LE. 0.D0 ) GOTO 41 CTHETA = PLLAB / SQRT( PTSQ(I) + PLLAB**2 ) IF ( CTHETA .LT. C(27) ) GOTO 41 CALL ADDANG( COSTHE,PHI, CTHETA,PHIPAR(I), SECPAR(3),SECPAR(4) ) IF ( SECPAR(3) .LT. C(29) ) GOTO 41 SECPAR(1) = NTYP(I) SECPAR(2) = GCM / PAMA(NTYP(I)) * ( PL(I) * BETACM + E(I) ) CALL TSTACK 41 CONTINUE MSMM = MSMM + NPI 999 CONTINUE RETURN END