GEANT321/fluka/rbkekv.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 29/03/94  15.41.44  by  S.Giani
  11 *-- Author :
  12 *$ CREATE RBKEKV.FOR
  13 *COPY RBKEKV
  14 *                                                                      *
  15 *=== rbkekv ===========================================================*
  16 *                                                                      *
  17       SUBROUTINE RBKEKV ( IT, EXSOP, TO, AMSS, TKIN, TSTRCK,
  18      &                    PSTRCK, ARECL, TKRECL, COD, SID )
  19
  20 #include "geant321/dblprc.inc"
  21 #include "geant321/dimpar.inc"
  22 #include "geant321/iounit.inc"
  23 *
  24 *----------------------------------------------------------------------*
  25 *     version by                     Alfredo Ferrari
  26 *                                    INFN - Milan
  27 *     last change 19 april 93 by     Alfredo Ferrari
  28 *                                    INFN - Milan
  29 *
  30 *     To be called from the high energy production
  31 *
  32 *     this is a subroutine of fluka to sample intranuclear cascade
  33 *     particles: it is based on the old Rakeka from J. Ranft
  34 *
  35 *     input variables:
  36 *        it    = type of the secondary requested; 1=proton, 2=neutron
  37 *        bbtar = atomic weight of the medium
  38 *
  39 *     output variables:
  40 *        tkin  = kinetic energy of the secondary in GeV before applying
  41 *                the nuclear well (and eventually the Coulomb barreer)
  42 *        tstrck= kinetic energy of the secondary in GeV
  43 *        pstrck= momentum of the secondary in GeV/c
  44 *        sid,cod = sine and cosine of the angle between
  45 *                  the directions of the primary
  46 *                  and secondary particles
  47 *
  48 *********************************************************************
  49 *
  50       PARAMETER ( PI   = PIPIPI )
  51       PARAMETER ( PIO2 = PI / 2.D+00 )
  52 #include "geant321/nucdat.inc"
  53       LOGICAL LINIT, LZEROI
  54       DIMENSION EXSOP (2), AMMOLD (2)
  55       DIMENSION PFINIT (50,2), TKINIT (50,2), TSINIT (50,2),
  56      &          TRINIT (50,2), NINI (2)
  57       REAL RNDM(4)
  58       SAVE AMMOLD, PFINIT, TKINIT, TSINIT, TRINIT, NINI
  59       DATA AMMOLD / 0.9382796D+00, 0.9395731D+00 /
  60 *
  61       LINIT = .FALSE.
  62       IF ( NINI (IT) .GT. 0 ) THEN
  63          PKFRMI = PFINIT (NINI(IT),IT)
  64          TKIN   = TKINIT (NINI(IT),IT)
  65          TSTRCK = TSINIT (NINI(IT),IT)
  66          TKRECL = TRINIT (NINI(IT),IT)
  67          NINI  (IT) = NINI (IT) - 1
  68          EXSOP (IT) = 0.D+00
  69          GO TO 3500
  70       ELSE
  71          GO TO 200
  72       END IF
  73       ENTRY RBKINI ( IT, LZEROI, EXSOP, TKIN, TSTRCK,
  74      &               PSTRCK, ARECL, TKRECL )
  75       LINIT = .TRUE.
  76       IF ( LZEROI ) THEN
  77          NINI (1) = 0
  78          NINI (2) = 0
  79          RETURN
  80       END IF
  81   200 CONTINUE
  82 *  In this version the low energy component has been suppressed, since
  83 *  they are now produced by the evaporation model, A. Ferrari.
  84 *  The parameters needed for the sampling have been already set in
  85 *  Corrin:
  86       EXSOP (IT) = 0.D+00
  87 *  Sample the Fermi momentum
  88       ESLPFF = ESLOPE (IT)
  89       EXUPFF = EXUPNU (IT)
  90       EKUPFF = EKUPNU (IT)
  91       EXDWFF = EXMNNU (IT)
  92       EKDWFF = EKMNNU (IT)
  93       ERCLFF = ERCLAV (IT)
  94       IRECNT = 0
  95 *  +-------------------------------------------------------------------*
  96 *  |
  97  100  CONTINUE
  98 *  |  Sample the Fermi momentum
  99          CALL GRNDM(RNDM,4)
 100          PKFRMI = PFRMMX (IT) * MAX ( RNDM (1), RNDM (2),
 101      &            RNDM (3) )
 102          PKFRSQ = PKFRMI * PKFRMI
 103          TKIN   = - ESLPFF * LOG ( EXDWFF - RNDM (4) * (
 104      &              EXDWFF - EXUPFF ) )
 105          TKRECL = 0.5D+00 * PKFRSQ / ( AMUC12 * ARECL ) * ( 1.D+00 -
 106      &            0.25D+00 * PKFRSQ / ( ARECL**2 * AMUCSQ ) )
 107          TSTRCK = SQRT ( AMNUSQ (IT) + PKFRSQ ) + TKIN - AMNUCL (IT)
 108      &          - V0WELL (IT) - TKRECL - EBNDNG (IT)
 109 *  |  +----------------------------------------------------------------*
 110 *  |  |
 111          IF ( TSTRCK .LE. VEFFNU (IT) - V0WELL (IT) ) THEN
 112             EXSOP (IT) = EXSOP (IT) + TKIN
 113             IRECNT = IRECNT + 1
 114 *  |  |  +-------------------------------------------------------------*
 115 *  |  |  |
 116             IF ( IRECNT .GT. 10 ) THEN
 117 *  |  |  |  +----------------------------------------------------------*
 118 *  |  |  |  |
 119                IF ( TKRECL - ERCLFF .GT. EKUPFF - EKDWFF .AND.
 120      &              IRECNT .LT. 20 ) THEN
 121                   ERCLFF = ERCLFF + TKRECL - ERCLFF
 122                   EKUPFF = EKUPFF + TKRECL - ERCLFF
 123                   EKDWFF = EKDWFF + TKRECL - ERCLFF
 124                   AHELP  = EXP ( - ( TKRECL - ERCLFF ) / ESLPFF )
 125                   EXUPFF = EXUPFF * AHELP
 126                   EXDWFF = EXDWFF * AHELP
 127 *  |  |  |  |
 128 *  |  |  |  +----------------------------------------------------------*
 129 *  |  |  |  |
 130                ELSE IF ( IRECNT .GT. 15 ) THEN
 131                   TKRECL = MAX ( TKRECL, ERCLFF )
 132                   ERCLFF = ERCLFF + TKRECL
 133                   EKUPFF = EKUPFF + TKRECL - ERCLFF
 134                   EKDWFF = EKDWFF + TKRECL - ERCLFF
 135                   AHELP  = EXP ( - TKRECL / ESLPFF )
 136                   EXUPFF = EXUPFF * AHELP
 137                   EXDWFF = EXDWFF * AHELP
 138                END IF
 139 *  |  |  |  |
 140 *  |  |  |  +----------------------------------------------------------*
 141             END IF
 142 *  |  |  |
 143 *  |  |  +-------------------------------------------------------------*
 144             GO TO 100
 145 *  +-<|--<--<--<--<--< go to resampling
 146          END IF
 147 *  |  |
 148 *  |  +----------------------------------------------------------------*
 149 *  +-------------------------------------------------------------------*
 150 *  |
 151       IF ( LINIT ) THEN
 152          NINI (IT) = NINI (IT) + 1
 153          PFINIT (NINI(IT),IT) = PKFRMI
 154          TKINIT (NINI(IT),IT) = TKIN
 155          TSINIT (NINI(IT),IT) = TSTRCK
 156          TRINIT (NINI(IT),IT) = TKRECL
 157          RETURN
 158       END IF
 159 *  |
 160 *  +-------------------------------------------------------------------*
 161  3500 CONTINUE
 162 *  | Masses have been updated
 163       PSTRCK = SQRT ( TSTRCK * ( TSTRCK + 2.D+00 * AMNUCL (IT) ) )
 164 *
 165 ********************* Sample the angle ********************************
 166 * Polar angle selection
 167       ADE=0.090D0*(1.D0+0.081D0*ATO1O3)/TKIN
 168       DEX=EXP(- PIO2 * PIO2 / ADE)
 169       AN1=(1.D0-DEX)*ADE/2.D0
 170       AN2=DEX*PIO2
 171       AN=AN1+AN2
 172       AN1=AN1/AN
 173       CALL GRNDM(RNDM,1)
 174       IF(RNDM(1).GT.AN1)  GO TO 3
 175     2 CONTINUE
 176       CALL GRNDM(RNDM,1)
 177       DE=SQRT(-ADE*LOG(1.D0-RNDM(1)*(1.D0-DEX)))
 178       IF(DE.GT.PIO2) GO TO 2
 179 C     WRITE(LUNOUT,*)IT,TO,AMSS,SQAMSS,T,P,DE
 180       COD = COS (DE)
 181       SID = SIN (DE)
 182       RETURN
 183     3 CONTINUE
 184       CALL GRNDM(RNDM,1)
 185       COD = - RNDM(1)
 186       SID = SQRT ( (1.D0 + COD ) * ( 1.D0 - COD ) )
 187 C     WRITE(LUNOUT,*)IT,TO,AMSS,SQAMSS,T,P,DE
 188       RETURN
 189       ENTRY RBKMIN (IT)
 190       NINI(IT) = NINI(IT)-1
 191       RETURN
 192       END