Bugfix in AliPoints2Memory
[u/mrichter/AliRoot.git] / GEANT321 / gphys / gdecay.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:21:23  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.21  by  S.Giani
11 *-- Author :
12       SUBROUTINE GDECAY
13 C.
14 C.    ******************************************************************
15 C.    *                                                                *
16 C.    *  Control routine for generation of particle decays.            *
17 C.    *                                                                *
18 C.    *  MODE(I)        I'th decay mode of current particle            *
19 C.    *  BRATIO(I)      branching ratio for I'th decay mode of         *
20 C.    *                 current particle.                              *
21 C.    *                                                                *
22 C.    *    ==>Called by : GHSTOP,GTHADR,GTNEUT,GTMUON                  *
23 C.    *       Author    G.Patrick *********                            *
24 C.    *                                                                *
25 C.    ******************************************************************
26 C.
27 #include "geant321/gcbank.inc"
28 #include "geant321/gctrak.inc"
29 #include "geant321/gconsp.inc"
30 #include "geant321/gcking.inc"
31 #include "geant321/gckine.inc"
32 #include "geant321/gcunit.inc"
33 #include "geant321/gcphys.inc"
34 #if defined(CERNLIB_USRJMP)
35 #include "geant321/gcjump.inc"
36 #endif
37       DIMENSION BAREA(7)
38       DIMENSION BETA(4)
39       DIMENSION BRATIO(6)
40       DIMENSION MODE(6)
41       DIMENSION NTYPE(3)
42       DIMENSION PCM(4,3)
43       DIMENSION XM(3)
44       DIMENSION RNDM(1)
45       LOGICAL ROTATE
46 C.
47 C.    ------------------------------------------------------------------
48 C.
49 C             Search for parent particle in decay list.
50 C
51 C
52       KCASE   = NAMEC(5)
53       NGKINE  = 0
54       IF(IDCAY.EQ.2) THEN
55          DESTEP = DESTEP+GETOT
56          ISTOP  = 2
57          GO TO 99
58       ENDIF
59       DMASS = AMASS
60       JPA = LQ(JPART-IPART)
61       JDK1 = LQ(JPA-1)
62       JDK2 = LQ(JPA-2)
63       IF (JDK1.LE.0)                               GO TO 90
64       IF (JDK2.LE.0)                               GO TO 90
65       DO 5 I=1,6
66          BRATIO(I)=Q(JDK1+I)
67          MODE(I)=IQ(JDK2+I)
68    5  CONTINUE
69 C
70 C             Generate branching ratio and select decay mode.
71 C
72       NBR      = 1
73       BAREA(1) = 0.
74       DO 10 I=2,7
75          BRADD    = BRATIO(I-1)
76          IF (BRADD.EQ.0.) GO TO 20
77          NBR      = NBR+1
78          BAREA(I) = BAREA(I-1)+BRADD
79   10  CONTINUE
80 C
81   20  CALL GRNDM(RNDM,1)
82       BRAND    = 100.*RNDM(1)
83       IF (BRAND.GE.BAREA(NBR)) GO TO 99
84       ID       = IABS((LOCATF(BAREA,NBR,BRAND)))
85 C
86 C             Unpack decay mode.
87 C
88       MXX      = MODE(ID)
89       NTYPE(1) = MOD(MXX,100)
90       NTYPE(2) = MOD(MXX/100,100)
91       JP1 = LQ(JPART-NTYPE(1))
92       JP2 = LQ(JPART-NTYPE(2))
93       XM(1) = Q(JP1+7)
94       XM(2) = Q(JP2+7)
95       IF (MXX.LT.10000)THEN
96 C
97 C             Two body decay.
98 C
99          NGKINE  = 2
100          IF (TLIFE.LT.1.E-15) THEN
101             XMTOT = XM(1)+XM(2)
102             DO 30 I=1,1000
103 C--  Create Lorentz distributed energy with FWHM HBAR/TLIFE.
104 C--  (via integral-transformation of Lorentz-distribution)
105 C--                 (M.Guckes)
106               CALL GRNDM(RNDM,1)
107               RMASS = DMASS
108      1                + 3.291086E-25/TLIFE * TAN(PI*(RNDM(1)-0.5))
109               IF (RMASS.GE.XMTOT) GO TO 40
110  30         CONTINUE
111             WRITE(CHMAIL,1000) IPART, NTYPE(1), NTYPE(2)
112             CALL GMAIL(0,0)
113             NGKINE=0
114             GO TO 99
115  40         DMASS = RMASS
116          END IF
117          CALL GDECA2(DMASS,XM(1),XM(2),PCM)
118       ELSE
119 C
120 C             Three body decay.
121 C
122          NTYPE(3) = MXX/10000
123          NGKINE  = 3
124          JP3 = LQ(JPART-NTYPE(3))
125          XM(3) = Q(JP3+7)
126          IF (TLIFE.LT.1.E-15) THEN
127             XMTOT = XM(1)+XM(2)+XM(3)
128             DO 31 I=1,1000
129 C--  Create Lorentz distributed energy with FWHM HBAR/TLIFE.
130 C--  (via integral-transformation of Lorentz-distribution)
131               CALL GRNDM(RNDM,1)
132               RMASS = DMASS
133      1                + 3.291086E-25/TLIFE * TAN(PI*(RNDM(1)-0.5))
134               IF (RMASS.GE.XMTOT) GO TO 41
135  31         CONTINUE
136             WRITE(CHMAIL,1000) IPART, NTYPE(1), NTYPE(2), NTYPE(3)
137             CALL GMAIL(0,0)
138             NGKINE=0
139             GO TO 99
140  41         DMASS = RMASS
141          END IF
142          CALL GDECA3(DMASS,XM(1),XM(2),XM(3),PCM)
143       ENDIF
144 C
145 C             LORENTZ boost into LAB system defined along parent vector
146 C             followed by rotation back into GEANT system.
147 C
148       P0       = VECT(7)
149       E0       = SQRT(P0*P0+DMASS*DMASS)
150       BETA(1)  = 0.
151       BETA(2)  = 0.
152       BETA(3)  = -P0/E0
153       BETA(4)  = E0/DMASS
154       CALL GFANG(VECT(4),COSTH,SINTH,COSPH,SINPH,ROTATE)
155 C
156       DO 60 K=1,NGKINE
157          IF (P0.LE.0.) THEN
158             DO 59 I = 1,3
159    59       GKIN(I,K) = PCM(I,K)
160          ELSE
161             CALL GLOREN (BETA, PCM(1,K), GKIN(1,K))
162          ENDIF
163          IF(ROTATE) CALL GDROT  (GKIN(1,K),COSTH,SINTH,COSPH,SINPH)
164          GKIN(4,K)=SQRT(GKIN(1,K)**2+GKIN(2,K)**2+GKIN(3,K)**2+XM(K)**2)
165          GKIN(5,K)=NTYPE(K)
166          TOFD(K)=0.
167          GPOS(1,K) = VECT(1)
168          GPOS(2,K) = VECT(2)
169          GPOS(3,K) = VECT(3)
170   60  CONTINUE
171       GO TO 99
172 C
173 C             No branching ratio defined. Call user routine
174 C
175 #if !defined(CERNLIB_USRJMP)
176   90  CALL GUDCAY
177 #endif
178 #if defined(CERNLIB_USRJMP)
179   90  CALL JUMPT0(JUDCAY)
180 #endif
181 C
182   99  RETURN
183  1000 FORMAT(' ***** GDECAY ERROR : Not enough energy available for ',
184      +       'decay of resonance',I3,' to',3I3,'; no decay.')
185       END