*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:20:01  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.44  by  S.Giani
*-- Author :
*$ CREATE THREPD.FOR
*COPY THREPD
*
*=== threpd ===========================================================*
*
      SUBROUTINE THREPD(UMO,ECM1,ECM2,ECM3,PCM1,PCM2,PCM3,COD1,COF1,
     *SIF1,COD2,COF2,SIF2,COD3,COF3,SIF3,AM1,AM2,AM3)
 
#include "geant321/dblprc.inc"
#include "geant321/dimpar.inc"
#include "geant321/iounit.inc"
*
*----------------------------------------------------------------------*
*  Threpd89: slight revision by A. Ferrari                             *
*----------------------------------------------------------------------*
*
      DIMENSION F(5),XX(5)
C***THREE PARTICLE DECAY IN THE CM - SYSTEM
      COMMON /FKGAMR/ REDU,AMO,AMM(15 )
      COMMON/FKDREI/UUMO,AAM1,AAM2,AAM3,S22,UMO2,AM11,AM22,AM33,S2SUP
     *,S2SAP(2)
      COMMON/FKPRUN/ISYS
      REAL RNDM(2)
      SAVE EPS
      DATA EPS/1.D-16/
      UMOO=UMO+UMO
C***S1, S2, S3 ARE THE INVARIANT MASSES OF THE PARTICLES 1, 2, 3
C***J. VON NEUMANN - RANDOM - SELECTION OF S2
C***CALCULATION OF THE MAXIMUM OF THE S2 - DISTRIBUTION
      UUMO=UMO
      AAM1=AM1
      AAM2=AM2
      AAM3=AM3
 
      GU=(AM2+AM3)**2
      GO=(UMO-AM1)**2
      UFAK=1.0000000000001D0
      IF (GU.GT.GO) UFAK=0.9999999999999D0
      OFAK=2.D0-UFAK
      GU=GU*UFAK
      GO=GO*OFAK
      DS2=(GO-GU)/99.D0
      AM11=AM1*AM1
      AM22=AM2*AM2
      AM33=AM3*AM3
      UMO2=UMO*UMO
      RHO2=0.D0
      S22=GU
      DO 124 I=1,100
         S21=S22
         S22=GU+(I-1.D0)*DS2
         RHO1=RHO2
         RHO2=XLAMB(S22,UMO2,AM11)*XLAMB(S22,AM22,AM33)/(S22+EPS)
         IF(RHO2.LT.RHO1) GO TO 125
  124 CONTINUE
  125 S2SUP=(S22-S21)*.5D0+S21
      SUPRHO=XLAMB(S2SUP,UMO2,AM11)*XLAMB(S2SUP,AM22,AM33)/(S2SUP+EPS)
      SUPRHO=SUPRHO*1.05D0
      XO=S21-DS2
      IF (GU.LT.GO.AND.XO.LT.GU) XO=GU
      IF (GU.GT.GO.AND.XO.GT.GU) XO=GU
      XX(1)=XO
      XX(3)=S22
      X1=(XO+S22)*0.5D0
      XX(2)=X1
      F(3)=RHO2
      F(1)=XLAMB(XO,UMO2,AM11)*XLAMB(XO,AM22,AM33)/(XO+EPS)
      F(2)=XLAMB(X1,UMO2,AM11)*XLAMB(X1,AM22,AM33)/(X1+EPS)
      DO 126 I=1,16
         X4=(XX(1)+XX(2))*0.5D0
         X5=(XX(2)+XX(3))*0.5D0
         F(4)=XLAMB(X4,UMO2,AM11)*XLAMB(X4,AM22,AM33)/(X4+EPS)
         F(5)=XLAMB(X5,UMO2,AM11)*XLAMB(X5,AM22,AM33)/(X5+EPS)
         XX(4)=X4
         XX(5)=X5
         DO 128 II=1,5
            IA=II
            DO 128 III=IA,5
               IF (F (II).GE.F (III)) GO TO 128
               FH=F(II)
               F(II)=F(III)
               F(III)=FH
               FH=XX(II)
               XX(II)=XX(III)
               XX(III)=FH
128      CONTINUE
         SUPRHO=F(1)
         S2SUP=XX(1)
         DO 129 II=1,3
            IA=II
            DO 129 III=IA,3
               IF (XX(II).GE.XX(III)) GO TO 129
               FH=F(II)
               F(II)=F(III)
               F(III)=FH
               FH=XX(II)
               XX(II)=XX(III)
               XX(III)=FH
