]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | *CMZ : 24/11/95 16.28.16 by S.Ravndal |
2 | *CMZ : 3.21/02 29/03/94 15.41.49 by S.Giani | |
3 | *-- Author : | |
4 | SUBROUTINE ERTRCH | |
5 | C. | |
6 | C. ****************************************************************** | |
7 | C. * * | |
8 | C. * Average charged track is extrapolated by one step * | |
9 | C. * * | |
10 | C. * ==>Called by : ERTRGO * | |
11 | C. * Original routine : GTHADR * | |
12 | C. * Authors M.Maire, E.Nagy ********* * | |
13 | C. * * | |
14 | C. ****************************************************************** | |
15 | C. | |
16 | #include "geant321/gcbank.inc" | |
17 | #include "geant321/gconst.inc" | |
18 | #include "geant321/gccuts.inc" | |
19 | #include "geant321/gcphys.inc" | |
20 | #include "geant321/gcjloc.inc" | |
21 | #include "geant321/gckine.inc" | |
22 | #include "geant321/gcmate.inc" | |
23 | #include "geant321/gcmulo.inc" | |
24 | #include "geant321/gctmed.inc" | |
25 | #include "geant321/gctrak.inc" | |
26 | #include "geant321/gcunit.inc" | |
27 | #include "geant321/ertrio.inc" | |
28 | #include "geant321/erwork.inc" | |
29 | ||
30 | #if defined(CERNLIB_SINGLE) | |
31 | PARAMETER (EPSMAC=1.E-11) | |
32 | #endif | |
33 | #if !defined(CERNLIB_SINGLE)&&!defined(CERNLIB_IBM) | |
34 | PARAMETER (EPSMAC=5.E-6) | |
35 | #endif | |
36 | #if !defined(CERNLIB_SINGLE)&&defined(CERNLIB_IBM) | |
37 | PARAMETER (EPSMAC=5.E-5) | |
38 | #endif | |
39 | #if !defined(CERNLIB_SINGLE) | |
40 | DOUBLE PRECISION GKR,DEMEAN,STOPP1,STOPP2,STOPMX,STOPRG,STOPC | |
41 | DOUBLE PRECISION EKIPR | |
42 | #endif | |
43 | REAL VNEXT(6) | |
44 | SAVE CFLD,CHARG2,RMASS,CUTPRO,IKCUT,STOPC | |
45 | C. | |
46 | C. ------------------------------------------------------------------ | |
47 | * | |
48 | * | |
49 | * *** Update local pointers if medium has changed | |
50 | * | |
51 | IF (IUPD.EQ.0) THEN | |
52 | IUPD = 1 | |
53 | CHARG2 = CHARGE*CHARGE | |
54 | IF (IPART.LE.3) THEN | |
55 | CUTEK = CUTELE | |
56 | RMASS = 1. | |
57 | JRANG = LQ(JMA-15) | |
58 | ELSE IF (IPART.LE.6) THEN | |
59 | CUTEK = CUTMUO | |
60 | RMASS = 1. | |
61 | JRANG = LQ(JMA-16) | |
62 | ELSE | |
63 | CUTEK = CUTHAD | |
64 | RMASS = PMASS/AMASS | |
65 | JRANG = LQ(JMA-16) + NEK1 | |
66 | ENDIF | |
67 | CUTPRO = MAX(CUTEK*RMASS,ELOW(1)) | |
68 | IKCUT = GEKA*LOG10(CUTPRO) + GEKB | |
69 | GKR = (CUTPRO - ELOW(IKCUT))/(ELOW(IKCUT+1) - ELOW(IKCUT)) | |
70 | STOPC = (1.-GKR)*Q(JRANG+IKCUT) + GKR*Q(JRANG+IKCUT+1) | |
71 | CFLD = 0. | |
72 | IF (FIELDM.NE.0.) CFLD = 3333.*DEGRAD*TMAXFD/ABS(FIELDM*CHARGE) | |
73 | ENDIF | |
74 | * | |
75 | * *** Compute current step size | |
76 | * | |
77 | STEP = BIG | |
78 | GEKRT1 = 1. - GEKRAT | |
79 | * | |
80 | * *** Step limitation due to energy loss (stopping range) ? | |
81 | * | |
82 | IF (ILOSS*DEEMAX.GT.0.) THEN | |
83 | STOPP1 = GEKRT1*Q(JRANG+IEKBIN) + GEKRAT*Q(JRANG+IEKBIN+1) | |
84 | STOPMX = (STOPP1 - STOPC)/(RMASS*CHARG2) | |
85 | EKF = (1. - BACKTR*DEEMAX)*GEKIN*RMASS | |
86 | IF (EKF.LT.ELOW(1)) EKF = ELOW(1) | |
87 | IF (EKF.GE.ELOW(NEK1)) EKF = ELOW(NEK1)*0.99 | |
88 | IKF=GEKA*LOG10(EKF)+GEKB | |
89 | GKR=(EKF-ELOW(IKF))/(ELOW(IKF+1)-ELOW(IKF)) | |
90 | STOPP2 = (1.-GKR)*Q(JRANG+IKF) + GKR*Q(JRANG+IKF+1) | |
91 | SLOSP = ABS (STOPP1 - STOPP2) | |
92 | STEP = SLOSP/(RMASS*CHARG2) | |
93 | ENDIF | |
94 | * | |
95 | * *** Step limitation due to energy loss in magnetic field ? | |
96 | * | |
97 | IF (IFIELD*FIELDM.NE.0.) THEN | |
98 | SFIELD = CFLD*VECT(7) | |
99 | IF (SFIELD.LT.STEP) STEP = SFIELD | |
100 | ENDIF | |
101 | * | |
102 | * *** Compute point where to store error matrix | |
103 | * | |
104 | LERST = 0 | |
105 | STEPER = BIG | |
106 | ASCL1 = BIG | |
107 | DO 20 IPR = 1,NEPRED | |
108 | STEPE = BIG | |
109 | IF (LELENG) STEPE = ERLENG(IPR) - SLENG | |
110 | IF (LEPLAN) THEN | |
111 | SCAL1 = 0. | |
112 | SCAL2 = 0. | |
113 | DO 18 I=1,3 | |
114 | SCAL1 = SCAL1 + ERPLO(I,4,IPR)*(ERPLO(I,3,IPR)-VECT(I)) | |
115 | SCAL2 = SCAL2 + ERPLO(I,4,IPR)*VECT(I+3) | |
116 | 18 CONTINUE | |
117 | STEPE = SCAL1/SCAL2 | |
118 | ENDIF | |
119 | IF (STEPE.LE.PREC) STEPE = BIG | |
120 | IF (STEPE.LT.STEPER) THEN | |
121 | STEPER = STEPE | |
122 | INLIST = IPR | |
123 | IF (LEPLAN) ASCL1 = ABS (SCAL1) | |
124 | ENDIF | |
125 | 20 CONTINUE | |
126 | IF (STEPER.LE.STEP) THEN | |
127 | STEP = STEPER | |
128 | LERST = 1 | |
129 | ENDIF | |
130 | * | |
131 | * *** Step limitation due to geometry ? | |
132 | * | |
133 | LNEXT = 0 | |
134 | IF (STEP.GE.0.95*SAFETY) THEN | |
135 | CALL GTNEXT | |
136 | IF (IGNEXT.NE.0) THEN | |
137 | STEP = SNEXT + PREC | |
138 | LNEXT = 1 | |
139 | IF ((STEPER-SNEXT).GT.(2*PREC)) LERST = 0 | |
140 | ENDIF | |
141 | ENDIF | |
142 | * | |
143 | * *** Linear transport when no field or very short step | |
144 | * | |
145 | IF (IFIELD.EQ.0.OR.STEP.LE.2*PREC) THEN | |
146 | IF (IGNEXT.NE.0) THEN | |
147 | DO 25 I = 1,3 | |
148 | VECTMP = VECT(I) +STEP*VECT(I+3) | |
149 | IF(VECTMP.EQ.VECT(I)) THEN | |
150 | * | |
151 | * *** Correct for machine precision | |
152 | * | |
153 | IF(VECT(I+3).NE.0.) THEN | |
154 | VECTMP = | |
155 | + VECT(I)+ABS(VECT(I))*SIGN(1.,VECT(I+3))*EPSMAC | |
156 | #if defined(DEBUG) | |
157 | WRITE(CHMAIL, 10000) | |
158 | CALL GMAIL(0,0) | |
159 | WRITE(CHMAIL, 10100) GEKIN, NUMED, STEP, SNEXT | |
160 | CALL GMAIL(0,0) | |
161 | 10000 FORMAT(' Boundary correction in ERTRCH: ', | |
162 | + ' GEKIN NUMED STEP SNEXT') | |
163 | 10100 FORMAT(31X,E10.3,1X,I10,1X,E10.3,1X,E10.3,1X) | |
164 | #endif | |
165 | ENDIF | |
166 | ENDIF | |
167 | VOUT(I) = VECTMP | |
168 | 25 CONTINUE | |
169 | INWVOL = 2 | |
170 | NMEC = NMEC +1 | |
171 | LMEC(NMEC) = 1 | |
172 | ELSE | |
173 | DO 30 I = 1,3 | |
174 | VOUT(I) = VECT(I) +STEP*VECT(I+3) | |
175 | 30 CONTINUE | |
176 | ENDIF | |
177 | DO 35 I = 4,6 | |
178 | VOUT(I) = VECT(I) | |
179 | 35 CONTINUE | |
180 | GOTO 74 | |
181 | END IF | |
182 | * | |
183 | * *** otherwise, swim particle in magnetic field | |
184 | * | |
185 | NMEC = NMEC +1 | |
186 | LMEC(NMEC) = 4 | |
187 | * | |
188 | 50 LERST = 0 | |
189 | LNEXT = 0 | |
190 | CALL GUSWIM (CHTR , STEP, VECT, VOUT) | |
191 | * | |
192 | * When near to boundary, take proper action (cut-step,crossing...) | |
193 | IF (STEP.GE.SAFETY) THEN | |
194 | INEAR = 0 | |
195 | IF (IGNEXT.NE.0) THEN | |
196 | DO 51 I = 1,3 | |
197 | VNEXT(I+3) = VECT(I+3) | |
198 | VNEXT(I) = VECT(I) +SNEXT*VECT(I+3) | |
199 | 51 CONTINUE | |
200 | DO 52 I = 1,3 | |
201 | IF (ABS(VOUT(I)-VNEXT(I)).GT.EPSIL) GO TO 55 | |
202 | 52 CONTINUE | |
203 | INEAR = 1 | |
204 | ENDIF | |
205 | * | |
206 | 55 CALL GINVOL (VOUT,ISAME) | |
207 | IF (ISAME.EQ.0) THEN | |
208 | IF ((INEAR.NE.0).OR.(STEP.LT.EPSIL)) THEN | |
209 | INWVOL = 2 | |
210 | NMEC = NMEC +1 | |
211 | LMEC(NMEC) = 1 | |
212 | LNEXT = 1 | |
213 | ELSE | |
214 | * Cut step | |
215 | STEP = 0.5*STEP | |
216 | IF (LMEC(NMEC).NE.24) THEN | |
217 | NMEC = NMEC +1 | |
218 | LMEC(NMEC) = 24 | |
219 | ENDIF | |
220 | GOTO 50 | |
221 | ENDIF | |
222 | ENDIF | |
223 | ENDIF | |
224 | * | |
225 | * | |
226 | * preset plane reached ? | |
227 | 74 CONTINUE | |
228 | IF ((LEPLAN).AND.(STEP.GE.ASCL1)) THEN | |
229 | SCAL3 = 0. | |
230 | DO 28 I=1,3 | |
231 | SCAL3=SCAL3+ERPLO(I,4,INLIST)*(ERPLO(I,3,INLIST)-VOUT(I)) | |
232 | 28 CONTINUE | |
233 | ASCL3 = ABS(SCAL3) | |
234 | SSCL1 = ASCL1/SCAL1 | |
235 | IF (SCAL3*SSCL1.LT. -PREC) THEN | |
236 | * Cut step | |
237 | STEP = STEP*(ASCL1/(ASCL1+ASCL3)) | |
238 | NMEC = NMEC +1 | |
239 | LMEC(NMEC) = 24 | |
240 | GOTO 50 | |
241 | ELSE | |
242 | IF(ASCL3.LE.PREC) LERST = 1 | |
243 | ENDIF | |
244 | ENDIF | |
245 | * | |
246 | DO 75 I=1,6 | |
247 | VECT(I) = VOUT(I) | |
248 | 75 CONTINUE | |
249 | * | |
250 | IF (LELENG.AND.(STEP.GE.STEPER)) LERST = 1 | |
251 | * | |
252 | SLENG = SLENG + STEP | |
253 | * | |
254 | * *** Now apply selected mechanisms | |
255 | * | |
256 | IF (LNEXT.EQ.1) THEN | |
257 | INWVOL = 2 | |
258 | NMEC = NMEC + 1 | |
259 | LMEC(NMEC) = 1 | |
260 | ENDIF | |
261 | * | |
262 | * *** apply energy loss : find the kinetic energy corresponding | |
263 | * to the new stopping range = stopmx -/+ step | |
264 | * (take care of the back tracking !) | |
265 | * | |
266 | IF (ILOSS*DEEMAX.GT.0) THEN | |
267 | NMEC = NMEC +1 | |
268 | LMEC(NMEC) = 3 | |
269 | CALL ERLAND (STEP,Z,A,DENS,VECT(7),GETOT,AMASS,DEDX2) | |
270 | DEDX2 = DEDX2*CHARG2*CHARG2 | |
271 | STOPRG = STOPP1 - BACKTR*STEP*RMASS*CHARG2 | |
272 | IKF = IEKBIN | |
273 | IF (BACKTR.LE.0.) THEN | |
274 | 95 IF (STOPRG.LT.Q(JRANG+IKF)) THEN | |
275 | IKF = IKF - 1 | |
276 | IF (IKF.GT.1) GO TO 95 | |
277 | ENDIF | |
278 | ELSE | |
279 | 96 IF (STOPRG.GE.Q(JRANG+IKF+1)) THEN | |
280 | IKF = IKF + 1 | |
281 | IF (IKF.LT.NEK1) GO TO 96 | |
282 | ENDIF | |
283 | ENDIF | |
284 | GKR = (STOPRG - Q(JRANG+IKF)) / (Q(JRANG+IKF+1) - Q(JRANG+IKF)) | |
285 | EKIPR = (1. -GKR)*ELOW(IKF) + GKR*ELOW(IKF+1) | |
286 | GEKINT = EKIPR/RMASS | |
287 | IF (GEKINT.GT.CUTEK) THEN | |
288 | DESTEP = ABS (GEKIN - GEKINT) | |
289 | GEKIN = GEKINT | |
290 | GETOT = GEKIN + AMASS | |
291 | VECT(7)= SQRT((GETOT+AMASS)*GEKIN) | |
292 | CALL GEKBIN | |
293 | ELSE | |
294 | DESTEP = GEKINT | |
295 | GEKIN = 0. | |
296 | GETOT = AMASS | |
297 | VECT(7)= 0. | |
298 | INWVOL = 0 | |
299 | ISTOP = 2 | |
300 | NMEC = NMEC + 1 | |
301 | LMEC(NMEC) = 30 | |
302 | ENDIF | |
303 | ENDIF | |
304 | * | |
305 | * *** Propagate error matrix | |
306 | * | |
307 | IF (.NOT. LEONLY) CALL ERPROP | |
308 | * | |
309 | * *** Store informations | |
310 | * | |
311 | IF(LERST.EQ.1) THEN | |
312 | NMEC = NMEC + 1 | |
313 | LMEC(NMEC) = 27 | |
314 | CALL ERSTOR | |
315 | ENDIF | |
316 | * | |
317 | END |