C------------ Decay --------------------------------------
C LAST DATE OF CHANGE 20.03.08
      SUBROUTINE DECAY(IHAD)
      COMMON/FINPAR/PXF(10000),PYF(10000),PZF(10000),HEF(10000),
     *AMF(10000),ICHF(10000),IBARF(10000),ANF(10000),NREF(10000)

      COMMON/DECAYC/ ZKNAME(533),NZK(533,3),WT(533)

      COMMON/PART/ANAME(180),AM(180),GA(180),TAU(180),ICH(180),IBAR(180)
     *,K1(180),K2(180)

      COMMON/METLSP/ IS,ITS(100000),CXS(100000),CYS(100000),CZS(100000),
     *ELS(100000),PLS(100000)

      COMMON/DREI/ TEST(12)

      COMMON/GAMRED/REDU,AMO,AMM(15)

      COMMON/PRINT/ISYS

      COMMON/LIMMAS/IDSTAB(180), SUMKM(533), AML(180), FI0ML(180)

C
      DIMENSION AMS(100000), WCHAN(20), NUMCH(20)
C
      REAL*8 ANF
      REAL*8 ZKNAME
      REAL*8 ANAME
C
      INTEGER VW
      VW=-1
      REDU=2.
      DO 10   I=1,IHAD
      ITS(I)=NREF(I)
      PLS(I)=SQRT(PXF(I)**2+PYF(I)**2+PZF(I)**2)
C
      CXS(I)=0.
      CYS(I)=0.
      CZS(I)=1.
C
      IF(PLS(I).NE.0.)CXS(I)=PXF(I)/PLS(I)
      IF(PLS(I).NE.0.)CYS(I)=PYF(I)/PLS(I)
      IF(PLS(I).NE.0.)CZS(I)=PZF(I)/PLS(I)
      ELS(I)=HEF(I)
      AMS(I)=AMF(I)
 10   CONTINUE
      IST=IHAD
      IR=0
 20   CONTINUE
C
C*****TEST STABLE OR UNSTABLE
C$$$$$ISTAB=1/2/3 MEAN ... STRONG + WEAK DECAY / ONLY STRONG DECAYS /
C*****STRONG DECAY + WEAK DECAYS FOR CHARMED PARTICLES AND TAU LEPTONS
C
      IDPART=ITS(IST)
C
      IF(IDSTAB(IDPART).EQ.1)   GO TO 60
C
      GO TO 70
C
C         STORING OF THE STABLE PARTICLES
C
 60   IR=IR+1
      IF(IR.LE.100000) GO TO 62
C
*      WRITE(ISYS,61)
* 61   FORMAT(2X,' NUMBER OF STAB. FINAL PART. IS GREATER THAN 100000'/,
*     ,2X,'   ONLY 100000 PARTICLES ARE STORED IN COMMON/FINPAR/'/
*     ,2X,'   THE OTHERS ARE LOSSED.      DECAY')
      GO TO 180
C
 62   CONTINUE
      NREF(IR)=ITS(IST)
      ITT=ITS(IST)
      AMF(IR)=AMS(IST)
      ANF(IR)=ANAME(ITT)
      ICHF(IR)=ICH(ITT)
      IBARF(IR)=IBAR(ITT)
      HEF(IR)=ELS(IST)
      PXF(IR)=CXS(IST)*PLS(IST)
      PYF(IR)=CYS(IST)*PLS(IST)
      PZF(IR)=CZS(IST)*PLS(IST)
      IST=IST-1
      IF(IST.GE.1) GO TO 20
      GO TO 180
C
C  SIMULATION OF DECAY OF PARTICLE WITH NUMBER --- IST
C
 70   IT=ITS(IST)
      GAM=ELS(IST)/AMS(IST)
      BGAM=PLS(IST)/AMS(IST)
      ECO=AMS(IST)
      KZ1=K1(IT)
      KZ2=K2(IT)
C
      SUMWT=0.
      NCHAN=0
      DO 80   I=KZ1, KZ2
      IF((SUMKM(I).GT.ECO).OR.(WT(I).EQ.0.)) GO TO 80
      NCHAN=NCHAN+1
      SUMWT=SUMWT+WT(I)
      WCHAN(NCHAN)=SUMWT
      NUMCH(NCHAN)=I
 80   CONTINUE
C
      DO 90   I=1,NCHAN
 90   WCHAN(I)=WCHAN(I)/SUMWT
