GEANT321/peanut/bimsel.F

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