1 #include "geant321/pilot.h"
2 *CMZ : 3.21/02 03/07/94 17.58.49 by S.Giani
6 C. ******************************************************************
8 C. * Heavy ion type track. Computes step size and propagates *
9 C. * particle through step. *
11 C. * The ionisation energy loss is calculated here (mean + *
13 C. * The fluctuations are the same for ILOSS=1,2,3 and *
14 C. * there is no fluctuation for ILOSS=4. *
16 C. * ==>Called by : GTRACK *
17 C. * Authors R.Brun, F.Bruyant, M.Maire, L.Urban *** *
19 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"
38 #if !defined(CERNLIB_OLD)
39 #include "geant321/gcvolu.inc"
40 #include "geant321/gcvdma.inc"
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
47 #if defined(CERNLIB_SINGLE)
48 PARAMETER (EPSMAC=1.E-11)
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)
55 SAVE RMASS,CUTPRO,IKCUT,STOPC,FACFLU,CHAR23
57 C. ------------------------------------------------------------------
59 * *** Particle below energy threshold ? short circuit
61 IF (GEKIN.LE.CUTHAD) GO TO 100
63 * *** Update local pointers if medium has changed
68 JRANG = LQ(JMA-16) + NEK1
69 JCOEF = LQ(JMA-18) + 3*NEK1
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
87 * *** Compute energy dependent parameters
90 BET2=GEKIN*GAMASS/(GETOT*GETOT)
92 W1=1.034-0.1777*EXP(-0.08114*CHARGE)
94 W3=121.4139*W2+0.0378*SIN(190.7165*W2)
95 CHARG1=CHARGE*(1.-W1*EXP(-W3))
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
102 IF(CHARG1.LT.0.) GOTO 100
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)
113 * *** Compute current step size
119 * ** Step limitation due to hadron interaction ?
122 #if !defined(CERNLIB_USRJMP)
125 #if defined(CERNLIB_USRJMP)
128 IF (SHADR.LT.STEP) THEN
129 IF (SHADR.LE.0.) SHADR = PREC
135 * ** Step limitation due to delta-ray production ?
136 * (Cannot be tabulated easily because dependent on AMASS)
140 IF (GEKIN.GT.DCUTM) THEN
141 TMAX = EMASS*GEKIN*GAMASS/(0.5*AMASS*AMASS+EMASS*GETOT)
142 IF (TMAX.GT.DCUTM) THEN
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
151 SDRAY = STEPDR*ZINTDR
152 IF (SDRAY.LE.STEP) THEN
166 * ** Step limitation due to energy loss (stopping range) ?
169 IF(GEKRAT.LT.THRESH) THEN
172 I1 = MIN(IEKBIN,NEKBIN-1)
176 XCOEF2 = Q(JCOEF+I1+1)
177 XCOEF3 = Q(JCOEF+I1+2)
179 STOPP = -XCOEF2+SIGN(ONE,XCOEF1)* SQRT(XCOEF2
180 + **2 -(XCOEF3-GEKIN*RMASS/XCOEF1))
182 STOPP = - (XCOEF3-GEKIN*RMASS)/XCOEF2
184 STOPMX = (STOPP - STOPC)/(RMASS*CHARG2)
185 IF (STOPMX.LT.MIN(STEP,STMIN)) THEN
193 EKF = (1. - DEEMAX)*GEKIN*RMASS
194 IF (EKF.LT.ELOW(1)) THEN
196 ELSEIF (EKF.GE.ELOW(NEK1)) THEN
197 EKF = ELOW(NEK1)*0.99
199 IKF=GEKA*LOG10(EKF)+GEKB
200 GKR=(EKF-ELOW(IKF))/(ELOW(IKF+1)-ELOW(IKF))
201 IF(GKR.LT.THRESH) THEN
204 IK1 = MIN(IKF,NEKBIN-1)
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-
214 SLOSP = - (YCOEF3-EKF)/YCOEF2
216 SLOSP = STOPP - SLOSP
217 SLOSS = MAX(STMIN, SLOSP/(RMASS*CHARG2) )
218 IF (SLOSS.LT.STEP) THEN
224 * ** Step limitation due to energy loss in magnetic field ?
226 IF (IFIELD.NE.0) THEN
227 SFIELD = CFLD*VECT(7)
228 SFIELD=MAX(SFIELD, STMIN)
229 IF (SFIELD.LT.STEP) THEN
235 * ** Step limitation due to multiple scattering ?
238 SMULS=MIN(2232.*RADL*((VECT(7)**2)/(GETOT*CHARG1))**2,10.*RADL)
239 SMULS = MAX(STMIN, SMULS )
240 IF (SMULS.LT.STEP) THEN
248 * ** Step limitation due to Cerenkov production ?
250 IF (IMCKOV.GT.0) THEN
252 STCKOV = MXPHOT/MAX(3.*DNDL,1E-10)
253 SMULS = MAX(STMIN, STCKOV)
254 IF (SMULS.LT.STEP) THEN
260 * ** Step limitation due to geometry ?
262 IF (STEP.GE.0.95*SAFETY) THEN
264 IF (IGNEXT.NE.0) THEN
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
281 * *** Linear transport when no field or very short step
283 IF (IFIELD.EQ.0.OR.STEP.LE.PREC) THEN
285 IF (IGNEXT.NE.0) THEN
287 VECTMP = VECT(I) +STEP*VECT(I+3)
288 IF(VECTMP.EQ.VECT(I)) THEN
290 * *** Correct for machine precision
292 IF(VECT(I+3).NE.0.) THEN
294 + VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC
296 IF(LMEC(NMEC).EQ.104) NMEC=NMEC-1
300 #if defined(CERNLIB_DEBUG)
303 WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT
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)
318 VECT(I) = VECT(I) +STEP*VECT(I+3)
323 * *** otherwise, swim particle in magnetic field
325 #if !defined(CERNLIB_OLD)
326 if(mycoun.gt.1.and.nfmany.gt.0.and.step.ge.safety)then
327 nlevel=manyle(nfmany)
329 names(i)=manyna(nfmany,i)
330 number(i)=manynu(nfmany,i)
332 call glvolu(nlevel,names,number,ier)
333 if(ier.ne.0)print *,'Fatal error in GLVOLU'
340 #if !defined(CERNLIB_USRJMP)
341 50 CALL GUSWIM (CHARG1, STEP, VECT, VOUT)
343 #if defined(CERNLIB_USRJMP)
344 50 CALL JUMPT4(JUSWIM, CHARG1, STEP, VECT, VOUT)
347 * ** When near to boundary, take proper action (cut-step,crossing...)
349 IF(STEP.GE.SAFETY)THEN
351 IF (IGNEXT.NE.0) THEN
353 VNEXT(I+3) = VECT(I+3)
354 VNEXT(I) = VECT(I) +SNEXT*VECT(I+3)
357 IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 80
362 80 CALL GINVOL (VOUT, ISAME)
364 IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN
371 IF (LMEC(NMEC).NE.24) THEN
386 * *** Correct the step due to multiple scattering
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
396 * *** Generate Cherenkov photons if required
404 * *** apply energy loss : find the kinetic energy corresponding
405 * to the new stopping range = stopmx - step
410 STOPRG = STOPP - STEP*RMASS*CHARG2
411 IF (STOPRG.LE.STOPC) THEN
415 IF(XCOEF1.NE.0.) THEN
416 EKIPR = XCOEF1*(XCOEF3+STOPRG*(2.*XCOEF2+STOPRG))
418 EKIPR = XCOEF2*STOPRG+XCOEF3
420 DEMEAN=GEKIN - EKIPR/RMASS
421 IF(DEMEAN.LE.5.*GEKIN*EPSMAC) THEN
422 DEMEAN=(GEKRT1*Q(JLOSS+IEKBIN)+GEKRAT*Q(JLOSS+IEKBIN+1))
426 * fluctuations : differ from that of 'ordinary' hadrons
428 IF (ILOSS.EQ.4.OR.IEKBIN.LE.IKCUT+1) THEN
432 * Charge exchange fluctuations + Gaussian 'Landau' fluctuations
433 * (it is the same for ILOSS=1,2,3 !)
435 SIGMA2=CNORM*CHARG1*(1.-CHARG1/CHARGE)
436 SIGMA2=MAX(SIGMA2,0.)
439 SIGMA2=SIGMA2+2.+TAM*(2.+TAM)
441 SIGMA2=FACFLU*CHARG2*STEP*SIGMA2
442 IF(SIGMA2.GT.0.0) THEN
448 * Check if we are in 'Gaussian' regime ...
450 CAPPA=(1.+TAM)/(TAM*(2.+TAM)*EMASS)
451 CAPPA=0.5*CAPPA**2*FACFLU*CHARG2*STEP
453 * ... if not , correct SIGMA !
455 IF( (CAPPA.LT.10.) .AND. (CAPPA.GT.0.0) ) THEN
456 SIGMA=SIGMA/(0.97+0.03*SQRT(10./CAPPA))
460 DEFLUC=SIGMA*SIN(TWOPI*RNDM(1))*SQRT(-2.*LOG(RNDM(2)))
464 * protection against negative destep
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
473 VECT(7)= SQRT((GETOT+AMASS)*GEKIN)
477 * *** Apply multiple scattering.
482 * check charge dependence ...........!!!!!!! (later..)
486 * *** Update time of flight
488 SUMLIF = SUMLIF -STEP*AMASS/VECT(7)
489 TOFG = TOFG +STEP*GETOT/(VECT(7)*CLIGHT)
490 IF (TOFG.GE.TOFMAX) THEN
497 * *** Update interaction probabilities
499 IF (IHADR.GT.0) ZINTHA = ZINTHA*(1.-STEP/SHADR)
500 IF (IDRAY.GT.0) ZINTDR = ZINTDR -STEP/STEPDR
504 * ** Special treatment for overstopped tracks
515 IF (IHADR.EQ.0) GO TO 999
518 * *** apply slected process if any
520 110 IF (IPROC.EQ.0) GO TO 999
524 * ** Hadron interaction ?
526 IF (IPROC.EQ.12) THEN
527 #if !defined(CERNLIB_USRJMP)
530 #if defined(CERNLIB_USRJMP)
533 * * Check time cut-off for decays at rest
534 IF (LMEC(NMEC).EQ.5) THEN
535 TOFG = TOFG +SUMLIF/CLIGHT
537 IF (TOFG.GE.TOFMAX) THEN
546 ELSE IF (IPROC.EQ.10) THEN