This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gtrak / gtmuon.F
1 #include "geant321/pilot.h"
2 *CMZ :  3.21/02 29/03/94  15.41.24  by  S.Giani
3 *-- Author :
4       SUBROUTINE GTMUON
5 C.
6 C.    ******************************************************************
7 C.    *                                                                *
8 C.    *   Muon track. Computes step size and propagates particle       *
9 C.    *    through step.                                               *
10 C.    *                                                                *
11 C.    *   ==>Called by : GTRACK                                        *
12 C.    *      Authors     R.Brun, F.Bruyant, M.Maire ********           *
13 C.    *                                                                *
14 C.    ******************************************************************
15 C.
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
48 C.
49 C.    ------------------------------------------------------------------
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)
217 10000 FORMAT(' Boundary correction in GTMUON: ',
218      +       '    GEKIN      NUMED       STEP      SNEXT')
219 10100 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
468 c----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