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