[u/mrichter/AliRoot.git] / GEANT321 / fluka / hadrin.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:20:00  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.44  by  S.Giani
*-- Author :
*$ CREATE HADRIN.FOR
*COPY HADRIN
*
*=== hadrin ===========================================================*
*
*
*----------------------------------------------------------------------*
*    hadrin89: slight modifications by A. Ferrari
*----------------------------------------------------------------------*
*
      SUBROUTINE HADRIN (N,PLAB,ELAB,CX,CY,CZ,ITTA)
 
#include "geant321/dblprc.inc"
#include "geant321/dimpar.inc"
#include "geant321/iounit.inc"
#include "geant321/metlsp.inc"
#include "geant321/finlsp.inc"
      PARAMETER ( AMPROT = 0.93827231D+00 )
*
#include "geant321/reac.inc"
#include "geant321/redver.inc"
#include "geant321/split.inc"
*
****  INTEGER*2ICH,IBAR,K1,K2,NZK,NRK,IEII,IKII,NURE
      COMMON / FKGAMR / REDU, AMO, AMM (15)
C
      COMMON / FKABLT / AM   (110), GA   (110), TAU  (110), ICH   (110),
     &                  IBAR (110), K1   (110), K2   (110)
      COMMON / FKCODS / COD1, COD2, COD3, COF1, COF2, COF3, SIF1, SIF2,
     &                  SIF3, ECM1, ECM2, ECM3, PCM1, PCM2, PCM3
      COMMON / FKRUN    / RUNTES, EFTES
*
      SAVE UMODA, ITPRF, NNN
      DIMENSION ITPRF(110)
      REAL RNDM(1)
      DATA NNN/0/
      DATA UMODA/0.D0/
      DATA ITPRF/-1,-1,5*1,-1,-1,1,1,1,-1,-1,-1,-1,6*1,-1,-1,-1,85*1/
C
C-----------------------------
C*** INPUT VARIABLES LIST:
C*** SAMPLING OF HADRON NUCLEON INTERACTION FOR (ABOUT) 0.1 LE PLAB LE 6
C*** GEV/C LABORATORY MOMENTUM REGION
C*** N    - PROJECTILE HADRON INDEX
C*** PLAB - LABORATORY MOMENTUM OF N (GEV/C)
C*** ELAB - LABORATORY ENERGY OF N (GEV)
C*** CX,CY,CZ - DIRECTION COSINES OF N IN THE LABORATORY SYSTEM
C*** ITTA - TARGET NUCLEON INDEX
C*** OUTPUT VARIABLES LIST OF PARTICLE CHARACTERISTICS IN /FINLSP/
C  IR COUNTS THE NUMBER OF PRODUCED PARTICLES
C*** ITR - PARTICLE INDEX, CXR,CYR,CZR - DIRECTION COSINES (LAB. SYST.)
C*** ELR,PLR LAB. ENERGY AND LAB. MOMENTUM OF THE SAMPLED PARTICLE
C*** RESPECT., UNITS (GEV/C AND GEV)
C----------------------------
      LOWP=0
      IF (ITPRF(  N  ).LT.0) GO TO 99999
      WRITE(LUNOUT,99998) N
*     STOP    Commented out A. Fasso' 1989
      IR=0
      RETURN
99998 FORMAT (3(5H ****/),
     *45H FALSE USE OF THE PARTICLE TYPE INDEX, VALUE ,I4,3(5H ****/))
99999 CONTINUE
      IATMPT=0
         IF (ABS(PLAB-5.D0).GE.4.99999D0) THEN
            WRITE(LUNOUT,99996) PLAB
*           STOP     Commented out A. Fasso' 1989
            IR=0
            RETURN
99996       FORMAT (3(5H ****/),64
     *H PROJECTILE HADRON MOMENTUM OUTSIDE OF THE ALLOWED REGION, PLAB=,
     *1E15.5/,3(5H ****/))
      END IF
      UMODAT=N*1.11111D0+ITTA*2.19291D0
      IF (UMODAT.NE.UMODA.OR.(AMPROT-AM(1)).GT.1.D-6)CALL CALUMO(N,ITTA)
      UMODA=UMODAT
 1009 CONTINUE
      IATMPT=0
      LOWP=LOWP+1
 1000 CONTINUE
      IMACH=0
      REDU=2.D0
      IW1=0
      IF (LOWP.GT.20) GO TO 8
      NNN=N
      IF (NNN.EQ.N) GO TO 4322
      RUNTES=0.D0
      EFTES=0.D0
 4322 CONTINUE
      IS=1
      IR=0
      IST=1
      NSTAB=25
      IF ( ITTA .EQ. 1 ) THEN
         IRE=NURE(N,1)
      ELSE
         IRE=NURE(N,2)
      END IF
