GEANT321/fluka/evevap.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:19:56  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 EVEVAP.FOR
  13 *COPY EVEVAP
  14 *
  15 *=== evevap ===========================================================*
  16 *
  17       SUBROUTINE EVEVAP ( WEE )
  18
  19 #include "geant321/dblprc.inc"
  20 #include "geant321/dimpar.inc"
  21 #include "geant321/iounit.inc"
  22 *
  23 *----------------------------------------------------------------------*
  24 *                                                                      *
  25 *  EVent EVAPoration: this routine is used to steer both the evapora-  *
  26 *  tion, the high energy fission, possibly a future fragmentation      *
  27 *  and the gamma deexcitation routines                                 *
  28 *                                                                      *
  29 *  Created  on  15  may  1991   by   Alfredo Ferrari & Paola Sala      *
  30 *                                             INFN - Milan             *
  31 *                                                                      *
  32 *  Last change  on 19-apr-93    By   Alfredo Ferrari, INFN - Milan     *
  33 *                                                                      *
  34 *----------------------------------------------------------------------*
  35 *
  36 #include "geant321/balanc.inc"
  37 #include "geant321/eva1.inc"
  38 #include "geant321/fheavy.inc"
  39 #include "geant321/finuc.inc"
  40 #include "geant321/hetc5.inc"
  41 #include "geant321/hetc7.inc"
  42 #include "geant321/hettp.inc"
  43 #include "geant321/higfis.inc"
  44 #include "geant321/labcos.inc"
  45 #include "geant321/nucdat.inc"
  46 #include "geant321/parevt.inc"
  47 #include "geant321/part.inc"
  48 #include "geant321/resnuc.inc"
  49 *
  50       PARAMETER ( AMUMEV = 1.D+03 * AMUAMU )
  51 *
  52       COMMON /FKEVNT/ LNUCRI, LHADRI
  53       LOGICAL LNUCRI, LHADRI
  54 *  The initial excitation energy, mass and charge of the nucleus are
  55 *  put into Ex, Apr, Zpr (common Hetc5)
  56       EX  = MAX ( 1000 * TVCMS, ANGLGB )
  57       APR = ANOW
  58       ZPR = ZNOW
  59 *  Reset the fission/fragmentation counter:
  60       NFISS = 0
  61 *  Ammres is the atomic mass of the residual nucleus
  62 *  Reset accumulators for the energy conservation check (they are only
  63 *  local)
  64       EOTEST = AMMRES + TVCMS + TVRECL
  65       ETEVAP = 0.D+00
  66 *  +-------------------------------------------------------------------*
  67 *  |  Set the variables recording the recoil direction of the residual
  68 *  |  nucleus:
  69       IF ( PTRES .GT. 0.D+00 ) THEN
  70          COSLBR (1) = PXRES / PTRES
  71          COSLBR (2) = PYRES / PTRES
  72          COSLBR (3) = PZRES / PTRES
  73 *  |
  74 *  +-------------------------------------------------------------------*
  75 *  |  It can happen for pion capture for example that ptres=0
  76 *  |  ( it is always 0 if no "direct" particle is emitted )
  77       ELSE
  78          COSLBR (1) = 0.D+00
  79          COSLBR (2) = 0.D+00
  80          COSLBR (3) = 1.D+00
  81       END IF
  82 *  |
  83 *  +-------------------------------------------------------------------*
  84 *  The call to getrig is useless, since we actually need no rotation
  85 *     CALL GETRIG ( ZERZER, ZERZER, ONEONE )
  86       EREC = 1.D+03 * TVRECL
  87       CALL FKERUP (0)
  88 *  +-------------------------------------------------------------------*
  89 *  |  Check for fission/fragmentation: if it occurred loop back on the
  90 *  |  fission fragments to possibly evaporate further particles:
  91       IF ( FISINH ) THEN
  92          LRNFSS = .TRUE.
  93          FISINH = .FALSE.
  94          JFISS  = 0
  95 *  |  +----------------------------------------------------------------*
  96 *  |  |  Update the partial counters of evaporated particles
  97          DO 40 J = 1,6
  98             NPARTF (J,JFISS) = NPART (J)
  99             HEVFIS (JFISS)   = HEVSUM
 100    40    CONTINUE
 101 *  |  |
 102 *  |  +----------------------------------------------------------------*
 103 *  |  +----------------------------------------------------------------*
 104 *  |  |  The following "do" is not structured as a do since Nfiss can
 105 *  |  |  be incremented during evaporation/fragmentation of the
 106 *  |  |  previously generated fragments
 107    50    CONTINUE
 108             JFISS  = JFISS + 1
 109             AMMRES = 1.D-03 * AMFIS (JFISS)
 110             PTRES  = 1.D-03 * PPFIS (JFISS)
 111             EREC = EKFIS (JFISS)
 112             APR  = AFIS  (JFISS)
 113             ZPR  = ZFIS  (JFISS)
 114             EX   = MAX ( UFIS (JFISS), ANGLGB )
 115             COSLBR (1) = COSLFF (1,JFISS)
 116             COSLBR (2) = COSLFF (2,JFISS)
 117             COSLBR (3) = COSLFF (3,JFISS)
 118 *  |  |  The call to getrig is useless, since we need no rotation
 119 *           CALL GETRIG ( ZERZER, ZERZER, ONEONE )
 120             CALL FKERUP (JFISS)
 121             ANOW = APR
 122             ZNOW = ZPR
 123             ICHLP = NINT (ZNOW)
 124             IBHLP = NINT (ANOW)
 125 *  |  |  +-------------------------------------------------------------*
 126 *  |  |  |  If we enter this branch the present fragment has been
 127 *  |  |  |  completely evaporated without further fragmentation and
 128 *  |  |  |  it is ready for the final gamma deexcitation and for
 129 *  |  |  |  residual nuclei scoring
 130             IF ( .NOT. FISINH .AND. IBHLP .GT. 0 ) THEN
 131                AMTFIS (JFISS) = ANOW * AMUMEV +  FKENER ( ANOW, ZNOW )
 132                UTFIS  (JFISS) = UU
 133                RECFIS (JFISS) = EREC
 134                PPTFIS (JFISS) = SQRT ( EREC * ( EREC + TWOTWO
 135      &                        * ( AMTFIS (JFISS) + UTFIS (JFISS) ) ) )
 136                ATFIS  (JFISS) = ANOW
 137                ZTFIS  (JFISS) = ZNOW
 138                COSLFF (1,JFISS) = COSLBR (1)
 139                COSLFF (2,JFISS) = COSLBR (2)
 140                COSLFF (3,JFISS) = COSLBR (3)
 141                ETEVAP = ETEVAP + 1.D-03 * ( EREC + AMTFIS (JFISS)
 142      &                + UTFIS (JFISS) )
 143 *  |  |  |
 144 *  |  |  +-------------------------------------------------------------*
 145 *  |  |  |  Fragment furtherly fragmented or completely evaporated into
 146 *  |  |  |  p,n,d,t,3-He and alphas
 147             ELSE
 148                FISINH = .FALSE.
 149                ATFIS  (JFISS) = ZERZER
 150                ZTFIS  (JFISS) = ZERZER
 151             END IF
 152 *  |  |  |
 153 *  |  |  +-------------------------------------------------------------*
 154 *  |  |  +-------------------------------------------------------------*
 155 *  |  |  |  Update the partial counters of evaporated particles
 156             DO 60 J = 1,6
 157                NPARTF (J,JFISS) = NPART (J)
 158                HEVFIS (JFISS)   = HEVSUM
 159    60       CONTINUE
 160 *  |  |  |
 161 *  |  |  +-------------------------------------------------------------*
 162          IF ( JFISS .LT. NFISS ) GO TO 50
 163 *  |  |
 164 *  |  +----------------------------------------------------------------*
 165          FISINH = .FALSE.
 166       END IF
 167 *  |
 168 *  +-------------------------------------------------------------------*
 169       IEVNEU = NPART (1)
 170       IEVPRO = NPART (2)
 171       IEVDEU = NPART (3)
 172       IEVTRI = NPART (4)
 173       IEV3HE = NPART (5)
 174       IEV4HE = NPART (6)
 175       IEVAPL = IEVNEU + IEVPRO
 176       IEVAPH = IEVDEU + IEVTRI + IEV3HE + IEV4HE
 177 *  +-------------------------------------------------------------------*
 178 *  |              Add to the secondary stack the evaporated neutrons
 179       DO 100 IP = 1, NPART (1)
 180          NP = NP + 1
 181          KPART (NP) = 8
 182          TKI   (NP) = 1.D-03 * EPART ( IP, 1 )
 183          WEI   (NP) = WEE
 184          CXR   (NP) = COSEVP ( 1, IP, 1 )
 185          CYR   (NP) = COSEVP ( 2, IP, 1 )
 186          CZR   (NP) = COSEVP ( 3, IP, 1 )
 187          PLR   (NP) = SQRT ( TKI (NP) * ( TKI (NP) + 2.D+00 * AM (8) ) )
 188          ETEVAP = ETEVAP + TKI (NP) + AMHEAV (1)
 189   100 CONTINUE
 190 *  |
 191 *  +-------------------------------------------------------------------*
 192
 193 *  +-------------------------------------------------------------------*
 194 *  |              Add to the secondary stack the evaporated protons
 195       DO 200 IP = 1, NPART (2)
 196          NP = NP + 1
 197          KPART (NP) = 1
 198          TKI   (NP) = 1.D-03 * EPART ( IP, 2 )
 199          WEI   (NP) = WEE
 200          CXR   (NP) = COSEVP ( 1, IP, 2 )
 201          CYR   (NP) = COSEVP ( 2, IP, 2 )
 202          CZR   (NP) = COSEVP ( 3, IP, 2 )
 203          PLR   (NP) = SQRT ( TKI (NP) * ( TKI (NP) + 2.D+00 * AM (1) ) )
 204          ETEVAP = ETEVAP + TKI (NP) + AMHEAV (2)
 205   200 CONTINUE
 206 *  |
 207 *  +-------------------------------------------------------------------*
 208
 209 *  +-------------------------------------------------------------------*
 210 *  |         Add to the heavy stack the other evaporated (if requested)
 211       IF ( LHEAVY ) THEN
 212          NPHEAV = 0
 213 *  |  +----------------------------------------------------------------*
 214 *  |  |  Loop over the particle types:
 215          DO 400 JP = 3, 6
 216 *  |  |  +-------------------------------------------------------------*
 217 *  |  |  |
 218             DO 300 IP = 1, NPART (JP)
 219                NPHEAV = NPHEAV + 1
 220                KHEAVY (NPHEAV) = JP
 221                TKHEAV (NPHEAV) = 1.D-03 * EPART ( IP, JP )
 222                WHEAVY (NPHEAV) = WEE
 223                CXHEAV (NPHEAV) = COSEVP ( 1, IP, JP )
 224                CYHEAV (NPHEAV) = COSEVP ( 2, IP, JP )
 225                CZHEAV (NPHEAV) = COSEVP ( 3, IP, JP )
 226                PHEAVY (NPHEAV) = SQRT ( ( TKHEAV (NPHEAV) + TWOTWO
 227      &                         * AMHEAV (JP) ) * TKHEAV (NPHEAV) )
 228                ETEVAP = ETEVAP + TKHEAV (NPHEAV) + AMHEAV (JP)
 229   300       CONTINUE
 230 *  |  |  |
 231 *  |  |  +-------------------------------------------------------------*
 232   400    CONTINUE
 233 *  |  |
 234 *  |  +----------------------------------------------------------------*
 235 *  |
 236 *  +-------------------------------------------------------------------*
 237 *  |
 238       ELSE
 239          NPHEAV = 0
 240          ETEVAP = ETEVAP + 1.D-03 * HEVSUM + IEVDEU * AMHEAV (3)
 241      &          + IEVTRI * AMHEAV (4)
 242      &          + IEV3HE * AMHEAV (5)
 243      &          + IEV4HE * AMHEAV (6)
 244       END IF
 245 *  |
 246 *  +-------------------------------------------------------------------*
 247 *  +-------------------------------------------------------------------*
 248 *  |  Fission and/or fragmentation occurred:
 249       IF ( LRNFSS ) THEN
 250          TVHEAV = 1.D-03 * HEVSUM
 251          IF ( ABS ( ETEVAP - EOTEST )/ EOTEST .GT. 1.D-07 ) THEN
 252             WRITE ( LUNOUT, * )
 253      &            ' Evevap_fis: failure in energy conservation!!',
 254      &                        ETEVAP, EOTEST
 255             WRITE ( LUNERR, * )
 256      &            ' Evevap_fis: failure in energy conservation!!',
 257      &                        ETEVAP, EOTEST
 258          END IF
 259          TVCHLP = ZERZER
 260          IDEHLP = 0
 261 *  |  +----------------------------------------------------------------*
 262 *  |  |  Loop on fission/fragmentation fragments
 263          DO 5000 JFISS = 1, NFISS
 264             ANOW  = ATFIS (JFISS)
 265             ZNOW  = ZTFIS (JFISS)
 266             IBRES = NINT ( ANOW )
 267             ICRES = NINT ( ZNOW )
 268 *  |  |  +-------------------------------------------------------------*
 269 *  |  |  |  Check the residual nucleus:
 270             IF ( IBRES .EQ. 0 ) THEN
 271                AMMRES = ZERZER
 272                TVCMS  = ZERZER
 273                TVRECL = ZERZER
 274                PTRES  = ZERZER
 275                PXRES  = ZERZER
 276                PYRES  = ZERZER
 277                PZRES  = ZERZER
 278                ERES   = ZERZER
 279 *  |  |  |
 280 *  |  |  +-------------------------------------------------------------*
 281 *  |  |  |  real fragment:
 282             ELSE
 283                AMMRES = 1.D-03 * AMTFIS (JFISS)
 284                TVCMS  = 1.D-03 * UTFIS  (JFISS)
 285                TVRECL = 1.D-03 * RECFIS (JFISS)
 286                PTRES  = 1.D-03 * PPTFIS (JFISS)
 287                PXRES  = PTRES * COSLFF (1,JFISS)
 288                PYRES  = PTRES * COSLFF (2,JFISS)
 289                PZRES  = PTRES * COSLFF (3,JFISS)
 290                ERES   = AMMRES + TVCMS + TVRECL
 291                EKRES  = TVRECL
 292             END IF
 293 *  |  |  |
 294 *  |  |  +-------------------------------------------------------------*
 295 *  |  |  +-------------------------------------------------------------*
 296 *  |  |  |  Check if the deexcitation module have to be called
 297             IF ( LDEEXG ) THEN
 298                IDEEXG = 0
 299                CALL EVDEEX ( WEE )
 300                IDEHLP = IDEHLP + IDEEXG
 301 *  |  |  |
 302 *  |  |  +-------------------------------------------------------------*
 303 *  |  |  |
 304             ELSE
 305                TVCHLP = TVCHLP + TVCMS
 306             END IF
 307 *  |  |  |
 308 *  |  |  +-------------------------------------------------------------*
 309 *  |  |  +-------------------------------------------------------------*
 310 *  |  |  |  Check if fission fragments have to be put on stack
 311             IF ( LHEAVY .AND. IBRES .GT. 0 ) THEN
 312                NPHEAV = NPHEAV + 1
 313                TKHEAV (NPHEAV) = EKRES
 314                PHEAVY (NPHEAV) = PTRES
 315                CXHEAV (NPHEAV) = PXRES / PTRES
 316                CYHEAV (NPHEAV) = PYRES / PTRES
 317                CZHEAV (NPHEAV) = PZRES / PTRES
 318                WHEAVY (NPHEAV) = WEE
 319                KHEAVY (NPHEAV) = 6 + JFISS
 320                AMHEAV (KHEAVY(NPHEAV)) = AMMRES
 321                IBHEAV (KHEAVY(NPHEAV)) = IBRES
 322                ICHEAV (KHEAVY(NPHEAV)) = ICRES
 323             END IF
 324 *  |  |  |
 325 *  |  |  +-------------------------------------------------------------*
 326             TVHEAV = TVHEAV + TVRECL
 327  5000    CONTINUE
 328 *  |  |
 329 *  |  +----------------------------------------------------------------*
 330          IDEEXG = IDEHLP
 331          TVCMS  = TVCHLP
 332          ANOW   = ZERZER
 333          ZNOW   = ZERZER
 334          IBRES  = 0
 335          ICRES  = 0
 336          AMMRES = ZERZER
 337          TVRECL = ZERZER
 338          PTRES  = ZERZER
 339          PXRES  = ZERZER
 340          PYRES  = ZERZER
 341          PZRES  = ZERZER
 342          ERES   = ZERZER
 343 *  |
 344 *  +-------------------------------------------------------------------*
 345 *  |  Normal evaporation:
 346       ELSE
 347          ANOW  = APR
 348          ZNOW  = ZPR
 349          IBRES = NINT ( ANOW )
 350          ICRES = NINT ( ZNOW )
 351 *  |  Ammres is the atomic mass of the residual nucleus
 352 *  |  +----------------------------------------------------------------*
 353 *  |  |  Check the residual nucleus:
 354          IF ( IBRES .EQ. 0 ) THEN
 355             AMMRES = ZERZER
 356             TVCMS  = ZERZER
 357             TVRECL = ZERZER
 358             PTRES  = ZERZER
 359             PXRES  = ZERZER
 360             PYRES  = ZERZER
 361             PZRES  = ZERZER
 362             ERES   = ZERZER
 363 *  |  |
 364 *  |  +----------------------------------------------------------------*
 365 *  |  |
 366          ELSE
 367             AMMRES = ANOW * AMUAMU + 1.D-03 * FKENER ( ANOW, ZNOW )
 368             TVCMS  = 1.D-03 * UU
 369             TVRECL = 1.D-03 * EREC
 370             PTRES  = SQRT ( TVRECL * ( TVRECL + 2.D+00 * ( AMMRES +
 371      &                      TVCMS ) ) )
 372             PXRES  = PTRES * COSLBR (1)
 373             PYRES  = PTRES * COSLBR (2)
 374             PZRES  = PTRES * COSLBR (3)
 375             ERES   = AMMRES + TVCMS + TVRECL
 376             EKRES  = TVRECL
 377          END IF
 378 *  |  |
 379 *  |  +----------------------------------------------------------------*
 380          TVHEAV = 1.D-03 * HEVSUM
 381          ETEVAP = ETEVAP + ERES
 382          IF ( ABS ( ETEVAP - EOTEST )/ EOTEST .GT. 1.D-07 ) THEN
 383             WRITE ( LUNOUT, * )
 384      &            ' Evevap: failure in energy conservation!!',
 385      &                        ETEVAP, EOTEST
 386             WRITE ( LUNERR, * )
 387      &            ' Evevap: failure in energy conservation!!',
 388      &                        ETEVAP, EOTEST
 389          END IF
 390 *  |   Check if the deexcitation module have to be called
 391          IF ( LDEEXG ) CALL EVDEEX ( WEE )
 392       END IF
 393 *  |
 394 *  +-------------------------------------------------------------------*
 395       RETURN
 396 *=== End of subroutine Evevap =========================================*
 397       END