source: trunk/MagicSoft/Simulation/Corsika/Mmcs/mudecy.f

Last change on this file 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: 4.8 KB
Line 
1 SUBROUTINE MUDECY
2
3C-----------------------------------------------------------------------
4C MU(ON) DEC(A)Y
5C
6C TREATES DECAY OF MUON INTO ELECTRON (INCLUDING POLARISATION EFFECTS)
7C INCLUDING NEUTRINOS, IF SELECTED
8C THIS SUBROUTINE IS CALLED FROM MUTRAC
9C
10C DESIGN : D. HECK IK3 FZK KARLSRUHE
11C-----------------------------------------------------------------------
12
13 IMPLICIT NONE
14*KEEP,CONST.
15 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
16 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
17*KEEP,GENER.
18 COMMON /GENER/ GEN,ALEVEL
19 DOUBLE PRECISION GEN,ALEVEL
20*KEEP,PAM.
21 COMMON /PAM/ PAMA,SIGNUM
22 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
23*KEEP,PARPAR.
24 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
25 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
26 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
27 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
28 INTEGER ITYPE,LEVL
29*KEEP,PARPAE.
30 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
31 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
32 * (CURPAR(4), PHI ), (CURPAR(5), H ),
33 * (CURPAR(6), T ), (CURPAR(7), X ),
34 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
35 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
36 * (CURPAR(12),ECM )
37*KEEP,POLAR.
38 COMMON /POLAR/ POLART,POLARF
39 DOUBLE PRECISION POLART,POLARF
40*KEEP,RANDPA.
41 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
42 DOUBLE PRECISION FAC,U1,U2
43 REAL RD(3000)
44 INTEGER ISEED(103,10),NSEQ
45 LOGICAL KNOR
46*KEEP,RUNPAR.
47 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
48 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
49 * MONIOU,MDEBUG,NUCNUC,
50 * CETAPE,
51 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
52 * N1STTR,MDBASE,
53 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
54 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
55 * ,GHEISH,GHESIG
56 COMMON /RUNPAC/ DSN,HOST,USER
57 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
58 REAL STEPFC
59 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
60 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
61 * N1STTR,MDBASE
62 INTEGER CETAPE
63 CHARACTER*79 DSN
64 CHARACTER*20 HOST,USER
65
66 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
67 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
68 * ,GHEISH,GHESIG
69*KEND.
70
71 DOUBLE PRECISION AUX2,COSDE,COSTH3,COS3CM,COS3C1,COS3C2,
72 * E3CM,GAMMA3,PHI3CM,PHI3C2,PHI31,
73 * P3CM,XI
74 INTEGER I
75C-----------------------------------------------------------------------
76
77 IF ( DEBUG ) WRITE(MDEBUG,444) (CURPAR(I),I=1,9)
78 444 FORMAT(' MUDECY: CURPAR=',1P,9E10.3)
79
80C MUON DECAYS INTO ELECTRON AND NEUTRINOS
81 XI = 2*ITYPE - 11
82C ELECTRON ENERGY SPECTRUM N(E) * DE = CONST * E**2 * (3/2*E0-E) * DE
83C IS GAINED BY THE REJECTION/REFLECTION METHOD
84 6 CALL RMMAR( RD,4,1 )
85 IF ( RD(1)**2*(3.-RD(1)*2.) .LT. RD(2) ) RD(1) = 1.-RD(1)
86 E3CM = PAMA(2) + RD(1) * ( C(8) - PAMA(2) )
87 IF ( E3CM .GT. 0.5D0*PAMA(5) ) GOTO 6
88 P3CM = SQRT( E3CM**2 - PAMA(2)**2 )
89C NOW DETERMINE COS3C1 AND PHI31 BY RANDOM SELECTION
90C WITH RESPECT TO THE POLARIZATION DIRECTION OF THE MUON IN THE MU CM
91C GIVEN BY POLART, POLARF
92 COSDE = 2.D0 * RD(4) - 1.D0
93 AUX2 = ( 1. - 2.*RD(1) ) / ( 3. - 2.*RD(1) )
94 IF ( ABS(AUX2) .GT. 1.D-2 ) THEN
95 COS3C1 = XI*(SQRT(1.D0-(2.D0*COSDE-AUX2)*AUX2) - 1.D0) / AUX2
96 ELSE
97 COS3C1 = -XI * COSDE
98 ENDIF
99 PHI31 = RD(3)*PI2
100
101C NOW ADD ELECTRON EMISSION ANGLE COS3C1 TO THE POLARISATION DIRECTION
102C TO GET THE DIRECTION (RELATIVE TO THE CORSIKA COORDINATE SYSTEM)
103 CALL ADDANG( POLART,POLARF, COS3C1,PHI31, COS3C2,PHI3C2 )
104C GET THE ELECTRON DIRECTION RELATIVE TO THE MUON LAB DIRECTION
105 CALL ADDANI( CURPAR(3),CURPAR(4), COS3C2,PHI3C2, COS3CM,PHI3CM )
106C LORENTZ TRANSFORMATION TO THE LAB SYSTEM
107 GAMMA3 = GAMMA * ( E3CM + BETA * P3CM * COS3CM ) / PAMA(2)
108 COSTH3 = MIN( 1.D0, GAMMA * (P3CM * COS3CM + BETA * E3CM) /
109 * (PAMA(2) * SQRT(GAMMA3**2 - 1.D0)) )
110 CALL ADDANG( CURPAR(3),CURPAR(4), COSTH3,PHI3CM,
111 * SECPAR(3),SECPAR(4) )
112 IF ( SECPAR(3) .GE. C(29) ) THEN
113 SECPAR(1) = ITYPE - 3
114 SECPAR(2) = GAMMA3
115 DO 10 I = 5,8
116 SECPAR(I) = CURPAR(I)
117 10 CONTINUE
118 SECPAR( 9) = GEN
119 SECPAR(10) = ALEVEL
120 CALL TSTACK
121 ENDIF
122 POLART = 0.D0
123 POLARF = 0.D0
124
125 RETURN
126 END
Note: See TracBrowser for help on using the repository browser.