GEANT321/peanut/nwisel.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:22:01  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 NWISEL.FOR
  13 *COPY NWISEL
  14 *
  15 *=== nwisel ===========================================================*
  16 *
  17       SUBROUTINE NWISEL ( IKN, LBCHCK )
  18
  19 #include "geant321/dblprc.inc"
  20 #include "geant321/dimpar.inc"
  21 #include "geant321/iounit.inc"
  22 *
  23 *----------------------------------------------------------------------*
  24 *----------------------------------------------------------------------*
  25 *
  26 #include "geant321/balanc.inc"
  27 #include "geant321/eva0.inc"
  28 #include "geant321/nucdat.inc"
  29 #include "geant321/nucgeo.inc"
  30 #include "geant321/parevt.inc"
  31 #include "geant321/parnuc.inc"
  32 #include "geant321/part.inc"
  33 #include "geant321/resnuc.inc"
  34 *
  35       PARAMETER ( FEFFEC = 1.518066780142162 D+00 )
  36       PARAMETER ( BETMAX = 0.4 D+00 )
  37 *
  38       REAL RNDM(1), RNGAUS, DUMNOR
  39 *
  40       LOGICAL LBCHCK, LFERMI, LFRMCK, LLMDBR
  41 *
  42       SAVE XBIMTR, YBIMTR, ZBIMTR, COSTHE, SINTHE, SIGMP0, EKENWI,
  43      &     SIGMN0, RHOBIM, RPRONU, RADPRP, RADPRN, DSKRED, RHRUSF,
  44      &     AUSFL , ZUSFL , PRCOLP, PRCOLN, KPROJ , IKPMX , LFRMCK
  45       SAVE SBHAL0, SBSKI0, SBCEN0, SBCEN1, SBSKI1, SBHAL1
  46 *
  47       LFRMCK = IKN .LT. 0
  48       IKPMX  = ABS (IKN)
  49       KPROJ  = KPNUCL (IKPMX)
  50       EKENWI = EKECON
  51       AUSFL  = IBTAR
  52       ZUSFL  = ICHTAR
  53       RHRUSF = 1.D+00
  54       BEPROJ = PNUCCO / ( EKENWI + AM (KPROJ) )
  55       RHOBIM = - AINFNT
  56 *
  57       IF ( KPROJ .EQ. 1 .OR. KPROJ .EQ. 8 ) THEN
  58          IPWELL = 1 + KPROJ / 8
  59          WLLRED = 1.D+00
  60       ELSE
  61          IPWELL = ITNCMX
  62          IF ( IBAR (KPROJ) .EQ. 0 ) THEN
  63             IF ( KPROJ .LE. 11 ) THEN
  64                WLLRED = 0.D+00
  65             ELSE
  66                WLLRED = POTMES
  67             END IF
  68          ELSE
  69             WLLRED = POTBAR
  70          END IF
  71       END IF
  72       IF ( IBAR (KPROJ) .NE. 0 ) THEN
  73          RPRONU = 1.D+00
  74       ELSE IF ( KPROJ .NE. 7 ) THEN
  75          RPRONU = 0.8164965809277260D+00
  76       ELSE
  77          RPRONU = 0.D+00
  78       END IF
  79       IF ( LBCHCK ) THEN
  80          LFERMI = .FALSE.
  81          EKESIG = EKENWI
  82          PPRSIG = PNUCCO
  83          CALL SIGFER ( KPTOIP (KPROJ), EKESIG, PPRSIG, LFERMI )
  84          PRCOLP = ZUSFL / AUSFL * SIGMAP
  85          PRCOLN = ( AUSFL - ZUSFL ) / AUSFL * SIGMAN
  86          SIGMAA = PRCOLP + PRCOLN
  87          PRCOLP = PRCOLP / SIGMAA
  88          PRCOLN = 1.D+00 - PRCOLP
  89       END IF
  90       IF ( RPRONU .GT. ANGLGB ) THEN
  91          IF ( LPARWV ) THEN
  92             LLMDBR = .TRUE.
  93             TMP102 = 1.D-02
  94             PDEBRO = MAX ( PNUCL (IKPMX), TMP102 )
  95             DEBRLM = 0.5D+00 * PLABRC / PDEBRO
  96             RADCOR = 1.D+00
  97          ELSE
  98             LLMDBR = .FALSE.
  99             DEBRLM = 0.D+00
 100             RADCOR = 1.D+00
 101          END IF
 102       ELSE
 103          RADCOR = 0.D+00
 104          LLMDBR = .FALSE.
 105          DEBRLM = 0.D+00
 106       END IF
 107       RADCO2 = RADCOR
 108       RADPRO = MIN ( TWOTWO * RMSPRO * RPRONU * RADCOR, SKGT16 )
 109       RADPRP = RADPRO
 110       RADPRN = RADPRO
 111       CXIMPC = PXNUCL (IKPMX) / PNUCL (IKPMX)
 112       CYIMPC = PYNUCL (IKPMX) / PNUCL (IKPMX)
 113       CZIMPC = PZNUCL (IKPMX) / PNUCL (IKPMX)
 114       COSTHE = ( CXIMPC * XSTNUC (IKPMX) + CYIMPC * YSTNUC (IKPMX)
 115      &         + CZIMPC * ZSTNUC (IKPMX) ) / RSTNUC (IKPMX)
 116       SINTHE = SQRT ( 1.D+00 - MIN ( COSTHE * COSTHE, ONEONE ) )
 117       BIMPTR = RSTNUC (IKPMX) * SINTHE
 118       IF ( LLMDBR ) THEN
 119          BIMOLD = BIMPTR
 120          CALL GRNDM(RNDM,1)
 121          BIMPTR = BIMPTR + ( 2.D+00 * RNDM (1) - 1.D+00 )
 122      &          * DEBRLM
 123          CALL GRANOR ( RNGAUS, DUMNOR)
 124          HEISEX = DEBRLM * RNGAUS
 125          DSTMIN = - RSTNUC (IKPMX) * COSTHE
 126          XBIMTR = XSTNUC (IKPMX) + DSTMIN * CXIMPC
 127          YBIMTR = YSTNUC (IKPMX) + DSTMIN * CYIMPC
 128          ZBIMTR = ZSTNUC (IKPMX) + DSTMIN * CZIMPC
 129          XPARAL = XSTNUC (IKPMX) - XBIMTR
 130          YPARAL = YSTNUC (IKPMX) - YBIMTR
 131          ZPARAL = ZSTNUC (IKPMX) - ZBIMTR
 132          IF ( BIMOLD .GT. ANGLGB ) THEN
 133             BIMOLD = BIMPTR / BIMOLD
 134             XBIMTR = XBIMTR * BIMOLD
 135             YBIMTR = YBIMTR * BIMOLD
 136             ZBIMTR = ZBIMTR * BIMOLD
 137          ELSE
 138             CALL RACO ( CXHLP, CYHLP, CZHLP )
 139             XBIMTR = BIMPTR * ( CYHLP * CZIMPC - CZHLP * CYIMPC )
 140             YBIMTR = BIMPTR * ( CZHLP * CXIMPC - CXHLP * CZIMPC )
 141             ZBIMTR = BIMPTR * ( CXHLP * CYIMPC - CYHLP * CXIMPC )
 142          END IF
 143          XSTNUC (IKPMX) = XPARAL + XBIMTR + HEISEX * CXIMPC
 144          YSTNUC (IKPMX) = YPARAL + YBIMTR + HEISEX * CYIMPC
 145          ZSTNUC (IKPMX) = ZPARAL + ZBIMTR + HEISEX * CZIMPC
 146          RSTNUC (IKPMX) = SQRT ( XSTNUC (IKPMX)**2 + YSTNUC (IKPMX)**2
 147      &                  + ZSTNUC (IKPMX)**2 )
 148          COSTHE = ( CXIMPC * XSTNUC (IKPMX) + CYIMPC * YSTNUC (IKPMX)
 149      &          + CZIMPC * ZSTNUC (IKPMX) ) / RSTNUC (IKPMX)
 150          BIMPTR = ABS (BIMPTR)
 151          SINTHE = BIMPTR / MAX ( RSTNUC (IKPMX), ANGLGB )
 152       ELSE
 153          DSTMIN = - RSTNUC (IKPMX) * COSTHE
 154          XBIMTR = XSTNUC (IKPMX) + DSTMIN * CXIMPC
 155          YBIMTR = YSTNUC (IKPMX) + DSTMIN * CYIMPC
 156          ZBIMTR = ZSTNUC (IKPMX) + DSTMIN * CZIMPC
 157       END IF
 158       IF ( COSTHE .GE. 0.D+00 .AND. RSTNUC (IKPMX) .GE. RADTOT
 159      &   + RADPRO ) GO TO 4500
 160 *
 161       SBUSED = 0.D+00
 162       IF ( BIMPTR .GT. RADTOT - RADPRO ) THEN
 163          BIMPCT = 0.5D+00 * ( RADTOT + BIMPTR - RADPRO )
 164          TRUFAC = BIMPCT / BIMPTR
 165          IF ( BIMPTR .GE. RADTOT ) THEN
 166             X1 = BIMPTR - RADTOT
 167             IF ( X1 .LE. RADPRO ) THEN
 168                ANGRED = ACOS ( 2.D+00 * X1 / ( RADPRO + X1 ) ) / PI
 169                X1 = X1 / ( R0PROT * RPRONU * RADCO2 )
 170                DSKRED = ( 0.5D+00 * X1 * X1 + X1 + 1.D+00 )
 171      &                * EXP (-X1) * ANGRED
 172             END IF
 173          ELSE
 174             X1 = RADPRO + BIMPTR - RADTOT
 175             ANGRED = ACOS ( 2.D+00 * X1 / ( RADPRO + X1 ) ) / PI
 176             X1 = X1 / ( R0PROT * RPRONU * RADCO2 )
 177             DSKRED = 1.D+00 - ( 0.5D+00 * X1 * X1 + X1 + 1.D+00 )
 178      &             * EXP (-X1) * ANGRED
 179          END IF
 180       ELSE
 181          BIMPCT = BIMPTR
 182          TRUFAC = 1.D+00
 183          DSKRED = 1.D+00
 184       END IF
 185       IF ( BIMPCT .LT. RADTOT ) THEN
 186          IF ( .NOT. LBCHCK ) THEN
 187             RHOSAV = RHOBIM
 188             RHOBIM = FRHONC ( BIMPCT )
 189             IF ( RHOBIM .EQ. RHOSAV ) GO TO 5500
 190             PFRBIM = FPFRNC ( RHOBIM, ITNCMX )
 191             EKFBIM = FEKFNC ( PFRBIM, ITNCMX )
 192             RHOHLP = FRHONC ( BIMPTR )
 193             PFRHLP = FPFRNC ( RHOHLP, IPWELL )
 194             PFRHLP = 0.5D+00 * PFRHLP * PFRHLP / AMNUSQ (IPWELL)
 195             VPRBIM = WLLRED * ( AMNUCL (IPWELL) * PFRHLP
 196      &             * ( 1.D+00 - 0.5D+00 * PFRHLP ) + BNDNUC )
 197             LFERMI = .TRUE.
 198             EKESIG = EKENWI
 199             PPRSIG = PNUCCO
 200             CALL SIGFER ( KPTOIP (KPROJ), EKESIG, PPRSIG, LFERMI )
 201             PRCOLP = ZUSFL / AUSFL * SIGMAP
 202             PRCOLN = ( AUSFL - ZUSFL ) / AUSFL * SIGMAN
 203             SIGMAA = PRCOLP + PRCOLN
 204             IF ( SIGMAA .GT. 0.D+00 ) THEN
 205                PRCOLP = PRCOLP / SIGMAA
 206                PRCOLN = 1.D+00 - PRCOLP
 207             ELSE
 208                PRCOLP = 0.D+00
 209                PRCOLN = 0.D+00
 210                XBIMPC = XBIMTR * TRUFAC
 211                YBIMPC = YBIMTR * TRUFAC
 212                ZBIMPC = ZBIMTR * TRUFAC
 213                GO TO 4500
 214             END IF
 215  5500       CONTINUE
 216          END IF
 217          XBIMPC = XBIMTR * TRUFAC
 218          YBIMPC = YBIMTR * TRUFAC
 219          ZBIMPC = ZBIMTR * TRUFAC
 220          IF ( COSTHE .GE. 0.D+00 ) THEN
 221             R0TRAJ = RSTNUC (IKPMX) * TRUFAC
 222             R1TRAJ = RADTOT
 223             IF ( R0TRAJ .GT. RADTOT ) GO TO 4500
 224          ELSE
 225             R0TRAJ = - MIN ( RSTNUC (IKPMX) * TRUFAC, RADTOT )
 226             R1TRAJ = RADTOT
 227          END IF
 228          CALL GRNDM(RNDM,1)
 229          ANMFP  = - LOG ( 1.D+00 - RNDM (1) ) / DSKRED
 230          IF ( BIMPCT .GT. RAD1O2 ) THEN
 231             SBTTSQ = 4.D+00 * ( RADTOT**2 - BIMPCT**2 ) * RHOBIM**2
 232             IF ( SBTTSQ .LE. ( ANMFP / SIGMAA )**2 ) GO TO 4500
 233          END IF
 234          CALL SBCOMP ( SBHAL0, SBSKI0, SBCEN0, SBCEN1, SBSKI1, SBHAL1 )
 235          SBTOT  = SBHAL0 + SBSKI0 + SBCEN0 + SBCEN1 + SBSKI1 + SBHAL1
 236          SBTOT  = RHRUSF * SBTOT
 237          IF ( SBTOT * SIGMAA .GT. ANMFP ) GO TO 5000
 238       END IF
 239  4500 CONTINUE
 240       BIMPCT = RADTOT + RADPRO + 0.0001D+00
 241       IF ( BIMPTR .LT. BIMPCT ) THEN
 242          DSTMIN = SQRT ( BIMPCT**2 - BIMPTR**2 )
 243          RIMPTR = BIMPCT
 244          XIMPTR = XBIMTR + CXIMPC * DSTMIN
 245          YIMPTR = YBIMTR + CYIMPC * DSTMIN
 246          ZIMPTR = ZBIMTR + CZIMPC * DSTMIN
 247       ELSE
 248          XIMPTR = XBIMTR
 249          YIMPTR = YBIMTR
 250          ZIMPTR = ZBIMTR
 251          RIMPTR = BIMPTR
 252       END IF
 253       CXIMPC = -XIMPTR / RIMPTR
 254       CYIMPC = -YIMPTR / RIMPTR
 255       CZIMPC = -ZIMPTR / RIMPTR
 256       RETURN
 257       ENTRY NWINXT ( LBCHCK )
 258  4800 CONTINUE
 259          SIGMAP = SIGMP0
 260          SIGMAN = SIGMN0
 261  4900    CONTINUE
 262            CALL GRNDM(RNDM,1)
 263            ANMFP = - LOG ( 1.D+00 - RNDM (1) ) / DSKRED
 264          IF ( SBRES * SIGMAA .LE. ANMFP ) GO TO 4500
 265  5000 CONTINUE
 266       SBUSED = SBUSED * RHRUSF + ANMFP / SIGMAA
 267       SBRES  = SBTOT  - SBUSED
 268       SBUSED = SBUSED / RHRUSF
 269       CALL GRNDM(RNDM,1)
 270       IF ( RNDM (1) .LT. PRCOLP ) THEN
 271          KNUCIM = 1
 272          ITFRMI = 1
 273       ELSE
 274          KNUCIM = 8
 275          ITFRMI = 2
 276       END IF
 277       CALL RSCOMP ( SBHAL0, SBSKI0, SBCEN0, SBCEN1, SBSKI1, SBHAL1 )
 278       RHOIMP = FRHONC ( ABS (RIMPCT) )
 279       PFRIMP = FPFRNC ( RHOIMP, ITFRMI )
 280       EKFIMP = FEKFNC ( PFRIMP, ITFRMI )
 281       RHOIMT = FRHONC ( ABS (RIMPTR) )
 282       PFRPRO = FPFRNC ( RHOIMT, IPWELL )
 283       EKFPRO = FEKFNC ( PFRPRO, IPWELL )
 284       VPRWLL = WLLRED * ( EKFPRO + BNDNUC )
 285       LELSTC = .TRUE.
 286       NTARGT = 1
 287       EKEWLL = EKENWI + VPRWLL
 288       EPSWLL = EKEWLL + AM (KPROJ)
 289       IF ( .NOT. LBCHCK ) THEN
 290          PPRWLL = SQRT ( EKEWLL * ( EKEWLL + 2.D+00 * AM (KPROJ) ) )
 291          PXPROJ = PPRWLL * CXIMPC
 292          PYPROJ = PPRWLL * CYIMPC
 293          PZPROJ = PPRWLL * CZIMPC
 294          PNFRMI = PFNCLV ( ITFRMI, .TRUE. )
 295          IF ( PNFRMI .LT. - 100.D+00 ) GO TO 4900
 296          CALL RACO ( PXFERM, PYFERM, PZFERM )
 297          PXFERM = PXFERM * PNFRMI
 298          PYFERM = PYFERM * PNFRMI
 299          PZFERM = PZFERM * PNFRMI
 300          EKFERM = SQRT ( PNFRMI**2 + AM (KNUCIM)**2 ) - AM (KNUCIM)
 301          ERES   = EKEWLL + AM (KPROJ) + AM (KNUCIM) + EKFERM
 302          PXRES  = PXPROJ + PXFERM
 303          PYRES  = PYPROJ + PYFERM
 304          PZRES  = PZPROJ + PZFERM
 305          PTRES2 = PXRES**2 + PYRES**2 + PZRES**2
 306          UMO2   = ERES**2 - PTRES2
 307          EKESIG = 0.5D+00 * ( UMO2 - AM (KPROJ)**2 - AM (KNUCIM)**2 )
 308      &          / AM (KNUCIM) - AM (KPROJ)
 309          EKFIMP = MAX ( EKFERM, EKFIMP )
 310       ELSE
 311          EKESIG = EKEWLL
 312       END IF
 313       PPRSIG = SQRT ( EKESIG * ( EKESIG + 2.D+00 * AM (KPROJ) ) )
 314       SIGMN0 = SIGMAN
 315       SIGMP0 = SIGMAP
 316       LFERMI = .FALSE.
 317       CALL SIGFER ( KPTOIP (KPROJ), EKESIG, PPRSIG, LFERMI )
 318       IF ( KNUCIM .EQ. 1 ) THEN
 319          SIGMAR = SIGMAP / SIGMP0
 320       ELSE
 321          SIGMAR = SIGMAN / SIGMN0
 322       END IF
 323       SIGMAR = MIN ( SIGMAR, ONEONE )
 324       CALL GRNDM(RNDM,1)
 325       RNDREJ = RNDM (1)
 326       IF ( RNDREJ .GE. SIGMAR ) GO TO 4800
 327       IF ( LBCHCK ) THEN
 328          ZITA   = 0.5D+00 * ( EKFIMP + EKFPRO ) / EKEWLL
 329          IF ( ZITA .LE. 0.5D+00 ) THEN
 330             PZITA = 1.D+00 - 1.4D+00 * ZITA
 331          ELSE
 332             PZITA = 1.D+00 - 1.4D+00 * ZITA + 0.4D+00 * ZITA
 333      &            * ( 2.D+00 - 1.D+00 / ZITA )**2.5D+00
 334          END IF
 335          RNDREJ = RNDREJ / SIGMAR
 336          IF ( RNDREJ .GE. PZITA ) GO TO 4800
 337       ELSE IF ( LFRMCK ) THEN
 338          EKFPRJ = EKFPRO + DEFNUC (IPWELL)
 339          EKFTAR = EKFIMP + DEFNUC (ITFRMI)
 340          IF ( EKEWLL + EKFERM .LE. EKFTAR + EKFPRJ ) GO TO 4800
 341          PCMSSQ = ( PXPROJ - 0.5D+00 * PXRES )**2
 342      &          + ( PYPROJ - 0.5D+00 * PYRES )**2
 343      &          + ( PZPROJ - 0.5D+00 * PZRES )**2
 344          RNDREJ = RNDREJ / SIGMAR
 345          PPROSQ = EKFPRJ * ( EKFPRJ + 2.D+00 * AM (KPROJ) )
 346          PTARSQ = EKFTAR * ( EKFTAR + 2.D+00 * AMNUCL (ITFRMI) )
 347          HLPHLP = PCMSSQ + 0.25D+00 * PTRES2
 348          HLPHL2 = PCMSSQ * PTRES2
 349          IF ( RNDREJ .GE. 0.5D+00 ) THEN
 350             RNDREJ = 2.D+00 * ( RNDREJ - 0.5D+00 )
 351             IF ( HLPHLP .LT. PTARSQ ) THEN
 352                COSQMX = 0.D+00
 353             ELSE
 354                COSQMX = ( HLPHLP - PTARSQ )**2 / HLPHL2
 355             END IF
 356             IF ( HLPHLP .GT. PPROSQ ) THEN
 357                COSQMN = 0.D+00
 358             ELSE
 359                COSQMN = ( PPROSQ - HLPHLP )**2 / HLPHL2
 360             END IF
 361          ELSE
 362             RNDREJ = 2.D+00 * RNDREJ
 363             IF ( HLPHLP .LT. PPROSQ ) THEN
 364                COSQMX = 0.D+00
 365             ELSE
 366                COSQMX = ( HLPHLP - PPROSQ )**2 / HLPHL2
 367             END IF
 368             IF ( HLPHLP .GT. PTARSQ ) THEN
 369                COSQMN = 0.D+00
 370             ELSE
 371                COSQMN = ( PTARSQ - HLPHLP )**2 / HLPHL2
 372             END IF
 373             COSQMX = ( 0.25D+00 * PTRES2 + PCMSSQ - PPROSQ )**2
 374      &             / ( PTRES2 * PCMSSQ )
 375          END IF
 376          RNDREJ = RNDREJ**2
 377          IF ( RNDREJ .GE. COSQMX .OR. RNDREJ .LT. COSQMN ) GO TO 4800
 378       END IF
 379       RETURN
 380 *=== End of subroutine Nwisel =========================================*
 381       END