source: trunk/MagicSoft/Simulation/Corsika/Mmcs/mubrem.f@ 1105

Last change on this file since 1105 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: 6.2 KB
Line 
1 SUBROUTINE MUBREM
2
3C-----------------------------------------------------------------------
4C MU(ON) BREM(SSTRAHLUNG)
5C
6C TREATES MUON BREMSSTRAHLUNG (WITHOUT POLARISATION EFFECTS)
7C IN ANALOGY WITH SUBROUTINE GBREMM FROM GEANT WRITTEN BY L. URBAN
8C EXPLANATIONS SEE CERN PROGRM LIBRARY LONG WRITEUP W5013
9C THIS SUBROUTINE IS CALLED FROM MUTRAC
10C
11C DESIGN : D. HECK IK3 FZK KARLSRUHE
12C-----------------------------------------------------------------------
13
14 IMPLICIT NONE
15*KEEP,CONST.
16 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
17 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
18*KEEP,ELABCT.
19 COMMON /ELABCT/ ELCUT
20 DOUBLE PRECISION ELCUT(4)
21*KEEP,GENER.
22 COMMON /GENER/ GEN,ALEVEL
23 DOUBLE PRECISION GEN,ALEVEL
24*KEEP,MUPART.
25 COMMON /MUPART/ AMUPAR,BCUT,CMUON,FMUBRM,FMUORG
26 DOUBLE PRECISION AMUPAR(14),BCUT,CMUON(11)
27 LOGICAL FMUBRM,FMUORG
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,POLAR.
46 COMMON /POLAR/ POLART,POLARF
47 DOUBLE PRECISION POLART,POLARF
48*KEEP,RANDPA.
49 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
50 DOUBLE PRECISION FAC,U1,U2
51 REAL RD(3000)
52 INTEGER ISEED(103,10),NSEQ
53 LOGICAL KNOR
54*KEEP,REST.
55 COMMON /REST/ CONTNE,TAR,LT
56 DOUBLE PRECISION CONTNE(3),TAR
57 INTEGER LT
58*KEEP,RUNPAR.
59 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
60 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
61 * MONIOU,MDEBUG,NUCNUC,
62 * CETAPE,
63 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
64 * N1STTR,MDBASE,
65 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
66 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
67 * ,GHEISH,GHESIG
68 COMMON /RUNPAC/ DSN,HOST,USER
69 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
70 REAL STEPFC
71 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
72 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
73 * N1STTR,MDBASE
74 INTEGER CETAPE
75 CHARACTER*79 DSN
76 CHARACTER*20 HOST,USER
77
78 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
79 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
80 * ,GHEISH,GHESIG
81*KEND.
82
83 DOUBLE PRECISION ALFA1,AUXIL,BETA1,COSTH3,COSTH4,CREJ,D,F1,
84 * EKIN,EMUON,PHI3,SCREJ,SINTH3,THETA3,U,UMAX,
85 * V,VC,VM,V1,W1,Z
86 INTEGER I
87 DATA ALFA1/0.625D0/
88C-----------------------------------------------------------------------
89
90 IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9)
91 444 FORMAT(' MUBREM: CURPAR=',1P,9E10.3)
92
93C COPY COORDINATES TO SECPAR
94 DO 11 I = 5,8
95 SECPAR(I) = CURPAR(I)
96 11 CONTINUE
97 SECPAR( 9) = GEN
98C TOTAL AND KINETIC ENERGY OF MUON
99 EMUON = PAMA(5) * GAMMA
100 EKIN = EMUON - PAMA(5)
101 IF ( EKIN .LE. BCUT ) THEN
102C MUON ENERGY IS TOO LOW TO PRODUCE BREMSSTRAHLUNG
103 SECPAR(2) = CURPAR(2)
104 GOTO 900
105 ENDIF
106 VC = BCUT/EMUON
107 VM = 1.D0 - CMUON(6+LT)/EMUON
108 IF ( VM .LE. 0.D0 ) THEN
109C MAXIMUM OF BREMSSTRAHLUNG SPECTRUM IS NEGATIVE, NO BREMSSTRAHLUNG
110 SECPAR(2) = CURPAR(2)
111 GOTO 900
112 ENDIF
113 CREJ = CMUON(3+LT)/EMUON
114
115 50 CALL RMMAR(RD,2,1)
116 V = VC*(VM/VC)**RD(1)
117 V1 = 1.D0 - V
118C COMPUTE REJECTION FUNCTION
119 F1 = CMUON(LT) - LOG(1.D0 + CREJ*V/V1)
120 SCREJ = (V1 + 0.75D0*V*V)*F1/CMUON(LT)
121 IF ( RD(2) .GT. SCREJ ) GOTO 50
122
123C PHOTON ENERGY
124 SECPAR(2) = EMUON * V
125
126C RADIATED GAMMA BELOW CUT? THEN REDUCE ENERGY OF MUON
127 IF ( SECPAR(2) .LE. ELCUT(4) ) THEN
128 GO TO 800
129 ENDIF
130
131C SET MATERIAL CONSTANTS CMUON(.) ACCORDING TO
132C TARGET INDEX LT (1=N, 2=O, 3=AR) WHICH HAS BEEN SET IN BOX2
133 IF ( LT .EQ. 1 ) THEN
134 Z = 7.D0
135 ELSEIF ( LT .EQ. 2 ) THEN
136 Z = 8.D0
137 ELSE
138 Z = 18.D0
139 ENDIF
140
141C GENERATE EMITTED PHOTON ANGLES WITH RESPECT TO MUON DIRECTION
142C PHI IS GENERATED ISOTROPICALLY AND THETA IS ASSIGNED A UNIVERSAL
143C ANGULAR DISTRIBUTION WITH D=D(Z,E,V)
144C THIS FUNCTION APPROXIMATES THE REAL DISTRIBUTION FUNCTION WHICH CAN
145C BE FOUND IN: YUNG-SU TSAI, REV. MOD. PHYS. 46(1974)815
146C +ERRATUM: REV. MOD. PHYS. 49(1977)421
147 D = 0.13D0 *(0.8D0 + 1.3D0/Z) * (100.D0 + 1.D0/EMUON) * (1.D0 + V)
148 W1 = 9.D0 / (9.D0 + D)
149 UMAX = EMUON * PI / PAMA(5)
15010 CALL RMMAR(RD,3,1)
151 IF ( RD(1) .LE. W1 ) THEN
152 BETA1 = ALFA1
153 ELSE
154 BETA1 = 3.D0 * ALFA1
155 ENDIF
156 U = -( LOG(RD(2) * RD(3)) ) / BETA1
157C CUT: THETA SHOULD BE .LE. PI !
158C THIS CONDITION DEPENDS ON E IN THE CASE OF D=CONST TOO!
159 IF ( U .GE. UMAX ) GOTO 10
160
161 THETA3 = U * PAMA(ITYPE) / EMUON
162 COSTH3 = COS( THETA3 )
163 SINTH3 = SIN( THETA3 )
164 CALL RMMAR(RD,1,1)
165
166 PHI3 = PI2 * RD(1)
167 CALL ADDANG( COSTHE,PHI, COSTH3,PHI3, SECPAR(3),SECPAR(4))
168 IF ( SECPAR(3) .GT. C(29) ) THEN
169C WRITE BREMSSTRAHLUNG PHOTON TO STACK
170 SECPAR( 1) = 1.D0
171 SECPAR(10) = H
172 CALL TSTACK
173 ENDIF
174
175C REDUCE ENERGY OF MUON
176 800 CONTINUE
177 EMUON = EMUON * V1
178 SECPAR(2) = EMUON/PAMA(5)
179
180 900 CONTINUE
181C WRITE MUON TO STACK
182 SECPAR( 1) = CURPAR(1)
183 SECPAR( 3) = CURPAR(3)
184 SECPAR( 4) = CURPAR(4)
185 SECPAR(10) = ALEVEL
186 CALL TSTACK
187
188 RETURN
189 END
Note: See TracBrowser for help on using the repository browser.