GEANT321/fluka/threpd.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:01  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.44  by  S.Giani
  11 *-- Author :
  12 *$ CREATE THREPD.FOR
  13 *COPY THREPD
  14 *
  15 *=== threpd ===========================================================*
  16 *
  17       SUBROUTINE THREPD(UMO,ECM1,ECM2,ECM3,PCM1,PCM2,PCM3,COD1,COF1,
  18      *SIF1,COD2,COF2,SIF2,COD3,COF3,SIF3,AM1,AM2,AM3)
  19
  20 #include "geant321/dblprc.inc"
  21 #include "geant321/dimpar.inc"
  22 #include "geant321/iounit.inc"
  23 *
  24 *----------------------------------------------------------------------*
  25 *  Threpd89: slight revision by A. Ferrari                             *
  26 *----------------------------------------------------------------------*
  27 *
  28       DIMENSION F(5),XX(5)
  29 C***THREE PARTICLE DECAY IN THE CM - SYSTEM
  30       COMMON /FKGAMR/ REDU,AMO,AMM(15 )
  31       COMMON/FKDREI/UUMO,AAM1,AAM2,AAM3,S22,UMO2,AM11,AM22,AM33,S2SUP
  32      *,S2SAP(2)
  33       COMMON/FKPRUN/ISYS
  34       REAL RNDM(2)
  35       SAVE EPS
  36       DATA EPS/1.D-16/
  37       UMOO=UMO+UMO
  38 C***S1, S2, S3 ARE THE INVARIANT MASSES OF THE PARTICLES 1, 2, 3
  39 C***J. VON NEUMANN - RANDOM - SELECTION OF S2
  40 C***CALCULATION OF THE MAXIMUM OF THE S2 - DISTRIBUTION
  41       UUMO=UMO
  42       AAM1=AM1
  43       AAM2=AM2
  44       AAM3=AM3
  45
  46       GU=(AM2+AM3)**2
  47       GO=(UMO-AM1)**2
  48       UFAK=1.0000000000001D0
  49       IF (GU.GT.GO) UFAK=0.9999999999999D0
  50       OFAK=2.D0-UFAK
  51       GU=GU*UFAK
  52       GO=GO*OFAK
  53       DS2=(GO-GU)/99.D0
  54       AM11=AM1*AM1
  55       AM22=AM2*AM2
  56       AM33=AM3*AM3
  57       UMO2=UMO*UMO
  58       RHO2=0.D0
  59       S22=GU
  60       DO 124 I=1,100
  61          S21=S22
  62          S22=GU+(I-1.D0)*DS2
  63          RHO1=RHO2
  64          RHO2=XLAMB(S22,UMO2,AM11)*XLAMB(S22,AM22,AM33)/(S22+EPS)
  65          IF(RHO2.LT.RHO1) GO TO 125
  66   124 CONTINUE
  67   125 S2SUP=(S22-S21)*.5D0+S21
  68       SUPRHO=XLAMB(S2SUP,UMO2,AM11)*XLAMB(S2SUP,AM22,AM33)/(S2SUP+EPS)
  69       SUPRHO=SUPRHO*1.05D0
  70       XO=S21-DS2
  71       IF (GU.LT.GO.AND.XO.LT.GU) XO=GU
  72       IF (GU.GT.GO.AND.XO.GT.GU) XO=GU
  73       XX(1)=XO
  74       XX(3)=S22
  75       X1=(XO+S22)*0.5D0
  76       XX(2)=X1
  77       F(3)=RHO2
  78       F(1)=XLAMB(XO,UMO2,AM11)*XLAMB(XO,AM22,AM33)/(XO+EPS)
  79       F(2)=XLAMB(X1,UMO2,AM11)*XLAMB(X1,AM22,AM33)/(X1+EPS)
  80       DO 126 I=1,16
  81          X4=(XX(1)+XX(2))*0.5D0
  82          X5=(XX(2)+XX(3))*0.5D0
  83          F(4)=XLAMB(X4,UMO2,AM11)*XLAMB(X4,AM22,AM33)/(X4+EPS)
  84          F(5)=XLAMB(X5,UMO2,AM11)*XLAMB(X5,AM22,AM33)/(X5+EPS)
  85          XX(4)=X4
  86          XX(5)=X5
  87          DO 128 II=1,5
  88             IA=II
  89             DO 128 III=IA,5
  90                IF (F (II).GE.F (III)) GO TO 128
  91                FH=F(II)
  92                F(II)=F(III)
  93                F(III)=FH
  94                FH=XX(II)
  95                XX(II)=XX(III)
  96                XX(III)=FH
  97 128      CONTINUE
  98          SUPRHO=F(1)
  99          S2SUP=XX(1)
 100          DO 129 II=1,3
 101             IA=II
 102             DO 129 III=IA,3
 103                IF (XX(II).GE.XX(III)) GO TO 129
 104                FH=F(II)
 105                F(II)=F(III)
 106                F(III)=FH
 107                FH=XX(II)
 108                XX(II)=XX(III)
 109                XX(III)=FH
 110 129      CONTINUE
 111 126   CONTINUE
 112       AM23=(AM2+AM3)**2
 113       ITH=0
 114       REDU=2.D0
 115     1 CONTINUE
 116       ITH=ITH+1
 117       IF (ITH.GT.200) REDU=-9.D0
 118       IF (ITH.GT.200) GO TO 400
 119       CALL GRNDM(RNDM,2)
 120       C=RNDM(1)
 121       S2=AM23+C*((UMO-AM1)**2-AM23)
 122       Y=RNDM(2)
 123       Y=Y*SUPRHO
 124       RHO=XLAMB(S2,UMO2,AM11)*XLAMB(S2,AM22,AM33)/S2
 125       IF(Y.GT.RHO) GO TO 1
 126 C***RANDOM SELECTION OF S3 AND CALCULATION OF S1
 127       CALL GRNDM(RNDM,1)
 128       S1=RNDM(1)
 129       S1=S1*RHO+AM11+AM22-(S2-UMO2+AM11)*(S2+AM22-AM33)/(2.D0*S2)-
 130      &RHO*.5D0
 131       S3=UMO2+AM11+AM22+AM33-S1-S2
 132       ECM1=(UMO2+AM11-S2)/UMOO
 133       ECM2=(UMO2+AM22-S3)/UMOO
 134       ECM3=(UMO2+AM33-S1)/UMOO
 135        PCM1=SQRT((ECM1+AM1)*(ECM1-AM1))
 136        PCM2=SQRT((ECM2+AM2)*(ECM2-AM2))
 137        PCM3=SQRT((ECM3+AM3)*(ECM3-AM3))
 138       CALL SFECFE(SFE,CFE)
 139 C***TH IS THE ANGLE BETWEEN PARTICLES 1 AND 2
 140 C***TH1, TH2 ARE THE ANGLES BETWEEN PARTICLES 1, 2 AND THE DIRECTION OF
 141       IF ((PCM1.LE.1.D-3).OR.(PCM2.LE.1.D-3)) GO TO 200
 142       COSTH=(ECM1*ECM2+0.5D0*(AM11+AM22-S1))/(PCM1*PCM2)
 143       GO TO 300
 144   200 CALL GRNDM(RNDM,1)
 145       UW=RNDM(1)
 146       COSTH=(UW-.5D0)*2.D0
 147  300  CONTINUE
 148       TMPONE = 0.9999999999999999D0
 149       IF(ABS(COSTH).GT.0.9999999999999999D0)
 150      &COSTH=SIGN(TMPONE,COSTH)
 151       IF (REDU.LT.1.D0) RETURN
 152       COSTH2=(PCM3*PCM3+PCM2*PCM2-PCM1*PCM1)/(2.D0*PCM2*PCM3)
 153       IF(ABS(COSTH2).GT.0.9999999999999999D0)
 154      &COSTH2=SIGN(TMPONE,COSTH2)
 155       SINTH2=SQRT(1.D0-COSTH2*COSTH2)
 156       SINTH1=COSTH2*SQRT(1.D0-COSTH*COSTH)-COSTH*SINTH2
 157       COSTH1=COSTH*COSTH2+SINTH2*SQRT(1.D0-COSTH*COSTH)
 158 C***RANDOM SELECTION OF THE SPHERICAL COORDINATES OF THE DIRECTION OF PA
 159 C***CFE, SFE ARE COS AND SIN OF THE ROTATION ANGLE OF THE SYSTEM 1, 2 AR
 160 C***THE DIRECTION OF PARTICLE 3
 161 C***CALCULATION OF THE SPHERICAL COORDINATES OF PARTICLES 1, 2
 162       CX11=-COSTH1
 163       CY11=SINTH1*CFE
 164       CZ11=SINTH1*SFE
 165       CX22=-COSTH2
 166       CY22=-SINTH2*CFE
 167       CZ22=-SINTH2*SFE
 168       CALL SFECFE(SIF3,COF3)
 169       CALL GRNDM(RNDM,1)
 170       COD3=RNDM(1)
 171       COD3=2.D0*COD3-1.D0
 172       SID3=SQRT(1.D0-COD3*COD3)
 173     2 FORMAT(5F20.15)
 174       COD1=CX11*COD3+CZ11*SID3
 175       IF((1.D0-COD1*COD1).LT.1.D-14)WRITE(ISYS,2)COD1,COF3,SID3,
 176      &CX11,CZ11
 177       SID1=SQRT(1.D0-COD1*COD1)
 178       COF1=(CX11*SID3*COF3-CY11*SIF3-CZ11*COD3*COF3)/SID1
 179       SIF1=(CX11*SID3*SIF3+CY11*COF3-CZ11*COD3*SIF3)/SID1
 180       COD2=CX22*COD3+CZ22*SID3
 181       SID2=SQRT(1.D0-COD2*COD2)
 182       COF2=(CX22*SID3*COF3-CY22*SIF3-CZ22*COD3*COF3)/SID2
 183       SIF2=(CX22*SID3*SIF3+CY22*COF3-CZ22*COD3*SIF3)/SID2
 184   400 RETURN
 185       END