]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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 | |
33 | #undef CERNLIB_GEANT321_GCLINK_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 | |
186 | #undef CERNLIB_GEANT321_GCLINK_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 |