Introduced M.Kowalski modifications for very short steps.
[u/mrichter/AliRoot.git] / GEANT321 / gphys / gfluct.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
24b41ee6 5* Revision 1.1.1.1 1999/05/18 15:55:20 fca
6* AliRoot sources
7*
fe4da5cc 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)
16C.
17C. ******************************************************************
18C. * *
19C. * Subroutine to decide which method is used to simulate *
20C. * the straggling around the mean energy loss. *
21C. * *
22C. * *
23C. * DNMIN: <---------->1<-------->30<--------->50<---------> *
24C. * *
25C. * LOSS=2 : *
26C. * STRA=0 <----------GLANDZ-------------------><--GLANDO--> *
27C. * LOSS=1,3: *
28C. * STRA=0 <---------------------GLANDZ--------------------> *
29C. * *
30C. * STRA=1 <-----------PAI---------------------><--GLANDZ--> *
31#if defined(CERNLIB_ASHO)
32C. * STRA=2 <---PAI----><---ASHO---><----PAI----><--GLANDZ--> *
33C. * *
34#endif
35C. * *
36C. * DNMIN : an estimation of the number of collisions *
37C. * with energy close to the ionization energy *
38C. * (see PHYS333) *
39C. * *
40C. * Input : DEMEAN (mean energy loss) *
41C. * Output : DE (energy loss in the current step) *
42C. * *
43C. * ==>Called by : GTELEC,GTMUON,GTHADR *
44C. * *
45C. ******************************************************************
46C.
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
24b41ee6 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
fe4da5cc 68 DIMENSION RNDM(2)
24b41ee6 69 DE2(DPOT,RAN)=(DPOT**XEXPO*(1-RAN)+EEND**XEXPO*RAN)**YEXPO
fe4da5cc 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)
24b41ee6 77 IF(ISTRA.EQ.0.AND.(ILOSS.EQ.1.OR.ILOSS.EQ.3)) THEN
fe4da5cc 78 CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,Q(JPROB+10))
24b41ee6 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*
fe4da5cc 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