This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / fluka / hadriv.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:19:57  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.43  by  S.Giani
11 *-- Author :
12 *$ CREATE HADRIV.FOR
13 *COPY HADRIV
14 *
15 *=== hadriv ===========================================================*
16 *
17 *
18 *----------------------------------------------------------------------*
19 *                                                                      *
20 *    Modified version of Hadrin created by Alfredo Ferrari, INFN-Milan *
21 *                                                                      *
22 *    Last change  on  20-jun-92  by  Alfredo Ferrari, INFN - MIlan     *
23 *                                                                      *
24 *    hadriv: this is a modified version of Hadrin, used by the Eventv  *
25 *            package for hadron-hadron interactions below 5 GeV        *
26 *                                                                      *
27 *----------------------------------------------------------------------*
28 *
29       SUBROUTINE HADRIV ( N, PLAB, ELAB, CX, CY, CZ, ITTA )
30  
31 #include "geant321/dblprc.inc"
32 #include "geant321/dimpar.inc"
33 #include "geant321/iounit.inc"
34 *
35 #include "geant321/finlsp.inc"
36 #include "geant321/hadflg.inc"
37 #include "geant321/metlsp.inc"
38 #include "geant321/reac.inc"
39 #include "geant321/redver.inc"
40 #include "geant321/split.inc"
41 *
42       COMMON / FKGAMR / REDU, AMO, AMM (15)
43 C
44       COMMON / FKABLT / AM   (110), GA   (110), TAU  (110), ICH   (110),
45      &                  IBAR (110), K1   (110), K2   (110)
46       COMMON / FKCODS / COD1, COD2, COD3, COF1, COF2, COF3, SIF1, SIF2,
47      &                  SIF3, ECM1, ECM2, ECM3, PCM1, PCM2, PCM3
48       COMMON / FKRUN    / RUNTES, EFTES
49 *
50       PARAMETER ( AMPROT = 0.93827231D+00 )
51       DIMENSION WCHANN (40), WCUMCH (0:40), IKIK (40)
52       DIMENSION ITPRF(110)
53       REAL RNDM(1)
54       LOGICAL LSWAP
55 *
56       SAVE IKIK,ITPRF
57       DATA NNN / 0 /
58       DATA IKIK /  1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15,
59      &            16, 17, 18 ,19, 20, 21, 22, 23, 24, 25, 26, 27, 28,
60      &            29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40 /
61       DATA ITPRF/-1,-1,5*1,-1,-1,1,1,1,-1,-1,-1,-1,6*1,-1,-1,-1,85*1/
62 C
63 C-----------------------------
64 C*** INPUT VARIABLES LIST:
65 C*** SAMPLING OF HADRON NUCLEON INTERACTION FOR (ABOUT) 0.1 LE PLAB LE 6
66 C*** GEV/C LABORATORY MOMENTUM REGION
67 C*** N    - PROJECTILE HADRON INDEX
68 C*** PLAB - LABORATORY MOMENTUM OF N (GEV/C)
69 C*** ELAB - LABORATORY ENERGY OF N (GEV)
70 C*** CX,CY,CZ - DIRECTION COSINES OF N IN THE LABORATORY SYSTEM
71 C*** ITTA - TARGET NUCLEON INDEX
72 C*** OUTPUT VARIABLES LIST OF PARTICLE CHARACTERISTICS IN /FINLSP/
73 C  IR COUNTS THE NUMBER OF PRODUCED PARTICLES
74 C*** ITR - PARTICLE INDEX, CXR,CYR,CZR - DIRECTION COSINES (LAB. SYST.)
75 C*** ELR,PLR LAB. ENERGY AND LAB. MOMENTUM OF THE SAMPLED PARTICLE
76 C*** RESPECT., UNITS (GEV/C AND GEV)
77 C----------------------------
78       LOWP=0
79       IF (ITPRF(  N  ).LT.0) GO TO 99999
80       WRITE(LUNOUT,99998) N
81 *     STOP    Commented out A. Fasso' 1989
82       IR=0
83       RETURN
84 99998 FORMAT (3(5H ****/),
85      *45H FALSE USE OF THE PARTICLE TYPE INDEX, VALUE ,I4,3(5H ****/))
86 99999 CONTINUE
87       IATMPT=0
88          IF (ABS(PLAB-5.D0).GE.4.99999D0) THEN
89             WRITE(LUNOUT,99996) PLAB
90 *           STOP     Commented out A. Fasso' 1989
91             IR=0
92             RETURN
93 99996       FORMAT (3(5H ****/),64
94      *H PROJECTILE HADRON MOMENTUM OUTSIDE OF THE ALLOWED REGION, PLAB=,
95      *1E15.5/,3(5H ****/))
96       END IF
97       INEWHD = N + 1000 * ITTA
98       IF ( INEWHD .NE. IOLDHD ) THEN
99          CALL CALUMV (N,ITTA)
100       ELSE IF ( (AMPROT-AM(1)) .GT. 1.D-6 ) THEN
101          CALL CALUMV (N,ITTA)
102          INEWHD = - INEWHD
103       END IF
104       IOLDHD = INEWHD
105  1009 CONTINUE
106       IATMPT=0
107       LOWP=LOWP+1
108  1000 CONTINUE
109       IMACH=0
110       REDU=2.D0
111       IW1=0
112       IF (LOWP.GT.20) GO TO 8
113       NNN=N
114 *  +-------------------------------------------------------------------*
115 *  |  The following condition is never verified
116       IF (NNN.NE.N) THEN
117          RUNTES=0.D0
118          EFTES=0.D0
119       END IF
120 *  |
121 *  +-------------------------------------------------------------------*
122       IS=1
123       IR=0
124       IST=1
125       NSTAB=25
126 *  +-------------------------------------------------------------------*
127 *  |  Select the reaction channel Ire: proton target
128       IF ( ITTA .EQ. 1 ) THEN
129          IRE = NURE(N,1)
130 *  |
131 *  +-------------------------------------------------------------------*
132 *  |  neutron target
133       ELSE
134          IRE = NURE(N,2)
135       END IF
136 *  |
137 *  +-------------------------------------------------------------------*
138 *  Elastic scattering index:
139       IELSCT = MIN (N,ITTA) + 1000 * MAX (N,ITTA)
140 C
141 C-----------------------------
142 C*** IE,AMT,ECM,SI DETERMINATION
143 C----------------------------
144       CALL FKSIGI(IRE,PLAB,N,IE,AMT,AMN,ECM,SI,ITTA)
145 *  +------------------------------------------------------------------*
146 *  |  Check if masses have been changed for annihilation treated with
147 *  |  pseudo masses
148       IF ( AMPROT - AM(1) .GT. 1.D-6 ) THEN
149          IANTH = 1
150          SI = 1.D0
151 *  |
152 *  +------------------------------------------------------------------*
153 *  |
154       ELSE
155          IANTH = -1
156       END IF
157 *  |
158 *  +------------------------------------------------------------------*
159       ECMMH=ECM
160 C
161 C-----------------------------
162 C    ENERGY INDEX
163 C  IRE CHARACTERIZES THE REACTION
164 C  IE IS THE ENERGY INDEX
165 C----------------------------
166       IF (SI.LT.1.D-10) GO TO 8
167 * The following condition should be always verified
168       IF (N .LE.NSTAB) GO TO 1
169       RUNTES=RUNTES+1.D0
170       IF (RUNTES.LT.20.D0) WRITE(LUNOUT,602)N
171  602   FORMAT(3H N=,I10,30H THE PROEKTILE IS A RESONANCE  )
172       IF(IBAR(N).EQ.1) N=8
173       IF(IBAR(N).EQ.-1)  N=9
174 *  **** Come here ( 1 continue ) every time we unsuccessfully selected
175 *       for more than 50 times the mass of the produced resonaces ****
176    1  CONTINUE
177       IMACH=IMACH+1
178       IF (IMACH.GT.10) GO TO 8
179       ECM =ECMMH
180       AMN2=AMN**2
181       AMT2=AMT**2
182 *  CMS energy of the projectile
183       ECMN=(ECM**2+AMN2-AMT2)/(2.D0*ECM)
184 *  It should never happen
185       IF(ECMN.LE.AMN) ECMN=AMN
186 *  CMS momentum of the projectile (and of the target)
187       PCMN=SQRT(ECMN**2-AMN2)
188       GAM=(ELAB+AMT)/ECM
189       BGAM=PLAB/ECM
190       IF (IANTH.GE.0) ECM=2.1D0
191 *  From this point starts the random choiche of the reaction channel:
192 *  it was extensively modified by A. Ferrari
193       IST=0
194 * Initial energy index for the reaction IRE (index 0)
195       IIEI=IEII(IRE)
196 * Number of energy intervals for the reaction IRE
197       IDWK=IEII(IRE+1)-IIEI
198 * Initial index for the exit channel weights
199       IIWK=IRII(IRE)
200 * Initial index (for 0) of the exit channels
201       IIKI=IKII(IRE)
202 * Number of exit channels of reaction ire
203       IKE =IKII(IRE+1)-IIKI
204 * *** Shrinkage to the considered energy region for the use of weights
205       HECM=ECM
206 * The following cards assure that Ecm =< Umax + DUmax for this reaction
207 * where:
208 *       Umax  = max cms energy at which data are tabulated
209 *       DUmax = width of the last interval for the tabulated data
210       HUMO=2.D0*UMO(IIEI+IDWK)-UMO(IIEI+IDWK-1)
211       IF (HUMO.LT.ECM) ECM=HUMO
212 * *** Interpolation preparation
213 * Cms energy of the upper limit of the considered energy interval
214       ECMO=UMO(IE)
215 * Cms energy of the lower limit of the considered energy interval
216       ECM1=UMO(IE-1)
217 * Width of the considered interval
218       DECM=ECMO-ECM1
219 * Width from actual value to the lower limit (note that if Ecm > Ecmo
220 * it can be larger than Decm but it is always less than 2xdecm for the
221 * above condition on Humo
222       DEC0=ECM -ECM1
223 * Set to 1 the default total weight
224       WACCUM = 1.D+00
225       WCUMCH (0) = 0.D+00
226       WCUMIE = 0.D+00
227       WCUMI0 = 0.D+00
228       WCSUM0 = 0.D+00
229       CALL GRNDM(RNDM,1)
230       RNDMIK = RNDM (1)
231       IOUT1  = NRK (1,IIKI+1)
232       IOUT2  = NRK (2,IIKI+1)
233       IELCHK = MIN ( IOUT1, IOUT2 ) + 1000 * MAX ( IOUT1, IOUT2 )
234 *  +-------------------------------------------------------------------*
235 *  |  Look for "inverse" reactions for which the elastic channel is not
236 *  |  the first: for these reactions we must exchange the weight of
237 *  |  the elastic channel with the charge exchange one at minimum, to
238 *  |  fulfill the detailed balance theorem
239       IF ( IELCHK .NE. IELSCT ) THEN
240          LSWAP  = .TRUE.
241          IELCHA = IKCHXG (IRE)
242          ICXCHA = 1
243          IKIK (1) = IELCHA
244          IKIK (IELCHA) = 1
245 *  |
246 *  +-------------------------------------------------------------------*
247 *  |  Loop to find the elastic channel
248       ELSE
249          LSWAP  = .FALSE.
250          IELCHA = 1
251          ICXCHA = IKCHXG (IRE)
252       END IF
253 *  |
254 *  +-------------------------------------------------------------------*
255 *  +-------------------------------------------------------------------*
256 *  |  Loop on the exit channels:
257 *  |                              Jkk = index for weights
258 *  |                              Ik  = index for reaction
259 *  |  usually they are the same, but for "inverse" reactions where the
260 *  |  elastic and the charge exchange channels are exchanged, for these
261 *  |  two channels they are crossed
262       DO 2000 JKK = 1, IKE
263 *  |  Ik : index of the exit channel under consideration
264          IK = IKIK (JKK)
265 *  |  Get the index for the weight of ikth channel at ieth energy (upper
266 *  |  limit of the interval in energy)
267          IWK = IIWK + (JKK-1) * IDWK + IE - IIEI
268 *  |  Cumulative Weight of channels 1...ik at energy Ie
269          WIEK   = WK (IWK)
270 *  |  +----------------------------------------------------------------*
271 *  |  |  Check if we are in the first interval: is so all the weights
272 *  |  |  are set to zero for Ie-1, for all channels, so set them to the
273 *  |  |  same value as for Ie to get the proper normalization to 1,
274 *  |  |  then possible thresholds are accounted for after
275          IF ( IE .LE. IIEI+2 ) THEN
276             WIEM1K = WIEK
277 *  |  |
278 *  |  +----------------------------------------------------------------*
279 *  |  |  Cumulative Weight of channels 1...ik at energy Ie-1
280          ELSE
281             WIEM1K = WK (IWK-1)
282          END IF
283 *  |  |
284 *  |  +----------------------------------------------------------------*
285 *  |  Compute the weight at Ie for this channel
286          WIEK   = WIEK   - WCUMIE
287 *  |  Compute the weight at Ie-1 for this channel
288          WIEM1K = WIEM1K - WCUMI0
289 *  |  +----------------------------------------------------------------*
290 *  |  |  This channel is not open at energies Ie and Ie-1
291          IF ( WIEM1K + WIEK .LE. ANGLGB ) THEN
292             WCHANN (JKK) = 0.D+00
293             WCUMCH (JKK) = WCUMCH (JKK-1)
294 *  |  |
295 *  |  +----------------------------------------------------------------*
296 *  |  |
297          ELSE
298 *  |  |  Set Eklim to a negative value to flag for Iefun it is a
299 *  |  |  cms energy and not a lab momentum: this is the cms threshold
300 *  |  |  for the exit channel ik
301             EKLIM = - THRESH (IIKI+IK)
302 *  |  |  Iefun returns the energy index of upper limit of the interval
303 *  |  |  containing the threshold
304             IELIM = IEFUN (EKLIM,IRE)
305             EKLIM = - EKLIM
306             WCHAN0 = WIEM1K + ( WIEK - WIEM1K ) * DEC0 / DECM
307             WCSUM0 = WCSUM0 + WCHAN0
308 *  |  |  +-------------------------------------------------------------*
309 *  |  |  |  Check if we are below threshold
310             IF ( ECM .LE. EKLIM ) THEN
311                WCHANN (JKK) = 0.D+00
312                WCUMCH (JKK) = WCUMCH (JKK-1)
313                WCUMIE = WCUMIE + WIEK
314                WCUMI0 = WCUMI0 + WIEM1K
315                RNDRED = WCHAN0
316                WACCUM = WACCUM - RNDRED
317 *  |  |  |
318 *  |  |  +-------------------------------------------------------------*
319 *  |  |  |  We are above threshold
320             ELSE
321                WCUMIE = WCUMIE + WIEK
322                WCUMI0 = WCUMI0 + WIEM1K
323                ECMD   = MAX ( EKLIM, ECM1 )
324                DEC    = ECM  - ECMD
325                DECC   = ECMO - ECMD
326                WCHANN (JKK) = WIEM1K + ( WIEK - WIEM1K ) * DEC / DECC
327 *  |  |  |  If we are beyond the last tabulated point and the xsec for
328 *  |  |  |  this channel is going down it can happen that Wchann < 0
329 *  |  |  |  set it to 0 and correct according to the usual formalism
330                WCHANN (JKK) = MAX ( WCHANN (JKK), ZERZER )
331                RNDRED = WCHAN0 - WCHANN (JKK)
332                WACCUM = WACCUM - RNDRED
333 *  |  |  |  +----------------------------------------------------------*
334 *  |  |  |  |  Ielflg check: first of all check if this channel
335 *  |  |  |  |  is the elastic one and if so if reduction must be applied
336 *  |  |  |  |  to
337                IF ( IK .EQ. IELCHA .AND. IELFLG .NE. 0 ) THEN
338 *  |  |  |  |  +-------------------------------------------------------*
339 *  |  |  |  |  |  Elastic scattering reduced due to collision in the
340 *  |  |  |  |  |  nucleus and not with a free proton/neutron
341                   IF ( IELFLG .LT. 0 ) THEN
342                      IF ( IBAR (N) .NE. 0 ) THEN
343                         REDUC = PLAB / PPAMXB
344                         REDUC = PAUMXB * REDUC / ( 1.D+00 + REDUC**2 )
345                      ELSE
346                         REDUC = PLAB / PPAMXM
347                         REDUC = PAUMXM * REDUC / ( 1.D+00 + REDUC**2 )
348                      END IF
349                      RNDRED = WCHANN (JKK) * ( 1.D+00 - REDUC )
350                      WCHANN (JKK) = WCHANN (JKK) - RNDRED
351 *  |  |  |  |  |
352 *  |  |  |  |  +-------------------------------------------------------*
353 *  |  |  |  |  |  Elastic scattering forbidden
354                   ELSE
355                      RNDRED = WCHANN (JKK)
356                      WCHANN (JKK) = 0.D+00
357                   END IF
358 *  |  |  |  |  |
359 *  |  |  |  |  +-------------------------------------------------------*
360                   WACCUM = WACCUM - RNDRED
361 *  |  |  |  |
362 *  |  |  |  +----------------------------------------------------------*
363 *  |  |  |  |  Icxflg check: first of all check if this channel
364 *  |  |  |  |  is the charge exchange one and if so if reduction must
365 *  |  |  |  |  be applied to
366                ELSE IF ( IK .EQ. ICXCHA .AND. ICXFLG .NE. 0 )THEN
367 *  |  |  |  |  +-------------------------------------------------------*
368 *  |  |  |  |  |  Charge exchange reduced due to collision in the
369 *  |  |  |  |  |  nucleus and not with a free proton/neutron
370                   IF ( ICXFLG .LT. 0 ) THEN
371                      IF ( IBAR (N) .NE. 0 ) THEN
372                         REDUC = PLAB / PPAMXB
373                         REDUC = PAUMXB * REDUC / ( 1.D+00 + REDUC**2 )
374                      ELSE
375                         REDUC = PLAB / PPAMXM
376                         REDUC = PAUMXM * REDUC / ( 1.D+00 + REDUC**2 )
377                      END IF
378                      RNDRED = WCHANN (JKK) * ( 1.D+00 - REDUC )
379                      WCHANN (JKK) = WCHANN (JKK) - RNDRED
380 *  |  |  |  |  |
381 *  |  |  |  |  +-------------------------------------------------------*
382 *  |  |  |  |  | Charge exchange forbidden
383                   ELSE
384                      RNDRED = WCHANN (JKK)
385                      WCHANN (JKK) = 0.D+00
386                   END IF
387 *  |  |  |  |  |
388 *  |  |  |  |  +-------------------------------------------------------*
389                   WACCUM = WACCUM - RNDRED
390                END IF
391 *  |  |  |  |
392 *  |  |  |  +----------------------------------------------------------*
393                WCUMCH (JKK) = WCUMCH (JKK-1) + WCHANN (JKK)
394             END IF
395 *  |  |  |
396 *  |  |  +-------------------------------------------------------------*
397             RNDMCH = RNDMIK * WACCUM
398 *  |  |  +-------------------------------------------------------------*
399 *  |  |  |  Check if a possible decrease/increase of this channel
400 *  |  |  |  opened one of the already examinated channels
401             IF ( RNDMCH .LT. WCUMCH (JKK-1) ) THEN
402 *  |  |  |  +----------------------------------------------------------*
403 *  |  |  |  |  Loop on the previous channels
404                DO 1900 JPK = 1, JKK - 1
405                   IF ( RNDMCH .LT. WCUMCH (JPK) ) THEN
406                      IK = IKIK (JPK)
407                      GO TO 2100
408                   END IF
409  1900          CONTINUE
410 *  |  |  |  |
411 *  |  |  |  +----------------------------------------------------------*
412 *  |  |  |
413 *  |  |  +-------------------------------------------------------------*
414 *  |  |  |  Check if this one is the right channel
415             ELSE IF ( RNDMCH .LT. WCUMCH (JKK) ) THEN
416                GO TO 2100
417             END IF
418 *  |  |  |
419 *  |  |  +-------------------------------------------------------------*
420          END IF
421 *  |  |
422 *  |  +----------------------------------------------------------------*
423  2000 CONTINUE
424 *  |
425 *  +-------------------------------------------------------------------*
426       IK = IELCHA
427       WRITE (LUNERR,*)' **** Hadriv: elastic channel selected when',
428      &                ' prohibited ! ****',N,ELAB,PLAB
429 *  Finally we selected channel ik
430  2100 CONTINUE
431 * *** Ik is the reaction channel ***
432 *  +-------------------------------------------------------------------*
433 *  |  Set to the default values the array Ikik
434       IF ( LSWAP ) THEN
435          IKIK (1) = 1
436          IKIK (IELCHA) = IELCHA
437       END IF
438 *  |
439 *  +-------------------------------------------------------------------*
440       INRK=IKII(IRE)+IK
441 *  First resonance to be created
442       IT1=NRK(1,INRK)
443 *  Second resonance to be created
444       IT2=NRK(2,INRK)
445       ECM=HECM
446       I1001 =0
447 *  +-------------------------------------------------------------------*
448 *  |  Rejection loop for the choiche of the resonance masses
449  1001 CONTINUE
450          IF (I1001.GT.50) GO TO 1
451 *  |  Selection of the resonance mass according to its width
452          AM1=AMGA(IT1)
453 *  |  Selection of the resonance mass according to its width
454          AM2=AMGA(IT2)
455          AMS=AM1+AM2
456          I1001=I1001+1
457       IF ( IT2*AMS .GT. IT2*ECM ) GO TO 1001
458 *  |--<--<--<--<--<--<  Loop back if m1+m2 > Ecm
459 *  +-------------------------------------------------------------------*
460       IT11=IT1
461       IT22=IT2
462       IF (IANTH.GE.0) ECM=ELAB+AMT+0.00000001D0
463       AM11=AM1
464       AM22=AM2
465 *  +-------------------------------------------------------------------*
466 *  |  Direct (single) resonances
467       IF (IT2.LE.0) THEN
468 *  | Random choice of decay channels of the direct resonance  it1
469          KZ1=K1(IT1)
470          IST=IST+1
471          IECO=0
472 *  |   Here was the mistake in the pseudo-masses treatment!!!!!
473 *        ECO=ECM
474          ECO=ECMMH
475          GAM=(ELAB+AMT)/ECO
476          BGAM=PLAB/ECO
477          CXS(1)=CX
478          CYS(1)=CY
479          CZS(1)=CZ
480          GO TO 310
481       END IF
482 *  |
483 *  +-------------------------------------------------------------------*
484 *  +-------------------------------------------------------------------*
485 *  |
486       CALL GRNDM(RNDM,1)
487       IF ( RNDM(1) .LT. 0.5D+00 ) THEN
488          IT1=IT22
489          IT2=IT11
490          AM1=AM22
491          AM2=AM11
492       END IF
493 *  |
494 *  +-------------------------------------------------------------------*
495 C
496 C-----------------------------
497 C   THE FIRST PARTICLE IS DEFINED TO BE THE FORWARD GOING ONE AT SMALL T
498       IBN=IBAR(N)
499       IB1=IBAR(IT1)
500       IT11=IT1
501       IT22=IT2
502       AM11=AM1
503       AM22=AM2
504 *  +-------------------------------------------------------------------*
505 *  |
506       IF(IB1.NE.IBN) THEN
507          IT1=IT22
508          IT2=IT11
509          AM1=AM22
510          AM2=AM11
511       END IF
512 *  |
513 *  +-------------------------------------------------------------------*
514 C-----------------------------
515 C***IT1,IT2 ARE THE CREATED PARTICLES
516 C***MOMENTA AND DIRECTION COSINA IN THE CM - SYSTEM
517 C------------------------
518       CALL TWOPAR(ECM1,ECM2,PCM1,PCM2,COD1,COD2,COF1,COF2,SIF1,SIF2,
519      *IT1,IT2,ECM,ECMN,PCMN,N,AM1,AM2)
520       IST=IST+1
521       ITS(IST)=IT1
522       AMM(IST)=AM1
523 C
524 C-----------------------------
525 C***TRANSFORMATION INTO LAB SYSTEM AND ROTATION
526 C----------------------------
527       CALL TRAFO(GAM,BGAM,CX,CY,CZ,COD1,COF1,SIF1,PCM1,ECM1,PLS(IST),
528      *CXS(IST),CYS(IST),CZS(IST),ELS(IST))
529       IST=IST+1
530       ITS(IST)=IT2
531       AMM(IST)=AM2
532       CALL TRAFO(GAM,BGAM,CX,CY,CZ,COD2,COF2,SIF2,PCM2,ECM2,PLS(IST),CXS
533      *(IST),CYS(IST),CZS(IST),ELS(IST))
534   200 CONTINUE
535 C
536 C-----------------------------
537 C***TEST   STABLE OR UNSTABLE
538 C----------------------------
539       IF(ITS(IST).GT.NSTAB) GO TO 300
540       IR=IR+1
541 C
542 C-----------------------------
543 C***IR IS THE NUMBER OF THE FINAL STABLE PARTICLE
544 C----------------------------
545       IF (REDU.LT.0.D0) GO TO 1009
546       ITR(IR)=ITS(IST)
547       PLR(IR)=PLS(IST)
548       CXR(IR)=CXS(IST)
549       CYR(IR)=CYS(IST)
550       CZR(IR)=CZS(IST)
551       ELR(IR)=ELS(IST)
552       IST=IST-1
553       IF(IST.GE.1) GO TO 200
554          GO TO 500
555   300 CONTINUE
556 C
557 C  RANDOM CHOICE OF DECAY CHANNELS
558 C----------------------------
559 C
560       IT=ITS(IST)
561       ECO=AMM(IST)
562       GAM=ELS(IST)/ECO
563       BGAM=PLS(IST)/ECO
564       IECO=0
565       KZ1=K1(IT)
566   310 CONTINUE
567       IECO=IECO+1
568       CALL GRNDM(RNDM,1)
569       VV=RNDM(1)
570       IIK=KZ1-1
571   301 CONTINUE
572          IIK=IIK+1
573       IF (VV.GT.WT(IIK)) GO TO 301
574 C
575 C  IIK IS THE DECAY CHANNEL
576 C----------------------------
577       IT1=NZK(IIK,1)
578       I310=0
579  1310 CONTINUE
580          I310=I310+1
581          AM1=AMGA(IT1)
582          IT2=NZK(IIK,2)
583          AM2=AMGA(IT2)
584          IF (IT2-1.LT.0) GO TO 110
585          IT3=NZK(IIK,3)
586          AM3=AMGA(IT3)
587          AMS=AM1+AM2+AM3
588 C
589 C  IF  IIK-KIN.LIM.GT.ACTUAL TOTAL CM-ENERGY, DO AGAIN RANDOM IIK-CHOICE
590 C----------------------------
591          IF (IECO.GT.10) THEN
592             IATMPT=IATMPT+1
593 * Note: we can go to 8 also for too many iterations
594             IF (IATMPT.GT.3) GO TO 8
595             GO TO 1000
596          END IF
597          IF (I310.GT.50) GO TO 310
598       IF (AMS.GT.ECO) GO TO 1310
599 C
600 C  FOR THE DECAY CHANNEL
601 C  IT1,IT2, IT3 ARE THE PRODUCED PARTICLES FROM  IT
602 C----------------------------
603       IF (REDU.LT.0.D0) GO TO 1009
604       ITWTHC=0
605       REDU=2.D0
606       IF(IT3.EQ.0) GO TO 400
607  4001 CONTINUE
608       ITWTH=1
609       CALL THREPD(ECO,ECM1,ECM2,ECM3,PCM1,PCM2,PCM3,COD1,COF1,SIF1,
610      *COD2,COF2,SIF2,COD3,COF3,SIF3,AM1,AM2,AM3)
611       GO TO 411
612   400 CONTINUE
613       CALL TWOPAD(ECO,ECM1,ECM2,PCM1,PCM2,COD1,COF1,SIF1,COD2,COF2,SIF2,
614      *AM1,AM2)
615       ITWTH=-1
616       IT3=0
617   411 CONTINUE
618       ITWTHC=ITWTHC+1
619       IF (REDU.GT.0.D0) GO TO 110
620       REDU=2.D0
621       IF (ITWTHC.GT.100) GO TO 1009
622       IF (ITWTH) 400,400,4001
623   110 CONTINUE
624       ITS(IST)=IT1
625       IF (IT2.LE.0) GO TO 305
626       ITS(IST+1)=IT2
627       ITS(IST+2)=IT3
628       RX=CXS(IST)
629       RY=CYS(IST)
630       RZ=CZS(IST)
631       AMM(IST)=AM1
632       CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD1,COF1,SIF1,PCM1,ECM1,
633      *PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
634       IST=IST+1
635       AMM(IST)=AM2
636       CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD2,COF2,SIF2,PCM2,ECM2,
637      *PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
638       IF (IT3.LE.0) GO TO 305
639       IST=IST+1
640       AMM(IST)=AM3
641       CALL TRAFO(GAM,BGAM,RX,RY,RZ,COD3,COF3,SIF3,PCM3,ECM3,
642      *PLS(IST),CXS(IST),CYS(IST),CZS(IST),ELS(IST))
643   305 CONTINUE
644       GO TO 200
645   500 CONTINUE
646   631 CONTINUE
647       RETURN
648     8 CONTINUE
649 C
650 C----------------------------
651 C
652 C   ZERO CROSS SECTION CASE
653 C
654 C----------------------------
655 C
656          IR=1
657          ITR(1)=N
658          CXR(1)=CX
659          CYR(1)=CY
660          CZR(1)=CZ
661          ELR(1)=ELAB
662          PLR(1)=PLAB
663       RETURN
664       END