Introduced M.Kowalski modifications for very short steps.
[u/mrichter/AliRoot.git] / 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