SUBROUTINE PIGEN1 C C********************************************************************* C DESIGN : D. HECK IK3 FZK KARLSRUHE C DATE : JUL 31, 1989 C********************************************************************* C THIS SUBROUTINE DESCRIBES THE PHOTONUCLEAR REACTION C GAMMA + NUCLEON -----> PION + NUCLEON C********************************************************************* DOUBLE PRECISION BETA,DUMMY,ENUCL,ESQ,E3CM,GAMMA DOUBLE PRECISION PEIG,PEOP,PT,PTRANS,P3CM,W0,W0I,W0S,W0SI *KEEP,ELABCT. COMMON /ELABCT/ ELCUT DOUBLE PRECISION ELCUT(4) *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 *KEND. DOUBLE PRECISION PI0MSQ COMMON/PION/PI0MSQ,PITHR,PICMAS,PI0MAS,AMASK0,AMASKC,AMASPR,AMASNT * *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,STACKE. COMMON/STACKE/ E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP DOUBLE PRECISION E(60),TIME(60) REAL X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60) INTEGER IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP *KEND. COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2 COMMON/ACLOCK/NCLOCK,JCLOCK C_____IF (NCLOCK.GT.JCLOCK) THEN C______WRITE(MDEBUG,* )' PIGEN1:NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP) C______CALL AUSGB2 C_____END IF IF(DEBUG)WRITE(MDEBUG,*)'PIGEN1: E=',E(NP) PEIG=E(NP) C*** NUMBERS AT THE VARIABLES MEAN : C*** 1 INCOMING GAMMA RAY C*** 2 HIT NUCLEON C*** 3 PRODUCED PION C*** 4 RECOILING NUCLEON C*** LOOK WHICH TYPE OF REACTION CALL RMMAR(RD,2,2) RNNO91=RD(1) RNNO92=RD(2) C*** 0.49923 IS THE FRACTION OF PROTONS IN AIR IF (RNNO91.LE.0.49923) THEN C *** HIT NUCLEON IS PROTON AMASS2=AMASPR C *** 33% CHANCE FOR CHARGE EXCHANGE IF (RNNO92.LE.0.333333) THEN C *** PI(+) + NEUTRON PRODUCED IQ(NP)=8 IQ(NP+1)=13 ELSE C *** PI(0) + PROTON PRODUCED IQ(NP)=7 IQ(NP+1)=14 END IF ELSE C *** HIT NUCLEON IS NEUTRON AMASS2=AMASNT C *** 33% CHANCE FOR CHARGE EXCHANGE IF (RNNO92.LE.0.333333) THEN C *** PI(-) + PROTON PRODUCED IQ(NP)=9 IQ(NP+1)=14 ELSE C *** PI(0) + NEUTRON PRODUCED IQ(NP)=7 IQ(NP+1)=13 END IF END IF AMAS2I=1./AMASS2 C*** NOTE: THE ENERGIES IN EGS ARE IN MEV, IN CORSIKA IN GEV AMASS3=PAMA(IQ(NP))*1.D3 AMASS4=PAMA(IQ(NP+1))*1.D3 C*** TOTAL LABORATORY ENERGY AND ITS INVERSE W0 =PEIG+AMASS2 W0I=1.D0/W0 C*** TOTAL.C.M. ENERGY AND INVERSE OF TOTAL C.M.ENERGY W0S = SQRT(AMASS2*(AMASS2+2.D0*PEIG)) W0SI=1.D0/W0S C*** THRESHOLD ENERGY ETH=0.5*((AMASS3+AMASS4)**2-AMASS2**2)*AMAS2I C*** BETA,GAMMA, ESQ, BRATIO, G3 ARE AUXILIARY QUANTITIES BETA=PEIG*W0I GAMMA=W0*W0SI ED =0.5*((AMASS3-AMASS4)**2-AMASS2**2)*AMAS2I ESQ = SQRT((PEIG-ETH)*(PEIG-ED)) BRATIO = PEIG/ESQ G3 = W0I*BRATIO*(PEIG-ETH+AMASS3*AMAS2I*(AMASS3+AMASS4)) C*** C.M. ENERGY OF PION E3CM=G3*AMASS2*GAMMA/BRATIO C*** C.M. PION MOMENTUM P3CM=AMASS2*W0SI*ESQ B3CM2=P3CM**2/(P3CM**2+AMASS3**2) B3CM=SQRT(B3CM2) C*** DETERMINE THETA IN C.M. SYSTEM BY CHANCE. IF (PEIG.LE.900.D0) THEN C *** PHOTON ENERGY IS BELOW 900 MEV 210 CONTINUE CALL RMMAR(RD,2,2) RNNO93=RD(1) RNNO94=RD(2) IF (IQ(NP).EQ.7) THEN C *** NEUTRAL PION EMITTED, TAKE PURE C *** DIPOLE RADIATION: W(COSTH) = 1-3/5*COSTH**2 COSTE3 = 2.*RNNO93-1. IF((RNNO94 .GT. 1.-0.6*COSTE3**2))GOTO 210 ELSE C *** CHARGED PION EMITTED, TAKE MODIFIED DIPOLE RADIATION C *** WITH ASYMMETRY TERM 1/(1-BETACM*COSTE3)**2 COSTE3 = 1./B3CM-1./(RNNO93*2.*B3CM2/(1.-B3CM2)+B3CM/(1.+B3CM)) IF((RNNO94*2.5 .GT. 1.+COSTE3*(-1.8+COSTE3*(.65+COSTE3*(.34 -.18 * *COSTE3 )))))GOTO 210 END IF ELSE IF(PEIG.LE.1300.D0) THEN C *** PHOTON ENERGY BETWEEN 900 AND 1300 MEV 220 CONTINUE CALL RMMAR(RD,2,2) RNNO93=RD(1) RNNO94=RD(2) IF (IQ(NP).EQ.7) THEN C *** NEUTRAL PION EMITTED, TAKE PURE QUADRUPOLE C *** RADIATION: W(COSTH) = 1+6*COSTH**2-5*COSTH**4 COSTE3 = 2.*RNNO93-1. IF((2.8*RNNO94 .GT. 1.+6.*COSTE3**2-5.*COSTE3**4))GOTO 220 ELSE C *** CHARGED PION EMITTED, TAKE MODIFIED QUADRUPOLE C *** RADIATION WITH ASYMMETRY TERM: 1/(1-BETACM*COSTE3)**2 COSTE3 = 1./B3CM-1./(RNNO93*2.*B3CM2/(1.-B3CM2)+B3CM/(1.+B3CM)) IF((13.2*RNNO94 .GT. 1.+COSTE3*(-2.18+COSTE3*(7.20+COSTE3*(-2.55 * +COSTE3*(-15.39+COSTE3*(6.36+COSTE3*(13.80-COSTE3*8.235)))))))) * GOTO 220 END IF ELSE C *** ABOVE 1300 MEV THE ANGULAR DISTRIBUTION IS DETERMINED C *** BY THE TRANSVERSE MOMENTUM OF THE PION PT=1.D3*PTRANS(DUMMY) COSTE3=SQRT(MAX(0.D0,P3CM**2-PT**2))/P3CM END IF C*** PRECISE ENERGY OUTGOING PION = PEOP PEOP =GAMMA*(E3CM+BETA*P3CM*COSTE3) C*** ENERGY OF OUTGOING PION IN STACK POSITION NP E(NP)=PEOP C*** MOMENTUM OF OUTGOING PION = AMOM3 C*** COSTHE AND SINTHE ARE ANGLES IN LAB SYSTEM FOR PARTICLE 3 (PION) C*** SEE SLAC-265, P. 52 AMOM3=SQRT(MAX(0.D0,PEOP**2-AMASS3**2)) IF (AMOM3.GT.0.) THEN COSTHE=(AMASS4**2-AMASS2**2-AMASS3**2+2.*PEOP*W0-2.*PEIG*AMASS2) * /(2.*PEIG* AMOM3) ELSE COSTHE=1. END IF SINTHE= SQRT(MAX(0.0,1.-COSTHE**2)) CALL UPHI(2,1) C*** TOTAL ENERGY OF RECOILING NUCLEON = ENUCL ENUCL=W0-PEOP IF ((ENUCL-AMASS4).GT.ELCUT(1)*1000.D0) THEN C *** RECOIL ENERGY IS TOO LARGE, MUST TREAT THE NUCLEON NP=NP+1 E(NP)=ENUCL C *** MOMENTUM OF RECOIL NUCLEON AMOM4=SQRT(ENUCL**2-AMASS4**2) C *** COSTHE AND SINTHE ARE ANGLES IN LAB SYSTEM FOR RECOIL NUCLEON C *** SEE SLAC-265, P. 52 COSTHE=(AMASS3**2-AMASS2**2-AMASS4**2+2.*ENUCL*W0-2.*PEIG*AMASS2) * / (2. * PEIG*AMOM4) SINTHE=-SQRT(MAX(0.0,1.-COSTHE**2)) CALL UPHI(3,2) IF (W(NP).GT.C(29)) THEN C *** ANGLE WITH RESPECT TO X AXIS IF (U(NP)**2+V(NP)**2.GT.3.E-38) THEN ANGLEX = -ATAN2(V(NP),U(NP)) ELSE ANGLEX = 0. END IF C *** ADD NUCLEON TO CORSIKA STACK SECPAR(1)=IQ(NP) SECPAR(2)=E(NP)/AMASS4 SECPAR(3)=W(NP) SECPAR(4)=ANGLEX SECPAR(5)=-Z(NP) SECPAR(6)=TIME(NP) SECPAR(7)=X(NP) SECPAR(8)=-Y(NP) SECPAR(11)=1.D0 SECPAR(12)=0.D0 CALL TSTOUT END IF C *** ELIMINATE NUCLEON FROM EGS-STACK NP=NP-1 END IF C*** END OF RECOIL NUCLEON TREATEMENT CASE RETURN END