GEANT321/peanut/bimsel.F

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