129      CONTINUE
126   CONTINUE
      AM23=(AM2+AM3)**2
      ITH=0
      REDU=2.D0
    1 CONTINUE
      ITH=ITH+1
      IF (ITH.GT.200) REDU=-9.D0
      IF (ITH.GT.200) GO TO 400
      CALL GRNDM(RNDM,2)
      C=RNDM(1)
      S2=AM23+C*((UMO-AM1)**2-AM23)
      Y=RNDM(2)
      Y=Y*SUPRHO
      RHO=XLAMB(S2,UMO2,AM11)*XLAMB(S2,AM22,AM33)/S2
      IF(Y.GT.RHO) GO TO 1
C***RANDOM SELECTION OF S3 AND CALCULATION OF S1
      CALL GRNDM(RNDM,1)
      S1=RNDM(1)
      S1=S1*RHO+AM11+AM22-(S2-UMO2+AM11)*(S2+AM22-AM33)/(2.D0*S2)-
     &RHO*.5D0
      S3=UMO2+AM11+AM22+AM33-S1-S2
      ECM1=(UMO2+AM11-S2)/UMOO
      ECM2=(UMO2+AM22-S3)/UMOO
      ECM3=(UMO2+AM33-S1)/UMOO
       PCM1=SQRT((ECM1+AM1)*(ECM1-AM1))
       PCM2=SQRT((ECM2+AM2)*(ECM2-AM2))
       PCM3=SQRT((ECM3+AM3)*(ECM3-AM3))
      CALL SFECFE(SFE,CFE)
C***TH IS THE ANGLE BETWEEN PARTICLES 1 AND 2
C***TH1, TH2 ARE THE ANGLES BETWEEN PARTICLES 1, 2 AND THE DIRECTION OF
      IF ((PCM1.LE.1.D-3).OR.(PCM2.LE.1.D-3)) GO TO 200
      COSTH=(ECM1*ECM2+0.5D0*(AM11+AM22-S1))/(PCM1*PCM2)
      GO TO 300
  200 CALL GRNDM(RNDM,1)
      UW=RNDM(1)
      COSTH=(UW-.5D0)*2.D0
 300  CONTINUE
      TMPONE = 0.9999999999999999D0
      IF(ABS(COSTH).GT.0.9999999999999999D0)
     &COSTH=SIGN(TMPONE,COSTH)
      IF (REDU.LT.1.D0) RETURN
      COSTH2=(PCM3*PCM3+PCM2*PCM2-PCM1*PCM1)/(2.D0*PCM2*PCM3)
      IF(ABS(COSTH2).GT.0.9999999999999999D0)
     &COSTH2=SIGN(TMPONE,COSTH2)
      SINTH2=SQRT(1.D0-COSTH2*COSTH2)
      SINTH1=COSTH2*SQRT(1.D0-COSTH*COSTH)-COSTH*SINTH2
      COSTH1=COSTH*COSTH2+SINTH2*SQRT(1.D0-COSTH*COSTH)
C***RANDOM SELECTION OF THE SPHERICAL COORDINATES OF THE DIRECTION OF PA
C***CFE, SFE ARE COS AND SIN OF THE ROTATION ANGLE OF THE SYSTEM 1, 2 AR
C***THE DIRECTION OF PARTICLE 3
C***CALCULATION OF THE SPHERICAL COORDINATES OF PARTICLES 1, 2
      CX11=-COSTH1
      CY11=SINTH1*CFE
      CZ11=SINTH1*SFE
      CX22=-COSTH2
      CY22=-SINTH2*CFE
      CZ22=-SINTH2*SFE
      CALL SFECFE(SIF3,COF3)
      CALL GRNDM(RNDM,1)
      COD3=RNDM(1)
      COD3=2.D0*COD3-1.D0
      SID3=SQRT(1.D0-COD3*COD3)
    2 FORMAT(5F20.15)
      COD1=CX11*COD3+CZ11*SID3
      IF((1.D0-COD1*COD1).LT.1.D-14)WRITE(ISYS,2)COD1,COF3,SID3,
     &CX11,CZ11
      SID1=SQRT(1.D0-COD1*COD1)
      COF1=(CX11*SID3*COF3-CY11*SIF3-CZ11*COD3*COF3)/SID1
      SIF1=(CX11*SID3*SIF3+CY11*COF3-CZ11*COD3*SIF3)/SID1
      COD2=CX22*COD3+CZ22*SID3
      SID2=SQRT(1.D0-COD2*COD2)
      COF2=(CX22*SID3*COF3-CY22*SIF3-CZ22*COD3*COF3)/SID2
      SIF2=(CX22*SID3*SIF3+CY22*COF3-CZ22*COD3*SIF3)/SID2
  400 RETURN
      END