C
      IREPIT=0
 100  CONTINUE
      IREPIT=IREPIT+1
C
C         CHOISE OF A DECAY CHANNEL
C
      IF(IREPIT.LE.100) GO TO 110
      ITTT=ITS(IST)
      IRRR=IR+1
*      PRINT 101,IRRR,AMS(IST),ANAME(ITTT),AM(ITTT),ICH(ITTT),IBAR(ITTT),
*     ,          ITTT
 101  FORMAT(' A PARTICLE STORED IN COMMON/FINPAR/ WITH NUMBER',
     , I6,' AND MASS ',E11.4,/1X,
     ,'DID NOT DECAY OR DECAY CHANNEL WAS NOT DEFINED'/,15X,
     ,'IT IS CHARACTERIZED BY'/,
     ,'  NAME  TAB. MASS     CHARGE BAR. NUMBER INDEX'/,
     ,1X,A6,1X,E11.4,3X,I4,5X,I4,6X,I4,'       DECAY'/)
      GO TO 60
 110  VV=RNDM1(VW)-1.E-17
C
      DO 120  I=1,NCHAN
      IF(VV.LE.WCHAN(I)) GO TO 130
 120  CONTINUE
 130  IIK=NUMCH(I)
C
C   IIK IS THE DECAY CHANNEL
C
      IT1=NZK(IIK,1)
      IT2=NZK(IIK,2)
      IT3=NZK(IIK,3)
C
C  IT1,IT2, IT3 ARE THE PRODUCED PARTICLES FROM IT
C
      IF (IT2-1. GE.0) GO TO 140
C
C         FOR ONE-PARTICLE DECAY
C
      AM1=ECO
      GO TO 160
C
 140  CONTINUE
      IF(IT3.NE.0) GO TO 150
C
C         FOR TWO-PARTICLES DECAY
C
      AM1=AMASS(IT1)
      AM2=AMASS(IT2)
      IF(AM1+AM2.GE.ECO) GO TO 100
      CALL TWOPAD(ECO,ECM1,ECM2,PCM1,PCM2,COD1,COF1,SIF1,COD2,COF2,SIF2,
     * AM1,AM2)
C
C     IF(IBAR(IT).EQ.0)                     GO TO 160
C
C     IF(.NOT.((IBAR(IT1).NE.0))) GO TO 160
C
C     AEPS=RNDM1(-1)
C     IF(RNDM1(-1).LE.0.5) THEN
C       COD1=SQRT(AEPS)
C     ELSE
C       COD1=-SQRT(AEPS)
C     ENDIF
C      COD1=0.
C      COD2=0.
C
C     COD2=-COD1
      GO TO 160
C
 150  CONTINUE
C
C        FOR THRE-PARTICLES DECAY
C
      AM1=AMASS(IT1)
      AM2=AMASS(IT2)
      AM3=AMASS(IT3)
      IF(AM1+AM2+AM3.GE.ECO) GO TO 100
      CALL THREPD(ECO,ECM1,ECM2,ECM3,PCM1,PCM2,PCM3,COD1,COF1,SIF1,
     *COD2,COF2,SIF2,COD3,COF3,SIF3,AM1,AM2,AM3)
