*CMZ : 24/11/95 16.28.16 by S.Ravndal *CMZ : 3.21/02 29/03/94 15.41.49 by S.Giani *-- Author : SUBROUTINE ERTRCH C. C. ****************************************************************** C. * * C. * Average charged track is extrapolated by one step * C. * * C. * ==>Called by : ERTRGO * C. * Original routine : GTHADR * C. * Authors M.Maire, E.Nagy ********* * C. * * C. ****************************************************************** C. #include "geant321/gcbank.inc" #include "geant321/gconst.inc" #include "geant321/gccuts.inc" #include "geant321/gcphys.inc" #include "geant321/gcjloc.inc" #include "geant321/gckine.inc" #include "geant321/gcmate.inc" #include "geant321/gcmulo.inc" #include "geant321/gctmed.inc" #include "geant321/gctrak.inc" #include "geant321/gcunit.inc" #include "geant321/ertrio.inc" #include "geant321/erwork.inc" #if defined(CERNLIB_SINGLE) PARAMETER (EPSMAC=1.E-11) #endif #if !defined(CERNLIB_SINGLE)&&!defined(CERNLIB_IBM) PARAMETER (EPSMAC=5.E-6) #endif #if !defined(CERNLIB_SINGLE)&&defined(CERNLIB_IBM) PARAMETER (EPSMAC=5.E-5) #endif #if !defined(CERNLIB_SINGLE) DOUBLE PRECISION GKR,DEMEAN,STOPP1,STOPP2,STOPMX,STOPRG,STOPC DOUBLE PRECISION EKIPR #endif REAL VNEXT(6) SAVE CFLD,CHARG2,RMASS,CUTPRO,IKCUT,STOPC C. C. ------------------------------------------------------------------ * * * *** Update local pointers if medium has changed * IF (IUPD.EQ.0) THEN IUPD = 1 CHARG2 = CHARGE*CHARGE IF (IPART.LE.3) THEN CUTEK = CUTELE RMASS = 1. JRANG = LQ(JMA-15) ELSE IF (IPART.LE.6) THEN CUTEK = CUTMUO RMASS = 1. JRANG = LQ(JMA-16) ELSE CUTEK = CUTHAD RMASS = PMASS/AMASS JRANG = LQ(JMA-16) + NEK1 ENDIF CUTPRO = MAX(CUTEK*RMASS,ELOW(1)) IKCUT = GEKA*LOG10(CUTPRO) + GEKB GKR = (CUTPRO - ELOW(IKCUT))/(ELOW(IKCUT+1) - ELOW(IKCUT)) STOPC = (1.-GKR)*Q(JRANG+IKCUT) + GKR*Q(JRANG+IKCUT+1) CFLD = 0. IF (FIELDM.NE.0.) CFLD = 3333.*DEGRAD*TMAXFD/ABS(FIELDM*CHARGE) ENDIF * * *** Compute current step size * STEP = BIG GEKRT1 = 1. - GEKRAT * * *** Step limitation due to energy loss (stopping range) ? * IF (ILOSS*DEEMAX.GT.0.) THEN STOPP1 = GEKRT1*Q(JRANG+IEKBIN) + GEKRAT*Q(JRANG+IEKBIN+1) STOPMX = (STOPP1 - STOPC)/(RMASS*CHARG2) EKF = (1. - BACKTR*DEEMAX)*GEKIN*RMASS IF (EKF.LT.ELOW(1)) EKF = ELOW(1) IF (EKF.GE.ELOW(NEK1)) EKF = ELOW(NEK1)*0.99 IKF=GEKA*LOG10(EKF)+GEKB GKR=(EKF-ELOW(IKF))/(ELOW(IKF+1)-ELOW(IKF)) STOPP2 = (1.-GKR)*Q(JRANG+IKF) + GKR*Q(JRANG+IKF+1) SLOSP = ABS (STOPP1 - STOPP2) STEP = SLOSP/(RMASS*CHARG2) ENDIF * * *** Step limitation due to energy loss in magnetic field ? * IF (IFIELD*FIELDM.NE.0.) THEN SFIELD = CFLD*VECT(7) IF (SFIELD.LT.STEP) STEP = SFIELD ENDIF * * *** Compute point where to store error matrix * LERST = 0 STEPER = BIG ASCL1 = BIG DO 20 IPR = 1,NEPRED STEPE = BIG IF (LELENG) STEPE = ERLENG(IPR) - SLENG IF (LEPLAN) THEN SCAL1 = 0. SCAL2 = 0. DO 18 I=1,3 SCAL1 = SCAL1 + ERPLO(I,4,IPR)*(ERPLO(I,3,IPR)-VECT(I)) SCAL2 = SCAL2 + ERPLO(I,4,IPR)*VECT(I+3) 18 CONTINUE STEPE = SCAL1/SCAL2 ENDIF IF (STEPE.LE.PREC) STEPE = BIG IF (STEPE.LT.STEPER) THEN STEPER = STEPE INLIST = IPR IF (LEPLAN) ASCL1 = ABS (SCAL1) ENDIF 20 CONTINUE IF (STEPER.LE.STEP) THEN STEP = STEPER LERST = 1 ENDIF * * *** Step limitation due to geometry ? * LNEXT = 0 IF (STEP.GE.0.95*SAFETY) THEN CALL GTNEXT IF (IGNEXT.NE.0) THEN STEP = SNEXT + PREC LNEXT = 1 IF ((STEPER-SNEXT).GT.(2*PREC)) LERST = 0 ENDIF ENDIF * * *** Linear transport when no field or very short step * IF (IFIELD.EQ.0.OR.STEP.LE.2*PREC) THEN IF (IGNEXT.NE.0) THEN DO 25 I = 1,3 VECTMP = VECT(I) +STEP*VECT(I+3) IF(VECTMP.EQ.VECT(I)) THEN * * *** Correct for machine precision * IF(VECT(I+3).NE.0.) THEN VECTMP = + VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC #if defined(DEBUG) WRITE(CHMAIL, 10000) CALL GMAIL(0,0) WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT CALL GMAIL(0,0) 10000 FORMAT(' Boundary correction in ERTRCH: ', + ' GEKIN NUMED STEP SNEXT') 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X) #endif ENDIF ENDIF VOUT(I) = VECTMP 25 CONTINUE INWVOL = 2 NMEC = NMEC +1 LMEC(NMEC) = 1 ELSE DO 30 I = 1,3 VOUT(I) = VECT(I) +STEP*VECT(I+3) 30 CONTINUE ENDIF DO 35 I = 4,6 VOUT(I) = VECT(I) 35 CONTINUE GOTO 74 END IF * * *** otherwise, swim particle in magnetic field * NMEC = NMEC +1 LMEC(NMEC) = 4 * 50 LERST = 0 LNEXT = 0 CALL GUSWIM (CHTR , STEP, VECT, VOUT) * * When near to boundary, take proper action (cut-step,crossing...) IF (STEP.GE.SAFETY) THEN INEAR = 0 IF (IGNEXT.NE.0) THEN DO 51 I = 1,3 VNEXT(I+3) = VECT(I+3) VNEXT(I) = VECT(I) +SNEXT*VECT(I+3) 51 CONTINUE DO 52 I = 1,3 IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 55 52 CONTINUE INEAR = 1 ENDIF * 55 CALL GINVOL (VOUT,ISAME) IF (ISAME.EQ.0) THEN IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN INWVOL = 2 NMEC = NMEC +1 LMEC(NMEC) = 1 LNEXT = 1 ELSE * Cut step STEP = 0.5*STEP IF (LMEC(NMEC).NE.24) THEN NMEC = NMEC +1 LMEC(NMEC) = 24 ENDIF GOTO 50 ENDIF ENDIF ENDIF * * * preset plane reached ? 74 CONTINUE IF ((LEPLAN).AND.(STEP.GE.ASCL1)) THEN SCAL3 = 0. DO 28 I=1,3 SCAL3=SCAL3+ERPLO(I,4,INLIST)*(ERPLO(I,3,INLIST)-VOUT(I)) 28 CONTINUE ASCL3 = ABS(SCAL3) SSCL1 = ASCL1/SCAL1 IF (SCAL3*SSCL1.LT. -PREC) THEN * Cut step STEP = STEP*(ASCL1/(ASCL1+ASCL3)) NMEC = NMEC +1 LMEC(NMEC) = 24 GOTO 50 ELSE IF(ASCL3.LE.PREC) LERST = 1 ENDIF ENDIF * DO 75 I=1,6 VECT(I) = VOUT(I) 75 CONTINUE * IF (LELENG.AND.(STEP.GE.STEPER)) LERST = 1 * SLENG = SLENG + STEP * * *** Now apply selected mechanisms * IF (LNEXT.EQ.1) THEN INWVOL = 2 NMEC = NMEC + 1 LMEC(NMEC) = 1 ENDIF * * *** apply energy loss : find the kinetic energy corresponding * to the new stopping range = stopmx -/+ step * (take care of the back tracking !) * IF (ILOSS*DEEMAX.GT.0) THEN NMEC = NMEC +1 LMEC(NMEC) = 3 CALL ERLAND (STEP,Z,A,DENS,VECT(7),GETOT,AMASS,DEDX2) DEDX2 = DEDX2*CHARG2*CHARG2 STOPRG = STOPP1 - BACKTR*STEP*RMASS*CHARG2 IKF = IEKBIN IF (BACKTR.LE.0.) THEN 95 IF (STOPRG.LT.Q(JRANG+IKF)) THEN IKF = IKF - 1 IF (IKF.GT.1) GO TO 95 ENDIF ELSE 96 IF (STOPRG.GE.Q(JRANG+IKF+1)) THEN IKF = IKF + 1 IF (IKF.LT.NEK1) GO TO 96 ENDIF ENDIF GKR = (STOPRG - Q(JRANG+IKF)) / (Q(JRANG+IKF+1) - Q(JRANG+IKF)) EKIPR = (1. -GKR)*ELOW(IKF) + GKR*ELOW(IKF+1) GEKINT = EKIPR/RMASS IF (GEKINT.GT.CUTEK) THEN DESTEP = ABS (GEKIN - GEKINT) GEKIN = GEKINT GETOT = GEKIN + AMASS VECT(7)= SQRT((GETOT+AMASS)*GEKIN) CALL GEKBIN ELSE DESTEP = GEKINT GEKIN = 0. GETOT = AMASS VECT(7)= 0. INWVOL = 0 ISTOP = 2 NMEC = NMEC + 1 LMEC(NMEC) = 30 ENDIF ENDIF * * *** Propagate error matrix * IF (.NOT. LEONLY) CALL ERPROP * * *** Store informations * IF(LERST.EQ.1) THEN NMEC = NMEC + 1 LMEC(NMEC) = 27 CALL ERSTOR ENDIF * END