This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gtrak / gthion.F
1 #include "geant321/pilot.h"
2 *CMZ :  3.21/02 03/07/94  17.58.49  by  S.Giani
3 *-- Author :
4       SUBROUTINE GTHION
5 C.
6 C.    ******************************************************************
7 C.    *                                                                *
8 C.    *   Heavy  ion  type track. Computes step size and propagates    *
9 C.    *    particle through step.                                      *
10 C.    *                                                                *
11 C.    *   The ionisation energy loss is calculated here (mean +        *
12 C.    *       fluctuations)                                            *
13 C.    *   The fluctuations are the same for ILOSS=1,2,3 and            *
14 C.    *       there is no fluctuation for ILOSS=4.                     *
15 C.    *                                                                *
16 C.    *   ==>Called by : GTRACK                                        *
17 C.    *       Authors    R.Brun, F.Bruyant, M.Maire, L.Urban ***       *
18 C.    *                                                                *
19 C.    ******************************************************************
20 C.
21 #include "geant321/gcbank.inc"
22 #include "geant321/gccuts.inc"
23 #include "geant321/gcjloc.inc"
24 #include "geant321/gckine.inc"
25 #include "geant321/gcking.inc"
26 #include "geant321/gcmate.inc"
27 #include "geant321/gcmulo.inc"
28 #include "geant321/gconsp.inc"
29 #include "geant321/gcphys.inc"
30 #include "geant321/gcstak.inc"
31 #include "geant321/gctmed.inc"
32 #include "geant321/gctrak.inc"
33 #include "geant321/gcunit.inc"
34 #if defined(CERNLIB_USRJMP)
35 #include "geant321/gcjump.inc"
36 #endif
37  
38 #if !defined(CERNLIB_OLD)
39 #include "geant321/gcvolu.inc"
40 #include "geant321/gcvdma.inc"
41 #endif
42 #if !defined(CERNLIB_SINGLE)
43       PARAMETER (EPSMAC=1.E-6)
44       DOUBLE PRECISION GKR,DEMEAN,STOPP,STOPMX,STOPRG,STOPC,EKIPR
45       DOUBLE PRECISION ONE,XCOEF1,XCOEF2,XCOEF3,YCOEF1,YCOEF2,YCOEF3
46 #endif
47 #if defined(CERNLIB_SINGLE)
48       PARAMETER (EPSMAC=1.E-11)
49 #endif
50       PARAMETER (THRESH=0.7,ONE=1)
51       PARAMETER (TWOTHR=2*ONE/3,AMU=0.9314943)
52       PARAMETER (DME=7.84572E-8,CNORM=2.5)
53       REAL VNEXT(6)
54       DIMENSION RNDM(2)
55       SAVE RMASS,CUTPRO,IKCUT,STOPC,FACFLU,CHAR23
56 C.
57 C.    ------------------------------------------------------------------
58 *
59 * *** Particle below energy threshold ? short circuit
60 *
61       IF (GEKIN.LE.CUTHAD) GO TO 100
62 *
63 * *** Update local pointers if medium has changed
64 *
65       IF (IUPD.EQ.0) THEN
66          IUPD  = 1
67          JLOSS = LQ(JMA-3)
68          JRANG = LQ(JMA-16) + NEK1
69          JCOEF = LQ(JMA-18) + 3*NEK1
70          RMASS  = PMASS/AMASS
71          CUTPRO = MAX(CUTHAD*RMASS,ELOW(1))
72          IKCUT = GEKA*LOG10(CUTPRO) + GEKB
73          GKR   = (CUTPRO - ELOW(IKCUT))/(ELOW(IKCUT+1) - ELOW(IKCUT))
74          STOPC = (1.-GKR)*Q(JRANG+IKCUT) + GKR*Q(JRANG+IKCUT+1)
75          FACFLU = DME*(Z*DENS/A)
76          CHAR23 = ONE/CHARGE**TWOTHR
77          IF(IMCKOV.EQ.1) THEN
78             JTCKOV = LQ(JTM-3)
79             JABSCO = LQ(JTCKOV-1)
80             JEFFIC = LQ(JTCKOV-2)
81             JINDEX = LQ(JTCKOV-3)
82             JCURIN = LQ(JTCKOV-4)
83             NPCKOV = Q(JTCKOV+1)
84          ENDIF
85       ENDIF
86 *
87 * *** Compute energy dependent parameters
88 *
89       GAMASS=GETOT+AMASS
90       BET2=GEKIN*GAMASS/(GETOT*GETOT)
91       BET=SQRT(BET2)
92       W1=1.034-0.1777*EXP(-0.08114*CHARGE)
93       W2=BET*CHAR23
94       W3=121.4139*W2+0.0378*SIN(190.7165*W2)
95       CHARG1=CHARGE*(1.-W1*EXP(-W3))
96 *
97 *              the effective charge  CHARG1
98 *            can be negative only for very low energy and
99 *     for CHARGE > 20 ( very low energy : T/A < 20 keV/nucleon)
100 *              in this case short circuit
101 *
102       IF(CHARG1.LT.0.) GOTO 100
103       CHARG2=CHARG1**2
104 *
105       OMCMOL=Q(JPROB+21)*CHARG2
106       CHCMOL=Q(JPROB+25)*ABS(CHARG1)
107       IF(FIELDM.NE.0.) THEN
108          CFLD=3333.*DEGRAD*TMAXFD/ABS(FIELDM*CHARG1)
109       ELSE
110          CFLD=BIG
111       ENDIF
112 *
113 * *** Compute current step size
114 *
115       STEP   = STEMAX
116       IPROC  = 103
117       GEKRT1 = 1. -GEKRAT
118 *
119 *  **   Step limitation due to hadron interaction ?
120 *
121       IF (IHADR.GT.0) THEN
122 #if !defined(CERNLIB_USRJMP)
123          CALL GUPHAD
124 #endif
125 #if defined(CERNLIB_USRJMP)
126          CALL JUMPT0(JUPHAD)
127 #endif
128          IF (SHADR.LT.STEP) THEN
129             IF (SHADR.LE.0.) SHADR = PREC
130             STEP  = SHADR
131             IPROC = 12
132          ENDIF
133       ENDIF
134 *
135 *  ** Step limitation due to delta-ray production ?
136 *       (Cannot be tabulated easily because dependent on AMASS)
137 *
138       IF (IDRAY.GT.0) THEN
139          STEPDR = BIG
140          IF (GEKIN.GT.DCUTM) THEN
141             TMAX   = EMASS*GEKIN*GAMASS/(0.5*AMASS*AMASS+EMASS*GETOT)
142             IF (TMAX.GT.DCUTM) THEN
143                Y    = DCUTM/TMAX
144                SIG  = (1.-Y+BET2*Y*LOG(Y))/DCUTM
145 *              extra term for spin 1/2
146                IF (AMASS.GT.0.9) SIG=SIG+0.5*(TMAX-DCUTM)/(GETOT*GETOT)
147                SIG = SIG*Q(JPROB+17)*CHARG2*EMASS/BET2
148 *
149                IF (SIG.GT.0.) THEN
150                   STEPDR = 1./SIG
151                   SDRAY  = STEPDR*ZINTDR
152                   IF (SDRAY.LE.STEP) THEN
153                      STEP  = SDRAY
154                      IPROC = 10
155                   ENDIF
156                ENDIF
157             ENDIF
158          ENDIF
159       ENDIF
160 *
161       IF (STEP.LE.0.) THEN
162          STEP  = 0.
163          GO TO 110
164       ENDIF
165 *
166 *  **   Step limitation due to energy loss (stopping range) ?
167 *
168       IF (ILOSL.GT.0) THEN
169          IF(GEKRAT.LT.THRESH) THEN
170             I1 = MAX(IEKBIN-1,1)
171          ELSE
172             I1 = MIN(IEKBIN,NEKBIN-1)
173          ENDIF
174          I1 = 3*(I1-1)+1
175          XCOEF1 = Q(JCOEF+I1)
176          XCOEF2 = Q(JCOEF+I1+1)
177          XCOEF3 = Q(JCOEF+I1+2)
178          IF(XCOEF1.NE.0) THEN
179             STOPP = -XCOEF2+SIGN(ONE,XCOEF1)* SQRT(XCOEF2
180      +      **2 -(XCOEF3-GEKIN*RMASS/XCOEF1))
181          ELSE
182             STOPP = - (XCOEF3-GEKIN*RMASS)/XCOEF2
183          ENDIF
184          STOPMX = (STOPP - STOPC)/(RMASS*CHARG2)
185          IF (STOPMX.LT.MIN(STEP,STMIN)) THEN
186             STEP = STOPMX
187             IPROC = 0
188             IF(STEP.LE.0.)THEN
189                GO TO 100
190             ENDIF
191             GO TO 10
192          ENDIF
193          EKF = (1. - DEEMAX)*GEKIN*RMASS
194          IF (EKF.LT.ELOW(1)) THEN
195             EKF = ELOW(1)
196          ELSEIF (EKF.GE.ELOW(NEK1)) THEN
197             EKF = ELOW(NEK1)*0.99
198          ENDIF
199          IKF=GEKA*LOG10(EKF)+GEKB
200          GKR=(EKF-ELOW(IKF))/(ELOW(IKF+1)-ELOW(IKF))
201          IF(GKR.LT.THRESH) THEN
202             IK1 = MAX(IKF-1,1)
203          ELSE
204             IK1 = MIN(IKF,NEKBIN-1)
205          ENDIF
206          IK1 = 3*(IK1-1)+1
207          YCOEF1=Q(JCOEF+IK1)
208          YCOEF2=Q(JCOEF+IK1+1)
209          YCOEF3=Q(JCOEF+IK1+2)
210          IF(YCOEF1.NE.0.) THEN
211             SLOSP = -YCOEF2+SIGN(ONE,YCOEF1)*SQRT(YCOEF2**2- (YCOEF3-
212      +      EKF/YCOEF1))
213          ELSE
214             SLOSP = - (YCOEF3-EKF)/YCOEF2
215          ENDIF
216          SLOSP = STOPP - SLOSP
217          SLOSS = MAX(STMIN, SLOSP/(RMASS*CHARG2) )
218          IF (SLOSS.LT.STEP) THEN
219             STEP = SLOSS
220             IPROC = 0
221          ENDIF
222       ENDIF
223 *
224 *  **   Step limitation due to energy loss in magnetic field ?
225 *
226       IF (IFIELD.NE.0) THEN
227          SFIELD = CFLD*VECT(7)
228          SFIELD=MAX(SFIELD, STMIN)
229          IF (SFIELD.LT.STEP) THEN
230             STEP  = SFIELD
231             IPROC = 0
232          ENDIF
233       ENDIF
234 *
235 *  **   Step limitation due to multiple scattering ?
236 *
237       IF (IMULL.GT.0) THEN
238          SMULS=MIN(2232.*RADL*((VECT(7)**2)/(GETOT*CHARG1))**2,10.*RADL)
239          SMULS  = MAX(STMIN, SMULS )
240          IF (SMULS.LT.STEP) THEN
241             STEP  = SMULS
242             IPROC = 0
243          ENDIF
244       ENDIF
245 *
246    10 CONTINUE
247 *
248 *  **   Step limitation due to Cerenkov production ?
249 *
250       IF (IMCKOV.GT.0) THEN
251          CALL GNCKOV
252          STCKOV = MXPHOT/MAX(3.*DNDL,1E-10)
253          SMULS  = MAX(STMIN, STCKOV)
254          IF (SMULS.LT.STEP) THEN
255             STEP  = STCKOV
256             IPROC = 0
257          ENDIF
258       ENDIF
259 *
260 *  **   Step limitation due to geometry ?
261 *
262       IF (STEP.GE.0.95*SAFETY) THEN
263          CALL GTNEXT
264          IF (IGNEXT.NE.0) THEN
265             STEP   = SNEXT + PREC
266             IPROC = 0
267          ENDIF
268 *
269 *        Update SAFETY in stack companions, if any
270          IF (IQ(JSTAK+3).NE.0) THEN
271             DO 20 IST = IQ(JSTAK+3),IQ(JSTAK+1)
272                JST    = JSTAK + 3 + (IST-1)*NWSTAK
273                Q(JST+11) = SAFETY
274    20       CONTINUE
275             IQ(JSTAK+3) = 0
276          ENDIF
277       ELSE
278          IQ(JSTAK+3) = 0
279       ENDIF
280 *
281 * *** Linear transport when no field or very short step
282 *
283       IF (IFIELD.EQ.0.OR.STEP.LE.PREC) THEN
284 *
285          IF (IGNEXT.NE.0) THEN
286             DO 30 I = 1,3
287                VECTMP  = VECT(I) +STEP*VECT(I+3)
288                IF(VECTMP.EQ.VECT(I)) THEN
289 *
290 * *** Correct for machine precision
291 *
292                   IF(VECT(I+3).NE.0.) THEN
293                      VECTMP =
294      +               VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC
295                      IF(NMEC.GT.0) THEN
296                         IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
297                      ENDIF
298                      NMEC=NMEC+1
299                      LMEC(NMEC)=104
300 #if defined(CERNLIB_DEBUG)
301                      WRITE(CHMAIL, 10000)
302                      CALL GMAIL(0,0)
303                      WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
304                      CALL GMAIL(0,0)
305 10000 FORMAT(' Boundary correction in GTHION: ',
306      +       '    GEKIN      NUMED       STEP      SNEXT')
307 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X)
308 #endif
309                   ENDIF
310                ENDIF
311                VECT(I) = VECTMP
312    30       CONTINUE
313             INWVOL = 2
314             NMEC = NMEC +1
315             LMEC(NMEC) = 1
316          ELSE
317             DO 40 I = 1,3
318                VECT(I)  = VECT(I) +STEP*VECT(I+3)
319    40       CONTINUE
320          ENDIF
321       ELSE
322 *
323 * ***   otherwise, swim particle in magnetic field
324 *
325 #if !defined(CERNLIB_OLD)
326          if(mycoun.gt.1.and.nfmany.gt.0.and.step.ge.safety)then
327            nlevel=manyle(nfmany)
328            do 99 i=1,nlevel
329              names(i)=manyna(nfmany,i)
330              number(i)=manynu(nfmany,i)
331  99        continue
332            call glvolu(nlevel,names,number,ier)
333            if(ier.ne.0)print *,'Fatal error in GLVOLU'
334            ingoto=0
335          endif
336 #endif
337          NMEC = NMEC +1
338          LMEC(NMEC) = 4
339 *
340 #if !defined(CERNLIB_USRJMP)
341    50    CALL GUSWIM (CHARG1, STEP, VECT, VOUT)
342 #endif
343 #if defined(CERNLIB_USRJMP)
344    50    CALL JUMPT4(JUSWIM, CHARG1, STEP, VECT, VOUT)
345 #endif
346 *
347 *  ** When near to boundary, take proper action (cut-step,crossing...)
348 *
349          IF(STEP.GE.SAFETY)THEN
350             INEAR = 0
351             IF (IGNEXT.NE.0) THEN
352                DO 60 I = 1,3
353                   VNEXT(I+3) = VECT(I+3)
354                   VNEXT(I) = VECT(I) +SNEXT*VECT(I+3)
355    60          CONTINUE
356                DO 70 I = 1,3
357                   IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 80
358    70          CONTINUE
359                INEAR = 1
360             ENDIF
361 *
362    80       CALL GINVOL (VOUT, ISAME)
363             IF (ISAME.EQ.0)THEN
364                IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN
365                   INWVOL = 2
366                   NMEC = NMEC +1
367                   LMEC(NMEC) = 1
368                ELSE
369 *              Cut step
370                   STEP = 0.5*STEP
371                   IF (LMEC(NMEC).NE.24) THEN
372                      NMEC = NMEC +1
373                      LMEC(NMEC) = 24
374                   ENDIF
375                   GO TO 50
376                ENDIF
377             ENDIF
378          ENDIF
379 *
380          DO 90 I = 1,6
381             VECT(I) = VOUT(I)
382    90    CONTINUE
383 *
384       ENDIF
385 *
386 * *** Correct the step due to multiple scattering
387       IF (IMULL.NE.0) THEN
388          STMULS = STEP
389          CORR=0.0001*CHARG2*(STEP/RADL)*(GETOT/(VECT(7)*VECT(7)))**2
390          IF (CORR.GT.0.25) CORR = 0.25
391          STEP  = (1.+CORR)*STEP
392       ENDIF
393 *
394       SLENG = SLENG + STEP
395 *
396 * *** Generate Cherenkov photons if required
397 *
398       IF(IMCKOV.EQ.1) THEN
399          CALL GGCKOV
400          NMEC=NMEC+1
401          LMEC(NMEC)=105
402       ENDIF
403 *
404 * *** apply energy loss : find the kinetic energy corresponding
405 *      to the new stopping range = stopmx - step
406 *
407       IF (ILOSL.NE.0) THEN
408          NMEC = NMEC +1
409          LMEC(NMEC) = 3
410          STOPRG = STOPP - STEP*RMASS*CHARG2
411          IF (STOPRG.LE.STOPC) THEN
412             STEP = STOPMX
413             GO TO 100
414          ENDIF
415          IF(XCOEF1.NE.0.) THEN
416             EKIPR = XCOEF1*(XCOEF3+STOPRG*(2.*XCOEF2+STOPRG))
417          ELSE
418             EKIPR = XCOEF2*STOPRG+XCOEF3
419          ENDIF
420          DEMEAN=GEKIN - EKIPR/RMASS
421          IF(DEMEAN.LE.5.*GEKIN*EPSMAC) THEN
422             DEMEAN=(GEKRT1*Q(JLOSS+IEKBIN)+GEKRAT*Q(JLOSS+IEKBIN+1))
423      +             *STEP*CHARG2
424          ENDIF
425 *
426 *        fluctuations : differ from that of 'ordinary' hadrons
427 *
428          IF (ILOSS.EQ.4.OR.IEKBIN.LE.IKCUT+1) THEN
429             DESTEP = DEMEAN
430          ELSE
431 *
432 *     Charge exchange fluctuations + Gaussian 'Landau' fluctuations
433 *           (it is the same for ILOSS=1,2,3 !)
434 *
435             SIGMA2=CNORM*CHARG1*(1.-CHARG1/CHARGE)
436             SIGMA2=MAX(SIGMA2,0.)
437             TA = RMASS*GEKIN
438             TAM=TA/AMU
439             SIGMA2=SIGMA2+2.+TAM*(2.+TAM)
440 *
441             SIGMA2=FACFLU*CHARG2*STEP*SIGMA2
442             IF(SIGMA2.GT.0.0) THEN
443                 SIGMA=SQRT(SIGMA2)
444             ELSE
445                 SIGMA= 0.0
446             END IF
447 *
448 *     Check if we are in 'Gaussian' regime ...
449 *
450             CAPPA=(1.+TAM)/(TAM*(2.+TAM)*EMASS)
451             CAPPA=0.5*CAPPA**2*FACFLU*CHARG2*STEP
452 *
453 *     ... if not , correct SIGMA !
454  
455             IF( (CAPPA.LT.10.) .AND. (CAPPA.GT.0.0) ) THEN
456                SIGMA=SIGMA/(0.97+0.03*SQRT(10./CAPPA))
457             ENDIF
458 *
459             CALL GRNDM(RNDM,2)
460             DEFLUC=SIGMA*SIN(TWOPI*RNDM(1))*SQRT(-2.*LOG(RNDM(2)))
461             DESTEP=DEMEAN+DEFLUC
462          ENDIF
463 *
464 *     protection against negative destep
465 *
466          IF(DESTEP.LT.0.) DESTEP=DEMEAN
467 *                          IF (DESTEP.LT.0.) DESTEP = 0.
468          GEKINT = GEKIN -DESTEP
469          IF (GEKINT.LE.(1.01*CUTHAD)) GO TO 100
470          DESTEL = DESTEP
471          GEKIN  = GEKINT
472          GETOT  = GEKIN +AMASS
473          VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
474          CALL GEKBIN
475       ENDIF
476 *
477 * *** Apply multiple scattering.
478 *
479       IF (IMULL.NE.0) THEN
480          NMEC = NMEC +1
481          LMEC(NMEC) = 2
482 *      check charge dependence ...........!!!!!!!  (later..)
483          CALL GMULTS
484       ENDIF
485 *
486 * *** Update time of flight
487 *
488       SUMLIF = SUMLIF -STEP*AMASS/VECT(7)
489       TOFG   = TOFG +STEP*GETOT/(VECT(7)*CLIGHT)
490       IF (TOFG.GE.TOFMAX) THEN
491          ISTOP = 4
492          NMEC  = NMEC +1
493          LMEC(NMEC) = 22
494          GO TO 999
495       ENDIF
496 *
497 * *** Update interaction probabilities
498 *
499       IF (IHADR.GT.0) ZINTHA = ZINTHA*(1.-STEP/SHADR)
500       IF (IDRAY.GT.0) ZINTDR = ZINTDR -STEP/STEPDR
501 *
502       GO TO 110
503 *
504 *  **   Special treatment for overstopped tracks
505 *
506   100 DESTEP = GEKIN
507       DESTEL = DESTEP
508       GEKIN  = 0.
509       GETOT  = AMASS
510       VECT(7)= 0.
511       INWVOL = 0
512       ISTOP  = 2
513       NMEC = NMEC + 1
514       LMEC(NMEC) = 30
515       IF (IHADR.EQ.0) GO TO 999
516       IPROC = 12
517 *
518 * *** apply slected process if any
519 *
520   110 IF (IPROC.EQ.0) GO TO 999
521       NMEC = NMEC +1
522       LMEC(NMEC) = IPROC
523 *
524 *  **   Hadron interaction ?
525 *
526       IF (IPROC.EQ.12) THEN
527 #if !defined(CERNLIB_USRJMP)
528          CALL GUHADR
529 #endif
530 #if defined(CERNLIB_USRJMP)
531          CALL JUMPT0(JUHADR)
532 #endif
533 *   *   Check time cut-off for decays at rest
534          IF (LMEC(NMEC).EQ.5) THEN
535             TOFG   = TOFG +SUMLIF/CLIGHT
536             SUMLIF = 0.
537             IF (TOFG.GE.TOFMAX) THEN
538                NGKINE = 0
539                ISTOP  = 4
540                LMEC(NMEC) = 22
541             ENDIF
542          ENDIF
543 *
544 *  **   Delta-ray ?
545 *
546       ELSE IF (IPROC.EQ.10) THEN
547          CALL GDRAY
548       ENDIF
549   999 END