This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / 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