      SUBROUTINE DECAY6(AM0,AM3,AM4,AM5,PARAMA,PARAMB,PARAMC,AMPMX,MODE)
 
C-----------------------------------------------------------------------
C  DECAY (INTO 3 PARTICLES)
C
C  TREATES DECAY INTO 3 PARTICLES; FULLY CONSERVING ENERGY AND MOMENTA
C  KINEMATIC RANGE PARAMETRISATION SEE PHYS. LETT. 204B (1988) 90-91
C  FOR LEPTONIC KAON DACAY: THE POLARIZATION OF THE MUON AND
C  THE NEUTRINO PRODUCTION IS INCLUDED.
C  THIS SUBROUTINE IS CALLED FROM ETADEC, KDECAY, AND PI0DEC
C  ARGUMENTS:
C   AM0    = MASS OF DECAYING PARTICLE
C   AM3, AM4, AM5 = MASSES OF RESULTING PARTICLES
C   PARAMA = DALITZ AMPLITUDE PARAMETER (SEE BELOW)
C   PARAMB = DALITZ AMPLITUDE PARAMETER (SEE BELOW)
C   PARAMC = DALITZ AMPLITUDE PARAMETER (SEE BELOW)
C   AMPMX  = MAXIMUM AMPLITUDE OF DALITZ PLOT
C   MODE   = 1  FOR DECAY KAON  ----> 3 PIONS
C          = 2  FOR DECAY ETA   ----> 3 PIONS OR 2 PIONS + GAMMA
C               FOR DECAY PI(0) ----> ELECTRON + POSITRON + GAMMA
C          = 3  FOR DECAY KAON  ----> PION + MUON + NEUTRINO
C          = 4  FOR DECAY KAON  ----> PION + ELECTRON + NEUTRINO
C
C  AMPLITUDE PARAMETERS PARAMA, PARAMB, PARAMC ARE DEPENDENT ON MODE:
C  FOR MODE=1: PARAMA = G      DALITZ AMPLITUDE PARAMETRISATION SEE
C              PARAMB = H      PHYS. LETT. 204B (1988) 181 - 193
C              PARAMC = K
C
C  FOR MODE=2: PARAMA = A      DALITZ AMPLITUDE PARAMETRISATION SEE
C              PARAMB = DUMMY  PHYS. LETT. 204B (1988) 173 - 175;
C              PARAMC = DUMMY  J.G. LAYTER ET.AL. PHYS.REV.D7(1973)2565
C
C  FOR MODE>2: PARAMA = LAMBDA-PLUS    DALITZ AMPLITUDE PARAMETRISATION
C              PARAMB = LAMBDA-ZERO    SEE PHYS. LETT. 204B (1988)
C              PARAMC = DUMMY          182 - 194
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,DECAY.
      COMMON /DECAY/   GAM345,COS345,PHI345
      DOUBLE PRECISION GAM345(3),COS345(3),PHI345(3)
*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,POLAR.
      COMMON /POLAR/   POLART,POLARF
      DOUBLE PRECISION POLART,POLARF
*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 ABYM,AMPLI,AMPMX,AM0,AM3,AM34I,AM34SQ,AM35SQ,
     *                 AM4,AM5,APARAL,APERPN,AUXA,AUXB,AUX1,AUX2,AUX2A,
     *                 AUX3,AUX4,AUX4A,AUX5,AUX6,
     *                 AUX7,AUX8,AUX10,AUX12,AUX14,BBYM,BOFQ,
     *                 CM0SQ,CM3SQ,CM3SQI,CM4SQ,CM5SQ,COSALF,COSBET,
     *                 COSFI4,COSOME,COSPHI,COSPSI,COS3CM,COS4CM,COS5CM,
     *                 DISCR,EPIPRM,E3CM,E3STAR,E4CM,E5CM,E5STAR,FACT,
     *                 GRLAMD,OMEGA,PA,PARAMA,PARAMB,PARAMC,PB,PC,PSI,
     *                 P3CM,P3SQ,P4CM,P4SQ,P5CM,P5SQ,ROOT1,ROOT2,
     *                 SINALF,SINBET,SINFI4,SINFI5,SINOMG,SINPHI,SINPSI,
     *                 SINT4,SINT4I,SINT5I,SIN3CM,S0,TBYMSS,XIT,XI0
      INTEGER          MODE
C-----------------------------------------------------------------------
 
      IF ( DEBUG ) WRITE(MDEBUG,444) AM0,AM3,AM4,AM5
  444 FORMAT(' DECAY6: AM0',1P,E10.3,' AM3',E10.3,' AM4',E10.3,
     *       ' AM5',E10.3)
 
