GEANT321/gphys/gprobi.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:21:33  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.23  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE GPROBI
  13 C.
  14 C.    ******************************************************************
  15 C.    *                                                                *
  16 C.    *  Initialise material constants used in the computation of      *
  17 C.    *  the probability for various interactions.                     *
  18 C.    *                                                                *
  19 C.    *    ==>Called by : GPHYSI                                       *
  20 C.    *       Authors    R.Brun, G.Patrick, L.Urban  *********         *
  21 C.    *                                                                *
  22 C.    ******************************************************************
  23 C.
  24 #include "geant321/gcbank.inc"
  25 #include "geant321/gconsp.inc"
  26 #include "geant321/gcmate.inc"
  27 #include "geant321/gcjloc.inc"
  28       DIMENSION EK(4),EL1(4),EL2(4)
  29       DATA EK / 0.66644E-8 , 0.22077E-9 ,-0.32552E-11, 0.18199E-13/
  30       DATA EL1/-0.29179E-9 , 0.87983E-10,-0.12589E-11, 0.69602E-14/
  31       DATA EL2/-0.68606E-9 , 0.10078E-9 ,-0.14496E-11, 0.78809E-14/
  32       DATA ALFA/7.29735E-3/
  33       DATA REL/0.2817938/
  34 C.
  35 C.    ------------------------------------------------------------------
  36 C.
  37       IF(Z.LT.1.) GOTO 999
  38       AEFF=A
  39       JPROB = LQ(JMA-4)
  40       IF(JMIXT.GT.0)THEN
  41          JMI1=LQ(JMIXT-1)
  42          AEFF=Q(JMI1+1)
  43       ENDIF
  44 C
  45 C             store constants for PAIR/BREMS routines
  46 C
  47       X      = (Z*ALFA)**2
  48       FC     = (( - 0.002 * X + 0.0083) * X - 0.0369) * X + 0.20206
  49       FC     = X * (FC + 1. / (1. + X))
  50       C1=Z**0.333333
  51       C2=LOG(C1)
  52       C3=LOG(183./C1)-FC
  53       C4=LOG(1440./(C1*C1))/C3
  54       Q(JPROB+1)=Z*(Z+C4)*C3/A
  55       Q(JPROB+2)=C1
  56       Q(JPROB+3)=C2
  57       Q(JPROB+4)=FC
  58 C
  59 C             constants for PHOTOEFFECT
  60 C
  61       Z2   = Z*Z
  62       EKZ  = Z2*(EK(1) +Z*(EK(2) +Z*(EK(3) +Z*EK(4))))
  63       EL1Z = Z2*(EL1(1)+Z*(EL1(2)+Z*(EL1(3)+Z*EL1(4))))
  64       EL2Z = Z2*(EL2(1)+Z*(EL2(2)+Z*(EL2(3)+Z*EL2(4))))
  65       Q(JPROB+5)=EKZ
  66       Q(JPROB+6)=EL1Z
  67       Q(JPROB+7)=EL2Z
  68 C
  69 C             Constants for Hadronic interactions
  70 C
  71       Q(JPROB+8)= 1000.*AEFF/(AVO*DENS)
  72 C
  73 C             Constants for electron/positron ionisation losses
  74 C             and S5 for one-photon annihilation
  75 C
  76       IF(JMIXT.LE.0)THEN
  77          POTI=16.E-9*Z**0.9
  78          S1=Z/A
  79          S5=Z**5/A*ALFA**4
  80       ELSE
  81          NLMAT=Q(JMA+11)
  82          NLM=IABS(NLMAT)
  83          S1=0.
  84          S2=0.
  85          S5=0.
  86          DO 10 J=1,NLM
  87             AJ=Q(JMIXT+J)
  88             ZJ=Q(JMIXT+NLM+J)
  89             WJ=Q(JMIXT+2*NLM+J)
  90             S1=S1+WJ*ZJ/AJ
  91             S2=S2+WJ*ZJ*LOG(ZJ)/AJ
  92             S5=S5+WJ*ZJ**5/AJ*ALFA**4
  93    10    CONTINUE
  94          POTI=16.E-9*EXP(0.9*S2/S1)
  95       ENDIF
  96       Q(JPROB+9) = POTI
  97       Q(JPROB+10) = LOG(POTI)
  98 C
  99       CON1=LOG(POTI/EMASS)
 100       CON2=DENS*S1
 101       CON3=1.+2.*LOG(POTI/(28.8E-9*SQRT(CON2)))
 102 C
 103 C             Condensed material ?
 104 C             (at present that means: DENS.GT.0.05 g/cm**3)
 105 C
 106       IF(DENS.GT.0.05)THEN
 107          IF(POTI.LT.1.E-7)THEN
 108             IF(CON3.LT.3.681)THEN
 109                CON4=0.2
 110             ELSE
 111                CON4=0.326*CON3-1.
 112             ENDIF
 113             CON5=2.
 114          ELSE
 115             IF(CON3.LT.5.215)THEN
 116                CON4=0.2
 117             ELSE
 118                CON4=0.326*CON3-1.5
 119             ENDIF
 120             CON5=3.
 121          ENDIF
 122       ELSE
 123 C
 124 C             Gas (T=0 C, P= 1 ATM)
 125 C             if T.NE. 0 C and/or P.NE. 1 ATM
 126 C             you have to modify the variable X
 127 C             X=>X+0.5*LOG((273+T C)/(273*P ATM))
 128 C             in the function GDRELE
 129 C             ------------------------
 130 C
 131          IF(CON3.LE.12.25)THEN
 132             IP=INT((CON3-10.)/0.5)+1
 133             IF(IP.LT.0) IP=0
 134             IF(IP.GT.4) IP=4
 135             CON4=1.6+0.1*FLOAT(IP)
 136             CON5=4.
 137          ELSE
 138             IF(CON3.LE.13.804)THEN
 139                CON4=2.
 140                CON5=5.
 141             ELSE
 142                CON4=0.326*CON3-2.5
 143                CON5=5.
 144             ENDIF
 145          ENDIF
 146       ENDIF
 147 C
 148       XA=CON3/4.606
 149       CON6=4.606*(XA-CON4)/(CON5-CON4)**3.
 150       Q(JPROB+11)=CON1
 151       Q(JPROB+12)=CON2
 152       Q(JPROB+13)=-CON3
 153       Q(JPROB+14)=CON4
 154       Q(JPROB+15)=CON5
 155       Q(JPROB+16)=CON6
 156 C
 157 C            constant for delta rays
 158 C            (the same constant is used in the Compton
 159 C              and Annihilation subroutines )
 160 C            and for one-photon annihilation
 161 C
 162       Q(JPROB+17)=AVO*TWOPI*REL*REL*DENS*S1
 163       Q(JPROB+18)=AVO*TWOPI*REL*REL*DENS*S5
 164 C
 165 C            Constants for Moliere scattering
 166 C
 167       IF(JMIXT.LE.0)THEN
 168          CALL GMOLI(A,Z,1.,1,DENS,Q(JPROB+21),Q(JPROB+25))
 169       ELSE
 170          CALL GMOLI(Q(JMIXT+1),Q(JMIXT+NLM+1),Q(JMIXT+2*NLM+1),
 171      +              NLM,DENS,Q(JPROB+21),Q(JPROB+25))
 172       ENDIF
 173 C
 174 C                Constants for muon bremsstrahlung
 175 C
 176       Q(JPROB+31)=LOG(189.*EMMU/(EMASS*C1))
 177       IF(Z.GT.10)Q(JPROB+31)=Q(JPROB+31)+LOG(0.666666/C1)
 178       SE         =SQRT(2.71828)
 179       Q(JPROB+32)=189.*SE*EMMU*EMMU/(2.*EMASS*C1)
 180       Q(JPROB+33)=0.75*SE*EMMU*C1
 181 C
 182   999 END