source: trunk/MagicSoft/Simulation/Corsika/Mmcs/electr.f@ 10083

Last change on this file since 10083 was 286, checked in by harald, 25 years ago
This is the start point for further developments of the Magic Monte Carlo Simulation written by Jose Carlos Gonzales. Now it is under control of one CVS repository for the whole collaboration. Everyone should use this CVS repository for further developments.
File size: 17.9 KB
Line 
1 SUBROUTINE ELECTR(IRCODE)
2C VERSION 4.00 -- 26 JAN 1986/1900
3C******************************************************************
4 DOUBLE PRECISION PEIE
5 COMMON/BOUNDS/ECUT(6),PCUT(6),VACDST
6*KEEP,BUFFS.
7 COMMON /BUFFS/ RUNH,RUNE,EVTH,EVTE,DATAB,LH
8 INTEGER MAXBUF,MAXLEN
9 PARAMETER (MAXBUF=39*7)
10 PARAMETER (MAXLEN=12)
11 REAL RUNH(MAXBUF),EVTH(MAXBUF),EVTE(MAXBUF),
12 * RUNE(MAXBUF),DATAB(MAXBUF)
13 INTEGER LH
14 CHARACTER*4 CRUNH,CRUNE,CEVTH,CEVTE
15 EQUIVALENCE (RUNH(1),CRUNH), (RUNE(1),CRUNE)
16 EQUIVALENCE (EVTH(1),CEVTH), (EVTE(1),CEVTE)
17*KEEP,CEREN1.
18 COMMON /CEREN1/ CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD,
19 * CERSIZ,LCERFI
20 DOUBLE PRECISION CERELE,CERHAD,ETADSN,WAVLGL,WAVLGU,CYIELD
21 REAL CERSIZ
22 LOGICAL LCERFI
23*KEND.
24 COMMON/ELECIN/EKELIM,ICOMP,EKE0,EKE1,CMFP0,CMFP1,RANGE0,RANGE1, XR
25 *0,TEFF0,BLCC,XCC,PICMP0(1),PICMP1(1),EICMP0(1),EICMP1(1),MPEEM(1),
26 * ESIG0(500),ESIG1(500),PSIG0(500),PSIG1(500),EDEDX0(500),EDEDX1(50
27 *0),PDEDX0(500),PDEDX1(500),EBR10(500),EBR11(500),PBR10(500),PBR11(
28 *500),PBR20(500),PBR21(500),TMXS0(500),TMXS1(500),CMFPE0(1),CMFPE1(
29 *1),CMFPP0(1),CMFPP1(1),ERANG0(1),ERANG1(1),PRANG0(1),PRANG1(1),CXC
30 *2E0(1),CXC2E1(1),CXC2P0(1),CXC2P1(1),CLXAE0(1),CLXAE1(1),CLXAP0(1)
31 *,CLXAP1(1), THR0(1,1),THR1(1,1),THR2(1,1),THRI0(1,1),THRI1(1,1),TH
32 *RI2(1,1),FSTEP(16),FSQR(16),MSMAP(200), VERT1(1000),VERT2(100,16),
33 *MSTEPS,JRMAX,MXV1, MXV2,NBLC,NRNTH,NRNTHI,BLC0,BLC1,RTHR0,RTHR1,RT
34 *HRI0,RTHRI1
35*KEEP,EPCONT.
36 COMMON/EPCONT/ EDEP,RATIO,TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,IDISC,
37 * IROLD,IRNEW,RHOFAC, EOLD,ENEW,EKE,ELKE,BETA2,GLE,
38 * TSCAT,IAUSFL
39 DOUBLE PRECISION EDEP,RATIO
40 REAL TSTEP,TUSTEP,USTEP,TVSTEP,VSTEP,RHOFAC,EOLD,ENEW,
41 * EKE,ELKE,BETA2,GLE,TSCAT
42 INTEGER IDISC,IROLD,IRNEW,IAUSFL(29)
43*KEND.
44 COMMON/GEOM/ZALTIT,BOUND(6),NEWOBS,OBSLVL(10)
45*KEEP,LONGI.
46 COMMON /LONGI/ APLONG,HLONG,PLONG,SPLONG,THSTEP,THSTPI,
47 * NSTEP,LLONGI,FLGFIT
48 DOUBLE PRECISION APLONG(0:1040,9),HLONG(0:1024),PLONG(0:1040,9),
49 * SPLONG(0:1040,9),THSTEP,THSTPI
50 INTEGER NSTEP
51 LOGICAL LLONGI,FLGFIT
52*KEEP,MAGNET.
53 COMMON /MAGNET/ BX,BZ,BVAL,BNORMC,BNORM,COSB,SINB,BLIMIT
54 DOUBLE PRECISION BX,BZ,BVAL,BNORMC
55 REAL BNORM,COSB,SINB,BLIMIT
56*KEND.
57 COMMON /MEDIA/ NMED, RLC,RLDU,RLDUI,RHO,MSGE,MGE,MSEKE,MEKE,MLEKE,
58 *MCMFP,MRANGE,IRAYLM,HBARO(6),HBAROI(6)
59 CHARACTER MEDIA*24
60 COMMON/MEDIAC/MEDIA
61 COMMON/MISC/KMPI,KMPO,DUNIT,NOSCAT,MED(6),RHOR(6),IRAYLR(6)
62 DOUBLE PRECISION PRRMMU
63 COMMON/MUON/PRRMMU,RMMU,RMMUT2
64*KEEP,OBSPAR.
65 COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
66 * THETPR,PHIPR,NOBSLV
67 DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
68 * THETAP,THETPR(2),PHIP,PHIPR(2)
69 INTEGER NOBSLV
70*KEEP,PARPAR.
71 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
72 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
73 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
74 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
75 INTEGER ITYPE,LEVL
76*KEND.
77 COMMON/PATHCM/NPTH,B0PTH,B1PTH,PTH0(6),PTH1(6),PTH2(6)
78*KEEP,RANDPA.
79 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
80 DOUBLE PRECISION FAC,U1,U2
81 REAL RD(3000)
82 INTEGER ISEED(103,10),NSEQ
83 LOGICAL KNOR
84*KEEP,REJECT.
85 COMMON /REJECT/ AVNREJ,
86 * ALTMIN,ANEXP,THICKA,THICKD,CUTLN,EONCUT,
87 * FNPRIM
88 DOUBLE PRECISION AVNREJ(10)
89 REAL ALTMIN(10),ANEXP(10),THICKA(10),THICKD(10),
90 * CUTLN,EONCUT
91 LOGICAL FNPRIM
92*KEEP,RUNPAR.
93 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
94 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
95 * MONIOU,MDEBUG,NUCNUC,
96 * CETAPE,
97 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
98 * N1STTR,MDBASE,
99 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
100 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
101 * ,GHEISH,GHESIG
102 COMMON /RUNPAC/ DSN,HOST,USER
103 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
104 REAL STEPFC
105 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
106 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
107 * N1STTR,MDBASE
108 INTEGER CETAPE
109 CHARACTER*79 DSN
110 CHARACTER*20 HOST,USER
111
112 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
113 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
114 * ,GHEISH,GHESIG
115*KEEP,STACKE.
116 COMMON/STACKE/ E,TIME,X,Y,Z,U,V,W,DNEAR,IQ,IGEN,IR,IOBS,LPCTE,NP
117 DOUBLE PRECISION E(60),TIME(60)
118 REAL X(60),Y(60),Z(60),U(60),V(60),W(60),DNEAR(60)
119 INTEGER IQ(60),IGEN(60),IR(60),IOBS(60),LPCTE(60),NP
120*KEND.
121 COMMON/THRESH/RMT2,RMSQ,ESCD2,AP,API,AE,UP,UE,TE,THMOLL
122 COMMON/UPHIIN/SINC0,SINC1,SIN0(20002),SIN1(20002)
123 COMMON/UPHIOT/THETA,SINTHE,COSTHE,SINPHI, COSPHI,PI,TWOPI,PI5D2
124 DOUBLE PRECISION PZERO,PRM,PRMT2,RMI,VC
125 COMMON/USEFUL/PZERO,PRM,PRMT2,RMI,VC,RM,MEDIUM,MEDOLD,IBLOBE,ICALL
126 COMMON/ACLOCK/NCLOCK,JCLOCK
127 DOUBLE PRECISION THICK
128 DATA NSTPCN/0/
129
130
131
132C_____NCLOCK = NCLOCK+1
133C_____IF (NCLOCK.GT.JCLOCK) THEN
134C______WRITE(MDEBUG,* )' ELECTR:NP=',NP,' IR=',IR(NP),' IOBS=',IOBS(NP)
135C______CALL AUSGB2
136C_____END IF
137 NEWOBS=IOBS(NP)
138 IRCODE=1
139 IROLD=IR(NP)
140 IRL=IR(NP)
141 MEDIUM=MED(IRL)
142380 CONTINUE
143381 CONTINUE
144 LELEC=5-2*IQ(NP)
145 PEIE=E(NP)
146 EIE=PEIE
147 IF((EIE.LE.ECUT(IRL)))GO TO 390
148 MEDIUM=MED(IRL)
149400 CONTINUE
150401 CONTINUE
151 IF (MEDIUM.NE.0) THEN
152 EKE=EIE-RM
153 ELKE=LOG(EKE)
154 CALL RMMAR(RNNE1,1,2)
155 IF ((RNNE1.EQ.0.0)) THEN
156 RNNE1=1.E-30
157 END IF
158 DEMFP=AMAX1(-ALOG(RNNE1),1.E-6)
159 LELKE=EKE1*ELKE+EKE0
160 IF (LELEC.LT.0) THEN
161 SIG0=ESIG1(LELKE)*ELKE+ESIG0(LELKE)
162 ELSE
163 SIG0=PSIG1(LELKE)*ELKE+PSIG0(LELKE)
164 END IF
165 END IF
166450 CONTINUE
167451 CONTINUE
168 IF (MEDIUM.EQ.0) THEN
169 TSTEP=VACDST
170 USTEP=TSTEP
171 TUSTEP=USTEP
172 ELSE
173 RHOFAC=RHOR(IRL)/RHO
174 RHOFI=1./RHOFAC
175 SIG=SIG0*RHOFAC
176 IF (SIG.LE.0.0) THEN
177 TSTEP=VACDST
178 ELSE
179 TSTEP=DEMFP/SIG
180 END IF
181 TMXS=TMXS1(LELKE)*ELKE+TMXS0(LELKE)
182 TMXS=MIN(TMXS,STEPFC*200.*TEFF0)
183 TMXS=TMXS*RHOFI
184 TUSTEP=MIN(TSTEP,TMXS)
185 IF (LELEC.LT.0) THEN
186 DEDX0=EDEDX1(LELKE)*ELKE+EDEDX0(LELKE)
187 ELSE
188 DEDX0=PDEDX1(LELKE)*ELKE+PDEDX0(LELKE)
189 END IF
190 DEDX=RHOFAC*MIN(DEDX0,(86.65-Z(NP)*8.E-6)*RLDUI)
191 RANGE=(EIE-ECUT(IRL)+0.001)/DEDX
192 BETA2=MAX(1.E-8,1.-RMSQ/(EIE*EIE))
193 BETA3=EIE*BETA2*0.094315
194 TSCAT=RLDU*BETA3*BETA3
195 TSCAT=TSCAT*RHOFI
196 TUSTEP=MIN(TUSTEP,0.3*TSCAT,RANGE)
197 RATIO=TUSTEP/TSCAT
198 USTEP=TUSTEP*(1.D0-RATIO)
199 USTEPU=USTEP
200 ALTEXP=EXP(-Z(NP)*HBAROI(IRL))
201 USTEP=USTEP*ALTEXP
202 DISC=W(NP)*USTEP*HBAROI(IRL)
203 IF (ABS(DISC).LT.0.065) THEN
204 USTEP=USTEP*(1.-0.5*DISC*(1.-0.6666667*DISC* (1.-0.75*DISC *
205 * (1.-0.8*DISC))))
206 ELSE IF(DISC.LE.-1.) THEN
207 USTEP=VACDST
208 ELSE
209 USTEP=USTEP/DISC*LOG(DISC+1.)
210 END IF
211 TUSTPC=USTEP/(1.D0-RATIO)
212 END IF
213 IRNEW=IR(NP)
214 IDISC=0
215 USTEP0=USTEP
216 USTEP=MIN(USTEP,BLIMIT*EIE)
217 IF((USTEP.GT.DNEAR(NP) ))CALL HOWFAR
218 IF((IDISC.GT.0))GO TO 420
219 IF (USTEP.LE.0.0) THEN
220 IF (USTEP.LT.-1.E-4) THEN
221 WRITE(KMPO,460)USTEP
222460 FORMAT(' ELECTR: NEGATIVE USTEP=',G20.10,' CM')
223 WRITE(KMPO,470)Z(NP),DNEAR(NP),IR(NP),IRNEW,W(NP)
224470 FORMAT (' Z=',G15.7, ' DNEAR=',G15.7,' IR=',I5, ' IRNEW=',I5,
225 * ' W=',G15.7)
226 NSTPCN=NSTPCN+1
227 IF (NSTPCN.GE.20) THEN
228 CALL AUSGB2
229 WRITE(KMPO,480) NSTPCN
230480 FORMAT (' ELECTR: PROGRAM STOPPED BECAUSE OF FREQUENT NEGA',
231 * 'TIVE USTEP, COUNTER = ',I5)
232 STOP
233 END IF
234 END IF
235 USTEP=0.
236 END IF
237 IF (USTEP.EQ.0.0.OR.MEDIUM.EQ.0) THEN
238 IF (USTEP.NE.0.0) THEN
239 VSTEP=USTEP
240 TVSTEP=VSTEP
241 EDEP=PZERO
242 TVSTPC=TVSTEP
243 ALPHA=VSTEP*LELEC*BNORM/EIE
244 TVSTPC=TVSTPC*(1.+0.04166667*ALPHA*ALPHA)
245 U0=U(NP)
246 V0=V(NP)
247 W0=W(NP)
248 FNORM=1.-0.5*(ALPHA*ALPHA)*(1.-0.75*(ALPHA*ALPHA))
249 F1SIN=(1.-FNORM)*SINB
250 F1COS=(1.-FNORM)*COSB
251 V1=V0*ALPHA*FNORM
252 U(NP)=U0*(1.-F1SIN*SINB)+W0*F1SIN*COSB+V1*SINB
253 V(NP)=FNORM*(V0-ALPHA*(U0*SINB-W0*COSB))
254 W(NP)=W0*(1.-F1COS*COSB)+U0*F1COS*SINB-V1*COSB
255 RADINV=1.5-0.5*(U(NP)**2+V(NP)**2+W(NP)**2)
256 U(NP)=U(NP)*RADINV
257 V(NP)=V(NP)*RADINV
258 W(NP)=W(NP)*RADINV
259 X(NP)=X(NP)+(VSTEP*0.5)*(U0+U(NP))
260 Y(NP)=Y(NP)+(VSTEP*0.5)*(V0+V(NP))
261 ZOLD =Z(NP)
262 Z(NP)=Z(NP)+(VSTEP*0.5)*(W0+W(NP))
263 TIME(NP)=TIME(NP)+TVSTPC*VC*E(NP)/SQRT((E(NP)-PRM)*(E(NP)+PRM
264 * ))
265
266
267C GENERATE CERENKOV PHOTONS
268 IF ( FNPRIM ) CALL CERENE(TVSTPC)
269C ADD ELECTRONS TO THE LONGITUDINAL DEVELOPMENT
270C FIND FIRST THE EQUIVALENT LEVELS
271 IF ( LLONGI ) THEN
272C IF STARTING POINT BELOW LOWEST LEVEL THEN DON'T CHECK
273 IF ( HLONG(NSTEP) .LE. -ZOLD ) THEN
274 LPCT1 = LPCTE(NP)
275C Z NEW IS PROBABLY ONLY LITTLE BELOW Z OLD, THEREFORE INCREMENTAL SEARCH
276 DO 6002 I1 = LPCT1,NSTEP
277 IF ( HLONG(I1) .LT. -Z(NP) ) GOTO 6003
278 6002 CONTINUE
279 I1 = NSTEP + 1
280 6003 CONTINUE
281 LPCT2 = I1 - 1
282 DO 485 I=LPCT1,LPCT2
283 PLONG(I,IQ(NP)) = PLONG(I,IQ(NP)) + 1.D0
284 485 CONTINUE
285C STORE END POINT AS POSSIBLE STARTPOINT OF NEXT TRACK
286 LPCTE(NP) = LPCT2 + 1
287 ENDIF
288 ENDIF
289 DNEAR(NP)=DNEAR(NP)-VSTEP
290 IROLD=IR(NP)
291 END IF
292 IR(NP)=IRNEW
293 IRL=IRNEW
294 MEDIUM=MED(IRL)
295 IF((EIE.LE.ECUT(IRL)))GO TO 390
296 IF (USTEP.NE.0.0) THEN
297 IF (NEWOBS.GT.IOBS(NP)) THEN
298 CALL AUSGAB
299 IOBS(NP)=NEWOBS
300 END IF
301 END IF
302 GO TO 401
303 END IF
304 VSTEP=USTEP
305 IF (USTEP.EQ.USTEP0) THEN
306 TVSTEP=TUSTEP
307 TVSTPC=TUSTPC
308 ELSE
309 VSTEPU=VSTEP
310 DISC=W(NP)*VSTEPU*HBAROI(IRL)
311 VSTEPU=VSTEPU/ALTEXP
312 IF (ABS(DISC).LT.0.16) THEN
313 VSTEPU=VSTEPU*(1.+.5*DISC*(1.+.33333333*DISC* (1.+0.25*DISC*
314 * (1.+0.2*DISC))))
315 ELSE
316 VSTEPU=VSTEPU/DISC*(EXP(DISC)-1.)
317 END IF
318 VSTP=VSTEPU/TSCAT
319 IPTH=B0PTH+B1PTH*VSTP
320 IF (IPTH.GT.NPTH) THEN
321 CALL AUSGB2
322 WRITE(KMPO,490) VSTP,IPTH,NPTH
323490 FORMAT (' ELECTR: OUT OF BOUNDS IPTH: VSTP,IPTH,NPTH=' , 1P ,
324 * G15.6,2I10)
325 STOP
326 END IF
327 PTH=PTH0(IPTH)+VSTP*(PTH1(IPTH)+VSTP*PTH2(IPTH))
328 TVSTEP=PTH*VSTEPU
329 TVSTPC=PTH*VSTEP
330 END IF
331 ALPHA=VSTEP*LELEC*BNORM/EIE
332 TVSTPC=TVSTPC*(1.+0.04166667*ALPHA*ALPHA)
333 DE=DEDX*TVSTEP
334 EDEP=DE
335 EKEF=EKE-DE
336 EOLD=EIE
337 ENEW=EOLD-DE
338 CALL MSCAT
339 U0=U(NP)
340 V0=V(NP)
341 W0=W(NP)
342 FNORM=1.-0.5*(ALPHA*ALPHA)*(1.-0.75*(ALPHA*ALPHA))
343 F1SIN=(1.-FNORM)*SINB
344 F1COS=(1.-FNORM)*COSB
345 V1=V0*ALPHA*FNORM
346 U(NP)=U0*(1.-F1SIN*SINB)+W0*F1SIN*COSB+V1*SINB
347 V(NP)=FNORM*(V0-ALPHA*(U0*SINB-W0*COSB))
348 W(NP)=W0*(1.-F1COS*COSB)+U0*F1COS*SINB-V1*COSB
349 RADINV=1.5-0.5*(U(NP)**2+V(NP)**2+W(NP)**2)
350 U(NP)=U(NP)*RADINV
351 V(NP)=V(NP)*RADINV
352 W(NP)=W(NP)*RADINV
353 X(NP)=X(NP)+(VSTEP*0.5)*(U0+U(NP))
354 Y(NP)=Y(NP)+(VSTEP*0.5)*(V0+V(NP))
355 ZOLD = Z(NP)
356 Z(NP)=Z(NP)+(VSTEP*0.5)*(W0+W(NP))
357 TIME(NP)=TIME(NP)+TVSTPC*VC*E(NP)/SQRT((E(NP)-PRM)*(E(NP)+PRM))
358
359
360C GENERATE CERENKOV PHOTONS
361 IF ( FNPRIM ) CALL CERENE(TVSTPC)
362C ADD ELECTRONS TO THE LONGITUDINAL DEVELOPMENT
363C FIND FIRST THE EQUIVALENT LEVELS
364 IF ( LLONGI ) THEN
365C IF STARTING POINT BELOW LOWEST LEVEL THEN DON'T CHECK
366 IF ( HLONG(NSTEP) .LE. -ZOLD ) THEN
367 LPCT1 = LPCTE(NP)
368C Z NEW IS PROBABLY ONLY LITTLE BELOW Z OLD, THEREFORE INCREMENTAL SEARCH
369 DO 6102 I1 = LPCT1,NSTEP
370 IF ( HLONG(I1) .LT. -Z(NP) ) GOTO 6103
371 6102 CONTINUE
372 I1 = NSTEP + 1
373 6103 CONTINUE
374 LPCT2 = I1 - 1
375 DO 495 I=LPCT1,LPCT2
376 PLONG(I,IQ(NP)) = PLONG(I,IQ(NP)) + 1.D0
377 495 CONTINUE
378 LPCTE(NP) = LPCT2 + 1
379 ENDIF
380 ENDIF
381 DNEAR(NP)=DNEAR(NP)-VSTEP
382 IROLD=IR(NP)
383 CALL RMMAR(RNNO38,1,2)
384 PHI=RNNO38*TWOPI
385 LPHI=SINC1*PHI+SINC0
386 SINPHI=SIN1(LPHI)*PHI+SIN0(LPHI)
387 CPHI=PI5D2-PHI
388 LCPHI=SINC1*CPHI+SINC0
389 COSPHI=SIN1(LCPHI)*CPHI+SIN0(LCPHI)
390 A=U(NP)
391 B=V(NP)
392 CC=W(NP)
393 SINPS2=A*A+B*B
394 IF (SINPS2.LT.1.0E-10) THEN
395 U(NP)=SINTHE*COSPHI
396 V(NP)=SINTHE*SINPHI
397 W(NP)=CC*COSTHE
398 ELSE
399 SINPSI=SQRT(SINPS2)
400 US=SINTHE*COSPHI
401 VS=SINTHE*SINPHI
402 SINDEL=B*(1./SINPSI)
403 COSDEL=A*(1./SINPSI)
404 U(NP)=CC*COSDEL*US-SINDEL*VS+A*COSTHE
405 V(NP)=CC*SINDEL*US+COSDEL*VS+B*COSTHE
406 W(NP)=-SINPSI*US+CC*COSTHE
407 END IF
408 RADINV=1.5-0.5*(U(NP)**2+V(NP)**2+W(NP)**2)
409 U(NP)=U(NP)*RADINV
410 V(NP)=V(NP)*RADINV
411 W(NP)=W(NP)*RADINV
412 PEIE=PEIE-EDEP
413 EIE=PEIE
414 E(NP)=PEIE
415 IF((EIE.LE.ECUT(IRL)))GO TO 390
416 MEDOLD=MEDIUM
417 IF (MEDIUM.NE.0) THEN
418 EKEOLD=EKE
419 EKE=EIE-RM
420 ELKE=LOG(EKE)
421 LELKE=EKE1*ELKE+EKE0
422 END IF
423 IF (IRNEW.NE.IROLD) THEN
424 IR(NP)=IRNEW
425 IRL=IRNEW
426 MEDIUM=MED(IRL)
427 END IF
428 IF((EIE.LE.ECUT(IRL)))GO TO 390
429 IF (NEWOBS.GT.IOBS(NP)) THEN
430 CALL AUSGAB
431 IOBS(NP)=NEWOBS
432 END IF
433 IF((IDISC.LT.0))GO TO 420
434 IF((MEDIUM.NE.MEDOLD))GO TO 401
435 DEMFP=MAX(0.,DEMFP-TVSTEP*SIG)
436 IF(((DEMFP.LT.1.E-6)))GO TO452
437 GO TO 451
438452 CONTINUE
439 IF (LELEC.LT.0) THEN
440 SIGF=ESIG1(LELKE)*ELKE+ESIG0(LELKE)
441 ELSE
442 SIGF=PSIG1(LELKE)*ELKE+PSIG0(LELKE)
443 END IF
444 CALL RMMAR(RFICT,1,2)
445 IF(((RFICT.LE.SIGF/SIG0)))GO TO402
446 GO TO 401
447402 CONTINUE
448 IF ( .NOT. FNPRIM ) THEN
449 X(1)=0.
450 Y(1)=0.
451 EVTH(5)=X(1)
452 EVTH(6)=-Y(1)
453 IF (FIX1I) THEN
454 Z(1)=-FIXHEI
455 NP=1
456 LPCTE(1)=MIN(NSTEP,INT(THICK(FIXHEI)*THSTPI)+1)
457 SITHET=SQRT(1.D0-SECPAR(3)**2)
458 U(1)=SITHET*COS(-SECPAR(4))
459 V(1)=SITHET*SIN(-SECPAR(4))
460 W(1)=SECPAR(3)
461 RADINV=1.5-0.5*(U(1)**2+V(1)**2+W(1)**2)
462 U(1)=U(1)*RADINV
463 V(1)=V(1)*RADINV
464 W(1)=W(1)*RADINV
465 END IF
466 EVTH(7)=-Z(1)
467 CALL TOBUF(EVTH,0)
468C OUTPUT OF EVENTHEADER TO THE CERENKOV FILE
469 IF (LCERFI) CALL TOBUFC(EVTH,0)
470 CALL COORIN(DBLE(-Z(1)))
471 TIME(1)=0.D0
472 FNPRIM =.TRUE.
473 IF (FPRINT) THEN
474 WRITE(KMPO,* )' FIRST INTERACTION AT ',EVTH(7)*0.01,' M'
475 END IF
476 END IF
477 IF (LELEC.LT.0) THEN
478 EBR1=EBR11(LELKE)*ELKE+EBR10(LELKE)
479 CALL RMMAR(RNNO24,1,2)
480 IF (RNNO24.LE.EBR1) THEN
481 GO TO 500
482 ELSE
483 IF (E(NP).LE.THMOLL) THEN
484 IF((EBR1.LE.0.0))GO TO 380
485 GO TO 500
486 END IF
487 CALL MOLLER
488 END IF
489 GO TO 380
490 END IF
491 PBR1=PBR11(LELKE)*ELKE+PBR10(LELKE)
492 CALL RMMAR(RNNO25,1,2)
493 IF((RNNO25.LT.PBR1))GO TO 500
494 PBR2=PBR21(LELKE)*ELKE+PBR20(LELKE)
495 IF (RNNO25.LT.PBR2) THEN
496 CALL BHABHA
497 ELSE
498 CALL ANNIH
499 GO TO 382
500 END IF
501 GO TO 381
502382 CONTINUE
503 RETURN
504500 CONTINUE
505 CALL BREMS
506 IF (IQ(NP).EQ.1) THEN
507 RETURN
508 ELSE
509 GO TO 380
510 END IF
511390 IF (EIE.GT.AE) THEN
512 IDR=1
513 IF (LELEC.LT.0) THEN
514 EDEP=PEIE-PRM
515 ELSE
516 EDEP=PEIE-PRM
517 END IF
518 ELSE
519 IDR=2
520 EDEP=PEIE-PRM
521 END IF
522 IF (LELEC.GT.0) THEN
523 IF (EDEP.LT.PEIE) THEN
524 CALL RMMAR(RD,2,2)
525 COSTHE=RD(1)
526 FLIP=RD(2)
527 IF((FLIP.LE.0.5))COSTHE=-COSTHE
528 SINTHE=SQRT(MAX(0.,1.0-COSTHE*COSTHE))
529 E(NP)=PRM
530 IQ(NP)=1
531 U(NP)=0.
532 V(NP)=0.
533 W(NP)=1.
534 CALL UPHI(2,1)
535 NP=NP+1
536 E(NP)=PRM
537 IQ(NP)=1
538 X(NP)=X(NP-1)
539 Y(NP)=Y(NP-1)
540 Z(NP)=Z(NP-1)
541 LPCTE(NP)=LPCTE(NP-1)
542 IR(NP)=IR(NP-1)
543 DNEAR(NP)=DNEAR(NP-1)
544 TIME(NP)=TIME(NP-1)
545 IGEN(NP)=IGEN(NP-1)
546 IOBS(NP)=IOBS(NP-1)
547 U(NP)=-U(NP-1)
548 V(NP)=-V(NP-1)
549 W(NP)=-W(NP-1)
550 RETURN
551 END IF
552 END IF
553 NP=NP-1
554 IRCODE=2
555 RETURN
556420 IF (LELEC.LT.0) THEN
557 EDEP=PEIE-PRM
558 ELSE
559 EDEP=PEIE+PRM
560 END IF
561 IRCODE=2
562 NP=NP-1
563 RETURN
564 END
Note: See TracBrowser for help on using the repository browser.