This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gtrak / gtelec.F
1 #include "geant321/pilot.h"
2 *CMZ :  3.21/04 22/02/95  16.08.31  by  S.Giani
3 *-- Author :
4       SUBROUTINE GTELEC
5 C.
6 C.    ******************************************************************
7 C.    *                                                                *
8 C.    *   Electron type track. Computes step size and propagates       *
9 C.    *    particle through step.                                      *
10 C.    *                                                                *
11 C.    *   ==> Called by : GTRACK                                       *
12 C.    *       Authors    R.Brun, F.Bruyant, M.Maire L.Urban ********   *
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,ZERO
41 #endif
42 #if defined(CERNLIB_SINGLE)
43       PARAMETER (EPSMAC=1.E-11)
44 #endif
45       PARAMETER (ONE=1,ZERO=0.)
46       REAL VNEXT(6)
47       SAVE IKCUT,STOPC
48 *>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
49       DIMENSION RNDM(3)
50       PARAMETER ( TLIM = 0.0002)
51       IABAN = NINT(DPHYS1)
52 *<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
53 C.
54 C.    ------------------------------------------------------------------
55 *
56 * *** Particle below energy threshold ? short circuit
57 *
58       IF (GEKIN.LE.CUTELE) GO TO 100
59 *
60 * *** Update local pointers if medium or particle code has changed
61       IF (IUPD.EQ.0) THEN
62          IUPD  = 1
63          JMULOF= LQ(JTM-1)
64          IF (CHARGE.LT.0.) THEN
65             JBREM = LQ(JMA-9)
66             JLOSS = LQ(JMA-1)
67             JDRAY = LQ(JMA-11)
68             JRANG = LQ(JMA-15)
69             JCOEF = LQ(JMA-17)
70          ELSE
71             JBREM = LQ(JMA-9)  +NEK1
72             JLOSS = LQ(JMA-1)  +NEK1
73             JDRAY = LQ(JMA-11) +NEK1
74             JRANG = LQ(JMA-15) +NEK1
75             JCOEF = LQ(JMA-17) +3*NEK1
76             JANNI = LQ(JMA-7)
77          ENDIF
78          IF(IMCKOV.EQ.1) THEN
79             JTCKOV = LQ(JTM-3)
80             JABSCO = LQ(JTCKOV-1)
81             JEFFIC = LQ(JTCKOV-2)
82             JINDEX = LQ(JTCKOV-3)
83             JCURIN = LQ(JTCKOV-4)
84             NPCKOV = Q(JTCKOV+1)
85          ENDIF
86          OMCMOL = Q(JPROB+21)
87          CHCMOL = Q(JPROB+25)
88          IKCUT  = Q(JMULOF+NEK1+1)
89          STOPC  = Q(JMULOF+NEK1+2)
90          IF(ISTRA.GT.0) THEN
91             JTSTRA = LQ(JMA-19)
92             JTSTCO = LQ(JTSTRA-1)
93             JTSTEN = LQ(JTSTRA-2)
94 #if defined(CERNLIB_ASHO)
95             IF(ISTRA.EQ.2) THEN
96                JTASHO = LQ(JMA-20)
97             ENDIF
98 #endif
99          ENDIF
100       ENDIF
101 *
102 * *** Compute current step size
103 *
104       STEP   = STEMAX
105       IPROC  = 103
106       GEKRT1 = 1. -GEKRAT
107 *
108 *  **   Step limitation due to bremsstrahlung ?
109 *
110       IF (IBREM.GT.0) THEN
111          STEPBR = GEKRT1*Q(JBREM+IEKBIN) +GEKRAT*Q(JBREM+IEKBIN+1)
112          SBREM  = STEPBR*ZINTBR
113          IF (SBREM.LT.STEP) THEN
114             STEP  = SBREM
115             IPROC = 9
116          ENDIF
117       ENDIF
118 *
119 *  **   Step limitation due to delta-ray production ?
120 *
121       IF (IDRAY.GT.0) THEN
122          STEPDR = GEKRT1*Q(JDRAY+IEKBIN) +GEKRAT*Q(JDRAY+IEKBIN+1)
123          SDRAY  = STEPDR*ZINTDR
124          IF (SDRAY.LT.STEP) THEN
125             STEP  = SDRAY
126             IPROC = 10
127          ENDIF
128       ENDIF
129 *
130 *  **   Step limitation due to annihilation ?
131 *
132       IF (CHARGE.GT.0.) THEN
133          IF (IANNI.GT.0) THEN
134             STEPAN = GEKRT1*Q(JANNI+IEKBIN) +GEKRAT*Q(JANNI+IEKBIN+1)
135             SANNI  = STEPAN*ZINTAN
136             IF (SANNI.LT.STEP) THEN
137                STEP  = SANNI
138                IPROC = 11
139             ENDIF
140          ENDIF
141       ENDIF
142 *
143       IF (STEP.LE.0.) THEN
144          STEP = 0.
145          GO TO 90
146       ENDIF
147 *
148 *  **   Step limitation due to energy-loss,multiple scattering
149 *             or magnetic field ?
150 *
151       IF (JMULOF.NE.0) THEN
152          SMULOF  = GEKRT1*Q(JMULOF+IEKBIN) +GEKRAT*Q(JMULOF+IEKBIN+1)
153          IF (SMULOF.LT.STEP) THEN
154             STEP  = SMULOF
155             IPROC = 0
156          ENDIF
157       ENDIF
158 *
159 *  **   Step limitation due to geometry ?
160 *
161       IF (STEP.GE.0.95*SAFETY) THEN
162          CALL GTNEXT
163          IF (IGNEXT.NE.0) THEN
164             STEP   = SNEXT + PREC
165             IPROC = 0
166          ENDIF
167 *
168 *        Update SAFETY in stack companions, if any
169          IF (IQ(JSTAK+3).NE.0) THEN
170             DO 10 IST = IQ(JSTAK+3),IQ(JSTAK+1)
171                JST    = JSTAK + 3 + (IST-1)*NWSTAK
172                Q(JST+11) = SAFETY
173    10       CONTINUE
174             IQ(JSTAK+3) = 0
175          ENDIF
176       ELSE
177          IQ(JSTAK+3) = 0
178       ENDIF
179 *
180       IF (IFIELD.EQ.0.OR.STEP.LE.PREC) THEN
181 *
182 * *** Linear transport when no field or very short step
183 *
184          IF (IGNEXT.NE.0) THEN
185 *
186 * *** Particle is supposed to cross the boundary during step
187 *
188             DO 20 I = 1,3
189                VECTMP  = VECT(I) +STEP*VECT(I+3)
190                IF(VECTMP.EQ.VECT(I)) THEN
191 *
192 * *** Correct for machine precision
193 *
194                   IF(VECT(I+3).NE.0.) THEN
195                      VECTMP =
196      +               VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC
197                      IF(NMEC.GT.0) THEN
198                         IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
199                      ENDIF
200                      NMEC=NMEC+1
201                      LMEC(NMEC)=104
202 #if defined(CERNLIB_DEBUG)
203                      WRITE(CHMAIL, 10000)
204                      CALL GMAIL(0,0)
205                      WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
206                      CALL GMAIL(0,0)
207 10000 FORMAT(' Boundary correction in GTELEC: ',
208      +       '    GEKIN      NUMED       STEP      SNEXT')
209 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
210 #endif
211                   ENDIF
212                ENDIF
213                VECT(I) = VECTMP
214    20       CONTINUE
215             INWVOL = 2
216             NMEC = NMEC +1
217             LMEC(NMEC) = 1
218          ELSE
219             DO 30 I = 1,3
220                VECT(I)  = VECT(I) +STEP*VECT(I+3)
221    30       CONTINUE
222          ENDIF
223       ELSE
224 *
225 * *** otherwise, swim particle in magnetic field
226 *
227 #if !defined(CERNLIB_OLD)
228          if(mycoun.gt.1.and.nfmany.gt.0.and.step.ge.safety)then
229            nlevel=manyle(nfmany)
230            do 99 i=1,nlevel
231              names(i)=manyna(nfmany,i)
232              number(i)=manynu(nfmany,i)
233  99        continue
234            call glvolu(nlevel,names,number,ier)
235            if(ier.ne.0)print *,'Fatal error in GLVOLU'
236            ingoto=0
237          endif
238 #endif
239          NMEC = NMEC +1
240          LMEC(NMEC) = 4
241 *
242 #if !defined(CERNLIB_USRJMP)
243    40    CALL GUSWIM (CHARGE, STEP, VECT, VOUT)
244 #endif
245 #if defined(CERNLIB_USRJMP)
246    40    CALL JUMPT4(JUSWIM,CHARGE, STEP, VECT, VOUT)
247 #endif
248 *
249 *  ** When near to boundary, take proper action (cut-step,crossing...)
250 *
251          IF(STEP.GE.SAFETY)THEN
252             INEAR = 0
253             IF (IGNEXT.NE.0) THEN
254                DO 50 I = 1,3
255                   VNEXT(I+3) = VECT(I+3)
256                   VNEXT(I) = VECT(I) +SNEXT*VECT(I+3)
257    50          CONTINUE
258                DO 60 I = 1,3
259                   IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 70
260    60          CONTINUE
261                INEAR = 1
262             ENDIF
263 *
264    70       CALL GINVOL (VOUT, ISAME)
265             IF (ISAME.EQ.0) THEN
266                IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN
267                   INWVOL = 2
268                   NMEC = NMEC +1
269                   LMEC(NMEC) = 1
270                ELSE
271 *              Cut step
272                   STEP = 0.5*STEP
273                   IF (LMEC(NMEC).NE.24) THEN
274                      NMEC = NMEC +1
275                      LMEC(NMEC) = 24
276                   ENDIF
277                   GO TO 40
278                ENDIF
279             ENDIF
280          ENDIF
281 *
282 *
283          DO 80 I = 1,6
284             VECT(I) = VOUT(I)
285    80    CONTINUE
286 *
287       ENDIF
288 *
289 *
290 * *** Correct the step due to multiple scattering
291       IF (IMULL.NE.0) THEN
292          STMULS = STEP
293          CORR=0.0001*(STEP/RADL)*(GETOT/(VECT(7)*VECT(7)))**2
294          IF (CORR.GT.0.25) CORR = 0.25
295          STEP  = (1.+CORR)*STEP
296       ENDIF
297 *
298       SLENG = SLENG + STEP
299 *
300 *  **  Apply synchrotron radiation if required
301 *
302       IF(ISYNC*IFIELD.NE.0) THEN
303          CALL GSYNC
304          NMEC = NMEC+1
305          LMEC(NMEC) = 108
306       ENDIF
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          IF(GEKRAT.LT.0.7) THEN
324             I1 = MAX(IEKBIN-1,1)
325          ELSE
326             I1 = MIN(IEKBIN,NEKBIN-1)
327          ENDIF
328          I1 = 3*(I1-1)+1
329          XCOEF1 = Q(JCOEF+I1)
330          XCOEF2 = Q(JCOEF+I1+1)
331          XCOEF3 = Q(JCOEF+I1+2)
332          IF(XCOEF1.NE.0.) THEN
333             STOPMX = -XCOEF2+SIGN(ONE,XCOEF1)*SQRT(XCOEF2**2 - (XCOEF3-
334      +      GEKIN/XCOEF1))
335          ELSE
336             STOPMX = - (XCOEF3-GEKIN)/XCOEF2
337          ENDIF
338 *
339          STOPRG = STOPMX - STEP
340          IF (STOPRG.LT.STOPC) THEN
341             STEP = MAX(STOPMX - STOPC,ZERO)
342             GO TO 100
343          ENDIF
344 *
345          LMEC(NMEC) = 3
346          IF(XCOEF1.NE.0.) THEN
347             DEMEAN=GEKIN-XCOEF1*(XCOEF3+STOPRG*(2.*XCOEF2+STOPRG))
348          ELSE
349             DEMEAN=GEKIN-XCOEF2*STOPRG-XCOEF3
350          ENDIF
351          IF(DEMEAN.LE.5.*GEKIN*EPSMAC) THEN
352             DEMEAN=(GEKRT1*Q(JLOSS+IEKBIN)+GEKRAT*Q(JLOSS+IEKBIN+1))
353      +     *STEP
354          ENDIF
355          IF (ILOSS.EQ.4.OR.IEKBIN.LE.IKCUT+1) THEN
356             DESTEP = DEMEAN
357          ELSE
358             DEMS = DEMEAN
359             CALL GFLUCT (DEMS,DESTEP)
360          ENDIF
361          DESTEP=MAX(DESTEP,0.)
362          GEKINT = GEKIN -DESTEP
363          IF (GEKINT.LE.(1.01*CUTELE)) GO TO 100
364          DESTEL = DESTEP
365          GEKIN  = GEKINT
366          GETOT  = GEKIN +AMASS
367          VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
368          CALL GEKBIN
369 *
370 *        IABAN = 1 does not distinguish between sensitive and
371 *        non-sensitive volumes, and can stop particles above the CUTE
372 *
373          IF(IABAN.EQ.1) THEN
374 *
375 *        STOP electron/positron if range < safety AND no brems
376 *
377            IF(STOPMX.LE.SAFETY) THEN
378              IF(SBREM.GT.STOPMX) THEN
379 *          + condition in case of Cherenkov generation:
380 *          STOP if E/p >refractive index (i.e. e+/e- is below threshold)
381                IF(IMCKOV.EQ.1) THEN
382                  THRIND=GETOT/VECT(7)
383 *****            IF(THRIND.GT.Q(JINDEX+1)) GOTO 100
384                  IF(THRIND.GT.Q(JINDEX+1)) GOTO 98 
385                ELSE
386                  GOTO 98 
387                ENDIF
388              ENDIF
389            ENDIF
390            STOPRG = STOPMX - STEP
391            IF (STOPRG.LT.STOPC) THEN
392             STEP = MAX(STOPMX - STOPC,ZERO)
393             GO TO 98 
394            ENDIF
395          END IF
396 *
397 *        IABAN = 2 distinguishes between sensitive and non-sensitive
398 *        volumes.
399 *        In sensitive volumes additional tests are applied before
400 *        the particle is stopped
401 *
402          IF(IABAN.EQ.2) THEN
403            IF(ISVOL.LE.0) THEN
404              IF(STOPMX.LE.SAFETY) THEN
405 *
406 *        test for brems and annihilation only
407 *
408                IF((IBREM.GT.0).AND.(SBREM.LE.STOPMX)) GOTO 97
409                IF((CHARGE.GT.0.).AND.(IANNI.GT.0).and.
410      +                                 (SANNI.LE.STOPMX)) GOTO 97
411                GOTO 98 
412              END IF
413            ELSE
414 *
415 *        sensitive volume ---> more tests !!
416 *        is energy below TLIM (=200 keV ) ?
417 *
418              IF(GEKIN.LE.TLIM) THEN
419 *
420 *     range of the particle is the overestimated stopping range here
421 *
422                 TOSTOP=STOPMX-STEP*DESTEP/DEMEAN-STOPC
423 *
424 *        does the track remain in the actual volume?
425 *
426                 IF(TOSTOP.LE.SAFETY) THEN
427 *
428 *        is there no delta ray, brems, annihilation ?
429 *
430                   IF((IDRAY.GT.0).AND.(SDRAY.LE.TOSTOP)) GOTO 97
431                   IF((IBREM.GT.0).AND.(SBREM.LE.TOSTOP)) GOTO 97
432                   IF((CHARGE.GT.0.).AND.(IANNI.GT.0).and.
433      +                                   (SANNI.LE.TOSTOP)) GOTO 97
434 *
435 *        extra condition in case of Cherenkov generation:
436 *
437                   IF(IMCKOV.EQ.1) THEN
438                      THRIND=GETOT/VECT(7)
439 *
440 *        continue only if e+/e- below threshold
441 *
442                      IF(THRIND.LT.Q(JINDEX+1)) GOTO 97
443                    ELSE
444 *
445 *        do not make transport if this estimated range negative...
446 *
447                      IF(TOSTOP.LE.0.) GOTO 98 
448 *
449 *        estimate final position/direction of the particle
450 *        from energy loss + multiple scattering
451 *       ( multiple scattering with path length = range of the particle)
452 *
453                      ALFA=0.18*Z
454                      ALFA1=ALFA+1.
455                      TGTH2=ALFA*ALFA1/(1.2+ALFA)
456                      S=TOSTOP*SQRT(1.+TGTH2)/ALFA1
457                      TGTH=SQRT(TGTH2)
458                      THET=ATAN(TGTH)
459                      IF(THET.Lt.0.) THET=PI-THET
460 *
461 *        correct direction
462 *
463                      CT=COS(THET)
464                      ST=SQRT(1.-CT*CT)
465                      CALL GRNDM(RNDM,1)
466                      PHI=TWOPI*RNDM(1)
467 *
468                      D1=ST*COS(PHI)
469                      D2=ST*SIN(PHI)
470                      D3=CT
471                      VMM=SQRT(VECT(4)*VECT(4)+VECT(5)*VECT(5))
472                      IF(VMM.NE.0.) THEN
473                        PD1=VECT(4)/VMM
474                        PD2=VECT(5)/VMM
475                        V4=PD1*VECT(6)*D1-PD2*D2+VECT(4)*D3
476                        V5=PD2*VECT(6)*D1+PD1*D2+VECT(5)*D3
477                        V6=-VMM*D1+VECT(6)*D3
478                      ELSE
479                        V4=D1
480                        V5=D2
481                        V6=D3
482                      ENDIF
483                      VP=1./SQRT(V4*V4+V5*V5+V6*V6)
484                      VECT(4)=V4*VP
485                      VECT(5)=V5*VP
486                      VECT(6)=V6*VP
487 *
488 *       transport particle ( assuming FIELD = 0. )
489 *
490                      DO 123 I=1,3
491                      VECT(I)=VECT(I)+S*VECT(I+3)
492 123                  CONTINUE
493 *
494 *       put back into GEKIN the original value in order to have
495 *       a correct DESTEP
496 *
497                   GOTO 98 
498                 END IF
499               END IF
500             END IF
501 *++++++++++++++++++++++++++++++++++++++++++++++++++++++
502           END IF
503 97        CONTINUE
504         ENDIF
505       ENDIF
506 *<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
507 *
508 * *** Apply multiple scattering.
509 *
510       IF (IMULL.NE.0) THEN
511          NMEC = NMEC +1
512          LMEC(NMEC) = 2
513          CALL GMULTS
514       ENDIF
515 *
516 * *** Update time of flight
517 *
518       TOFG = TOFG +STEP*GETOT/(VECT(7)*CLIGHT)
519       IF (TOFG.GE.TOFMAX) THEN
520          ISTOP = 4
521          NMEC  = NMEC +1
522          LMEC(NMEC) = 22
523          GO TO 999
524       ENDIF
525 *
526 * *** Update interaction probabilities
527 *
528       IF (IBREM.GT.0)    ZINTBR = ZINTBR -STEP/STEPBR
529       IF (IDRAY.GT.0)    ZINTDR = ZINTDR -STEP/STEPDR
530       IF (CHARGE.GT.0.) THEN
531          IF (IANNI.GT.0) ZINTAN = ZINTAN -STEP/STEPAN
532       ENDIF
533 *
534 * *** Apply the selected process if any
535 *
536    90 IF (IPROC.EQ.0) GO TO 999
537       NMEC = NMEC +1
538       LMEC(NMEC) = IPROC
539 *
540 *  **   Bremsstrahlung ?
541 *
542       IF (IPROC.EQ.9) THEN
543          CALL GBREME
544 *
545 *  **   Delta ray ?
546 *
547       ELSE IF (IPROC.EQ.10) THEN
548 *
549        IF((IPART.EQ.2).OR.((IPART.EQ.3).AND.(GEKIN.GT.2.*DCUTE))) THEN
550          CALL GDRAY
551        ELSE
552          GOTO 98 
553        ENDIF
554 *
555 *  **   Positron annihilation ?
556 *
557       ELSE IF (IPROC.EQ.11) THEN
558          CALL GANNI
559  
560       ENDIF
561       GO TO 999
562 *
563 * *** Special treatment for overstopped tracks
564 *
565   98  GEKIN=GEKIN+DESTEP
566   100 DESTEP = GEKIN
567       DESTEL = DESTEP
568       GEKIN  = 0.
569       GETOT  = AMASS
570       VECT(7)= 0.
571       INWVOL = 0
572       NMEC   = NMEC +1
573       LMEC(NMEC) = 30
574       IF ((CHARGE.LT.0.).OR.(IANNI.EQ.0)) THEN
575          ISTOP = 2
576       ELSE
577          NMEC = NMEC +1
578          LMEC(NMEC) = 11
579          CALL GANNIR
580       ENDIF
581   999 IF(NGPHOT.GT.0) THEN
582          IF(ITCKOV.EQ.2.AND.ISTOP.EQ.0) THEN
583 *
584 *  The electron has produced Cerenkov photons and it is still alive
585 *  we put it in the stack and we let the photons to be tracked
586             NGKINE = NGKINE+1
587             GKIN(1,NGKINE) = VECT(4)*VECT(7)
588             GKIN(2,NGKINE) = VECT(5)*VECT(7)
589             GKIN(3,NGKINE) = VECT(6)*VECT(7)
590             GKIN(4,NGKINE) = GETOT
591             GKIN(5,NGKINE) = IPART
592             TOFD(NGKINE) = 0.
593             ISTOP = 1
594 c----put position as well
595             GPOS(1,NGKINE)=VECT(1)
596             GPOS(2,NGKINE)=VECT(2)
597             GPOS(3,NGKINE)=VECT(3)
598          ENDIF
599       ENDIF
600       END
601