source: trunk/MagicSoft/Simulation/Corsika/Mmcs/pigen1.f@ 4787

      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