source: trunk/MagicSoft/Simulation/Corsika/Mmcs/decay1.f@ 10099

      SUBROUTINE DECAY1( M0,M3,M4 )
 
C-----------------------------------------------------------------------
C  DECAY (INTO TWO PARTICLES)
C
C  TWO PARTICLE DECAY WITH FULL KINEMATIC; ENERGY AND MOMENTA CONSERVED
C  THIS SUBROUTINE IS CALLED FROM KDECAY, RESDEC, AND STRDEC
C  ARGUMENTS:
C   M0     = TYPE OF DECAYING PARTICLE
C   M3     = TYPE OF FIRST PRODUCT PARTICLE
C   M4     = TYPE OF SECOND PRODUCT PARTICLE
C
C  DESIGN  : D. HECK    IK3  FZK KARLSRUHE
C-----------------------------------------------------------------------
 
      IMPLICIT NONE
*KEEP,CONST.
      COMMON /CONST/   PI,PI2,OB3,TB3,ENEPER
      DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
*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
*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,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
*KEND.
 
      DOUBLE PRECISION AUX1,AUX2,AUX2A,AUX3,AUX4,COSTCM,COSTH3,COSTH4,
     *                 GAMMA3,GAMMA4,PHI4,WORK1,WORK2
      INTEGER          I,M0,M3,M4
C-----------------------------------------------------------------------
 
      IF ( DEBUG ) WRITE(MDEBUG,444) BETA,M0,M3,M4
  444 FORMAT(' DECAY1: BETA,M0,M3,M4=',1P,E10.3,3I5)
 
C  PARTICLE COORDINATES 5..10 ARE COPIED INTO SECPAR IN CALLING PROGRAM
C  CALCULATE AUXILIARY QUANTITIES
      AUX1   = ( ( PAMA(M0)**2 + PAMA(M3)**2 - PAMA(M4)**2 )
     *            / (2.D0*PAMA(M0)) )**2  - PAMA(M3)**2
      AUX2   = 1.D0 + AUX1 / PAMA(M3)**2
      AUX2A  = SQRT(AUX2)
      AUX3   = SQRT( 1.D0 - 1.D0 / AUX2 )
 
      WORK1  = GAMMA * AUX2A
      WORK2  = AUX3 * BETA * WORK1
 
C  DETERMINE POLAR ANGLE IN CM SYSTEM
      CALL RMMAR( RD,2,1 )
      COSTCM = 2.D0 * RD(1) - 1.D0
      GAMMA3 = WORK1 + WORK2 * COSTCM
C  SECOND PRODUCT PARTICLE WITH NONVANISHING REST MASS
      IF ( PAMA(M4) .NE. 0.D0 ) THEN
        GAMMA4 = (PAMA(M0)*GAMMA - PAMA(M3)*GAMMA3) / PAMA(M4)
        AUX4   = (PAMA(M0)**2 + PAMA(M4)**2 - PAMA(M3)**2 )
     *            / (2.D0*PAMA(M0)*PAMA(M4))
        COSTH4 = MIN( 1.D0, (GAMMA*GAMMA4 - AUX4)
     *                     / (BETA * GAMMA * SQRT(GAMMA4**2 - 1.D0)) )
      ELSE
C  SECOND PRODUCT PARTICLE IS GAMMA; THEN GAMMA4 IS THE ENERGY
        GAMMA4 = PAMA(M0)*GAMMA - PAMA(M3)*GAMMA3
        COSTH4 = MIN( 1.D0, (BETA - COSTCM)/(1.D0 - BETA*COSTCM) )
      ENDIF
      PHI4 = RD(2)*PI2
      CALL ADDANG( COSTHE,PHI, COSTH4,PHI4, SECPAR(3),SECPAR(4) )
      IF ( SECPAR(3) .GT. C(29) ) THEN
        SECPAR(1) = M4
        SECPAR(2) = GAMMA4
        IF ( DEBUG ) WRITE(MDEBUG,445) (SECPAR(I),I=1,9)
  445   FORMAT(' DECAY1: SECPAR=',1P,9E10.3)
        CALL TSTACK
      ENDIF
C  FIRST PRODUCT PARTICLE
      COSTH3 = MIN( 1.D0, (GAMMA * GAMMA3 - AUX2A)
     *                   / (BETA * GAMMA * SQRT(GAMMA3**2 - 1.D0)) )
      CALL ADDANG( COSTHE,PHI, COSTH3,PHI4+PI, SECPAR(3),SECPAR(4) )
      IF ( SECPAR(3) .GT. C(29) ) THEN
        SECPAR(1) = M3
        SECPAR(2) = GAMMA3
        IF ( DEBUG ) WRITE(MDEBUG,445) (SECPAR(I),I=1,9)
        CALL TSTACK
      ENDIF
 
      RETURN
      END