SUBROUTINE CGHEI C----------------------------------------------------------------------- C C(ORSIKA) GHE(ISHA) I(NTERFACE) C C MAIN STEERING ROUTINE FOR HADRON PACKAGE GHEISHA *** C THIS SUBROUTINE IS CALLED FROM NUCINT C C ORIGIN : F.CARMINATI, H.FESEFELDT (ROUTINE GHESIG) C REDESIGN: P. GABRIEL IK1 FZK KARLSRUHE C----------------------------------------------------------------------- *KEEP,CGCOMP. PARAMETER (KK=3) COMMON/CGCOMP/ ACOMP,ZCOMP,WCOMP REAL ACOMP(KK),ZCOMP(KK),WCOMP(KK) *KEEP,ELABCT. COMMON /ELABCT/ ELCUT DOUBLE PRECISION ELCUT(4) *KEEP,ELADPM. COMMON /ELADPM/ ELMEAN,ELMEAA,IELDPM,IELDPA DOUBLE PRECISION ELMEAN(37),ELMEAA(37) INTEGER IELDPM(37,13),IELDPA(37,13) *KEEP,ELASTY. COMMON /ELASTY/ ELAST,IELIS,IELHM,IELNU,IELPI DOUBLE PRECISION ELAST INTEGER IELIS(20),IELHM(20),IELNU(20),IELPI(20) *KEEP,GENER. COMMON /GENER/ GEN,ALEVEL DOUBLE PRECISION GEN,ALEVEL *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,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 *KEND. DOUBLE PRECISION ELASTI,ELABOR,PLX,PLY,PLZ,PLSQ,PLTOT,RMASSK COMMON/GSECTI/ AIEL(20),AIIN(20),AIFI(20),AICA(20),ALAM,K0FLAG INTEGER K0FLAG REAL AIEL,AIIN,AIFI,AICA,ALAM C --- GHEISHA COMMONS --- PARAMETER (MXGKGH=100) PARAMETER (MXGKPV=MXGKGH) COMMON /VECUTY/ PV(10,MXGKPV) COMMON/CONSTS/ PI,TWPI,PIBTW,MP,MPI,MMU,MEL,MKCH,MK0,SMP,SMPI, $ SMU,CT,CTKCH,CTK0, $ ML0,MSP,MS0,MSM,MX0,MXM,CTL0,CTSP,CTSM,CTX0,CTXM, $ RMASS(35),RCHARG(35) REAL MP,MPI,MMU,MEL,MKCH,MK0, * ML0,MSP,MS0,MSM,MX0,MXM PARAMETER (MXEVEN=12*MXGKGH) COMMON/EVENT / NSIZE,NCUR,NEXT,NTOT,EVE(MXEVEN) COMMON/PRNTFL/INBCD,NEWBCD,INBIN,NEWBIN,NPEVT,NEVTP,LPRT,NPRT(10) LOGICAL LPRT,NPRT C --- "NEVENT" CHANGED TO "KEVENT" IN COMMON /CURPAR/ DUE TO CLASH --- C --- WITH VARIABLE "NEVENT" IN GEANT COMMON --- PARAMETER (MXGKCU=MXGKGH) COMMON /CURPAR/ WEIGHT(10),DDELTN,IFILE,IRUN,NEVT,KEVENT,SHFLAG, $ ITHST,ITTOT,ITLST,IFRND,TOFCUT,CMOM(5),CENG(5), $ RS,S,ENP(10),NP,NM,NN,NR,NO,NZ,IPA(MXGKCU), $ ATNO2,ZNO2 C --- "IPART" CHANGED TO "KPART" IN COMMON /RESULT/ DUE TO CLASH --- C --- WITH VARIABLE "IPART" IN GEANT COMMON --- COMMON /RESULT/ XEND,YEND,ZEND,RCA,RCE,AMAS,NCH,TOF,PX,PY,PZ, $ USERW,INTCT,P,EN,EK,AMASQ,DELTN,ITK,NTK,KPART,IND, $ LCALO,ICEL,SINL,COSL,SINP,COSP, $ XOLD,YOLD,ZOLD,POLD,PXOLD,PYOLD,PZOLD, $ XSCAT,YSCAT,ZSCAT,PSCAT,PXSCAT,PYSCAT,PZSCAT REAL NCH,INTCT C --- "ABSL(21)" CHANGED TO "ABSLTH(21)" IN COMMON /MAT/ DUE TO CLASH --- C --- WITH VARIABLE "ABSL" IN GEANT COMMON --- COMMON /MAT/ LMAT, $ DEN(21),RADLTH(21),ATNO(21),ZNO(21),ABSLTH(21), $ CDEN(21),MDEN(21),X0DEN(21),X1DEN(21),RION(21), $ MATID(21),MATID1(21,24),PARMAT(21,10), $ IFRAT,IFRAC(21),FRAC1(21,10),DEN1(21,10), $ ATNO1(21,10),ZNO1(21,10) DIMENSION IPELOS(35) REAL EMAX,EEESQ SAVE IDEOL DIMENSION RNDM(1) C --- DIMENSION STMTS. FOR GEANT/GHEISHA PARTICLE CODE CONVERSIONS --- C --- KIPART(I)=GHEISHA CODE CORRESPONDING TO GEANT CODE I --- C --- IKPART(I)=GEANT CODE CORRESPONDING TO GHEISHA CODE I --- DIMENSION KIPART(48),IKPART(35) C --- DATA STMTS. FOR GEANT/GHEISHA PARTICLE CODE CONVERSIONS --- C --- KIPART(I)=GHEISHA CODE CORRESPONDING TO GEANT CODE I --- C --- IKPART(I)=GEANT CODE CORRESPONDING TO GHEISHA CODE I --- DATA KIPART/ $ 1, 3, 4, 2, 5, 6, 8, 7, $ 9, 12, 10, 13, 16, 14, 15, 11, $ 35, 18, 20, 21, 22, 26, 27, 33, $ 17, 19, 23, 24, 25, 28, 29, 34, $ 35, 35, 35, 35, 35, 35, 35, 35, $ 35, 35, 35, 35, 30, 31, 32, 35/ DATA IKPART/ $ 1, 4, 2, 3, 5, 6, 8, 7, $ 9, 11, 16, 10, 12, 14, 15, 13, $ 25, 18, 26, 19, 20, 21, 27, 28, $ 29, 22, 23, 30, 31, 45, 46, 47, $ 24, 32, 48/ C --- DENOTE STABLE PARTICLES ACCORDING TO GHEISHA CODE --- C --- STABLE : GAMMA, NEUTRINO, ELECTRON, PROTON AND HEAVY FRAGMENTS --- C --- WHEN STOPPING THESE PARTICLES ONLY LOOSE THEIR KINETIC ENERGY --- DATA IPELOS/ $ 1, 1, 0, 1, 0, 0, 0, 0, $ 0, 0, 0, 0, 0, 1, 0, 0, $ 0, 0, 0, 0, 0, 0, 0, 0, $ 0, 0, 0, 0, 0, 1, 1, 1, $ 0, 0, 1/ C --- LOWERBOUND OF KINETIC ENERGY BIN IN N CROSS-SECTION TABLES --- DATA TEKLOW /0.0001/ C --- KINETIC ENERGY TO SWITCH FROM "CASN" TO "GNSLWD" FOR N CASCADE --- DATA SWTEKN /0.05/ DATA IDEOL/0/ C----------------------------------------------------------------------- IF ( DEBUG ) WRITE(MDEBUG,*) 'CGHEI :' C --- DEFINE PARTICLE TYPE IF ( ITYPE .LE. 48 ) THEN IPART = ITYPE ELSEIF ( ITYPE .EQ. 201 ) THEN IPART = 45 ELSEIF ( ITYPE .EQ. 301 ) THEN IPART = 46 ELSEIF ( ITYPE .EQ. 402 ) THEN IPART = 47 ELSE WRITE (MONIOU,7795) ITYPE 7795 FORMAT (//,' *CGHEI* ILLEGAL PARTICLE TYPE OCCURS =',I5) IPART = 48 ENDIF NETEST=IKPART(KPART) IF ( NETEST .EQ. IPART ) GO TO 9004 WRITE(MONIOU,8881) IPART,KPART 8881 FORMAT(' *CGHEI* IPART,KPART = ',2(I3,1X)/ $ ' *CGHEI* ======> PARTICLE TYPES DO NOT MATCH <=======') STOP 9004 CONTINUE KPART=KIPART(IPART) KKPART=KPART C --- TRANSPORT THE TRACK NUMBER TO GHEISHA AND INITIALISE SOME NUMBERS C --- NTK=ITRA ITRA = CURRENT TRACK NUMBER IN GEANT (GCKINE) NTK=0 INTCT=0.0 NEXT=1 NTOT=0 INT=0 TOF=0.0 C --- STORE COORDINATES FOR SECONDARIES AND RESET ITYPE SECPAR(1) = 0. DO 7001 LK = 5, 8 SECPAR(LK) = CURPAR(LK) 7001 CONTINUE C --- FILL RESULT COMMON FOR THIS TRACK WITH CORSIKA VALUES --- AMAS=RMASS(KPART) NCH=RCHARG(KPART) 107 XEND = CURPAR(7) YEND = CURPAR(8) ZEND = CURPAR(5) SINL = -CURPAR(3) PHI = CURPAR(4) USERW=0.0 AMASQ=AMAS*AMAS EN = CURPAR(2) * ABS(AMAS) EK = ABS ( EN - ABS(AMAS) ) ENOLD=EN EMAX = 0. P = SQRT ( EN*EN - AMASQ ) ELABOR = EN SINP = SIN(PHI) COSP = COS(PHI) COSL = SQRT ( ABS(1.-SINL**2) ) PX = COSL * COSP PY = COSL * SINP PZ = SINL C --- SET GHEISHA INDEX FOR THE CURRENT MEDIUM ALWAYS TO 1 --- IND=1 C --- TRANSFER GLOBAL MATERIAL CONSTANTS FOR CURRENT MEDIUM --- C --- DETAILED DATA FOR COMPOUNDS IS OBTAINED VIA ROUTINE COMPO --- ATNO(IND+1) = 14.56 ZNO(IND+1) = 7.265 DEN(IND+1) = 0.0 RADLTH(IND+1)= 0.0 ABSLTH(IND+1)= 0.0 C --- SETUP PARMAT FOR PHYSICS STEERING --- PARMAT(IND+1,10)=0.0 5 CONTINUE C --- INDICATE LIGHT (<= PI) AND HEAVY PARTICLES (HISTORICALLY) --- C --- CALIM CODE --- J=2 TEST=RMASS(7)-0.001 IF (ABS(AMAS) .LT. TEST) J=1 C *** DIVISION INTO VARIOUS INTERACTION CHANNELS DENOTED BY "INT" *** C THE CONVENTION FOR "INT" IS THE FOLLOWING C INT = -1 REACTION CROSS SECTIONS NOT YET TABULATED/PROGRAMMED C = 0 NO INTERACTION C = 1 ELEASTIC SCATTERING C = 2 INELASTIC SCATTERING C = 3 NUCLEAR FISSION WITH INELEASTIC SCATTERING C = 4 NEUTRON CAPTURE C INT = 3, 4 SHOULD BE DELETED FOR AIR TARGET C --- INTACT CODE --- ALAM1=0.0 CALL GRNDM(RNDM,1) RAT=RNDM(1)*ALAM ATNO2 = 14.56 ZNO2 = 7.265 DO 6 K=1,KK ATNO2 = ACOMP(K) ZNO2 = ZCOMP(K) C --- TRY FOR ELASTIC SCATTERING --- INT=1 ALAM1=ALAM1+AIEL(K) IF (RAT .LT. ALAM1) GO TO 8 C --- TRY FOR INELASTIC SCATTERING --- INT=2 ALAM1=ALAM1+AIIN(K) IF (RAT .LT. ALAM1) GO TO 8 C --- TRY FOR NEUTRON CAPTURE --- INT=4 ALAM1=ALAM1+AICA(K) IF (RAT .LT. ALAM1) GO TO 8 6 CONTINUE C --- NO REACTION SELECTED ==> ELASTIC SCATTERING --- INT=1 C *** TAKE ACTION ACCORDING TO SELECTED REACTION CHANNEL *** C --- FOLLOWING CODE IS A TRANSLATION OF "CALIM" INTO GEANT JARGON --- 8 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,1001) INT 1001 FORMAT(' *CGHEI* INTERACTION TYPE CHOSEN INT = ',I3) IF (INT .NE. 4) GO TO 10 C --- NEUTRON CAPTURE --- IF (NPRT(9)) WRITE(MDEBUG,2000) 2000 FORMAT(' *CGHEI* ROUTINE CAPTUR WILL BE CALLED') CALL CAPTUR(NOPT) GO TO 40 10 CONTINUE C --- ELASTIC AND INELASTIC SCATTERING --- PV(1,MXGKPV)=P*PX PV(2,MXGKPV)=P*PY PV(3,MXGKPV)=P*PZ PV(4,MXGKPV)=EN PV(5,MXGKPV)=AMAS PV(6,MXGKPV)=NCH PV(7,MXGKPV)=TOF PV(8,MXGKPV)=KPART PV(9,MXGKPV)=0. PV(10,MXGKPV)=USERW C --- ADDITIONAL PARAMETERS TO SIMULATE FERMI MOTION AND EVAPORATION --- DO 111 JENP=1,10 ENP(JENP)=0. 111 CONTINUE ENP(5)=EK ENP(6)=EN ENP(7)=P IF (INT .NE. 1) GO TO 12 C *** ELASTIC SCATTERING PROCESSES *** C --- ONLY NUCLEAR INTERACTIONS FOR HEAVY FRAGMENTS --- IF ((KPART .GE. 30) .AND. (KPART .LE. 32)) GO TO 35 C --- NORMAL ELASTIC SCATTERING FOR LIGHT MEDIA --- IF (ATNO2 .LT. 1.5) GO TO 35 C --- COHERENT ELASTIC SCATTERING FOR HEAVY MEDIA --- IF (NPRT(9)) WRITE(MDEBUG,2002) 2002 FORMAT(' *CGHEI* ROUTINE COSCAT WILL BE CALLED') CALL COSCAT GO TO 40 C *** NON-ELASTIC SCATTERING PROCESSES *** 12 CONTINUE C --- ONLY NUCLEAR INTERACTIONS FOR HEAVY FRAGMENTS --- IF ((KPART .GE. 30) .AND. (KPART .LE. 32)) GO TO 35 C *** USE SOMETIMES NUCLEAR REACTION ROUTINE "NUCREC" FOR LOW ENERGY *** C *** PROTON AND NEUTRON SCATTERING *** CALL GRNDM(RNDM,1) TEST1=RNDM(1) TEST2=4.5*(EK-0.01) IF ((KPART .EQ. 14) .AND. (TEST1 .GT. TEST2)) GO TO 85 IF ((KPART .EQ. 16) .AND. (TEST1 .GT. TEST2)) GO TO 86 C *** FERMI MOTION AND EVAPORATION *** TKIN=CINEMA(EK) PV(9,MXGKPV) = TKIN ENP(5)=EK+TKIN C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES --- IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW ENP(6)=ENP(5)+ABS(AMAS) ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS) ENP(7)=SQRT(ABS(ENP(7))) TKIN=FERMI(ENP(5)) ENP(5)=ENP(5)+TKIN C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES --- IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW ENP(6)=ENP(5)+ABS(AMAS) ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS) ENP(7)=SQRT(ABS(ENP(7))) TKIN=EXNU(ENP(5)) ENP(5)=ENP(5)-TKIN C --- CHECK FOR LOWERBOUND OF EKIN IN CROSS-SECTION TABLES --- IF (ENP(5) .LE. TEKLOW) ENP(5)=TEKLOW ENP(6)=ENP(5)+ABS(AMAS) ENP(7)=(ENP(6)-AMAS)*(ENP(6)+AMAS) ENP(7)=SQRT(ABS(ENP(7))) C *** IN CASE OF ENERGY ABOVE CUT-OFF LET THE PARTICLE CASCADE *** IF ( ENP(5) .GT. ELCUT(1)) GO TO 35 C --- SECOND CHANCE FOR ANTI-BARYONS DUE TO POSSIBLE ANNIHILATION --- IF ((AMAS .GE. 0.0) .OR. (KPART .LE. 14)) GO TO 13 ANNI=1.3*P IF (ANNI .GT. 0.4) ANNI=0.4 CALL GRNDM(RNDM,1) TEST=RNDM(1) IF (TEST .GT. ANNI) GO TO 35 C *** PARTICLE WITH ENERGY BELOW CUT-OFF *** C --- ==> ONLY NUCLEAR EVAPORATION AND QUASI-ELASTIC SCATTERING --- 13 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,1002) KPART,EK,EN,P,ENP(5),ENP(6),ENP(7) 1002 FORMAT(' *CGHEI* ENERGY BELOW CUT-OFF FOR GHEISHA PARTICLE ',I3/ $ ' EK,EN,P,ENP(5),ENP(6),ENP(7) = ',6(G12.5,1X)) IF ((KPART .NE. 14) .AND. (KPART .NE. 16)) GO TO 14 IF (KPART .EQ. 16) GO TO 86 C --- SLOW PROTON --- 85 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2003) EK,KPART 2003 FORMAT(' *CGHEI* ROUTINE NUCREC WILL BE CALLED', $ ' EK = ',G12.5,' GEV KPART = ',I3) CALL NUCREC(NOPT,2) IF (NOPT .NE. 0) GO TO 50 IF (NPRT(9)) WRITE(MDEBUG,2004)EK,KPART 2004 FORMAT(' *CGHEI* ROUTINE COSCAT WILL BE CALLED', $ ' EK = ',G12.5,' GEV KPART = ',I3) CALL COSCAT GO TO 40 C --- SLOW NEUTRON --- 86 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2015) NUCFLG=0 CALL GNSLWD(NUCFLG,INT,NFL,TEKLOW) IF (NUCFLG .NE. 0) GO TO 50 GO TO 40 C --- OTHER SLOW PARTICLES --- 14 CONTINUE IPA(1)=KPART C --- DECIDE FOR PROTON OR NEUTRON TARGET --- IPA(2)=16 CALL GRNDM(RNDM,1) TEST1=RNDM(1) TEST2=ZNO2/ATNO2 IF (TEST1 .LT. TEST2) IPA(2)=14 AVERN=0.0 NFL=1 IF (IPA(2) .EQ. 16) NFL=2 IPPP=KPART IF (NPRT(9)) WRITE(MDEBUG,2005) 2005 FORMAT(' *CGHEI* ROUTINE TWOB WILL BE CALLED') CALL TWOB(IPPP,NFL,AVERN) GOTO 40 C --- INITIALISATION OF CASCADE QUANTITIES --- 35 CONTINUE C *** CASCADE GENERATION *** C --- CALCULATE FINAL STATE MULTIPLICITY AND LONGITUDINAL AND --- C --- TRANSVERSE MOMENTUM DISTRIBUTIONS --- C --- FIXED PARTICLE TYPE TO STEER THE CASCADE --- KKPART=KPART C --- NO CASCADE FOR LEPTONS --- IF (KKPART .LE. 6) GO TO 9999 C *** WHAT TO DO WITH "NEW PARTICLES" FOR GHEISHA ?????? *** C --- RETURN FOR THE TIME BEING --- IF (KKPART .GE. 35) GO TO 9999 C --- CASCADE OF HEAVY FRAGMENTS IF ((KKPART .GE. 30) .AND. (KKPART .LE. 32)) GO TO 390 C --- INITIALIZE THE IPA ARRAY --- CALL VZERO(IPA(1),MXGKCU) C --- CASCADE OF OMEGA - AND OMEGA - BAR --- IF (KKPART .EQ. 33) GO TO 330 IF (KKPART .EQ. 34) GO TO 331 NVEPAR=KKPART-17 IF (NVEPAR .LE. 0) GO TO 15 GO TO (318,319,320,321,322,323,324,325,326,327,328,329),NVEPAR 15 CONTINUE NVEPAR=KKPART-6 GO TO (307,308,309,310,311,312,313,314,315,316,317,318),NVEPAR C --- PI+ CASCADE --- 307 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2006) 2006 FORMAT(' *CGHEI* ROUTINE CASPIP WILL BE CALLED') CALL CASPIP(J,INT,NFL) GO TO 40 C --- PI0 ==> NO CASCADE --- 308 CONTINUE GO TO 40 C --- PI- CASCADE --- 309 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2007) 2007 FORMAT(' *CGHEI* ROUTINE CASPIM WILL BE CALLED') CALL CASPIM(J,INT,NFL) GO TO 40 C --- K+ CASCADE --- 310 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2008) 2008 FORMAT(' *CGHEI* ROUTINE CASKP WILL BE CALLED') CALL CASKP(J,INT,NFL) GO TO 40 C --- K0 CASCADE --- 311 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2009) 2009 FORMAT(' *CGHEI* ROUTINE CASK0 WILL BE CALLED') CALL CASK0(J,INT,NFL) GO TO 40 C --- K0 BAR CASCADE --- 312 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2010) 2010 FORMAT(' *CGHEI* ROUTINE CASK0B WILL BE CALLED') CALL CASK0B(J,INT,NFL) GO TO 40 C --- K- CASCADE --- 313 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2011) 2011 FORMAT(' *CGHEI* ROUTINE CASKM WILL BE CALLED') CALL CASKM(J,INT,NFL) GO TO 40 C --- PROTON CASCADE --- 314 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2012) 2012 FORMAT(' *CGHEI* ROUTINE CASP WILL BE CALLED') CALL CASP(J,INT,NFL) GO TO 40 C --- PROTON BAR CASCADE --- 315 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2013) 2013 FORMAT(' *CGHEI* ROUTINE CASPB WILL BE CALLED') CALL CASPB(J,INT,NFL) GO TO 40 C --- NEUTRON CASCADE --- 316 CONTINUE NUCFLG=0 IF (EK .GT. SWTEKN) THEN CALL CASN(J,INT,NFL) IF (NPRT(9)) WRITE(MDEBUG,2014) 2014 FORMAT(' *CGHEI* ROUTINE CASN WILL BE CALLED') ELSE CALL GNSLWD(NUCFLG,INT,NFL,TEKLOW) IF (NPRT(9)) WRITE(MDEBUG,2015) 2015 FORMAT(' *CGHEI* ROUTINE GNSLWD WILL BE CALLED') ENDIF IF (NUCFLG .NE. 0) GO TO 50 GO TO 40 C --- NEUTRON BAR CASCADE --- 317 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2016) 2016 FORMAT(' *CGHEI* ROUTINE CASNB WILL BE CALLED') CALL CASNB(J,INT,NFL) GO TO 40 C --- LAMBDA CASCADE --- 318 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2017) 2017 FORMAT(' *CGHEI* ROUTINE CASL0 WILL BE CALLED') CALL CASL0(J,INT,NFL) GO TO 40 C --- LAMBDA BAR CASCADE --- 319 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2018) 2018 FORMAT(' *CGHEI* ROUTINE CASAL0 WILL BE CALLED') CALL CASAL0(J,INT,NFL) GO TO 40 C --- SIGMA + CASCADE --- 320 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2019) 2019 FORMAT(' *CGHEI* ROUTINE CASSP WILL BE CALLED') CALL CASSP(J,INT,NFL) GO TO 40 C --- SIGMA 0 ==> NO CASCADE --- 321 CONTINUE GO TO 40 C --- SIGMA - CASCADE --- 322 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2020) 2020 FORMAT(' *CGHEI* ROUTINE CASSM WILL BE CALLED') CALL CASSM(J,INT,NFL) GO TO 40 C --- SIGMA + BAR CASCADE --- 323 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2021) 2021 FORMAT(' *CGHEI* ROUTINE CASASP WILL BE CALLED') CALL CASASP(J,INT,NFL) GO TO 40 C --- SIGMA 0 BAR ==> NO CASCADE --- 324 CONTINUE GO TO 40 C --- SIGMA - BAR CASCADE --- 325 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2022) 2022 FORMAT(' *CGHEI* ROUTINE CASASM WILL BE CALLED') CALL CASASM(J,INT,NFL) GO TO 40 C --- XI 0 CASCADE --- 326 CONTINUE IF (NPRT(9)) PRINT 2023 2023 FORMAT(' *CGHEI* ROUTINE CASX0 WILL BE CALLED') CALL CASX0(J,INT,NFL) GO TO 40 C --- XI - CASCADE --- 327 CONTINUE IF (NPRT(9)) PRINT 2024 2024 FORMAT(' *CGHEI* ROUTINE CASXM WILL BE CALLED') CALL CASXM(J,INT,NFL) GO TO 40 C --- XI 0 BAR CASCADE --- 328 CONTINUE IF (NPRT(9)) PRINT 2025 2025 FORMAT(' *CGHEI* ROUTINE CASAX0 WILL BE CALLED') CALL CASAX0(J,INT,NFL) GO TO 40 C --- XI - BAR CASCADE --- 329 CONTINUE IF (NPRT(9)) PRINT 2026 2026 FORMAT(' *CGHEI* ROUTINE CASAXM WILL BE CALLED') CALL CASAXM(J,INT,NFL) GO TO 40 C --- OMEGA - CASCADE --- 330 CONTINUE IF (NPRT(9)) PRINT 2027 2027 FORMAT(' *CGHEI* ROUTINE CASOM WILL BE CALLED') CALL CASOM(J,INT,NFL) GO TO 40 C --- OMEGA - BAR CASCADE --- 331 CONTINUE IF (NPRT(9)) PRINT 2028 2028 FORMAT(' *CGHEI* ROUTINE CASAOM WILL BE CALLED') CALL CASAOM(J,INT,NFL) GO TO 40 C --- HEAVY FRAGMENT CASCADE --- 390 CONTINUE IF (NPRT(9)) WRITE(MDEBUG,2090) 2090 FORMAT(' *CGHEI* ROUTINE CASFRG WILL BE CALLED') NUCFLG=0 CALL CASFRG(NUCFLG,INT,NFL) IF (NUCFLG .NE. 0) GO TO 50 C *** CHECK WHETHER THERE ARE NEW PARTICLES GENERATED *** 40 CONTINUE IF ((NTOT .NE. 0) .OR. (KKPART .NE. KPART)) GO TO 50 50 CONTINUE NVEDUM=KIPART(IPART) IF (NPRT(9)) WRITE(MDEBUG,1004)NTOT,IPART,KPART,KKPART,NVEDUM 1004 FORMAT(' *CGHEI* SEC. GEN. NTOT,IPART,KPART,KKPART,KIPART = ', $ 5(I3,1X)) C --- INITIAL PARTICLE TYPE HAS BEEN CHANGED ==> PUT NEW TYPE ON --- C --- THE TEMPORARY STACK --- C --- MAKE CHOICE BETWEEN K0 LONG / K0 SHORT --- IF ((KPART .NE. 11) .AND. (KPART .NE. 12)) GO TO 52 CALL GRNDM(RNDM,1) KPART=11.5+RNDM(1) 52 CONTINUE C --- IN CASE THE NEW PARTICLE IS A NEUTRINO ==> FORGET IT --- IF (KPART .EQ. 2) GO TO 60 C --- PUT CURRENT GHEISHA PARTICLE ON THE CORSIKA STACK C --- ( IF SURVIVING ANGLE CUT ! ) NGKINE = 1 SECPAR(3) = -PZ C --- CALCULATE ELASTICITY IF ( EN .GT. EMAX ) THEN EMAX = EN ENDIF IF ( SECPAR(3) .GT. C(29) ) THEN ITY=IKPART(KPART) IF ( ITY .LT. 45 ) THEN SECPAR(1) = DBLE(ITY) ELSEIF ( ITY .EQ. 45 ) THEN SECPAR(1) = 201.D0 ELSEIF ( ITY .EQ. 46 ) THEN SECPAR(1) = 301.D0 ELSEIF ( ITY .EQ. 47 ) THEN SECPAR(1) = 402.D0 ENDIF IF ( ABS(AMAS) .LT. 1.E-9 ) THEN SECPAR(2) = EN ELSE SECPAR(2) = DBLE(EN) / DBLE(ABS(AMAS)) ENDIF IF ( PX .NE. 0. .OR. PY .NE. 0. ) THEN SECPAR(4) = ATAN2 ( DBLE(PY) , DBLE(PX) ) ELSE SECPAR(4) = 0.D0 ENDIF CALL TSTACK ENDIF C *** CHECK WHETHER SECONDARIES HAVE BEEN GENERATED AND COPY THEM *** C *** ALSO ON THE GEANT STACK *** 60 CONTINUE C --- ALL QUANTITIES ARE TAKEN FROM THE GHEISHA STACK WHERE THE --- C --- CONVENTION IS THE FOLLOWING --- C C EVE(INDEX+ 1)= X C EVE(INDEX+ 2)= Y C EVE(INDEX+ 3)= Z C EVE(INDEX+ 4)= NCAL C EVE(INDEX+ 5)= NCELL C EVE(INDEX+ 6)= MASS C EVE(INDEX+ 7)= CHARGE C EVE(INDEX+ 8)= TOF C EVE(INDEX+ 9)= PX C EVE(INDEX+10)= PY C EVE(INDEX+11)= PZ C EVE(INDEX+12)= TYPE IF (NTOT .LE. 0) GO TO 9999 C --- ONE OR MORE SECONDARIES HAVE BEEN GENERATED --- DO 61 L=1,NTOT INDEX=(L-1)*12 JND=EVE(INDEX+12) C --- MAKE CHOICE BETWEEN K0 LONG / K0 SHORT --- IF ((JND .NE. 11) .AND. (JND .NE. 12)) GO TO 63 CALL GRNDM(RNDM,1) JND=11.5+RNDM(1) C --- FORGET ABOUT NEUTRINOS --- 63 CONTINUE IF (JND .EQ. 2) GO TO 61 C --- SWITCH TO CORSIKA QUANTITIES --- ITY=IKPART(JND) IF (NPRT(9)) WRITE(MDEBUG,1006) ITY,NGKINE,L,(EVE(INDEX+J),J=1,12) 1006 FORMAT(' *CGHEI* GEANT PART. ',I3,' ALSO PUT ONTO STACK AT', $ ' POS. ',I3/ $ ' EVE(',I2,') = ',(' ',10G12.5)) PLX = EVE(INDEX+9) PLY = EVE(INDEX+10) PLZ = EVE(INDEX+11) PLSQ = PLX**2 + PLY**2 + PLZ**2 PLTOT = SQRT (PLSQ) RMASSK = ABS(RMASS(JND)) C FIND HIGHEST ENERGY PARTICLE FOR ELASTICITY EEESQ = PLSQ + RMASSK**2 IF ( EEESQ .GT. EMAX**2 ) THEN EMAX = SQRT(EEESQ) ENDIF C --- APPLY ANGLE CUT AND C --- ADD PARTICLE TO THE CORSIKA STACK (RESTRICTED TO 100) --- IF ( PLTOT .LE. 1.D-10 ) GOTO 61 SECPAR(3) = -PLZ / PLTOT IF ( SECPAR(3) .LE. C(29) ) GOTO 61 IF (NGKINE .GE. MXGKGH) GO TO 9999 NGKINE=NGKINE+1 IF ( ITY .LT. 45 ) THEN SECPAR(1) = DBLE(ITY) ELSEIF ( ITY .EQ. 45 ) THEN SECPAR(1) = 201.D0 ELSEIF ( ITY .EQ. 46 ) THEN SECPAR(1) = 301.D0 ELSEIF ( ITY .EQ. 47 ) THEN SECPAR(1) = 402.D0 ELSE SECPAR(1) = 0.D0 WRITE(MONIOU,*) '*CGHEI* ILLEGAL PARTICLE TYPE' ENDIF IF ( RMASSK .LT. 1.D-9 ) THEN SECPAR(2) = PLTOT ELSE SECPAR(2) = SQRT (PLSQ+RMASSK**2) / RMASSK ENDIF IF ( PLX .NE. 0.D0 .OR. PLY .NE. 0.D0 ) THEN SECPAR(4) = ATAN2 ( PLY,PLX ) ELSE SECPAR(4) = 0.D0 ENDIF CALL TSTACK 61 CONTINUE C --- COUNTER FOR ENERGY-MULTIPLICITY MATRIX MSMM = MSMM + NTOT C --- FILL ELASTICITY IN MATRICES ELASTI = EMAX/ENOLD MELL = MIN ( 1.D0+10.D0* MAX( 0.D0, ELASTI ) , 11.D0 ) MEN = MIN ( 4.D0+ 3.D0*LOG10(MAX( .1D0, EKINL )), 37.D0 ) IELDPM(MEN,MELL) = IELDPM(MEN,MELL) + 1 IELDPA(MEN,MELL) = IELDPA(MEN,MELL) + 1 IF ( ELASTI .LT. 1. ) THEN ELMEAN(MEN) = ELMEAN(MEN) + ELASTI ELMEAA(MEN) = ELMEAA(MEN) + ELASTI ENDIF 9999 CONTINUE C --- LIMIT THE VALUE OF NGKINE IN CASE OF OVERFLOW --- NGKINE=MIN(NGKINE,MXGKGH) RETURN END