GEANT321/gphys/glandz.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/04 22/02/95  10.38.36  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GLANDZ(Z,STEP,P,E,DEDX,DE,POTI,POTIL)
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *  Energy straggling using a Monte Carlo model.                  *
  17 C.    *  It can be used with or without delta ray generation.          *
  18 C.    *                                                                *
  19 C.    *   It is a NEW VERSION of the model , which reproduces          *
  20 C.    *   the experimental data rather well.                           *
  21 C.    *                                                                *
  22 C.    *  Input : STEP  = current step-length (cm)                      *
  23 C.    *                                                                *
  24 C.    *  Output: DE    = actual energy loss (Gev)                      *
  25 C.    *                 ( NOT the fluctuation DE/DX-<DE/DX> !)         *
  26 C.    *                                                                *
  27 C.    *     ==> Called by : GTELEC,GTHADR,GTMUON                       *
  28 C.    *                                                                *
  29 C.    *  Author      : L.Urban                                         *
  30 C.    *  Date        : 28.04.1988       Last update :  1.02.90         *
  31 C.    *                                                                *
  32 C.    ******************************************************************
  33 C.
  34 #include "geant321/gconsp.inc"
  35 #include "geant321/gcflag.inc"
  36 #include "geant321/gckine.inc"
  37 #include "geant321/gccuts.inc"
  38       PARAMETER (MAXRND=100)
  39       DIMENSION RNDM(MAXRND),APOIS(3),NPOIS(3),FPOIS(3)
  40 #if !defined(CERNLIB_SINGLE)
  41       DOUBLE PRECISION E1,E2,P2,B2,BG2,TM,DEC,W,WW,WWW,RR
  42       DOUBLE PRECISION X,AK,ALFA,EA,SA,AKL,REALK,FPOIS
  43       DOUBLE PRECISION ONE,HALF,ZERO
  44       PARAMETER (HMXINT=2**30)
  45 #endif
  46 #if defined(CERNLIB_SINGLE)
  47       PARAMETER (HMXINT=2**45)
  48 #endif
  49       PARAMETER (ONE=1,HALF=ONE/2,ZERO=0)
  50       PARAMETER (RCD=0.40, RCD1=1.-RCD, PROBLM=0.01)
  51       PARAMETER( C1=4, C2=16)
  52 *
  53 *     ------------------------------------------------------------------
  54 *
  55 *     POTI=1.6E-8*Z**0.9
  56 *     POTIL=LOG(POTI)
  57 *
  58       IF(Z.GT.2.) THEN
  59          F2=2./Z
  60       ELSE
  61          F2=0.
  62       ENDIF
  63       F1=1.-F2
  64 *
  65       E2 = 1.E-8*Z*Z
  66       E2L= LOG(E2)
  67       E1L= (POTIL-F2*E2L)/F1
  68       E1 = EXP(E1L)
  69 *
  70 *
  71       P2=P*P
  72       B2=P2/(E*E)
  73       BG2=P2/(AMASS*AMASS)
  74       IF(ITRTYP.EQ.2) THEN
  75          TM=P2/(E+EMASS)
  76          IF(CHARGE.LT.0.) TM=TM/2.
  77          TM=TM-POTI
  78          IF(TM.GT.DCUTE) TM=DCUTE
  79       ELSE
  80          TM=EMASS*P2/(0.5*AMASS*AMASS+EMASS*E)
  81          TM=TM-POTI
  82          IF(TM.GT.DCUTM) TM=DCUTM
  83       ENDIF
  84 *
  85 *
  86 * *** Protection against negative TM    ---------------------
  87 *     TM can be negative only for heavy particles with  a very
  88 *     low kinetic energy (e.g. for proton with T  100-300 kev)
  89       TM=MAX(TM,ZERO)
  90 *
  91       W  = TM+POTI
  92       WW = W/POTI
  93       WWW= 2.*EMASS*BG2
  94       WL = LOG(WWW)
  95       CSB=STEP*RCD1*DEDX/(WL-POTIL-B2)
  96       APOIS(1)=CSB*F1*(WL-E1L-B2)/E1
  97       APOIS(2)=CSB*F2*(WL-E2L-B2)/E2
  98 *
  99       IF(TM.GT.0.) THEN
 100          APOIS(3)=RCD*DEDX*STEP*TM/(POTI*W*LOG(WW))
 101       ELSE
 102          APOIS(1)=APOIS(1)/RCD1
 103          APOIS(2)=APOIS(2)/RCD1
 104          APOIS(3)=0.
 105       ENDIF
 106 *
 107 *    calculate the probability of the zero energy loss
 108 *
 109       APSUM=APOIS(1)+APOIS(2)+APOIS(3)
 110       IF(APSUM.LT.0.) THEN
 111          PROB=1.
 112       ELSEIF(APSUM.LT.50.) THEN
 113          PROB=EXP(-APSUM)
 114       ELSE
 115          PROB=0.
 116       ENDIF
 117 *
 118 *
 119 *      do it differently if prob > problm  <====================
 120       IF(PROB.GT.PROBLM) THEN
 121          E0=1.E-8
 122          EMEAN=DEDX*STEP
 123          IF(TM.LE.0.) THEN
 124 *      excitation only ....
 125             APOIS(1)=EMEAN/E0
 126 *
 127             CALL GPOISS(APOIS,NPOIS,1)
 128             FPOIS(1)=NPOIS(1)
 129             DE=FPOIS(1)*E0
 130 *
 131          ELSE
 132 *         ionization only ....
 133             EM=TM+E0
 134             APOIS(1)=EMEAN*(EM-E0)/(EM*E0*LOG(EM/E0))
 135             CALL GPOISS(APOIS,NPOIS,1)
 136             NN=NPOIS(1)
 137             DE=0.
 138 *
 139             IF(NN.GT.0) THEN
 140                RCORR=1.
 141                IF(NN.GT.MAXRND) THEN
 142                   RCORR=FLOAT(NN)/MAXRND
 143                   NN=MAXRND
 144 *
 145                ENDIF
 146                W=(EM-E0)/EM
 147                CALL GRNDM(RNDM,NN)
 148                DO 10 I=1,NN
 149                   DE=DE+E0/(1.-W*RNDM(I))
 150    10          CONTINUE
 151                DE=RCORR*DE
 152 *
 153             ENDIF
 154          ENDIF
 155          GOTO 999
 156       ENDIF
 157 *
 158       IF(MAX(APOIS(1),APOIS(2),APOIS(3)).LT.HMXINT) THEN
 159          CALL GPOISS(APOIS,NPOIS,3)
 160          FPOIS(1)=NPOIS(1)
 161          FPOIS(2)=NPOIS(2)
 162          FPOIS(3)=NPOIS(3)
 163       ELSE
 164          DO 20 JPOIS=1, 3
 165             IF(APOIS(JPOIS).LT.HMXINT) THEN
 166                CALL GPOISS(APOIS(JPOIS),NPOIS(JPOIS),1)
 167                FPOIS(JPOIS)=NPOIS(JPOIS)
 168             ELSE
 169                CALL GRNDM(RNDM,2)
 170                FPOIS(JPOIS)=ABS(SQRT(-2.*LOG(RNDM(1)*ONE)
 171      +         *APOIS(JPOIS))*SIN(TWOPI*RNDM(2)*ONE)+APOIS(JPOIS))
 172             ENDIF
 173    20    CONTINUE
 174       ENDIF
 175
 176 *
 177 *          Now we have all three numbers in REAL/DOUBLE
 178 *          variables. REALK is actually an INTEGER that now may
 179 *          exceed the machine representation limit for integers.
 180 *
 181       DE=FPOIS(1)*E1+FPOIS(2)*E2
 182 *
 183 *     smear to avoid peaks in the energy loss (note: E1<<E2)
 184 *
 185       IF(DE.GT.0.) THEN
 186          CALL GRNDM(RNDM,1)
 187          DE=DE+E1*(1.-2.*RNDM(1))
 188       ENDIF
 189 *
 190       ALFA=1.
 191       REALK=0
 192       DEC=0.
 193 *
 194       ANC=FPOIS(3)
 195       IF(ANC.GE.C2) THEN
 196          R=ANC/(C2+ANC)
 197          AN=ANC*R
 198          SN=C1*R
 199          CALL GRNDM(RNDM,2)
 200          RR=SQRT(-2.*LOG(RNDM(1)))
 201          PHI=TWOPI*RNDM(2)
 202          X=RR*COS(PHI)
 203          AK=AN+SN*X
 204          ALFA=WW*(C2+ANC)/(C2*WW+ANC)
 205          EA=AK*POTI*ALFA*LOG(ALFA)/(ALFA-1.)
 206          SA=SQRT(ABS(AK*ALFA*POTI*POTI-EA*EA/AK))
 207          AKL=(EA-C1*SA)/POTI
 208          IF(AK.LE.AKL) THEN
 209             X=RR*SIN(PHI)
 210             DEC=EA+SA*X
 211             REALK=AK+HALF-MOD(AK+HALF,ONE)
 212          ELSE
 213             ALFA=1.
 214          ENDIF
 215       ENDIF
 216       NN=NINT(FPOIS(3)-REALK)
 217       IF(NN.GT.MAXRND) THEN
 218          W=1.-ALFA/WW
 219          WW=POTI*ALFA
 220 *
 221 *     Here we take a gaussian distribution to avoid loosing
 222 *     too much time in computing
 223 *
 224          AVERAG=-LOG(1.-W)/W
 225          SIGMA =SQRT(NN*(1./(1.-W)-AVERAG**2))
 226          CALL GRNDM(RNDM,2)
 227          DEC   = DEC+WW*(NN*AVERAG+SIGMA*SQRT(-2.*LOG(RNDM(1)))*
 228      +           SIN(TWOPI*RNDM(2)))
 229          DE=DE+DEC
 230       ELSEIF(NN.GT.0) THEN
 231          W=1.-ALFA/WW
 232          WW=POTI*ALFA
 233          CALL GRNDM(RNDM,NN)
 234          DO 30 I=1,NN
 235             DEC=DEC+WW/(1.-W*RNDM(I))
 236    30    CONTINUE
 237          DE=DE+DEC
 238       ENDIF
 239 *
 240   999 END