source: trunk/MagicSoft/Simulation/Corsika/Mmcs/hmeson.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: 5.7 KB
Line 
1 SUBROUTINE HMESON( E,AMASS,ASMASS )
2
3C-----------------------------------------------------------------------
4C H(EAVY) MESON
5C
6C HANDLES PION INITIATED HEAVY MESON AND ITS DECAY IN UP TO 3 PIONS
7C HEAVY MESON EMITTED FORWARD
8C THIS SUBROUTINE IS CALLED FROM BOX67 AND BOX69
9C ARGUMENTS:
10C E = AVAILABLE ENERGY IN CM
11C AMASS = MASS OF HEAVY MESON
12C ASMASS = MASS TO BE LEFT OVER FOR OTHER PARTICLES
13C-----------------------------------------------------------------------
14
15 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
16*KEEP,BAL.
17 COMMON /BAL/ EBAL
18 DOUBLE PRECISION EBAL(10)
19*KEEP,CONST.
20 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
21 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
22*KEEP,ELASTY.
23 COMMON /ELASTY/ ELAST,IELIS,IELHM,IELNU,IELPI
24 DOUBLE PRECISION ELAST
25 INTEGER IELIS(20),IELHM(20),IELNU(20),IELPI(20)
26*KEEP,MULT.
27 COMMON /MULT/ EKINL,MSMM,MULTMA,MULTOT
28 DOUBLE PRECISION EKINL
29 INTEGER MSMM,MULTMA(37,13),MULTOT(37,13)
30*KEEP,PAM.
31 COMMON /PAM/ PAMA,SIGNUM
32 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
33*KEEP,PARPAR.
34 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
35 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
36 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
37 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
38 INTEGER ITYPE,LEVL
39*KEEP,PARPAE.
40 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
41 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
42 * (CURPAR(4), PHI ), (CURPAR(5), H ),
43 * (CURPAR(6), T ), (CURPAR(7), X ),
44 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
45 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
46 * (CURPAR(12),ECM )
47*KEEP,RANDPA.
48 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
49 DOUBLE PRECISION FAC,U1,U2
50 REAL RD(3000)
51 INTEGER ISEED(103,10),NSEQ
52 LOGICAL KNOR
53*KEEP,RUNPAR.
54 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
55 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
56 * MONIOU,MDEBUG,NUCNUC,
57 * CETAPE,
58 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
59 * N1STTR,MDBASE,
60 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
61 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
62 * ,GHEISH,GHESIG
63 COMMON /RUNPAC/ DSN,HOST,USER
64 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
65 REAL STEPFC
66 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
67 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
68 * N1STTR,MDBASE
69 INTEGER CETAPE
70 CHARACTER*79 DSN
71 CHARACTER*20 HOST,USER
72
73 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
74 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
75 * ,GHEISH,GHESIG
76*KEEP,VKIN.
77 COMMON /VKIN/ BETACM
78 DOUBLE PRECISION BETACM
79*KEND.
80
81C-----------------------------------------------------------------------
82
83 IF ( DEBUG ) WRITE(MDEBUG,*) 'HMESON: E,AMASS,ASMASS=',
84 * SNGL(E),SNGL(AMASS),SNGL(ASMASS)
85
86 IPI = 0
87 EDHM = 0.D0
88 PACC = 0.D0
89 W = 0.6D0
90
91C GAMMA AND BETA OF HEAVY MESON IN CM AND LAB
92C E > AMASS + ASMASS TO KEEP GHMCM > 1.
93 GHMCM = ( E**2+AMASS**2-ASMASS**2 ) / ( 2.D0*E*AMASS )
94 BHMCM = SQRT(GHMCM**2 - 1.D0) / GHMCM
95 GHMLAB = GCM * GHMCM * (1.D0 + BETACM * BHMCM)
96 BHMLAB = SQRT(GHMLAB**2 - 1.D0) / GHMLAB
97
98C DECAY OF HEAVY MESON
99 7 CONTINUE
100 IPI = IPI + 1
101C CHOSE TRANSVERSE MOMENTUM RANDOMLY
102 PTPI = PTRANS(DUMMY)
103C CHOSE LONGITUDINAL MOMENTUM RANDOMLY
104 IF ( IPI .LT. 3 ) THEN
105 P = PCL(C(40),W)
106 ELSE
107 P2 = RESTE**2 - PAMA(8)**2 - PTPI**2
108 P = SQRT(MAX( P2, 0.D0 ))
109 ENDIF
110
111 PTPI = PTRANS(DUMMY)
112 GPIHM = SQRT( P**2 / PAMA(8)**2 + 1.D0 )
113 BPIHM = SQRT( GPIHM**2-1.D0 ) / GPIHM
114 EDHM = EDHM + SQRT( PAMA(8)**2 + P**2 + PTPI**2 )
115 RESTE = AMASS - EDHM
116
117C FOR FIRST 2 PARTICLES CHOSE RANDOMLY WHETHER FORWARD OR BACKWARD
118C FOR 3. PARTICLE DECIDE ACCORDING TO ACCULMULATED P
119 CALL RMMAR( RD,3,1 )
120 IF ( IPI .EQ. 3 ) THEN
121 IF ( PACC .LE. 0.D0 ) THEN
122 RD(1) = 0.
123 ELSE
124 RD(1) = 1.
125 ENDIF
126 ENDIF
127
128 IF ( RD(1) .GT. 0.5 ) THEN
129C BACKWARD PION
130 GPILAB = GHMLAB*GPIHM*(1.D0-BHMLAB*BPIHM)
131 PACC = PACC - P
132 ELSE
133C FORWARD PION
134 GPILAB = GHMLAB*GPIHM*(1.D0+BHMLAB*BPIHM)
135 PACC = PACC + P
136 ENDIF
137C CORRECTIVE ACTION IF GPILAB LESS OR EQUAL TO 1.0
138 GPILAB = MAX( GPILAB, 1.D0 )
139
140C GET NEW DIRECTION
141 PLLAB2 = PAMA(8)**2 *(GPILAB**2 - 1.D0)
142 CTHETA = SQRT( PLLAB2 / (PTPI**2+PLLAB2) )
143 IF ( CTHETA .GE. C(27) ) THEN
144 CALL ADDANG( COSTHE,PHI, CTHETA,RD(2)*PI2, SECPAR(3),SECPAR(4) )
145 IF ( SECPAR(3) .GE. C(29) ) THEN
146 SECPAR(2) = GPILAB
147
148C RANDOM CHARGE ASSIGNMENT FOR PIONS
149 IF ( RD(3) .LE. OB3 ) THEN
150 SECPAR(1) = 7.D0
151 ELSEIF ( RD(3) .LE. TB3 ) THEN
152 SECPAR(1) = 8.D0
153 ELSE
154 SECPAR(1) = 9.D0
155 ENDIF
156
157 DO 4 J = 5,8
158 SECPAR(J) = CURPAR(J)
159 4 CONTINUE
160
161 CALL TSTACK
162 ENDIF
163 ENDIF
164 IF ( IPI .LT. 3 .AND. RESTE .GT. PAMA(8) ) GOTO 7
165
166C STATISTICS ON ENERGY BALANCE, MULTIPLICITY AND ELASTICITY
167 EBAL(4) = EBAL(4) + RESTE
168 MSMM = MSMM + IPI
169
170C INELASTICITY STATISTICS
171 IN = 1.D0 + SECPAR(2) / GAMMA * 20.D0
172 IN = MIN( IN, 20 )
173 IELHM(IN) = IELHM(IN) + 1
174
175 RETURN
176 END
Note: See TracBrowser for help on using the repository browser.