GEANT321/fluka/corrin.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:19:54  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.42  by  S.Giani
  11 *-- Author :
  12 *$ CREATE CORRIN.FOR
  13 *COPY CORRIN
  14 *                                                                      *
  15 *=== corrin ===========================================================*
  16 *                                                                      *
  17       SUBROUTINE CORRIN ( ZZTAR, BBTAR, KPROJ, PPROJ, EKE )
  18
  19 #include "geant321/dblprc.inc"
  20 #include "geant321/dimpar.inc"
  21 #include "geant321/iounit.inc"
  22 *
  23 *----------------------------------------------------------------------*
  24 *                                                                      *
  25 *     Created on  10  may  1990    by    Alfredo Ferrari & Paola Sala  *
  26 *                                                   Infn - Milan       *
  27 *                                                                      *
  28 *     Last change on 23-mar-93     by    Alfredo Ferrari               *
  29 *                                                                      *
  30 *    This version has been developed by A. Ferrari starting from the   *
  31 *    the one by J. M. Zazula: it includes a few differences mainly in  *
  32 *    the correlation since now we try to get the excitation energy     *
  33 *    from a correct energy balance                                     *
  34 *----------------------------------------------------------------------*
  35 *
  36 #include "geant321/balanc.inc"
  37 #include "geant321/corinc.inc"
  38 #include "geant321/nucdat.inc"
  39 #include "geant321/parevt.inc"
  40 #include "geant321/part.inc"
  41 #include "geant321/resnuc.inc"
  42       PARAMETER ( ALFMAX = 5.00D+00 )
  43       PARAMETER ( FINCFR = 1.00D+00 )
  44       PARAMETER ( TMNOLD = 0.02D+00 )
  45       COMMON / FKNUCO / HELP (2), HHLP (2), FTVTH (2), FINCX (2),
  46      &                  EKPOLD (2), BBOLD, ZZOLD, SQROLD, ASEASQ,
  47      &                  FSPRED, FEX0RD
  48       DIMENSION KPA (26)
  49       DIMENSION FRA(2,110), ESLOLD (2), AVMULT (2), EXDUMM (2),
  50      &          V0WSAV (2), VEWSAV (2)
  51       REAL RNDM(3)
  52       LOGICAL LSTOP, LNUCNW, LLLIN, LLPOW
  53       SAVE FRA, FRAC, FRAC0, FRAMAX, KPA, IJJ, LSTOP, V0WSAV, VEWSAV
  54       DATA KPA/1,1,3,3,3,3,3,2,2,3,3,3,3,3,3,3,2,2,3,1,1,2,3,3,3,3/
  55 *
  56 *  Reduction factors for intran. cascade energy, taken from Alsmiller
  57 *  Incoming baryons
  58       DATA (FRA(1,I), I=1,107) / .048D0,.076D0,0.10D0,.12D0,
  59      * .14D0,.16D0,.17D0,.19D0,.21D0,.22D0,.24D0,
  60      1 .25D0,.26D0,.28D0,.29D0,.30D0,.315D0,.33D0,.34D0,
  61      * .35D0,.36D0,.38D0,.39D0,.40D0,.41D0,.42D0,
  62      2 .43D0,.44D0,.46D0,.47D0,.48D0,.49D0,.50D0,.51D0,
  63      * .52D0,.53D0,.54D0,.55D0,.56D0,.57D0,.58D0,.59D0,
  64      3 .60D0,.61D0,.62D0,.63D0,.635D0,.64D0,.65D0,.66D0,
  65      * .67D0,.675D0,.68D0,.69D0,.70D0,.71D0,.715D0,
  66      4 .72D0,.725D0,.73D0,.735D0,.74D0,.75D0,.755D0,
  67      * .76D0,.767D0,.77D0,.77D0,.775D0,.78D0,.783D0,
  68      5 .786D0,.79D0,.795D0,.80D0,.805D0,.81D0,.812D0,
  69      * .815D0,.82D0,.822D0,.824D0,.825D0,.825D0,.83D0,
  70      6 .832D0,.834D0,.836D0,.838D0,.84D0,.843D0,.836D0,
  71      * .849D0,.852D0,.855D0,.856D0,.857D0,.858D0,
  72      7 .859D0,.860D0,.862D0,.864D0,.866D0,.868D0,.870D0,
  73      * .872D0,.874D0/
  74 *  Incoming mesons
  75       DATA (FRA(2,I), I=1,63) / .048D0,.076D0,0.10D0,.12D0,
  76      * .14D0,.16D0,.17D0,.19D0,.21D0,.22D0,.24D0,
  77      1 .25D0,.262D0,.274D0,.285D0,.295D0,.305D0,.315D0,
  78      * .327D0,.340D0,.345D0,.350D0,.360D0,
  79      2 .370D0,.380D0,.385D0,.390D0,.398D0,.406D0,.415D0,
  80      * .420D0,.425D0,.430D0,.435D0,.440D0,.446D0,
  81      3 .452D0,.458D0,.464D0,.470D0,.474D0,.478D0,.482D0,
  82      * .486D0,.490D0,.495D0,.500D0,.505D0,.510D0,
  83      4 .515D0,.519D0,.523D0,.526D0,.520D0,.533D0,.537D0,
  84      * .540D0,.543D0,.547D0,.550D0,.555D0,.559D0,
  85      5 .5625D0 /
  86 *
  87       IBPROJ = IBAR (KPROJ)
  88       IJJ    = KPA  (KPROJ)
  89       ANUAV  = 1.D+00
  90       ANCOLL = 1.D+00
  91       ANUSEA = 1.D+00
  92       NSEA   = 0
  93       LSTOP = .FALSE.
  94 *  +-------------------------------------------------------------------*
  95 *  |                                        Incoming baryons
  96       IF ( IBPROJ .NE. 0 ) THEN
  97 *  |  +----------------------------------------------------------------*
  98 *  |  |
  99          IF ( IBTAR .LE. 107 ) THEN
 100             FRAC = FRA ( 1, IBTAR )
 101 *  |  |
 102 *  |  +----------------------------------------------------------------*
 103 *  |  |
 104          ELSE IF (IBTAR .LE. 206) THEN
 105             FRAC = 0.739D+00 + 0.00126D+00 * BBTAR
 106 *  |  |
 107 *  |  +----------------------------------------------------------------*
 108 *  |  |
 109          ELSE
 110             FRAC = 1.D+00
 111          END IF
 112 *  |  |
 113 *  |  +----------------------------------------------------------------*
 114 *  |
 115 *  +-------------------------------------------------------------------*
 116 *  |                                         Incoming mesons
 117       ELSE
 118 *  |  +----------------------------------------------------------------*
 119 *  |  |
 120          IF ( IBTAR .LE. 63 ) THEN
 121             FRAC = FRA ( 2, IBTAR ) * ( 1.D+00 + 0.1333D+00 * MAX (
 122      &             0.D+00, BBTAR - 35.D+00 ) / 28.D+00 )
 123 *  |  |
 124 *  |  +----------------------------------------------------------------*
 125 *  |  |
 126          ELSE IF ( IBTAR .LE. 107 ) THEN
 127             FRAC = 0.85D+00 * FRA ( 1, IBTAR )
 128 *  |  |
 129 *  |  +----------------------------------------------------------------*
 130 *  |  |
 131          ELSE IF ( IBTAR .LE. 206 ) THEN
 132             FRAC = 0.6279D+00 + 0.001077D+00 * BBTAR
 133 *  |  |
 134 *  |  +----------------------------------------------------------------*
 135 *  |  |
 136          ELSE
 137             FRAC = 0.85D+00
 138          END IF
 139 *  |  |
 140 *  |  +----------------------------------------------------------------*
 141       END IF
 142 *  |
 143 *  +-------------------------------------------------------------------*
 144 *  +-------------------------------------------------------------------*
 145 *  |
 146       IF ( BBTAR .NE. BBOLD .OR. ZZTAR .NE. ZZOLD ) THEN
 147          LNUCNW = .TRUE.
 148 *  |  Supply the fraction of the total kinetic energy to be
 149 *  |  used for intranuclear cascade nucleons
 150          SQRAMS = SQRT ( BBTAR )
 151          ATO1O3 = BBTAR**0.3333333333333333D+00
 152 *        ZTO1O3 = ZZTAR**0.3333333333333333D+00
 153          HKAP   = BBTAR**2 / ( ZZTAR**2 + ( BBTAR - ZZTAR )**2 )
 154          HHLP (1) = ( HKAP * ZZTAR )**0.3333333333333333D+00 / ATO1O3
 155          HHLP (2) = ( HKAP * ( BBTAR - ZZTAR ) )
 156      &              **0.3333333333333333D+00 / ATO1O3
 157          RDSNUC = R0NUCL * ATO1O3
 158          RDSCOU = RCCOUL * ATO1O3
 159          FLKCOU = DOST ( 1, ZZTAR )
 160          VEFFNU (1) = COULPR * FLKCOU * ZZTAR / RDSCOU
 161          VEFFNU (2) = 0.D+00
 162          AMRCAV = MAX ( 1.D+00, ( BBTAR - 2.D+00 ) ) * AMUC12
 163          AMRCSQ = AMRCAV * AMRCAV
 164 *  |  +----------------------------------------------------------------*
 165 *  |  |
 166          DO 3000 I = 1, 2
 167             PFRMMX (I) = HHLP (I) * APFRMX
 168             P2HELP     = PFRMMX (I)**2
 169             EFRMMX (I) = SQRT ( P2HELP + AMNUSQ (I) ) - AMNUCL (I)
 170             EFRMAV (I) = 0.3D+00 * P2HELP / AMNUCL (I) * ( 1.D+00
 171      &                 - P2HELP / ( 5.6D+00 * AMNUSQ (I) ) )
 172             ERCLAV (I) = 0.3D+00 * P2HELP / AMRCAV * ( 1.D+00
 173      &                 - P2HELP / ( 5.6D+00 * AMRCSQ ) )
 174             V0WELL (I) = EFRMMX (I) + EBNDNG (I)
 175             VEFFNU (I) = VEFFNU (I) + V0WELL (I)
 176             V0WSAV (I) = V0WELL (I)
 177             VEWSAV (I) = VEFFNU (I)
 178             ERCLMX     = 0.5D+00 * P2HELP / AMRCAV * ( 1.D+00
 179      &                 - 0.25D+00 * P2HELP / AMRCSQ )
 180             EKMNNU (I) = VEFFNU (I) + EBNDNG (I)
 181      &                 - EFRMMX (I) + ERCLMX
 182             EKMXNU (I) = VEFFNU (I) + EBNDNG (I)
 183             EKMNAV (I) = VEFFNU (I) + EBNDNG (I)
 184      &                 - EFRMAV (I) + ERCLAV (I)
 185             ESWELL (I) = EBNDNG (I) + V0WELL (I)
 186      &                 - EFRMAV (I) + ERCLAV (I)
 187             FTVTH  (I) = ( EKMNNU (I) / EFRMMX (I) )**3
 188             FTVTH  (I) = 0.4D+00 * SQRT  ( FTVTH (I) )
 189 3000     CONTINUE
 190 *  |  |
 191 *  |  +----------------------------------------------------------------*
 192          BBOLD  = BBTAR
 193          ZZOLD  = ZZTAR
 194          SQROLD = SQRAMS
 195 *  |
 196 *  +-------------------------------------------------------------------*
 197 *  |
 198       ELSE
 199          V0WSAV (1) = V0WELL (1)
 200          VEWSAV (1) = VEFFNU (1)
 201          V0WSAV (2) = V0WELL (2)
 202          VEWSAV (2) = VEFFNU (2)
 203          LNUCNW = .FALSE.
 204       END IF
 205 *  |
 206 *  +-------------------------------------------------------------------*
 207       FRAC0 = FRAC
 208       RETURN
 209       ENTRY CORSTP ( EKEFF )
 210          LSTOP = .TRUE.
 211       ENTRY CORRNC ( EKEFF )
 212       WEIGH1 = MIN ( 1.D+00, EKEFF / 4.D+00 )
 213       FRAC   = ( WEIGH1 * 2.5D+00 + ( 1.D+00 - WEIGH1 ) ) * FRAC0
 214       FRAMAX = 2.5D+00
 215 *  The following is a reduction factor to bring Hannes's parametri-
 216 *  zations for the average energy carried by cascade nucleons in
 217 *  better agreement with experimental data
 218       FSPRED = FSPRD0 + ( 1.D+00 - FSPRD0 ) * MAX ( 8.D+00 - EKEFF,
 219      &         0.D+00 ) / 8.D+00
 220       FSPRED = FSPRED * MIN ( EKEFF / 0.8D+00, 1.D+00 )
 221       TMPFSP = 1.25D+00 * FSPRD0
 222       FSPRED = MIN ( FSPRED, TMPFSP )
 223       AVMULT (1) = FRAC * BNKEKA ( 1, EKEFF, BBOLD, SQROLD )
 224       AVMULT (2) = FRAC * BNKEKA ( 2, EKEFF, BBOLD, SQROLD )
 225       EKUPNU (1) = BEKEKA ( 2, EKEFF, BBOLD, SQROLD )
 226       EKUPNU (2) = BEKEKA ( 3, EKEFF, BBOLD, SQROLD )
 227       EINCP  = FRAC * EKUPNU (1)
 228       EINCN  = FRAC * EKUPNU (2)
 229       EKMAX  = EKEFF - EBNDAV
 230       EKUPN0 = MAX ( EKUPNU (1), EKUPNU (2) )
 231 *  +-------------------------------------------------------------------*
 232 *  |
 233       DO 4000 I = 1, 2
 234          FINCUP (I) = BKEKA ( I, EKEFF, BBOLD )
 235          RATRAT  =  MIN ( 1.D+00, 0.5D+00 * EKMAX / FINCUP (I) )
 236          TMPFIN  = 0.5D+00 * EKMAX
 237          FINCUP (I) = MIN ( FINCUP (I), TMPFIN )
 238          EKPOLD (I) = EKUPNU (I) * FRAC * RATRAT
 239          EKUPNU (I) = FINCUP (I)
 240          TMPEKU     = 1.5D+00 * EKMXNU (I)
 241          EKUPNU (I) = MAX ( FRAMAX * EKUPNU (I), TMPEKU )
 242          TMPEKU     = 0.7D+00 * EKEFF
 243          EKUPNU (I) = MIN ( EKUPNU (I), TMPEKU )
 244          ESLOLD (I) = MAX ( FINCUP (I), TMNOLD )
 245          FINCX  (I) = ( ESLOLD (I) + ESWELL (I) ) / ESLOLD (I) * FSPRED
 246          ESLOPE (I) = MIN ( ESLOLD (I) * FINCX (I), EKMAX )
 247          AHELP      = EXP ( - EKPOLD (I) / ESLOLD (I) )
 248          EKNOLD     = ESLOLD (I) - EKPOLD (I) * AHELP /
 249      &              ( 1.D+00 - AHELP )
 250          EXMNAV (I) = EXP ( - EKMNAV (I) / ESLOPE (I) )
 251          EXMNNU (I) = EXP ( - EKMNNU (I) / ESLOPE (I) )
 252          EXUPNU (I) = EXP ( - EKUPNU (I) / ESLOPE (I) )
 253          EKINAV (I) = ESLOPE (I) + ( EKMNAV (I) * EXMNAV (I) -
 254      &                EKUPNU (I) * EXUPNU (I) ) / ( EXMNAV (I) -
 255      &                EXUPNU (I) )
 256          FINCUP (I) = AVMULT (I) * EKINAV (I) / EKPOLD (I)
 257 *        FINCUP (I) = EKINAV (I) / EKNOLD
 258 4000  CONTINUE
 259 *  |
 260 *  +-------------------------------------------------------------------*
 261       EINCP = EKPOLD (1)
 262       EINCN = EKPOLD (2)
 263       IF ( LSTOP ) RETURN
 264 *  Power sampling
 265       LLPOW  = .TRUE.
 266       LLLIN  = .FALSE.
 267 *  Take into account that a fraction 1/A^2/3 is lost because of
 268 *  peripheral collisions
 269       PERCOR = ATO1O3 * ATO1O3
 270       PERCOR = PERCOR / ( PERCOR - 1.D+00 )
 271       EINCP0 = EINCP  * FINCUP (1) * PERCOR
 272       EINCN0 = EINCN  * FINCUP (2) * PERCOR
 273       EINCT  = EINCP0 + EINCN0
 274 * For the other distr,:
 275       EINCMX = EKEFF - EBNDAV - ( AV0WEL - AEFRMA ) / ( 1.D+00 +
 276      &         4.D+00 * ( AV0WEL - AEFRMA ) / EKEFF )
 277       EIUSE  = 0.5D+00 * ( EINCMX + EKMAX )
 278       AHNORM = 0.D+00
 279       EINCM0 = - AINFNT
 280 *  +-------------------------------------------------------------------*
 281 *  |  Compute alfa:
 282       IF ( EINCMX .GT. 2.1D+00 * EINCT ) THEN
 283          EINFIN = MIN ( 0.95D+00 * EINCMX, ( ALFMAX + 2.D+00 ) * EINCT )
 284          ALFA   = EINFIN / EINCT - 2.D+00
 285 *  |
 286 *  +-------------------------------------------------------------------*
 287 *  |  Energy interval too small, sample E uniformly
 288       ELSE
 289          LLPOW  = .FALSE.
 290       END IF
 291 *  |
 292 *  +-------------------------------------------------------------------*
 293       EINCMN = 0.D+00
 294 *  +-------------------------------------------------------------------*
 295 *  | Energy is so small that we cannot produce cascade nucleons,
 296 *  | sample an excitation energy and return
 297       IF ( EINCMX .LE. EINCMN + 0.25D+00 * EBNDAV ) THEN
 298          ISAMPL = 50
 299 *  |
 300 *  +-------------------------------------------------------------------*
 301 *  |
 302       ELSE
 303          ISAMPL = 0
 304       END IF
 305 *  |
 306 *  +-------------------------------------------------------------------*
 307 *  +-------------------------------------------------------------------*
 308 *  |  Now make the linear sampling for Eincp, Eincn, on [0,Eicnmx]
 309 *  |  according to P(E)dE=A(Efin-E)^a dE
 310 5000  CONTINUE
 311          ISAMPL = ISAMPL + 1
 312          CALL GRNDM(RNDM,1)
 313          IF ( LLPOW ) THEN
 314             EIHELP = EINFIN * ( 1.D+00 - RNDM (1)**(1.D+00
 315      &             /(ALFA+1.D+00)) )
 316          ELSE IF ( LLLIN ) THEN
 317             EINCMX = EINCM0
 318             EIHELP = EINFIN - SQRT ( EINFIN**2 - 2.D+00 * RNDM (1)
 319      &             / AHNORM )
 320          ELSE
 321             EIHELP = EINCMX * RNDM (1)
 322          END IF
 323 *  |  *** end linear sampling ***
 324          EINCP  = EINCP0 * EIHELP / EINCT
 325          EINCN  = EINCN0 * EIHELP / EINCT
 326          AHELP  = EKMNNU (1) / ESLOPE (1)
 327          AHELP  = AHELP * FTVTH (1) * EINCP * ( 1.D+00 + AHELP *
 328      &            0.3333333333333333D+00 ) / ( 1.D+00 - EXUPNU (1) )
 329          TVGRE0 = AHELP
 330          AHELP  = EKMNNU (2) / ESLOPE (2)
 331          AHELP  = AHELP * FTVTH (2) * EINCN * ( 1.D+00 + AHELP *
 332      &            0.3333333333333333D+00 ) / ( 1.D+00 - EXUPNU (2) )
 333          TVGRE0 = TVGRE0 + AHELP
 334 *  |  +----------------------------------------------------------------*
 335 *  |  |  Energy is so small that we cannot produce cascade nucleons,
 336 *  |  |  sample an excitation energy and return
 337          IF ( ISAMPL .GE. 50 ) THEN
 338             CALL GRNDM(RNDM,3)
 339             PRNDM  = MAX ( RNDM (1), RNDM (2), RNDM (3) )
 340             TVGRE0 = MAX ( AV0WEL - PRNDM**2 * AEFRMX - EBNDAV, ZERZER )
 341             EINCP  = 0.D+00
 342             EINCN  = 0.D+00
 343             TVGREY = 0.D+00
 344             RETURN
 345          END IF
 346 *  |  |
 347 *  |  +----------------------------------------------------------------*
 348 *  | Changed:
 349          IF ( EIHELP + TVGRE0 .GE. EIUSE )  GO TO 5000
 350 *  |
 351 *  +--<--<--<--<--<  Resampling!
 352 *  |
 353 *  +-------------------------------------------------------------------*
 354       PCR    = EIHELP / EINCT
 355       FRAINC = FRAC * PCR
 356       TVGREY = 0.D+00
 357       AGREYP = EINCP / EKINAV (1)
 358       AGREYN = EINCN / EKINAV (2)
 359 * ==== Discretize the distribution !!! ==== *
 360       CALL GRNDM(RNDM,2)
 361       NGREYP = INT ( AGREYP )
 362       IF ( RNDM(1) .LT. AGREYP - NGREYP ) NGREYP = NGREYP + 1
 363       NGREYN = INT ( AGREYN )
 364       IF ( RNDM(2) .LT. AGREYN - NGREYN ) NGREYN = NGREYN + 1
 365 * ====
 366       AGREYT = AGREYP + AGREYN
 367       NGREYT = NGREYP + NGREYN
 368       PPROCS = AGREYP / AGREYT
 369       NGREYP = 0
 370       NGREYN = 0
 371       EINCP  = 0.D+00
 372       EINCN  = 0.D+00
 373       IF ( NGREYT .LE. 0 ) GO TO 9000
 374       ARDUMM = MAX ( 1.D+00, BBOLD - 2.D+00 )
 375       CALL RBKINI ( 1, .TRUE., EXDUMM, TKDUMM, TSDUMM,
 376      &               PSDUMM, ARDUMM, TRDUMM )
 377       V0SAV1 = V0WELL (1)
 378       V0SAV2 = V0WELL (2)
 379       VESAV1 = VEFFNU (1)
 380       VESAV2 = VEFFNU (2)
 381       V0WELL (1) = V0WSAV (1)
 382       V0WELL (2) = V0WSAV (2)
 383       VEFFNU (1) = VEWSAV (1)
 384       VEFFNU (2) = VEWSAV (2)
 385       IRETRY = 0
 386       DO 8000 I = 1, NGREYT
 387          CALL GRNDM(RNDM,1)
 388          RNDMPR = RNDM (1)
 389          ARDUMM = MAX ( 1.D+00, ARDUMM - 1.D+00 )
 390  7500    CONTINUE
 391          IF ( RNDMPR .LT. PPROCS ) THEN
 392             NGREYP = NGREYP + 1
 393             CALL RBKINI ( 1, .FALSE., EXDUMM, TKDUMM, TSDUMM,
 394      &                    PSDUMM, ARDUMM, TRDUMM )
 395             EINCP = EINCP + TKDUMM
 396             IF ( EINCP + EINCN .GT. EINCMX .AND. IRETRY .LT. 5 ) THEN
 397                EINCP = EINCP - TKDUMM
 398                CALL RBKMIN (1)
 399                IRETRY = IRETRY + 1
 400                GO TO 7500
 401             END IF
 402          ELSE
 403             NGREYN = NGREYN + 1
 404             CALL RBKINI ( 2, .FALSE., EXDUMM, TKDUMM, TSDUMM,
 405      &                    PSDUMM, ARDUMM, TRDUMM )
 406             EINCN = EINCN + TKDUMM
 407             IF ( EINCP + EINCN .GT. EINCMX .AND. IRETRY .LT. 5 ) THEN
 408                EINCN = EINCN - TKDUMM
 409                CALL RBKMIN (2)
 410                IRETRY = IRETRY + 1
 411                GO TO 7500
 412             END IF
 413          END IF
 414  8000 CONTINUE
 415       V0WELL (1) = V0SAV1
 416       V0WELL (2) = V0SAV2
 417       VEFFNU (1) = VESAV1
 418       VEFFNU (2) = VESAV2
 419  9000 CONTINUE
 420       EINCT  = EINCP  + EINCN
 421       IF ( EINCT + TVGRE0 .GT. EIUSE ) THEN
 422          TVGRE0 = 0.D+00
 423          EIUSE = MIN ( ONEONE, EINCMX / EINCT )
 424          EINCP = EINCP * EIUSE
 425          EINCN = EINCN * EIUSE
 426          EINCT = EINCP + EINCN
 427       END IF
 428       RETURN
 429 *=== End of subroutine corrin =========================================*
 430       END