SUBROUTINE NKG( ENERN ) C----------------------------------------------------------------------- C N(ISHIMURA) K(AMATA) G(REISEN) C C CALCULATES ELECTROMAGNETIC COMPONENT OF SHOWERS USING THE ANALYTIC C NKG FORMULAS, INCLUDING ELECTRON ENERGY THRESHOLD ELCUT(3) C SEE J.N. CAPDEVIELLE, 22ND ICRC, DUBLIN 1991, CONTRIB. HE 3.5.10 C THIS SUBROUTINE IS CALLED FROM EM C ARGUMENTS: C ENERN = ENERGY OF ELECTRON/PHOTON GENERATING A SUBSHOWER C NEGATIVE FOR SUBSHOWERS TO BE SUBTRACTED AFTER C PHOTONUCLEAR REACTION C----------------------------------------------------------------------- IMPLICIT NONE *KEEP,CONST. COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER *KEEP,ELABCT. COMMON /ELABCT/ ELCUT DOUBLE PRECISION ELCUT(4) *KEEP,NKGI. COMMON /NKGI/ SEL,SELLG,STH,ZEL,ZELLG,ZSL,DIST, * DISX,DISY,DISXY,DISYX,DLAX,DLAY,DLAXY,DLAYX, * OBSATI,RADNKG,RMOL,TLEV,TLEVCM,IALT DOUBLE PRECISION SEL(10),SELLG(10),STH(10),ZEL(10),ZELLG(10), * ZSL(10),DIST(10), * DISX(-10:10),DISY(-10:10), * DISXY(-10:10,2),DISYX(-10:10,2), * DLAX (-10:10,2),DLAY (-10:10,2), * DLAXY(-10:10,2),DLAYX(-10:10,2), * OBSATI(2),RADNKG,RMOL(2),TLEV(10),TLEVCM(10) INTEGER IALT(2) *KEEP,NKGS. COMMON /NKGS/ CZX,CZY,CZXY,CZYX,SAH,SL,ZNE DOUBLE PRECISION CZX(-10:10,2),CZY(-10:10,2),CZXY(-10:10,2), * CZYX(-10:10,2),SAH(10),SL(10),ZNE(10) *KEEP,OBSPAR. COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP, * THETPR,PHIPR,NOBSLV DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10), * THETAP,THETPR(2),PHIP,PHIPR(2) INTEGER NOBSLV *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,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 AE,AS,ASE,AUXIL,BS,CCP,CPC,CPCP,CPH,CSGA, * DE,DISTL,ECRI,ECR1,ECR2,ENERN,GAM,GRCUT, * G1,G2,G3,S,SC1,SC2,SIGNE,SM,SMRM, * SQRZ1I,SQZC1I,SQZC2I,SS2,SS45,TEX,THICK,THICKP, * XMOL,XNE,XS,X0,YM,YS,ZC1,ZC2,ZG1,ZG2,ZG3,Z1 INTEGER IL,IOL,M EXTERNAL GAM,THICK C X0 IS RADIATON LENGTH IN AIR (G/CM**2) C (SEE ALSO MIKOCKI ET AL. J.PHYS.G.:NUCL.PART.PHYS. 17 (1991) 1303 ) C GRCUT IS GREISEN CUT OFF, ECRI IS CRITICAL ENERGY IN AIR C ECR2 IS 0.4 * ECRI DATA X0 / 37.1D0 /, GRCUT / 0.1D0 /, ECRI / 0.082D0 / DATA ECR2 / 0.0328D0 / C----------------------------------------------------------------------- IF (DEBUG) WRITE(MDEBUG,*)'NKG : ',SNGL(SECPAR(1)),SNGL(ENERN) C CHECK WETHER SUBSHOWER IS SUBTRACTED IF ( ENERN .GE. 0.D0 ) THEN SIGNE = +1.D0 ELSE ENERN = -ENERN SIGNE = -1.D0 ENDIF C ENERGY CUT OFF IN GREISEN FORMULA C (EM PARTICLE BELOW THIS CUT CAN NOT PRODUCE A SHOWER) IF ( ENERN .LT. GRCUT ) RETURN C DON'T CALCULATE NKG FOR BACKWARD GOING PARTICLES IF ( SECPAR(3) .LE. 0.D0 ) RETURN C DON'T CALCULATE NKG IF PARTICLE BELOW THE LOWEST OBSERVATION LEVEL IF ( SECPAR(5) .LT. OBSATI(1) ) RETURN Z1 = LOG(ENERN / ECRI) SQRZ1I = 1.D0 / SQRT(Z1) C THIS CUT IS ONLY IMPORTANT FOR ELCUT > .0672 ECR1 = ECR2 + ELCUT(3) IF ( ENERN .LT. ECR1 ) RETURN ZC1 = LOG(ENERN / ECR1) SQZC1I = 1.D0 / SQRT(ZC1) C LOG(ENERN/ECR2) IS LOG(ENERN / ECRI) - LOG(0.4) ZC2 = Z1 + 0.916290732D0 SQZC2I = 1.D0 / SQRT(ZC2) THICKP = THICK(SECPAR(5)) C LOOP OVER LEVELS DO 14 IL = 1,IALT(1) C DISREGARD LEVELS ABOVE THE PARTICLE IF ( TLEVCM(IL) .GT. SECPAR(5) ) GOTO 14 C DISTANCE IN G/CM**2 .... (ALONG PHOTON-AXIS) IN RADIATION LENGTHS XMOL = (TLEV(IL) - THICKP) / ( X0 * SECPAR(3) ) C CORRECT DEPTH FOR SUBSHOWERS TO BE SUBTRACTED BY 9/7 IF ( SIGNE .LT. 0.D0 ) XMOL = XMOL + 1.285714286D0 C XMOL IS DEPTH IN RADIATION LENGTHS IF ( XMOL .GT. 60.D0 .OR. XMOL .LT. 1.D0 ) GOTO 14 C S IS AGE PARAMETER S = 3.D0 * XMOL / (XMOL + 2.D0 * Z1) IF ( S .LE. 0.2D0 ) GOTO 14 SC1 = 3.D0 * XMOL / (XMOL + 2.D0 * ZC1) SC2 = 3.D0 * XMOL / (XMOL + 2.D0 * ZC2) C ELECTRON NUMBER AT OBSERVATION LEVEL CPH = .31D0 * EXP( XMOL * (1.D0 - 1.5D0 * LOG(S) ) ) * SQRZ1I CPC = EXP( XMOL * ( 1.D0 - 1.5D0 * LOG(SC1) ) ) * SQZC1I CCP = EXP( XMOL * ( 1.D0 - 1.5D0 * LOG(SC2) ) ) * SQZC2I CPCP = SIGNE * CPH * CPC / CCP C INTERMEDIATE FACTORS FOR LATERAL DISTRIBUTION AND AGE PARAMETER AE = 4.D0 * EXP( 0.915D0 * (S - 1.D0) ) / S DE = ( 1.D0 + S ) / ( 1.15D0 + 0.15D0 * S ) ASE = AE**DE ZG3 = GAM( (S + 2.D0) * DE ) IF ( ZG3 .LE. 0.D0 ) GOTO 14 ZG1 = GAM(S * DE) ZG2 = GAM( (S + 1.D0) * DE ) AUXIL = 4.D0 / (S * ASE) XNE = CPCP * ( ZG2 + AUXIL * ZG3 ) / ( ASE * (ZG1 + AUXIL*ZG2) ) C SUM OF N_E AT FIXED LEVEL ZNE(IL) = ZNE(IL) + XNE SL(IL) = SL(IL) + CPCP C CALCULATE THE ELECTRON LATERAL DISTRIBUTION FOR THE 2 SELECTED C OBSERVATION LEVELS IF ( IL .EQ. IALT(1) ) THEN IOL = 1 ELSEIF ( IL .EQ. IALT(2) ) THEN IOL = 2 ELSE GOTO 14 ENDIF C CALCULATION OF LATERAL ELECTRON DISTRIBUTION IF ( SC1 .GE. 2.25D0 ) GOTO 14 G1 = GAM(4.5D0 - SC1) G2 = GAM(SC1) G3 = GAM(4.5D0 - 2.D0 * SC1) C DISTANCE IN CM BETWEEN PHOTON INITIATION AND OBSERVATION (VERTICAL) DISTL = SECPAR(5) - TLEVCM(IL) C MODULATION BY AGE PARAMETER FOLLOWING LAGUTIN & UCHAIKIN C (AGE PARAMETER LIES BETWEEN 0.2 AND 2.25) SM = 0.78D0 - 0.21D0 * SC1 SMRM = 1.D0 / ( SM * RMOL(IOL) ) CSGA = CPCP * SMRM**2 * G1 / ( PI2 * G2 * G3 ) SS2 = SC1 - 2.D0 SS45 = SC1 - 4.5D0 AS = SIN( SECPAR(4) ) BS = COS( SECPAR(4) ) TEX = DISTL * SQRT( 1.D0 - SECPAR(3)**2 ) / SECPAR(3) C DISTANCE TO THE CENTER OF THE CASCADE (IN CM) XS = SECPAR(7) + TEX * BS - XOFF(NOBSLV+1-IOL) YS = SECPAR(8) + TEX * AS - YOFF(NOBSLV+1-IOL) C NKG-FORMULA C LOOP OVER ALL LATERAL DISTANCES GETTING THE DENSITY IN MOLIERE UNITS DO 171 M = -10,10 IF ( M .EQ. 0 ) GOTO 171 C X DIRECTION YM = SMRM * MAX( SQRT((DISX(M)-XS)**2 + YS**2), 1.D0 ) CZX (M,IOL) = CZX (M,IOL) + CSGA * YM**SS2 * (YM+1.D0)**SS45 C Y DIRECTION YM = SMRM * MAX( SQRT(XS**2 + (DISY(M)-YS)**2), 1.D0 ) CZY (M,IOL) = CZY (M,IOL) + CSGA * YM**SS2 * (YM+1.D0)**SS45 C XY DIRECTION YM = SMRM * * MAX( SQRT((DISXY(M,1)-XS)**2 + (DISXY(M,2)-YS)**2), 1.D0 ) CZXY(M,IOL) = CZXY(M,IOL) + CSGA * YM**SS2 * (YM+1.D0)**SS45 C YX DIRECTION YM = SMRM * * MAX( SQRT((DISYX(M,1)-XS)**2 + (DISYX(M,2)-YS)**2), 1.D0 ) CZYX(M,IOL) = CZYX(M,IOL) + CSGA * YM**SS2 * (YM+1.D0)**SS45 171 CONTINUE 14 CONTINUE RETURN END