GEANT321/gphys/glando.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:25  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 GLANDO(IMODE,STEP,Z,A,RHO,P,E,XMASS,DE,IFLAG)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *                                                                *
  17 C.    *  GLANDO                                                        *
  18 C.    *  ------                                                        *
  19 C.    *                                                                *
  20 C.    *  Energy straggling using Gaussian, Landau & Vavilov theories.  *
  21 C.    *                                                                *
  22 C.    *  Input                                                         *
  23 C.    *  -----                                                         *
  24 C.    *  IMODE  =  3   Landau(RANLAN)  sampling                        *
  25 C.    *         =  4   Landau(GENLAN)  sampling                        *
  26 C.    *         =  5   Vavilov(DINVAV) sampling                        *
  27 C.    *         =  6   Gaussian sampling                               *
  28 C.    *         =  2   Automatic selection of relevant distribution    *
  29 C.    *                                                                *
  30 C.    *  STEP   =  current step-length (cm)                            *
  31 C.    *                                                                *
  32 C.    *  Output                                                        *
  33 C.    *  ------                                                        *
  34 C.    *  DE     =  DE/DX - <DE/DX>     (GeV)                           *
  35 C.    *                                                                *
  36 C.    *  IFLAG  =  3   Landau(RANLAN) sampling used                    *
  37 C.    *         =  4   Landau(GENLAN) sampling used                    *
  38 C.    *         =  5   Vavilov(GVAVIV) sampling used                   *
  39 C.    *         =  6   Gaussian sampling used                          *
  40 C.    *                                                                *
  41 C.    *  Warning                                                       *
  42 C.    *  -------                                                       *
  43 C.    *  Only Landau sampling should be used since this has been well  *
  44 C.    *  tested whereas both Vavilov and Gaussian sampling are being   *
  45 C.    *  developed.                                                    *
  46 C.    *                                                                *
  47 C.    *  Author      : G.N. Patrick                                    *
  48 C.    *  Date        : 03.05.1985                                      *
  49 C.    *  Last update : 09.09.1985                                      *
  50 C.    *                                                                *
  51 C.    ******************************************************************
  52 C.
  53       PARAMETER (EULER=0.577215)
  54       PARAMETER (P1=.60715,P2=.67794,P3=.52382E-1,P4=.94753,
  55      +           P5=.74442,P6=1.1934)
  56 #include "geant321/gconsp.inc"
  57       DIMENSION RNDM(2)
  58       FLAND(X) = P1+P6*X+(P2+P3*X)*EXP(P4*X+P5)
  59 C.
  60 C.    ------------------------------------------------------------------
  61 C.
  62       IF(STEP.LT.1.E-7)THEN
  63          DE=0.
  64          IFLAG=0
  65          RETURN
  66       ENDIF
  67 C
  68 C     Calculate xi factor (KeV).
  69 C
  70       BETA   = P/E
  71       GAMMA  = E/XMASS
  72       XI     = (153.5*Z*STEP*RHO)/(A*BETA*BETA)
  73 C
  74 C     Maximum energy transfer to atomic electron (KeV).
  75 C
  76       ETA    = BETA*GAMMA
  77       ETASQ  = ETA*ETA
  78       RATIO  = EMASS/XMASS
  79       F1     = 2.*EMASS*ETASQ
  80       F2     = 1.+2.*RATIO*GAMMA+RATIO*RATIO
  81       EMAX   = F1*1.E+6/F2
  82 C
  83 C     Calculate Kappa significance ratio.
  84 C
  85       CAPPA  = XI/EMAX
  86 C
  87 C     Choose correct function if IMODE set to automatic selection.
  88 C
  89       IMODEI = IMODE
  90       IF (IMODEI.EQ.2)    THEN
  91          IF (CAPPA.LT.0.01)                              IMODEI = 4
  92 *****    IF (CAPPA.GE.0.01.AND.CAPPA.LE.10.)             IMODEI = 3
  93          IF (CAPPA.GE.0.01.AND.CAPPA.LE.10.)             IMODEI = 5
  94          IF (CAPPA.GT.10.)                               IMODEI = 6
  95       ENDIF
  96 C
  97 C     +---------------------------------------------------------------+
  98 C     I Sample lambda variable from Kolbig/Schorr Landau distribution I
  99 C     +---------------------------------------------------------------+
 100 C
 101    10 CONTINUE
 102       IF (IMODEI.EQ.3) THEN
 103          XMEAN = -BETA**2-LOG(XI/EMAX)+EULER-1.
 104          XLAMX = FLAND(XMEAN)
 105    20    CALL GRNDM(RNDM,1)
 106          IF( RNDM(1) .GT. 0.980 ) GO TO 20
 107          XLAMB = GLANDR(RNDM(1))
 108          IF(XLAMB.GT.XLAMX) GO TO 20
 109          IFLAG  = 3
 110 C
 111 C            +----------------------------------------------+
 112 C            I  Sample lambda variable from James & Hancock I
 113 C            I  Landau distribution                         I
 114 C            +----------------------------------------------+
 115 C
 116
 117       ELSEIF (IMODEI.EQ.4) THEN
 118          XMEAN = -BETA**2-LOG(XI/EMAX)+EULER-1.
 119          XLAMX = FLAND(XMEAN)
 120   25     CALL GLANDG(XLAMB)
 121          IF(XLAMB.GT.XLAMX) GO TO 25
 122          IFLAG = 4
 123 C
 124 C            +---------------------------------------------------------+
 125 C            I Sample lambda variable (Landau not Vavilov) from        I
 126 C            I Rotondi&Montagna&Kolbig Vavilov distribution            I
 127 C            +---------------------------------------------------------+
 128 C
 129       ELSEIF (IMODEI.EQ.5) THEN
 130 C
 131 C            Keep within computation bounds.
 132 C
 133          RKA = CAPPA
 134          BE2 = BETA*BETA
 135          CALL GRNDM(RNDM,1)
 136          XLAMB = GVAVIV(RKA,BE2,RNDM(1))
 137          IFLAG = 5
 138 C
 139 C     +-----------------------------------+
 140 C     I Sample from Gaussian distribution I
 141 C     +-----------------------------------+
 142 C
 143       ELSEIF (IMODEI.EQ.6) THEN
 144          SIGMA  = XI*EMAX*(1.-(BETA*BETA/2.))
 145          SIGMA  = SQRT(SIGMA)
 146    30    CALL GRNDM(RNDM,2)
 147          IF(RNDM(1).LE.0.)GO TO 30
 148          F1     = -2.*LOG(RNDM(1))
 149          DEKEV  = SIGMA*SQRT(F1)*COS(2.*PI*RNDM(2))
 150          IFLAG  = 6
 151          GOTO 40
 152       ENDIF
 153 C
 154 C     Calculate DE/DX - <DE/DX>
 155 C
 156       DEKEV  = XI*(XLAMB+BETA*BETA+LOG(XI/EMAX)-EULER+1.)
 157 C
 158    40 DE=DEKEV*1.E-6
 159       END