C  CALCULATE AUXILIARY QUANTITIES
      CM0SQ = AM0**2
      CM3SQ = AM3**2
      CM4SQ = AM4**2
      CM5SQ = AM5**2
      AUX1  = (AM3 + AM4)**2
      AUX2A = (AM0 - AM5)**2
      AUX2  = AUX2A - AUX1
      AUX3  = (AM3 + AM5)**2
      AUX4A = (AM0 - AM4)**2
      AUX4  = AUX4A - AUX3
      AUX5  = CM3SQ - CM4SQ
      AUX6  = CM0SQ - CM5SQ
      AUX7  = 0.5D0 / AM0
      IF     ( MODE .EQ. 1 ) THEN
        AUX8   = (AM0 - AM3)**2
        S0     = OB3 * ( CM0SQ + CM3SQ + CM4SQ + CM5SQ )
        AUX10  = 1.D0 / PAMA(8)**2
      ELSEIF ( MODE .EQ. 2 ) THEN
        AUX14  = 1.D0 / (AM0 - AM3 - AM4 - AM5)
      ELSEIF ( MODE .EQ. 3  .OR.  MODE .EQ. 4 ) THEN
        CM3SQI = 1.D0 / CM3SQ
        AUX12  = (CM0SQ + CM3SQ - CM4SQ) * AUX7
C  XI0 IS XI(0); GRLAMD IS GREAT LAMBDA
        XI0    = ( CM0SQ - CM3SQ) * CM3SQI * (PARAMB - PARAMA)
        GRLAMD = -XI0 * PARAMA
      ELSE
        WRITE(MONIOU,*) 'DECAY6: UNEXPECTED MODE =',MODE
        RETURN
      ENDIF
 
 100  CALL RMMAR( RD,3,1 )
C  ARE INVARIANT MASS SQUARES INSIDE BOUNDARY OF DALITZ PLOT?
      AM34SQ = AUX2 * RD(1) + AUX1
      AM35SQ = AUX4 * RD(2) + AUX3
      AM34I  = 0.5D0 / SQRT(AM34SQ)
      E3STAR = (AUX5 + AM34SQ) * AM34I
      E5STAR = (AUX6 - AM34SQ) * AM34I
      ROOT1  = SQRT(E3STAR**2 - CM3SQ )
      ROOT2  = SQRT(E5STAR**2 - CM5SQ )
      DISCR  = AM35SQ - (E3STAR + E5STAR)**2
C  REJECT RANDOM NUMBERS, IF OUTSIDE KINEMATIC BOUNDARY OF DALITZ PLOT
      IF ( DISCR .GT. -(ROOT1 - ROOT2)**2 ) GOTO 100
      IF ( DISCR .LT. -(ROOT1 + ROOT2)**2 ) GOTO 100
C  E3CM, E4CM, E5CM ARE ENERGIES IN THE C. M. SYSTEM
      E4CM   = (CM0SQ + CM4SQ - AM35SQ) * AUX7
      E5CM   = (CM0SQ + CM5SQ - AM34SQ) * AUX7
      E3CM   = AM0 - E4CM - E5CM
 
      IF     ( MODE .EQ. 1 ) THEN
        FACT  = AUX10 * (AUX2A - 2.D0*AM0*(E5CM-AM5) - S0)
C  AMPLITUDE OF SQUARED MATRIX ELEMENT (SEE PHYS. LETT. B204 (1988) 181)
        AMPLI = 1.D0 + PARAMA*FACT + PARAMB*FACT**2 + PARAMC*( AUX10
     *           * ( AUX4A -AUX8 -2.D0*(E4CM-AM4-E3CM+AM3)*AM0 ) )**2
 
      ELSEIF ( MODE .EQ. 2 ) THEN
C  AMPLITUDE OF SQUARED MATRIX ELEMENT (SEE PHYS. LETT. B204 (1988) 173)
C  REF: J. G. LAYTER ET AL., PHYS. REV. D7 (1973) 2565
        AMPLI = 1.D0 + PARAMA * ( 3.D0 * (E5CM - AM5) * AUX14 - 1.D0 )
 
      ELSE
C  EPIPRM IS (ENERGY OF PION)PRIMED
        EPIPRM = AUX12 - E3CM
C  PA, PB, AND PC ARE THE A, B, AND C PARAMETERS
        PA     = AM0 * ( 2.D0 * E4CM * E5CM - AM0 * EPIPRM )
     *            + CM4SQ * ( 0.25D0 * EPIPRM - E5CM )
        PB     = CM4SQ * ( E5CM - 0.5D0 * EPIPRM )
        PC     = CM4SQ * EPIPRM * 0.25D0
C  TBYMSS IS T DIVIDED BY MASS SQUARE OF PION
        TBYMSS = (CM0SQ + CM3SQ - 2.D0 * AM0 * E3CM) * CM3SQI
C  XIT IS XI(T)
        XIT   = XI0 + GRLAMD*TBYMSS
