1 *CMZ :          02/02/99  18.09.13  by  Federico Carminati
2 *-- Author :    Federico Carminati   02/02/99
3       SUBROUTINE GFLUCT(DEMEAN,DE)
4 C.
5 C.    ******************************************************************
6 C.    *                                                                *
7 C.    *   Subroutine to decide which method is used to simulate        *
8 C.    *   the straggling around the mean energy loss.                  *
9 C.    *                                                                *
10 C.    *                                                                *
11 C.    *   DNMIN:  <---------->1<-------->30<--------->50<--------->    *
12 C.    *                                                                *
13 C.    *   LOSS=2  :                                                    *
14 C.    *   STRA=0  <----------GLANDZ-------------------><--GLANDO-->    *
15 C.    *   LOSS=1,3:                                                    *
16 C.    *   STRA=0  <---------------------GLANDZ-------------------->    *
17 C.    *                                                                *
18 C.    *   STRA=1  <-----------PAI---------------------><--GLANDZ-->    *
19 C.    *                                                                *
20 C.    *   DNMIN :  an estimation of the number of collisions           *
21 C.    *            with energy close to the ionization energy          *
22 C.    *            (see PHYS333)                                       *
23 C.    *                                                                *
24 C.    *   Input  : DEMEAN (mean energy loss)                           *
25 C.    *   Output : DE   (energy loss in the current step)              *
26 C.    *                                                                *
27 C.    *    ==>Called by : GTELEC,GTMUON,GTHADR                         *
28 C.    *                                                                *
29 C.    ******************************************************************
30 C.
31 #include "geant321/pilot.h"
32 #undef CERNLIB_GEANT321_GCBANK_INC
34 #include "geant321/gcbank.inc"
35 #undef CERNLIB_GEANT321_GCJLOC_INC
36 #include "geant321/gcjloc.inc"
37 #undef CERNLIB_GEANT321_GCONSP_INC
38 #include "geant321/gconsp.inc"
39 #undef CERNLIB_GEANT321_GCMATE_INC
40 #include "geant321/gcmate.inc"
41 #undef CERNLIB_GEANT321_GCCUTS_INC
42 #include "geant321/gccuts.inc"
43 #undef CERNLIB_GEANT321_GCKINE_INC
44 #include "geant321/gckine.inc"
45 #undef CERNLIB_GEANT321_GCMULO_INC
46 #include "geant321/gcmulo.inc"
47 #undef CERNLIB_GEANT321_GCPHYS_INC
48 #include "geant321/gcphys.inc"
49 #undef CERNLIB_GEANT321_GCTRAK_INC
50 #include "geant321/gctrak.inc"
51 *KEND.
52       PARAMETER (EULER=0.577215,GAM1=EULER-1)
53       PARAMETER (P1=.60715,P2=.67794,P3=.52382E-1,P4=.94753,
54      +           P5=.74442,P6=1.1934)
55       PARAMETER (DGEV=0.153536 E-3, DNLIM=50)
56 * These parameters are needed by M.Kowalski's fluctuation algorithm
57       PARAMETER (FPOT=20.77E-9, EEND=10E-6, EEXPO=2.2)
58       PARAMETER (XEXPO=-EEXPO+1, YEXPO=1/XEXPO)
59 * These parameters are needed by M.Kowalski's fluctuation algorithm
60       DIMENSION RNDM(2)
61       DE2(DPOT,RAN)=(DPOT**XEXPO*(1-RAN)+EEND**XEXPO*RAN)**YEXPO
62       FLAND(X) = P1+P6*X+(P2+P3*X)*EXP(P4*X+P5)
63       IF(STEP.LE.0) THEN
64          DE=DEMEAN
65       ELSE
66          DEDX = DEMEAN/STEP
67          POTI=Q(JPROB+9)
68          IF(ISTRA.EQ.0.AND.(ILOSS.EQ.1.OR.ILOSS.EQ.3)) THEN
69             CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,Q(JPROB+10))
70          ELSEIF (ILOSS.EQ.5) THEN
71 * This is Marek Kowalski's fluctuation algorithm, it works only when
72 * the step size has been limited to one ionisation on average
73             CALL GRNDM(RNDM,1)
74             DE=DE2(FPOT,RNDM(1))
75 *
76          ELSE
77 * *** mean ionization potential (GeV)
78 *        POTI=16E-9*Z**0.9
79             GAMMA = GETOT/AMASS
80             BETA = VECT(7)/GETOT
81             BET2 = BETA**2
82 * ***    low energy transfer
83             XI = DGEV*CHARGE**2*STEP*DENS*Z/(A*BET2)
84 * ***    regime
85 * ***    ISTRA = 1 --> PAI + URBAN
86             DNMIN = MIN(XI,DEMEAN)/POTI
87             IF (ISTRA.EQ.0) THEN
88                IF(DNMIN.GE.DNLIM) THEN
89 *  Energy straggling using Gaussian, Landau & Vavilov theories.
90 *  STEP   =  current step-length (cm)
91 *  DELAND =  DE/DX - <DE/DX>     (GeV)
92 *  Author      : G.N. Patrick
93                   IF(STEP.LT.1.E-7)THEN
94                      DELAND=0.
95                   ELSE
96 *     Maximum energy transfer to atomic electron (GeV).
97                      ETA = BETA*GAMMA
98                      RATIO = EMASS/AMASS
99 *     Calculate Kappa significance ratio.
100 *                 EMAX=(2*EMASS*ETA**2)/(1+2*RATIO*GAMMA+RATIO**2)
101 *                 CAPPA = XI/EMAX
102                      CAPPA = XI*(1+2*RATIO*GAMMA+RATIO**2)/(2*EMASS*
103      +               ETA**2)
104                      IF (CAPPA.GE.10.) THEN
105 *     +-----------------------------------+
106 *     I Sample from Gaussian distribution I
107 *     +-----------------------------------+
108                         SIGMA = XI*SQRT((1.-0.5*BET2)/CAPPA)
109                         CALL GRNDM(RNDM,2)
110                         F1 = -2.*LOG(RNDM(1))
111                         DELAND = SIGMA*SQRT(F1)*COS(TWOPI*RNDM(2))
112                      ELSE
113                         XMEAN = -BET2-LOG(CAPPA)+GAM1
114                         IF (CAPPA.LT.0.01) THEN
115                            XLAMX = FLAND(XMEAN)
116 *     +---------------------------------------------------------------+
117 *     I Sample lambda variable from Kolbig/Schorr Landau distribution I
118 *     +---------------------------------------------------------------+
119 *  10                   CALL GRNDM(RNDM,1)
120 *                       IF( RNDM(1) .GT. 0.980 ) GO TO 10
121 *                       XLAMB = GLANDR(RNDM(1))
122 *     +---------------------------------------------------------------+
123 *     I Sample lambda variable from James/Hancock Landau distribution I
124 *     +---------------------------------------------------------------+
125    10                      CALL GLANDG(XLAMB)
126                            IF(XLAMB.GT.XLAMX) GO TO 10
127                         ELSE
128 *            +---------------------------------------------------------+
129 *            I Sample lambda variable (Landau not Vavilov) from        I
130 *            I Rotondi&Montagna&Kolbig Vavilov distribution            I
131 *            +---------------------------------------------------------+
132                            CALL GRNDM(RNDM,1)
133                            XLAMB = GVAVIV(CAPPA,BET2,RNDM(1))
134                         ENDIF
135 *     Calculate DE/DX - <DE/DX>
136                         DELAND = XI*(XLAMB-XMEAN)
137                      ENDIF
138                   ENDIF
139                   DE = DEMEAN + DELAND
140                ELSE
141                   CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,
142      +            Q(JPROB+ 10))
143                ENDIF
144             ELSE IF (ISTRA.LE.2) THEN
145                IF(DNMIN.GE.DNLIM) THEN
146                   CALL GLANDZ(Z,STEP,VECT(7),GETOT,DEDX,DE,POTI,
147      +            Q(JPROB+ 10))
148                ELSE
149                   NMEC = NMEC+1
150                   LMEC(NMEC)=109
151                   DE = GSTREN(GAMMA,DCUTE,STEP)
152 * ***   Add brem losses to ionisation
153                   IF(ITRTYP.EQ.2) THEN
154                      JBASE = LQ(JMA-1)+2*NEK1+IEKBIN
155                      DE = DE +(1.-GEKRAT)*Q(JBASE)+GEKRAT*Q(JBASE+1)
156                   ELSEIF(ITRTYP.EQ.5) THEN
157                      JBASE = LQ(JMA-2)+NEK1+IEKBIN
158                      DE = DE +(1.-GEKRAT)*Q(JBASE)+GEKRAT*Q(JBASE+1)
159                   ENDIF
160                ENDIF
161             ENDIF
162          ENDIF
163       ENDIF
164       END
165 *CMZ :          16/02/99  19.24.38  by  Federico Carminati
166 *CMZ :  2.03/01 28/08/98  09.33.11  by  Federico Carminati
167 *CMZ :  2.01/00 18/06/98  17.34.13  by  Federico Carminati
168 *CMZ :  2.00/05 28/05/98  19.04.13  by  Federico Carminati
169 *-- Author :
170       SUBROUTINE GTREVE
171 C.
172 C.    ******************************************************************
173 C.    *                                                                *
174 C.    *    SUBR. GTREVE                                                *
175 C.    *                                                                *
176 C.    *   Controls tracking of all particles belonging to the current  *
177 C.    *    event.                                                      *
178 C.    *                                                                *
179 C.    *   Called by : GUTREV, called by GTRIG                          *
180 C.    *   Authors   : R.Brun, F.Bruyant                                *
181 C.    *                                                                *
182 C.    ******************************************************************
183 C.
184 #include "geant321/pilot.h"
185 #undef CERNLIB_GEANT321_GCBANK_INC
187 #include "geant321/gcbank.inc"
188 #undef CERNLIB_GEANT321_GCFLAG_INC
189 #include "geant321/gcflag.inc"
190 #undef CERNLIB_GEANT321_GCKINE_INC
191 #include "geant321/gckine.inc"
192 #undef CERNLIB_GEANT321_GCNUM_INC
193 #include "geant321/gcnum.inc"
194 #undef CERNLIB_GEANT321_GCSTAK_INC
195 #include "geant321/gcstak.inc"
196 #undef CERNLIB_GEANT321_GCTMED_INC
197 #include "geant321/gctmed.inc"
198 #undef CERNLIB_GEANT321_GCTRAK_INC
199 #include "geant321/gctrak.inc"
200 #undef CERNLIB_GEANT321_GCUNIT_INC
201 #include "geant321/gcunit.inc"
202 #include "sckine.inc"
203 *KEND.
204 *
205       REAL UBUF(2)
206       EQUIVALENCE (UBUF(1),WS(1))
207       LOGICAL   BTEST
208       DIMENSION PMOM(3),VPOS(3)
209 C.
210 C.    ------------------------------------------------------------------
211       NTMSTO = 0
212       NSTMAX = 0
213       NALIVE = 0
214 *         Kick start the creation of the vertex
215       VPOS(1)=0
216       VPOS(2)=0
217       VPOS(3)=0
218       PMOM(1)=0
219       PMOM(2)=0
220       PMOM(3)=0
221       IPART=1
222       CALL GSVERT(VPOS,0,0,UBUF,0,NVTX)
223       CALL GSKINE(PMOM,IPART,NVTX,UBUF,0,NT)
224 *
225       MTRACK=-999
226  10   MTROLD=MTRACK
227       CALL RXGTRAK(MTRACK,IPART,PMOM,E,VPOS,TTOF)
228       IF(MTROLD.LT.0) THEN
229          MPRIMA=MTRACK
230       ENDIF
231       IF(MTRACK.LE.MPRIMA) THEN
232          IF(ISWIT(4).GT.0.AND.MTRACK.GT.0) THEN
233             IF(ISWIT(4).EQ.1.OR.MOD(MTRACK,ISWIT(4)).EQ.0) THEN
234                WRITE(CHMAIL,10200) MTRACK
235                CALL GMAIL(0,0)
236             ENDIF
237          ENDIF
238          IF(MTROLD.GT.0) THEN
239 C --- Output root hits tree only for each primary MTRACK
240             CALL RXOUTH
241          ENDIF
242       ENDIF
243       IF(MTRACK.LE.0) GOTO 999
244 * Set the vertex
245       JV=LQ(JVERTX-1)
246       Q(JV + 1) = VPOS(1)
247       Q(JV + 2) = VPOS(2)
248       Q(JV + 3) = VPOS(3)
249       Q(JV + 4) = TTOF
250       Q(JV + 5) = 0
251       Q(JV + 6) = 0
252 * Set the track
253       JK=LQ(JKINE-1)
254       Q(JK + 1) = PMOM(1)
255       Q(JK + 2) = PMOM(2)
256       Q(JK + 3) = PMOM(3)
257       Q(JK + 4) = E
258       Q(JK + 5) = IPART
259       Q(JK + 6) = 1
260 * Now transport
261 C      CALL GPVERT(0)
262 C      CALL GPKINE(0)
263 * Normal Gtreve code
264       NV = NVERTX
265       DO 40  IV = 1,NV
266 * ***   For each vertex in turn ..
267          JV = LQ(JVERTX-IV)
268          NT = Q(JV+7)
269          IF (NT.LE.0) GO TO 40
270          TOFG   = Q(JV+4)
271          SAFETY = 0.
272          IF (NJTMAX.GT.0) THEN
273             CALL GMEDIA (Q(JV+1), NUMED)
274             IF (NUMED.EQ.0) THEN
275                WRITE (CHMAIL, 10000) (Q(JV+I), I=1,3)
276                CALL GMAIL (0, 0)
277                GO TO 40
278             ENDIF
279             CALL GLSKLT
280          ENDIF
281 *  **   Loop over tracks attached to current vertex
282          DO 20  IT = 1,NT
283             JV   = LQ(JVERTX-IV)
284             ITRA = Q(JV+7+IT)
285             IF (BTEST(IQ(LQ(JKINE-ITRA)),0)) GO TO 20
286             CALL GFKINE (ITRA, VERT, PVERT, IPART, IVERT, UBUF, NWBUF)
287             IF (IVERT.NE.IV) THEN
288                WRITE (CHMAIL, 10100) IV, IVERT
289                CALL GMAIL (0, 0)
290                GO TO 999
291             ENDIF
292 *   *      Store current track parameters in stack JSTAK
293             CALL GSSTAK (2)
294    20    CONTINUE
295 *  **   Start tracking phase
296    30    IF (NALIVE.NE.0) THEN
297             NALIVE = NALIVE -1
298 *   *      Pick-up next track in stack JSTAK, if any
299 *   *         Initialize tracking parameters
300             CALL GLTRAC
301             IF (NUMED.EQ.0) GO TO 30
302 *   *       Resume tracking
303             CALL GUTRAK
304             IF (IEOTRI.NE.0) GO TO 999
305             GO TO 30
306          ENDIF
307 *
308    40 CONTINUE
309       GOTO 10
310 *
311 10000 FORMAT (' GTREVE : Vertex outside setup, XYZ=',3G12.4)
312 10100 FORMAT (' GTREVE : Abnormal track/vertex connection',2I8)
313 10200 FORMAT (' GTREVE : Transporting primary track No ',I10)
314 *                                                             END GTREVE
315   999 END