GEANT321/fluka/nudisv.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:19:58  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 25/07/94  15.04.20  by  S.Giani
  11 *-- Author :
  12 *$ CREATE NUDISV.FOR
  13 *COPY NUDISV
  14 *
  15 *=== nudisv ===========================================================*
  16 *
  17       INTEGER FUNCTION NUDISV ( ANU, KB, EXPLAM, ASEASQ, APOWER, PRZERO)
  18
  19 #include "geant321/dblprc.inc"
  20 #include "geant321/dimpar.inc"
  21 #include "geant321/iounit.inc"
  22 *
  23 *----------------------------------------------------------------------*
  24 *                                                                      *
  25 *     Last change  on  05-may-92    by  Alfredo Ferrari, INFN-Milan    *
  26 *                                                                      *
  27 *     Nudist91:                                                        *
  28 *      multiplicity distribution of chains produced in hadron-nucleus  *
  29 *      collisions                                                      *
  30 *                                                                      *
  31 *      Input variables:                                                *
  32 *                       anu    = average number of collisions          *
  33 *                       kb     = baryon number of the projectile       *
  34 *                       explam = exp ( - (anu-1) )                     *
  35 *      Output variables:                                               *
  36 *                       nudisv = number of high energy collisions      *
  37 *----------------------------------------------------------------------*
  38 *
  39       PARAMETER ( INMAX  = 30 )
  40       PARAMETER ( IAMAX  = 13 )
  41       PARAMETER ( ANUMAX = 6.0D+00 )
  42       PARAMETER ( ANUMED = 0.5D+00 * ( 3.75D+00 + ANUMAX ) )
  43 *
  44 C      AB,(KB.NE.0) BARYON NUCLEUS
  45 C
  46       DIMENSION AB (INMAX,IAMAX), ANUB (IAMAX), ANUM (9),
  47      &          ANSESQ (IAMAX), ANPOWR (IAMAX), NULAST (IAMAX)
  48       DIMENSION AAB (IAMAX), DIST (0:INMAX), IEX (3,5)
  49       REAL RNDM(1)
  50       LOGICAL LFIRST
  51       SAVE LFIRST, ANUB, ANUM, AB, ANSESQ, NULAST, ANPOWR
  52       DATA ANUB/1.D0,1.1D0,1.25D0,1.45D0,1.65D0,1.74D0,2.09D0,2.66D0,
  53      &          3.0D0,3.1D0,3.75D0,ANUMED,ANUMAX/
  54       DATA ANUM/1.D0,1.45D0,1.52D0,1.77D0,2.19D0,2.42D0,2.47D0,
  55      *2.90D0,5.D0/
  56       DATA AB  /   1.D+00, 29*0.D+00, 90*0.D+00,
  57      & 0.5771D+00, 0.2565D+00, 0.1159D+00, 0.0397D+00,
  58      & 0.0095D+00, 0.0013D+00, 0.0001D+00, 23*0.D+00,
  59      & 0.5457D+00, 0.2548D+00, 0.1303D+00, 0.0515D+00,
  60      & 0.0146D+00, 0.0028D+00, 0.0003D+00, 23*0.D+00,
  61      & 0.4462D+00, 0.2445D+00, 0.1599D+00, 0.0919D+00, 0.0406D+00,
  62      & 0.0132D+00, 0.0032D+00, 0.0005D+00, 0.0001D+00, 21*0.D+00,
  63      & 0.3387D+00, 0.2105D+00, 0.1670D+00, 0.1291D+00, 0.0849D+00,
  64      & 0.0441D+00, 0.0180D+00, 0.0059D+00, 0.0014D+00, 0.0003D+00,
  65      & 0.0001D+00, 19*0.D+00,
  66      & 0.2959D+00, 0.1885D+00, 0.1673D+00, 0.1344D+00, 0.1044D+00,
  67      & 0.0653D+00, 0.0336D+00, 0.0144D+00, 0.0046D+00, 0.0012D+00,
  68      & 0.0003D+00, 19*0.D+00,
  69      & 0.2835D+00, 0.1818D+00, 0.1572D+00, 0.1357D+00, 0.1082D+00,
  70      & 0.0715D+00, 0.0380D+00, 0.0161D+00, 0.0057D+00, 0.0018D+00,
  71      & 0.0004D+00, 0.0001D+00, 18*0.D+00,
  72      & 0.2247D+00, 0.1481D+00, 0.1342D+00, 0.1328D+00, 0.1211D+00,
  73      & 0.0984D+00, 0.0700D+00, 0.0394D+00, 0.0198D+00, 0.0077D+00,
  74      & 0.0028D+00, 0.0008D+00, 0.0002D+00, 17*0.D+00,  60*0.D+00/
  75       DATA IEX / 1, 3, 6, 1, 2, 3, 2, 4, 5, 10, 12, 11, 11, 13, 12 /
  76       DATA LFIRST / .TRUE. /
  77 *
  78 *  +-------------------------------------------------------------------*
  79 *  |      First call: perform the initialization
  80       IF ( LFIRST ) THEN
  81          LFIRST = .FALSE.
  82          IPOWER = NINT (-APOWER)
  83 *  |  +----------------------------------------------------------------*
  84 *  |  |
  85          DO 100 IA = 1, IAMAX
  86             AAB (IA) = 0.D+00
  87 *  |  |  +-------------------------------------------------------------*
  88 *  |  |  |          This loop computes the normalization factor
  89             DO 80 IN = 1, INMAX
  90                AAB (IA) = AAB (IA) + AB (IN,IA)
  91    80       CONTINUE
  92 *  |  |  |
  93 *  |  |  +-------------------------------------------------------------*
  94 *  |  |  +-------------------------------------------------------------*
  95 *  |  |  |
  96             IF ( AAB (IA) .GT. 0.D+00 ) THEN
  97                ANUAVE = 0.D+00
  98 *  |  |  |  +----------------------------------------------------------*
  99 *  |  |  |  |
 100                DO 90 IN = 1, INMAX
 101                   AB (IN,IA) = AB (IN,IA) / AAB (IA)
 102                   ANUAVE = ANUAVE + IN * AB (IN,IA)
 103    90          CONTINUE
 104 *  |  |  |  |
 105 *  |  |  |  +----------------------------------------------------------*
 106                ANUB (IA) = ANUAVE
 107             END IF
 108 *  |  |  |
 109 *  |  |  +-------------------------------------------------------------*
 110   100    CONTINUE
 111 *  |  |
 112 *  |  +----------------------------------------------------------------*
 113 *  |  +----------------------------------------------------------------*
 114 *  |  |
 115          DO 200 IE = 1,5
 116             IA  = IEX (2,IE)
 117             IR1 = IEX (1,IE)
 118             IR2 = IEX (3,IE)
 119             ANUEN1 = 0.D+00
 120             ANUEN2 = 0.D+00
 121             AAB (IA) = 0.D+00
 122             AEXTR1 = 0.5D+00 * MIN ( 1, IR1 - 1 )
 123             AEXTR2 = 0.5D+00 * MIN ( 1, IR2 - 1 )
 124 *  |  |  +-------------------------------------------------------------*
 125 *  |  |  |          This loop fills the ab(i,ia) array,
 126 *  |  |  |          extrapolating from the data for <nu>s,
 127 *  |  |  |          anub (ir1,r2)
 128             DO 130 IN = 1, INMAX
 129                ANUBE1 = ANUEN1
 130                ANUBE2 = ANUEN2
 131                ANUEN1 = ( AEXTR1 + IN  ) * ANUB (IA)
 132      &                / ANUB (IR1)
 133                ANUEN2 = ( AEXTR2 + IN  ) * ANUB (IA)
 134      &                / ANUB (IR2)
 135                NABE1  = INT (ANUBE1)
 136                NABE2  = INT (ANUBE2)
 137                NAEN1  = NINT (ANUEN1)
 138                NAEN2  = NINT (ANUEN2)
 139                IF ( AB (IN,IR1) .GT. 0.D+00 ) THEN
 140                   DO 110 INN = MAX (NABE1,1), NAEN1
 141                      IF (INN .GT. 1) THEN
 142                         ANUBE0 = 0.5D+00 + (INN-1)
 143                         IF (INN .GT. INMAX ) GO TO 110
 144                      ELSE
 145                         ANUBE0 = 0.D+00
 146                      END IF
 147                      ANUEN0 = 0.5D+00 + INN
 148                      ANUBEG = MAX ( ANUBE1, ANUBE0 )
 149                      ANUEND = MIN ( ANUEN1, ANUEN0 )
 150                      WEIGHT = AB (IN,IR1) * MAX ( ZERZER, ( ANUEND
 151      &                      - ANUBEG ) / ( ANUEN1 - ANUBE1 ) )
 152                      AB (INN,IA) = AB (INN,IA) + WEIGHT
 153                      AAB (IA) = AAB (IA) + WEIGHT
 154   110             CONTINUE
 155                END IF
 156                IF ( AB (IN,IR2) .GT. 0.D+00 ) THEN
 157                   DO 120 INN = MAX (NABE2,1), NAEN2
 158                      IF (INN .GT. 1) THEN
 159                         ANUBE0 = 0.5D+00 + (INN-1)
 160                         IF (INN .GT. INMAX ) GO TO 120
 161                      ELSE
 162                         ANUBE0 = 0.D+00
 163                      END IF
 164                      ANUEN0 = 0.5D+00 + INN
 165                      ANUBEG = MAX ( ANUBE2, ANUBE0 )
 166                      ANUEND = MIN ( ANUEN2, ANUEN0 )
 167                      WEIGHT = AB (IN,IR2) * MAX ( ZERZER, ( ANUEND
 168      &                      - ANUBEG ) / ( ANUEN2 - ANUBE2 ) )
 169                      AB (INN,IA) = AB (INN,IA) + WEIGHT
 170                      AAB (IA) = AAB (IA) + WEIGHT
 171   120             CONTINUE
 172                END IF
 173   130       CONTINUE
 174 *  |  |  |
 175 *  |  |  +-------------------------------------------------------------*
 176 *  |  |  +-------------------------------------------------------------*
 177 *  |  |  |
 178             DO 140 IN =1, INMAX
 179                AB (IN,IA) = AB (IN,IA) / AAB (IA)
 180   140       CONTINUE
 181 *  |  |  |
 182 *  |  |  +-------------------------------------------------------------*
 183   200    CONTINUE
 184 *  |  |
 185 *  |  +----------------------------------------------------------------*
 186 *  |  +----------------------------------------------------------------*
 187 *  |  |             This loop simply creates a cumulative distribution
 188          DO 12 IA = 1, IAMAX
 189             AAB (IA) = 0.D+00
 190             NULAST (IA) = INMAX
 191 *  |  |  +-------------------------------------------------------------*
 192 *  |  |  |          This loop computes the normalization factor
 193             DO 18 IN = INMAX, 1, -1
 194                AAB (IA) = AAB (IA) + AB (IN,IA)
 195                IF ( AB (IN,IA) .LE. 0.D+00 ) NULAST (IA) = IN - 1
 196    18       CONTINUE
 197 *  |  |  |
 198 *  |  |  +-------------------------------------------------------------*
 199             ANUAVE = AB (1,IA)
 200             ANSESQ (IA) = 0.D+00
 201             AB   (1,IA) = AB (1,IA) / AAB (IA)
 202 *  |  |  +-------------------------------------------------------------*
 203 *  |  |  |          Create the cumulative distribution
 204             DO 13 IN = 2, INMAX
 205                AB (IN,IA) = AB (IN,IA) / AAB (IA) + AB (IN-1,IA)
 206                ANUAVE = ANUAVE + IN * ( AB (IN,IA)
 207      &                - AB (IN-1,IA) )
 208                ANSESQ (IA) = ANSESQ (IA) + (IN-1)*(IN-1)
 209      &                     * ( AB (IN,IA) - AB (IN-1,IA) )
 210    13       CONTINUE
 211 *  |  |  |
 212 *  |  |  +-------------------------------------------------------------*
 213             ANUB   (IA) = ANUAVE
 214 *  |  |  +-------------------------------------------------------------*
 215 *  |  |  |
 216             IF ( IPOWER .LT. 7 ) THEN
 217                ANPOWR (IA) = AB (1,IA)
 218 *  |  |  |
 219 *  |  |  +-------------------------------------------------------------*
 220 *  |  |  |
 221             ELSE IF ( IPOWER .LT. 11 ) THEN
 222                ANPOWR (IA) = AB (1,IA) * FPOWER ( IPOWER, 1, ANUAVE )
 223 *  |  |  |
 224 *  |  |  +-------------------------------------------------------------*
 225 *  |  |  |
 226             ELSE
 227                ANPOWR (IA) = AB (1,IA) * FPOWER ( IPOWER, 1, ANUAVE )
 228             END IF
 229 *  |  |  |
 230 *  |  |  +-------------------------------------------------------------*
 231 *  |  |  +-------------------------------------------------------------*
 232 *  |  |  |          Compute < nu**y(nu) >:
 233             DO 11 IN = 2, INMAX
 234                IF ( IPOWER .LT. 7 ) THEN
 235                   IF ( APOWER .GT. 0.D+00 ) THEN
 236                      DPOWER = APOWER
 237                   ELSE
 238                      DPOWER = FPOWER ( IPOWER, IN, ANUAVE )
 239                   END IF
 240                   ANPOWR (IA) = ANPOWR (IA) + IN**DPOWER
 241      &                        * ( AB (IN,IA) - AB (IN-1,IA) )
 242                ELSE IF ( IPOWER .LT. 11 ) THEN
 243                   ANPOWR (IA) = ANPOWR (IA) + IN
 244      &                        * FPOWER ( IPOWER, IN, ANUAVE )
 245      &                        * ( AB (IN,IA) - AB (IN-1,IA) )
 246                ELSE
 247                   ANPOWR (IA) = ANPOWR (IA) + IN
 248      &                        **FPOWER ( IPOWER, IN, ANUAVE )
 249      &                        * ( AB (IN,IA) - AB (IN-1,IA) )
 250                END IF
 251    11       CONTINUE
 252 *  |  |  |
 253 *  |  |  +-------------------------------------------------------------*
 254    12    CONTINUE
 255 *  |  |
 256 *  |  +----------------------------------------------------------------*
 257 *  |  +----------------------------------------------------------------*
 258 *  |  |
 259          DO 17 IA = 1, IAMAX
 260             ANUAC = AB (1,IA)
 261             ANUSQ = AB (1,IA)
 262             ANUTH = AB (1,IA)
 263 *  |  |  +-------------------------------------------------------------*
 264 *  |  |  |
 265             DO 15 I = 2, NULAST (IA)
 266                ANUAC = ANUAC + I     * ( AB (I,IA) - AB(I-1,IA) )
 267                ANUSQ = ANUSQ + I*I   * ( AB (I,IA) - AB(I-1,IA) )
 268                ANUTH = ANUTH + I*I*I * ( AB (I,IA) - AB(I-1,IA) )
 269    15       CONTINUE
 270 *  |  |  |
 271 *  |  |  +-------------------------------------------------------------*
 272 *  |  |  +-------------------------------------------------------------*
 273 *  |  |  |
 274             IF ( IA .GT. 1 ) THEN
 275                ANUSE2 = ( ANUSQ - 2.D+00 * ANUB (IA) + 1.D+00 ) /
 276      &                ( ANUB (IA) - 1.D+00 )**2
 277                ANUSE3 = ( ANUTH - 3.D+00 * ANUSQ + 3.D+00 * ANUB (IA)
 278      &                - 1.D+00 ) / ( ANUB (IA) - 1.D+00 )**3
 279 *  |  |  |
 280 *  |  |  +-------------------------------------------------------------*
 281 *  |  |  |
 282             ELSE
 283                ANUSE2 = 1.D+00
 284                ANUSE3 = 1.D+00
 285             END IF
 286 *  |  |  |
 287 *  |  |  +-------------------------------------------------------------*
 288    17    CONTINUE
 289 *  |  |
 290 *  |  +----------------------------------------------------------------*
 291          NUDISV = -1
 292          RETURN
 293       END IF
 294 *  |
 295 *  +-------------------------------------------------------------------*
 296 *  +-------------------------------------------------------------------*
 297 *  |                Select the proper distributions for interpolations
 298       DO 40 I = 1, IAMAX
 299          IF (ANU .LT. ANUB(I)) GO TO 41
 300    40 CONTINUE
 301 *  |
 302 *  +-------------------------------------------------------------------*
 303       I=IAMAX + 1
 304    41 CONTINUE
 305       NANU  = I - 1
 306       IF (NANU.GE.IAMAX) GO TO 50
 307       IF (NANU.LE.0) NANU = 1
 308       NANUN = NANU + 1
 309       WEIGH1 = ( ANU - ANUB (NANU)  ) / ( ANUB (NANUN) - ANUB (NANU) )
 310       WEIGH2 = ( ANUB (NANUN) - ANU ) / ( ANUB (NANUN) - ANUB (NANU) )
 311       DIST (0) = 0.D+00
 312       ASEASQ = WEIGH1 * ANSESQ (NANUN) + WEIGH2 * ANSESQ (NANU)
 313       APOWER = WEIGH1 * ANPOWR (NANUN) + WEIGH2 * ANPOWR (NANU)
 314 *  +-------------------------------------------------------------------*
 315 *  |   Create the proper cumulative distribution by interpolation
 316 *  |   ( note that since weigh1 + weigh2 = 1 it will be automatically
 317 *  |   normalized )
 318       DO 20 IN = 1, NULAST (NANUN)
 319           DIST (IN) = WEIGH1 * AB (IN,NANUN) + WEIGH2 * AB (IN,NANU)
 320           WEIGHT = (IN-1) * ( DIST (IN) - DIST (IN-1) )
 321    20 CONTINUE
 322 *  |
 323 *  +-------------------------------------------------------------------*
 324       PRZERO = DIST (1)
 325       CALL GRNDM(RNDM,1)
 326       X=RNDM(1)
 327 *  +-------------------------------------------------------------------*
 328 *  |   Compute the proper nu from the cumulative distribution
 329       DO 30 IN = 1, NULAST (NANUN)
 330          IF ( X .LE. DIST (IN) ) GO TO 31
 331    30 CONTINUE
 332 *  |
 333 *  +-------------------------------------------------------------------*
 334       IN = NULAST (NANUN)
 335    31 CONTINUE
 336       NUDISV = IN
 337       RETURN
 338 *  Come here for <nu> larger than 8
 339    50 CONTINUE
 340       APOWER = -1.D+00
 341       ASEASQ = -1.D+00
 342       NUD    = KPOIS (EXPLAM)
 343       NUDISV = NUD + 1
 344       PRZERO = EXPLAM
 345       WRITE (LUNERR,*)' *** Nudisv called with <nu> >= 8 !!',REAL(ANU),
 346      &                ' ***'
 347       RETURN
 348 *=== end of function Nudisv ===========================================*
 349       END