GEANT321/peanut/peanut.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:22:02  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.46  by  S.Giani
  11 *-- Author :
  12 *$ CREATE PEANUT.FOR
  13 *COPY PEANUT
  14 *
  15 *=== peanut ===========================================================*
  16 *
  17       SUBROUTINE PEANUT ( KPROJ, EKE, PPROJ, TXX, TYY, TZZ, WEE )
  18
  19 #include "geant321/dblprc.inc"
  20 #include "geant321/dimpar.inc"
  21 #include "geant321/iounit.inc"
  22 *
  23 *----------------------------------------------------------------------*
  24 *----------------------------------------------------------------------*
  25 *
  26 *
  27 #include "geant321/balanc.inc"
  28 #include "geant321/eva0.inc"
  29 #include "geant321/fheavy.inc"
  30 #include "geant321/finuc.inc"
  31 #include "geant321/higfis.inc"
  32 #include "geant321/nucdat.inc"
  33 #include "geant321/nucgeo.inc"
  34 #include "geant321/parevt.inc"
  35 #include "geant321/parnuc.inc"
  36 #include "geant321/part.inc"
  37 #include "geant321/resnuc.inc"
  38 *
  39 *
  40       COMMON / FKCOSP / C1ST (3), C2ND (3), LEMISS
  41       LOGICAL LEMISS
  42       COMMON / FKCMCY / NPCYCL (20,2), IEVT, LOUT
  43 *
  44       COMMON / FKPLOC / IABCOU
  45       LOGICAL LBCHCK, LBIMPC, LTRPPD, LASSOR, LEXIT, LPRCYC, LEXPLC,
  46      &        LNWINT
  47       DIMENSION IPTYPE (39)
  48       REAL RNDM(2)
  49       SAVE IPTYPE
  50       DATA IPTYPE /  1,  2,  0,  0,  0,  0,  0,  1,  2,  0,  0,  0,
  51      &               3,  3,  4,  4,  5,  5,  0,  5,  5,  5,  3,  4,
  52      &               4,  0,  0,  0,  0,  0,  6,  6,  6,  7,  8,  7,
  53      &               8,  9, 10 /
  54 *
  55 *
  56       IEVPRE = IEVPRE + 1
  57       NUSCIN = 0
  58       IABCOU = 0
  59       IF ( EKE .GT. 2.D+00 * GAMMIN ) THEN
  60          EOTEST = ETTOT
  61          PTORI  = PPROJ
  62          PXORI  = PTORI * TXX
  63          PYORI  = PTORI * TYY
  64          PZORI  = PTORI * TZZ
  65       ELSE
  66          EOTEST = ETTOT - EKE
  67          ETTOT  = EOTEST
  68          PTORI  = 0.D+00
  69          PXTTOT = 0.D+00
  70          PYTTOT = 0.D+00
  71          PZTTOT = 0.D+00
  72          PTTOT  = 0.D+00
  73          PXORI  = 0.D+00
  74          PYORI  = 0.D+00
  75          PZORI  = 0.D+00
  76       END IF
  77       ETEPS  = 1.D-10 * ETTOT
  78       ICHTOT = ICHTAR + ICH  (KPROJ)
  79       IBTOT  = IBTAR  + IBAR (KPROJ)
  80       IBNUCL = IBTAR
  81       IBORI  = IBAR   (KPROJ)
  82       IPTORI = IPTYPE (KPTOIP(KPROJ))
  83       KPORI  = KPROJ
  84       EKORI  = EKE
  85       ZZTAR  = ICHTAR
  86       BBTAR  = IBTAR
  87       IF ( ICH (KPROJ) .NE. 0 .AND. EKE .GT. 2.D+00 * GAMMIN ) THEN
  88          FLKCOU = DOST ( 1, ZZTAR )
  89          CCOUL  = DOST ( 3, ZZTAR )
  90          CLMBBR = ICH (KPROJ) * COULBH * ZZTAR / RMASS (IBTAR)
  91          IF ( CLMBBR .GT. 0.9D+00 * EKE ) THEN
  92             TMPEKE = 0.98D+00 * EKE
  93             CLMHLP = MIN ( CLMBBR * FLKCOU, TMPEKE )
  94             CLMBBR = MIN  ( CLMBBR, EKE )
  95             WEIGH1 = 10.0D+00 * ( CLMBBR / EKE - 0.9D+00 )
  96             CLMBBR = WEIGH1 * CLMHLP + ( 1.D+00 - WEIGH1 ) * CLMBBR
  97          END IF
  98          BFCLMB = SQRT ( 1.D+00 - CLMBBR / EKE )
  99          RDCLMB = ICH (KPROJ) * COULPR * ZZTAR / CLMBBR
 100       ELSE
 101          CLMBBR = 0.D+00
 102          BFCLMB = 1.D+00
 103          RDCLMB = AINFNT
 104       END IF
 105       IBRES = IBTOT
 106       ICRES = ICHTOT
 107       BBRES = IBRES
 108       ZZRES = ICRES
 109       AMMRES = BBRES * AMUAMU + 0.001D+00 * FKENER ( BBRES, ZZRES )
 110       AMNRES = AMMRES - ZZRES * AMELEC + ELBNDE ( ICRES )
 111       NPROT = 0
 112       NNEUT = 0
 113       NHOLE = 0
 114       AVEBIN = ( ( BBTAR - ZZTAR ) * AMNUCL (2) + ZZTAR * AMNUCL (1)
 115      &       - AMNTAR ) / BBTAR
 116       AMMHLP = ( BBTAR - 1.D+00 ) * AMUAMU + 1.D-03
 117      &       * FKENER ( BBTAR - ONEONE, ZZTAR - ONEONE )
 118       AMNHLP = AMMHLP - ( ZZTAR - 1.D+00 ) * AMELEC + ELBNDE (ICHTAR-1)
 119       BNENRG (1) = AMNHLP + AMNUCL (1) - AMNTAR
 120       AMMHLP = ( BBTAR - 1.D+00 ) * AMUAMU + 1.D-03
 121      &       * FKENER ( BBTAR - ONEONE, ZZTAR )
 122       AMNHLP = AMMHLP - ZZTAR * AMELEC + ELBNDE (ICHTAR)
 123       BNENRG (2) = AMNHLP + AMNUCL (2) - AMNTAR
 124
 125       BNENRG (3) = 0.5D+00 * ( BNENRG (1) + BNENRG (2) )
 126       RHORED = 1.D+00
 127       NPNUC  = 0
 128       NNUCTS = 0
 129       NHLEXP = 0
 130       JNUCTS = 0
 131       ACOLL  = ANOW
 132       ZCOLL  = ZNOW
 133       IF ( .NOT. LPREEX ) THEN
 134          IF ( KPROJ .EQ. 1 .OR. KPROJ .EQ. 8 ) THEN
 135             LEXPLC = EKE .GT. 0.250D+00
 136          ELSE IF ( KPROJ .EQ. 14 .AND. PTTOT .LE. 0.D+00 ) THEN
 137             LEXPLC = .TRUE.
 138          ELSE
 139             STOP 'LEXPLC'
 140          END IF
 141       ELSE
 142          LEXPLC = .TRUE.
 143       END IF
 144       IF ( LEXPLC .AND. EKE .GT. IBAR(KPROJ) * EKEEXP ) THEN
 145          ICYCL  = 0
 146          IREINT = 0
 147          LPRCYC = .TRUE.
 148          LBCHCK = .FALSE.
 149          LBIMPC = .TRUE.
 150          LELSTC = .FALSE.
 151          LABRST = .FALSE.
 152          LABSRP = .FALSE.
 153          LINELS = .FALSE.
 154          LCHEXC = .FALSE.
 155          RHOEXP = 0.D+00
 156          EKFEXP = 0.D+00
 157          EKFREI = 0.D+00
 158          RHOREI = 0.D+00
 159          KPRIN  = KPROJ
 160          KRFLIN = 0
 161          ERECLD = 0.D+00
 162          BNPREV = 0.D+00
 163          EKECON = EKE
 164          PNUCCO = PPROJ
 165          CALL BIMSEL ( KPROJ, TXX, TYY, TZZ, LBCHCK )
 166          WLLPRO = WLLRED
 167          BNPROJ = WLLRED * BNDNUC
 168          RHOMEM = 0.5D+00 * ( RHOIMP + RHOIMT )
 169          EKFMEM = 0.5D+00 * ( EKFIMP + EKFPRO )
 170          IAAFT  = IBRES - IBAR (KPROJ)
 171          IZAFT  = ICRES - ICH  (KPROJ)
 172          CALL EEXLVL ( IAAFT, IZAFT, EEXDEL, EEXMIN, EEXDUM )
 173          EEXANY = EEXDEL
 174          IF ( ALPHAL .LE. 0.D+00 ) THEN
 175             DEFPRO = 0.D+00
 176             DEFNEU = 0.D+00
 177          ELSE
 178             DEFPRO = 0.D+00
 179             DEFNEU = 0.D+00
 180          END IF
 181          DEFNUC (1) = DEFPRO
 182          DEFNUC (2) = DEFNEU
 183   100    CONTINUE
 184          IF ( LELSTC ) THEN
 185             CALL NUCNUC ( IKPMX , KRFLIN, WEE   , ERECMN, LBIMPC,
 186      &                    LBCHCK, ICYCL , NHOLE , NPROT , NNEUT ,
 187      &                    LEXIT , LNWINT )
 188          ELSE IF ( KPRIN .EQ. 14 ) THEN
 189             CALL PIOABS ( IKPMX , KRFLIN, WEE   , ERECMN, LBIMPC,
 190      &                    LBCHCK, ICYCL , NHOLE , NPROT , NNEUT ,
 191      &                    LEXIT , LNWINT )
 192          ELSE
 193             STOP 'Int_kind'
 194          END IF
 195          IF ( LNWINT ) GO TO 100
 196          BBRES = IBRES
 197          ZZRES = ICRES
 198          IF ( .NOT. LEXIT ) THEN
 199             LPRCYC = .FALSE.
 200          ELSE
 201             BNPREV = BNPREV + BNDUSE
 202          END IF
 203   200    CONTINUE
 204          LELSTC = .FALSE.
 205          LABRST = .FALSE.
 206          LABSRP = .FALSE.
 207          LINELS = .FALSE.
 208          LCHEXC = .FALSE.
 209          GAMMAX = 0.D+00
 210          IREFMN = 10000
 211          IKPMX  = 0
 212          IBCHCK = 0
 213          ICCHCK = 0
 214          IBNUCL = 0
 215          ICNUCL = 0
 216          ILIVE  = 0
 217          LTRPPD = .FALSE.
 218          DO 300 KP = 1, NPNUC
 219             IF ( KPNUCL (KP) .LE. 0 ) GO TO 300
 220             ILIVE = ILIVE + 1
 221             KPNUC = KPNUCL (KP)
 222             IPTNUC = IPTYPE (KPTOIP(KPNUC))
 223             IF ( IPTNUC .EQ. 1 ) THEN
 224                BNDNU0 = BNENRG (1+KPNUC/8)
 225                WLLRE0 = POTBAR
 226             ELSE
 227                IF ( IBAR (KPNUC) .NE. 0 ) THEN
 228                   WLLRE0 = POTBAR
 229                   BNDNU0 = BNENRG (3)
 230                ELSE IF ( KPNUC .LE. 11 ) THEN
 231                   WLLRE0 = 0.D+00
 232                   BNDNU0 = 0.D+00
 233                ELSE
 234                   WLLRE0 = POTMES
 235                   BNDNU0 = BNENRG (3)
 236                END IF
 237             END IF
 238             IF ( EKFNUC (KP) .GT. -100.D+00 ) THEN
 239                GAMMA = ( ENNUC (KP) - WLLRE0 * ( EKFNUC (KP) + BNDNU0 )
 240      &               ) / AM (KPNUC)
 241             ELSE
 242                IF ( AM (KPNUC) .LE. ANGLGB ) THEN
 243                   GAMMA = AINFNT
 244                ELSE
 245                   GAMMA = ENNUC (KP) / AM (KPNUC)
 246                END IF
 247             END IF
 248             IF ( IBAR (KPNUC) .GT. 0 ) THEN
 249                IBNUCL = IBNUCL + IBAR (KPNUC)
 250                ICNUCL = ICNUCL + ICH  (KPNUC)
 251             END IF
 252             IBCHCK = IBCHCK + IBAR (KPNUC)
 253             ICCHCK = ICCHCK + ICH  (KPNUC)
 254             IF ( KRFNUC (KP) .LT. IREFMN ) THEN
 255                IREFMN = KRFNUC (KP)
 256                GAMMAX = GAMMA
 257                IKPMX  = KP
 258                WLLRED = WLLRE0
 259                BNDNUC = BNDNU0
 260             ELSE IF ( KRFNUC (KP) .EQ. IREFMN ) THEN
 261                IF ( GAMMA .GT. GAMMAX ) THEN
 262                   GAMMAX = GAMMA
 263                   IKPMX  = KP
 264                   WLLRED = WLLRE0
 265                   BNDNUC = BNDNU0
 266                END IF
 267             END IF
 268   300    CONTINUE
 269          IBNUCL = IBRES - IBNUCL - NPROT - NNEUT
 270          ICNUCL = ICRES - ICNUCL - NPROT
 271          ACOLL  = IBNUCL
 272          ZCOLL  = ICNUCL
 273          RHORED = ACOLL / BBTAR
 274          IF ( IKPMX .LE. 0 ) THEN
 275             IBCKC = IBTOT  - IBINTR - IBNUCR
 276             ICCKC = ICHTOT - ICINTR - ICNUCR
 277             IF ( IBCKC .NE. IBRES .OR. ICCKC .NE. ICRES ) THEN
 278                ICRES = ICCKC
 279                IBRES = IBCKC
 280             END IF
 281             NEXPEM = NP-NP0
 282             DO 450 IJJ = 1,IGREYN
 283                NPCYCL (IJJ,1) = 0
 284   450       CONTINUE
 285             DO 460 IJJ = 1,IGREYP
 286                NPCYCL (IJJ,2) = 0
 287   460       CONTINUE
 288             BBRES = IBRES
 289             ZZRES = ICRES
 290             ANOW  = BBRES
 291             ZNOW  = ZZRES
 292             AMMRES = BBRES * AMUAMU + 0.001D+00 * FKENER ( BBRES, ZZRES)
 293             AMNRES = AMMRES - ZZRES * AMELEC + ELBNDE ( ICRES )
 294             PXRES = PXTTOT - PXNUCR - PXINTR
 295             PYRES = PYTTOT - PYNUCR - PYINTR
 296             PZRES = PZTTOT - PZNUCR - PZINTR
 297             PTRES2= PXRES**2 + PYRES**2 + PZRES**2
 298             PTRES = SQRT ( PTRES2 )
 299             ERES  = ETTOT - EINTR - ENUCR
 300             UMO2  = ( ERES - PTRES ) * ( ERES + PTRES )
 301             IF ( UMO2 .LT. ONEMNS*AMNRES**2 ) THEN
 302                UMO = SQRT (UMO2)
 303                WRITE ( LUNOUT,* )' 2:UMO,AMNRES',UMO,AMNRES
 304                GO TO 530
 305             ELSE IF ( UMO2 - AMNRES*AMNRES .LT. AMNRES*TVEPSI ) THEN
 306                UMO = SQRT (UMO2)
 307                GO TO 530
 308             END IF
 309             IF ( ICYCL .NE. NHOLE - IABCOU ) THEN
 310                WRITE (LUNOUT,*)' *** KPORI, ICYCL, NHOLE, IABCOU',
 311      &                               KPORI, ICYCL, NHOLE, IABCOU
 312                ICYCL = NHOLE - IABCOU
 313             END IF
 314             NPTOT  = NPROT + NNEUT
 315             NHLEXP = NHOLE
 316             IF ( .NOT. LPRCYC .OR. NPTOT .LE. 0 .OR. NNUCTS .GT. 0 )
 317      &         THEN
 318                LEMISS = .FALSE.
 319                NHOLE  = NHOLE + 1
 320                IF ( EKFREI .GT. ANGLGB ) THEN
 321                   RHOIMP = ( RHOEXP + RHOREI ) / NHOLE
 322                   EKFIMP = ( EKFEXP + EKFREI ) / NHOLE
 323                ELSE
 324                   RHOIMP = ( RHOEXP + RHOAVE ) / NHOLE
 325                   EKFIMP = ( EKFEXP + EKFAVE ) / NHOLE
 326                END IF
 327                ANPROT = NPROT
 328                ANNEUT = NNEUT
 329                ACOLL  = ACOLL - 1.D+00
 330                IF ( NPTOT .GT. 0 ) THEN
 331                   PNPROT = ( ZNOW - ANPROT ) * ( 3.D+00 * ANNEUT
 332      &                   + ANPROT ) / ( ANPROT * ( ZNOW - ANPROT
 333      &                   + 3.D+00 * ( ANOW - ANNEUT - ZNOW ) ) + ANNEUT
 334      &                   * ( 3.D+00 * ( ZNOW
 335      &                   - ANPROT ) + ANOW - ANNEUT - ZNOW ) )
 336                ELSE
 337                   PNPROT = ZNOW / ANOW
 338                END IF
 339                CALL GRNDM(RNDM,1)
 340                IF ( RNDM (1) .LT. PNPROT ) THEN
 341                   NPROT = NPROT + 1
 342                   ZCOLL = ZCOLL - 1.D+00
 343                ELSE
 344                   NNEUT = NNEUT + 1
 345                END IF
 346             ELSE
 347                RHOIMP = RHOEXP / NHLEXP
 348                EKFIMP = EKFEXP / NHLEXP
 349                IF ( ICYCL .EQ. 1 ) THEN
 350                   IF ( NP .LE. NP0 ) THEN
 351                      LEMISS = .FALSE.
 352                      ICYCL  = ICYCL - 1
 353                   ELSE
 354                      LEMISS = .TRUE.
 355                   END IF
 356                ELSE
 357                   LEMISS = .FALSE.
 358                   ICYCL  = ICYCL - 1
 359                END IF
 360             END IF
 361             GO TO 500
 362          END IF
 363          KPNUC  = KPNUCL (IKPMX)
 364          IPTNUC = IPTYPE (KPTOIP(KPNUC))
 365          ERECMN = MAX ( ERECLD, ERECMN )
 366          ERECLD = ERECMN
 367          ERECMN = ERECMN / ( ICYCL + IGREYP + IGREYN )
 368          IAAFT  = IBRES - IBAR (KPNUC)
 369          IZAFT  = ICRES - ICH  (KPNUC)
 370          CALL EEXLVL ( IAAFT, IZAFT, EEXDEL, EEXMIN, EEXDUM )
 371          IF ( NP .EQ. NP0 .AND. KPNUC .EQ. KPROJ ) THEN
 372             EEXANY = EEXDEL
 373          ELSE
 374             EEXANY = 0.D+00
 375          END IF
 376          AAFT = BBRES - IBAR (KPNUC)
 377          ZAFT = ZZRES - ICH  (KPNUC)
 378          AMMAFT = AAFT * AMUAMU + 0.001D+00 * FKENER ( AAFT, ZAFT )
 379          AMNAFT = AMMAFT - ZAFT * AMELEC + ELBNDE ( NINT (ZAFT) )
 380          IF ( WLLRED .GT. ANGLGB ) THEN
 381             IF ( EKFNUC (IKPMX) .GT. -100.D+00 ) THEN
 382                BNDGEN = IBAR (KPNUC) * AM (KPNUC) + AMNAFT - AMNRES
 383                IF ( NP .EQ. NP0 .AND. IPTNUC .EQ. IPTORI ) THEN
 384                   BNDUSE = BNPROJ + AMNAFT - AMNTAR + AM (KPNUC)
 385      &                   - AM (KPROJ)
 386                ELSE IF ( NP .EQ. NP0 ) THEN
 387                   IF ( IPTNUC .EQ. 1 ) THEN
 388                      BNDUSE = AMNAFT - AMNTAR + AM (KPNUC)
 389                      BNDUSE = MAX ( BNDUSE, ZERZER )
 390                   ELSE
 391                      BNDUSE = WLLRED * BNDNUC
 392                   END IF
 393                ELSE IF ( NUSCIN .EQ. 1 ) THEN
 394                   AMEMIT = 0.D+00
 395                   DO 430 KP = NP0+1, NP
 396                      IPTPAR = IPTYPE (KPTOIP(KPART(KP)))
 397                      IF ( IPTPAR .EQ. 1 .OR. IPTPAR .EQ. IPTORI )
 398      &               AMEMIT = AMEMIT + AM (KPTOIP(KPART(KP)))
 399   430             CONTINUE
 400                   IF ( IPTNUC .EQ. IPTORI ) THEN
 401                      BNTRUE = AMNAFT + AMEMIT + AM (KPNUC) - AMNTAR
 402      &                      - AM (KPROJ)
 403                      BNDUSE = BNPROJ + BNTRUE - BNPREV
 404                      BNDUSE = MAX ( BNDUSE, ZERZER )
 405                   ELSE IF ( IPTNUC .EQ. 1 ) THEN
 406                      BNDUSE = AMNAFT + AMEMIT + AM (KPNUC) - AMNTAR
 407      &                      - BNPREV
 408                      BNDUSE = MAX ( BNDUSE, ZERZER )
 409                   ELSE
 410                      BNDUSE = WLLRED * BNDNUC
 411                   END IF
 412                ELSE
 413                   BNDUSE = WLLRED * MAX ( BNDGEN, ZERZER )
 414                END IF
 415                EKFPRE = EKFNUC (IKPMX)
 416                VWELL0 = WLLRED * EKFNUC (IKPMX) + BNDUSE + ERECMN
 417                ENNUC  (IKPMX) = ENNUC (IKPMX) - VWELL0
 418                EKFNUC (IKPMX) = WLLRED * BNDNUC - BNDUSE - ERECMN
 419             ELSE
 420                BNDGEN = IBAR (KPNUC) * AM (KPNUC) + AMNAFT - AMNRES
 421                IF ( NP .EQ. NP0 .AND. KPNUC .EQ. KPORI ) THEN
 422                   BNDUSE = BNPROJ
 423                ELSE
 424                   BNDUSE = WLLRED * MAX ( BNDGEN, ZERZER )
 425                END IF
 426                EKFPRE = EKFAVE
 427                RHNUCL (IKPMX) = RHOAVE
 428                VWELL0 = BNDUSE - WLLRED * BNDNUC + ERECMN
 429                ENNUC  (IKPMX) = ENNUC (IKPMX) - VWELL0
 430                EKFNUC (IKPMX) = - VWELL0
 431             END IF
 432          ELSE
 433             VWELL0 = 0.D+00
 434             EKFNUC (IKPMX) = 0.D+00
 435             EKFPRE = 0.D+00
 436             BNDUSE = 0.D+00
 437          END IF
 438          EKNNUC = ENNUC (IKPMX) - AM (KPNUC)
 439          EKECON = EKNNUC - EKFNUC (IKPMX)
 440          IF ( ICH (KPNUC) .GT. 0 ) THEN
 441             FLKCOU = DOST ( 1, ZAFT )
 442             ETHRES = FLKCOU * ICH (KPNUC) * COULBH * ZAFT
 443      &             / RMASS ( NINT (AAFT) )
 444             IF ( EKNNUC .GT. ETHRES ) THEN
 445                LASSOR = .FALSE.
 446                FREJE  = 1.D+00 - ( ETHRES / EKNNUC )**3
 447                CALL GRNDM(RNDM,1)
 448                IF ( RNDM (1) .GE. FREJE ) LTRPPD = .TRUE.
 449             ELSE
 450                LASSOR = .TRUE.
 451             END IF
 452          ELSE
 453             ETHRES = 0.D+00
 454             IF ( EKNNUC .GT. ETHRES ) THEN
 455                LASSOR = .FALSE.
 456             ELSE
 457                LASSOR = .TRUE.
 458             END IF
 459          END IF
 460          IF ( LASSOR .OR. LTRPPD ) THEN
 461             IF ( KPNUC .EQ. 1 .OR. KPNUC .EQ. 8 ) THEN
 462                KPNUCL (IKPMX) = - KPNUCL (IKPMX)
 463                ENNUC (IKPMX) = ENNUC (IKPMX) - EKFNUC (IKPMX)
 464                EKFNUC (IKPMX) = -AINFNT
 465                NPROT  = NPROT + ICH (KPNUC)
 466                NNEUT  = NNEUT + 1 - ICH (KPNUC)
 467                IREINT = IREINT + 1
 468                LPRCYC = .FALSE.
 469                GO TO 200
 470             ELSE
 471                LTRPPD = .TRUE.
 472                STOP 'KPNUCL_TRAPPED'
 473             END IF
 474          ELSE
 475             LTRPPD = .FALSE.
 476          END IF
 477          IKPNWI = IKPMX
 478          IF ( KPNUC .EQ. 1 .OR. KPNUC .EQ. 8 ) THEN
 479             ETHMNM = ETHRES + EKEMNM
 480             ETHREI = MAX ( EKREXP, ETHMNM )
 481             IF ( EKECON .LT. ETHMNM ) THEN
 482                NPROT = NPROT + ICH (KPNUC)
 483                NNEUT = NNEUT + 1 - ICH (KPNUC)
 484                ENNUC  (IKPMX) = EKECON + AM (KPNUC)
 485                KPNUCL (IKPMX) = - KPNUCL (IKPMX)
 486                NNUCTS = NNUCTS + 1
 487                INUCTS (NNUCTS) = IKPMX
 488                JNUCTS = NUSCIN
 489                EKFNUC (IKPMX) = EKFPRE
 490                ENNUC  (IKPMX) = ENNUC (IKPMX) - BNDUSE + BNDGEN
 491                RSTNUC (IKPMX) = BNDGEN
 492                LPRCYC = .TRUE.
 493                GO TO 200
 494             ELSE IF ( EKECON .LT. ETHREI ) THEN
 495                LBCHCK = .FALSE.
 496                IKPNWI = - IKPMX
 497             ELSE
 498                LBCHCK = .FALSE.
 499             END IF
 500          ELSE
 501             LBCHCK = .FALSE.
 502          END IF
 503          PNUCCO = SQRT ( EKECON * ( EKECON + 2.D+00 * AM (KPNUC) ) )
 504          CALL NWISEL ( IKPNWI, LBCHCK )
 505   350    CONTINUE
 506          IF ( BIMPCT .GT. RADTOT .AND. .NOT. LTRPPD ) THEN
 507             KPRIN = KPNUC
 508             IF ( EKNNUC .NE. EKECON ) PNUCCO = SQRT ( EKNNUC * ( EKNNUC
 509      &         + 2.D+00 * AM (KPNUC) ) )
 510             IF ( ABS ( PNUCL (IKPMX) - PNUCCO ) .GT. ANGLGB * PNUCCO )
 511      &           CALL PHDSET ( IKPMX )
 512             IBRES = IBRES - IBAR (KPNUC)
 513             ICRES = ICRES - ICH  (KPNUC)
 514             BBRES = IBRES
 515             ZZRES = ICRES
 516             AMMRES = AMMAFT
 517             AMNRES = AMNAFT
 518             CALL UMOFIN ( IKPMX, BBRES, ZZRES, LTRPPD )
 519             EKNNUC = ENNUC (IKPMX) - AM (KPNUC)
 520             IF ( LTRPPD ) GO TO 350
 521             NP = NP + 1
 522             TKI   (NP) = ENNUC  (IKPMX) - AM (KPNUC)
 523             KPART (NP) = KPNUC
 524             PLR   (NP) = PNUCL  (IKPMX)
 525             CXR   (NP) = PXNUCL (IKPMX) / PLR (NP)
 526             CYR   (NP) = PYNUCL (IKPMX) / PLR (NP)
 527             CZR   (NP) = PZNUCL (IKPMX) / PLR (NP)
 528             WEI   (NP) = WEE
 529             KPNUCL (IKPMX) = 0
 530             IF ( KPNUC .EQ. 1 .OR. KPNUC .EQ. 8 ) THEN
 531                IGREYP = IGREYP + ICH (KPNUC)
 532                IGREYN = IGREYN + 1 - ICH (KPNUC)
 533                PXINTR = PXINTR + PXNUCL (IKPMX)
 534                PYINTR = PYINTR + PYNUCL (IKPMX)
 535                PZINTR = PZINTR + PZNUCL (IKPMX)
 536                EINTR  = EINTR  + ENNUC  (IKPMX)
 537                IBINTR = IBINTR + IBAR   (KPNUC)
 538                ICINTR = ICINTR + ICH    (KPNUC)
 539             ELSE
 540                IOTHER = IOTHER + 1
 541                PXNUCR = PXNUCR + PXNUCL (IKPMX)
 542                PYNUCR = PYNUCR + PYNUCL (IKPMX)
 543                PZNUCR = PZNUCR + PZNUCL (IKPMX)
 544                ENUCR  = ENUCR  + ENNUC  (IKPMX)
 545                IBNUCR = IBNUCR + IBAR   (KPNUC)
 546                ICNUCR = ICNUCR + ICH    (KPNUC)
 547             END IF
 548             BNPREV = BNPREV + BNDUSE
 549             IF ( IREINT .LE. 0 ) THEN
 550                LPRCYC = .TRUE.
 551             ELSE
 552                LPRCYC = .FALSE.
 553             END IF
 554             GO TO 200
 555          ELSE IF ( BIMPCT .GT. RADTOT ) THEN
 556             KRFNUC (IKPMX) = KRFNUC (IKPMX) + 1
 557             CALL GRNDM(RNDM,1)
 558             SINTHE = RNDM ( 1 )
 559             COSTHE = SQRT ( 1.D+00 - SINTHE )
 560             SINTHE = SQRT ( SINTHE )
 561   400       CONTINUE
 562                CALL GRNDM(RNDM,2)
 563                RPHI1 = 2.D+00 * RNDM (1) - 1.D+00
 564                RPHI2 = 2.D+00 * RNDM (2) - 1.D+00
 565                RPHI12 = RPHI1 * RPHI1
 566                RPHI22 = RPHI2 * RPHI2
 567                RSQ = RPHI12 + RPHI22
 568             IF ( RSQ .GT. 1.D+00 ) GO TO 400
 569             SINPHI = 2.D+00 * RPHI1 * RPHI2 / RSQ
 570             COSPHI = ( RPHI12 - RPHI22 ) / RSQ
 571             SINT02 = CXIMPC**2 + CYIMPC**2
 572             IF ( SINT02 .LT. ANGLSQ ) THEN
 573                PXNUCL (IKPMX) = COSPHI * SINTHE * PNUCCO
 574                PYNUCL (IKPMX) = SINPHI * SINTHE * PNUCCO
 575                PZNUCL (IKPMX) = CZIMPC * COSTHE * PNUCCO
 576             ELSE
 577                SINTH0 = SQRT ( SINT02 )
 578                UPRIME = SINTHE * COSPHI
 579                VPRIME = SINTHE * SINPHI
 580                COSPH0 = CXIMPC / SINTH0
 581                SINPH0 = CYIMPC / SINTH0
 582                PXNUCL (IKPMX) = ( UPRIME * COSPH0 * CZIMPC - VPRIME
 583      &                        * SINPH0 + COSTHE * CXIMPC ) * PNUCCO
 584                PYNUCL (IKPMX) = ( UPRIME * SINPH0 * CZIMPC + VPRIME
 585      &                        * COSPH0 + COSTHE * CYIMPC ) * PNUCCO
 586                PZNUCL (IKPMX) = ( - UPRIME * SINTH0 + COSTHE * CZIMPC )
 587      &                        * PNUCCO
 588             END IF
 589             PNUCL  (IKPMX) = PNUCCO
 590             XSTNUC (IKPMX) = XIMPTR
 591             YSTNUC (IKPMX) = YIMPTR
 592             ZSTNUC (IKPMX) = ZIMPTR
 593             RSTNUC (IKPMX) = ABS (RIMPTR)
 594             ENNUC (IKPMX) = EKECON + AM (KPNUC)
 595             EKFNUC (IKPMX) = -AINFNT
 596             GO TO 200
 597          ELSE
 598             IF ( ( KPNUC .EQ. 1 .OR. KPNUC .EQ. 8 ) .AND. EKECON .LE.
 599      &           ETHREI ) THEN
 600                KPNUCL (IKPMX) = - KPNUCL (IKPMX)
 601                ENNUC (IKPMX) = EKECON + AM (KPNUC)
 602                NPROT = NPROT + ICH (KPNUC)
 603                NNEUT = NNEUT + 1 - ICH (KPNUC)
 604                IF ( .NOT. LBCHCK .AND. IKPNWI .GT. 0 ) THEN
 605                   LPRCYC = .TRUE.
 606                   NNUCTS = NNUCTS + 1
 607                   INUCTS (NNUCTS) = IKPMX
 608                   JNUCTS = NUSCIN
 609                   EKFNUC (IKPMX) = EKFPRE
 610                   ENNUC  (IKPMX) = ENNUC (IKPMX) - BNDUSE + BNDGEN
 611                   RSTNUC (IKPMX) = BNDGEN
 612                   LPRCYC = .TRUE.
 613                ELSE
 614                   EKFNUC (IKPMX) = -AINFNT
 615                   IREINT = IREINT + 1
 616                   IF ( EKFREI .LT. ANGLGB ) THEN
 617                      EKFREI = 0.5D+00 * ( EKFIMP + EKFPRO )
 618                      RHOREI = 0.5D+00 * ( RHOIMP + RHOIMT )
 619                   END IF
 620                   LPRCYC = .FALSE.
 621                END IF
 622                GO TO 200
 623             END IF
 624             LBIMPC = .FALSE.
 625             KPRIN  = KPNUC
 626             KPNUCL (IKPMX) = 0
 627             KRFLIN = KRFNUC (IKPMX)
 628             IF ( EKNNUC .NE. EKECON ) PNUCCO = SQRT ( EKNNUC * ( EKNNUC
 629      &         + 2.D+00 * AM (KPRIN) ) )
 630             CXIMPC = PXNUCL (IKPMX) / PNUCL (IKPMX)
 631             CYIMPC = PYNUCL (IKPMX) / PNUCL (IKPMX)
 632             CZIMPC = PZNUCL (IKPMX) / PNUCL (IKPMX)
 633             XSTNUC (IKPMX) = XIMPTR
 634             YSTNUC (IKPMX) = YIMPTR
 635             ZSTNUC (IKPMX) = ZIMPTR
 636             RSTNUC (IKPMX) = ABS (RIMPTR)
 637             GO TO 100
 638          END IF
 639       END IF
 640       NEXPEM = 0
 641       IF ( LGDHPR ) THEN
 642          LBCHCK = .TRUE.
 643          EKECON = EKE
 644          PNUCCO = PPROJ
 645          CALL BIMSEL ( KPROJ, TXX, TYY, TZZ, LBCHCK )
 646          LELSTC = .FALSE.
 647          RHOIMP = 0.5D+00 * ( RHOIMP + RHOIMT )
 648          EKFIMP = 0.5D+00 * ( EKFIMP + EKFPRO )
 649          RHOMEM = RHOIMP
 650          EKFMEM = EKFIMP
 651       END IF
 652 *
 653       ANOW  = BBRES
 654       ZNOW  = ZZRES
 655       PXRES = PXORI
 656       PYRES = PYORI
 657       PZRES = PZORI
 658       PTRES = PTORI
 659       ERES  = EKE  + AM (KPROJ) + AMNTAR
 660       IF ( LGDHPR ) THEN
 661          IF ( KPROJ .EQ. 1 ) THEN
 662             NPROT = NPROT + 1
 663          ELSE IF ( KPROJ .EQ. 8 ) THEN
 664             NNEUT  = NNEUT + 1
 665          END IF
 666          ACOLL = BBTAR - 1.D+00
 667          IF ( KNUCIM .EQ. 1 ) THEN
 668             NPROT = NPROT + 1
 669             ZCOLL = ZZTAR - 1.D+00
 670          ELSE
 671             NNEUT = NNEUT + 1
 672             ZCOLL = ZZTAR
 673          END IF
 674          NHOLE = NHOLE + 1
 675       ELSE
 676          ACOLL = BBTAR - 1.D+00
 677          IF ( KPROJ .EQ. 1 ) THEN
 678             NPROT = NPROT + 1
 679             PRPONP = ZNOW / ( 3.D+00 * ANOW - 2.D+00 * ZNOW )
 680             CALL GRNDM(RNDM,1)
 681             IF ( RNDM (1) .LT. PRPONP ) THEN
 682                NPROT = NPROT + 1
 683                ZCOLL = ZZTAR - 1.D+00
 684                IPRTYP = KPROJ * 10 + 1
 685             ELSE
 686                NNEUT = NNEUT + 1
 687                ZCOLL = ZZTAR
 688                IPRTYP = KPROJ * 10 + 8
 689             END IF
 690             NHOLE = NHOLE + 1
 691          ELSE IF ( KPROJ .EQ. 8 ) THEN
 692             NNEUT  = NNEUT + 1
 693             PRNONP = 3.D+00 * ZNOW / ( 2.D+00 * ZNOW + ANOW )
 694             CALL GRNDM(RNDM,1)
 695             IF ( RNDM (1) .LT. PRNONP ) THEN
 696                NPROT = NPROT + 1
 697                ZCOLL = ZZTAR  - 1.D+00
 698                IPRTYP = KPROJ * 10 + 1
 699             ELSE
 700                NNEUT = NNEUT + 1
 701                ZCOLL = ZZTAR
 702                IPRTYP = KPROJ * 10 + 8
 703             END IF
 704             NHOLE = NHOLE + 1
 705          ELSE
 706             STOP 'KPROJ'
 707          END IF
 708       END IF
 709       ICYCL  = 0
 710       LEMISS = .FALSE.
 711   500 CONTINUE
 712       CALL PREPRE ( WEE, NNEUT, NPROT, NHOLE, ICYCL )
 713   530 CONTINUE
 714       IF ( IBRES .GT. 0 ) THEN
 715          EKR0   = ERES - AMNRES
 716          ATTNUM = ELBNDE (ICHTAR) - ELBNDE (ICRES) - EKR0 * ( AMMRES
 717      &          - AMNRES ) / AMMRES
 718          ERES   = ERES + AMMTAR - AMNTAR - ( ZZTAR - ZNOW ) * AMELEC
 719      &          + ATTNUM
 720          EKRES  = ERES - AMMRES
 721       ELSE
 722          AMMRES = 0.D+00
 723          AMNRES = 0.D+00
 724          ERES   = 0.D+00
 725          EKR0   = 0.D+00
 726          EKRES  = 0.D+00
 727          TVTENT = 0.D+00
 728          GO TO 600
 729       END IF
 730       IF ( EKRES .LE. 0.D+00 ) THEN
 731          WRITE ( LUNERR,* )' Peanut: negative kinetic energy for',
 732      &                     ' the residual nucleus!!',ICRES,IBRES,
 733      &                       REAL (EKRES)
 734          IF ( EKRES .LT. -3.D-3 ) THEN
 735             LRESMP = .TRUE.
 736             RETURN
 737          END IF
 738          EKRES  = 0.D+00
 739          TVRECL = 0.D+00
 740          AMSTAR = AMMRES
 741          TVCMS  = 0.D+00
 742          PTRES2 = 0.D+00
 743          PXRES  = 0.D+00
 744          PYRES  = 0.D+00
 745          PZRES  = 0.D+00
 746          PTRES  = 0.D+00
 747       ELSE
 748          PTRES2 = PTRES * PTRES
 749          AMSTAR = ( ERES - PTRES ) * ( ERES + PTRES )
 750          IF ( AMSTAR .GE. AMMRES**2 ) THEN
 751             AMSTAR = SQRT ( AMSTAR )
 752             TVCMS  = AMSTAR - AMMRES
 753          ELSE IF ( AMMRES**2 - AMSTAR .LT. 2.D+00 * AMSTAR * TVEPSI
 754      &             ) THEN
 755             AMSTAR = AMMRES
 756             ERES   = SQRT ( AMSTAR**2 + PTRES**2 )
 757             TVCMS  = 0.D+00
 758          ELSE IF ( AMSTAR .LE. 0.D+00 ) THEN
 759             WRITE ( LUNERR,* )' Peanut: immaginary mass for',
 760      &                        ' the residual nucleus!!',ICRES,IBRES,
 761      &                          REAL (AMSTAR)
 762             LRESMP = .TRUE.
 763             RETURN
 764          ELSE
 765             AMSTAR = SQRT ( AMSTAR )
 766             IF ( AMMRES - AMSTAR .LT. TVEPSI ) THEN
 767                AMSTAR = AMMRES
 768                TVCMS  = 0.D+00
 769                TVRECL = ERES - AMSTAR
 770                GO TO 550
 771             END IF
 772             WRITE ( LUNERR,* )' Peanut: negative excitation energy for',
 773      &                        ' the residual nucleus!!',ICRES,IBRES,
 774      &                          REAL ( AMSTAR - AMMRES )
 775             IF ( AMSTAR - AMMRES .LT. -3.D-3 ) THEN
 776                LRESMP = .TRUE.
 777                RETURN
 778             END IF
 779             AMSTAR = AMMRES
 780             TVCMS  = 0.D+00
 781             HELP   = SQRT ( ( ERES - AMMRES ) * ( ERES + AMMRES ) )
 782      &             / PTRES
 783             PXRES = PXRES * HELP
 784             PYRES = PYRES * HELP
 785             PZRES = PZRES * HELP
 786             PTRES = PTRES * HELP
 787          END IF
 788          TVRECL = ERES - AMSTAR
 789       END IF
 790   550 CONTINUE
 791       IF ( TVRECL .LT. 0.D+00 ) THEN
 792          TVRECL = 0.D+00
 793       END IF
 794       TV     = 0.D+00
 795       EKRES  = TVRECL
 796   600 CONTINUE
 797       EOTEST = EOTEST - EINTR - ENUCR - EKR0 - AMNRES
 798       IF ( ABS (EOTEST) .GT. ETEPS ) THEN
 799          WRITE (LUNERR,*)' Peanut: eotest failure',EOTEST
 800          LRESMP = .TRUE.
 801          RETURN
 802       END IF
 803       IF ( IBRES .EQ. 0 ) RETURN
 804       EOTEST = ETTOT + AMMTAR - AMNTAR + ATTNUM
 805       IF ( KPROJ .EQ. 1 ) THEN
 806          EOTEST = EOTEST + AMHEAV (2) - AM (1)
 807       ELSE IF ( KPROJ .EQ. 8 ) THEN
 808          EOTEST = EOTEST + AMHEAV (1) - AM (8)
 809       ELSE
 810          EOTEST = EOTEST + ICH(KPROJ) * AMELEC
 811       END IF
 812       IF ( LEVPRT ) THEN
 813          CALL EVEVAP ( WEE )
 814          IF ( LRESMP ) RETURN
 815       ELSE
 816          TVHEAV = 0.D+00
 817       END IF
 818       DO 1000 KP = NP0+1,NP
 819          IF ( KPART (KP) .EQ. 1 ) THEN
 820             EOTEST = EOTEST - TKI (KP) - AMHEAV (2)
 821             IBTOT  = IBTOT  - 1
 822             ICHTOT = ICHTOT - 1
 823          ELSE IF ( KPART (KP) .EQ. 8 ) THEN
 824             EOTEST = EOTEST - TKI (KP) - AMHEAV (1)
 825             IBTOT  = IBTOT  - 1
 826          ELSE
 827             EOTEST = EOTEST - TKI (KP) - AM (KPART(KP))
 828             IBTOT  = IBTOT  - IBAR (KPART(KP))
 829             ICHTOT = ICHTOT - ICH  (KPART(KP))
 830          END IF
 831  1000 CONTINUE
 832       EOTEST = EOTEST - TVHEAV - IEVDEU * AMHEAV (3)
 833      &       - IEVTRI * AMHEAV (4)
 834      &       - IEV3HE * AMHEAV (5)
 835      &       - IEV4HE * AMHEAV (6)
 836      &       - AMMRES - TVRECL
 837       IBTOT  = IBTOT  - IEVDEU * 2 - IEVTRI * 3 - IEV3HE * 3
 838      &       - IEV4HE * 4
 839       ICHTOT = ICHTOT - IEVDEU - IEVTRI - IEV3HE * 2
 840      &       - IEV4HE * 2
 841       IF ( LRNFSS ) THEN
 842          IF ( LHEAVY ) THEN
 843             DO 2000 JP = 1, NPHEAV
 844                IF ( KHEAVY (JP) .GT. 6 ) THEN
 845                   EOTEST = EOTEST - AMHEAV (JP)
 846                   IBTOT  = IBTOT  - IBHEAV (KHEAVY(JP))
 847                   ICHTOT = ICHTOT - ICHEAV (KHEAVY(JP))
 848                END IF
 849  2000       CONTINUE
 850          ELSE
 851             DO 2100 JFISS = 1, NFISS
 852                IBHLP = NINT (ATFIS(JFISS))
 853                IF ( IBHLP .GT. 0 ) THEN
 854                   ICHLP  = NINT (ZTFIS(JFISS))
 855                   EOTEST = EOTEST - 1.D-03 * AMTFIS (JFISS)
 856                   IBTOT  = IBTOT  - IBHLP
 857                   ICHTOT = ICHTOT - ICHLP
 858                END IF
 859  2100       CONTINUE
 860          END IF
 861       END IF
 862       IF ( ABS (EOTEST) .GT. 1.D+3 * ETEPS ) THEN
 863          WRITE (LUNERR,*)
 864      &   ' Peanut failure!!, Eotest,Ammres,Tvrecl,Ibres,Icres',
 865      &                       EOTEST,AMMRES,TVRECL,IBRES,ICRES
 866       END IF
 867       IF ( IBTOT .NE. IBRES .OR. ICHTOT .NE. ICRES ) THEN
 868          WRITE (LUNERR,*)
 869      &   ' Peanut failure!!, Ichtot, Icres, Ibtot, Ibres',
 870      &                       ICHTOT, ICRES, IBTOT, IBRES
 871       END IF
 872 *=== End of subroutine peanut =========================================*
 873       RETURN
 874       END