C
C-----------------------------
C*** IE,AMT,ECM,SI DETERMINATION
C----------------------------
      CALL FKSIGI(IRE,PLAB,N,IE,AMT,AMN,ECM,SI,ITTA)
      IF ( AMPROT - AM(1) .GT. 1.D-6 ) THEN
         IANTH = 1
         SI = 1.D0
      ELSE
         IANTH = -1
      END IF
      ECMMH=ECM
C
C-----------------------------
C    ENERGY INDEX
C  IRE CHARACTERIZES THE REACTION
C  IE IS THE ENERGY INDEX
C----------------------------
      IF (SI.LT.1.D-6) GO TO 8
      IF (N .LE.NSTAB) GO TO 1
      RUNTES=RUNTES+1.D0
      IF (RUNTES.LT.20.D0) WRITE(LUNOUT,602)N
 602   FORMAT(3H N=,I10,30H THE PROEKTILE IS A RESONANCE  )
      IF(IBAR(N).EQ.1) N=8
      IF(IBAR(N).EQ.-1)  N=9
   1  CONTINUE
      IMACH=IMACH+1
      IF (IMACH.GT.10) GO TO 8
      ECM =ECMMH
      AMN2=AMN**2
      AMT2=AMT**2
      ECMN=(ECM**2+AMN2-AMT2)/(2.D0*ECM)
      IF(ECMN.LE.AMN) ECMN=AMN
      PCMN=SQRT(ECMN**2-AMN2)
      GAM=(ELAB+AMT)/ECM
      BGAM=PLAB/ECM
      IF (IANTH.GE.0) ECM=2.1D0
C
C-----------------------------
C*** RANDOM CHOICE OF REACTION CHANNEL
C----------------------------
      IST=0
      CALL GRNDM(RNDM,1)
      VV = RNDM (1)
C
C-----------------------------
C***  PLACE REDUCED VERSION
C----------------------------
      IIEI=IEII(IRE)
      IDWK=IEII(IRE+1)-IIEI
      IIWK=IRII(IRE)
      IIKI=IKII(IRE)
C
C-----------------------------
C***  SHRINKAGE TO THE CONSIDERED ENERGY REGION FOR THE USE OF WEIGHTS
C----------------------------
      HECM=ECM
* The following cards assure that Ecm =< Umax + DUmax for this reaction
* where:
*       Umax  = max cms energy at which data are tabulated
*       DUmax = width of the last interval for the tabulated data
      HUMO=2.D0*UMO(IIEI+IDWK)-UMO(IIEI+IDWK-1)
      IF (HUMO.LT.ECM) ECM=HUMO
C
C-----------------------------
C*** INTERPOLATION PREPARATION
C----------------------------
* Cms energy of the upper limit of the considered energy interval
      ECMO=UMO(IE)
* Cms energy of the lower limit of the considered energy interval
      ECM1=UMO(IE-1)
* Width of the considered interval
      DECM=ECMO-ECM1
