source: trunk/MagicSoft/Simulation/Corsika/Mmcs/mmolie.f@ 526

Last change on this file since 526 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: 7.3 KB
Line 
1 SUBROUTINE MMOLIE(OMEGA,DENS,VSCAT)
2
3C-----------------------------------------------------------------------
4C M(UON) MOLIE(RE MULTIPLE SCATTERING)
5C
6C TREATES MOLIERE MULTIPLE SCATTERING FOR MUONS
7C CORRECTED FOR FINITE ANGLE SCATTERING
8C THIS SUBROUTINE IS IN ANALOGY WITH SUBROUTINE GMOLIE
9C (AUTHOR: M.S.DIXIT, NRCC, OTTAWA) OF GEANT321
10C SEE CERN PROGRAM LIBRARY LONG WRITEUP W5013
11C THIS SUBROUTINE IS CALLED FROM UPDATE
12C ARGUMENTS:
13C OMEGA = NUMBER OF SCATTERINGS FOR THE STEP
14C DENS = LOCAL DENSITY
15C VSCAT = SCATTERING ANGLE
16C
17C REDESIGN: D. HECK IK3 FZK KARLSRUHE
18C-----------------------------------------------------------------------
19
20 IMPLICIT NONE
21*KEEP,CONST.
22 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
23 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
24*KEEP,MUMULT.
25 COMMON /MUMULT/ CHC,OMC,FMOLI
26 DOUBLE PRECISION CHC,OMC
27 LOGICAL FMOLI
28*KEEP,PAM.
29 COMMON /PAM/ PAMA,SIGNUM
30 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
31*KEEP,PARPAR.
32 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
33 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
34 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
35 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
36 INTEGER ITYPE,LEVL
37*KEEP,PARPAE.
38 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
39 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
40 * (CURPAR(4), PHI ), (CURPAR(5), H ),
41 * (CURPAR(6), T ), (CURPAR(7), X ),
42 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
43 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
44 * (CURPAR(12),ECM )
45*KEEP,RANDPA.
46 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
47 DOUBLE PRECISION FAC,U1,U2
48 REAL RD(3000)
49 INTEGER ISEED(103,10),NSEQ
50 LOGICAL KNOR
51*KEEP,RUNPAR.
52 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
53 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
54 * MONIOU,MDEBUG,NUCNUC,
55 * CETAPE,
56 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
57 * N1STTR,MDBASE,
58 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
59 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
60 * ,GHEISH,GHESIG
61 COMMON /RUNPAC/ DSN,HOST,USER
62 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
63 REAL STEPFC
64 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
65 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
66 * N1STTR,MDBASE
67 INTEGER CETAPE
68 CHARACTER*79 DSN
69 CHARACTER*20 HOST,USER
70
71 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
72 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
73 * ,GHEISH,GHESIG
74*KEND.
75
76 DOUBLE PRECISION TINT(40),B,BINV,CHIC,CNST,DB,DENS,OMEGA,SINTH,
77 * TEST,TMP,VSCAT
78 REAL ARG(4),F0I(40),F1I(40),F2I(40),
79 * THRED(40),VAL(4),DIN(3),F,THRI,XINT
80 INTEGER IER,JA,L,M,NA,NA3,NA3M,NMAX
81 DATA THRED/ 0.00, 0.10, 0.20, 0.30
82 + , 0.40, 0.50, 0.60, 0.70
83 + , 0.80, 0.90, 1.00, 1.10
84 + , 1.20, 1.30, 1.40, 1.50
85 + , 1.60, 1.70, 1.80, 1.90
86 + , 2.00, 2.20, 2.40, 2.60
87 + , 2.80, 3.00, 3.20, 3.40
88 + , 3.60, 3.80, 4.00, 5.00
89 + , 6.00, 7.00, 8.00, 9.00
90 + , 10.00,11.00,12.00,13.00 /
91 DATA F0I/
92 + 0.000000E+00 ,0.995016E-02 ,0.392106E-01 ,0.860688E-01
93 + ,0.147856E+00 ,0.221199E+00 ,0.302324E+00 ,0.387374E+00
94 + ,0.472708E+00 ,0.555142E+00 ,0.632121E+00 ,0.701803E+00
95 + ,0.763072E+00 ,0.815480E+00 ,0.859142E+00 ,0.894601E+00
96 + ,0.922695E+00 ,0.944424E+00 ,0.960836E+00 ,0.972948E+00
97 + ,0.981684E+00 ,0.992093E+00 ,0.996849E+00 ,0.998841E+00
98 + ,0.999606E+00 ,0.999877E+00 ,0.999964E+00 ,0.999990E+00
99 + ,0.999998E+00 ,0.999999E+00 ,0.100000E+01 ,0.100000E+01
100 + ,0.100000E+01 ,0.100000E+01 ,0.100000E+01 ,0.100000E+01
101 + ,1.000000E+00 ,1.000000E+00 ,1.000000E+00 ,1.000000E+00 /
102 DATA F1I/
103 + 0.000000E+00, 0.414985E-02, 0.154894E-01, 0.310312E-01
104 + , 0.464438E-01, 0.569008E-01, 0.580763E-01, 0.468264E-01
105 + , 0.217924E-01,-0.163419E-01,-0.651205E-01,-0.120503E+00
106 + ,-0.178272E+00,-0.233580E+00,-0.282442E+00,-0.321901E+00
107 + ,-0.350115E+00,-0.366534E+00,-0.371831E+00,-0.367378E+00
108 + ,-0.354994E+00,-0.314803E+00,-0.266539E+00,-0.220551E+00
109 + ,-0.181546E+00,-0.150427E+00,-0.126404E+00,-0.107830E+00
110 + ,-0.933106E-01,-0.817375E-01,-0.723389E-01,-0.436650E-01
111 + ,-0.294700E-01,-0.212940E-01,-0.161406E-01,-0.126604E-01
112 + ,-0.102042E-01,-0.840465E-02,-0.704261E-02,-0.598886E-02/
113 DATA F2I/
114 + 0.000000 , 0.121500E-01, 0.454999E-01, 0.913000E-01
115 + , 0.137300E+00, 0.171400E+00, 0.183900E+00, 0.170300E+00
116 + , 0.132200E+00, 0.763000E-01, 0.126500E-01,-0.473500E-01
117 + ,-0.936000E-01,-0.119750E+00,-0.123450E+00,-0.106300E+00
118 + ,-0.732800E-01,-0.312400E-01, 0.128450E-01, 0.528800E-01
119 + , 0.844100E-01, 0.114710E+00, 0.106200E+00, 0.765830E-01
120 + , 0.435800E-01, 0.173950E-01, 0.695001E-03,-0.809500E-02
121 + ,-0.117355E-01,-0.125449E-01,-0.120280E-01,-0.686530E-02
122 + ,-0.385275E-02,-0.231115E-02,-0.147056E-02,-0.982480E-03
123 + ,-0.682440E-03,-0.489715E-03,-0.361190E-03,-0.272582E-03/
124C-----------------------------------------------------------------------
125
126 IF ( DEBUG ) WRITE(MDEBUG,*)'MMOLIE: OMEGA=',SNGL(OMEGA),
127 * ' DENS=',SNGL(DENS)
128
129C COMPUTE VSCAT ANGLE FROM MOLIERE DISTRIBUTION
130 VSCAT = 0.D0
131 IF ( OMEGA .LE. ENEPER ) RETURN
132 CNST = LOG(OMEGA)
133 B = 5.D0
134 DO 10 L = 1,10
135 IF ( ABS(B) .LT. 1.D-10 ) THEN
136 B = 1.D-10
137 ENDIF
138 DB = - (B - LOG(ABS(B)) - CNST)/(1.D0 - 1.D0/B)
139 B = B + DB
140 IF ( ABS(DB) .LE. 0.0001D0 ) GOTO 20
141 10 CONTINUE
142 RETURN
143 20 CONTINUE
144 IF ( B .LE. 0.D0 ) RETURN
145C CHC IS DEFINED DIFFERENTLY FROM GEANT
146 CHIC = CHC*SQRT(CHI)/(PAMA(5)*GAMMA*BETA**2)
147 BINV = 1.D0/B
148 TINT(1) = 0.D0
149 DO 30 JA = 2,4
150 TINT(JA) = F0I(JA) + ( F1I(JA) + F2I(JA)*BINV ) * BINV
151 30 CONTINUE
152 NMAX = 4
153 40 CONTINUE
154 CALL RMMAR(RD,2,1)
155 XINT = RD(2)
156 DO 50 NA = 3,40
157 IF ( NA .GT. NMAX ) THEN
158 TINT(NA) = F0I(NA) + ( F1I(NA) + F2I(NA)*BINV ) * BINV
159 NMAX = NA
160 ENDIF
161 IF ( XINT .LE. TINT(NA-1) ) GOTO 60
162 50 CONTINUE
163 IF ( XINT .LE. TINT(40) ) THEN
164 NA = 40
165 GOTO 60
166 ELSE
167 TMP = 1.D0 - ( 1.D0 - B*(1.D0-XINT) )**5
168 IF ( TMP .LE. 0.D0 ) GOTO 40
169 THRI = 5.D0 / TMP
170 GOTO 80
171 ENDIF
172 60 CONTINUE
173 NA = MAX(NA-1,3)
174 NA3 = NA-3
175 DO 70 M = 1,4
176 NA3M = NA3 + M
177 ARG(M) = TINT(NA3M)
178 VAL(M) = THRED(NA3M)**2
179 70 CONTINUE
180 F = THRED(NA) * .02D0
181 CALL MMOL4(THRI,XINT,VAL,ARG,F,IER)
182 80 CONTINUE
183 VSCAT = CHIC * SQRT( ABS(B*THRI) )
184 IF ( VSCAT .GT. PI ) GOTO 40
185 SINTH = SIN(VSCAT)
186 TEST = VSCAT * (RD(1))**2
187 IF ( TEST .GT. SINTH ) GOTO 40
188
189 RETURN
190 END
Note: See TracBrowser for help on using the repository browser.