#if HAVE_CONFIG_H #include "config.h" #endif #if __FLUKAINFN__ *-- Author : D. HECK IK FZK KARLSRUHE 18/03/2003 C======================================================================= DOUBLE PRECISION FUNCTION FLRNDM() C----------------------------------------------------------------------- C FL(UKA) R(A)ND(O)M (GENERATOR) C SEE SUBROUT. RMMARD C WE USE HERE A SIMPLIFIED FORM OF RMMARD WITH JSEQ=1, LENV=1. C THIS FUNCTON IS CALLED FROM FLUKA ROUTINES. C----------------------------------------------------------------------- IMPLICIT NONE #define __RANMA3INC__ #define __RANMA4INC__ #if __CONEX__ #define __CONEXINC__ #endif #include "corsika.h" #if __CONEX__ #include "conex.h" #endif SAVE C----------------------------------------------------------------------- JSEQ = 1 #if __CONEX__ IF ( FINCNX ) JSEQ = lseq #endif UNI = U(I97(JSEQ),JSEQ) - U(J97(JSEQ),JSEQ) IF ( UNI .LT. 0.D0 ) UNI = UNI + 1.D0 U(I97(JSEQ),JSEQ) = UNI I97(JSEQ) = I97(JSEQ) - 1 IF ( I97(JSEQ) .EQ. 0 ) I97(JSEQ) = 97 J97(JSEQ) = J97(JSEQ) - 1 IF ( J97(JSEQ) .EQ. 0 ) J97(JSEQ) = 97 C(JSEQ) = C(JSEQ) - CD IF ( C(JSEQ) .LT. 0.D0 ) C(JSEQ) = C(JSEQ) + CM UNI = UNI - C(JSEQ) IF ( UNI .LT. 0.D0 ) UNI = UNI + 1.D0 C AN EXACT ZERO HERE IS VERY UNLIKELY, BUT LET''S BE SAFE. IF ( UNI .EQ. 0.D0 ) UNI = TWOM48 FLRNDM = UNI NTOT(JSEQ) = NTOT(JSEQ) + 1 IF ( NTOT(JSEQ) .GE. MODCNS ) THEN NTOT2(JSEQ) = NTOT2(JSEQ) + 1 NTOT(JSEQ) = NTOT(JSEQ) - MODCNS ENDIF RETURN END *-- Author : D. HECK IK FZK KARLSRUHE 22/10/2002 C======================================================================= BLOCK DATA FLUDAT C----------------------------------------------------------------------- C FLU(KA) DAT(A) C C SETS PARTICLE CODE TABLES FOR CONVERIONS BETWEEN CORSIKA AND FLUKA C----------------------------------------------------------------------- IMPLICIT NONE #define __FLULININC__ #include "corsika.h" C ICTABL CONVERTS CORSIKA PARTICLES INTO FLUKA PARTICLES C FIRST TABLE ONLY IF CHARMED PARTICLES CAN BE TREATED C DATA ICFTABL/ C * 7, 4, 3, 0, 10, 11, 23, 13, 14, 12, ! 10 C * 15, 16, 8, 1, 2, 19, 0, 17, 21, 22, ! 20 C * 20, 34, 36, 38, 9, 18, 31, 32, 33, 34, ! 30 C * 37, 39, 8*0, C * 10*0, ! 50 C * 10*0, C * 0, 0, 0, 0, 0, 5, 6, 27, 28, 0, ! 70 C * 10*0, C * 10*0, C * 10*0, !100 C * 10*0, C * 0, 0, 0, 0, 0, 47, 45, 46, 48, 49, !120 C * 50, 0, 0, 0, 0, 0, 0, 0, 0, 0, !130 C * 41, 42, 43, 44, 0, 0, 51, 52, 53, 0, !140 C * 0, 0, 54, 55, 56, 0, 0, 0, 57, 58, !150 C * 59, 0, 0, 0, 60, 61, 62, 0, 0, 0, !160 C * 40*0/ C MODIFIED TABLE AS CHARMED PARTICLES CANNOT BE TREATED BY FLUKA C REPLACE CHARMED QUARK BY STRANGE QUARK DATA ICFTABL/ * 7, 4, 3, 0, 10, 11, 23, 13, 14, 12, ! 10 * 15, 16, 8, 1, 2, 19, 0, 17, 21, 22, ! 20 * 20, 34, 36, 38, 9, 18, 31, 32, 33, 34, ! 30 * 37, 39, 8*0, * 10*0, ! 50 * 10*0, * 0, 0, 0, 0, 0, 5, 6, 27, 28, 0, ! 70 * 10*0, * 10*0, * 10*0, !100 * 10*0, * 0, 0, 0, 0, 0, 16, 25, 24, 15, 0, !120 * 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, !130 * 41, 42, 43, 44, 0, 0, 17, 34, 36, 21, !140 * 22, 20, 34, 36, 38, 0, 0, 0, 18, 35, !150 * 37, 31, 32, 33, 35, 37, 39, 0, 0, 0, !160 * 40*0/ C IFCTABL CONVERTS FLUKA PARTICLES INTO CORSIKA PARTICLES DATA IFCTABL/ * 402, 302, 301, 201, 0, 0, 0, * 14, 15, 3, 2, 66, 67, 1, 13, 25, 5, * 6, 10, 8, 9, 11, 12, 18, 26, 16, 21, * 19, 20, 7, 0, 0, 0, 68, 69, 0, 0, * 27, 28, 29, 22, 30, 23, 31, 24, 32, 0, * 131, 132, 133, 134, 117, 118, 116, 119, 120, 121, * 137, 138, 139, 143, 144, 145, 149, 150, 151, 155, * 156, 157, 0, 0, 36*0/ END *-- Author : D. HECK IK FZK KARLSRUHE 22/10/2002 *-- New version: A. FERRARI IAP KARLSRUHE 28/10/2022 C======================================================================= SUBROUTINE FLUINI C----------------------------------------------------------------------- C FLU(KA) INI(TIALIZATION) C C INITIALIZES FLUKA 2021.2 C THIS SUBROUTINE IS CALLED FROM START. C----------------------------------------------------------------------- #if __LINUX__ INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' #else INCLUDE 'DBLPRC' INCLUDE 'DIMPAR' INCLUDE 'IOUNIT' #endif INTEGER MXCKMT, MXCKEL, NMATFL, NELMFL, MTFLKA, IZELFL, & MXELFL, IFLXYZ LOGICAL LPRINT DOUBLE PRECISION WFELFL, PPTMAX, EF2DP3, DF2DP3 CHARACTER CRVRCK*8 * Maximum number of Fluka materials (change as you wish) PARAMETER ( MXCKMT = 5 ) * Cumulative maximum number of elements in the Fluka materials PARAMETER ( MXCKEL = 10 ) DIMENSION NELMFL (MXCKMT), MTFLKA (MXCKMT), IZELFL (MXCKEL), & WFELFL (MXCKEL) #define __AIRINC__ #define __RUNPARINC__ #include "corsika.h" #if !__GFORTRAN__ SAVE #endif C----------------------------------------------------------------------- IF ( DEBUG ) WRITE(MDEBUG,*) 'FLUINI:' * Define three materials NMATFL = 3 * Air: NELMFL (1) = 3 IZELFL (1) = 7 IZELFL (2) = 8 IZELFL (3) =18 * Weight fractions: WFELFL (1) = 0.92561D-03 WFELFL (2) = 0.28361D-03 WFELFL (3) = 0.15776D-04 * Iron NELMFL (2) = 1 IZELFL (4) = 26 * Weight fraction: WFELFL (4) = 1.D+00 * Water: NELMFL (3) = 2 IZELFL (5) = 1 IZELFL (6) = 8 * Weight fractions: WFELFL (5) = 0.11190D+00 WFELFL (6) = 0.88810D+00 * Total number of components for all materials: MXELFL = 6 * Maximum momentum (GeV/c) for Fluka initializations PPTMAX = 1.D+11 * Transition energy between Fluka (Peanut) and Dpmjet3 for hA * interactions (0 leaves at the Fluka default ~ 10 TeV) EF2DP3 = 0.D+00 * Smearing (GeV/c) of the transition energy between Fluka (Peanut) * and Dpmjet3 for hA interactions (-1 leaves at the Fluka default, * 0 is accepted, as well as positive values) DF2DP3 =-1.D+00 * Iflxyz must be consistent in all calls, to Stpxyz, Sgmxyz, Evtxyz * Iflxyz = 1 -> only inelastic * Iflxyz = 10 -> only elastic * Iflxyz = 11 -> inelastic + elastic * Iflxyz =100 -> only emd * Iflxyz =101 -> inelastic + emd * Iflxyz =110 -> elastic + emd * Iflxyz =111 -> inelastic + elastic + emd IFLXYZ = 11 LPRINT = .TRUE. * Simple check to verify the `right' Fluka is linked: ICHECK = NINT ( AMPRMU * 1.D+12 - 1.0072D+12 ) WRITE ( CRVRCK, '(I8)' ) ICHECK * End check CALL STPXYZ ( NMATFL, NELMFL, IZELFL, WFELFL, MXELFL, & PPTMAX, EF2DP3, DF2DP3, IFLXYZ, LPRINT, & MTFLKA, CRVRCK ) * Simple check to verify the `right' Fluka is linked: IF ( CRVRCK .NE. CRVRFL ) & CALL FLABRT ( 'CORSIKA INI', 'MISMATCHED FLUKA VERSION' ) * End check * Initialization of rm48 is done in start IF ( DEBUG ) WRITE(MDEBUG,*)'FLUINI: INITIALIZATION DONE' RETURN END *-- Author : D. HECK IK FZK KARLSRUHE 22/10/2002 *-- New version: A. FERRARI IAP KARLSRUHE 28/10/2022 C======================================================================= SUBROUTINE FLULNK C----------------------------------------------------------------------- C FLU(KA) L(I)NK ROUTINE C C LINKING SUBROUTINE TO FLUKA 2021.2 C THIS SUBROUTINE IS CALLED FROM SDPM. C----------------------------------------------------------------------- #if __LINUX__ INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' * INCLUDE '(EVTFLG)' INCLUDE '(PAPROP)' INCLUDE '(RESNUC)' #else INCLUDE 'DBLPRC' INCLUDE 'DIMPAR' INCLUDE 'IOUNIT' * INCLUDE 'EVTFLG' INCLUDE 'PAPROP' INCLUDE 'RESNUC' #endif #define __FLULININC__ #define __INTERINC__ #define __PAMINC__ #define __PARPARINC__ #define __PARPAEINC__ #define __RANDPAINC__ #define __RESTINC__ #define __RUNPARINC__ #define __SIGMINC__ #include "corsika.h" DOUBLE PRECISION ETOTAL, CUMSGI, CUMSGE, CUMSGM DIMENSION CUMSGI (0:NELEMX), CUMSGE (0:NELEMX), CUMSGM (0:NELEMX) INTEGER IRAND(3),L #if !__GFORTRAN__ SAVE #endif C----------------------------------------------------------------------- IF ( DEBUG ) WRITE(MDEBUG,*) 'FLULNK: ' C CONVERT PARTICLE TYPE TO FLUKA KPROJ = ICFTABL(ITYPE) IF ( KPROJ .EQ. 0 ) THEN WRITE(MONIOU,*) 'FLULNK: ILLEGAL PARTICLE TYPE (CORSIKA):', * ITYPE RETURN ENDIF * Calculate kinetic energy ekin IF ( PAMA(ITYPE) .NE. 0.D0 ) THEN ETOTAL = CURPAR(1) * PAMA(ITYPE) ELSE ETOTAL = CURPAR(1) ENDIF ELAB = ETOTAL EKIN = ETOTAL - PAMA(ITYPE) * Calculate the momentum pproj consistent with Fluka masses PPROJ = SQRT( EKIN * ( EKIN + TWOTWO * AM (KPROJ) ) ) IF ( DEBUG ) THEN WRITE(MDEBUG,*) 'FLULNK: KPROJ,EKIN=',KPROJ,EKIN C RANDOM GENERATOR STATUS (SEQUENCE L=1) AT BEGINNING OF EVENT CALL RMMAQD( ISEED(1,1),1,'R' ) IRAND(1) = ISEED(1,1) C NUMBER OF CALLS IRAND(2) = ISEED(2,1) C NUMBER OF BILLIONS IRAND(3) = ISEED(3,1) WRITE(MDEBUG,158) (IRAND(J),J=1,3) 158 FORMAT(' FLULNK: RANDOM NUMBER GENERATOR AT BEGIN:' * ,' SEQUENCE= 1 SEED= ',I9,' CALLS=',I9, * ' BILLIONS=',I9) ENDIF * Use fluka material number 3 (= Air) MMAT = 3 * Use the Fluka-internally used direction cosines: TXX = CURPAR (3) TYY = CURPAR (4) TZZ = CURPAR (2) WEE = ONEONE IJ = KPROJ POO = PPROJ * Iflxyz must be consistent in all calls, to Stpxyz, Sgmxyz, Evtxyz * Iflxyz = 1 -> only inelastic * Iflxyz = 10 -> only elastic * Iflxyz = 11 -> inelastic + elastic * Iflxyz =100 -> only emd * Iflxyz =101 -> inelastic + emd * Iflxyz =110 -> elastic + emd * Iflxyz =111 -> inelastic + elastic + emd IFLXYZ = 11 * Output variables: * * Cumsgi(i) = cumulative macroscopic inelastic cross section * (cm^2/g x rho) for the i_th element of the mmat * material * Cumsge(i) = cumulative macroscopic elastic cross section * (cm^2/g x rho) for the i_th element of the mmat * material * Cumsgm(i) = cumulative macroscopic emd cross section * (cm^2/g x rho) for the i_th element of the mmat * material CALL EVTXYZ ( IJ , MMAT , EKIN , POO , TXX , TYY , & TZZ , IFLXYZ, CUMSGI, CUMSGE, CUMSGM ) IF ( DEBUG ) THEN WRITE(MDEBUG,*) 'FLULNK: A_TARGET=',IBTAR, * ' Z_TARGET=',ICHTAR IF ( LINEVT ) THEN WRITE(MDEBUG,*)' INELASTIC INTERACTION' ELSE IF ( LELEVT ) THEN WRITE(MDEBUG,*)' ELASTIC INTERACTION' ELSE IF ( LELDIS ) THEN WRITE(MDEBUG,*)' EMD INTERACTION' END IF END IF IF ( ICHTAR .EQ. 7 ) THEN TAR = 14.D0 ELSEIF ( ICHTAR .EQ. 8 ) THEN TAR = 16.D0 ELSEIF ( ICHTAR .EQ. 18 ) THEN TAR = 40.D0 ENDIF IF ( DEBUG ) WRITE(MDEBUG,*) 'FLULNK: NOW STORE PARTICLES' * Store the resulting particles (to be changed for Corsika8) CALL FLUSTR RETURN END *-- Author : D. HECK IK FZK KARLSRUHE 22/10/2002 *-- New version: A. FERRARI IAP KARLSRUHE 28/10/2022 C======================================================================= SUBROUTINE FLUSIG( EKIN0,PLAB ) C----------------------------------------------------------------------- C FLU(KA) SIG(MA) C C DETERMINES THE INTERACTION CROSS SECTION USED BY FLUKA 2021.2 C THIS SUBROUTINE IS CALLED FROM BOX2. C ARGUMENTS: C EKIN0 = KINETIC ENERGY OF PARTICLE (GEV) C PLAB = MOMENTUM OF PARTICLE (GEV) C----------------------------------------------------------------------- #if __LINUX__ INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(PAPROP)' #else INCLUDE 'DBLPRC' INCLUDE 'DIMPAR' INCLUDE 'IOUNIT' INCLUDE 'PAPROP' #endif #define __AIRINC__ #define __FLULININC__ #define __PAMINC__ #define __PARPARINC__ #define __RUNPARINC__ #define __SIGMINC__ #include "corsika.h" DOUBLE PRECISION EKIN0,EKIN,PLAB INTEGER MMAT #if !__GFORTRAN__ SAVE #endif C----------------------------------------------------------------------- IF ( DEBUG ) WRITE(MDEBUG,*) 'FLUSIG: EKIN=',EKIN0 EKIN = EKIN0 * Convert particle type to Fluka KPROJ = ICFTABL(ITYPE) IF ( KPROJ .EQ. 0 ) THEN WRITE(MONIOU,*) 'FLUSIG: ILLEGAL PARTICLE TYPE (CORSIKA):', * ITYPE SIGAIR = 0.D0 RETURN ENDIF * PPROJ = SQRT( EKIN * (EKIN + PAMA(ITYPE)*2.D0) ) * PPROJ = PLAB * Be sure the momentum is consistent with Fluka masses PPROJ = SQRT ( EKIN * ( EKIN + TWOTWO * AM (KPROJ) ) ) * Material number 3 is Air MMAT = 3 * Get the cross section: sigrea is microscopic cross section (mb) * Iflxyz must be consistent in all calls, to Stpxyz, Sgmxyz, Evtxyz * Iflxyz = 1 -> only inelastic * Iflxyz = 10 -> only elastic * Iflxyz = 11 -> inelastic + elastic * Iflxyz =100 -> only emd * Iflxyz =101 -> inelastic + emd * Iflxyz =110 -> elastic + emd * Iflxyz =111 -> inelastic + elastic + emd IFLXYZ = 11 SIGREA = SGMXYZ ( KPROJ, MMAT, EKIN, PPROJ, IFLXYZ ) SIGMA = 0.D0 SIGAIR = SIGREA IF ( DEBUG ) WRITE(MDEBUG,*) 'FLUSIG: SIGAIR=',SIGAIR RETURN END *-- Author : D. HECK IK FZK KARLSRUHE 22/10/2002 C======================================================================= SUBROUTINE FLUSTR C----------------------------------------------------------------------- C FLU(KA) ST(O)R(E) C C STORES THE SECONDARY PARTICLES OF FLUKA INTO CORSIKA. C FOR FLUKA 2021.2 C THIS SUBROUTINE IS CALLED FROM FLULNK. C----------------------------------------------------------------------- #if __LINUX__ INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' INCLUDE '(FLKCMP)' INCLUDE '(FHEAVY)' INCLUDE '(GENSTK)' INCLUDE '(FLKMAT)' INCLUDE '(NUCDAT)' INCLUDE '(NUCGEO)' INCLUDE '(PAPROP)' INCLUDE '(PAREVT)' INCLUDE '(PART2)' INCLUDE '(PHNCCM)' INCLUDE '(RESNUC)' #else INCLUDE 'DBLPRC' INCLUDE 'DIMPAR' INCLUDE 'IOUNIT' INCLUDE 'FLKCMP' INCLUDE 'FHEAVY' INCLUDE 'GENSTK' INCLUDE 'FLKMAT' INCLUDE 'NUCDAT' INCLUDE 'NUCGEO' INCLUDE 'PAPROP' INCLUDE 'PAREVT' INCLUDE 'PART2' INCLUDE 'PHNCCM' INCLUDE 'RESNUC' #endif #define __ELADPMINC__ #define __ELASTYINC__ #define __FLULININC__ #define __INTERINC__ #define __ISTAINC__ #define __LONGIINC__ #define __MULTINC__ #define __PAMINC__ #define __PARPARINC__ #define __PARPAEINC__ #define __RANDPAINC__ #define __RESTINC__ #define __RUNPARINC__ #define __SIGMINC__ #if __AUGERHIST__ || __EHISTORY__ #define __GENERINC__ #endif #if __AUGERHIST__ || __COASTUSERLIB__ #define __OBSPARINC__ #endif #include "corsika.h" DOUBLE PRECISION ELASTI,EMAX,ENP,ETOT,FAC1,FAC2,PTOTRES INTEGER I,J,KODCRS,LL #if __EHISTORY__ INTEGER IK #endif #if __AUGERHIST__ DOUBLE PRECISION EDEP,THICKLOC,THICK INTEGER II EXTERNAL THICK #endif #if __COASTUSERLIB__ c definition of the COAST crs::CInteraction class COMMON/coastInteraction/coastX, coastY, coastZ, & coastE, coastCX, coastEl, coastProjId, coastTargId, & coastT double precision coastX, coastY, coastZ double precision coastE, coastCX, coastEl double precision coastT integer coastProjId, coastTargId #endif #if !__GFORTRAN__ SAVE #endif C----------------------------------------------------------------------- IF ( DEBUG ) WRITE(MDEBUG,*) 'FLUSTR: ELAB=',ELAB IF ( DEBUG ) THEN WRITE(MDEBUG,*) 'FLUSTR:',NP,' HADRONIC SECONDARIES' WRITE(MDEBUG,*) 'NO TYPE EKIN CX CY CZ' DO I = 1, NP WRITE(MDEBUG,1010) I,KPART(I),TKI(I),CXR(I),CYR(I),CZR(I) 1010 FORMAT(I4,I4,(1P,E10.2),0P,3(1X,F8.5)) ENDDO ENDIF EMAX = 0.D0 ELASTI = 1.D0 ETOT = 0.D0 #if __EHISTORY__ C COPY PARTICLE INFORMATION, LATER TO BECOME GRANDMOTHER PARTICLE DO IK = 0, 6 SECPAR(28+IK) = CURPAR(IK) ENDDO C STORE GENERATION COUNTER OF MOTHER SECPAR(35) = GEN C STORE MASS PENETRATION BEFORE INTERACTION SECPAR(36) = CURPAR(9) #if __THIN__ SECPAR(37) = CURPAR(13) #endif #endif C LOOP OVER ALL PARTICLES IN COMMON /GENSTK/ DO 1001 J = 1, NP C TAKE K0SHORT OR K0LONG AT RANDOM INSTEAD OF K0 OR K0_BAR IF ( KPART(J) .EQ. 24 .OR. KPART(J) .EQ. 25 ) THEN CALL RMMARD( RD,1,1 ) IF ( RD(1) .GE. 0.5D0 ) THEN KPART(J) = 19 ELSE KPART(J) = 12 ENDIF ENDIF C CONVERT FLUKA PARTICLE CODE INTO CORSIKA PARTICLE CODE KODCRS = IFCTABL(KPART(J)) C CHECK PARTICLE CODE IF ( KODCRS .EQ. 0 ) THEN WRITE(MONIOU,*) 'FLUSTR: UNKNOWN PARTICLE NR.',J, * ' WITH FLUKA CODE =', KPART(J) GOTO 1001 ENDIF #if !__NEUTRINO__ C SKIP NEUTRINOS IF ( ( KODCRS .GE. 66 .AND. KODCRS .LE. 69 ) #if __CHARM__ || __TAULEP__ * .OR. KODCRS .EQ. 133 .OR. KODCRS .EQ. 134 #endif * ) THEN IF ( LLONGI ) THEN C FILL NEUTRINO ENERGY INTO LONGITUDINAL ENERGY DEPOSIT #if __THIN__ DLONG(LHEIGH,8) = DLONG(LHEIGH,8) + TKI(J) * WEIGHT #else DLONG(LHEIGH,8) = DLONG(LHEIGH,8) + TKI(J) #endif ENDIF GOTO 1001 ENDIF #endif C CALCULATE GAMMA FACTOR RSP. ENERGY IF ( PAMA(KODCRS) .NE. 0.D0 ) THEN ENP = TKI(J) + PAMA(KODCRS) SECPAR(1) = ENP / PAMA(KODCRS) ELSE ENP = TKI(J) SECPAR(1) = TKI(J) ENDIF ETOT = ETOT + ENP IF ( KODCRS .GT. 1 .AND. KODCRS .LE. 65 ) THEN IF ( ENP .GT. EMAX ) THEN EMAX = ENP C CALCULATE ELASTICITY FROM MOST ENERGETIC PARTICLE (LEADER) ELASTI = EMAX / ELAB ENDIF ENDIF C COUNTER FOR ENERGY-MULTIPLICITY MATRIX MSMM = MSMM + 1 SECPAR(0) = KODCRS C CALCULATE EMISSION ANGLES SECPAR(2) = CZR(J) SECPAR(3) = CXR(J) SECPAR(4) = CYR(J) #if __UPWARD__ IF ( SECPAR(2) .GE. C(29) ) THEN #else IF ( SECPAR(2) .GT. C(29) ) THEN #endif #if __EHISTORY__ C COPY PARTICLE INFORMATION, LATER TO BECOME MOTHER PARTICLE DO IK = 0, 8 SECPAR(17+IK) = SECPAR(IK) ENDDO #if __THIN__ SECPAR(26) = SECPAR(13) #endif #endif CALL TSTACK ELSE IF ( LLONGI ) THEN C ADD ENERGY TO LONGITUDINAL ENERGY DEPOSIT IF ( KODCRS .EQ. 1 ) THEN #if __THIN__ DLONG(LHEIGH,11) = DLONG(LHEIGH,11) + SECPAR(1)*WEIGHT ELSEIF ( KODCRS .EQ. 2 ) THEN DLONG(LHEIGH,13) = DLONG(LHEIGH,13) * + (SECPAR(1)+1.D0)*PAMA(2)*WEIGHT ELSEIF ( KODCRS .EQ. 3 ) THEN DLONG(LHEIGH,13) = DLONG(LHEIGH,13) * + (SECPAR(1)-1.D0)*PAMA(2)*WEIGHT ELSEIF ( KODCRS .EQ. 5 .OR. KODCRS .EQ. 6 ) THEN DLONG(LHEIGH,15) = DLONG(LHEIGH,15) * + SECPAR(1)*PAMA(5)*WEIGHT ELSE IF ( KODCRS .EQ. 8 .OR. KODCRS .EQ. 9 .OR. * KODCRS .EQ. 11 .OR. KODCRS .EQ. 12 ) THEN FAC1 = 0.25D0 FAC2 = 0.75D0 ELSEIF ( KODCRS .EQ. 10 .OR. KODCRS .EQ. 16 ) THEN FAC1 = 0.5D0 FAC2 = 0.5D0 ELSE FAC1 = 1.D0 FAC2 = 0.D0 ENDIF C ADD TO THE HADRON ENERGY DEPOSIT DLONG(LHEIGH,17) = DLONG(LHEIGH,17) * + ( SECPAR(1)*PAMA(KODCRS)-RESTMS(KODCRS) )*WEIGHT*FAC1 C ADD TO THE NEUTRINO DEPOSIT DLONG(LHEIGH,18) = DLONG(LHEIGH,18) * + ( SECPAR(1)*PAMA(KODCRS)-RESTMS(KODCRS) )*WEIGHT*FAC2 #else DLONG(LHEIGH,11) = DLONG(LHEIGH,11) + SECPAR(1) ELSEIF ( KODCRS .EQ. 2 ) THEN DLONG(LHEIGH,13) = DLONG(LHEIGH,13) * + (SECPAR(1)+1.D0)*PAMA(2) ELSEIF ( KODCRS .EQ. 3 ) THEN DLONG(LHEIGH,13) = DLONG(LHEIGH,13) * + (SECPAR(1)-1.D0)*PAMA(2) ELSEIF ( KODCRS .EQ. 5 .OR. KODCRS .EQ. 6 ) THEN DLONG(LHEIGH,15) = DLONG(LHEIGH,15) * + SECPAR(1) * PAMA(5) ELSE IF ( KODCRS .EQ. 8 .OR. KODCRS .EQ. 9 .OR. * KODCRS .EQ. 11 .OR. KODCRS .EQ. 12 ) THEN FAC1 = 0.25D0 FAC2 = 0.75D0 ELSEIF ( KODCRS .EQ. 10 .OR. KODCRS .EQ. 16 ) THEN FAC1 = 0.5D0 FAC2 = 0.5D0 ELSE FAC1 = 1.D0 FAC2 = 0.D0 ENDIF C ADD TO THE HADRON ENERGY DEPOSIT DLONG(LHEIGH,17) = DLONG(LHEIGH,17) * + ( SECPAR(1)*PAMA(KODCRS)-RESTMS(KODCRS) )*FAC1 C ADD TO THE NEUTRINO DEPOSIT DLONG(LHEIGH,18) = DLONG(LHEIGH,18) * + ( SECPAR(1)*PAMA(KODCRS)-RESTMS(KODCRS) )*FAC2 #endif ENDIF ENDIF #if __AUGERHIST__ THICKLOC = THICK( H ) DO LL = 1, NOBSLV IF ( THICKLOC .GE. THCKOB(LL) .AND. * THICKLOC .LT. THCKOB(LL)+SAMPTH ) THEN C THICKH IS WITHIN 1 G/CM^2 BELOW OBSLEV(LL) C BRING THE ENERGY BELOW ANGULAR CUT TO THE HISTO OF LEVEL LL DO II = 0, 8 OUTPAR(II) = SECPAR(II) ENDDO OUTPAR( 9) = GEN OUTPAR(10) = ALEVEL OUTPAR(13) = WEIGHT IF ( KODCRS .EQ. 1 ) THEN EDEP = OUTPAR(1) * WEIGHT ELSE EDEP = ( OUTPAR(1) * PAMA(KODCRS) * - RESTMS(KODCRS) ) * WEIGHT ENDIF IF ( KODCRS .EQ. 2. .OR. KODCRS .EQ. 3 ) * OUTPAR(1) = OUTPAR(1) * PAMA(2) C WE HAVE ANGULAR CUT CALL AUGERDEPFIL( EDEP,LL,1 ) ELSEIF ( THICKLOC .LT. THCKOB(LL) ) THEN GOTO 113 ENDIF ENDDO 113 CONTINUE #endif ENDIF C COUNTERS FOR FIRST INTERACTION IF ( FIRSTI ) THEN IF ( SECPAR(0).EQ. 7.D0 .OR. SECPAR(0).EQ. 8.D0 * .OR. SECPAR(0).EQ. 9.D0 ) THEN IFINPI = IFINPI + 1 ELSEIF ( SECPAR(0).EQ.13.D0 .OR. SECPAR(0).EQ.14.D0 * .OR. SECPAR(0).EQ.15.D0 .OR. SECPAR(0).EQ.25.D0 ) THEN IFINNU = IFINNU + 1 ELSEIF ( SECPAR(0).EQ.10.D0 .OR. SECPAR(0).EQ.11.D0 * .OR. SECPAR(0).EQ.12.D0 .OR. SECPAR(0).EQ.16.D0 ) THEN IFINKA = IFINKA + 1 ELSEIF ( SECPAR(0).EQ.17.D0 ) THEN IFINET = IFINET + 1 ELSEIF ((SECPAR(0).GE.18.D0 .AND. SECPAR(0).LE.24.D0) * .OR. (SECPAR(0).GE.26.D0 .AND. SECPAR(0).LE.32.D0)) THEN IFINHY = IFINHY + 1 ELSEIF ( SECPAR(0).GE.51.D0 .AND. SECPAR(0).LE.53.D0 ) THEN IFINRHO = IFINRHO + 1 ELSE IFINOT = IFINOT + 1 ENDIF ENDIF 1001 CONTINUE C END OF ORDINARY PARTICLE LOOP C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C NOW TREAT THE FRAGMENTS OF TARGET C CALCULATE INTERACTING TARGET NUCLEONS (INCLUDING EVAPORATION) IWOUNT = 0 IF ( DEBUG ) WRITE(MDEBUG,*) 'FLUSTR: ',NPHEAV,' HEAVY PARTICLES' IF ( NPHEAV .GT. 0 ) THEN IF ( DEBUG ) THEN WRITE(MDEBUG,*) * ' J KHY IBH ICH TKHEAV CX CY CZ' DO J = 1, NPHEAV WRITE(MDEBUG,1020) J,KHEAVY(J),IBHEAV(J),ICHEAV(J), * TKHEAV(J),CXHEAV(J),CYHEAV(J),CZHEAV(J) 1020 FORMAT(1H ,I2,3I4,1P,4(E12.5,1X),0P) ENDDO ENDIF DO J = 1, NPHEAV IF ( KHEAVY(J) .EQ. 1 ) THEN C NEUTRON KODCRS = 13 IWOUNT = IWOUNT + 1 ELSEIF ( KHEAVY(J) .EQ. 2 ) THEN C PROTON KODCRS = 14 IWOUNT = IWOUNT + 1 ELSEIF ( KHEAVY(J) .EQ. 3 ) THEN C DEUTERON KODCRS = 201 IWOUNT = IWOUNT + 2 ELSEIF ( KHEAVY(J) .EQ. 4 ) THEN C TRITON KODCRS = 301 IWOUNT = IWOUNT + 3 ELSEIF ( KHEAVY(J) .EQ. 5 ) THEN C 3-HELIUM KODCRS = 302 IWOUNT = IWOUNT + 3 ELSEIF ( KHEAVY(J) .EQ. 6 ) THEN C 4-HELIUM KODCRS = 402 IWOUNT = IWOUNT + 4 ELSEIF ( KHEAVY(J) .GE. 7 .AND. KHEAVY(J) .LE. 12 ) THEN C SPECIFY HEAVY FRAGMENT BY NUCLEON AND CHARGE NUMBER KODCRS = IBHEAV(KHEAVY(J))*100 + ICHEAV(KHEAVY(J)) C REPLACE SINGLE NUCLEONS BY THE ORDINARY HADRONIC CODE IF ( KODCRS .EQ. 100 ) THEN KODCRS = 13 ELSEIF ( KODCRS .EQ. 101 ) THEN KODCRS = 14 C SKIP ILLEGAL FRAGMENT ELSEIF ( KODCRS .EQ. 0 ) THEN WRITE(MONIOU,*) * 'FLUSTR: WRONG FRAGMENT: KHEAVY,IBHEAV,ICHEAV=', * KHEAVY(J),IBHEAV(KHEAVY(J)),ICHEAV(KHEAVY(J)) GOTO 1030 ENDIF IWOUNT = IWOUNT + IBHEAV(KHEAVY(J)) ELSE WRITE(MONIOU,*) 'FLUSTR: WRONG FRAGMENT = ',KHEAVY(J) GOTO 1030 ENDIF ENP = TKHEAV(J) + PAMA(KODCRS) SECPAR(0) = KODCRS SECPAR(1) = ENP / PAMA(KODCRS) ETOT = ETOT + TKHEAV(J) IF ( KODCRS .GT. 1 .AND. KODCRS .LE. 65 ) THEN IF ( ENP .GT. EMAX ) THEN EMAX = ENP C CALCULATE ELASTICITY FROM MOST ENERGETIC PARTICLE (LEADER) ELASTI = EMAX / ELAB ENDIF ENDIF C CALCULATE EMISSION ANGLES SECPAR(2) = CZHEAV(J) SECPAR(3) = CXHEAV(J) SECPAR(4) = CYHEAV(J) #if __UPWARD__ IF ( SECPAR(2) .GE. C(29) ) THEN #else IF ( SECPAR(2) .GT. C(29) ) THEN #endif #if __EHISTORY__ C COPY PARTICLE INFORMATION, LATER TO BECOME MOTHER PARTICLE DO IK = 0, 8 SECPAR(17+IK) = SECPAR(IK) ENDDO #if __THIN__ SECPAR(26) = SECPAR(13) #endif #endif CALL TSTACK ELSE IF ( LLONGI ) THEN C ADD ENERGY TO LONGITUDINAL ENERGY DEPOSIT #if __THIN__ DLONG(LHEIGH,17) = DLONG(LHEIGH,17) * + ( SECPAR(1)*PAMA(KODCRS)-RESTMS(KODCRS) )*WEIGHT #else DLONG(LHEIGH,17) = DLONG(LHEIGH,17) * + SECPAR(1)*PAMA(KODCRS)-RESTMS(KODCRS) #endif ENDIF #if __AUGERHIST__ THICKLOC = THICK( H ) DO LL = 1, NOBSLV IF ( THICKLOC .GE. THCKOB(LL) .AND. * THICKLOC .LT. THCKOB(LL)+SAMPTH ) THEN C THICKH IS WITHIN 1 G/CM^2 BELOW OBSLEV(LL) C BRING THE ENERGY BELOW ANGULAR CUT TO THE HISTO OF LEVEL LL DO II = 0, 8 OUTPAR(II) = SECPAR(II) ENDDO OUTPAR( 9) = GEN OUTPAR(10) = ALEVEL OUTPAR(13) = WEIGHT EDEP = ( OUTPAR(1) * PAMA(KODCRS) * - RESTMS(KODCRS) ) * WEIGHT C WE HAVE ANGULAR CUT CALL AUGERDEPFIL( EDEP,LL,1 ) ELSEIF ( THICKLOC .LT. THCKOB(LL) ) THEN GOTO 111 ENDIF ENDDO 111 CONTINUE #endif ENDIF C COUNTERS FOR FIRST INTERACTION IF ( FIRSTI ) THEN IF ( SECPAR(0).EQ. 7.D0 .OR. SECPAR(0).EQ. 8.D0 * .OR. SECPAR(0).EQ. 9.D0 ) THEN IFINPI = IFINPI + 1 ELSEIF ( SECPAR(0).EQ.13.D0 .OR. SECPAR(0).EQ.14.D0 * .OR. SECPAR(0).EQ.15.D0 .OR. SECPAR(0).EQ.25.D0 ) THEN IFINNU = IFINNU + 1 ELSEIF ( SECPAR(0).EQ.10.D0 .OR. SECPAR(0).EQ.11.D0 * .OR. SECPAR(0).EQ.12.D0 .OR. SECPAR(0).EQ.16.D0 ) THEN IFINKA = IFINKA + 1 ELSEIF ( SECPAR(0).EQ.17.D0 ) THEN IFINET = IFINET + 1 ELSEIF ((SECPAR(0).GE.18.D0 .AND. SECPAR(0).LE.24.D0) * .OR. (SECPAR(0).GE.26.D0 .AND. SECPAR(0).LE.32.D0)) THEN IFINHY = IFINHY + 1 ELSEIF ( SECPAR(0).GE.51.D0 .AND. SECPAR(0).LE.53.D0 ) THEN IFINRHO = IFINRHO + 1 ELSE IFINOT = IFINOT + 1 ENDIF ENDIF 1030 CONTINUE ENDDO C CALCULATE INTERACTING PROJECTILE NUCLEONS (INCLUDING EVAPORATION) IWOUNP = 1 C END OF HEAVY PARTICLE LOOP ENDIF C - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - C NOW TREAT THE REMAINING TARGET NUCLEUS IF ( DEBUG ) WRITE(MDEBUG,*) 'FLUSTR: REMAINING NUCLEUS' IF ( DEBUG ) WRITE(MDEBUG,*) ' MASS=',IBRES,' CHARGE=',ICRES, * ' ETOT=',EKRES KODCRS = IBRES*100 + ICRES IF ( KODCRS .EQ. 0 ) THEN GOTO 112 C REPLACE SINGLE NUCEONS BY THE ORDINARY HADRONIC CODE ELSEIF ( KODCRS .EQ. 100 ) THEN KODCRS = 13 ELSEIF ( KODCRS .EQ. 101 ) THEN KODCRS = 14 ENDIF ETOT = ETOT + EKRES SECPAR(0) = KODCRS SECPAR(1) = (EKRES+PAMA(KODCRS)) / PAMA(KODCRS) C CALCULATE RECOIL ANGLES PTOTRES = SQRT( PXRES**2 + PYRES**2 + PZRES**2 ) IF ( PTOTRES .EQ. 0.D0 ) GOTO 112 SECPAR(2) = PZRES / PTOTRES #if __UPWARD__ IF ( SECPAR(2) .GE. C(29) ) THEN #else IF ( SECPAR(2) .GT. C(29) ) THEN #endif SECPAR(3) = PXRES / PTOTRES SECPAR(4) = PYRES / PTOTRES #if __EHISTORY__ C COPY PARTICLE INFORMATION, LATER TO BECOME MOTHER PARTICLE DO IK = 0, 8 SECPAR(17+IK) = SECPAR(IK) ENDDO #if __THIN__ SECPAR(26) = SECPAR(13) #endif #endif CALL TSTACK ELSE IF ( LLONGI ) THEN C ADD ENERGY TO LONGITUDINAL ENERGY DEPOSIT #if __THIN__ DLONG(LHEIGH,17) = DLONG(LHEIGH,17) * + ( SECPAR(1)*PAMA(KODCRS) - RESTMS(KODCRS) )*WEIGHT #else DLONG(LHEIGH,17) = DLONG(LHEIGH,17) * + ( SECPAR(1)*PAMA(KODCRS) - RESTMS(KODCRS) ) #endif ENDIF #if __AUGERHIST__ THICKLOC = THICK( H ) DO LL = 1, NOBSLV IF ( THICKLOC .GE. THCKOB(LL) .AND. * THICKLOC .LT. THCKOB(LL)+SAMPTH ) THEN C THICKH IS WITHIN 1 G/CM^2 BELOW OBSLEV(LL) C BRING THE ENERGY BELOW ANGULAR CUT TO THE HISTO OF LEVEL LL DO II = 2, 8 OUTPAR(II) = CURPAR(II) ENDDO OUTPAR( 0) = KODCRS OUTPAR( 1) = ERES / PAMA(KODCRS) OUTPAR( 9) = GEN OUTPAR(10) = ALEVEL OUTPAR(13) = WEIGHT EDEP = EKRES * WEIGHT C WE HAVE ENERGY CUT CALL AUGERDEPFIL( EDEP,LL,2 ) ELSEIF ( THICKLOC .LT. THCKOB(LL) ) THEN GOTO 110 ENDIF ENDDO 110 CONTINUE #endif ENDIF 112 CONTINUE C FILL ELASTICITY IN MATRICES MEL = MIN ( 1.D0+10.D0* MAX( 0.D0, ELASTI ) , 11.D0 ) MEN = MIN ( 4.D0+ 3.D0*LOG10(MAX( .1D0, EKINL )), 40.D0 ) #if __THIN__ IELDPM(MEN,MEL) = IELDPM(MEN,MEL) + NINT( WEIGHT ) IELDPA(MEN,MEL) = IELDPA(MEN,MEL) + NINT( WEIGHT ) IF ( ELASTI .LT. 1.D0 ) THEN ELMEAN(MEN) = ELMEAN(MEN) + ELASTI * WEIGHT ELMEAA(MEN) = ELMEAA(MEN) + ELASTI * WEIGHT #else IELDPM(MEN,MEL) = IELDPM(MEN,MEL) + 1 IELDPA(MEN,MEL) = IELDPA(MEN,MEL) + 1 IF ( ELASTI .LT. 1.D0 ) THEN ELMEAN(MEN) = ELMEAN(MEN) + ELASTI ELMEAA(MEN) = ELMEAA(MEN) + ELASTI #endif ENDIF #if __COASTUSERLIB__ coastProjId = nint(curpar(0)) coastTargId = nint(tar) coastX = curpar(7) coastY = curpar(8) #if __CURVED__ coastZ = curpar(14) #else coastX = coastX - XOFF(NOBSLV) coastY = coastY - YOFF(NOBSLV) coastZ = curpar(5) #endif coastT = curpar(6) coastE = pama(coastProjId)*curpar(1) coastCX = sigair coastEl = elasti call interaction(coastX) #endif IF ( FIRSTI ) THEN TARG1I = NINT( TAR ) SIG1I = SIGAIR ELAST = ELASTI FIRSTI = .FALSE. C RANDOM GENERATOR STATUS (SEQUENCE L=1) AT END OF EVENT LL = 1 CALL RMMAQD( ISEED(1,LL),LL,'R' ) C SEED ISEED1I(1) = ISEED(1,LL) C NUMBER OF CALLS ISEED1I(2) = ISEED(2,LL) C NUMBER OF BILLIONS ISEED1I(3) = ISEED(3,LL) ENDIF IF ( DEBUG ) WRITE(MDEBUG,*) 'FLUSTR: ETOT=',ETOT RETURN END *-- Author : Alfredo Ferrari, INFN - Milan 24/10/2002 C======================================================================= SUBROUTINE ZEREMF C----------------------------------------------------------------------- C THIS DUMMY SUBROUTINE IS NECESSARY TO OVERRIDE A FLUKA SUBROUTINE C WITH IDENTICAL NAME WHICH OTHERWISE WOULD ERASE SOME CORSIKA COMMONS. C THIS SUBROUTINE IS CALLED FROM ZEROIN. C----------------------------------------------------------------------- #if __LINUX__ INCLUDE '(DBLPRC)' INCLUDE '(DIMPAR)' INCLUDE '(IOUNIT)' #else INCLUDE 'DBLPRC' INCLUDE 'DIMPAR' INCLUDE 'IOUNIT' #endif C----------------------------------------------------------------------- RETURN END #endif