source: trunk/MagicSoft/Simulation/Corsika/Mmcs/decay1.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 DECAY1( M0,M3,M4 )
2
3C-----------------------------------------------------------------------
4C DECAY (INTO TWO PARTICLES)
5C
6C TWO PARTICLE DECAY WITH FULL KINEMATIC; ENERGY AND MOMENTA CONSERVED
7C THIS SUBROUTINE IS CALLED FROM KDECAY, RESDEC, AND STRDEC
8C ARGUMENTS:
9C M0 = TYPE OF DECAYING PARTICLE
10C M3 = TYPE OF FIRST PRODUCT PARTICLE
11C M4 = TYPE OF SECOND PRODUCT PARTICLE
12C
13C DESIGN : D. HECK IK3 FZK KARLSRUHE
14C-----------------------------------------------------------------------
15
16 IMPLICIT NONE
17*KEEP,CONST.
18 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
19 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
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,RANDPA.
38 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
39 DOUBLE PRECISION FAC,U1,U2
40 REAL RD(3000)
41 INTEGER ISEED(103,10),NSEQ
42 LOGICAL KNOR
43*KEEP,RUNPAR.
44 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
45 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
46 * MONIOU,MDEBUG,NUCNUC,
47 * CETAPE,
48 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
49 * N1STTR,MDBASE,
50 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
51 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
52 * ,GHEISH,GHESIG
53 COMMON /RUNPAC/ DSN,HOST,USER
54 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
55 REAL STEPFC
56 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
57 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
58 * N1STTR,MDBASE
59 INTEGER CETAPE
60 CHARACTER*79 DSN
61 CHARACTER*20 HOST,USER
62
63 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
64 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
65 * ,GHEISH,GHESIG
66*KEND.
67
68 DOUBLE PRECISION AUX1,AUX2,AUX2A,AUX3,AUX4,COSTCM,COSTH3,COSTH4,
69 * GAMMA3,GAMMA4,PHI4,WORK1,WORK2
70 INTEGER I,M0,M3,M4
71C-----------------------------------------------------------------------
72
73 IF ( DEBUG ) WRITE(MDEBUG,444) BETA,M0,M3,M4
74 444 FORMAT(' DECAY1: BETA,M0,M3,M4=',1P,E10.3,3I5)
75
76C PARTICLE COORDINATES 5..10 ARE COPIED INTO SECPAR IN CALLING PROGRAM
77C CALCULATE AUXILIARY QUANTITIES
78 AUX1 = ( ( PAMA(M0)**2 + PAMA(M3)**2 - PAMA(M4)**2 )
79 * / (2.D0*PAMA(M0)) )**2 - PAMA(M3)**2
80 AUX2 = 1.D0 + AUX1 / PAMA(M3)**2
81 AUX2A = SQRT(AUX2)
82 AUX3 = SQRT( 1.D0 - 1.D0 / AUX2 )
83
84 WORK1 = GAMMA * AUX2A
85 WORK2 = AUX3 * BETA * WORK1
86
87C DETERMINE POLAR ANGLE IN CM SYSTEM
88 CALL RMMAR( RD,2,1 )
89 COSTCM = 2.D0 * RD(1) - 1.D0
90 GAMMA3 = WORK1 + WORK2 * COSTCM
91C SECOND PRODUCT PARTICLE WITH NONVANISHING REST MASS
92 IF ( PAMA(M4) .NE. 0.D0 ) THEN
93 GAMMA4 = (PAMA(M0)*GAMMA - PAMA(M3)*GAMMA3) / PAMA(M4)
94 AUX4 = (PAMA(M0)**2 + PAMA(M4)**2 - PAMA(M3)**2 )
95 * / (2.D0*PAMA(M0)*PAMA(M4))
96 COSTH4 = MIN( 1.D0, (GAMMA*GAMMA4 - AUX4)
97 * / (BETA * GAMMA * SQRT(GAMMA4**2 - 1.D0)) )
98 ELSE
99C SECOND PRODUCT PARTICLE IS GAMMA; THEN GAMMA4 IS THE ENERGY
100 GAMMA4 = PAMA(M0)*GAMMA - PAMA(M3)*GAMMA3
101 COSTH4 = MIN( 1.D0, (BETA - COSTCM)/(1.D0 - BETA*COSTCM) )
102 ENDIF
103 PHI4 = RD(2)*PI2
104 CALL ADDANG( COSTHE,PHI, COSTH4,PHI4, SECPAR(3),SECPAR(4) )
105 IF ( SECPAR(3) .GT. C(29) ) THEN
106 SECPAR(1) = M4
107 SECPAR(2) = GAMMA4
108 IF ( DEBUG ) WRITE(MDEBUG,445) (SECPAR(I),I=1,9)
109 445 FORMAT(' DECAY1: SECPAR=',1P,9E10.3)
110 CALL TSTACK
111 ENDIF
112C FIRST PRODUCT PARTICLE
113 COSTH3 = MIN( 1.D0, (GAMMA * GAMMA3 - AUX2A)
114 * / (BETA * GAMMA * SQRT(GAMMA3**2 - 1.D0)) )
115 CALL ADDANG( COSTHE,PHI, COSTH3,PHI4+PI, SECPAR(3),SECPAR(4) )
116 IF ( SECPAR(3) .GT. C(29) ) THEN
117 SECPAR(1) = M3
118 SECPAR(2) = GAMMA3
119 IF ( DEBUG ) WRITE(MDEBUG,445) (SECPAR(I),I=1,9)
120 CALL TSTACK
121 ENDIF
122
123 RETURN
124 END
Note: See TracBrowser for help on using the repository browser.