GEANT321/gphys/gfluct.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:24  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 GFLUCT(DEMEAN,DE)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *   Subroutine to decide which method is used to simulate        *
  17 C.    *   the straggling around the mean energy loss.                  *
  18 C.    *                                                                *
  19 C.    *                                                                *
  20 C.    *   DNMIN:  <---------->1<-------->30<--------->50<--------->    *
  21 C.    *                                                                *
  22 C.    *   LOSS=2  :                                                    *
  23 C.    *   STRA=0  <----------GLANDZ-------------------><--GLANDO-->    *
  24 C.    *   LOSS=1,3:                                                    *
  25 C.    *   STRA=0  <---------------------GLANDZ-------------------->    *
  26 C.    *                                                                *
  27 C.    *   STRA=1  <-----------PAI---------------------><--GLANDZ-->    *
  28 #if defined(CERNLIB_ASHO)
  29 C.    *   STRA=2  <---PAI----><---ASHO---><----PAI----><--GLANDZ-->    *
  30 C.    *                                                                *
  31 #endif
  32 C.    *                                                                *
  33 C.    *   DNMIN :  an estimation of the number of collisions           *
  34 C.    *            with energy close to the ionization energy          *
  35 C.    *            (see PHYS333)                                       *
  36 C.    *                                                                *
  37 C.    *   Input  : DEMEAN (mean energy loss)                           *
  38 C.    *   Output : DE   (energy loss in the current step)              *
  39 C.    *                                                                *
  40 C.    *    ==>Called by : GTELEC,GTMUON,GTHADR                         *
  41 C.    *                                                                *
  42 C.    ******************************************************************
  43 C.
  44 #include "geant321/gcbank.inc"
  45 #include "geant321/gcjloc.inc"
  46 #include "geant321/gconsp.inc"
  47 #include "geant321/gcmate.inc"
  48 #include "geant321/gccuts.inc"
  49 #include "geant321/gckine.inc"
  50 #include "geant321/gcmulo.inc"
  51 #include "geant321/gcphys.inc"
  52 #include "geant321/gctrak.inc"
  53 *
  54       PARAMETER (EULER=0.577215,GAM1=EULER-1)
  55       PARAMETER (P1=.60715,P2=.67794,P3=.52382E-1,P4=.94753,
  56      +           P5=.74442,P6=1.1934)
  57       PARAMETER (DGEV=0.153536 E-3, DNLIM=50)
  58 #if defined(CERNLIB_ASHO)
  59       PARAMETER (ASHMIN=1,ASHMAX=30)
  60 #endif
  61       DIMENSION RNDM(2)
  62       FLAND(X) = P1+P6*X+(P2+P3*X)*EXP(P4*X+P5)
  63 *
  64       IF(STEP.LE.0) THEN
  65          DE=DEMEAN
  66       ELSE
  67          DEDX = DEMEAN/STEP
  68          POTI=Q(JPROB+9)
  69          IF(ISTRA.EQ.0.AND.ILOSS.NE.2) THEN
  70             CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,Q(JPROB+10))
  71          ELSE
  72 *
  73 * *** mean ionization potential (GeV)
  74 *        POTI=16E-9*Z**0.9
  75 *
  76             GAMMA = GETOT/AMASS
  77             BETA = VECT(7)/GETOT
  78             BET2 = BETA**2
  79 *
  80 * ***    low energy transfer
  81             XI = DGEV*CHARGE**2*STEP*DENS*Z/(A*BET2)
  82 *
  83 * ***    regime
  84 * ***    ISTRA = 1 --> PAI + URBAN
  85 #if defined(CERNLIB_ASHO)
  86 * ***    ISTRA = 2 --> PAI + URBAN + ASHO
  87 #endif
  88             DNMIN = MIN(XI,DEMEAN)/POTI
  89 *
  90             IF (ISTRA.EQ.0) THEN
  91                IF(DNMIN.GE.DNLIM) THEN
  92 *
  93 *  Energy straggling using Gaussian, Landau & Vavilov theories.
  94 *
  95 *  STEP   =  current step-length (cm)
  96 *
  97 *  DELAND =  DE/DX - <DE/DX>     (GeV)
  98 *
  99 *  Author      : G.N. Patrick
 100 *
 101                   IF(STEP.LT.1.E-7)THEN
 102                      DELAND=0.
 103                   ELSE
 104 *
 105 *     Maximum energy transfer to atomic electron (GeV).
 106                      ETA = BETA*GAMMA
 107                      RATIO = EMASS/AMASS
 108 *
 109 *     Calculate Kappa significance ratio.
 110 *                 EMAX=(2*EMASS*ETA**2)/(1+2*RATIO*GAMMA+RATIO**2)
 111 *                 CAPPA = XI/EMAX
 112                      CAPPA = XI*(1+2*RATIO*GAMMA+RATIO**2)/(2*EMASS*
 113      +               ETA**2)
 114                      IF (CAPPA.GE.10.) THEN
 115 *
 116 *     +-----------------------------------+
 117 *     I Sample from Gaussian distribution I
 118 *     +-----------------------------------+
 119                         SIGMA = XI*SQRT((1.-0.5*BET2)/CAPPA)
 120                         CALL GRNDM(RNDM,2)
 121                         F1 = -2.*LOG(RNDM(1))
 122                         DELAND = SIGMA*SQRT(F1)*COS(TWOPI*RNDM(2))
 123                      ELSE
 124                         XMEAN = -BET2-LOG(CAPPA)+GAM1
 125                         IF (CAPPA.LT.0.01) THEN
 126                            XLAMX = FLAND(XMEAN)
 127 *     +---------------------------------------------------------------+
 128 *     I Sample lambda variable from Kolbig/Schorr Landau distribution I
 129 *     +---------------------------------------------------------------+
 130 *  10                   CALL GRNDM(RNDM,1)
 131 *                       IF( RNDM(1) .GT. 0.980 ) GO TO 10
 132 *                       XLAMB = GLANDR(RNDM(1))
 133 *     +---------------------------------------------------------------+
 134 *     I Sample lambda variable from James/Hancock Landau distribution I
 135 *     +---------------------------------------------------------------+
 136    10                      CALL GLANDG(XLAMB)
 137                            IF(XLAMB.GT.XLAMX) GO TO 10
 138                         ELSE
 139 *            +---------------------------------------------------------+
 140 *            I Sample lambda variable (Landau not Vavilov) from        I
 141 *            I Rotondi&Montagna&Kolbig Vavilov distribution            I
 142 *            +---------------------------------------------------------+
 143                            CALL GRNDM(RNDM,1)
 144                            XLAMB = GVAVIV(CAPPA,BET2,RNDM(1))
 145                         ENDIF
 146 *
 147 *     Calculate DE/DX - <DE/DX>
 148                         DELAND = XI*(XLAMB-XMEAN)
 149                      ENDIF
 150                   ENDIF
 151                   DE = DEMEAN + DELAND
 152                ELSE
 153                   CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,
 154      +            Q(JPROB+ 10))
 155                ENDIF
 156             ELSE IF (ISTRA.LE.2) THEN
 157                IF(DNMIN.GE.DNLIM) THEN
 158                   CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,
 159      +            Q(JPROB+ 10))
 160                ELSE
 161                   NMEC = NMEC+1
 162                   LMEC(NMEC)=109
 163 #if defined(CERNLIB_ASHO)
 164                   IF (DNMIN.GE.ASHMIN.AND.DNMIN.LT.ASHMAX .AND.ISTRA.EQ
 165      +            .2) THEN
 166                      CALL GASHO(VECT(7),AMASS,STEP,DE)
 167                   ELSE
 168                      DE = GSTREN(GAMMA,DCUTE,STEP)
 169                   ENDIF
 170 #endif
 171 #if !defined(CERNLIB_ASHO)
 172                   DE = GSTREN(GAMMA,DCUTE,STEP)
 173 #endif
 174 *
 175 * ***   Add brem losses to ionisation
 176                   IF(ITRTYP.EQ.2) THEN
 177                      JBASE = LQ(JMA-1)+2*NEK1+IEKBIN
 178                      DE = DE +(1.-GEKRAT)*Q(JBASE)+GEKRAT*Q(JBASE+1)
 179                   ELSEIF(ITRTYP.EQ.5) THEN
 180                      JBASE = LQ(JMA-2)+NEK1+IEKBIN
 181                      DE = DE +(1.-GEKRAT)*Q(JBASE)+GEKRAT*Q(JBASE+1)
 182                   ENDIF
 183                ENDIF
 184             ENDIF
 185          ENDIF
 186       ENDIF
 187 *
 188       END