C  AMPLITUDE OF SQUARED MATRIX ELEMENT (PHYS. LETT. B204 (1988) 183)
        AMPLI = (1.D0 + PARAMA*TBYMSS)**2 * ( PA + XIT*PB + XIT**2 *PC )
      ENDIF
 
C  REJECT RANDOM NUMBERS, IF RD(3) IS LARGER THAN DALITZ PLOT AMPLITUDE
      IF ( RD(3)*AMPMX .GT. AMPLI ) GOTO 100
 
      IF (DEBUG) WRITE(MDEBUG,*)'DECAY6: E3CM,E4CM,E5CM=',
     *   SNGL(E3CM),SNGL(E4CM),SNGL(E5CM)
C  P3CM, P4CM, P5CM ARE MOMENTA IN THE C.M. SYSTEM
C  P3SQ, P4SQ, P5SQ ARE SQUARED MOMENTA IN C.M. SYSTEM
      P5SQ   = E5CM**2 - CM5SQ
      P5CM   = SQRT(P5SQ)
      P4SQ   = E4CM**2 - CM4SQ
      P4CM   = SQRT(P4SQ)
      P3SQ   = E3CM**2 - CM3SQ
      P3CM   = SQRT(P3SQ)
C  ANGLE ALFA AND BETA ARE BETWEEN PARTICLE 3 AND 4 RSP. 3 AND 5
      COSALF = (P5SQ - P3SQ - P4SQ) / (2.D0 * P3CM * P4CM)
      SINALF = -SQRT(1.D0 - COSALF**2)
      COSBET = (P4SQ - P3SQ - P5SQ) / (2.D0 * P3CM * P5CM)
      SINBET = SQRT(1.D0 - COSBET**2)
C  NOW SELECT RANDOM NUMBERS FOR THREE INDEPENDENT ANGLES IN CM-SYSTEM
C  COS3CM AND PHI ARE ANGLES OF PARTICLE 3 RELATIVE TO DECAYING PARTICLE
      CALL RMMAR( RD,3,1 )
      COS3CM = 2.D0*RD(1) - 1.D0
      SIN3CM = SQRT(1.D0 - COS3CM**2)
      PHI345(1) = PI2 * RD(2)
      COSPHI = COS( PHI345(1) )
      SINPHI = SIN( PHI345(1) )
C  ANGLE PSI GIVES ROTATION OF PLANE (3,4,5) RELATIVE TO PLANE (1,3)
      PSI    = PI2 * RD(3)
      COSPSI = COS(PSI)
      SINPSI = SIN(PSI)
C  CALCULATE ALL NEEDED POLAR AND AZIMUTHAL ANGLES IN THE CM-SYSTEM
      COS4CM = COS3CM * COSALF - SIN3CM * COSPSI * SINALF
      IF ( ABS(COS4CM) .LT. 1.D0 ) THEN
         SINT4     = SQRT(1.D0 - COS4CM**2)
         SINT4I    = 1.D0 / SINT4
         AUXA      = COS3CM * COSPSI * SINALF + SIN3CM * COSALF
         COSFI4    = (COSPHI*AUXA-SINPHI*SINPSI*SINALF) * SINT4I
         PHI345(2) = ACOS( MAX( -1.D0, MIN( 1.D0, COSFI4 ) ) )
         SINFI4    = (SINPHI*AUXA+COSPHI*SINPSI*SINALF) * SINT4I
         IF ( SINFI4 .LE. 0.D0 ) PHI345(2) = PI2 - PHI345(2)
      ELSE
         PHI345(2) = 0.D0
      ENDIF
C  CALCULATE GAMMA FACTORS AND POLAR ANGLES IN LABORATORY SYSTEM
      GAM345(1) = GAMMA * (E3CM + BETA * P3CM * COS3CM) / AM3
      COS345(1) = MIN( 1.D0, (BETA * E3CM + P3CM * COS3CM) * GAMMA
     *                      / ( AM3 * SQRT(GAM345(1)**2 - 1.D0) ) )
      GAM345(2) = GAMMA * (E4CM + BETA * P4CM * COS4CM) / AM4
      COS345(2) = MIN( 1.D0, (BETA * E4CM + P4CM * COS4CM) * GAMMA
     *                       / ( AM4 * SQRT(GAM345(2)**2 - 1.D0) ) )
