source: trunk/MagicSoft/Simulation/Corsika/Mmcs/decay6.f@ 6525

Last change on this file since 6525 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: 13.3 KB
Line 
1 SUBROUTINE DECAY6(AM0,AM3,AM4,AM5,PARAMA,PARAMB,PARAMC,AMPMX,MODE)
2
3C-----------------------------------------------------------------------
4C DECAY (INTO 3 PARTICLES)
5C
6C TREATES DECAY INTO 3 PARTICLES; FULLY CONSERVING ENERGY AND MOMENTA
7C KINEMATIC RANGE PARAMETRISATION SEE PHYS. LETT. 204B (1988) 90-91
8C FOR LEPTONIC KAON DACAY: THE POLARIZATION OF THE MUON AND
9C THE NEUTRINO PRODUCTION IS INCLUDED.
10C THIS SUBROUTINE IS CALLED FROM ETADEC, KDECAY, AND PI0DEC
11C ARGUMENTS:
12C AM0 = MASS OF DECAYING PARTICLE
13C AM3, AM4, AM5 = MASSES OF RESULTING PARTICLES
14C PARAMA = DALITZ AMPLITUDE PARAMETER (SEE BELOW)
15C PARAMB = DALITZ AMPLITUDE PARAMETER (SEE BELOW)
16C PARAMC = DALITZ AMPLITUDE PARAMETER (SEE BELOW)
17C AMPMX = MAXIMUM AMPLITUDE OF DALITZ PLOT
18C MODE = 1 FOR DECAY KAON ----> 3 PIONS
19C = 2 FOR DECAY ETA ----> 3 PIONS OR 2 PIONS + GAMMA
20C FOR DECAY PI(0) ----> ELECTRON + POSITRON + GAMMA
21C = 3 FOR DECAY KAON ----> PION + MUON + NEUTRINO
22C = 4 FOR DECAY KAON ----> PION + ELECTRON + NEUTRINO
23C
24C AMPLITUDE PARAMETERS PARAMA, PARAMB, PARAMC ARE DEPENDENT ON MODE:
25C FOR MODE=1: PARAMA = G DALITZ AMPLITUDE PARAMETRISATION SEE
26C PARAMB = H PHYS. LETT. 204B (1988) 181 - 193
27C PARAMC = K
28C
29C FOR MODE=2: PARAMA = A DALITZ AMPLITUDE PARAMETRISATION SEE
30C PARAMB = DUMMY PHYS. LETT. 204B (1988) 173 - 175;
31C PARAMC = DUMMY J.G. LAYTER ET.AL. PHYS.REV.D7(1973)2565
32C
33C FOR MODE>2: PARAMA = LAMBDA-PLUS DALITZ AMPLITUDE PARAMETRISATION
34C PARAMB = LAMBDA-ZERO SEE PHYS. LETT. 204B (1988)
35C PARAMC = DUMMY 182 - 194
36C
37C DESIGN : D. HECK IK3 FZK KARLSRUHE
38C-----------------------------------------------------------------------
39
40 IMPLICIT NONE
41*KEEP,CONST.
42 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
43 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
44*KEEP,DECAY.
45 COMMON /DECAY/ GAM345,COS345,PHI345
46 DOUBLE PRECISION GAM345(3),COS345(3),PHI345(3)
47*KEEP,PAM.
48 COMMON /PAM/ PAMA,SIGNUM
49 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
50*KEEP,PARPAR.
51 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
52 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
53 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
54 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
55 INTEGER ITYPE,LEVL
56*KEEP,PARPAE.
57 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
58 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
59 * (CURPAR(4), PHI ), (CURPAR(5), H ),
60 * (CURPAR(6), T ), (CURPAR(7), X ),
61 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
62 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
63 * (CURPAR(12),ECM )
64*KEEP,POLAR.
65 COMMON /POLAR/ POLART,POLARF
66 DOUBLE PRECISION POLART,POLARF
67*KEEP,RANDPA.
68 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
69 DOUBLE PRECISION FAC,U1,U2
70 REAL RD(3000)
71 INTEGER ISEED(103,10),NSEQ
72 LOGICAL KNOR
73*KEEP,RUNPAR.
74 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
75 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
76 * MONIOU,MDEBUG,NUCNUC,
77 * CETAPE,
78 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
79 * N1STTR,MDBASE,
80 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
81 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
82 * ,GHEISH,GHESIG
83 COMMON /RUNPAC/ DSN,HOST,USER
84 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
85 REAL STEPFC
86 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
87 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
88 * N1STTR,MDBASE
89 INTEGER CETAPE
90 CHARACTER*79 DSN
91 CHARACTER*20 HOST,USER
92
93 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
94 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
95 * ,GHEISH,GHESIG
96*KEND.
97
98 DOUBLE PRECISION ABYM,AMPLI,AMPMX,AM0,AM3,AM34I,AM34SQ,AM35SQ,
99 * AM4,AM5,APARAL,APERPN,AUXA,AUXB,AUX1,AUX2,AUX2A,
100 * AUX3,AUX4,AUX4A,AUX5,AUX6,
101 * AUX7,AUX8,AUX10,AUX12,AUX14,BBYM,BOFQ,
102 * CM0SQ,CM3SQ,CM3SQI,CM4SQ,CM5SQ,COSALF,COSBET,
103 * COSFI4,COSOME,COSPHI,COSPSI,COS3CM,COS4CM,COS5CM,
104 * DISCR,EPIPRM,E3CM,E3STAR,E4CM,E5CM,E5STAR,FACT,
105 * GRLAMD,OMEGA,PA,PARAMA,PARAMB,PARAMC,PB,PC,PSI,
106 * P3CM,P3SQ,P4CM,P4SQ,P5CM,P5SQ,ROOT1,ROOT2,
107 * SINALF,SINBET,SINFI4,SINFI5,SINOMG,SINPHI,SINPSI,
108 * SINT4,SINT4I,SINT5I,SIN3CM,S0,TBYMSS,XIT,XI0
109 INTEGER MODE
110C-----------------------------------------------------------------------
111
112 IF ( DEBUG ) WRITE(MDEBUG,444) AM0,AM3,AM4,AM5
113 444 FORMAT(' DECAY6: AM0',1P,E10.3,' AM3',E10.3,' AM4',E10.3,
114 * ' AM5',E10.3)
115
116C CALCULATE AUXILIARY QUANTITIES
117 CM0SQ = AM0**2
118 CM3SQ = AM3**2
119 CM4SQ = AM4**2
120 CM5SQ = AM5**2
121 AUX1 = (AM3 + AM4)**2
122 AUX2A = (AM0 - AM5)**2
123 AUX2 = AUX2A - AUX1
124 AUX3 = (AM3 + AM5)**2
125 AUX4A = (AM0 - AM4)**2
126 AUX4 = AUX4A - AUX3
127 AUX5 = CM3SQ - CM4SQ
128 AUX6 = CM0SQ - CM5SQ
129 AUX7 = 0.5D0 / AM0
130 IF ( MODE .EQ. 1 ) THEN
131 AUX8 = (AM0 - AM3)**2
132 S0 = OB3 * ( CM0SQ + CM3SQ + CM4SQ + CM5SQ )
133 AUX10 = 1.D0 / PAMA(8)**2
134 ELSEIF ( MODE .EQ. 2 ) THEN
135 AUX14 = 1.D0 / (AM0 - AM3 - AM4 - AM5)
136 ELSEIF ( MODE .EQ. 3 .OR. MODE .EQ. 4 ) THEN
137 CM3SQI = 1.D0 / CM3SQ
138 AUX12 = (CM0SQ + CM3SQ - CM4SQ) * AUX7
139C XI0 IS XI(0); GRLAMD IS GREAT LAMBDA
140 XI0 = ( CM0SQ - CM3SQ) * CM3SQI * (PARAMB - PARAMA)
141 GRLAMD = -XI0 * PARAMA
142 ELSE
143 WRITE(MONIOU,*) 'DECAY6: UNEXPECTED MODE =',MODE
144 RETURN
145 ENDIF
146
147 100 CALL RMMAR( RD,3,1 )
148C ARE INVARIANT MASS SQUARES INSIDE BOUNDARY OF DALITZ PLOT?
149 AM34SQ = AUX2 * RD(1) + AUX1
150 AM35SQ = AUX4 * RD(2) + AUX3
151 AM34I = 0.5D0 / SQRT(AM34SQ)
152 E3STAR = (AUX5 + AM34SQ) * AM34I
153 E5STAR = (AUX6 - AM34SQ) * AM34I
154 ROOT1 = SQRT(E3STAR**2 - CM3SQ )
155 ROOT2 = SQRT(E5STAR**2 - CM5SQ )
156 DISCR = AM35SQ - (E3STAR + E5STAR)**2
157C REJECT RANDOM NUMBERS, IF OUTSIDE KINEMATIC BOUNDARY OF DALITZ PLOT
158 IF ( DISCR .GT. -(ROOT1 - ROOT2)**2 ) GOTO 100
159 IF ( DISCR .LT. -(ROOT1 + ROOT2)**2 ) GOTO 100
160C E3CM, E4CM, E5CM ARE ENERGIES IN THE C. M. SYSTEM
161 E4CM = (CM0SQ + CM4SQ - AM35SQ) * AUX7
162 E5CM = (CM0SQ + CM5SQ - AM34SQ) * AUX7
163 E3CM = AM0 - E4CM - E5CM
164
165 IF ( MODE .EQ. 1 ) THEN
166 FACT = AUX10 * (AUX2A - 2.D0*AM0*(E5CM-AM5) - S0)
167C AMPLITUDE OF SQUARED MATRIX ELEMENT (SEE PHYS. LETT. B204 (1988) 181)
168 AMPLI = 1.D0 + PARAMA*FACT + PARAMB*FACT**2 + PARAMC*( AUX10
169 * * ( AUX4A -AUX8 -2.D0*(E4CM-AM4-E3CM+AM3)*AM0 ) )**2
170
171 ELSEIF ( MODE .EQ. 2 ) THEN
172C AMPLITUDE OF SQUARED MATRIX ELEMENT (SEE PHYS. LETT. B204 (1988) 173)
173C REF: J. G. LAYTER ET AL., PHYS. REV. D7 (1973) 2565
174 AMPLI = 1.D0 + PARAMA * ( 3.D0 * (E5CM - AM5) * AUX14 - 1.D0 )
175
176 ELSE
177C EPIPRM IS (ENERGY OF PION)PRIMED
178 EPIPRM = AUX12 - E3CM
179C PA, PB, AND PC ARE THE A, B, AND C PARAMETERS
180 PA = AM0 * ( 2.D0 * E4CM * E5CM - AM0 * EPIPRM )
181 * + CM4SQ * ( 0.25D0 * EPIPRM - E5CM )
182 PB = CM4SQ * ( E5CM - 0.5D0 * EPIPRM )
183 PC = CM4SQ * EPIPRM * 0.25D0
184C TBYMSS IS T DIVIDED BY MASS SQUARE OF PION
185 TBYMSS = (CM0SQ + CM3SQ - 2.D0 * AM0 * E3CM) * CM3SQI
186C XIT IS XI(T)
187 XIT = XI0 + GRLAMD*TBYMSS
188C AMPLITUDE OF SQUARED MATRIX ELEMENT (PHYS. LETT. B204 (1988) 183)
189 AMPLI = (1.D0 + PARAMA*TBYMSS)**2 * ( PA + XIT*PB + XIT**2 *PC )
190 ENDIF
191
192C REJECT RANDOM NUMBERS, IF RD(3) IS LARGER THAN DALITZ PLOT AMPLITUDE
193 IF ( RD(3)*AMPMX .GT. AMPLI ) GOTO 100
194
195 IF (DEBUG) WRITE(MDEBUG,*)'DECAY6: E3CM,E4CM,E5CM=',
196 * SNGL(E3CM),SNGL(E4CM),SNGL(E5CM)
197C P3CM, P4CM, P5CM ARE MOMENTA IN THE C.M. SYSTEM
198C P3SQ, P4SQ, P5SQ ARE SQUARED MOMENTA IN C.M. SYSTEM
199 P5SQ = E5CM**2 - CM5SQ
200 P5CM = SQRT(P5SQ)
201 P4SQ = E4CM**2 - CM4SQ
202 P4CM = SQRT(P4SQ)
203 P3SQ = E3CM**2 - CM3SQ
204 P3CM = SQRT(P3SQ)
205C ANGLE ALFA AND BETA ARE BETWEEN PARTICLE 3 AND 4 RSP. 3 AND 5
206 COSALF = (P5SQ - P3SQ - P4SQ) / (2.D0 * P3CM * P4CM)
207 SINALF = -SQRT(1.D0 - COSALF**2)
208 COSBET = (P4SQ - P3SQ - P5SQ) / (2.D0 * P3CM * P5CM)
209 SINBET = SQRT(1.D0 - COSBET**2)
210C NOW SELECT RANDOM NUMBERS FOR THREE INDEPENDENT ANGLES IN CM-SYSTEM
211C COS3CM AND PHI ARE ANGLES OF PARTICLE 3 RELATIVE TO DECAYING PARTICLE
212 CALL RMMAR( RD,3,1 )
213 COS3CM = 2.D0*RD(1) - 1.D0
214 SIN3CM = SQRT(1.D0 - COS3CM**2)
215 PHI345(1) = PI2 * RD(2)
216 COSPHI = COS( PHI345(1) )
217 SINPHI = SIN( PHI345(1) )
218C ANGLE PSI GIVES ROTATION OF PLANE (3,4,5) RELATIVE TO PLANE (1,3)
219 PSI = PI2 * RD(3)
220 COSPSI = COS(PSI)
221 SINPSI = SIN(PSI)
222C CALCULATE ALL NEEDED POLAR AND AZIMUTHAL ANGLES IN THE CM-SYSTEM
223 COS4CM = COS3CM * COSALF - SIN3CM * COSPSI * SINALF
224 IF ( ABS(COS4CM) .LT. 1.D0 ) THEN
225 SINT4 = SQRT(1.D0 - COS4CM**2)
226 SINT4I = 1.D0 / SINT4
227 AUXA = COS3CM * COSPSI * SINALF + SIN3CM * COSALF
228 COSFI4 = (COSPHI*AUXA-SINPHI*SINPSI*SINALF) * SINT4I
229 PHI345(2) = ACOS( MAX( -1.D0, MIN( 1.D0, COSFI4 ) ) )
230 SINFI4 = (SINPHI*AUXA+COSPHI*SINPSI*SINALF) * SINT4I
231 IF ( SINFI4 .LE. 0.D0 ) PHI345(2) = PI2 - PHI345(2)
232 ELSE
233 PHI345(2) = 0.D0
234 ENDIF
235C CALCULATE GAMMA FACTORS AND POLAR ANGLES IN LABORATORY SYSTEM
236 GAM345(1) = GAMMA * (E3CM + BETA * P3CM * COS3CM) / AM3
237 COS345(1) = MIN( 1.D0, (BETA * E3CM + P3CM * COS3CM) * GAMMA
238 * / ( AM3 * SQRT(GAM345(1)**2 - 1.D0) ) )
239 GAM345(2) = GAMMA * (E4CM + BETA * P4CM * COS4CM) / AM4
240 COS345(2) = MIN( 1.D0, (BETA * E4CM + P4CM * COS4CM) * GAMMA
241 * / ( AM4 * SQRT(GAM345(2)**2 - 1.D0) ) )
242C CALCULATE PARAMETERS OF PARTICLE 5, IF NEEDED
243 IF ( MODE .LE. 2 ) THEN
244 COS5CM = COS3CM * COSBET - SIN3CM * COSPSI * SINBET
245 IF ( ABS(COS5CM) .LT. 1.D0 ) THEN
246 SINT5I = 1.D0 / SQRT(1.D0 - COS5CM**2)
247 AUXB = COS3CM * COSPSI * SINBET + SIN3CM * COSBET
248 PHI345(3) = ACOS((COSPHI*AUXB-SINPHI*SINPSI*SINBET) * SINT5I)
249 SINFI5 = (SINPHI*AUXB+COSPHI*SINPSI*SINBET) * SINT5I
250 IF ( SINFI5 .LE. 0.D0 ) PHI345(3) = PI2 - PHI345(3)
251 ELSE
252 PHI345(3) = 0.D0
253 ENDIF
254 IF ( AM5 .NE. 0.D0 ) THEN
255 GAM345(3) = GAMMA * (E5CM + BETA * P5CM * COS5CM) / AM5
256 COS345(3) = MIN( 1.D0, (BETA * E5CM + P5CM * COS5CM) * GAMMA
257 * / ( AM5 * SQRT(GAM345(3)**2 - 1.D0) ) )
258 ELSE
259C IF PARTICLE 5 IS GAMMA RAY OR NEUTRINO, THEN GAM345(3) IS THE ENERGY
260 GAM345(3) = GAMMA * (E5CM + BETA * P5CM * COS5CM)
261 COS345(3) = MIN( 1.D0, (BETA * E5CM + P5CM * COS5CM) * GAMMA
262 * / GAM345(3) )
263 ENDIF
264 ENDIF
265
266 IF ( MODE .EQ. 3 ) THEN
267C CALCULATION OF MUON POLARIZATION. WE FOLLOW THE DESCRIPTION OF
268C L. JAUNEAU, IN: METHODS IN SUBNUCLEAR PHYSICS, VOL. 3, M. NIKOLIC ED.
269C (GORDON + BREACH, NEW YORK, 1969), P. 123
270C SEE ALSO: L.M. CHOUNET ET AL., PHYS. REP. 4 (1972) 199, APPENDIX 1.
271C SEE ALSO: N. CABBIBO, A. MAKSYMOWICZ, PHYS. LETT. 9 (1964) 352
272C (CORRECTIONS IN: PHYS. LETT. 11 (1964) 360; 14 (1965) 72)
273C WE DEFINE BOFQ (READ: B OF Q), WHICH IS -B(Q**2)*4
274 BOFQ = 1.D0 - XIT
275C ABYM AND BBYM (READ A BY M; B BY M) ARE THE QUANTITIES A/M AND B/M
276 ABYM = AM0 * ( BOFQ * EPIPRM - 2.D0 * E5CM )
277 BBYM = CM0SQ + 0.25D0 * CM4SQ * BOFQ**2 - BOFQ * AM0 * E4CM
278C NOW CALCULATE THE COMPONENTS APARAL (PARALLEL TO MU DIRECTION) AND
279C APERPN (PERPENDICULAR TO MU DIRECTION) USING QUANTITIES DEFINED IN
280C KAON REST SYSTEM. NOTE OUR DEFINITION OF SINALF (ALWAYS WITH NEGATIVE
281C SIGN) OPPOSITE TO CABBIBO'S SIN(PSI) AND JAUNEAU'S SIN(THETA)
282 APARAL = -P3CM*AM4*BBYM*COSALF - P4CM * ( AM0*ABYM - BBYM *
283 * ( P3CM*SINALF*(E4CM-AM4)/P4CM + AM0 - E3CM ) )
284 APERPN = P3CM*AM4*BBYM*SINALF
285C NOW NORMALIZE THE PARALLEL COMPONENT OF POLARIZATION; POLART IS
286C COSINE OF THE ANGLE BETWEEN MUON MOMENTUM AND POLARISATION
287 POLART = APARAL / SQRT(APARAL**2 + APERPN**2)
288C THE POLARIZATION VECTOR LIES IN THE PLANE OF MOMENTA (PION,MUON).
289C OMEGA IS THE ANGLE BY WHICH THE DECAY PLANE (PION,MUON) IS ROTATET
290C AROUND THE DIRECTION OF MUON RELATIVE TO THE PLANE (KAON,MUON)
291 IF ( ABS(COS4CM) .LT. 1.D0 .AND. SINALF .NE. 0.D0 ) THEN
292 COSOME = (COS4CM*COSALF - COS3CM)*SINT4I/SINALF
293 OMEGA = ACOS( MAX( -1.D0, MIN( 1.D0, COSOME ) ) )
294 IF ( SINFI4 .NE. 0.D0 ) THEN
295 SINOMG = ( COSFI4 * ( COSALF - COS3CM*COS4CM ) * SINT4I
296 * - SIN3CM * COSPHI ) / (SINALF*SINFI4)
297 IF ( SINOMG .LT. 0.D0 ) OMEGA = PI2 - OMEGA
298 ENDIF
299 ELSE
300 OMEGA = 0.D0
301 ENDIF
302 POLARF = OMEGA
303 ENDIF
304
305 RETURN
306 END
Note: See TracBrowser for help on using the repository browser.