SUBROUTINE PIGEN2 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 + PION + NUCLEON C********************************************************************* DOUBLE PRECISION BETA,DUMMY,ECM,ENUCL,GAMMA,PEIG,PTRANS *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/UPHIIN/SINC0,SINC1,SIN0(20002),SIN1(20002) COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2 COMMON/ACLOCK/NCLOCK,JCLOCK C_____IF (NCLOCK.GT.JCLOCK) THEN C______WRITE(MDEBUG,* )' PIGEN2:NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP) C______CALL AUSGB2 C_____END IF IF(DEBUG)WRITE(MDEBUG,*)'PIGEN2: E=',E(NP) PEIG=E(NP) C*** NUMBERS AT THE VARIABLES MEAN : C*** 1 INCOMING GAMMA RAY C*** 2 HIT NUCLEON C*** 3 FIRST PRODUCED PION C*** 4 SECOND PRODUCED PION C*** 5 RECOILING NUCLEON CALL RMMAR(RD,2,2) RNNO81=RD(1) RNNO82=RD(2) C*** LOOK WHICH TYPE OF REACTION C*** 0.49923 IS THE FRACTION OF PROTONS IN AIR IF (RNNO81.LE.0.49923) THEN C *** HIT NUCLEON IS PROTON AMASS2=AMASPR C *** BRANCHING FOR COLLISION WITH PROTON IF (RNNO82.LE.0.3) THEN C *** PI(0) + PI(0) + PROTON IQ(NP)= 7 IQ(NP+1)= 7 IQ(NP+2)= 14 ELSE IF(RNNO82.LE.0.6) THEN C *** PI(+) + PI(-) + PROTON IQ(NP)= 8 IQ(NP+1)= 9 IQ(NP+2)= 14 ELSE C *** PI(+) + PI(0) + NEUTRON IQ(NP)= 8 IQ(NP+1)= 7 IQ(NP+2)= 13 END IF ELSE C *** HIT NUCLEON IS NEUTRON C *** BRANCHING FOR COLLISION WITH NEUTRON AMASS2=AMASNT IF (RNNO82.LE.0.3) THEN C *** PI(0) + PI(0) + NEUTRON IQ(NP)= 7 IQ(NP+1)= 7 IQ(NP+2)= 13 ELSE IF(RNNO82.LE.0.6) THEN C *** PI(+) + PI(-) + NEUTRON IQ(NP)= 8 IQ(NP+1)= 9 IQ(NP+2)= 13 ELSE C *** PI(-) + PI(0) + PROTON IQ(NP)= 9 IQ(NP+1)= 7 IQ(NP+2)= 14 END IF END IF C*** CALCULATE AUXILIARY PARAMETERS ECM=SQRT(AMASS2*(AMASS2+2.D0*PEIG)) 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 AMASS5=PAMA(IQ(NP+2))*1.D3 AUX1=(AMASS3+AMASS4)**2 AUX2A=(ECM - AMASS5)**2 AUX2=AUX2A-AUX1 AUX3=(AMASS3+AMASS5)**2 AUX4A=(ECM - AMASS4)**2 AUX4=AUX4A-AUX3 AUX5=AMASS3**2-AMASS4**2 AUX6=ECM**2-AMASS5**2 AUX7=0.5/ECM AUX8=(ECM - AMASS3)**2 BETA=PEIG/(AMASS2+PEIG) GAMMA=2.*(PEIG+AMASS2)*AUX7 230 CONTINUE CALL RMMAR(RD,2,2) RNNO84=RD(1) RNNO85=RD(2) C*** ARE INVARIANT MASS SQUARES INSIDE BOUNDARY OF DALITZ PLOT? AM34SQ=AUX2*RNNO84+AUX1 AM35SQ=AUX4*RNNO85+AUX3 AM34I=0.5/SQRT(AM34SQ) E3STAR=(AUX5+AM34SQ)*AM34I E5STAR=(AUX6-AM34SQ)*AM34I ROOT1=SQRT(E3STAR**2-AMASS3**2) ROOT2=SQRT(E5STAR**2-AMASS5**2) C*** REJECT RANDOM NUMBERS, IF NOT INSIDE KINEMATIC BOUNDARY DISCR=AM35SQ-(E3STAR+E5STAR)**2 IF((DISCR.GT.-(ROOT1-ROOT2)**2))GOTO 230 IF((DISCR.LT.-(ROOT1+ROOT2)**2))GOTO 230 C*** E3CM,E4CM,E5CM ARE ENERGIES IN C.M. SYSTEM E4CM=(ECM**2+AMASS4**2-AM35SQ)*AUX7 E5CM=(ECM**2+AMASS5**2-AM34SQ)*AUX7 C*** NOW TAKE PION WITH HIGHEST ENERGY AS PARTICLE 3 E3CM=ECM-E4CM-E5CM IF (E4CM.GT.E3CM) THEN C *** INTERCHANGE PARTICLE 3 AND 4 HELP=E3CM E3CM=E4CM E4CM=HELP HELP=AMASS3 AMASS3=AMASS4 AMASS4=HELP IHELP=IQ(NP) IQ(NP)=IQ(NP+1) IQ(NP+1)=IHELP END IF C*** P3CM,P4CM,P5CM ARE MOMENTA IN C.M. SYSTEM C*** P3SQ,P4SQ,P5SQ ARE SQUARED MOMENTA IN C.M. SYSTEM P3SQ=E3CM**2-AMASS3**2 P3CM=SQRT(MAX(0.,P3SQ)) P4SQ=E4CM**2-AMASS4**2 P4CM=SQRT(MAX(0.,P4SQ)) P5SQ=E5CM**2-AMASS5**2 P5CM=SQRT(MAX(0.,P5SQ)) COSA=(P5SQ-P3SQ-P4SQ)/(2.*P3CM*P4CM) SINA=-SQRT(MAX(0.,1.-COSA**2)) COSB=(P4SQ-P3SQ-P5SQ)/(2.*P3CM*P5CM) SINB= SQRT(MAX(0.,1.-COSB**2)) C*** NOW SELECT THE THREE INDEPENDENT ANGLES IN C.M. SYSTEM PT3=1.D3*PTRANS(DUMMY) SIN3CM=MIN(1.,PT3/P3CM) COS3CM=SQRT(1.-SIN3CM**2) CALL RMMAR(RNNO86,1,2) PSI=TWOPI*RNNO86 LPSI=SINC1*PSI+SINC0 SINPSI=SIN1(LPSI)*PSI+SIN0(LPSI) CPSI=PI5D2-PSI LCPSI=SINC1*CPSI+SINC0 COSPSI=SIN1(LCPSI)*CPSI+SIN0(LCPSI) C*** THIRD INDEPENDENT ANGLE PHI IS CHOOSEN LATER IN SUBROUTINE UPHI C*** NOW MAKE LORENTZ-TRANSFORMATION FOR PARTICLE 3 (PION) E(NP)=GAMMA*(E3CM+BETA*P3CM*COS3CM) C*** COSTHE AND SINTHE ARE ANGLES IN LAB SYSTEM FOR PARTICLE 3 (PION) COSTHE= MIN((BETA*E3CM+P3CM*COS3CM)*GAMMA/ SQRT(MAX(0.D0,E(NP)**2 *-AMASS3**2)),1.D0) SINTHE=SQRT(MAX(0.0,1.-COSTHE**2)) C*** SINPHI AND COSPHI ARE NOW SET IN SUBROUTINE UPHI CALL UPHI(2,1) SINFI3=SINPHI COSFI3=COSPHI C*** NOW MAKE LORENTZ-TRANSFORMATION FOR PARTICLE 4 = PION COS4CM=COS3CM*COSA-SIN3CM*COSPSI*SINA NP=NP+1 E(NP)=GAMMA*(E4CM+BETA*P4CM*COS4CM) SINT4=SQRT(MAX(0.,1.-COS4CM**2)) IF (SINT4.NE.0.) THEN SINT4I =1./SINT4 AUXA=COS3CM*COSPSI*SINA+SIN3CM*COSA C *** COSPHI AND SINPHI ARE IN LAB SYSTEM FOR PARTICLE 4 (PION) COSPHI=(COSFI3*AUXA-SINFI3*SINPSI*SINA)*SINT4I SINPHI=(SINFI3*AUXA+COSFI3*SINPSI*SINA)*SINT4I ELSE COSPHI=0. SINPHI=1. END IF C*** COSTHE AND SINTHE ARE IN LAB SYSTEM FOR PARTICLE 4 (PION) COSTHE= MIN((BETA*E4CM+P4CM*COS4CM)*GAMMA/ SQRT(MAX(0.D0,E(NP)**2 *-AMASS4**2)),1.D0) SINTHE=SQRT(MAX(0.0,1.-COSTHE**2)) CALL UPHI(3,2) C*** NOW MAKE LORENTZ-TRANSFORMATION FOR PARTICLE 5 = RECOIL NUCLEON COS5CM=COS3CM*COSB-SIN3CM*COSPSI*SINB ENUCL=GAMMA*(E5CM+BETA*P5CM*COS5CM) IF ((ENUCL-AMASS5).GT.ELCUT(1)*1000.D0) THEN C *** RECOIL NUCLEON IS ABOVE THRESHOLD AND MUST BE TREATED NP=NP+1 E(NP)=ENUCL SINT5=SQRT(MAX(0.,1.-COS5CM**2)) IF (SINT5.NE.0.) THEN SINT5I =1./SINT5 AUXB=COS3CM*COSPSI*SINB+SIN3CM*COSB C *** COSPHI AND SINPHI ARE IN LAB SYSTEM FOR PART. 5 (NUCLEON) COSPHI=(COSFI3*AUXB-SINFI3*SINPSI*SINB)*SINT5I SINPHI=(SINFI3*AUXB+COSFI3*SINPSI*SINB)*SINT5I ELSE COSPHI=0. SINPHI=1. END IF C *** COSTHE AND SINTHE ARE IN LAB SYSTEM FOR PARTICLE 5 (NUCLEON) COSTHE=MIN((BETA*E5CM+P5CM*COS5CM)*GAMMA/SQRT(ENUCL**2-AMASS5**2) * , 1.D0) 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)/AMASS5 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