* Width from actual value to the upper limit (note that if Ecm > Ecmo
* it can be negative but its | | is always less than decm for the
* above condition on Humo
      DEC=ECMO-ECM
C
C-----------------------------
C*** RANDOM LOOP
C----------------------------
* Ik : index of the exit channel
      IK=0
      VFW=VV
      WKK=0.D+00
      WICOR=0.D+00
  111 CONTINUE
      IK=IK+1
* Save the weight accumulated up to now
      VFWO=WKK
* Get the index for the weight of ikth channel at ieth energy (upper
* limit of the interval in energy)
      IWK=IIWK+(IK-1)*IDWK+IE-IIEI
* Cumulative Weight of channels 1...ik at energy Ie
      WOK=WK(IWK)
* Difference for the cumulative weight of channels 1...ik at Ie and Ie-1
      WDK=WOK-WK(IWK-1)
* This card is not clear at all, it should zeroes the weight difference
* if we are in the first interval (that means take only the weights
* at the upper boundary of the interval)
      IF (PLAB.LT.PLABF(IIEI+2)) WDK=0.D+00
C
C-----------------------------
C*** TESTVARIABLE WICO/WICOR: IF CHANNEL IK HAS THE SAME WEIGHTS LIKE IK
C    GO TO NEXT CHANNEL, BECAUSE WKK((IK))-WKK((IK-1))=0, IK CAN NOT
C    CONTRIBUTE
C----------------------------
      WICO=WOK*1.23459876D0+WDK*1.735218469D0
      IF (WICO.EQ.WICOR) GO TO 111
* This card zeroes the weight difference if we are beyond the last
* interval
      IF (UMO(IIEI+IDWK).LT.HECM) WDK=0.D+00
* Save wico
      WICOR=WICO
C
C-----------------------------
C*** INTERPOLATION IN CHANNEL WEIGHTS
C----------------------------
* Set Eklim to a negative value to flag for Iefun it is a
* cms energy and not a lab momentum: this is the cms threshold
* for the exit channel ik
      EKLIM=-THRESH(IIKI+IK)
* Iefun returns the energy index of upper limit of the interval
* containing the threshold
      IELIM=IEFUN(EKLIM,IRE)
* Compute the difference between the upper limit and the
* threshold
      DELIM=UMO(IELIM)+EKLIM+ANGLSQ
* Dete is twice the difference between the actual energy value and
* the average value between Ecmo and the threshold
      DETE=(ECM-(ECMO-EKLIM)*.5D0)*2.D0
      IF ( DELIM*DELIM-DETE*DETE .GT. 0.D+00 ) THEN
         DECC=DELIM
      ELSE
         DECC=DECM
      END IF
      WKK=WOK-WDK*DEC/(DECC+ANGLSQ)
C
C-----------------------------
C*** RANDOM CHOICE
C----------------------------
C
      IF (VV.GT.WKK) GO TO 111
C
C***IK IS THE REACTION CHANNEL
C----------------------------
      INRK=IKII(IRE)+IK
      ECM=HECM
      I1001 =0
C
 1001 CONTINUE
      IT1=NRK(1,INRK)
      AM1=AMGA(IT1)
      IT2=NRK(2,INRK)
      AM2=AMGA(IT2)
      AMS=AM1+AM2
      I1001=I1001+1
      IF (I1001.GT.50) GO TO 1
C
      IF (IT2*AMS.GT.IT2*ECM) GO TO 1001
      IT11=IT1
      IT22=IT2
      IF (IANTH.GE.0) ECM=ELAB+AMT+0.00000001D0
      AM11=AM1
      AM22=AM2
      IF (IT2.GT.0) GO TO 401
C
C-----------------------------
C  INCLUSION OF DIRECT RESONANCES
C  RANDOM CHOICE OF DECAY CHANNELS OF THE DIRECT RESONANCE  IT1
C------------------------
      KZ1=K1(IT1)
      IST=IST+1
      IECO=0
*   Here was the mistake in the pseudo-masses treatment!!!!!
*     ECO=ECM
      ECO=ECMMH
      GAM=(ELAB+AMT)/ECO
      BGAM=PLAB/ECO
      CXS(1)=CX
      CYS(1)=CY
      CZS(1)=CZ
      GO TO 310
  401 CONTINUE
      CALL GRNDM(RNDM,1)
      IF ( RNDM(1) .LT. 0.5D0) GO TO 902
      IT1=IT22
      IT2=IT11
      AM1=AM22
      AM2=AM11
  902   CONTINUE
C
C-----------------------------
C   THE FIRST PARTICLE IS DEFINED TO BE THE FORWARD GOING ONE AT SMALL T
      IBN=IBAR(N)
      IB1=IBAR(IT1)
      IT11=IT1
      IT22=IT2
      AM11=AM1
      AM22=AM2
      IF(IB1.EQ.IBN) GO TO 901
      IT1=IT22
      IT2=IT11
      AM1=AM22
      AM2=AM11
 901  CONTINUE
C-----------------------------
C***IT1,IT2 ARE THE CREATED PARTICLES
C***MOMENTA AND DIRECTION COSINA IN THE CM - SYSTEM
C------------------------
      CALL TWOPAR(ECM1,ECM2,PCM1,PCM2,COD1,COD2,COF1,COF2,SIF1,SIF2,
     *IT1,IT2,ECM,ECMN,PCMN,N,AM1,AM2)
      IST=IST+1
      ITS(IST)=IT1
      AMM(IST)=AM1
C
C-----------------------------
C***TRANSFORMATION INTO LAB SYSTEM AND ROTATION
C----------------------------
      CALL TRAFO(GAM,BGAM,CX,CY,CZ,COD1,COF1,SIF1,PCM1,ECM1,PLS(IST),
     *CXS(IST),CYS(IST),CZS(IST),ELS(IST))
      IST=IST+1
      ITS(IST)=IT2
      AMM(IST)=AM2
      CALL TRAFO(GAM,BGAM,CX,CY,CZ,COD2,COF2,SIF2,PCM2,ECM2,PLS(IST),CXS
     *(IST),CYS(IST),CZS(IST),ELS(IST))
  200 CONTINUE
C
C-----------------------------
C***TEST   STABLE OR UNSTABLE
C----------------------------
      IF(ITS(IST).GT.NSTAB) GO TO 300
      IR=IR+1
C
C-----------------------------
C***IR IS THE NUMBER OF THE FINAL STABLE PARTICLE
C----------------------------
      IF (REDU.LT.0.D0) GO TO 1009
      ITR(IR)=ITS(IST)
       PLR(IR)=PLS(IST)
      CXR(IR)=CXS(IST)
      CYR(IR)=CYS(IST)
      CZR(IR)=CZS(IST)
      ELR(IR)=ELS(IST)
      IST=IST-1
      IF(IST.GE.1) GO TO 200
         GO TO 500
  300 CONTINUE
C
C  RANDOM CHOICE OF DECAY CHANNELS
C----------------------------
C
      IT=ITS(IST)
      ECO=AMM(IST)
      GAM=ELS(IST)/ECO
      BGAM=PLS(IST)/ECO
      IECO=0
      KZ1=K1(IT)
  310 CONTINUE
      IECO=IECO+1
      CALL GRNDM(RNDM,1)
      VV=RNDM(1)
      IIK=KZ1-1
  301 CONTINUE
      IIK=IIK+1
      IF (VV.GT.WT(IIK)) GO TO 301
C
C  IIK IS THE DECAY CHANNEL
C----------------------------
      IT1=NZK(IIK,1)
      I310=0
 1310 CONTINUE
      I310=I310+1
      AM1=AMGA(IT1)
      IT2=NZK(IIK,2)
      AM2=AMGA(IT2)
      IF (IT2-1.LT.0) GO TO 110
      IT3=NZK(IIK,3)
      AM3=AMGA(IT3)
      AMS=AM1+AM2+AM3
C
C  IF  IIK-KIN.LIM.GT.ACTUAL TOTAL CM-ENERGY, DO AGAIN RANDOM IIK-CHOICE
C----------------------------
      IF (IECO.LE.10) GO TO 1002
      IATMPT=IATMPT+1
* Note: we can go to 8 also for too many iterations
      IF (IATMPT.GT.3) GO TO 8
      GO TO 1000
 1002 CONTINUE
      IF (I310.GT.50) GO TO 310
      IF (AMS.GT.ECO) GO TO 1310
C
C  FOR THE DECAY CHANNEL
C  IT1,IT2, IT3 ARE THE PRODUCED PARTICLES FROM  IT
C----------------------------
      IF (REDU.LT.0.D0) GO TO 1009
      ITWTHC=0
      REDU=2.D0
      IF(IT3.EQ.0) GO TO 400
 4001 CONTINUE
      ITWTH=1
      CALL THREPD(ECO,ECM1,ECM2,ECM3,PCM1,PCM2,PCM3,COD1,COF1,SIF1,
     *COD2,COF2,SIF2,COD3,COF3,SIF3,AM1,AM2,AM3)
      GO TO 411
  400 CONTINUE
      CALL TWOPAD(ECO,ECM1,ECM2,PCM1,PCM2,COD1,COF1,SIF1,COD2,COF2,SIF2,
     *AM1,AM2)
      ITWTH=-1
      IT3=0
  411 CONTINUE
      ITWTHC=ITWTHC+1
      IF (REDU.GT.0.D0) GO TO 110
      REDU=2.D0
      IF (ITWTHC.GT.100) GO TO 1009
      IF (ITWTH) 400,400,4001
  110 CONTINUE
      ITS(IST  )=IT1
      IF (IT2-1.LT.0) GO TO 305
      ITS(IST+1)  =IT2
      ITS(IST+2)=IT3
      RX=CXS(IST)
      RY=CYS(IST)
      RZ=CZS(IST)
      AMM(IST)=AM1
      CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD1,COF1,SIF1,PCM1,ECM1,
     *PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
      IST=IST+1
      AMM(IST)=AM2
      CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD2,COF2,SIF2,PCM2,ECM2,
     *PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
      IF (IT3.LE.0) GO TO 305
      IST=IST+1
      AMM(IST)=AM3
      CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD3,COF3,SIF3,PCM3,ECM3,
     *PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
  305 CONTINUE
      GO TO 200
  500 CONTINUE
  631 CONTINUE
      RETURN
    8 CONTINUE
C
C----------------------------
C
C   ZERO CROSS SECTION CASE
C
C----------------------------
C
         IR=1
         ITR(1)=N
         CXR(1)=CX
         CYR(1)=CY
         CZR(1)=CZ
         ELR(1)=ELAB
         PLR(1)=PLAB
      RETURN
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1995/10/24 10:20:00 cernlib
	6	* Geant
	7	*
	8	*
	9	#include "geant321/pilot.h"
	10	*CMZ : 3.21/02 29/03/94 15.41.44 by S.Giani
	11	*-- Author :
	12	*$ CREATE HADRIN.FOR
	13	*COPY HADRIN
	14	*
	15	=== hadrin ===========================================================
	16	*
	17	*
	18	----------------------------------------------------------------------
	19	* hadrin89: slight modifications by A. Ferrari
	20	----------------------------------------------------------------------
	21	*
	22	SUBROUTINE HADRIN (N,PLAB,ELAB,CX,CY,CZ,ITTA)
	23
	24	#include "geant321/dblprc.inc"
	25	#include "geant321/dimpar.inc"
	26	#include "geant321/iounit.inc"
	27	#include "geant321/metlsp.inc"
	28	#include "geant321/finlsp.inc"
	29	PARAMETER ( AMPROT = 0.93827231D+00 )
	30	*
	31	#include "geant321/reac.inc"
	32	#include "geant321/redver.inc"
	33	#include "geant321/split.inc"
	34	*
	35	**** INTEGER*2ICH,IBAR,K1,K2,NZK,NRK,IEII,IKII,NURE
	36	COMMON / FKGAMR / REDU, AMO, AMM (15)
	37	C
	38	COMMON / FKABLT / AM (110), GA (110), TAU (110), ICH (110),
	39	& IBAR (110), K1 (110), K2 (110)
	40	COMMON / FKCODS / COD1, COD2, COD3, COF1, COF2, COF3, SIF1, SIF2,
	41	& SIF3, ECM1, ECM2, ECM3, PCM1, PCM2, PCM3
	42	COMMON / FKRUN / RUNTES, EFTES
	43	*
	44	SAVE UMODA, ITPRF, NNN
	45	DIMENSION ITPRF(110)
	46	REAL RNDM(1)
	47	DATA NNN/0/
	48	DATA UMODA/0.D0/
	49	DATA ITPRF/-1,-1,51,-1,-1,1,1,1,-1,-1,-1,-1,61,-1,-1,-1,85*1/
	50	C
	51	C-----------------------------
	52	C*** INPUT VARIABLES LIST:
	53	C*** SAMPLING OF HADRON NUCLEON INTERACTION FOR (ABOUT) 0.1 LE PLAB LE 6
	54	C*** GEV/C LABORATORY MOMENTUM REGION
	55	C*** N - PROJECTILE HADRON INDEX
	56	C*** PLAB - LABORATORY MOMENTUM OF N (GEV/C)
	57	C*** ELAB - LABORATORY ENERGY OF N (GEV)
	58	C*** CX,CY,CZ - DIRECTION COSINES OF N IN THE LABORATORY SYSTEM
	59	C*** ITTA - TARGET NUCLEON INDEX
	60	C*** OUTPUT VARIABLES LIST OF PARTICLE CHARACTERISTICS IN /FINLSP/
	61	C IR COUNTS THE NUMBER OF PRODUCED PARTICLES
	62	C*** ITR - PARTICLE INDEX, CXR,CYR,CZR - DIRECTION COSINES (LAB. SYST.)
	63	C*** ELR,PLR LAB. ENERGY AND LAB. MOMENTUM OF THE SAMPLED PARTICLE
	64	C*** RESPECT., UNITS (GEV/C AND GEV)
65	C----------------------------
66	LOWP=0
67	IF (ITPRF( N ).LT.0) GO TO 99999
68	WRITE(LUNOUT,99998) N
69	* STOP Commented out A. Fasso' 1989
70	IR=0
71	RETURN
72	99998 FORMAT (3(5H ****/),
73	45H FALSE USE OF THE PARTICLE TYPE INDEX, VALUE ,I4,3(5H ***/))
74	99999 CONTINUE
75	IATMPT=0
76	IF (ABS(PLAB-5.D0).GE.4.99999D0) THEN
77	WRITE(LUNOUT,99996) PLAB
78	* STOP Commented out A. Fasso' 1989
79	IR=0
80	RETURN
81	99996 FORMAT (3(5H ****/),64
82	*H PROJECTILE HADRON MOMENTUM OUTSIDE OF THE ALLOWED REGION, PLAB=,
83	1E15.5/,3(5H ***/))
84	END IF
85	UMODAT=N1.11111D0+ITTA2.19291D0
86	IF (UMODAT.NE.UMODA.OR.(AMPROT-AM(1)).GT.1.D-6)CALL CALUMO(N,ITTA)
87	UMODA=UMODAT
88	1009 CONTINUE
89	IATMPT=0
90	LOWP=LOWP+1
91	1000 CONTINUE
92	IMACH=0
93	REDU=2.D0
94	IW1=0
95	IF (LOWP.GT.20) GO TO 8
96	NNN=N
97	IF (NNN.EQ.N) GO TO 4322
98	RUNTES=0.D0
99	EFTES=0.D0
100	4322 CONTINUE
101	IS=1
102	IR=0
103	IST=1
104	NSTAB=25
105	IF ( ITTA .EQ. 1 ) THEN
106	IRE=NURE(N,1)
107	ELSE
108	IRE=NURE(N,2)
109	END IF
110	C
111	C-----------------------------
112	C*** IE,AMT,ECM,SI DETERMINATION
113	C----------------------------
114	CALL FKSIGI(IRE,PLAB,N,IE,AMT,AMN,ECM,SI,ITTA)
115	IF ( AMPROT - AM(1) .GT. 1.D-6 ) THEN
116	IANTH = 1
117	SI = 1.D0
118	ELSE
119	IANTH = -1
120	END IF
121	ECMMH=ECM
122	C
123	C-----------------------------
124	C ENERGY INDEX
125	C IRE CHARACTERIZES THE REACTION
126	C IE IS THE ENERGY INDEX
127	C----------------------------
128	IF (SI.LT.1.D-6) GO TO 8
129	IF (N .LE.NSTAB) GO TO 1
130	RUNTES=RUNTES+1.D0
131	IF (RUNTES.LT.20.D0) WRITE(LUNOUT,602)N
132	602 FORMAT(3H N=,I10,30H THE PROEKTILE IS A RESONANCE )
133	IF(IBAR(N).EQ.1) N=8
134	IF(IBAR(N).EQ.-1) N=9
135	1 CONTINUE
136	IMACH=IMACH+1
137	IF (IMACH.GT.10) GO TO 8
138	ECM =ECMMH
139	AMN2=AMN**2
140	AMT2=AMT**2
141	ECMN=(ECM*2+AMN2-AMT2)/(2.D0ECM)
142	IF(ECMN.LE.AMN) ECMN=AMN
143	PCMN=SQRT(ECMN**2-AMN2)
144	GAM=(ELAB+AMT)/ECM
145	BGAM=PLAB/ECM
146	IF (IANTH.GE.0) ECM=2.1D0
147	C
148	C-----------------------------
149	C*** RANDOM CHOICE OF REACTION CHANNEL
150	C----------------------------
151	IST=0
152	CALL GRNDM(RNDM,1)
153	VV = RNDM (1)
154	C
155	C-----------------------------
156	C*** PLACE REDUCED VERSION
157	C----------------------------
158	IIEI=IEII(IRE)
159	IDWK=IEII(IRE+1)-IIEI
160	IIWK=IRII(IRE)
161	IIKI=IKII(IRE)
162	C
163	C-----------------------------
164	C*** SHRINKAGE TO THE CONSIDERED ENERGY REGION FOR THE USE OF WEIGHTS
165	C----------------------------
166	HECM=ECM
167	* The following cards assure that Ecm =< Umax + DUmax for this reaction
168	* where:
169	* Umax = max cms energy at which data are tabulated
170	* DUmax = width of the last interval for the tabulated data
171	HUMO=2.D0*UMO(IIEI+IDWK)-UMO(IIEI+IDWK-1)
172	IF (HUMO.LT.ECM) ECM=HUMO
173	C
174	C-----------------------------
175	C*** INTERPOLATION PREPARATION
176	C----------------------------
177	* Cms energy of the upper limit of the considered energy interval
178	ECMO=UMO(IE)
179	* Cms energy of the lower limit of the considered energy interval
180	ECM1=UMO(IE-1)
181	* Width of the considered interval
182	DECM=ECMO-ECM1
183	* Width from actual value to the upper limit (note that if Ecm > Ecmo
184	* it can be negative but its \| \| is always less than decm for the
185	* above condition on Humo
186	DEC=ECMO-ECM
187	C
188	C-----------------------------
189	C*** RANDOM LOOP
190	C----------------------------
191	* Ik : index of the exit channel
192	IK=0
193	VFW=VV
194	WKK=0.D+00
195	WICOR=0.D+00
196	111 CONTINUE
197	IK=IK+1
198	* Save the weight accumulated up to now
199	VFWO=WKK
200	* Get the index for the weight of ikth channel at ieth energy (upper
201	* limit of the interval in energy)
202	IWK=IIWK+(IK-1)*IDWK+IE-IIEI
203	* Cumulative Weight of channels 1...ik at energy Ie
204	WOK=WK(IWK)
205	* Difference for the cumulative weight of channels 1...ik at Ie and Ie-1
206	WDK=WOK-WK(IWK-1)
207	* This card is not clear at all, it should zeroes the weight difference
208	* if we are in the first interval (that means take only the weights
209	* at the upper boundary of the interval)
210	IF (PLAB.LT.PLABF(IIEI+2)) WDK=0.D+00
211	C
212	C-----------------------------
213	C*** TESTVARIABLE WICO/WICOR: IF CHANNEL IK HAS THE SAME WEIGHTS LIKE IK
214	C GO TO NEXT CHANNEL, BECAUSE WKK((IK))-WKK((IK-1))=0, IK CAN NOT
215	C CONTRIBUTE
216	C----------------------------
217	WICO=WOK1.23459876D0+WDK1.735218469D0
218	IF (WICO.EQ.WICOR) GO TO 111
219	* This card zeroes the weight difference if we are beyond the last
220	* interval
221	IF (UMO(IIEI+IDWK).LT.HECM) WDK=0.D+00
222	* Save wico
223	WICOR=WICO
224	C
225	C-----------------------------
226	C*** INTERPOLATION IN CHANNEL WEIGHTS
227	C----------------------------
228	* Set Eklim to a negative value to flag for Iefun it is a
229	* cms energy and not a lab momentum: this is the cms threshold
230	* for the exit channel ik
231	EKLIM=-THRESH(IIKI+IK)
232	* Iefun returns the energy index of upper limit of the interval
233	* containing the threshold
234	IELIM=IEFUN(EKLIM,IRE)
235	* Compute the difference between the upper limit and the
236	* threshold
237	DELIM=UMO(IELIM)+EKLIM+ANGLSQ
238	* Dete is twice the difference between the actual energy value and
239	* the average value between Ecmo and the threshold
240	DETE=(ECM-(ECMO-EKLIM).5D0)2.D0
241	IF ( DELIMDELIM-DETEDETE .GT. 0.D+00 ) THEN
242	DECC=DELIM
243	ELSE
244	DECC=DECM
245	END IF
246	WKK=WOK-WDK*DEC/(DECC+ANGLSQ)
247	C
248	C-----------------------------
249	C*** RANDOM CHOICE
250	C----------------------------
251	C
252	IF (VV.GT.WKK) GO TO 111
253	C
254	C***IK IS THE REACTION CHANNEL
255	C----------------------------
256	INRK=IKII(IRE)+IK
257	ECM=HECM
258	I1001 =0
259	C
260	1001 CONTINUE
261	IT1=NRK(1,INRK)
262	AM1=AMGA(IT1)
263	IT2=NRK(2,INRK)
264	AM2=AMGA(IT2)
265	AMS=AM1+AM2
266	I1001=I1001+1
267	IF (I1001.GT.50) GO TO 1
268	C
269	IF (IT2AMS.GT.IT2ECM) GO TO 1001
270	IT11=IT1
271	IT22=IT2
272	IF (IANTH.GE.0) ECM=ELAB+AMT+0.00000001D0
273	AM11=AM1
274	AM22=AM2
275	IF (IT2.GT.0) GO TO 401
276	C
277	C-----------------------------
278	C INCLUSION OF DIRECT RESONANCES
279	C RANDOM CHOICE OF DECAY CHANNELS OF THE DIRECT RESONANCE IT1
280	C------------------------
281	KZ1=K1(IT1)
282	IST=IST+1
283	IECO=0
284	* Here was the mistake in the pseudo-masses treatment!!!!!
285	* ECO=ECM
286	ECO=ECMMH
287	GAM=(ELAB+AMT)/ECO
288	BGAM=PLAB/ECO
289	CXS(1)=CX
290	CYS(1)=CY
291	CZS(1)=CZ
292	GO TO 310
293	401 CONTINUE
294	CALL GRNDM(RNDM,1)
295	IF ( RNDM(1) .LT. 0.5D0) GO TO 902
296	IT1=IT22
297	IT2=IT11
298	AM1=AM22
299	AM2=AM11
300	902 CONTINUE
301	C
302	C-----------------------------
303	C THE FIRST PARTICLE IS DEFINED TO BE THE FORWARD GOING ONE AT SMALL T
304	IBN=IBAR(N)
305	IB1=IBAR(IT1)
306	IT11=IT1
307	IT22=IT2
308	AM11=AM1
309	AM22=AM2
310	IF(IB1.EQ.IBN) GO TO 901
311	IT1=IT22
312	IT2=IT11
313	AM1=AM22
314	AM2=AM11
315	901 CONTINUE
316	C-----------------------------
317	C***IT1,IT2 ARE THE CREATED PARTICLES
318	C***MOMENTA AND DIRECTION COSINA IN THE CM - SYSTEM
319	C------------------------
320	CALL TWOPAR(ECM1,ECM2,PCM1,PCM2,COD1,COD2,COF1,COF2,SIF1,SIF2,
321	*IT1,IT2,ECM,ECMN,PCMN,N,AM1,AM2)
322	IST=IST+1
323	ITS(IST)=IT1
324	AMM(IST)=AM1
325	C
326	C-----------------------------
327	C***TRANSFORMATION INTO LAB SYSTEM AND ROTATION
328	C----------------------------
329	CALL TRAFO(GAM,BGAM,CX,CY,CZ,COD1,COF1,SIF1,PCM1,ECM1,PLS(IST),
330	*CXS(IST),CYS(IST),CZS(IST),ELS(IST))
331	IST=IST+1
332	ITS(IST)=IT2
333	AMM(IST)=AM2
334	CALL TRAFO(GAM,BGAM,CX,CY,CZ,COD2,COF2,SIF2,PCM2,ECM2,PLS(IST),CXS
335	*(IST),CYS(IST),CZS(IST),ELS(IST))
336	200 CONTINUE
337	C
338	C-----------------------------
339	C***TEST STABLE OR UNSTABLE
340	C----------------------------
341	IF(ITS(IST).GT.NSTAB) GO TO 300
342	IR=IR+1
343	C
344	C-----------------------------
345	C***IR IS THE NUMBER OF THE FINAL STABLE PARTICLE
346	C----------------------------
347	IF (REDU.LT.0.D0) GO TO 1009
348	ITR(IR)=ITS(IST)
349	PLR(IR)=PLS(IST)
350	CXR(IR)=CXS(IST)
351	CYR(IR)=CYS(IST)
352	CZR(IR)=CZS(IST)
353	ELR(IR)=ELS(IST)
354	IST=IST-1
355	IF(IST.GE.1) GO TO 200
356	GO TO 500
357	300 CONTINUE
358	C
359	C RANDOM CHOICE OF DECAY CHANNELS
360	C----------------------------
361	C
362	IT=ITS(IST)
363	ECO=AMM(IST)
364	GAM=ELS(IST)/ECO
365	BGAM=PLS(IST)/ECO
366	IECO=0
367	KZ1=K1(IT)
368	310 CONTINUE
369	IECO=IECO+1
370	CALL GRNDM(RNDM,1)
371	VV=RNDM(1)
372	IIK=KZ1-1
373	301 CONTINUE
374	IIK=IIK+1
375	IF (VV.GT.WT(IIK)) GO TO 301
376	C
377	C IIK IS THE DECAY CHANNEL
378	C----------------------------
379	IT1=NZK(IIK,1)
380	I310=0
381	1310 CONTINUE
382	I310=I310+1
383	AM1=AMGA(IT1)
384	IT2=NZK(IIK,2)
385	AM2=AMGA(IT2)
386	IF (IT2-1.LT.0) GO TO 110
387	IT3=NZK(IIK,3)
388	AM3=AMGA(IT3)
389	AMS=AM1+AM2+AM3
390	C
391	C IF IIK-KIN.LIM.GT.ACTUAL TOTAL CM-ENERGY, DO AGAIN RANDOM IIK-CHOICE
392	C----------------------------
393	IF (IECO.LE.10) GO TO 1002
394	IATMPT=IATMPT+1
395	* Note: we can go to 8 also for too many iterations
396	IF (IATMPT.GT.3) GO TO 8
397	GO TO 1000
398	1002 CONTINUE
399	IF (I310.GT.50) GO TO 310
400	IF (AMS.GT.ECO) GO TO 1310
401	C
402	C FOR THE DECAY CHANNEL
403	C IT1,IT2, IT3 ARE THE PRODUCED PARTICLES FROM IT
404	C----------------------------
405	IF (REDU.LT.0.D0) GO TO 1009
406	ITWTHC=0
407	REDU=2.D0
408	IF(IT3.EQ.0) GO TO 400
409	4001 CONTINUE
410	ITWTH=1
411	CALL THREPD(ECO,ECM1,ECM2,ECM3,PCM1,PCM2,PCM3,COD1,COF1,SIF1,
412	*COD2,COF2,SIF2,COD3,COF3,SIF3,AM1,AM2,AM3)
413	GO TO 411
414	400 CONTINUE
415	CALL TWOPAD(ECO,ECM1,ECM2,PCM1,PCM2,COD1,COF1,SIF1,COD2,COF2,SIF2,
416	*AM1,AM2)
417	ITWTH=-1
418	IT3=0
419	411 CONTINUE
420	ITWTHC=ITWTHC+1
421	IF (REDU.GT.0.D0) GO TO 110
422	REDU=2.D0
423	IF (ITWTHC.GT.100) GO TO 1009
424	IF (ITWTH) 400,400,4001
425	110 CONTINUE
426	ITS(IST )=IT1
427	IF (IT2-1.LT.0) GO TO 305
428	ITS(IST+1) =IT2
429	ITS(IST+2)=IT3
430	RX=CXS(IST)
431	RY=CYS(IST)
432	RZ=CZS(IST)
433	AMM(IST)=AM1
434	CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD1,COF1,SIF1,PCM1,ECM1,
435	*PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
436	IST=IST+1
437	AMM(IST)=AM2
438	CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD2,COF2,SIF2,PCM2,ECM2,
439	*PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
440	IF (IT3.LE.0) GO TO 305
441	IST=IST+1
442	AMM(IST)=AM3
443	CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD3,COF3,SIF3,PCM3,ECM3,
444	*PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
445	305 CONTINUE
446	GO TO 200
447	500 CONTINUE
448	631 CONTINUE
449	RETURN
450	8 CONTINUE
451	C
452	C----------------------------
453	C
454	C ZERO CROSS SECTION CASE
455	C
456	C----------------------------
457	C
458	IR=1
459	ITR(1)=N
460	CXR(1)=CX
461	CYR(1)=CY
462	CZR(1)=CZ
463	ELR(1)=ELAB
464	PLR(1)=PLAB
465	RETURN
466	END