C
 160  CONTINUE
      ITS(IST)=IT1
      AMS(IST)=AM1
      IF (IT2-1.LT.0) GO TO 170
      ITS(IST+1)  =IT2
      RX=CXS(IST)
      RY=CYS(IST)
      RZ=CZS(IST)
      CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD1,COF1,SIF1,PCM1,ECM1,
     *PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
      IST=IST+1
      AMS(IST)=AM2
      CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD2,COF2,SIF2,PCM2,ECM2,
     *PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
      IF (IT3.LE.0) GO TO 170
      IST=IST+1
      ITS(IST)=IT3
      AMS(IST)=AM3
      CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD3,COF3,SIF3,PCM3,ECM3,
     *PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
 170  CONTINUE
      GO TO 20
 180  CONTINUE
      IHAD=IR
      RETURN
      END

      SUBROUTINE TWOPAD(UMO,ECM1,ECM2,PCM1,PCM2,COD1,COF1,SIF1,
     *COD2,COF2,SIF2,AM1,AM2)
      COMMON/GAMRED/REDU,AMO,AMM(15)

      INTEGER X
      X=-1
      UMO2=UMO**2
      AM11=AM1**2
      AM22=AM2**2
      ECM1=(UMO2+AM11-AM22)/(2.*UMO)
      ECM2=UMO-ECM1
      WAU=ECM1**2-AM11
      IF(WAU.LT.0)REDU=REDU-5.
      WAU=ABS(WAU)
      PCM1=SQRT(WAU)
      PCM2=PCM1
      CALL SFECFE(SIF1,COF1)
      COD1=2.*RNDM1(X)-1.
      COD2=-COD1
      COF2=-COF1
      SIF2=-SIF1
      RETURN
      END

      SUBROUTINE THREPD (UMO,ECM1,ECM2,ECM3,PCM1,PCM2,PCM3,COD1,COF1,
     *SIF1,COD2,COF2,SIF2,COD3,COF3,SIF3,AM1,AM2,AM3)
      DIMENSION F(5),XX(5)
      COMMON/GAMRED/REDU,AMO,AMM(15)

      COMMON/DREI/UUMO,AAM1,AAM2,AAM3,S22,UMO2,AM11,AM22,AM33,S2SUP,
     *S2SAP(2)

      COMMON/PRINT/ISYS

      INTEGER V
      DATA EPS/1.E-14/
      SAVE EPS
      V=-1
      UMOO=UMO+UMO
      UUMO=UMO
      AAM1=AM1
      AAM2=AM2
      AAM3=AM3
      GU=(AM2+AM3)**2
      GO=(UMO-AM1)**2
      UFAK=1.00001
      IF(GU.GT.GO) UFAK=0.99999
      OFAK=2.-UFAK
      GU=GU*UFAK
      GO=GO*OFAK
      DS2=(GO-GU)/99.
      AM11=AM1*AM1
      AM22=AM2*AM2
      AM33=AM3*AM3
      UMO2=UMO*UMO
      RHO2=0
      S22=GU
      DO 124 I=1,100
      S21=S22
      S22=GU+(I-1)*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)/2.+S21
      SUPRHO=XLAMB(S2SUP,UMO2,AM11)*XLAMB(S2SUP,AM22,AM33)/(S2SUP+EPS)
      SUPRHO=SUPRHO*1.05
      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.5
      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.5
      X5=(XX(2)+XX(3))*0.5
      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.
    1 CONTINUE
      ITH=ITH+1
      IF(ITH.GT.200)REDU=-9
      IF(ITH.GT.200) GO TO 400
      C=RNDM1(V)
      S2=AM23+C*((UMO-AM1)**2-AM23)
      Y=RNDM1(V)
      Y=Y*SUPRHO
      RHO=XLAMB(S2,UMO2,AM11)*XLAMB(S2,AM22,AM33)/S2
      IF(Y.GT.RHO) GO TO 1
      S1=RNDM1(V)
      S1=S1*RHO+AM11+AM22-(S2-UMO2+AM11)*(S2+AM22-AM33)/(2.*S2)-RHO/2.
      S3=UMO2+AM11+AM22+AM33-S1-S2
      ECM1=(UMO2+AM11-S2)/UMOO+1.E-6
      ECM2=(UMO2+AM22-S3)/UMOO+1.E-6
      ECM3=(UMO2+AM33-S1)/UMOO+1.E-6
      PCM1=SQRT((ECM1+AM1)*(ECM1-AM1))
      PCM2=SQRT((ECM2+AM2)*(ECM2-AM2))
      PCM3=SQRT((ECM3+AM3)*(ECM3-AM3))
      CALL SFECFE(SFE,CFE)
      IF((PCM1.LE.1.E-3).OR.(PCM2.LE.1.E-3)) GO TO 200
      COSTH=(ECM1*ECM2+0.5*(AM11+AM22-S1))/(PCM1*PCM2)
      GO TO 300
  200 UW=RNDM1(V)
      COSTH=(UW-0.5)*2.
  300 CONTINUE
      IF(ABS(COSTH).GT.0.99999) COSTH=SIGN(0.99999,COSTH)
      IF(REDU.LT.1.) RETURN
      COSTH2=(PCM3*PCM3+PCM2*PCM2-PCM1*PCM1)/(2.*PCM2*PCM3)
      IF(ABS(COSTH2).GT.0.99999) COSTH2=SIGN(0.99999,COSTH2)
      SINTH2=SQRT(1.-COSTH2**2)
      SINTH1=COSTH2*SQRT(1.-COSTH**2)-COSTH*SINTH2
      COSTH1=COSTH*COSTH2+SINTH2*SQRT(1.-COSTH**2)
      CX11=-COSTH1
      CY11=SINTH1*CFE
      CZ11=SINTH1*SFE
      CX22=-COSTH2
      CY22=-SINTH2*CFE
      CZ22=-SINTH2*SFE
      CALL SFECFE(SIF3,COF3)
      COD3=RNDM1(V)
      COD3=2.*COD3-1.
      SID3=SQRT(1.-COD3**2)
    2 FORMAT(5F20.15)
      COD1=CX11*COD3+CZ11*SID3
