      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
