source: trunk/MagicSoft/Simulation/Corsika/Mmcs/vhmeso.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: 7.1 KB
Line 
1 SUBROUTINE VHMESO( E,AMASS,ASMASS )
2
3C-----------------------------------------------------------------------
4C (STRANGE) H(EAVY) MESO(N)
5C
6C HANDLES KAON INITIATED HEAVY MESON AND ITS DECAY INTO 1 KAON AND 2 PI
7C STRANGE HEAVY MESON EMITTED FORWARD
8C THIS SUBROUTINE IS CALLED FROM BOX72 AND BOX74
9C ARGUMENTS:
10C E = AVAILABLE ENERGY IN CM
11C AMASS = MASS OF STRANGE 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,KAONS.
23 COMMON /KAONS/ CKA
24 DOUBLE PRECISION CKA(80)
25*KEEP,MULT.
26 COMMON /MULT/ EKINL,MSMM,MULTMA,MULTOT
27 DOUBLE PRECISION EKINL
28 INTEGER MSMM,MULTMA(37,13),MULTOT(37,13)
29*KEEP,PAM.
30 COMMON /PAM/ PAMA,SIGNUM
31 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
32*KEEP,PARPAR.
33 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
34 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
35 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
36 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
37 INTEGER ITYPE,LEVL
38*KEEP,PARPAE.
39 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
40 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
41 * (CURPAR(4), PHI ), (CURPAR(5), H ),
42 * (CURPAR(6), T ), (CURPAR(7), X ),
43 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
44 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
45 * (CURPAR(12),ECM )
46*KEEP,RANDPA.
47 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
48 DOUBLE PRECISION FAC,U1,U2
49 REAL RD(3000)
50 INTEGER ISEED(103,10),NSEQ
51 LOGICAL KNOR
52*KEEP,RUNPAR.
53 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
54 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
55 * MONIOU,MDEBUG,NUCNUC,
56 * CETAPE,
57 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
58 * N1STTR,MDBASE,
59 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
60 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
61 * ,GHEISH,GHESIG
62 COMMON /RUNPAC/ DSN,HOST,USER
63 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
64 REAL STEPFC
65 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
66 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
67 * N1STTR,MDBASE
68 INTEGER CETAPE
69 CHARACTER*79 DSN
70 CHARACTER*20 HOST,USER
71
72 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
73 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
74 * ,GHEISH,GHESIG
75*KEEP,VKIN.
76 COMMON /VKIN/ BETACM
77 DOUBLE PRECISION BETACM
78*KEND.
79
80C-----------------------------------------------------------------------
81
82 IF ( DEBUG ) WRITE(MDEBUG,201) E,AMASS,ASMASS
83 201 FORMAT(' VHMESO: E,AMASS,ASMASS=',1P,3E10.4)
84
85 EDVH = 0.D0
86 W = 0.6D0
87
88C GAMMA AND BETA OF HEAVY MESON IN CM AND LAB
89C E > AMASS + ASMASS TO KEEP GHMCM > 1.
90 GVHCM = (E**2 + AMASS**2 - ASMASS**2) / (2.D0 * E * AMASS)
91 BVHCM = SQRT(GVHCM**2 - 1.D0) / GVHCM
92 GVHLAB = GCM * GVHCM * (1.D0+BETACM*BVHCM)
93 BVHLAB = SQRT(GVHLAB**2 - 1.D0) / GVHLAB
94
95C DECAY OF HEAVY MESON
96
97C KAON PART
98C CHOSE LONGITUDINAL MOMENTUM RANDOMLY FROM EXPONENTIAL DISTRIBUTION
99 P = PCL(CKA(2),W)
100 PT = PTRANS(DUMMY)
101C CALCULATE REST OF ENERGY FOR OTHER PARTICLES
102 EDVH = SQRT( PAMA(ITYPE)**2 + P**2 + PT**2 )
103 RESTE = AMASS - EDVH
104C GAMMA AND BETA OF KAON
105 GKAVH = SQRT( P**2 / PAMA(ITYPE)**2 + 1.D0 )
106 BKAVH = SQRT(GKAVH**2 - 1.D0) / GKAVH
107C KAON FORWARD OR BACKWARD ?
108 CALL RMMAR( RD,1,1 )
109 IF ( RD(1) .LT. 0.5 ) THEN
110 GKALAB = GVHLAB * GKAVH * (1.D0+BKAVH*BVHLAB)
111 PACC = P
112 ELSE
113 GKALAB = GVHLAB * GKAVH * (1.D0-BKAVH*BVHLAB)
114 PACC = -P
115 ENDIF
116C GET NEW DIRECTION
117 PLLAB2 = PAMA(ITYPE)**2 * (GKALAB**2 - 1.D0)
118 PLLAB2 = MAX( 1.D-6, PLLAB2 )
119 CTHETA = SQRT( PLLAB2 / (PT**2+PLLAB2) )
120 IF ( CTHETA .GE. C(27) ) THEN
121 CALL RMMAR( RD,1,1 )
122 CALL ADDANG( COSTHE,PHI, CTHETA,RD(1)*PI2, SECPAR(3),SECPAR(4) )
123 IF ( SECPAR(3) .GE. C(29) ) THEN
124 SECPAR(2) = GKALAB
125C CHARGE ASSIGNMENT
126 SECPAR(1) = CURPAR(1)
127 DO 8 J = 5,8
128 SECPAR(J) = CURPAR(J)
129 8 CONTINUE
130 CALL TSTACK
131 ENDIF
132 ENDIF
133
134C PION PART
135
136C WHAT DECAY MODE ?
137C NDEC=1 : K,PI0,PI0
138C NDEC=2 : K,PI+,PI-
139C NDEC=3 : K,PI-,PI+
140 CALL RMMAR( RD,2,1 )
141 IF ( RD(1) .LT. 0.5 ) THEN
142 NDEC = 1
143 ELSE
144 IF ( RD(2) .LT. 0.5 ) THEN
145 NDEC = 2
146 ELSE
147 NDEC = 3
148 ENDIF
149 ENDIF
150
151 IPI = 0
152 2 CONTINUE
153 IPI = IPI + 1
154
155C TRANSVERS MOMENTUM
156 PT = PTRANS(DUMMY)
157 IF ( IPI .EQ. 1 ) THEN
158C LONGITUDINAL MOMENTUM IS SELECTED FROM EXPONENTIAL DISTRIBUTION
159 P = PCL(C(42),W)
160C CHARGE ASSIGNMENT
161 IF ( NDEC .EQ. 1 ) THEN
162 SECPAR(1) = 7.D0
163 PMA = PAMA(7)
164 ELSEIF ( NDEC .EQ. 2 ) THEN
165 PMA = PAMA(8)
166 SECPAR(1) = 8.D0
167 ELSE
168 PMA = PAMA(9)
169 SECPAR(1) = 9.D0
170 ENDIF
171 ELSE
172C LONGITUDINAL MOMENTUM AS ENERGY IS LEFT
173 P2 = RESTE**2 - PMA**2 - PT**2
174 P = SQRT(MAX( P2, 0.D0 ))
175 IF ( NDEC .EQ. 1 ) THEN
176 SECPAR(1) = 7.D0
177 ELSEIF ( NDEC .EQ. 2 ) THEN
178 SECPAR(1) = 9.D0
179 ELSE
180 SECPAR(1) = 8.D0
181 ENDIF
182 ENDIF
183C REST OF ENERGY FOR OTHER PARTICLES
184 EDVH = EDVH + SQRT( PMA**2 + P**2 + PT**2 )
185 RESTE = AMASS - EDVH
186
187C GAMMA AND BETA OF PION
188 GPIVH = SQRT(P**2/PMA**2 + 1.D0)
189 BEPIVH = SQRT(GPIVH**2 - 1.D0) / GPIVH
190C FOR FIRST PION CHOSE RANDOMLY WHETHER FORWARD OR BACKWARD
191C FOR SECOND PION DECIDE ACCORDING TO ACCUMULATED P
192 CALL RMMAR( RD,2,1 )
193 IF ( IPI .EQ. 2 ) THEN
194 IF ( PACC .LT. 0.D0 ) THEN
195 RD(1) = 0.
196 ELSE
197 RD(1) = 1.
198 ENDIF
199 ENDIF
200 IF ( RD(1) .LT. 0.5 ) THEN
201C BACKWARD PION
202 GPILAB = GVHLAB * GPIVH * (1.D0-BEPIVH*BVHLAB)
203 PACC = PACC - P
204 ELSE
205C FORWARD PION
206 GPILAB = GVHLAB * GPIVH * (1.D0+BEPIVH*BVHLAB)
207 PACC = PACC + P
208 ENDIF
209
210 PLLAB2 = PMA**2 * (GPILAB**2 - 1.D0)
211 PLLAB2 = MAX( 1.D-6, PLLAB2 )
212 CTHETA = SQRT( PLLAB2 / (PT**2+PLLAB2) )
213 IF ( CTHETA .GE. C(27) ) THEN
214 CALL ADDANG( COSTHE,PHI, CTHETA,RD(2)*PI2, SECPAR(3),SECPAR(4) )
215 IF ( SECPAR(3) .GE. C(29) ) THEN
216 SECPAR(2) = GPILAB
217 DO 11 J = 5,8
218 SECPAR(J) = CURPAR(J)
219 11 CONTINUE
220 CALL TSTACK
221 ENDIF
222 ENDIF
223
224 IF ( IPI .LT. 2 .AND. RESTE .GT. PMA ) GOTO 2
225
226 EBAL(5) = EBAL(5) + RESTE
227 MSMM = MSMM + IPI
228
229 RETURN
230 END
231
Note: See TracBrowser for help on using the repository browser.