*      IF((1.-COD1*COD1).LT.1.E-14) WRITE(ISYS,2) COD1,COF3,SID3,CX11,CZ1
*     *1
      SID1=SQRT(1.-COD1**2)                                  ! Uzhi ???
      if(SID1.ge.1.0e-9) then
      COF1=(CX11*SID3*COF3-CY11*SIF3-CZ11*COD3*COF3)/SID1
      SIF1=(CX11*SID3*SIF3+CY11*COF3-CZ11*COD3*SIF3)/SID1
      else
        COF1=1.
        SIF1=0.
      endif
      COD2=CX22*COD3+CZ22*SID3
      SID2=SQRT(1.-COD2**2)
      if(SID2.ge.1.0e-9) then
      COF2=(CX22*SID3*COF3-CY22*SIF3-CZ22*COD3*COF3)/SID2
      SIF2=(CX22*SID3*SIF3+CY22*COF3-CZ22*COD3*SIF3)/SID2
      else
        COF2=1.
        SIF2=0.
      endif
  400 RETURN
      END

      FUNCTION XLAMB(X,Y,Z)
      COMMON/IDGB/IDGB,IDG

      COMMON/GAMRED/REDU,AMO,AMM(15)

      COMMON/DREI/TEST(12)

      COMMON/PRINT/ISYS

      YZ=Y-Z
      XLAMB=X*X-2.*X*(Y+Z)+YZ**2
      XLAM=XLAMB
      IF(IDGB.LE.0) GO TO 11
      IF(XLAM .GT.1.E-12) GO TO 11
*      WRITE(ISYS,12)
*      WRITE(ISYS,10) XLAM,X,Y,Z,TEST
*      WRITE(ISYS,13)
   12 FORMAT(/,10X,12H XLAMB PRINT)
   13 FORMAT(10X,60(1H*))
   10 FORMAT(4E20.8,5HXLAMB,/,12F10.5)
   11 CONTINUE
      IF(XLAM .LE.0) XLAMB=ABS(XLAMB)
      XLAMB=SQRT(XLAMB)
      RETURN
      END

      SUBROUTINE TRAFO(GAM,BGAM,CX,CY,CZ,COD,COF,SIF,P,ECM,
     1PL,CXL,CYL,CZL,EL)
C    LORENTZ TRANSFORMATION INTO THE LAB -SYSTEM
      SID=SQRT(1.-COD*COD)
      PLX=P*SID*COF
      PLY=P*SID*SIF
      PCMZ=P*COD
      PLZ=GAM*PCMZ+BGAM*ECM
      PL=SQRT(PLX*PLX+PLY*PLY+PLZ*PLZ)
      EL=GAM*ECM+BGAM*PCMZ
C    ROTATION INTO THE ORIGINAL DIRECTION
      COZ=1.
      IF(PL.GE.1.0E-6) COZ=PLZ/PL
      SIZ=SQRT(1.-COZ**2)
      CALL TRANS(CX,CY,CZ,COZ,SIZ,SIF,COF,CXL,CYL,CZL)
      RETURN
      END

      SUBROUTINE TRANS(XO,YO,ZO,CDE,SDE,SFE,CFE,X,Y,Z)
      IF (ABS(XO)-0.0001) 1,1,2
 1    IF(ABS(YO)-0.0001) 3,3,2
 3    CONTINUE
      X=SDE*CFE
      Y=SDE*SFE
      Z=CDE
      RETURN
 2    CONTINUE
      XI=SDE*CFE
      YI=SDE*SFE
      ZI=CDE
      A=SQRT(XO**2+YO**2)
      X=-YO*XI/A-ZO*XO*YI/A+XO*ZI
      Y=XO*XI/A-ZO*YO*YI/A+YO*ZI
      Z=A*YI+ZO*ZI
      RETURN
      END