GEANT321/fluka/pmprab.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 PMPRAB.FOR
  13 *COPY PMPRAB
  14 *
  15 *=== pmprab ===========================================================*
  16 *
  17       SUBROUTINE PMPRAB ( KPROJ, EKIN, PPROJ, TXX, TYY, TZZ, WEE )
  18
  19 #include "geant321/dblprc.inc"
  20 #include "geant321/dimpar.inc"
  21 #include "geant321/iounit.inc"
  22 *
  23 *----------------------------------------------------------------------*
  24 *                                                                      *
  25 *     Created on 22 september 1991 by    Alfredo Ferrari & Paola Sala  *
  26 *                                                   Infn - Milan       *
  27 *                                                                      *
  28 *     Last change on 22-sep-91     by    Alfredo Ferrari               *
  29 *                                                                      *
  30 *                                                                      *
  31 *----------------------------------------------------------------------*
  32 *
  33 #include "geant321/balanc.inc"
  34 #include "geant321/finuc.inc"
  35 #include "geant321/nucdat.inc"
  36 #include "geant321/nucgeo.inc"
  37 #include "geant321/parevt.inc"
  38 #include "geant321/parnuc.inc"
  39 #include "geant321/part.inc"
  40 #include "geant321/resnuc.inc"
  41 *
  42       REAL RNDM(1)
  43 *
  44       IF ( KPROJ .NE. 14 .OR. EKIN .GT. 2.D+00 * GAMMIN .OR. IBTAR .NE.
  45      &     1 .OR. ICHTAR .NE. 1 ) THEN
  46          WRITE (LUNOUT,*)' **** Pmprab: kproj,ekin,ibtar,ichtar',
  47      &                     KPROJ,EKIN,IBTAR,ICHTAR
  48          WRITE (LUNERR,*)' **** Pmprab: kproj,ekin,ibtar,ichtar',
  49      &                     KPROJ,EKIN,IBTAR,ICHTAR
  50       END IF
  51       PXRES  = PXTTOT
  52       PYRES  = PYTTOT
  53       PZRES  = PZTTOT
  54       PTRES  = PTTOT
  55       CALL GRNDM(RNDM,1)
  56       RNDPAN = RNDM (1)
  57       IF ( RNDPAN .GE. 1.D+00 / PNFRAT ) THEN
  58          ERES   = EKIN + AM (KPROJ) + EKFERM + AM (1)
  59          UMO2   = ( ERES - PTRES ) * ( ERES + PTRES )
  60          UMO    = SQRT (UMO2)
  61          GAMCM = ERES  / UMO
  62          ETAX  = PXRES / UMO
  63          ETAY  = PYRES / UMO
  64          ETAZ  = PZRES / UMO
  65          ECMSNU = 0.5D+00 * ( UMO2 + AMNUSQ (2) ) / UMO
  66          PCMS   = UMO - ECMSNU
  67          CALL RACO ( PCMSX, PCMSY, PCMSZ )
  68          PCMSX = PCMS * PCMSX
  69          PCMSY = PCMS * PCMSY
  70          PCMSZ = PCMS * PCMSZ
  71          ETAPCM = ETAX * PCMSX + ETAY * PCMSY + ETAZ * PCMSZ
  72          PHELP  = PCMS + ETAPCM / ( GAMCM + 1.D+00 )
  73          ENPHOT = GAMCM * PCMS + ETAPCM
  74          PXHELP = PCMSX + ETAX * PHELP
  75          PYHELP = PCMSY + ETAY * PHELP
  76          PZHELP = PCMSZ + ETAZ * PHELP
  77          PXRES = PXRES - PXHELP
  78          PYRES = PYRES - PYHELP
  79          PZRES = PZRES - PZHELP
  80          ERES  = ERES  - ENPHOT
  81          NP = NP + 1
  82          TKI   (NP) = ENPHOT
  83          KPART (NP) = 7
  84          PLR   (NP) = ENPHOT
  85          CXR   (NP) = PXHELP / PLR (NP)
  86          CYR   (NP) = PYHELP / PLR (NP)
  87          CZR   (NP) = PZHELP / PLR (NP)
  88          WEI   (NP) = WEE
  89          IOTHER = IOTHER + 1
  90          PXNUCR = PXNUCR + PXHELP
  91          PYNUCR = PYNUCR + PYHELP
  92          PZNUCR = PZNUCR + PZHELP
  93          ENUCR  = ENUCR  + TKI (NP)
  94          IBNUCR = IBNUCR + IBAR (KPART(NP))
  95          ICNUCR = ICNUCR + ICH  (KPART(NP))
  96          ETAPCM = - ETAPCM
  97          PHELP  = ECMSNU + ETAPCM / ( GAMCM + 1.D+00 )
  98          ENNEU  = GAMCM * ECMSNU + ETAPCM
  99          PXHELP = -PCMSX + ETAX * PHELP
 100          PYHELP = -PCMSY + ETAY * PHELP
 101          PZHELP = -PCMSZ + ETAZ * PHELP
 102          NP = NP + 1
 103          TKI   (NP) = ENNEU - AM (8)
 104          KPART (NP) = 8
 105          PLR (NP) = SQRT ( PXHELP**2 + PYHELP**2 + PZHELP**2 )
 106          CXR   (NP) = PXHELP / PLR (NP)
 107          CYR   (NP) = PYHELP / PLR (NP)
 108          CZR   (NP) = PZHELP / PLR (NP)
 109          WEI   (NP) = WEE
 110          ERES  = ERES  - ENNEU
 111          PXRES = PXRES - PXHELP
 112          PYRES = PYRES - PYHELP
 113          PZRES = PZRES - PZHELP
 114          IBRES = 0
 115          ICRES = 0
 116          PTRES = 0.D+00
 117          ANOW  = 0.D+00
 118          ZNOW  = 0.D+00
 119       ELSE
 120          ERES   = EKIN + AM (KPROJ) + EKFERM + AM (1)
 121          UMO2   = ( ERES - PTRES ) * ( ERES + PTRES )
 122          UMO    = SQRT (UMO2)
 123          GAMCM = ERES  / UMO
 124          ETAX  = PXRES / UMO
 125          ETAY  = PYRES / UMO
 126          ETAZ  = PZRES / UMO
 127          ECMSNU = 0.5D+00 * ( UMO2 + AM (8)* AM (8) - AM (23) * AM (23)
 128      &          ) / UMO
 129          ECMSP0 = UMO - ECMSNU
 130          PCMS = SQRT ( ( ECMSP0 - AM (23) ) * ( ECMSP0 + AM (23) ) )
 131          CALL RACO ( PCMSX, PCMSY, PCMSZ )
 132          PCMSX = PCMS * PCMSX
 133          PCMSY = PCMS * PCMSY
 134          PCMSZ = PCMS * PCMSZ
 135          ETAPCM = ETAX * PCMSX + ETAY * PCMSY + ETAZ * PCMSZ
 136          PHELP  = ECMSP0 + ETAPCM / ( GAMCM + 1.D+00 )
 137          ENPIO0 = GAMCM * ECMSP0 + ETAPCM
 138          PXHELP = PCMSX + ETAX * PHELP
 139          PYHELP = PCMSY + ETAY * PHELP
 140          PZHELP = PCMSZ + ETAZ * PHELP
 141          PXRES = PXRES - PXHELP
 142          PYRES = PYRES - PYHELP
 143          PZRES = PZRES - PZHELP
 144          ERES  = ERES  - ENPIO0
 145          NP = NP + 1
 146          TKI   (NP) = ENPIO0 - AM (23)
 147          KPART (NP) = 23
 148          PLR (NP) = SQRT ( PXHELP**2 + PYHELP**2 + PZHELP**2 )
 149          CXR   (NP) = PXHELP / PLR (NP)
 150          CYR   (NP) = PYHELP / PLR (NP)
 151          CZR   (NP) = PZHELP / PLR (NP)
 152          WEI   (NP) = WEE
 153          IOTHER = IOTHER + 1
 154          PXNUCR = PXNUCR + PXHELP
 155          PYNUCR = PYNUCR + PYHELP
 156          PZNUCR = PZNUCR + PZHELP
 157          ENUCR  = ENUCR  + TKI (NP)
 158          IBNUCR = IBNUCR + IBAR (KPART(NP))
 159          ICNUCR = ICNUCR + ICH  (KPART(NP))
 160          ETAPCM = - ETAPCM
 161          PHELP  = ECMSNU + ETAPCM / ( GAMCM + 1.D+00 )
 162          ENNEU  = GAMCM * ECMSNU + ETAPCM
 163          PXHELP = -PCMSX + ETAX * PHELP
 164          PYHELP = -PCMSY + ETAY * PHELP
 165          PZHELP = -PCMSZ + ETAZ * PHELP
 166          NP = NP + 1
 167          TKI   (NP) = ENNEU - AM (8)
 168          KPART (NP) = 8
 169          PLR (NP) = SQRT ( PXHELP**2 + PYHELP**2 + PZHELP**2 )
 170          CXR   (NP) = PXHELP / PLR (NP)
 171          CYR   (NP) = PYHELP / PLR (NP)
 172          CZR   (NP) = PZHELP / PLR (NP)
 173          WEI   (NP) = WEE
 174          ERES  = ERES  - ENNEU
 175          PXRES = PXRES - PXHELP
 176          PYRES = PYRES - PYHELP
 177          PZRES = PZRES - PZHELP
 178          IBRES = 0
 179          ICRES = 0
 180          PTRES = 0.D+00
 181          ANOW  = 0.D+00
 182          ZNOW  = 0.D+00
 183       END IF
 184       RETURN
 185 *=== End of subroutine PMPRAB =========================================*
 186       END