C  CALCULATE PARAMETERS OF PARTICLE 5, IF NEEDED
      IF ( MODE .LE. 2 ) THEN
        COS5CM      = COS3CM * COSBET - SIN3CM * COSPSI * SINBET
        IF ( ABS(COS5CM) .LT. 1.D0 ) THEN
          SINT5I    = 1.D0 / SQRT(1.D0 - COS5CM**2)
          AUXB      = COS3CM * COSPSI * SINBET + SIN3CM * COSBET
          PHI345(3) = ACOS((COSPHI*AUXB-SINPHI*SINPSI*SINBET) * SINT5I)
          SINFI5    =      (SINPHI*AUXB+COSPHI*SINPSI*SINBET) * SINT5I
          IF ( SINFI5 .LE. 0.D0 ) PHI345(3) = PI2 - PHI345(3)
        ELSE
          PHI345(3) = 0.D0
        ENDIF
        IF ( AM5 .NE. 0.D0 ) THEN
          GAM345(3) = GAMMA * (E5CM + BETA * P5CM * COS5CM) / AM5
          COS345(3) = MIN( 1.D0, (BETA * E5CM + P5CM * COS5CM) * GAMMA
     *                           / ( AM5 * SQRT(GAM345(3)**2 - 1.D0) ) )
        ELSE
C  IF PARTICLE 5 IS GAMMA RAY OR NEUTRINO, THEN GAM345(3) IS THE ENERGY
          GAM345(3) = GAMMA * (E5CM + BETA * P5CM * COS5CM)
          COS345(3) = MIN( 1.D0, (BETA * E5CM + P5CM * COS5CM) * GAMMA
     *                           / GAM345(3) )
        ENDIF
      ENDIF
 
      IF ( MODE .EQ. 3 ) THEN
C  CALCULATION OF MUON POLARIZATION. WE FOLLOW THE DESCRIPTION OF
C  L. JAUNEAU, IN: METHODS IN SUBNUCLEAR PHYSICS, VOL. 3, M. NIKOLIC ED.
C  (GORDON + BREACH, NEW YORK, 1969), P. 123
C  SEE ALSO: L.M. CHOUNET ET AL., PHYS. REP. 4 (1972) 199, APPENDIX 1.
C  SEE ALSO: N. CABBIBO, A. MAKSYMOWICZ, PHYS. LETT. 9 (1964) 352
C  (CORRECTIONS IN: PHYS. LETT. 11 (1964) 360; 14 (1965) 72)
C  WE DEFINE  BOFQ (READ: B OF Q), WHICH IS -B(Q**2)*4
        BOFQ = 1.D0 - XIT
C  ABYM AND BBYM (READ A BY M; B BY M) ARE THE QUANTITIES A/M AND B/M
        ABYM = AM0 * ( BOFQ * EPIPRM - 2.D0 * E5CM )
        BBYM = CM0SQ + 0.25D0 * CM4SQ * BOFQ**2 - BOFQ * AM0 * E4CM
C  NOW CALCULATE THE COMPONENTS APARAL (PARALLEL TO MU DIRECTION) AND
C  APERPN (PERPENDICULAR TO MU DIRECTION) USING QUANTITIES DEFINED IN
C  KAON REST SYSTEM. NOTE OUR DEFINITION OF SINALF (ALWAYS WITH NEGATIVE
C  SIGN) OPPOSITE TO CABBIBO'S SIN(PSI) AND JAUNEAU'S SIN(THETA)
        APARAL = -P3CM*AM4*BBYM*COSALF - P4CM * ( AM0*ABYM - BBYM *
     *           ( P3CM*SINALF*(E4CM-AM4)/P4CM + AM0 - E3CM ) )
        APERPN =  P3CM*AM4*BBYM*SINALF
C  NOW NORMALIZE THE PARALLEL COMPONENT OF POLARIZATION; POLART IS
C  COSINE OF THE ANGLE BETWEEN MUON MOMENTUM AND POLARISATION
        POLART = APARAL / SQRT(APARAL**2 + APERPN**2)
C  THE POLARIZATION VECTOR LIES IN THE PLANE OF MOMENTA (PION,MUON).
C  OMEGA IS THE ANGLE BY WHICH THE DECAY PLANE (PION,MUON) IS ROTATET
C  AROUND THE DIRECTION OF MUON RELATIVE TO THE PLANE (KAON,MUON)
        IF ( ABS(COS4CM) .LT. 1.D0  .AND.  SINALF .NE. 0.D0 ) THEN
          COSOME = (COS4CM*COSALF - COS3CM)*SINT4I/SINALF
          OMEGA  = ACOS( MAX( -1.D0, MIN( 1.D0, COSOME ) ) )
          IF ( SINFI4 .NE. 0.D0 ) THEN
            SINOMG = ( COSFI4 * ( COSALF - COS3CM*COS4CM ) * SINT4I
     *                 - SIN3CM * COSPHI ) / (SINALF*SINFI4)
            IF ( SINOMG .LT. 0.D0 ) OMEGA = PI2 - OMEGA
          ENDIF
        ELSE
          OMEGA = 0.D0
        ENDIF
        POLARF  = OMEGA
      ENDIF
 
      RETURN
      END
