GEANT321/gphys/gfluct.F

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