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