This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / fluka / ferevv.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.22.04  by  S.Giani
11 *-- Author :
12       SUBROUTINE FEREVV ( NHAD, KP, KT, PM, EKM, TXX, TYY, TZZ, ATEMP,
13      &                    ZTEMP, IVVFLG )
14 #include "geant321/dblprc.inc"
15 #include "geant321/dimpar.inc"
16 #include "geant321/iounit.inc"
17 *
18 *----------------------------------------------------------------------*
19 *                                                                      *
20 *  Ferevt90 by A. Ferrari                                              *
21 *                                                                      *
22 *      Last change  on  15-aug-93   by   Alfredo Ferrari, INFN - Milan *
23 *                                                                      *
24 *      Ferevt calculates hadron nucleon collisions                     *
25 *      including the Fermi momentum of the target:                     *
26 *      now there are two entries, one for valence collisions (or for   *
27 *      sea collisions with one of the last two nucleons) and one       *
28 *      (fersea) for sea collisions.                                    *
29 *                                                                      *
30 *      Nhad = number of final hadrons                                  *
31 *      Kp and Kt are indices of the projectile and target nucleons     *
32 *      Pm   = momentum of the projectile (for ferevv entry)            *
33 *      Ekm  = kinetic energy of the projectile (for both entries)      *
34 *      Txx, Tyy, Tzz = direction cosines of the incident projectile,   *
35 *                      THEY ARE CHANGED in the routine for sea intera- *
36 *                      ctions                                          *
37 *      Atemp, Ztemp  = mass and atomic number of the residula nucleus  *
38 *                      after the "use" of the Kt nucleon               *
39 *      Kprim = index of the real projectile (only for Fersea)          *
40 *      Eprim = energy of the real projectile after the emission of     *
41 *              the virtual meson (only for Fersea)                     *
42 *      Pprim = momentum of the real projectile after the emission of   *
43 *              the virtual meson (only for Fersea)                     *
44 *      Eprold= energy of the real projectile before the emission of    *
45 *              the virtual meson (only for Fersea)                     *
46 *      Pprold= momentum of the real projectile before the emission of  *
47 *              the virtual meson (only for Fersea)                     *
48 *                                                                      *
49 *----------------------------------------------------------------------*
50 *
51 #include "geant321/balanc.inc"
52 #include "geant321/depnuc.inc"
53 #include "geant321/hadpar.inc"
54 #include "geant321/inpdat2.inc"
55 #include "geant321/nucdat.inc"
56 #include "geant321/parevt.inc"
57 #include "geant321/part2.inc"
58 #include "geant321/resnuc.inc"
59       COMMON /FKPRIN/ IPRI, INIT
60       LOGICAL LSEACL, LSMPAN
61       DIMENSION PTHRSH (NALLWP), IJNUCR (NALLWP)
62       COMMON /FKEVNT/ LNUCRI, LHADRI
63       LOGICAL LNUCRI, LHADRI
64       REAL FRNDM(3)
65       REAL FRUUN(1)
66       SAVE PTHRSH, IJNUCR
67       DATA PTHRSH / 16*5.D+00,2*2.5D+00,5.D+00,3*2.5D+00,8*5.D+00,
68      &              9*2.5D+00 /
69       DATA IJNUCR / 16*1,2*0,1,3*0,8*1,9*0 /
70       IVFLAG = IVVFLG
71       AMPROJ = AAM (KP)
72       AMTAR  = AAM (KT)
73       PSPENT = PM
74       ESPENT = EKM + AMPROJ
75       PSPTOT = PSPENT
76       ESPTOT = ESPENT
77       LSEACL = .FALSE.
78       IF ( AMPROJ .LT. 1.D+00 ) THEN
79          AFACT = 2.8D+00
80       ELSE
81          AFACT = 1.8D+00
82       END IF
83       URMIN2 = ( AMPROJ + AMTAR )**2 + ( 1.2D+00 + DBLE (IVFLAG)
84      &       * AFACT ) * ETHSEA * AMTAR
85       GO TO 100
86          ENTRY FERSEA ( NHAD, KP, KT, EKM, TXX, TYY, TZZ, ATEMP, ZTEMP,
87      &                  KPRIM, EPRIM, PPRIM, EPROLD, PPROLD )
88    50    CONTINUE
89          LSEACL = .TRUE.
90          IVFLAG = 0
91          AMPROJ = AAM (KP)
92          AMTAR  = AAM (KT)
93          PSPTOT = PPROLD
94          ESPTOT = EPROLD
95          PSPENT = PPROLD - PPRIM
96          ESPENT = EKM + AMPROJ
97          IF ( KP .EQ. 26 .AND. KT .EQ. 8 ) THEN
98             URMIN2 = ( AAM (13) + AAM (56) )**2 + 0.2D+00 * ETHSEA
99      &             * AAM (56)
100          ELSE IF ( KP .EQ. 23 .AND. KT .EQ. 1 ) THEN
101             URMIN2 = ( AAM (14) + AAM (53) )**2 + 0.2D+00 * ETHSEA
102      &             * AAM (53)
103          ELSE
104             URMIN2 = ( AMPROJ + AMTAR )**2 + 1.2D+00 * ETHSEA * AMTAR
105          END IF
106  100  CONTINUE
107       ITJ = MIN ( 2, KT )
108       B1SAVE = B1BAMJ
109       B2SAVE = B2BAMJ
110       IF ( LLASTN ) THEN
111          IF ( LLAST1 ) THEN
112             PXNUC = PXXINC
113             PYNUC = PYYINC
114             PZNUC = PZZINC
115             EKNUC = EKINC
116          ELSE
117             LLAST1 = .TRUE.
118             EKNUC  = EKLAST
119             PXNUC  = PXLAST
120             PYNUC  = PYLAST
121             PZNUC  = PZLAST
122          END IF
123          EM    = EKM   + AMPROJ
124          EFER  = EKNUC + AMTAR
125          EKFER = EFER  - AMNUCL (ITJ) + EBNDNG (ITJ)
126          EFRM  = EFRM  + EKFER
127          PXFRM = PXFRM + PXNUC
128          PYFRM = PYFRM + PYNUC
129          PZFRM = PZFRM + PZNUC
130          ECHCK  = EM + EFER
131          PXCHCK = PXNUC + PM * TXX
132          PYCHCK = PYNUC + PM * TYY
133          PZCHCK = PZNUC + PM * TZZ
134          UMO  = SQRT ( ECHCK**2 - PXCHCK**2 - PYCHCK**2 - PZCHCK**2 )
135       ELSE
136          CALL GRNDM(FRNDM,3)
137          P2  = MAX ( FRNDM (1), FRNDM (2), FRNDM (3) )
138          P2  = PFRMMX (ITJ) * P2
139          P2SQ   = P2 * P2
140          AMTEMP = ATEMP * AMUC12
141          AMTMSQ = AMTEMP * AMTEMP
142          EKRECL = 0.5D+00 * P2SQ / AMTEMP * ( 1.D+00 - 0.25D+00 * P2SQ
143      &          / AMTMSQ )
144          EKREC0 = EKRECL
145          CALL POLI   (POLC, POLS)
146          CALL SFECFE (SFE , CFE )
147          CXTA = POLS * CFE
148          CYTA = POLS * SFE
149          CZTA = POLC
150          PXNUC = P2 * CXTA
151          PYNUC = P2 * CYTA
152          PZNUC = P2 * CZTA
153          PXFRM = PXFRM + PXNUC
154          PYFRM = PYFRM + PYNUC
155          PZFRM = PZFRM + PZNUC
156          EFER  = SQRT ( AMNUSQ (ITJ) + PXNUC**2 + PYNUC**2 + PZNUC**2 )
157          FRSURV = ANOW /  DBLE ( IBTAR )
158          IF ( .NOT. LSEACL .OR. FRSURV .LT. 0.66D+00 .OR. ATEMP .LT.
159      &        40.D+00 ) THEN
160             IATEMP = NINT (ATEMP)
161             IZTEMP = NINT (ZTEMP)
162             AMMRES = AMUAMU * ATEMP + 1.D-03 * FKENER ( ATEMP, ZTEMP )
163             AMMRE2 = AMMRES * AMMRES
164             EKRECL = SQRT ( P2SQ + AMMRE2 ) - AMMRES
165             EKREC0 = EKRECL
166             ELEFT  = ETTOT  - EINTR  - EUZ - ESPTOT - EFER + EKRECL
167      &             + V0WELL (ITJ) + EBNDNG (ITJ)
168             PXLAST = PXTTOT - PXINTR - PUX - PSPTOT * TXX - PXNUC
169             PYLAST = PYTTOT - PYINTR - PUY - PSPTOT * TYY - PYNUC
170             PZLAST = PZTTOT - PZINTR - PUZ - PSPTOT * TZZ - PZNUC
171             PPLAS2 = PXLAST**2 + PYLAST**2 + PZLAST**2
172             IUMO = 0
173             DELEFT = AMMTAR - AMNTAR - ( DBLE (ICHTAR) - ZTEMP )
174      &             * AMELEC
175  200        CONTINUE
176                EELEFT = ELEFT + DELEFT
177                UMO2 = EELEFT**2 - PPLAS2
178                DELTU2 = AMMRE2 - UMO2
179                IF ( DELTU2 .GT. 0.D+00 ) THEN
180                   IUMO = IUMO + 1
181                   IF ( LSEACL ) THEN
182                      DELTAE = 0.5D+00 * DELTU2 / EELEFT
183                      IF ( IUMO .GT. 2 ) THEN
184                      ELSE IF ( DELTAE .LT. 2.D+00 * EKREC0 ) THEN
185                         EKRECL = EKRECL + DELTAE
186                         ELEFT  = ELEFT  + DELTAE
187                         GO TO 200
188                      ELSE
189                         EKRECL = EKRECL + EKREC0
190                         ELEFT  = ELEFT  + EKREC0
191                         EELEFT = ELEFT + DELEFT
192                         UMO2 = EELEFT**2 - PPLAS2
193                         DELTU2 = AMMRE2 - UMO2
194                         PFDTPL = PXNUC * PXLAST + PYNUC * PYLAST
195      &                         + PZNUC * PZLAST
196                         PPLAST = SQRT ( PPLAS2 )
197                         DELTPR = 0.51D+00 * DELTU2 / ( PPLAST - EELEFT
198      &                         * PFDTPL / ( EFER * PPLAST ) )
199                         DELTPR = SIGN ( MIN ( ABS ( DELTPR ), HLFHLF
200      &                         * P2 ), DELTPR ) / PPLAST
201                         PXXINC = PXLAST * DELTPR
202                         PYYINC = PYLAST * DELTPR
203                         PZZINC = PZLAST * DELTPR
204                         PXLAST = PXLAST - PXXINC
205                         PYLAST = PYLAST - PYYINC
206                         PZLAST = PZLAST - PZZINC
207                         PXFRM = PXFRM + PXXINC
208                         PYFRM = PYFRM + PYYINC
209                         PZFRM = PZFRM + PZZINC
210                         PXNUC = PXNUC + PXXINC
211                         PYNUC = PYNUC + PYYINC
212                         PZNUC = PZNUC + PZZINC
213                         EFER0 = EFER
214                         EFER  = SQRT ( PXNUC**2 + PYNUC**2 +
215      &                                 PZNUC**2 + AMNUSQ  (ITJ) )
216                         DELTAE = EFER0 - EFER
217                         IF ( DELTAE .GT. 0.D+00 ) THEN
218                            ELEFT = ELEFT + DELTAE
219                         ELSE
220                            EFER = EFER0
221                         END IF
222                         PPLAS2 = PXLAST**2 + PYLAST**2 + PZLAST**2
223                         DELTAE = EKREC0
224                         GO TO 200
225                      END IF
226                   ELSE
227                      DELTAE = 0.5D+00 * DELTU2 / EELEFT
228                      IF ( IUMO .GT. 3 ) THEN
229                         WRITE ( LUNOUT, * )' Ferevv: valence call,',
230      &                  ' impossible to get',
231      &                  ' enough invariant mass for the residual',
232      &                  ' nucleus!!',IATEMP,IZTEMP,AMMRE2,UMO2,UMO,
233      &                    DELTAE
234                         WRITE ( LUNERR, * )' Ferevv: valence call,',
235      &                  ' impossible to get',
236      &                  ' enough invariant mass for the residual',
237      &                  ' nucleus!!',IATEMP,IZTEMP,AMMRE2,UMO2,UMO,
238      &                    DELTAE
239                      ELSE IF ( DELTAE .LT. 3.D+00 * EKREC0 ) THEN
240                         EKRECL = EKRECL + DELTAE
241                         ELEFT  = ELEFT  + DELTAE
242                         GO TO 200
243                      ELSE
244                         EKRECL = EKRECL + 2.D+00 * EKREC0
245                         ELEFT  = ELEFT  + 2.D+00 * EKREC0
246                         EELEFT = ELEFT + DELEFT
247                         UMO2 = EELEFT**2 - PPLAS2
248                         DELTU2 = AMMRE2 - UMO2
249                         PPDTPL = PM * ( PXLAST * TXX + PYLAST * TYY
250      &                         + PZLAST * TZZ )
251                         PPLAST = SQRT ( PPLAS2 )
252                         DELTPR = 0.51D+00 * DELTU2 / ( PPLAST - EELEFT
253      &                         * PPDTPL / ( ( EKM + AMPROJ ) * PPLAST
254      &                           ) )
255                         TMPPM  = 0.3D+00 * PM
256                         TMPPPL = 0.8D+00 * PPLAST
257                         DELTPR = SIGN ( MIN ( ABS ( DELTPR ), TMPPM
258      &                          , TMPPPL ), DELTPR )
259      &                         / PPLAST
260                         PXXINC = PXLAST * DELTPR
261                         PYYINC = PYLAST * DELTPR
262                         PZZINC = PZLAST * DELTPR
263                         PXLAST = PXLAST - PXXINC
264                         PYLAST = PYLAST - PYYINC
265                         PZLAST = PZLAST - PZZINC
266                         PXFRM = PXFRM + PXXINC
267                         PYFRM = PYFRM + PYYINC
268                         PZFRM = PZFRM + PZZINC
269                         PXXINC = PM * TXX + PXXINC
270                         PYYINC = PM * TYY + PYYINC
271                         PZZINC = PM * TZZ + PZZINC
272                         PM  = SQRT ( PXXINC**2 + PYYINC**2
273      &                             + PZZINC**2 )
274                         TXX = PXXINC / PM
275                         TYY = PYYINC / PM
276                         TZZ = PZZINC / PM
277                         DELTAE = EKM
278                         EKM = SQRT ( PM * PM + AMPROJ * AMPROJ )
279      &                      - AMPROJ
280                         DELTAE = DELTAE - EKM
281                         ELEFT  = ELEFT  + DELTAE
282                         EFRM   = EFRM   - DELTAE
283                         PPLAS2 = PXLAST**2 + PYLAST**2 + PZLAST**2
284                         PSPENT = PM
285                         ESPENT = EKM + AMPROJ
286                         PSPTOT = PSPENT
287                         ESPTOT = ESPENT
288                         GO TO 200
289                      END IF
290                   END IF
291                   UMO = ( ESPENT + EFER - V0WELL (ITJ) - EKRECL
292      &                  - EBNDNG (ITJ) )**2
293      &                - ( TXX * PSPENT + PXNUC )**2
294      &                - ( TYY * PSPENT + PYNUC )**2
295      &                - ( TZZ * PSPENT + PZNUC )**2
296                   IF ( UMO .LT. URMIN2 ) THEN
297                      WRITE ( LUNOUT, * )' Ferevv: impossible to get',
298      &                      ' enough invariant mass for interaction',
299      &                        IATEMP,IZTEMP,AMMRE2,UMO2,UMO
300                      WRITE ( LUNERR, * )' Ferevv: impossible to get',
301      &                      ' enough invariant mass for interaction',
302      &                        IATEMP,IZTEMP,AMMRE2,UMO2,UMO
303                   END IF
304                END IF
305             IF ( IUMO .GT. 0 .AND. .NOT. LSEACL .AND. DELTU2 .LT.
306      &           0.D+00 ) THEN
307                DELTAE = SQRT ( DELTU2 + EELEFT**2 )
308      &                - EELEFT
309                DELTAE = DELTAE + 10.D+00 * ANGLGB * AMMRES
310                EKRECL = EKRECL + DELTAE
311                ELEFT  = ELEFT  + DELTAE
312                EELEFT = ELEFT + DELEFT
313                UMO2 = EELEFT**2 - PPLAS2
314                DELTU2 = AMMRE2 - UMO2
315             END IF
316          ELSE
317             AMTEMP = ATEMP * AMUC12
318             AMTMSQ = AMTEMP * AMTEMP
319             EKRECL = 0.5D+00 * P2SQ / AMTEMP * ( 1.D+00 - 0.25D+00
320      &             * P2SQ / AMTMSQ )
321             EKREC0 = EKRECL
322          END IF
323          EKFER = EFER - AMNUCL (ITJ)
324          IF ( LSEACL ) THEN
325             AMPRI2 = EPRIM * EPRIM - PPRIM * PPRIM
326             PXCHCK = PXNUC + PPROLD * TXX
327             PYCHCK = PYNUC + PPROLD * TYY
328             PZCHCK = PZNUC + PPROLD * TZZ
329             P2CHCK = PXCHCK**2 + PYCHCK**2 + PZCHCK**2
330             ECHCK  = EPROLD + EFER
331             UMO2   = ECHCK**2 - P2CHCK
332             IF ( UMO2 .LT. 0.D+00 ) THEN
333                LRESMP = .TRUE.
334                WRITE (LUNERR,*)' FEREVV: SEA INT. UMO2 < 0 ',UMO2
335                RETURN
336             END IF
337             UMO    = SQRT ( UMO2 )
338             RMIN2  = URMIN2 / UMO2
339             GAMCM = ECHCK  / UMO
340             ETAX  = PXCHCK / UMO
341             ETAY  = PYCHCK / UMO
342             ETAZ  = PZCHCK / UMO
343             ETACM = SQRT ( ETAX**2 + ETAY**2 + ETAZ**2 )
344             CXCMS = ETAX / ETACM
345             CYCMS = ETAY / ETACM
346             CZCMS = ETAZ / ETACM
347             CALL SFECFE ( SFE, CFE )
348             IF ( ABS (CZCMS) .GT. 1.D-04 ) THEN
349                CXXTR  = - SFE * CZCMS
350                CYYTR  = CFE * CZCMS
351                CZZTR  = CXCMS * SFE - CYCMS * CFE
352             ELSE IF ( ABS (CYCMS) .GT. 1.D-04 ) THEN
353                CXXTR  = CYCMS * CFE
354                CYYTR  = CZCMS * SFE - CXCMS * CFE
355                CZZTR  = - SFE * CYCMS
356             ELSE
357                CXXTR  = CYCMS * SFE - CZCMS * CFE
358                CYYTR  = - SFE * CXCMS
359                CZZTR  = CFE * CXCMS
360             END IF
361             TXXOLD = TXX
362             TYYOLD = TYY
363             TZZOLD = TZZ
364             PCMSMX = PPRIM - ETACM * ( EPRIM - ETACM * PPRIM / ( GAMCM
365      &             + 1.D+00 ) )
366             RMAX2  = 1.D+00 + AMPRI2 / UMO2 - 2.D+00 * ( EPRIM -
367      &               ETACM * PCMSMX ) / ECHCK
368             IF ( RMAX2 .GT. RMIN2 ) THEN
369                PTRANS = 0.D+00
370             ELSE
371                PTRANS = 0.D+00
372             END IF
373  
374 1220        CONTINUE
375                PLONG2 = ( PPRIM - PTRANS ) * ( PPRIM + PTRANS )
376                PLONGI = SQRT ( PLONG2 )
377                PXLAST = PTRANS * CXXTR + PLONGI * CXCMS
378                PYLAST = PTRANS * CYYTR + PLONGI * CYCMS
379                PZLAST = PTRANS * CZZTR + PLONGI * CZCMS
380                PCMSLN = PLONGI - ETACM * ( EPRIM - ETACM * PLONGI /
381      &                ( GAMCM + 1.D+00 ) )
382                RURM2  = 1.D+00 + AMPRI2 / UMO2 - 2.D+00 * ( EPRIM -
383      &                  ETACM * PCMSLN ) / ECHCK
384                TXX = PXLAST / PPRIM
385                TYY = PYLAST / PPRIM
386                TZZ = PZLAST / PPRIM
387                IF ( RURM2 .LE. RMIN2 ) THEN
388                   PTRANS = 0.5D+00 * PTRANS
389                   WRITE ( LUNERR, * )' Ferevv: R2 < Rmin2 for Pt',
390      &                                 RURM2, RMIN2, URMIN2
391                   IF ( PTRANS .GT. ANGLGB ) GO TO 1220
392                END IF
393             PCMSX  = PCMSLN * CXCMS + PTRANS * CXXTR
394             PCMSY  = PCMSLN * CYCMS + PTRANS * CYYTR
395             PCMSZ  = PCMSLN * CZCMS + PTRANS * CZZTR
396             ERCMS  = 0.5D+00 * ( UMO2 * ( 1.D+00 + RURM2 ) - AMPRI2 )
397      &             / UMO
398             ETAPCM = PCMSLN * ETACM
399             ECHCK  = GAMCM * ERCMS - ETAPCM
400             PHELP  = - ETAPCM / (GAMCM + 1.D+00) + ERCMS
401             PXCHCK = - PCMSX + ETAX * PHELP
402             PYCHCK = - PCMSY + ETAY * PHELP
403             PZCHCK = - PCMSZ + ETAZ * PHELP
404             ECHCK = ECHCK - V0WELL (ITJ) - EKRECL - EBNDNG (ITJ)
405             UMO   = ECHCK**2 - PXCHCK**2 - PYCHCK**2
406      &            - PZCHCK**2
407             IF ( UMO .LT. 0.D+00 ) THEN
408                LRESMP = .TRUE.
409                RETURN
410             END IF
411             UMO   = SQRT ( UMO )
412          ELSE
413             EM = EKM + AMPROJ - V0WELL (ITJ) - EKRECL - EBNDNG (ITJ)
414             ECHCK  = EM + EFER
415             PXCHCK = PXNUC + PM * TXX
416             PYCHCK = PYNUC + PM * TYY
417             PZCHCK = PZNUC + PM * TZZ
418             UMO    = ECHCK**2 - PXCHCK**2 - PYCHCK**2
419      &             - PZCHCK**2
420             IF ( UMO .LT. 0.D+00 ) THEN
421                WRITE (LUNOUT,*)' *** Ferevv: Umo < 0 ',UMO
422                WRITE (LUNERR,*)' *** Ferevv: Umo < 0 ',UMO
423                LRESMP = .TRUE.
424                NHAD = 0
425                RETURN
426             END IF
427             UMO    = SQRT ( UMO )
428          END IF
429          EFRM  = EFRM  + EKFER - V0WELL (ITJ) - EKRECL
430          TVEUZ = TVEUZ + V0WELL (ITJ) - EKFER + EKRECL
431       END IF
432       UMO2   = UMO * UMO
433       EPROJX = HLFHLF * ( UMO2 - AMPROJ**2 - AMTAR**2 ) / AMTAR
434       IF ( EPROJX .LE. AMPROJ ) THEN
435          WRITE (LUNOUT,*)' Ferevv: Eprojx < Amproj after kin. sel. !',
436      &                     EPROJX, AMPROJ, LSEACL
437          WRITE (LUNERR,*)' Ferevv: Eprojx < Amproj after kin. sel. !',
438      &                     EPROJX, AMPROJ, LSEACL
439          PXFRM = PXFRM - PXNUC
440          PYFRM = PYFRM - PYNUC
441          PZFRM = PZFRM - PZNUC
442          EFRM  = EFRM  - EKFER + V0WELL (ITJ) + EKRECL
443          TVEUZ = TVEUZ - V0WELL (ITJ) + EKFER - EKRECL
444          IF ( LSEACL ) GO TO 50
445          LRESMP = .TRUE.
446          NHAD   = 0
447          RETURN
448       END IF
449       PPROJX = SQRT ( ( EPROJX - AMPROJ ) * ( EPROJX + AMPROJ ) )
450       ETOTX  = EPROJX + AMTAR
451       PTOSCA = PXCHCK * TXX + PYCHCK * TYY + PZCHCK * TZZ
452       PXTART = PXCHCK - PTOSCA * TXX
453       PYTART = PYCHCK - PTOSCA * TYY
454       PZTART = PZCHCK - PTOSCA * TZZ
455       PTRASQ = PXTART**2 + PYTART**2 + PZTART**2
456       AMTRSQ = AMTAR**2  + PTRASQ
457       UMOTR2 = UMO2 + PTRASQ
458       UMOTR  = SQRT (UMOTR2)
459       PPARSQ = ECHCK**2 - UMOTR2
460       PPARTT = SQRT (PPARSQ)
461       GAMCMS = ECHCK  / UMOTR
462       ETACMS = PPARTT / UMOTR
463       EPRCMS = HLFHLF * ( UMOTR2 + AMPROJ**2 - AMTRSQ ) / UMOTR
464       PPRCMS = SQRT ( ( EPRCMS - AMPROJ ) * ( EPRCMS + AMPROJ ) )
465       EPRLAB = GAMCMS * EPRCMS + ETACMS * PPRCMS
466       ETRLAB = ECHCK  - EPRLAB
467       PPRLAB = SQRT ( ( EPRLAB - AMPROJ ) * ( EPRLAB + AMPROJ ) )
468       PXTARG = PXCHCK - PPRLAB * TXX
469       PYTARG = PYCHCK - PPRLAB * TYY
470       PZTARG = PZCHCK - PPRLAB * TZZ
471       GAM    = ETRLAB / AMTAR
472       BGX    = PXTARG / AMTAR
473       BGY    = PYTARG / AMTAR
474       BGZ    = PZTARG / AMTAR
475       PPHELP = ( BGX * TXX + BGY * TYY + BGZ * TZZ ) * PPRLAB
476       ETAPCM = EPRLAB - PPHELP / ( GAM + ONEONE )
477       PXPROJ = PPRLAB * TXX - BGX * ETAPCM
478       PYPROJ = PPRLAB * TYY - BGY * ETAPCM
479       PZPROJ = PPRLAB * TZZ - BGZ * ETAPCM
480       UUOLD  = PXPROJ / PPROJX
481       VVOLD  = PYPROJ / PPROJX
482       WWOLD  = PZPROJ / PPROJX
483       SINT02 = UUOLD**2 + VVOLD**2
484       IF ( SINT02 .LE. ANGLSQ ) THEN
485          LSMPAN = .TRUE.
486          SINTH0 = ZERZER
487          COSPH0 = ONEONE
488          SINPH0 = ZERZER
489       ELSE
490          LSMPAN = .FALSE.
491          SINTH0 = SQRT (SINT02)
492          COSPH0 = UUOLD / SINTH0
493          SINPH0 = VVOLD / SINTH0
494       END IF
495       PLABS  = PPROJX
496       ELABS  = EPROJX
497       IF ( LSEACL .OR. .NOT. LDIFFR (KPTOIP(KP)) .OR. PLABS .LE. PTHDFF
498      &    ) THEN
499          FRUUN(1) = ONEONE
500       ELSE
501          IF ( RN1GSC .GE. ZERZER ) THEN
502             CALL GRNDM(FRUUN,1)
503             IF ( FRUUN(1) .LT. HLFHLF ) THEN
504                FRUUN(1) = RN1GSC
505             ELSE
506                FRUUN(1) = RN2GSC
507             END IF
508          ELSE
509             CALL GRNDM(FRUUN,1)
510          END IF
511       END IF
512  
513       IF ( UMO * UMO .LT. URMIN2 ) THEN
514          IF ( URMIN2 - UMO2 .LT. 0.1D+00 * URMIN2 .AND. IIBAR (KP) .LT.
515      &      0 ) GO TO 1550
516          IF ( .NOT. LSEACL .AND. PLABS .GT. 2.D+00 ) GO TO 1550
517          WRITE ( LUNOUT,* )' Ferevv: Umo2 < Urmin2 !!',UMO*UMO,URMIN2,
518      &                       KP,KT,AMPROJ,AMTAR
519          WRITE ( LUNERR,* )' Ferevv: Umo2 < Urmin2 !!',UMO*UMO,URMIN2,
520      &                       KP,KT,AMPROJ,AMTAR
521          NHAD = 2
522          NREH(1) = KT
523          PXH(1) = 0.D+00
524          PYH(1) = 0.D+00
525          PZH(1) = 0.D+00
526          HEPH(1) = AMTAR
527          AMH (1) = AMTAR
528          IBARH (1) = IIBAR (KT)
529          ICHH  (1) = IICH  (KT)
530          ANH   (1) = ANAME (KT)
531          IF ( LSEACL .OR. IVFLAG .EQ. 0 ) THEN
532             NREH(2) = 23
533          ELSE
534             NREH(2) = KP
535          END IF
536          PXH(2) = 0.D+00
537          PYH(2) = 0.D+00
538          PZH(2) = PLABS
539          HEPH(2) = ELABS
540          AMH (2) = AMPROJ
541          IBARH (2) = IIBAR (NREH(2))
542          ICHH  (2) = IICH  (NREH(2))
543          ANH   (2) = ANAME (NREH(2))
544          GO TO 1122
545       END IF
546  1550 CONTINUE
547       IF ( FRUUN(1) .LE. FRDIFF ) THEN
548          CALL DIFEVV ( NHAD, KP, KT, PLABS, ELABS, UMO )
549          LEVDIF = .TRUE.
550       ELSE IF ( .NOT. LSEACL .AND. PLABS .LT. DBLE(IJNUCR(KPTOIP(KP)))
551      &        * 0.8D+00 * PTHRSH (KPTOIP(KP)) ) THEN
552          IF ( KP .EQ. 26 ) THEN
553             KPP = 23
554          ELSE
555             KPP = KP
556          END IF
557          CALL HEVHIN ( NHAD, KPP, KT, PLABS, ELABS, UMO )
558          LHADRI = .TRUE.
559       ELSE
560          IF ( .NOT. LSEACL ) THEN
561             LEVDIF = .FALSE.
562             LHADRI = .FALSE.
563          END IF
564          CALL HADEVV ( NHAD, KP, KT, PLABS, ELABS, UMO )
565       END IF
566 1122  CONTINUE
567       DO 2000 I=1,NHAD
568          IF ( LSMPAN ) THEN
569             PZH (I) = WWOLD * PZH (I)
570          ELSE
571             PLRX = PXH (I) * COSPH0 * WWOLD - PYH (I) * SINPH0
572      &           + PZH (I) * UUOLD
573             PLRY = PXH (I) * SINPH0 * WWOLD + PYH (I) * COSPH0
574      &           + PZH (I) * VVOLD
575             PLRZ = - PXH (I) * SINTH0 + PZH (I) * WWOLD
576             PXH (I) = PLRX
577             PYH (I) = PLRY
578             PZH (I) = PLRZ
579          END IF
580          CALL ALTRA ( GAM, BGX, BGY, BGZ, PXH(I), PYH(I), PZH(I),
581      &                HEPH(I), PLR, PLRX, PLRY, PLRZ, ELR )
582          PXH(I) = PLRX
583          PYH(I) = PLRY
584          PZH(I) = PLRZ
585          HEPH(I) = ELR
586 2000  CONTINUE
587  
588       B1BAMJ = B1SAVE
589       B2BAMJ = B2SAVE
590       RETURN
591       END