[u/mrichter/AliRoot.git] / GEANT321 / gphys / gdecay.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:21:23  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.21  by  S.Giani
*-- Author :
      SUBROUTINE GDECAY
C.
C.    ******************************************************************
C.    *                                                                *
C.    *  Control routine for generation of particle decays.            *
C.    *                                                                *
C.    *  MODE(I)        I'th decay mode of current particle            *
C.    *  BRATIO(I)      branching ratio for I'th decay mode of         *
C.    *                 current particle.                              *
C.    *                                                                *
C.    *    ==>Called by : GHSTOP,GTHADR,GTNEUT,GTMUON                  *
C.    *       Author    G.Patrick *********                            *
C.    *                                                                *
C.    ******************************************************************
C.
#include "geant321/gcbank.inc"
#include "geant321/gctrak.inc"
#include "geant321/gconsp.inc"
#include "geant321/gcking.inc"
#include "geant321/gckine.inc"
#include "geant321/gcunit.inc"
#include "geant321/gcphys.inc"
#if defined(CERNLIB_USRJMP)
#include "geant321/gcjump.inc"
#endif
      DIMENSION BAREA(7)
      DIMENSION BETA(4)
      DIMENSION BRATIO(6)
      DIMENSION MODE(6)
      DIMENSION NTYPE(3)
      DIMENSION PCM(4,3)
      DIMENSION XM(3)
      DIMENSION RNDM(1)
      LOGICAL ROTATE
C.
C.    ------------------------------------------------------------------
C.
C             Search for parent particle in decay list.
C
C
      KCASE   = NAMEC(5)
      NGKINE  = 0
      IF(IDCAY.EQ.2) THEN
         DESTEP = DESTEP+GETOT
         ISTOP  = 2
         GO TO 99
      ENDIF
      DMASS = AMASS
      JPA = LQ(JPART-IPART)
      JDK1 = LQ(JPA-1)
      JDK2 = LQ(JPA-2)
      IF (JDK1.LE.0)                               GO TO 90
      IF (JDK2.LE.0)                               GO TO 90
      DO 5 I=1,6
         BRATIO(I)=Q(JDK1+I)
         MODE(I)=IQ(JDK2+I)
   5  CONTINUE
C
C             Generate branching ratio and select decay mode.
C
      NBR      = 1
      BAREA(1) = 0.
      DO 10 I=2,7
         BRADD    = BRATIO(I-1)
         IF (BRADD.EQ.0.) GO TO 20
         NBR      = NBR+1
         BAREA(I) = BAREA(I-1)+BRADD
  10  CONTINUE
C
  20  CALL GRNDM(RNDM,1)
      BRAND    = 100.*RNDM(1)
      IF (BRAND.GE.BAREA(NBR)) GO TO 99
      ID       = IABS((LOCATF(BAREA,NBR,BRAND)))
C
C             Unpack decay mode.
C
      MXX      = MODE(ID)
      NTYPE(1) = MOD(MXX,100)
      NTYPE(2) = MOD(MXX/100,100)
      JP1 = LQ(JPART-NTYPE(1))
      JP2 = LQ(JPART-NTYPE(2))
      XM(1) = Q(JP1+7)
      XM(2) = Q(JP2+7)
      IF (MXX.LT.10000)THEN
C
C             Two body decay.
C
         NGKINE  = 2
         IF (TLIFE.LT.1.E-15) THEN
            XMTOT = XM(1)+XM(2)
            DO 30 I=1,1000
C--  Create Lorentz distributed energy with FWHM HBAR/TLIFE.
C--  (via integral-transformation of Lorentz-distribution)
C--                 (M.Guckes)
              CALL GRNDM(RNDM,1)
              RMASS = DMASS
     1                + 3.291086E-25/TLIFE * TAN(PI*(RNDM(1)-0.5))
              IF (RMASS.GE.XMTOT) GO TO 40
 30         CONTINUE
            WRITE(CHMAIL,1000) IPART, NTYPE(1), NTYPE(2)
            CALL GMAIL(0,0)
            NGKINE=0
            GO TO 99
 40         DMASS = RMASS
         END IF
         CALL GDECA2(DMASS,XM(1),XM(2),PCM)
      ELSE
C
C             Three body decay.
C
         NTYPE(3) = MXX/10000
         NGKINE  = 3
         JP3 = LQ(JPART-NTYPE(3))
         XM(3) = Q(JP3+7)
         IF (TLIFE.LT.1.E-15) THEN
            XMTOT = XM(1)+XM(2)+XM(3)
            DO 31 I=1,1000
C--  Create Lorentz distributed energy with FWHM HBAR/TLIFE.
C--  (via integral-transformation of Lorentz-distribution)
              CALL GRNDM(RNDM,1)
              RMASS = DMASS
     1                + 3.291086E-25/TLIFE * TAN(PI*(RNDM(1)-0.5))
              IF (RMASS.GE.XMTOT) GO TO 41
 31         CONTINUE
            WRITE(CHMAIL,1000) IPART, NTYPE(1), NTYPE(2), NTYPE(3)
            CALL GMAIL(0,0)
            NGKINE=0
            GO TO 99
 41         DMASS = RMASS
         END IF
         CALL GDECA3(DMASS,XM(1),XM(2),XM(3),PCM)
      ENDIF
C
C             LORENTZ boost into LAB system defined along parent vector
C             followed by rotation back into GEANT system.
C
      P0       = VECT(7)
      E0       = SQRT(P0*P0+DMASS*DMASS)
      BETA(1)  = 0.
      BETA(2)  = 0.
      BETA(3)  = -P0/E0
      BETA(4)  = E0/DMASS
      CALL GFANG(VECT(4),COSTH,SINTH,COSPH,SINPH,ROTATE)
C
      DO 60 K=1,NGKINE
         IF (P0.LE.0.) THEN
            DO 59 I = 1,3
   59       GKIN(I,K) = PCM(I,K)
         ELSE
            CALL GLOREN (BETA, PCM(1,K), GKIN(1,K))
         ENDIF
         IF(ROTATE) CALL GDROT  (GKIN(1,K),COSTH,SINTH,COSPH,SINPH)
         GKIN(4,K)=SQRT(GKIN(1,K)**2+GKIN(2,K)**2+GKIN(3,K)**2+XM(K)**2)
         GKIN(5,K)=NTYPE(K)
         TOFD(K)=0.
         GPOS(1,K) = VECT(1)
         GPOS(2,K) = VECT(2)
         GPOS(3,K) = VECT(3)
  60  CONTINUE
      GO TO 99
C
C             No branching ratio defined. Call user routine
C
#if !defined(CERNLIB_USRJMP)
  90  CALL GUDCAY
#endif
#if defined(CERNLIB_USRJMP)
  90  CALL JUMPT0(JUDCAY)
#endif
C
  99  RETURN
 1000 FORMAT(' ***** GDECAY ERROR : Not enough energy available for ',
     +       'decay of resonance',I3,' to',3I3,'; no decay.')
      END
Commit	Line	Data
fe4da5cc	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(P0P0+DMASSDMASS)
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