GEANT321/fluka/ferhav.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 02/07/94  17.46.10  by  S.Giani
  11 *-- Author :
  12       SUBROUTINE FERHAV ( KP, EPROJ, PPROJ, TXX, TYY, TZZ )
  13
  14 #include "geant321/dblprc.inc"
  15 #include "geant321/dimpar.inc"
  16 #include "geant321/iounit.inc"
  17 #include "geant321/balanc.inc"
  18 #include "geant321/finlsp.inc"
  19 #include "geant321/hadflg.inc"
  20 #include "geant321/nucdat.inc"
  21 #include "geant321/qquark.inc"
  22 #include "geant321/part3.inc"
  23 #include "geant321/resnuc.inc"
  24       COMMON /FKABLT/ AM(110), GAA(110), TAU(110), ICH(110), IBAR(110),
  25      &                K1(110), K2(110)
  26       COMMON / FKNUCF / DELEFT, EKRECL, V0EXTR, ITTA, ITJ, LVMASS
  27       LOGICAL LVMASS, LSMPAN
  28       COMMON / FKEVNT / LNUCRI, LHADRI
  29       LOGICAL LNUCRI, LHADRI
  30       REAL FRNDM(3)
  31       SAVE ONEDUM, ZERDUM
  32       DATA ONEDUM / 1.D+00 /
  33       DATA ZERDUM / 0.D+00 /
  34       AMPROJ = AM (KP)
  35       AMTAR  = AM (ITTA)
  36       ECHCK  = EPROJ + EFRM - V0WELL (ITJ) - EKRECL - EBNDNG (ITJ)
  37       PXCHCK = PXFRM + PPROJ * TXX
  38       PYCHCK = PYFRM + PPROJ * TYY
  39       PZCHCK = PZFRM + PPROJ * TZZ
  40       UMIN2  = ( AMPROJ + AMTAR )**2
  41       P2CHCK = PXCHCK**2 + PYCHCK**2 + PZCHCK**2
  42       UMO2   = ECHCK**2  - P2CHCK
  43       IF ( ABS ( UMO2 - UMIN2 ) .GT. TWOTWO * ANGLGB * UMO2 ) THEN
  44          EPROJX = HLFHLF * ( UMO2 - AMPROJ**2 - AMTAR**2 ) / AMTAR
  45          IF ( EPROJX .LT. AMPROJ ) THEN
  46             WRITE (LUNERR,*)' Ferhav: trouble with pseudo-masses!!',
  47      &      EPROJX,AMPROJ,LVMASS
  48             EPROJX = AMPROJ
  49             PPROJX = ZERZER
  50             LRESMP = .TRUE.
  51             RETURN
  52          ELSE
  53             PPROJX = SQRT ( ( EPROJX - AMPROJ ) * ( EPROJX + AMPROJ ) )
  54          END IF
  55          ETOTX  = EPROJX + AMTAR
  56          AMTRMX = HLFHLF * ( UMO2 - AMPROJ**2 - AMTAR**2 ) / AMPROJ
  57          AMSQMX = AMTRMX**2
  58          PTOSCA = PXCHCK * TXX + PYCHCK * TYY + PZCHCK * TZZ
  59          PXTART = PXCHCK - PTOSCA * TXX
  60          PYTART = PYCHCK - PTOSCA * TYY
  61          PZTART = PZCHCK - PTOSCA * TZZ
  62          PTRASQ = PXTART**2 + PYTART**2 + PZTART**2
  63          AMTRSQ = AMTAR**2  + PTRASQ
  64          IF ( AMTRSQ .GT. AMSQMX ) THEN
  65             PPCHCK = SQRT (P2CHCK)
  66             PTOOLD = PTOSCA
  67             PTRASQ = ( AMTRMX - AMTAR ) * ( AMTRMX + AMTAR )
  68 ******            AMTRSQ = AMTRMX + PTRASQ
  69 ******            PTOSCA = SIGN (ONEONE,PTOOLD) * SQRT ( P2CHCK - PTRASQ )
  70             AMTRSQ = AMSQMX
  71             PTOSCA = SQRT ( P2CHCK - PTRASQ )
  72             ALPTUU = ( PTOSCA * SQRT ( ( PPCHCK - PTOOLD ) * ( PPCHCK
  73      &             + PTOOLD ) / ( PPCHCK - PTOSCA ) / ( PPCHCK + PTOSCA
  74      &             ) ) - PTOOLD ) / PPCHCK
  75             FNORM  = SQRT ( ONEONE + ALPTUU**2 + TWOTWO * ALPTUU
  76      &             * PTOOLD / PPCHCK )
  77             ALPTUU = ALPTUU / PPCHCK
  78             TXXX   = ( TXX + ALPTUU * PXCHCK ) / FNORM
  79             TYYY   = ( TYY + ALPTUU * PYCHCK ) / FNORM
  80             TZZZ   = ( TZZ + ALPTUU * PZCHCK ) / FNORM
  81             UMOTR2 = UMO2 + PTRASQ
  82             UMOTR  = SQRT (UMOTR2)
  83             AMTRAN = AMTRMX
  84             PPARSQ = PTOSCA**2
  85             PPARTT = PTOSCA
  86             GAMCMS = ECHCK  / UMOTR
  87             ETACMS = PPARTT / UMOTR
  88             EPRCMS = AMPROJ
  89             PPRCMS = ZERZER
  90          ELSE
  91             TXXX = TXX
  92             TYYY = TYY
  93             TZZZ = TZZ
  94             PTOOLD = PTOSCA
  95             UMOTR2 = UMO2 + PTRASQ
  96             UMOTR  = SQRT (UMOTR2)
  97             PPARTT = PTOSCA
  98             GAMCMS = ECHCK  / UMOTR
  99             ETACMS = PPARTT / UMOTR
 100             EPRCMS = HLFHLF * ( UMOTR2 + AMPROJ**2 - AMTRSQ ) / UMOTR
 101             PPRCMS = SQRT ( ( EPRCMS - AMPROJ ) * ( EPRCMS + AMPROJ ) )
 102          END IF
 103          EPRLAB = GAMCMS * EPRCMS + ETACMS * PPRCMS
 104          ETRLAB = ECHCK  - EPRLAB
 105          PPRLAB = SQRT ( ( EPRLAB - AMPROJ ) * ( EPRLAB + AMPROJ ) )
 106          PXTARG = PXCHCK - PPRLAB * TXXX
 107          PYTARG = PYCHCK - PPRLAB * TYYY
 108          PZTARG = PZCHCK - PPRLAB * TZZZ
 109          GAM    = ETRLAB / AMTAR
 110          BGX    = PXTARG / AMTAR
 111          BGY    = PYTARG / AMTAR
 112          BGZ    = PZTARG / AMTAR
 113          PPHELP = ( BGX * TXXX + BGY * TYYY + BGZ * TZZZ ) * PPRLAB
 114          ETAPCM = EPRLAB - PPHELP / ( GAM + ONEONE )
 115          PXPROJ = PPRLAB * TXXX - BGX * ETAPCM
 116          PYPROJ = PPRLAB * TYYY - BGY * ETAPCM
 117          PZPROJ = PPRLAB * TZZZ - BGZ * ETAPCM
 118          UUOLD  = PXPROJ / PPROJX
 119          VVOLD  = PYPROJ / PPROJX
 120          WWOLD  = PZPROJ / PPROJX
 121          SINT02 = UUOLD**2 + VVOLD**2
 122          IF ( SINT02 .LE. ANGLSQ ) THEN
 123             LSMPAN = .TRUE.
 124             SINTH0 = ZERZER
 125             COSPH0 = ONEONE
 126             SINPH0 = ZERZER
 127          ELSE
 128             LSMPAN = .FALSE.
 129             SINTH0 = SQRT (SINT02)
 130             COSPH0 = UUOLD / SINTH0
 131             SINPH0 = VVOLD / SINTH0
 132          END IF
 133       ELSE
 134          UMO2   = UMIN2
 135          EPROJX = AMPROJ
 136          PPROJX = ZERZER
 137          ETOTX  = AMPROJ + AMTAR
 138          LSMPAN = .FALSE.
 139          CALL POLI   ( COSTH0, SINTH0 )
 140          CALL SFECFE ( SINPH0, COSPH0 )
 141          UUOLD  = SINTH0 * COSPH0
 142          VVOLD  = SINTH0 * SINPH0
 143          WWOLD  = COSTH0
 144          AAFACT = ECHCK  + ETOTX
 145          BBFACT = PPROJX - PZCHCK
 146          DDENOM = ETOTX * AAFACT - PPROJX * BBFACT
 147          GAM = ( ECHCK * AAFACT + PPROJX * BBFACT ) / DDENOM
 148          BGZ = - BBFACT * AAFACT / DDENOM
 149          BGX = PXCHCK * ( GAM + ONEONE ) / AAFACT
 150          BGY = PYCHCK * ( GAM + ONEONE ) / AAFACT
 151       END IF
 152       PLABS  = PPROJX
 153       ELABS  = EPROJX
 154       IF ( PLABS .LT. 1.D-04 ) THEN
 155          WRITE (LUNERR,*)' Ferhav: kp,plabs,elabs,pprox,y,z,pfrmix,y,z'
 156      &   ,KP,PLABS,ELABS,PPROJ*TXX,PPROJ*TYY,PPROJ*TZZ,PXFRM,PYFRM,
 157      &    PZFRM
 158          WRITE (LUNERR,*)'   Lvmass,Am(kp),Eproj:',LVMASS,AM(KP),EPROJ
 159       ELSE IF ( PLABS .GT. 1.D+01 ) THEN
 160          WRITE (LUNERR,*)' Ferhav: kp,plabs,elabs,pprox,y,z,pfrmix,y,z'
 161      &   ,KP,PLABS,ELABS,PPROJ*TXX,PPROJ*TYY,PPROJ*TZZ,PXFRM,PYFRM,
 162      &    PZFRM
 163          WRITE (LUNERR,*)'   Lvmass,Am(kp),Eproj:',LVMASS,AM(KP),EPROJ
 164       END IF
 165       ISSU = 0
 166       DO 100 IQ = 1,3
 167          ISSU = ISSU + MQUARK (IQ,KP) / 3
 168   100 CONTINUE
 169       IF ( LVMASS ) THEN
 170          LHADRI = .TRUE.
 171          CALL HADRIN ( KP, PLABS, ELABS, ZERDUM, ZERDUM, ONEDUM, ITTA )
 172          IOLDHD = 0
 173       ELSE IF ( PLABS .GT. 7.D+00 ) THEN
 174          LHADRI = .FALSE.
 175          CALL HINHEV ( KP, PLABS, ELABS, ITTA )
 176       ELSE
 177          LHADRI = .TRUE.
 178          CALL HADRIV ( KP, PLABS, ELABS, ZERDUM, ZERDUM, ONEDUM, ITTA )
 179       END IF
 180       DO 2000 I=1,IR
 181          ECMS  = ELR (I)
 182          PCMSX = PLR (I) * CXR (I)
 183          PCMSY = PLR (I) * CYR (I)
 184          PCMSZ = PLR (I) * CZR (I)
 185          IF ( LSMPAN ) THEN
 186             PCMSX = PLR (I) * CXR (I)
 187             PCMSY = PLR (I) * CYR (I)
 188             PCMSZ = WWOLD * PLR (I) * CZR (I)
 189          ELSE
 190             PLRX = CXR (I) * COSPH0 * WWOLD - CYR (I) * SINPH0
 191      &           + CZR (I) * UUOLD
 192             PLRY = CXR (I) * SINPH0 * WWOLD + CYR (I) * COSPH0
 193      &           + CZR (I) * VVOLD
 194             PLRZ = - CXR (I) * SINTH0 + CZR (I) * WWOLD
 195             PCMSX = PLRX * PLR (I)
 196             PCMSY = PLRY * PLR (I)
 197             PCMSZ = PLRZ * PLR (I)
 198          END IF
 199          CALL ALTRA ( GAM, BGX, BGY, BGZ, PCMSX, PCMSY, PCMSZ,
 200      &                ECMS, PLR (I), PLRX, PLRY, PLRZ, ELR (I) )
 201          CXR (I) = PLRX / PLR (I)
 202          CYR (I) = PLRY / PLR (I)
 203          CZR (I) = PLRZ / PLR (I)
 204          DO 200 IQ = 1,3
 205             ISSU = ISSU - MQUARK (IQ,KPTOIP(ITR(I))) / 3
 206   200    CONTINUE
 207 2000  CONTINUE
 208       IF ( ISSU .NE. 0 ) THEN
 209          WRITE (LUNOUT,*)' *** Strangeness non conservation in Hadriv',
 210      &                     ISSU,KP,ITTA,' ***'
 211          WRITE (LUNERR,*)' *** Strangeness non conservation in Hadriv',
 212      &                     ISSU,KP,ITTA,' ***'
 213          LRESMP = .TRUE.
 214       END IF
 215       V0WELL (ITJ) = V0WELL (ITJ) - V0EXTR
 216       RETURN
 217       ENTRY FERSET
 218       FERM = PFRMMX (ITJ)
 219       CALL GRNDM(FRNDM,3)
 220       P2 = MAX ( FRNDM (1), FRNDM (2), FRNDM (3) )
 221       IF ( IBTAR .LE. 1 ) THEN
 222          FERM = ZERZER
 223       END IF
 224       P2=FERM*P2
 225       P2SQ   = P2 * P2
 226       IATEMP = IBTAR - 1
 227       ATEMP  = DBLE ( IBTAR ) - ONEONE
 228       IF ( ITJ .EQ. 1 ) THEN
 229          IZTEMP = ICHTAR - 1
 230       ELSE
 231          IZTEMP = ICHTAR
 232       END IF
 233       ZTEMP =  DBLE ( IZTEMP )
 234       DELCTR = ( DBLE (ICHTAR) - ZTEMP ) * AMELEC
 235       DELEFT = AMMTAR - AMNTAR - DELCTR
 236       AMMRES = AMUAMU * ATEMP + 1.D-03 * FKENER ( ATEMP, ZTEMP )
 237       AMNRES = AMMRES - ZTEMP * AMELEC + ELBNDE ( IZTEMP )
 238       AMTMSQ = AMMRES * AMMRES
 239       EKRECL = SQRT ( AMTMSQ + P2SQ ) - AMMRES
 240       CALL POLI ( POLC, POLS )
 241       CALL COSI ( SFE,  CFE )
 242       PXFRM = CFE * POLS * P2
 243       PYFRM = SFE * POLS * P2
 244       PZFRM = POLC * P2
 245       EFRM  = SQRT ( AMNUSQ (ITJ) + P2SQ )
 246       EKFER  = EFRM  - AMNUCL (ITJ)
 247       TVEUZ  = V0WELL (ITJ) - EFRM + EBNDNG (ITJ) + AMMTAR - AMMRES
 248      &       - DELCTR
 249       IF ( TVEUZ .LT. ZERZER ) THEN
 250          V0EXTR = - TVEUZ + TENTEN * TVEPSI
 251          TVEUZ  = TVEUZ  + V0EXTR
 252          V0WELL (ITJ) = V0WELL (ITJ) + V0EXTR
 253       ELSE
 254          V0EXTR = ZERZER
 255       END IF
 256       RETURN
 257       END