GEANT321/fluka/evdeex.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:19:55  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 EVDEEX.FOR
  13 *COPY EVDEEX
  14 *
  15 *=== evdeex ===========================================================*
  16 *
  17       SUBROUTINE EVDEEX ( WEE )
  18
  19 #include "geant321/dblprc.inc"
  20 #include "geant321/dimpar.inc"
  21 #include "geant321/iounit.inc"
  22 *
  23 *----------------------------------------------------------------------*
  24 *                                                                      *
  25 *  Evdeex :  created on 5-10-1990 by A. Ferrari & P. Sala, INFN Milan  *
  26 *                                                                      *
  27 *    Last change  on  28-jan-93   by   Alfredo Ferrari, INFN-Milan     *
  28 *                                                                      *
  29 *    This routine provides a simple model for sampling nuclear deexci- *
  30 *    tation gammas following the evaporation step                      *
  31 *                                                                      *
  32 *----------------------------------------------------------------------*
  33 *
  34 #include "geant321/balanc.inc"
  35 #include "geant321/eva0.inc"
  36 #include "geant321/finuc.inc"
  37 #include "geant321/nucdat.inc"
  38 #include "geant321/parevt.inc"
  39 #include "geant321/resnuc.inc"
  40 *
  41 *----------------------------------------------------------------------*
  42 *  Entering the routine we have:                                       *
  43 *            Ammres is the atomic mass of the residual nucleus         *
  44 *            Ibres (Anow) its mass number                              *
  45 *            Icres (Znow) its atomic number                            *
  46 *            Eres         its total energy                             *
  47 *            Ptres        its momentum                                 *
  48 *            Pxres        x-component of the momentum                  *
  49 *            Pyres        y-component of the momentum                  *
  50 *            Pzres        z-component of the momentum                  *
  51 *            Tvrecl       kinetic energy                               *
  52 *            Tvcms        excitation energy                            *
  53 *----------------------------------------------------------------------*
  54 *
  55       PARAMETER ( C0M1E1 = 0.306 D+00 )
  56       PARAMETER ( C0E2E1 = 7.1   D-01 )
  57       PARAMETER ( HNDFE1 = 5.0   D-03 )
  58       PARAMETER ( HNDFM1 = 5.0   D-02 )
  59       PARAMETER ( HNDFE2 = 1.0   D+01 )
  60 *
  61       DIMENSION COSGAM (3)
  62       REAL RNDM(2)
  63 *
  64       IDEEXG = 0
  65       IF ( IBRES .LE. 0 .OR. TVCMS .LE. GAMMIN ) THEN
  66          TVCMS = 0.D+00
  67          RETURN
  68       END IF
  69       ECHCK  = ERES - AMMRES
  70       ECHCK0 = ECHCK
  71       IBHELP = IBRES / 2
  72       IF ( IBHELP * 2 .LT. IBRES ) THEN
  73          CALL EEXLVL ( IBRES, ICRES, DELTA, EEX2ND, EEXDUM )
  74          IPAR  = 1
  75          JPAR  = 1
  76       ELSE
  77          ICHELP = ICRES / 2
  78          JPAR   = 0
  79          IF ( ICHELP * 2 .LT. ICRES ) THEN
  80             CALL EEXLVL ( IBRES, ICRES, DELTA, EEX1ST, EEXDUM )
  81             IPAR  = 1
  82          ELSE
  83             CALL EEXLVL ( IBRES, ICRES, EEX1ST, DELTA, EEXDUM )
  84             IPAR  = 2
  85          END IF
  86       END IF
  87       DELPAI = MIN ( EEXDUM, DELTA )
  88       RNUCL  = R0NUCL * RMASS (IBRES)
  89       AINERM = 0.24D+00 * AMMRES * RNUCL * RNUCL
  90       ROTEN0 = 0.5D+00 * PLABRC * PLABRC / AINERM
  91       ENMIN  = MAX ( DELTA, TWOTWO * ROTEN0 )
  92       RNMASS = AMMRES + TVCMS
  93       UMO = RNMASS
  94       GAMCM  = ERES  / RNMASS
  95       ETAX   = PXRES / RNMASS
  96       ETAY   = PYRES / RNMASS
  97       ETAZ   = PZRES / RNMASS
  98       IF ( TVCMS .LE. ENMIN .OR. IBRES .LE. 4 ) GO TO 3000
  99       AHELP  = HNDFE1 * RMASS (IBRES) * RMASS (IBRES)
 100       CFM1E1 = C0M1E1 * HNDFM1 / AHELP
 101       CFE2E1 = C0E2E1 * HNDFE2 / AHELP
 102  1000 CONTINUE
 103          UMEV   = 1.D+03 * TVCMS
 104          ASMALL = GETA ( UMEV  , ICRES , IBRES-ICRES, ILVMOD, ISDUM,
 105      &                   AOGMAX, AOGMIN )
 106          TEMPSQ = 1.D-03 * ( TVCMS - DELPAI ) / ASMALL
 107          TEMPER = SQRT ( TEMPSQ )
 108          HHH    = TVCMS / TEMPER
 109          HHHSQ  = HHH * HHH
 110          IF(HHH.GT.88) THEN
 111             AHELP=0.0
 112          ELSE
 113             AHELP  = EXP ( - HHH )
 114          ENDIF
 115          BRE1M1 = 6.D+00 - AHELP * ( HHHSQ * HHH + 3.D+00 * HHHSQ
 116      &          + 6.D+00 * HHH )
 117          BRE2   = 20.D+00 * BRE1M1 - AHELP * ( HHHSQ * HHHSQ * HHH
 118      &          + 5.D+00 * HHHSQ * HHHSQ )
 119          BRE1M1 = BRE1M1 * ( 1.D+00 + CFM1E1 )
 120          BRE2   = BRE2 * TEMPSQ * CFE2E1
 121          CALL GRNDM(RNDM,1)
 122          IF ( RNDM (1) .LT. BRE1M1 / ( BRE1M1 + BRE2 ) ) THEN
 123             LMULT = 1
 124          ELSE
 125             LMULT = 2
 126          END IF
 127          LEXPN = 2 * LMULT + 1
 128          DDLEXP = LEXPN
 129          XDISMX = MIN ( DDLEXP, HHH )
 130          IF(XDISMX.GT.88) THEN
 131             DISMX=0.
 132          ELSE
 133             DISMX  = XDISMX ** LEXPN * EXP ( - XDISMX )
 134          ENDIF
 135  2000    CONTINUE
 136             CALL GRNDM(RNDM,2)
 137             XTENT = RNDM (1) * HHH
 138             IF(XTENT.GT.88) THEN
 139                FREJE = 0.
 140             ELSE
 141                FREJE = XTENT ** LEXPN * EXP ( - XTENT ) / DISMX
 142             ENDIF
 143          IF ( RNDM (2) .GE. FREJE ) GO TO 2000
 144          ENERG0 = XTENT * TEMPER
 145          TVCMS  = TVCMS - ENERG0
 146          RNMASS = AMMRES + TVCMS
 147          CALL RACO ( COSGAM (1), COSGAM (2), COSGAM (3) )
 148          ERNCM = 0.5D+00 * ( UMO * UMO + RNMASS * RNMASS ) / UMO
 149          EEGCM = UMO - ERNCM
 150          PCMS  = EEGCM
 151          PCMSX = PCMS * COSGAM (1)
 152          PCMSY = PCMS * COSGAM (2)
 153          PCMSZ = PCMS * COSGAM (3)
 154          ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
 155          ENERG = GAMCM * EEGCM + ETAPCM
 156          PHELP = ETAPCM / ( GAMCM + 1.D+00 ) + EEGCM
 157          PLBGX = PCMSX + ETAX * PHELP
 158          PLBGY = PCMSY + ETAY * PHELP
 159          PLBGZ = PCMSZ + ETAZ * PHELP
 160          ERES  = GAMCM * ERNCM - ETAPCM
 161          EKRES = ERES - RNMASS
 162          PHELP = - ETAPCM / ( GAMCM + 1.D+00 ) + ERNCM
 163          PXRES = - PCMSX + ETAX * PHELP
 164          PYRES = - PCMSY + ETAY * PHELP
 165          PZRES = - PCMSZ + ETAZ * PHELP
 166          ECHCK  = ECHCK - ENERG
 167          IDEEXG = IDEEXG + 1
 168          NP     = NP + 1
 169          WEI   (NP) = WEE
 170          KPART (NP) = 7
 171          TKI   (NP) = ENERG
 172          PLR   (NP) = ENERG
 173          CXR   (NP) = PLBGX / ENERG
 174          CYR   (NP) = PLBGY / ENERG
 175          CZR   (NP) = PLBGZ / ENERG
 176          GAMCM  = ERES  / RNMASS
 177          ETAX   = PXRES / RNMASS
 178          ETAY   = PYRES / RNMASS
 179          ETAZ   = PZRES / RNMASS
 180          UMO    = RNMASS
 181       IF ( TVCMS .GT. ENMIN ) GO TO 1000
 182       IF ( NP .GE. MXP ) THEN
 183          WRITE ( LUNOUT, * )' **** Finuc overflow in Evdeex,',
 184      &                      ' program stopped ****'
 185          WRITE ( LUNERR, * )' **** Finuc overflow in Evdeex,',
 186      &                      ' program stopped ****'
 187          STOP
 188       END IF
 189  3000 CONTINUE
 190          IF ( TVCMS .LE. ( 6 + 2 * JPAR ) * ROTEN0 ) THEN
 191             ENERG0 = TVCMS
 192             TVCMS  = 0.D+00
 193             CALL RACO ( COSGAM (1), COSGAM (2), COSGAM (3) )
 194             ERNCM = 0.5D+00 * ( UMO * UMO + AMMRES * AMMRES ) / UMO
 195             EEGCM = UMO - ERNCM
 196             PCMS  = EEGCM
 197             PCMSX = PCMS * COSGAM (1)
 198             PCMSY = PCMS * COSGAM (2)
 199             PCMSZ = PCMS * COSGAM (3)
 200             ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
 201             ENERG = GAMCM * EEGCM + ETAPCM
 202             PHELP = ETAPCM / ( GAMCM + 1.D+00 ) + EEGCM
 203             PLBGX = PCMSX + ETAX * PHELP
 204             PLBGY = PCMSY + ETAY * PHELP
 205             PLBGZ = PCMSZ + ETAZ * PHELP
 206             ERES  = GAMCM * ERNCM - ETAPCM
 207             EKRES = ERES - AMMRES
 208             PHELP = - ETAPCM / ( GAMCM + 1.D+00 ) + ERNCM
 209             PXRES = - PCMSX + ETAX * PHELP
 210             PYRES = - PCMSY + ETAY * PHELP
 211             PZRES = - PCMSZ + ETAZ * PHELP
 212             ECHCK  = ECHCK - ENERG
 213             IDEEXG = IDEEXG + 1
 214             NP     = NP + 1
 215             WEI   (NP) = WEE
 216             KPART (NP) = 7
 217             TKI   (NP) = ENERG
 218             PLR   (NP) = ENERG
 219             CXR   (NP) = PLBGX / ENERG
 220             CYR   (NP) = PLBGY / ENERG
 221             CZR   (NP) = PLBGZ / ENERG
 222          ELSE
 223             AIAMOM = SQRT ( 0.25D+00 + TVCMS / ROTEN0
 224      &             + 0.75D+00 * JPAR ) - 0.25D+00
 225             IF ( IPAR .EQ. 1 ) THEN
 226                IF ( JPAR .EQ. 1 ) THEN
 227                   JAMIN = 2 * INT ( AIAMOM ) + 1
 228                   IAMIN = JAMIN / 2
 229                   IF ( MOD ( IAMIN, 2 ) .GT. 0 ) THEN
 230                      EGROUN = 3.75D+00 * ROTEN0
 231                   ELSE
 232                      EGROUN = 0.75D+00 * ROTEN0
 233                   END IF
 234                ELSE
 235                   IAMIN = INT ( AIAMOM )
 236                   JAMIN = 2 * IAMIN
 237                   IF ( MOD ( IAMIN, 2 ) .GT. 0 ) THEN
 238                      EGROUN = 2.D+00 * ROTEN0
 239                   ELSE
 240                      EGROUN = 0.D+00
 241                   END IF
 242                END IF
 243             ELSE
 244                IAMIN  = INT ( AIAMOM ) / IPAR
 245                IAMIN  = IAMIN * IPAR
 246                JAMIN  = 2 * IAMIN
 247                EGROUN = 0.D+00
 248             END IF
 249             DELTAE = TVCMS + EGROUN - 0.25D+00 * ROTEN0 * JAMIN *
 250      &             ( JAMIN + 2 )
 251             DO 4000 JAMOM = JAMIN, 4, -4
 252                ENERG0 = ROTEN0 * ( 2 * JAMOM - 2 )
 253      &                + DELTAE
 254                DELTAE = 0.D+00
 255                TVCMS  = TVCMS - ENERG0
 256                RNMASS = AMMRES + TVCMS
 257                CALL RACO ( COSGAM (1), COSGAM (2), COSGAM (3) )
 258                ERNCM = 0.5D+00 * ( UMO * UMO + RNMASS * RNMASS ) / UMO
 259                EEGCM = UMO - ERNCM
 260                PCMS  = EEGCM
 261                PCMSX = PCMS * COSGAM (1)
 262                PCMSY = PCMS * COSGAM (2)
 263                PCMSZ = PCMS * COSGAM (3)
 264                ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
 265                ENERG = GAMCM * EEGCM + ETAPCM
 266                PHELP = ETAPCM / ( GAMCM + 1.D+00 ) + EEGCM
 267                PLBGX = PCMSX + ETAX * PHELP
 268                PLBGY = PCMSY + ETAY * PHELP
 269                PLBGZ = PCMSZ + ETAZ * PHELP
 270                ERES  = GAMCM * ERNCM - ETAPCM
 271                EKRES = ERES - RNMASS
 272                PHELP = - ETAPCM / ( GAMCM + 1.D+00 ) + ERNCM
 273                PXRES = - PCMSX + ETAX * PHELP
 274                PYRES = - PCMSY + ETAY * PHELP
 275                PZRES = - PCMSZ + ETAZ * PHELP
 276                ECHCK  = ECHCK - ENERG
 277                IDEEXG = IDEEXG + 1
 278                NP     = NP + 1
 279                WEI   (NP) = WEE
 280                KPART (NP) = 7
 281                TKI   (NP) = ENERG
 282                PLR   (NP) = ENERG
 283                CXR   (NP) = PLBGX / ENERG
 284                CYR   (NP) = PLBGY / ENERG
 285                CZR   (NP) = PLBGZ / ENERG
 286                GAMCM  = ERES / RNMASS
 287                ETAX   = PXRES / RNMASS
 288                ETAY   = PYRES / RNMASS
 289                ETAZ   = PZRES / RNMASS
 290                UMO    = RNMASS
 291  4000       CONTINUE
 292          END IF
 293       ECHCK  = ECHCK - EKRES
 294       IF ( ABS ( ECHCK ) .GT. 1.D-7 * ECHCK0 ) THEN
 295          WRITE ( LUNERR, * )' **** No energy conservation in Evdeex',
 296      &                      ' ****', ECHCK, IDEEXG
 297       END IF
 298       TVCMS = 0.D+00
 299       IF ( NP .GT. MXP ) THEN
 300          WRITE ( LUNOUT, * )' **** Finuc overflow in Evdeex,',
 301      &                      ' program stopped ****'
 302          WRITE ( LUNERR, * )' **** Finuc overflow in Evdeex,',
 303      &                      ' program stopped ****'
 304          STOP
 305       END IF
 306       TVRECL = EKRES
 307       P2RES = PXRES * PXRES + PYRES * PYRES + PZRES * PZRES
 308       PTRES = SQRT (P2RES)
 309 *=== End of subroutine Evdeex =========================================*
 310       RETURN
 311       END