GEANT321/peanut/bimsel.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1999/05/18 15:55:22  fca
   6 * AliRoot sources
   7 *
   8 * Revision 1.2  1997/10/17 10:00:03  mclareni
   9 * Remove SAVE statement for NT
  10 *
  11 * Revision 1.1.1.1  1995/10/24 10:22:00  cernlib
  12 * Geant
  13 *
  14 *
  15 #include "geant321/pilot.h"
  16 *CMZ :  3.21/02 29/03/94  15.41.45  by  S.Giani
  17 *-- Author :
  18 *$ CREATE BIMSEL.FOR
  19 *COPY BIMSEL
  20 *
  21 *=== bimsel ===========================================================*
  22 *
  23       SUBROUTINE BIMSEL ( JPROJ, TXX, TYY, TZZ, LBCHCK )
  24
  25 #include "geant321/dblprc.inc"
  26 #include "geant321/dimpar.inc"
  27 #include "geant321/iounit.inc"
  28 *
  29 *----------------------------------------------------------------------*
  30 *----------------------------------------------------------------------*
  31 *
  32 #include "geant321/balanc.inc"
  33 #include "geant321/eva0.inc"
  34 #include "geant321/nucdat.inc"
  35 #include "geant321/nucgeo.inc"
  36 #include "geant321/part.inc"
  37 #include "geant321/parevt.inc"
  38 #include "geant321/resnuc.inc"
  39 *
  40       PARAMETER ( FEFFEC = 1.518066780142162 D+00 )
  41       PARAMETER ( BETMAX = 0.4 D+00 )
  42 *
  43       REAL RNDM(2)
  44 *
  45       LOGICAL LBCHCK, LFERMI, LLMDBR
  46 *
  47       SAVE ABTAR , ZZTAR , SIGMP0, SIGMN0,
  48      &     AMNHLP, RHOBIM, RPRONU, RADPRP, RADPRN, DSKRED, RHRUSF,
  49      &     AUSFL , ZUSFL , BNDSAV, RADHLP, BFCHLP, BIMCLM, PRCOLP,
  50      &     PRCOLN, IBTOLD, ICTOLD, KPROJ , NTRIAL, ITFRMI
  51       SAVE SBHAL0, SBSKI0, SBCEN0, SBCEN1, SBSKI1, SBHAL1
  52 #ifndef CERNLIB_WINNT
  53 *      SAVE
  54 #endif
  55       DATA IBTOLD, ICTOLD / 2*0 /
  56 *
  57       KPROJ  = JPROJ
  58       AUSFL  = IBTAR
  59       ZUSFL  = ICHTAR
  60       RHRUSF = 1.D+00
  61       BEPROJ = PNUCCO / ( EKECON + AM (KPROJ) )
  62       CXIMPC = TXX
  63       CYIMPC = TYY
  64       CZIMPC = TZZ
  65       NTRIAL = 0
  66       RHOBIM = - AINFNT
  67       IF ( KPROJ .EQ. 1 .OR. KPROJ .EQ. 8 ) THEN
  68          IPWELL = 1 + KPROJ / 8
  69          WLLRED = 1.D+00
  70          BNDNUC = BNENRG (IPWELL)
  71       ELSE
  72          IPWELL = 0
  73          IF ( IBAR (KPROJ) .EQ. 0 ) THEN
  74             IF ( KPROJ .LE. 11 ) THEN
  75                WLLRED = 0.D+00
  76                BNDNUC = 0.D+00
  77             ELSE
  78                WLLRED = POTMES
  79                BNDNUC = BNENRG (3)
  80             END IF
  81          ELSE
  82             WLLRED = POTBAR
  83             BNDNUC = BNENRG (3)
  84          END IF
  85       END IF
  86       BNDSAV = BNDNUC
  87       IF ( IBAR (KPROJ) .NE. 0 ) THEN
  88          RPRONU = 1.D+00
  89       ELSE IF ( KPROJ .NE. 7 ) THEN
  90          RPRONU = 0.8164965809277260D+00
  91       ELSE
  92          RPRONU = 0.D+00
  93       END IF
  94       IF ( LBCHCK ) THEN
  95          LFERMI = .FALSE.
  96          EKESIG = EKECON
  97          PPRSIG = PNUCCO
  98          CALL SIGFER ( KPTOIP (KPROJ), EKESIG, PPRSIG, LFERMI )
  99          PRCOLP = ZUSFL / AUSFL * SIGMAP
 100          PRCOLN = ( AUSFL - ZUSFL ) / AUSFL * SIGMAN
 101          SIGMAA = PRCOLP + PRCOLN
 102          PRCOLP = PRCOLP / SIGMAA
 103          PRCOLN = 1.D+00 - PRCOLP
 104       END IF
 105       IF ( RPRONU .GT. ANGLGB ) THEN
 106          IF ( LPARWV ) THEN
 107             LLMDBR = .TRUE.
 108             TMP102 = 1.D-02
 109             PDEBRO = MAX ( PNUCCO, TMP102 )
 110             ALMBAR = PLABRC / PDEBRO
 111             DEBRLM = 0.5D+00 * ALMBAR
 112             RADCOR = SQRT ( (RPRONU * RMSPRO)**2 + ALMBAR**2 )
 113      &             / ( RMSPRO * RPRONU )
 114          ELSE
 115             PDEBRO = ( EKECON + BNDNUC ) * ( EKECON + BNDNUC + 2.D+00
 116      &             * AM (KPROJ) )
 117             LLMDBR = .FALSE.
 118             DEBRLM = 0.D+00
 119             ALMBAR = 0.D+00
 120             LLLMAX = -1
 121             RADCOR = SQRT ( (RPRONU * RMSPRO)**2 + PLABRC**2 / PDEBRO
 122      &             ) / ( RMSPRO * RPRONU )
 123          END IF
 124       ELSE
 125          RADCOR = 0.D+00
 126          LLMDBR = .FALSE.
 127          DEBRLM = 0.D+00
 128       END IF
 129       RADCO2 = RADCOR
 130       RADPRO = MIN ( TWOTWO * RMSPRO * RPRONU * RADCOR, SKGT16 )
 131       RADPRP = RADPRO
 132       RADPRN = RADPRO
 133       IF ( IBTAR .NE. IBTOLD .OR. ICHTAR .NE. ICTOLD ) THEN
 134          IBTOLD = IBTAR
 135          ICTOLD = ICHTAR
 136          ABTAR  = IBTAR
 137          ZZTAR  = ICHTAR
 138          AR1O3  = RMASS (IBTAR)
 139          AMNHLP = 0.5D+00 * ( AMNUCL (1) + AMNUCL (2) )
 140          HKAP   = ABTAR**2 / ( ZZTAR**2 + ( ABTAR - ZZTAR )**2 )
 141          HHLP (1) = ( HKAP * ZZTAR )**0.3333333333333333D+00 / AR1O3
 142          HHLP (2) = ( HKAP * ( ABTAR - ZZTAR ) )
 143      &            **0.3333333333333333D+00 / AR1O3
 144          RHOCEN = RHOTAB (IBTAR)
 145          ALPHAL = ALPTAB (IBTAR)
 146          RADIU0 = RADTAB (IBTAR)
 147          SKINDP = SKITAB (IBTAR)
 148          HALODP = HALTAB (IBTAR)
 149          RADIU1 = RADIU0 + SKINDP
 150          RADTOT = RADIU1 + HALODP
 151          RHOCOR = ONEMNS * RHOCEN
 152          RHOSKN = ALPHAL * RHOCEN
 153          PFRCEN (1) = HHLP (1) * PFRTAB (IBTAR)
 154          PFRCEN (2) = HHLP (2) * PFRTAB (IBTAR)
 155          RHOAVE = RHATAB (IBTAR)
 156          PFRAVE = PFATAB (IBTAR)
 157          EKFAVE = EKATAB (IBTAR)
 158          OMALHL = 1.D+00 - ALPHAL
 159          RAD1O2 = RADIU0 + 0.5D+00 * SKINDP / OMALHL
 160          SKNEFF = SKINDP / OMALHL
 161          RADSKN = RADIU0 + SKNEFF
 162          EKFCEN (1) = SQRT ( AMNUSQ (1) + PFRCEN (1)**2 ) - AMNUCL (1)
 163          EKFCEN (2) = SQRT ( AMNUSQ (2) + PFRCEN (2)**2 ) - AMNUCL (2)
 164          IF ( PFRCEN (1) .GT. PFRCEN (2) ) THEN
 165             ITNCMX = 1
 166          ELSE
 167             ITNCMX = 2
 168          END IF
 169          CALL NCLVST ( IBTAR, ICHTAR )
 170       END IF
 171       IF ( IPWELL .LE. 0 ) IPWELL = ITNCMX
 172       CALL NCLVIN
 173       IF ( EKECON .LT. 2.D+00 * GAMMIN ) THEN
 174          EKECON = 0.D+00
 175          PNUCCO = 0.D+00
 176          LABRST = .TRUE.
 177          RADPRO = 0.D+00
 178          RADSIG = ( RADTOT + DEBRLM ) * BFCLMB
 179          RADMAX = RADTOT
 180          LLLMAX = -1
 181          OPACTY = 2.D+00
 182          CALL RSTSEL (KPROJ)
 183          RETURN
 184       END IF
 185       RADMAX = RADTOT + RADPRO
 186       BIMCLM = RDCLMB * BFCLMB
 187       IF ( LLMDBR ) THEN
 188          RADHLP = RADMAX
 189          IF ( RADHLP .LE. RDCLMB ) THEN
 190             BFCMAX = BFCLMB
 191             BFCHLP = 0.5D+00 * CLMBBR * RDCLMB / EKECON
 192          ELSE
 193             BFCHLP = 0.5D+00 * CLMBBR * RDCLMB / EKECON
 194             BFCMAX = SQRT ( 1.D+00 - CLMBBR * RDCLMB / EKECON / RADHLP )
 195          END IF
 196          BIMMAX = RADHLP * BFCMAX
 197          LLLMAX = INT ( BIMMAX / ALMBAR )
 198          RADSIG = ALMBAR * ( LLLMAX + 1.D+00 )
 199          SIGGEO = PI * RADSIG * RADSIG
 200       ELSE
 201          RADHLP = RADTOT + RADPRO + DEBRLM
 202          IF ( RADHLP .LE. RDCLMB ) THEN
 203             BFCMAX = BFCLMB
 204          ELSE
 205             BFCHLP = 0.5D+00 * CLMBBR * RDCLMB / EKECON
 206             BFCMAX = SQRT ( 1.D+00 - CLMBBR * RDCLMB / EKECON / RADHLP )
 207          END IF
 208          RADSIG = RADHLP * BFCMAX
 209       END IF
 210       R0TRAJ = - RADTOT
 211       R1TRAJ = - R0TRAJ
 212  4200  CONTINUE
 213           CALL GRNDM(RNDM,2)
 214           RPHI1 = 2.D+00 * RNDM (1) - 1.D+00
 215           RPHI2 = 2.D+00 * RNDM (2) - 1.D+00
 216           RPHI12 = RPHI1 * RPHI1
 217           RPHI22 = RPHI2 * RPHI2
 218           RSQ = RPHI12 + RPHI22
 219        IF ( RSQ .GT. 1.D+00 ) GO TO 4200
 220       SINPHI = 2.D+00 * RPHI1 * RPHI2 / RSQ
 221       COSPHI = ( RPHI12 - RPHI22 ) / RSQ
 222       SINT02 = 1.D+00 - CZIMPC * CZIMPC
 223       IF ( SINT02 .LT. ANGLSQ ) THEN
 224          UBIMPC = COSPHI
 225          VBIMPC = SINPHI
 226          WBIMPC = 0.D+00
 227       ELSE
 228          SINTH0 = SQRT ( SINT02 )
 229          SINPH0 = CYIMPC / SINTH0
 230          COSPH0 = CXIMPC / SINTH0
 231          UBIMPC = COSPHI * COSPH0 * CZIMPC - SINPHI * SINPH0
 232          VBIMPC = COSPHI * SINPH0 * CZIMPC + SINPHI * COSPH0
 233          WBIMPC = - COSPHI * SINTH0
 234       END IF
 235       GO TO 4500
 236       ENTRY BIMNXT ( LBCHCK )
 237          IF ( EKECON .LT. 2.D+00 * GAMMIN ) THEN
 238             LABRST = .TRUE.
 239             CALL RSTNXT
 240             RETURN
 241          END IF
 242  4300 CONTINUE
 243          BNDNUC = BNDSAV
 244          SIGMAP = SIGMP0
 245          SIGMAN = SIGMN0
 246  4400    CONTINUE
 247          CALL GRNDM(RNDM,1)
 248          ANMFP = - LOG ( 1.D+00 - RNDM (1) ) / DSKRED
 249          IF ( SBRES * SIGMAA .GT. ANMFP ) THEN
 250             GO TO 6000
 251          END IF
 252  4500 CONTINUE
 253  5000 CONTINUE
 254          SBUSED = 0.D+00
 255          NTRIAL = NTRIAL + 1
 256          IF ( LLMDBR ) THEN
 257             ALLMAX = LLLMAX + 1.D+00
 258             CALL GRNDM(RNDM,2)
 259             RNDLLL = ALLMAX * MAX ( RNDM (1), RNDM (2) )
 260             LLLACT = INT (RNDLLL)
 261             BIMPTR = RNDLLL * ALMBAR
 262             BIMPTR = ABS (BIMPTR)
 263             IF ( BIMPTR .LE. BIMCLM ) THEN
 264                BFCEFF = BFCLMB
 265             ELSE
 266                HLPHLP = BFCHLP / BIMPTR
 267                BFCEFF = 1.D+00 / ( HLPHLP + SQRT ( HLPHLP * HLPHLP
 268      &                + 1.D+00 ) )
 269             END IF
 270             BIMPTR = BIMPTR / BFCEFF
 271             IF ( BIMPTR .GT. RADHLP ) GO TO 5000
 272          ELSE
 273             CALL GRNDM(RNDM,2)
 274             BIMPTR = RADSIG * MAX ( RNDM (1), RNDM (2) )
 275             IF ( BIMPTR .LE. BIMCLM ) THEN
 276                BFCEFF = BFCLMB
 277             ELSE
 278                HLPHLP = BFCHLP / BIMPTR
 279                BFCEFF = 1.D+00 / ( HLPHLP + SQRT ( HLPHLP * HLPHLP
 280      &                + 1.D+00 ) )
 281             END IF
 282             BIMPTR = BIMPTR / BFCEFF
 283          END IF
 284          BIMMEM = BIMPTR
 285          IF ( BIMPTR .GT. RADTOT - RADPRO ) THEN
 286             BIMPCT = 0.5D+00 * ( RADTOT + BIMPTR - RADPRO )
 287             IF ( BIMPTR .GE. RADTOT ) THEN
 288                X1 = BIMPTR - RADTOT
 289                ANGRED = ACOS ( 2.D+00 * X1 / ( RADPRO + X1 ) ) / PI
 290                X1 = X1 / ( R0PROT * RPRONU * RADCO2 )
 291                DSKRED = ( 0.5D+00 * X1 * X1 + X1 + 1.D+00 ) * EXP (-X1)
 292      &                * ANGRED
 293             ELSE
 294                X1 = RADPRO + BIMPTR - RADTOT
 295                ANGRED = ACOS ( 2.D+00 * X1 / ( RADPRO + X1 ) ) / PI
 296                X1 = X1 / ( R0PROT * RPRONU * RADCO2 )
 297                DSKRED = 1.D+00 - ( 0.5D+00 * X1 * X1 + X1 + 1.D+00 )
 298      &                * EXP (-X1) * ANGRED
 299             END IF
 300          ELSE
 301             DSKRED = 1.D+00
 302             BIMPCT = BIMPTR
 303          END IF
 304          IF ( .NOT. LBCHCK ) THEN
 305             RHOSAV = RHOBIM
 306             RHOBIM = FRHONC ( BIMPCT )
 307             IF ( RHOBIM .EQ. RHOSAV ) GO TO 5500
 308             PFRBIM = FPFRNC ( RHOBIM, ITNCMX )
 309             EKFBIM = FEKFNC ( PFRBIM, ITNCMX )
 310             RHOHLP = FRHONC ( BIMPTR )
 311             PFRHLP = FPFRNC ( RHOHLP, IPWELL )
 312             PFRHLP = 0.5D+00 * PFRHLP * PFRHLP / AMNUSQ (IPWELL)
 313             IF ( BIMPTR .GT. RADTOT ) BNDNUC = BNDNUC * ( 1.D+00
 314      &         - ( BIMPTR - RADTOT ) / ( RADHLP - RADTOT ) )
 315             VPRBIM = WLLRED * ( AMNUCL (IPWELL) * PFRHLP
 316      &             * ( 1.D+00 - 0.5D+00 * PFRHLP ) + BNDNUC )
 317             LFERMI = .TRUE.
 318             EKESIG = EKECON
 319             PPRSIG = PNUCCO
 320             CALL SIGFER ( KPTOIP (KPROJ), EKESIG, PPRSIG, LFERMI )
 321             PRCOLP = ZUSFL / AUSFL * SIGMAP
 322             PRCOLN = ( AUSFL - ZUSFL ) / AUSFL * SIGMAN
 323             SIGMAA = PRCOLP + PRCOLN
 324             PRCOLP = PRCOLP / SIGMAA
 325             PRCOLN = 1.D+00 - PRCOLP
 326  5500       CONTINUE
 327          ELSE
 328             RHOBIM = FRHONC ( BIMPCT )
 329          END IF
 330          XBIMPC = UBIMPC * BIMPCT
 331          YBIMPC = VBIMPC * BIMPCT
 332          ZBIMPC = WBIMPC * BIMPCT
 333          CALL GRNDM(RNDM,1)
 334          ANMFP  = - LOG ( 1.D+00 - RNDM (1) ) / DSKRED
 335          IF ( BIMPCT .GT. RAD1O2 ) THEN
 336             SBTTSQ = 4.D+00 * ( RADTOT**2 - BIMPCT**2 ) * RHOBIM**2
 337             IF ( SBTTSQ .LE. ( ANMFP / SIGMAA )**2 ) GO TO 5000
 338          END IF
 339          CALL SBCOMP ( SBHAL0, SBSKI0, SBCEN0, SBCEN1, SBSKI1, SBHAL1 )
 340          SBTOT  = SBHAL0 + SBSKI0 + SBCEN0 + SBCEN1 + SBSKI1 + SBHAL1
 341          SBTOT  = RHRUSF * SBTOT
 342       IF ( SBTOT * SIGMAA .LE. ANMFP ) GO TO 5000
 343  6000 CONTINUE
 344       SBUSED = SBUSED * RHRUSF + ANMFP / SIGMAA
 345       SBRES  = SBTOT  - SBUSED
 346       SBUSED = SBUSED / RHRUSF
 347       LELSTC = .TRUE.
 348       NTARGT = 1
 349       CALL GRNDM(RNDM,1)
 350       IF ( RNDM (1) .LT. PRCOLP ) THEN
 351          KNUCIM = 1
 352          ITFRMI = 1
 353       ELSE
 354          KNUCIM = 8
 355          ITFRMI = 2
 356       END IF
 357       IPRTYP = KPROJ * 10 + KNUCIM
 358       CALL RSCOMP ( SBHAL0, SBSKI0, SBCEN0, SBCEN1, SBSKI1, SBHAL1 )
 359       BNDNUC = BNDSAV
 360       RHOIMP = FRHONC ( ABS (RIMPCT) )
 361       PFRIMP = FPFRNC ( RHOIMP, ITFRMI )
 362       EKFIMP = FEKFNC ( PFRIMP, ITFRMI )
 363       RIMHLP = ABS (RIMPTR)
 364       RHOIMT = FRHONC ( RIMHLP )
 365       PFRPRO = FPFRNC ( RHOIMT, IPWELL )
 366       EKFPRO = FEKFNC ( PFRPRO, IPWELL )
 367       IF ( RIMHLP .GT. RADTOT ) BNDNUC = BNDNUC * ( 1.D+00 - ( RIMHLP
 368      &                                 - RADTOT ) / ( RADHLP - RADTOT ))
 369       VPRWLL = WLLRED * ( EKFPRO + BNDNUC )
 370       EKEWLL = EKECON + VPRWLL
 371       EPSWLL = EKEWLL + AM (KPROJ)
 372       IF ( .NOT. LBCHCK ) THEN
 373          PPRWLL = SQRT ( EKEWLL * ( EKEWLL + 2.D+00 * AM (KPROJ) ) )
 374          CALL PHDWLL ( UBIMPC, VBIMPC, WBIMPC )
 375          PNFRMI = PFNCLV ( ITFRMI, .TRUE. )
 376          IF ( PNFRMI .LT. -100.D+00 ) GO TO 4400
 377          CALL RACO ( PXFERM, PYFERM, PZFERM )
 378          PXFERM = PXFERM * PNFRMI
 379          PYFERM = PYFERM * PNFRMI
 380          PZFERM = PZFERM * PNFRMI
 381          ERES   = EKEWLL + AM (KPROJ) + AM (KNUCIM) + EKFERM
 382          PXRES  = PXPROJ + PXFERM
 383          PYRES  = PYPROJ + PYFERM
 384          PZRES  = PZPROJ + PZFERM
 385          PTRES2 = PXRES**2 + PYRES**2 + PZRES**2
 386          UMO2   = ERES**2 - PTRES2
 387          EKESIG = 0.5D+00 * ( UMO2 - AM (KPROJ)**2 - AM (KNUCIM)**2 )
 388      &          / AM (KNUCIM) - AM (KPROJ)
 389          EKFIMP = MAX ( EKFERM, EKFIMP )
 390       ELSE
 391          EKESIG = EPSWLL - AM (KPROJ)
 392       END IF
 393       PPRSIG = SQRT ( EKESIG * ( EKESIG + 2.D+00 * AM (KPROJ) ) )
 394       SIGMN0 = SIGMAN
 395       SIGMP0 = SIGMAP
 396       LFERMI = .FALSE.
 397       CALL SIGFER ( KPTOIP (KPROJ), EKESIG, PPRSIG, LFERMI )
 398       IF ( KNUCIM .EQ. 1 ) THEN
 399          SIGMAR = SIGMAP / SIGMP0
 400       ELSE
 401          SIGMAR = SIGMAN / SIGMN0
 402       END IF
 403       SIGMAR = MIN ( SIGMAR, ONEONE )
 404       CALL GRNDM(RNDM,1)
 405       RNDREJ = RNDM(1)
 406       IF ( RNDREJ .GE. SIGMAR ) GO TO 4300
 407       IF ( LBCHCK ) THEN
 408          ZITA   = 0.5D+00 * ( EKFIMP + EKFPRO ) / EKEWLL
 409          IF ( ZITA .LE. 0.5D+00 ) THEN
 410             PZITA = 1.D+00 - 1.4D+00 * ZITA
 411          ELSE
 412             PZITA = 1.D+00 - 1.4D+00 * ZITA + 0.4D+00 * ZITA
 413      &            * ( 2.D+00 - 1.D+00 / ZITA )**2.5D+00
 414          END IF
 415          RNDREJ = RNDREJ / SIGMAR
 416          IF ( RNDREJ .GE. PZITA ) GO TO 4300
 417       ELSE
 418          PZITA = 1.D+00
 419       END IF
 420       OPACTY = 1.D+00 / NTRIAL
 421       RETURN
 422 *=== End of subroutine Bimsel =========================================*
 423       END