GEANT321/gphys/gdeca3.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:23  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.21  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GDECA3(XM0,XM1,XM2,XM3,PCM)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *   This routine simulates three-body decays in center-of-mass.  *
  17 C.    *   Written by Jurgen Schukraft and Vincent Hedberg for R807.    *
  18 C.    *   Adapted for Geant3 by Sverker Johansson (Aug-1987).          *
  19 C.    *   (Quite extensive modifications; comments etc. are not always *
  20 C.    *    consistent with the code anymore. )                         *
  21 C.    *                                                                *
  22 C.    *   As it is, the decay is done according to pure phase-space.   *
  23 C.    *   There is code in the routine to properly treat K decays (V-A)*
  24 C.    *   but it is not used. (Requires an additional parameter        *
  25 C.    *   in the call, to tell which type of decay.)                   *
  26 C.    *                                                                *
  27 C.    *    ==>Called by : GDECAY                                       *
  28 C.    *       Author    S. Johansson  *********                        *
  29 C.    *                                                                *
  30 C.    ******************************************************************
  31 C.
  32 ***************************************************************
  33 ******* original name :                     *******************
  34 ******SUBROUTINE  D E C 3 B (ID,IP,K1,K2,K3)*******************
  35 ***************************************************************
  36 C--------------------------------------------------------
  37 C      O B S O L E T E   C O M M E N T   ! ! !
  38 C     THIS ROUTINE SIMULATES THE THREE-BODY DECAYS OF:
  39 C--------------------------------------------------------
  40 C     ID = 1 : ETA   INTO  PI0 + PI0 + PI0     PHASE SPACE
  41 C     ID = 2 : OMEGA INTO  PI0 + PI- + PI+        - " -
  42 C     ID = 3 : ETAPR INTO  PI0 + PI0 + ETA        - " -
  43 C     ID = 4 : K+/-  INTO  E+/- + PI0 + NU        (V-A)
  44 C     ID = 5 : K0L   INTO  E+/- + PI-/+ + NU      - " -
  45 C     ID = 6 : D     INTO  E   +  K  + NU      PHASE SPACE
  46 C     ID = 7 : D     INTO  E   +  K* + NU         - " -
  47 **      (** or any other decay; input is now masses
  48 **          rather than process ID   S.J.
  49 **          But the ID parameter could someday be useful
  50 **          since it controls the non-phase-space code. ***)
  51 C--------------------------------------------------------
  52 C     THE OUTPUT COMMON BLOCK /LAB/ PLAB(20,J) IS FILLED
  53 *       (** Not a common block anymore.  S.J. **)
  54 C     ACCORDING TO THE ORDER OF THE SECONDARY PARTICLES
  55 C     AS THEY APPEARE IN THE COMMENTS ABOVE. FOR EXAMPLE,
  56 C     IN THE CASE OF ID=6 (D-MESON DECAY)
  57 C     J = K1 IN PLAB(20,J) CORRESPONDS TO ELECTRON,
  58 C     J = K2  -- " -- " -- " -- " --  K-MESON,
  59 C     J = K3  -- " -- " -- " -- " --  NEUTRINO
  60 C--------------------------------------------------------
  61       DIMENSION  P(20,20)
  62       DIMENSION  HM(4,7)
  63       DIMENSION  PCM(4,3)
  64       DIMENSION RNDM(3)
  65 *
  66 ****************************************************
  67 *   Translation from new to old input parameters :
  68 *      ( S.J. )
  69 ****************************************************
  70       ID = 1
  71       HM(1,ID) = XM1
  72       HM(2,ID) = XM2
  73       HM(3,ID) = XM3
  74       HM(4,ID) = XM0
  75       K1 = 1
  76       K2 = 2
  77       K3 = 3
  78       PI2 = 2. * 3.141592654
  79 ***********************
  80 *  Original code :
  81 C--------------------------------------------------------
  82 C     SIMULATION OF SECONDARY PARTICLE MOMENTA
  83 C     IN THE REST FRAME OF PARENT PARTICLE
  84 C--------------------------------------------------------
  85       T0=HM(4,ID)-HM(1,ID)-HM(2,ID)-HM(3,ID)
  86   10  CALL GRNDM(RNDM,2)
  87       G1=RNDM(1)
  88       G2=RNDM(2)
  89       GMI=MIN(G1,G2)
  90       GMA=MAX(G1,G2)
  91       TE=GMI*T0
  92       TM=(1.-GMA)*T0
  93       TN=(GMA-GMI)*T0
  94       PA1=SQRT(TE**2+2.*HM(1,ID)*TE)
  95       PA2=SQRT(TM**2+2.*HM(2,ID)*TM)
  96       PA3=SQRT(TN**2+2.*HM(3,ID)*TN)
  97 C--------------------------------------------------------
  98 C     CHECK OF THE MOMENTUM CONSERVATION
  99 C--------------------------------------------------------
 100       PMMM=MAX(PA1,PA2)
 101       PMAX=MAX(PA3,PMMM)
 102       SP=PA1+PA2+PA3
 103       IF(PMAX.GE.SP-PMAX) GOTO 10
 104       P(4,K1)=TE+HM(1,ID)
 105       P(4,K2)=TM+HM(2,ID)
 106       P(4,K3)=TN+HM(3,ID)
 107 C--------------------------------------------------------
 108 C     ACCOUNTING FOR MATRIX ELEMENT (FOR ID=4 AND 5)
 109 *     (** Calculates correct (V-A) decay for kaons.       **)
 110 *     (** Now inactive.  Change ID to activate it. (S.J)  **)
 111 C--------------------------------------------------------
 112       IF(ID.NE.4.AND.ID.NE.5) GOTO 20
 113       EEM=((HM(4,ID)-HM(2,ID))**2+HM(1,ID)**2-
 114      *      HM(3,ID)**2)/(2.*(HM(4,ID)-HM(2,ID)))
 115       HMMAX=EEM*(HM(4,ID)-HM(2,ID)-EEM)+
 116      *      0.5*(EEM**2-HM(1,ID)**2)
 117       HMA=P(4,K1)*P(4,K3)+0.25*(PA1**2+PA3**2-PA2**2)
 118       CALL GRNDM(RNDM,1)
 119       GG=RNDM(1)
 120       IF(HMA.LT.GG*HMMAX) GOTO 10
 121 C--------------------------------------------------------
 122 C     CALCULATIONS OF MOMENTA
 123 C--------------------------------------------------------
 124 20    CONTINUE
 125       CALL GRNDM(RNDM,3)
 126       CTE=2.*RNDM(1)-1.
 127       STE=SQRT(ABS(1.-CTE**2))
 128       FE =PI2*RNDM(2)
 129       CFE=COS(FE)
 130       SFE=SIN(FE)
 131       P(1,K1)=PA1*STE*CFE
 132       P(2,K1)=PA1*STE*SFE
 133       P(3,K1)=PA1*CTE
 134       CTEN=(PA2**2-PA1**2-PA3**2)/(2.*PA1*PA3)
 135       STEN=SQRT(ABS(1.-CTEN**2))
 136       FEN =PI2*RNDM(3)
 137       CFEN=COS(FEN)
 138       SFEN=SIN(FEN)
 139       P(1,K3)=PA3*(STEN*CFEN*CTE*CFE-STEN*SFEN*SFE+CTEN*STE*CFE)
 140       P(2,K3)=PA3*(STEN*CFEN*CTE*SFE+STEN*SFEN*CFE+CTEN*STE*SFE)
 141       P(3,K3)=PA3*(-STEN*CFEN*STE+CTEN*CTE)
 142       P(1,K2)=-P(1,K1)-P(1,K3)
 143       P(2,K2)=-P(2,K1)-P(2,K3)
 144       P(3,K2)=-P(3,K1)-P(3,K3)
 145 C-------------------------------------------------------
 146 ***   Re-translation from old to new output arrays : ***
 147 ***     (S.J.)                                       ***
 148 ***--------------------------------------------------***
 149       DO 37 I=1,4
 150         DO 38 J=1,3
 151           PCM(I,J) = P(I,J)
 152  38     CONTINUE
 153  37   CONTINUE
 154       END