]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/gtrak/gtmuon.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gtrak / gtmuon.F
CommitLineData
fe4da5cc 1#include "geant321/pilot.h"
2*CMZ : 3.21/02 29/03/94 15.41.24 by S.Giani
3*-- Author :
4 SUBROUTINE GTMUON
5C.
6C. ******************************************************************
7C. * *
8C. * Muon track. Computes step size and propagates particle *
9C. * through step. *
10C. * *
11C. * ==>Called by : GTRACK *
12C. * Authors R.Brun, F.Bruyant, M.Maire ******** *
13C. * *
14C. ******************************************************************
15C.
16#include "geant321/gcbank.inc"
17#include "geant321/gccuts.inc"
18#include "geant321/gcjloc.inc"
19#include "geant321/gckine.inc"
20#include "geant321/gcmate.inc"
21#include "geant321/gcmulo.inc"
22#include "geant321/gconsp.inc"
23#include "geant321/gcphys.inc"
24#include "geant321/gcstak.inc"
25#include "geant321/gctmed.inc"
26#include "geant321/gctrak.inc"
27#include "geant321/gcunit.inc"
28#include "geant321/gcking.inc"
29#if defined(CERNLIB_USRJMP)
30#include "geant321/gcjump.inc"
31#endif
32
33#if !defined(CERNLIB_OLD)
34#include "geant321/gcvolu.inc"
35#include "geant321/gcvdma.inc"
36#endif
37#if !defined(CERNLIB_SINGLE)
38 PARAMETER (EPSMAC=1.E-6)
39 DOUBLE PRECISION DEMEAN,STOPRG,STOPMX,STOPC
40 DOUBLE PRECISION ONE,XCOEF1,XCOEF2,XCOEF3
41#endif
42#if defined(CERNLIB_SINGLE)
43 PARAMETER (EPSMAC=1.E-11)
44#endif
45 PARAMETER (ONE=1)
46 REAL VNEXT(6)
47 SAVE IKCUT,STOPC
48C.
49C. ------------------------------------------------------------------
50*
51* *** Particle below energy threshold ? short circuit
52*
53 IF (GEKIN.LE.CUTMUO) GO TO 100
54*
55* *** Update local pointers if medium has changed
56 IF (IUPD.EQ.0) THEN
57 IUPD = 1
58 JLOSS = LQ(JMA-2)
59 JBREM = LQ(JMA-9)
60 JPAIR = LQ(JMA-10)
61 JDRAY = LQ(JMA-11)
62 JMUNU = LQ(JMA-14)
63 JRANG = LQ(JMA-16)
64 JCOEF = LQ(JMA-18)
65 JMULOF= LQ(JTM-2)
66 IF(IMCKOV.EQ.1) THEN
67 JTCKOV = LQ(JTM-3)
68 JABSCO = LQ(JTCKOV-1)
69 JEFFIC = LQ(JTCKOV-2)
70 JINDEX = LQ(JTCKOV-3)
71 JCURIN = LQ(JTCKOV-4)
72 NPCKOV = Q(JTCKOV+1)
73 ENDIF
74 OMCMOL= Q(JPROB+21)
75 CHCMOL= Q(JPROB+25)
76 IKCUT = Q(JMULOF+NEK1+1)
77 STOPC = Q(JMULOF+NEK1+2)
78 IF(ISTRA.GT.0) THEN
79 JTSTRA = LQ(JMA-19)
80 JTSTCO = LQ(JTSTRA-1)
81 JTSTEN = LQ(JTSTRA-2)
82#if defined(CERNLIB_ASHO)
83 IF(ISTRA.EQ.2) THEN
84 JTASHO = LQ(JMA-20)
85 ENDIF
86#endif
87 ENDIF
88 ENDIF
89*
90* *** Compute current step size
91*
92 STEP = STEMAX
93 IPROC = 103
94 GEKRT1 = 1. -GEKRAT
95 IEK1 = IEKBIN+NEK1
96 IEK2 = IEKBIN+2*NEK1
97*
98* ** Step limitation due to bremsstrahlung ?
99*
100 IF (IBREM.GT.0) THEN
101 STEPBR = GEKRT1*Q(JBREM+IEK2) +GEKRAT*Q(JBREM+IEK2+1)
102 SBREM = STEPBR*ZINTBR
103 IF (SBREM.LT.STEP) THEN
104 STEP = SBREM
105 IPROC = 9
106 ENDIF
107 ENDIF
108*
109* ** Step limitation due to pair production ?
110*
111 IF (IPAIR.GT.0) THEN
112 STEPPA = GEKRT1*Q(JPAIR+IEK1) +GEKRAT*Q(JPAIR+IEK1+1)
113 SPAIR = STEPPA*ZINTPA
114 IF (SPAIR.LT.STEP) THEN
115 STEP = SPAIR
116 IPROC = 6
117 ENDIF
118 ENDIF
119*
120* ** Step limitation due to decay ?
121*
122 IF (IDCAY.NE.0) THEN
123 SDCAY = SUMLIF*VECT(7)/AMASS
124 IF (SDCAY.LT.STEP) THEN
125 STEP = SDCAY
126 IPROC = 5
127 ENDIF
128 ENDIF
129*
130* ** Step limitation due to delta-ray ?
131*
132 IF (IDRAY.GT.0) THEN
133 STEPDR = GEKRT1*Q(JDRAY+IEK2) +GEKRAT*Q(JDRAY+IEK2+1)
134 SDRAY = STEPDR*ZINTDR
135 IF (SDRAY.LT.STEP) THEN
136 STEP = SDRAY
137 IPROC = 10
138 ENDIF
139 ENDIF
140*
141* ** Step limitation due to nuclear interaction ?
142*
143 IF (IMUNU.GT.0) THEN
144 IF(GEKIN.GE.5.)THEN
145 STEPMU = GEKRT1*Q(JMUNU+IEKBIN) +GEKRAT*Q(JMUNU+IEKBIN+1)
146 SMUNU = STEPMU*ZINTMU
147 IF (SMUNU.LT.STEP) THEN
148 STEP = SMUNU
149 IPROC = 21
150 ENDIF
151 ELSE
152 STEPMU = BIG
153 ENDIF
154 ENDIF
155*
156 IF (STEP.LE.0.) THEN
157 STEP = 0.
158 GO TO 90
159 ENDIF
160*
161* ** Step limitation due to energy-loss,multiple scattering
162* or magnetic field ?
163*
164 IF (JMULOF.NE.0) THEN
165 SMULOF = GEKRT1*Q(JMULOF+IEKBIN) +GEKRAT*Q(JMULOF+IEKBIN+1)
166 IF (SMULOF.LT.STEP) THEN
167 STEP = SMULOF
168 IPROC = 0
169 ENDIF
170 ENDIF
171*
172* ** Step limitation due to geometry ?
173*
174 IF (STEP.GE.0.95*SAFETY) THEN
175 CALL GTNEXT
176 IF (IGNEXT.NE.0) THEN
177 STEP = SNEXT + PREC
178 IPROC = 0
179 ENDIF
180*
181* Update SAFETY in stack companions, if any
182 IF (IQ(JSTAK+3).NE.0) THEN
183 DO 10 IST = IQ(JSTAK+3),IQ(JSTAK+1)
184 JST = JSTAK + 3 + (IST-1)*NWSTAK
185 Q(JST+11) = SAFETY
186 10 CONTINUE
187 IQ(JSTAK+3) = 0
188 ENDIF
189 ELSE
190 IQ(JSTAK+3) = 0
191 ENDIF
192*
193* *** Linear transport when no field or very short step
194*
195 IF (IFIELD.EQ.0.OR.STEP.LE.PREC) THEN
196*
197 IF (IGNEXT.NE.0) THEN
198 DO 20 I = 1,3
199 VECTMP = VECT(I) +STEP*VECT(I+3)
200 IF(VECTMP.EQ.VECT(I)) THEN
201*
202* *** Correct for machine precision
203*
204 IF(VECT(I+3).NE.0.) THEN
205 VECTMP =
206 + VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC
207 IF(NMEC.GT.0) THEN
208 IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
209 ENDIF
210 NMEC=NMEC+1
211 LMEC(NMEC)=104
212#if defined(CERNLIB_DEBUG)
213 WRITE(CHMAIL, 10000)
214 CALL GMAIL(0,0)
215 WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
216 CALL GMAIL(0,0)
21710000 FORMAT(' Boundary correction in GTMUON: ',
218 + ' GEKIN NUMED STEP SNEXT')
21910100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
220#endif
221 ENDIF
222 ENDIF
223 VECT(I) = VECTMP
224 20 CONTINUE
225 INWVOL = 2
226 NMEC = NMEC +1
227 LMEC(NMEC) = 1
228 ELSE
229 DO 30 I = 1,3
230 VECT(I) = VECT(I) +STEP*VECT(I+3)
231 30 CONTINUE
232 ENDIF
233 ELSE
234*
235* *** otherwise, swim particle in magnetic field
236*
237#if !defined(CERNLIB_OLD)
238 if(mycoun.gt.1.and.nfmany.gt.0.and.step.ge.safety)then
239 nlevel=manyle(nfmany)
240 do 99 i=1,nlevel
241 names(i)=manyna(nfmany,i)
242 number(i)=manynu(nfmany,i)
243 99 continue
244 call glvolu(nlevel,names,number,ier)
245 if(ier.ne.0)print *,'Fatal error in GLVOLU'
246 ingoto=0
247 endif
248#endif
249 NMEC = NMEC +1
250 LMEC(NMEC) = 4
251*
252#if !defined(CERNLIB_USRJMP)
253 40 CALL GUSWIM (CHARGE, STEP, VECT, VOUT)
254#endif
255#if defined(CERNLIB_USRJMP)
256 40 CALL JUMPT4(JUSWIM, CHARGE, STEP, VECT, VOUT)
257#endif
258*
259* ** When near to boundary, take proper action (cut-step,crossing...)
260*
261 IF(STEP.GE.SAFETY)THEN
262 INEAR = 0
263 IF (IGNEXT.NE.0) THEN
264 DO 50 I = 1,3
265 VNEXT(I+3) = VECT(I+3)
266 VNEXT(I) = VECT(I) +SNEXT*VECT(I+3)
267 50 CONTINUE
268 DO 60 I = 1,3
269 IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 70
270 60 CONTINUE
271 INEAR = 1
272 ENDIF
273*
274 70 CALL GINVOL (VOUT, ISAME)
275 IF (ISAME.EQ.0)THEN
276 IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN
277 INWVOL = 2
278 NMEC = NMEC +1
279 LMEC(NMEC) = 1
280 ELSE
281* Cut step
282 STEP = 0.5*STEP
283 IF (LMEC(NMEC).NE.24) THEN
284 NMEC = NMEC +1
285 LMEC(NMEC) = 24
286 ENDIF
287 GO TO 40
288 ENDIF
289 ENDIF
290 ENDIF
291*
292 DO 80 I = 1,6
293 VECT(I) = VOUT(I)
294 80 CONTINUE
295*
296 ENDIF
297*
298* *** Correct the step due to multiple scattering
299 IF (IMULL.NE.0) THEN
300 STMULS = STEP
301 CORR=0.0001*(STEP/RADL)*(GETOT/(VECT(7)*VECT(7)))**2
302 IF (CORR.GT.0.25) CORR = 0.25
303 STEP = (1.+CORR)*STEP
304 ENDIF
305*
306 SLENG = SLENG + STEP
307*
308* *** Generate Cherenkov photons if required
309*
310 IF(IMCKOV.EQ.1) THEN
311 CALL GGCKOV
312 IF(NGPHOT.NE.0) THEN
313 NMEC=NMEC+1
314 LMEC(NMEC)=105
315 ENDIF
316 ENDIF
317*
318* *** apply energy loss : find the kinetic energy corresponding
319* to the new stopping range = stopmx - step
320*
321 IF (ILOSL.NE.0) THEN
322 NMEC = NMEC +1
323 LMEC(NMEC) = 3
324 IF(GEKRAT.LT.0.7) THEN
325 I1 = MAX(IEKBIN-1,1)
326 ELSE
327 I1 = MIN(IEKBIN,NEKBIN-1)
328 ENDIF
329 I1 = 3*(I1-1)+1
330 XCOEF1 = Q(JCOEF+I1)
331 XCOEF2 = Q(JCOEF+I1+1)
332 XCOEF3 = Q(JCOEF+I1+2)
333 IF(XCOEF1.NE.0.) THEN
334 STOPMX = -XCOEF2+SIGN(ONE,XCOEF1)*SQRT(XCOEF2**2 - (XCOEF3-
335 + GEKIN/XCOEF1))
336 ELSE
337 STOPMX = - (XCOEF3-GEKIN)/XCOEF2
338 ENDIF
339 STOPRG = STOPMX - STEP
340 IF (STOPRG.LT.STOPC) THEN
341 STEP = STOPMX - STOPC
342 GO TO 100
343 ENDIF
344*
345 IF(XCOEF1.NE.0.) THEN
346 DEMEAN=GEKIN-XCOEF1*(XCOEF3+STOPRG*(2.*XCOEF2+STOPRG))
347 ELSE
348 DEMEAN=GEKIN-XCOEF2*STOPRG-XCOEF3
349 ENDIF
350 IF(DEMEAN.LE.5.*GEKIN*EPSMAC) THEN
351 DEMEAN=(GEKRT1*Q(JLOSS+IEKBIN)+GEKRAT*Q(JLOSS+IEKBIN+1))
352 + *STEP
353 ENDIF
354 IF (ILOSS.EQ.4.OR.IEKBIN.LE.IKCUT+1) THEN
355 DESTEP = DEMEAN
356 ELSE
357 DEMS = DEMEAN
358 CALL GFLUCT(DEMS,DESTEP)
359 ENDIF
360 IF (DESTEP.LT.0.) DESTEP = 0.
361 GEKINT = GEKIN -DESTEP
362 IF (GEKINT.LE.(1.01*CUTMUO)) GO TO 100
363 DESTEL = DESTEP
364 GEKIN = GEKINT
365 GETOT = GEKIN +AMASS
366 VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
367 CALL GEKBIN
368 ENDIF
369*
370* *** Apply multiple scattering.
371*
372 IF (IMULL.NE.0) THEN
373 NMEC = NMEC +1
374 LMEC(NMEC) = 2
375 CALL GMULTS
376 ENDIF
377*
378* *** Update time of flight
379*
380 SUMLIF = SUMLIF -STEP*AMASS/VECT(7)
381 TOFG = TOFG +STEP*GETOT/(VECT(7)*CLIGHT)
382 IF (TOFG.GE.TOFMAX) THEN
383 ISTOP = 4
384 NMEC = NMEC +1
385 LMEC(NMEC) = 22
386 GO TO 999
387 ENDIF
388*
389* *** Update interaction probabilities
390*
391 IF (IBREM.GT.0) ZINTBR = ZINTBR -STEP/STEPBR
392 IF (IPAIR.GT.0) ZINTPA = ZINTPA -STEP/STEPPA
393 IF (IDRAY.GT.0) ZINTDR = ZINTDR -STEP/STEPDR
394 IF (IMUNU.GT.0) ZINTMU = ZINTMU -STEP/STEPMU
395*
396* *** otherwise, apply the selected process if any
397*
398 90 IF (IPROC.EQ.0) GO TO 999
399 NMEC = NMEC +1
400 LMEC(NMEC) = IPROC
401*
402* ** Bremsstrahlung ?
403*
404 IF (IPROC.EQ.9) THEN
405 CALL GBREMM
406*
407* ** Pair production ?
408*
409 ELSE IF (IPROC.EQ.6) THEN
410 CALL GPAIRM
411*
412* ** Decay ?
413*
414 ELSE IF (IPROC.EQ.5) THEN
415 ISTOP = 1
416 CALL GDECAY
417*
418* ** Delta-ray ?
419*
420 ELSE IF (IPROC.EQ.10) THEN
421 CALL GDRAY
422*
423* ** Nuclear interaction ?
424*
425 ELSE IF (IPROC.EQ.21) THEN
426 CALL GMUNU
427 ENDIF
428 GO TO 999
429*
430* *** Special treatment for overstopped tracks
431*
432 100 DESTEP = GEKIN
433 DESTEL = DESTEP
434 GEKIN = 0.
435 GETOT = AMASS
436 VECT(7)= 0.
437 INWVOL = 0
438 NMEC = NMEC +1
439 LMEC(NMEC) = 30
440 IF (IDCAY.EQ.0) THEN
441 ISTOP = 2
442 ELSE
443 NMEC = NMEC +1
444 TOFG = TOFG +SUMLIF/CLIGHT
445 SUMLIF = 0.
446 IF (TOFG.GE.TOFMAX) THEN
447 ISTOP = 4
448 LMEC(NMEC) = 22
449 GO TO 999
450 ENDIF
451 LMEC(NMEC) = 5
452 ISTOP = 1
453 CALL GDECAY
454 ENDIF
455 999 IF(NGPHOT.GT.0) THEN
456 IF(ITCKOV.EQ.2.AND.ISTOP.EQ.0) THEN
457*
458* The muon has produced Cerenkov photons and it is still alive
459* we put it in the stack and we let the photons to be tracked
460 NGKINE = NGKINE+1
461 GKIN(1,NGKINE) = VECT(4)*VECT(7)
462 GKIN(2,NGKINE) = VECT(5)*VECT(7)
463 GKIN(3,NGKINE) = VECT(6)*VECT(7)
464 GKIN(4,NGKINE) = GETOT
465 GKIN(5,NGKINE) = IPART
466 TOFD(NGKINE) = 0.
467 ISTOP = 1
468c----put position as well
469 GPOS(1,NGKINE)=VECT(1)
470 GPOS(2,NGKINE)=VECT(2)
471 GPOS(3,NGKINE)=VECT(3)
472 ENDIF
473 ENDIF
474 END