source: branches/start/MagicSoft/Simulation/Corsika/Mmcs/mutrac.f@ 10113

Last change on this file since 10113 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: 8.3 KB
Line 
1 SUBROUTINE MUTRAC
2
3C-----------------------------------------------------------------------
4C MU(ON) TRAC(KING)
5C
6C TRACKS THE MUON REGARDING MAX. STEP LENGTH FOR MULTIPLE SCATTERING
7C CHECKS PASSAGE THROUGH OBSERVATION LEVELS
8C IRET1=1 KILLS PARTICLE
9C IRET2=1 PARTICLE HAS BEEN CUTTED IN UPDATE
10C THIS SUBROUTINE IS CALLED FROM BOX3
11C
12C DESIGN : D. HECK IK3 FZK KARLSRUHE
13C-----------------------------------------------------------------------
14
15 IMPLICIT NONE
16*KEEP,GENER.
17 COMMON /GENER/ GEN,ALEVEL
18 DOUBLE PRECISION GEN,ALEVEL
19*KEEP,IRET.
20 COMMON /IRET/ IRET1,IRET2
21 INTEGER IRET1,IRET2
22*KEEP,LONGI.
23 COMMON /LONGI/ APLONG,HLONG,PLONG,SPLONG,THSTEP,THSTPI,
24 * NSTEP,LLONGI,FLGFIT
25 DOUBLE PRECISION APLONG(0:1040,9),HLONG(0:1024),PLONG(0:1040,9),
26 * SPLONG(0:1040,9),THSTEP,THSTPI
27 INTEGER NSTEP
28 LOGICAL LLONGI,FLGFIT
29*KEEP,MUPART.
30 COMMON /MUPART/ AMUPAR,BCUT,CMUON,FMUBRM,FMUORG
31 DOUBLE PRECISION AMUPAR(14),BCUT,CMUON(11)
32 LOGICAL FMUBRM,FMUORG
33*KEEP,NPARTI.
34 COMMON /NPARTI/ NPARTO,MUOND
35 DOUBLE PRECISION NPARTO(10,25),NPHOTO(10),NPOSIT(10),NELECT(10),
36 * NNU(10),NMUP(10),NMUM(10),NPI0(10),NPIP(10),
37 * NPIM(10),NK0L(10),NKPL(10),NKMI(10),NNEUTR(10),
38 * NPROTO(10),NPROTB(10),NK0S(10),NHYP(10),
39 * NNEUTB(10),NDEUT(10),NTRIT(10),NALPHA(10),
40 * NOTHER(10),MUOND
41 EQUIVALENCE (NPARTO(1, 1),NPHOTO(1)), (NPARTO(1, 2),NPOSIT(1)),
42 * (NPARTO(1, 3),NELECT(1)), (NPARTO(1, 4),NNU(1)) ,
43 * (NPARTO(1, 5),NMUP(1)) , (NPARTO(1, 6),NMUM(1)) ,
44 * (NPARTO(1, 7),NPI0(1)) , (NPARTO(1, 8),NPIP(1)) ,
45 * (NPARTO(1, 9),NPIM(1)) , (NPARTO(1,10),NK0L(1)) ,
46 * (NPARTO(1,11),NKPL(1)) , (NPARTO(1,12),NKMI(1)) ,
47 * (NPARTO(1,13),NNEUTR(1)), (NPARTO(1,14),NPROTO(1)),
48 * (NPARTO(1,15),NPROTB(1)), (NPARTO(1,16),NK0S(1)) ,
49 * (NPARTO(1,18),NHYP(1)) , (NPARTO(1,19),NDEUT(1)) ,
50 * (NPARTO(1,20),NTRIT(1)) , (NPARTO(1,21),NALPHA(1)),
51 * (NPARTO(1,22),NOTHER(1)), (NPARTO(1,25),NNEUTB(1))
52*KEEP,OBSPAR.
53 COMMON /OBSPAR/ OBSLEV,THCKOB,XOFF,YOFF,THETAP,PHIP,
54 * THETPR,PHIPR,NOBSLV
55 DOUBLE PRECISION OBSLEV(10),THCKOB(10),XOFF(10),YOFF(10),
56 * THETAP,THETPR(2),PHIP,PHIPR(2)
57 INTEGER NOBSLV
58*KEEP,PAM.
59 COMMON /PAM/ PAMA,SIGNUM
60 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
61*KEEP,PARPAR.
62 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
63 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
64 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
65 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
66 INTEGER ITYPE,LEVL
67*KEEP,PARPAE.
68 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
69 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
70 * (CURPAR(4), PHI ), (CURPAR(5), H ),
71 * (CURPAR(6), T ), (CURPAR(7), X ),
72 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
73 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
74 * (CURPAR(12),ECM )
75*KEEP,RUNPAR.
76 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
77 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
78 * MONIOU,MDEBUG,NUCNUC,
79 * CETAPE,
80 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
81 * N1STTR,MDBASE,
82 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
83 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
84 * ,GHEISH,GHESIG
85 COMMON /RUNPAC/ DSN,HOST,USER
86 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
87 REAL STEPFC
88 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
89 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
90 * N1STTR,MDBASE
91 INTEGER CETAPE
92 CHARACTER*79 DSN
93 CHARACTER*20 HOST,USER
94
95 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
96 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
97 * ,GHEISH,GHESIG
98*KEND.
99
100 DOUBLE PRECISION CHITOT,HEIGH,HNEW,PROPAR(8),THCKHN
101 INTEGER I,IRET3,J,L,LPCT1,LPCT2
102 LOGICAL FSCAT
103 EXTERNAL HEIGH
104C-----------------------------------------------------------------------
105
106 IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9)
107 444 FORMAT(' MUTRAC: CURPAR=',1P,9E10.3)
108
109C THE PLACE OF NEXT INTERACTION WAS DETERMINED IN BOX2
110C KEEP TOTAL STEP LENGTH UNTIL DECAY OR INTERACTION OCCURS
111 CHITOT = CHI
112
113 10 CONTINUE
114
115C CALCULATE MAX STEP SIZE (10 RAD. LENGTH) FOR MULTIPLE SCATTERING
116 CHI = MIN( C(20), CHITOT )
117 IF ( CHI .EQ. CHITOT ) THEN
118 FSCAT = .FALSE.
119 IF (DEBUG) WRITE(MDEBUG,*)'MUTRAC: CHI=',SNGL(CHI)
120 ELSE
121 FSCAT = .TRUE.
122 IF (DEBUG) WRITE(MDEBUG,*)'MUTRAC: C(20)=',SNGL(C(20))
123 ENDIF
124
125
126C CALCULATE HIGHT DIFFERENCE IN CM FROM GIVEN CHI IN G/CM**2
127 THCKHN = THICKH + COSTHE * CHI
128 HNEW = HEIGH(THCKHN)
129 IF (DEBUG) WRITE(MDEBUG,*)'MUTRAC: THICKH,THCKHN,HNEW=',
130 * SNGL(THICKH),SNGL(THCKHN),SNGL(HNEW)
131C UPDATE MUON TO INTERACTION POINT (IF IT REACHES SO FAR)
132C AND STORE COORDINATES IN PROPAR
133 CALL UPDATE( HNEW, THCKHN, 0 )
134 IF ( DEBUG ) THEN
135 WRITE(MDEBUG,455) IRET1,IRET2
136 455 FORMAT(' MUTRAC: IRET1..2=',2I5)
137 IF ( IRET2 .EQ. 0 ) WRITE(MDEBUG,454) (OUTPAR(I),I=1,8)
138 454 FORMAT(' MUTRAC: OUTPAR=',1P,8E10.3)
139 ENDIF
140C STORE MUON FOR FURTHER TREATMENT
141 IF ( IRET2 .EQ. 0 ) THEN
142 DO 3 I = 1,8
143 PROPAR(I) = OUTPAR(I)
144 3 CONTINUE
145 IRET3 = 0
146 ELSE
147C MUON CUTTED AT INTERACTION POINT; IT MAY HOWEVER PASS SOME OF THE
148C OBSERVATION LEVELS
149 IRET3 = 1
150 ENDIF
151
152C HERE THE ENDPOINT OF THE CURRENT TRACKING STEP IS WELL DEFINED.
153C THE MUON IS TRACKED FROM THICKH DOWN TO THICKHN
154C COUNT THE MUONS FOR THE LONGITUDINAL DEVELOPMENT
155 IF ( LLONGI ) THEN
156 LPCT1 = INT(THICKH*THSTPI + 1.D0)
157 LPCT2 = INT(THCKHN*THSTPI)
158 LPCT2 = MIN(NSTEP,LPCT2)
159 IF ( ITYPE .EQ. 6 ) THEN
160 DO 5003 L = LPCT1,LPCT2
161 PLONG(L,4) = PLONG(L,4) + 1.D0
162 5003 CONTINUE
163 ELSEIF ( ITYPE .EQ. 5 ) THEN
164 DO 5013 L = LPCT1,LPCT2
165 PLONG(L,5) = PLONG(L,5) + 1.D0
166 5013 CONTINUE
167 ENDIF
168 ENDIF
169
170C CHECK OBSERVATION LEVEL PASSAGE AND UPDATE MUON COORDINATES
171 DO 1 J = 1,NOBSLV
172 IF ( HNEW .GT. OBSLEV(J) ) GOTO 2
173 IF ( H .LT. OBSLEV(J) ) GOTO 1
174C REMEMBER NUMBER OF LEVEL FOR OUTPUT
175 LEVL = J
176 CALL UPDATE( OBSLEV(J), THCKOB(J), J )
177 IF (DEBUG) WRITE(MDEBUG,456) J,IRET1,IRET2
178 456 FORMAT(' MUTRAC: OBSLEV=',I5,' IRET1,2=',2I5)
179
180C IF MUON IS NOT CUTTED, BRING IT TO OUTPUT
181 IF ( IRET2 .EQ. 0 ) THEN
182 CALL OUTPUT
183 ENDIF
184 1 CONTINUE
185
186C KILL MUON AS IT DECAYS OR INTERACTS BELOW LOWEST OBSLEVEL
187 IRET1 = 1
188 FMUORG = .FALSE.
189 RETURN
190
191C MUON SCATTERS, DECAYS OR INTERACTS BEFORE PASSING OBSLEVEL
192 2 CONTINUE
193
194 IF ( IRET3 .NE. 0 ) THEN
195C ELIMINATE MUON IF BELOW CUTS
196 IRET1 = 1
197 FMUORG = .FALSE.
198 RETURN
199 ENDIF
200C MUON IS NOW UPDATED TO POINT OF INTERACTION
201 DO 5 J = 1,8
202 CURPAR(J) = PROPAR(J)
203 5 CONTINUE
204 BETA = SQRT( GAMMA**2 - 1.D0 ) / GAMMA
205 IF ( FSCAT ) THEN
206C MUON HAS MADE MULTIPLE SCATTERING AND MUST NOW BE TRACKED FURTHER ON
207 CHITOT = CHITOT - C(20)
208 IF ( CHITOT .GT. 0.D0 ) THEN
209 THICKH = THCKHN
210 IF ( DEBUG ) WRITE(MDEBUG,457) (CURPAR(I),I=1,9)
211 457 FORMAT(' MUTRAC: SCATTER',1P,9E10.3)
212 GOTO 10
213 ENDIF
214 ENDIF
215C MUONS HAVE TO DECAY IMMEDIATELY OR TO UNDERGO BREMSSTR./PAIRPR.
216 IF ( FDECAY ) THEN
217 ALEVEL = H
218 CALL MUDECY
219 MUOND = MUOND + 1.D0
220 FMUORG = .FALSE.
221C MUDECY WRITES EM-PARTICLE TO STACK
222 ELSE
223 IF ( FMUBRM ) THEN
224 CALL MUBREM
225 ELSE
226 CALL MUPRPR
227 ENDIF
228C MUBREM AND MUPRPR WRITE EM-PARTICLES AND MUON TO STACK
229 ENDIF
230 IRET1 = 1
231
232 RETURN
233 END
Note: See TracBrowser for help on using the repository browser.