SUBROUTINE NUCINT C----------------------------------------------------------------------- C NUC(LEAR) INT(ERACTION) C C SELECTS TYPE OF INTERACTION PROCESS FOR ISOBAR MODEL, ACCORDING TO ECM C ISOBAR MASSES INDEPENDENT OF RESPECTIVE ENERGY RANGES C HEAVY PRIMARIES AND STRANGE BARYONS INCLUDED C THIS SUBROUTINE IS CALLED FROM MAIN C----------------------------------------------------------------------- IMPLICIT NONE *KEEP,AIR. COMMON /AIR/ COMPOS,PROBTA,AVERAW,AVOGAD DOUBLE PRECISION COMPOS(3),PROBTA(3),AVERAW,AVOGAD *KEEP,CONST. COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER *KEEP,GENER. COMMON /GENER/ GEN,ALEVEL DOUBLE PRECISION GEN,ALEVEL *KEEP,IRET. COMMON /IRET/ IRET1,IRET2 INTEGER IRET1,IRET2 *KEEP,KAONS. COMMON /KAONS/ CKA DOUBLE PRECISION CKA(80) *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,POLAR. COMMON /POLAR/ POLART,POLARF DOUBLE PRECISION POLART,POLARF *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,RANGE. COMMON /RANGE/ CC DOUBLE PRECISION CC(20) *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,SIGM. COMMON /SIGM/ SIGMA,SIGANN,SIGAIR,FRACTN,FRCTNO DOUBLE PRECISION SIGMA,SIGANN,SIGAIR,FRACTN,FRCTNO *KEEP,STATI. COMMON /STATI/ SABIN,SBBIN,INBIN,IPBIN,IKBIN,IHBIN DOUBLE PRECISION SABIN(37),SBBIN(37) INTEGER INBIN(37),IPBIN(37),IKBIN(37),IHBIN(37) *KEEP,VKIN. COMMON /VKIN/ BETACM DOUBLE PRECISION BETACM *KEND. DOUBLE PRECISION BETA3,COSMU,COSTCM,COSTH3,GAMMA3, * PHIMU,PHI3,WORK1,WORK2 INTEGER I,IGO,KJ C----------------------------------------------------------------------- IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9) 444 FORMAT(' NUCINT: CURPAR=',1P,9E10.3) C SET GENERATION AND LEVEL OF LAST INTERACTION SECPAR( 9) = GEN SECPAR(10) = ALEVEL C RESET POLARIZATION, NOT USED FOR PARTICLES OTHER THAN MUONS YET SECPAR(11) = 0.D0 SECPAR(12) = 0.D0 C CALCULATE KIN. ENERGY BIN EKINL = PAMA(ITYPE) * ( GAMMA - 1.D0 ) IF ( EKINL .GE. .1D0 ) THEN KJ = INT( MIN( 37.D0, 4.D0 + 3.D0*LOG10(EKINL) ) ) ELSE KJ = 1 ENDIF C----------------------------------------------------------------------- C PION INCIDENT IF ( ITYPE .EQ. 8 .OR. ITYPE .EQ. 9 ) THEN IF ( DEBUG ) WRITE(MDEBUG,*) 'NUCINT: PION EKINL=',SNGL(EKINL) IPBIN(KJ) = IPBIN(KJ) + 1 C DECAY OR INTERACTION FOR PIONS ? IF ( FDECAY ) THEN DO 10 I = 5,8 SECPAR(I) = CURPAR(I) 10 CONTINUE C DECAY PI(+,-) ----> MU(+,-) + (ANTI)-NEUTRINO(MU) WORK1 = C(48) * GAMMA WORK2 = C(49) * BETA * WORK1 CALL RMMAR( RD,2,1 ) COSTCM = 2.D0 * RD(1) - 1.D0 GAMMA3 = WORK1 + COSTCM * WORK2 BETA3 = SQRT( 1.D0 - 1.D0 / GAMMA3**2 ) COSTH3 = MIN( 1.D0, ( GAMMA * GAMMA3 - C(48) ) * /( BETA * GAMMA * BETA3 * GAMMA3 ) ) PHI3 = PI2 * RD(2) C MUON / NEUTRINO IS DROPPED CALL ADDANG( COSTHE,PHI, COSTH3,PHI3, COSMU,PHIMU ) IF ( COSMU .GT. C(29) ) THEN C DIRECTION OF PION IN THE MUON CM SYSTEM (= DIRECTION OF POLARIZATION) C SEE: G. BARR ET AL., PHYS. REV. D39 (1989) 3532, EQ. 5 C POLART IS COS OF ANGLE BETWEEN PION AND LABORATORY IN THE MU CM C POLARF IS ANGLE PHI AROUND THE LAB DIRECTION IN THE MU CM C POLART, POLARF ARE WITH RESPECT TO THE MU DIRECTION IN THE LAB SYSTEM POLART = ( 2.D0*PAMA(8)*GAMMA*C(7)/(PAMA(5)*GAMMA3) * - C(7) - 1.D0 ) / ( BETA3 * (1.D0 - C(7)) ) POLARF = PHI3 - PI C PION DIRECTION IS DIRECTION OF POLARIZATION FOR PI+, OPPOSITE FOR PI- IF ( ITYPE .EQ. 9 ) THEN POLART = -POLART POLARF = POLARF + PI ENDIF C GET THE POLARIZATION DIRECTION IN THE MU CM RELATIVE TO THE CORSIKA C COORDINATE SYSTEM CALL ADDANG( COSMU,PHIMU, POLART,POLARF, POLART,POLARF ) C MUON IS WRITTEN TO STACK SECPAR( 1) = CURPAR(1) - 3.D0 SECPAR( 2) = GAMMA3 SECPAR( 3) = COSMU SECPAR( 4) = PHIMU SECPAR(11) = POLART SECPAR(12) = POLARF CALL TSTACK SECPAR(11) = 0.D0 SECPAR(12) = 0.D0 ENDIF IRET1 = 1 RETURN ENDIF C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C PION INTERACTS C CALCULATE GAMMA, BETA AND ENERGY IN CENTER OF MASS ECM = SQRT( C(45) * GAMMA + C(46) ) GCM = (PAMA(ITYPE) * GAMMA + PAMA(14)) / ECM BETACM = SQRT( 1.D0 - 1.D0 / GCM**2 ) C LOW ENERGY HADRONIC INTERACTIONS C USE GHEISHA IF THE CROSS SECTION HAS BEEN CALCULATED FOR GHEISHA IF ( GHEISH ) THEN IF ( GHESIG ) THEN CALL CGHEI ELSE CALL SDPM ENDIF ELSE C DETERMINE TYPE OF INTERACTION FOR PIONS IF ( ECM .GT. CC(8) ) THEN C DUAL PARTON MODEL CALL SDPM ELSEIF ( ECM .GT. CC(7) ) THEN C HEAVY ISOBAR + HEAVY MESON CALL BOX69 ELSEIF ( ECM .GT. CC(6) ) THEN CALL RMMAR( RD,1,1 ) IF ( RD(1) .LE. 0.5 ) THEN C HEAVY ISOBAR + PION CALL BOX68 ELSE C HEAVY MESON + NUCLEON CALL BOX67 ENDIF ELSEIF ( ECM .GT. CC(5) ) THEN C LIGHT ISOBAR + PION CALL BOX66 ELSE C ELASTIC SCATTERING CALL BOX65 ENDIF ENDIF C----------------------------------------------------------------------- C NUCLEON OR ANTINUCLEON INCIDENT ELSEIF ( ITYPE .EQ. 13 .OR. ITYPE .EQ. 14 .OR. * ITYPE .EQ. 15 .OR. ITYPE .EQ. 25 ) THEN C CALCULATE GAMMA, BETA AND ENERGY IN CENTER OF MASS GCM = SQRT( GAMMA * 0.5D0 + 0.5D0 ) ECM = PAMA(ITYPE) * GCM * 2.D0 BETACM = SQRT( 1.D0 - 1.D0 / GCM**2 ) IF ( DEBUG ) WRITE(MDEBUG,*) 'NUCINT: NUCL EKINL=',SNGL(EKINL) INBIN(KJ) = INBIN(KJ) + 1 C LOW ENERGY HADRONIC INTERACTIONS C USE GHEISHA IF THE CROSS SECTION HAS BEEN CALCULATED FOR GHEISHA IF ( GHEISH ) THEN IF ( GHESIG ) THEN CALL CGHEI ELSE CALL SDPM ENDIF ELSE C DETERMINE TYPE OF INTERACTION FOR NUCLEONS AND ANTINUCLEONS IF ( ECM .GT. CC(4) ) THEN C DUAL PARTON MODEL CALL SDPM ELSEIF ( ECM .GT. CC(3) ) THEN C USE THE INTERACTION ROUTINES OF PKF GRIEDER C 2 HEAVY ISOBARS AND ANNIHILATION CALL BOX63 ELSEIF ( ECM .GT. CC(2) ) THEN C 1 HEAVY ISOBAR + NUCLEON AND ANNIHILATION CALL BOX62 ELSEIF ( ECM .GT. CC(1) ) THEN C 1 LIGHT ISOBAR + NUCLEON AND ANNIHILATION CALL BOX61 ELSE C ELASTIC SCATTERING AND ANNIHILATION CALL BOX60 ENDIF ENDIF C----------------------------------------------------------------------- C KAON INCIDENT ELSEIF ( ITYPE .EQ. 11 .OR. ITYPE .EQ. 12 .OR. * ITYPE .EQ. 10 .OR. ITYPE .EQ. 16 ) THEN IF ( DEBUG ) WRITE(MDEBUG,*) 'NUCINT: KAON EKINL=',SNGL(EKINL) IKBIN(KJ) = IKBIN(KJ) + 1 C DECAY OR INTERACTION FOR KAONS ? IF ( FDECAY ) THEN C KAON DECAYS. DETERMINE DECAY MODE FOR KAONS AND SET LIFE TIME IF ( ITYPE .EQ. 10 ) THEN C K(0,L)-MESON IGO = 4 ELSEIF ( ITYPE .EQ. 11 ) THEN C K(+)-MESON IGO = 1 ELSEIF ( ITYPE .EQ. 12 ) THEN C K(-)-MESON IGO = 2 ELSE C K(0,S)-MESON IGO = 3 ENDIF CALL KDECAY( IGO ) RETURN ELSE C- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C KAON INTERACTS C CALCULATE GAMMA, BETA AND ENERGY IN CENTER OF MASS ECM = SQRT( CKA(13) * GAMMA + CKA(14) ) GCM = ( PAMA(ITYPE) * GAMMA + PAMA(14) ) / ECM BETACM = SQRT( 1.D0 - 1.D0 / GCM**2 ) C LOW ENERGY HADRONIC INTERACTIONS C USE GHEISHA IF THE CROSS SECTION HAS BEEN CALCULATED FOR GHEISHA IF ( GHEISH ) THEN IF ( GHESIG ) THEN CALL CGHEI ELSE CALL SDPM ENDIF ELSE C KAON INTERACTS. DETERMINE TYPE OF INTERACTION FOR KAONS IF ( ECM .GT. CC(12) ) THEN C DUAL PARTON MODEL CALL SDPM C USE THE INTERACTION ROUTINES OF PKF GRIEDER ELSEIF ( ECM .GT. CC(11) ) THEN C HEAVY ISOBAR + STRANGE MESON CALL BOX74 ELSEIF ( ECM .GT. CC(10) ) THEN CALL RMMAR( RD,1,1 ) IF ( RD(1) .GT. CKA(21) ) THEN C HEAVY ISOBAR + KAON CALL BOX73 ELSE C STRANGE MESON + NUCLEON CALL BOX72 ENDIF ELSEIF ( ECM .GT. CC(9) ) THEN C LIGHT ISOBAR + KAON CALL BOX71 ELSE C ELASTIC SCATTERING CALL BOX70 ENDIF ENDIF ENDIF C----------------------------------------------------------------------- C STRANGE BARYON (LAMDA, SIGMA) INCIDENT ELSEIF ( (ITYPE .GE. 18 .AND. ITYPE .LE. 24) .OR. * (ITYPE .GE. 26 .AND. ITYPE .LE. 32) ) THEN IF ( DEBUG ) WRITE(MDEBUG,*) 'NUCINT: SBAR EKINL=',SNGL(EKINL) IHBIN(KJ) = IHBIN(KJ) + 1 C DECAY OR INTERACTION FOR STRANGE BARYONS? IF ( FDECAY ) THEN CALL STRDEC RETURN ENDIF C CALCULATE GAMMA, BETA AND ENERGY IN CENTER OF MASS ECM = SQRT( 2.D0 * PAMA(ITYPE) * PAMA(14) * GAMMA * + PAMA(ITYPE)**2 + PAMA(14)**2 ) GCM = ( PAMA(ITYPE) * GAMMA + PAMA(14)) / ECM BETACM = SQRT( 1.D0 - 1.D0 / GCM**2 ) C LOW ENERGY HADRONIC INTERACTIONS C USE GHEISHA IF THE CROSS SECTION HAS BEEN CALCULATED FOR GHEISHA IF ( GHEISH ) THEN IF ( GHESIG ) THEN CALL CGHEI ELSE C VENUS TREATS STRANGE BARYONS CALL SDPM ENDIF ELSE IF ( ECM .GT. CC(4) ) THEN CALL SDPM ELSE C USE THE INTERACTION ROUTINES OF PKF GRIEDER C ELASTIC SCATTERING CALL BOX60 ENDIF ENDIF C----------------------------------------------------------------------- C HEAVY PRIMARY INCIDENT ELSEIF ( ITYPE .GT. 100 ) THEN IF (DEBUG) WRITE(MDEBUG,*) 'NUCINT: HEAVY PRIMARY EKINL=', * SNGL(EKINL) C USE GHEISHA IF THE CROSS SECTION HAS BEEN CALCULATED FOR GHEISHA C (THIS MIGHT BE THE CASE FOR DEUTERONS, TRITONS AND ALPHAS) IF ( GHEISH ) THEN IF ( GHESIG ) THEN CALL CGHEI ELSE CALL SDPM ENDIF ELSE C TREAT HEAVY PRIMARY IN SDPM CALL SDPM ENDIF C----------------------------------------------------------------------- C ILLEGAL PARTICLE ELSE WRITE(MONIOU,*) 'NUCINT: ILLEGAL PARTICLE = ',ITYPE STOP ENDIF C----------------------------------------------------------------------- C KILL PARTICLE IRET1 = 1 RETURN END