This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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