[u/mrichter/AliRoot.git] / GEANT321 / fluka / difevv.F

*
* $Id$
*
* $Log$
* Revision 1.1.1.1  1995/10/24 10:19:54  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 29/03/94  15.41.42  by  S.Giani
*-- Author :
*$ CREATE DIFEVV.FOR
*COPY DIFEVV
*
*=== difevv ===========================================================*
*
      SUBROUTINE DIFEVV ( NHAD, KPROJ, KTARG, PPROJ, EPROJ, UMO )
 
#include "geant321/dblprc.inc"
#include "geant321/dimpar.inc"
#include "geant321/iounit.inc"
*----------------------------------------------------------------------*
C
C          DIFFRACTIVE HADRON -HADRON COLLISIONS
C     GENERATE HADRON PRODUCTION EVENT IN  KPROJ - KTARG  COLLISION
C     WITH  LAB PROJECTILE MOMENTUM  PPROJ
C
C********************************************************************
C
#include "geant321/auxpar.inc"
#include "geant321/balanc.inc"
#include "geant321/cmsres.inc"
#include "geant321/finpar.inc"
#include "geant321/hadpar.inc"
#include "geant321/inpdat2.inc"
#include "geant321/part.inc"
#include "geant321/qquark.inc"
      COMMON /FKPRIN/ IPRI, INIT
      REAL RNDM(2)
C
C*******************************************************************"
C
C     KINEMATICS
C
C********************************************************************
C
      IPRI = 0
      AMPROJ = AM(KPROJ)
      AMTAR  = AM(KTARG)
* The following are the Lorentz parameters to come from the system
* (projectile + target) rest frame to the starting one, which is the
* one where the target is at rest and the projectile is moving
* along the +z direction with Pproj: from now down to 600 continue
* we are working in the system rest frame !!!
      GAMCM = (EPROJ+AMTAR)/UMO
      BGCM  = PPROJ/UMO
C
*or      IF(IPRI.EQ.1) WRITE(LUNOUT,101)KPROJ,KTARG,PPROJ,AMPROJ,AMTAR,
*or     *EPROJ,UMO,GAMCM,BGCM
*or101   FORMAT(2I5,10F11.5)
C
C
      IBPROJ = IBAR(KPROJ)
      IBTARG = IBAR(KTARG)
C
C
C=====================================================================
C
C     SAMPLE X-VALUES OF QUARK-ANTIQUARK PAIRS
C
C======================================================================
      IF ( KPROJ .GT. 2 ) THEN
         UNOSEA = 5.D+00
      ELSE
         UNOSEA = 3.D+00
      END IF
* Come here if we need to resample xsea and xasea!!
  211 CONTINUE
      TMP005 = 0.05D+00
      XSEA  = BETARN(TMP005,UNOSEA)
      XASEA = BETARN(TMP005,UNOSEA)
      XPIO  = XSEA+XASEA
      IF (XPIO .GE. 1.D+00) GO TO 211
      XHAD = 1.D+00 - XPIO
      CALL GRNDM(RNDM,1)
      ISAM = 2.D+00 * RNDM(1) + 1.D+00
*or      IF (IPRI.EQ.1) WRITE(LUNOUT,371)XSEA,XASEA,XPIO,XHAD,ISAM
*or  371 FORMAT (' XSEA,XASEA,XPIO,XHAD,ISAM',4F10.5,I10)
C=====================================================================
C
C      CALL EQUIVALENT PIO HADRON COLLISIONS
C
C=====================================================================
      GO TO (250,260),ISAM
*  +-------------------------------------------------------------------*
*  | Target excited !!!!
  250 CONTINUE
         LPRDIF = .TRUE.
C=======================================================================
C     PROJECTILE MOVING WITH XHAD, TARGET EXCITED
C=======================================================================
         IIDIF = 1
         AMCH  = SQRT (XPIO) * UMO
         BITBIT = 0.5D+00
*  |  +----------------------------------------------------------------*
*  |  | The following condition roughly guarantees 1 GeV for the total
*  |  | energy of the pseudo-pion!
         IF ( AMCH .LE. AMTAR + BITBIT ) THEN
            AMCH = AMTAR + BITBIT
            XPIO = ( AMCH / UMO )**2
            XHAD = 1.D+00 - XPIO
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
*  | The following instructions make the division of the invariant
*  | mass of the system between the two particles, the two resulting
*  | energies being Eh1s and Eh2ex and the momentum Ph1s: the two
*  | particles are the original projectile and the excited target,
*  | ("mass" = amch)
         EH1S = (UMO**2 + AMPROJ**2 - AMCH**2) / (2.D0*UMO)
         IF (EH1S .LE. AMPROJ) GO TO 211
         EH2EX = UMO - EH1S
         PH1S  = SQRT (EH1S**2 - AMPROJ**2)
C*** AND  INT. CHAIN TRANSVERSE MOMENTA
         B3SAVE = B3BAMJ
         CALL GRNDM(RNDM,2)
         ES  = -2.D0/(B3BAMJ**2)*LOG(RNDM(1)*RNDM(2))
         HPS = SQRT(ES*ES+2.D0*ES*AMTAR)
         CALL SFECFE(SFE,CFE)
         PTXCH1 = HPS * CFE
         PTYCH1 = HPS * SFE
6171  CONTINUE
         PX1 = PTXCH1
         PY1 = PTYCH1
         ACH = PH1S**2 - PX1**2 - PY1**2
         IF (ACH .LE. 0.D+00) THEN
            PTXCH1 = 0.75D+00 * PTXCH1
            PTYCH1 = 0.75D+00 * PTYCH1
            GO TO 6171
         END IF
         PZ1  = SQRT (ACH)
*  | Now transform back the excited target to the lab system
         ECHCK  = GAMCM * EH2EX - BGCM * PZ1
         PXCHCK = - PX1
         PYCHCK = - PY1
         PZCHCK = - GAMCM * PZ1 + BGCM * EH2EX
*  | Now ..chck are the kinematical variables of the excited target
*  | in the lab frame, the invariant mass is always Amch
         CALL GRNDM(RNDM,1)
         IF (RNDM (1) .LE. 0.5D+00 ) THEN
            KPIO = 26
         ELSE
            KPIO = 23
         END IF
         AMPIO = AM (KPIO)
         EPIOL = 0.5D+00 * ( AMCH**2 - AMPIO**2 - AMTAR**2 ) / AMTAR
         PPIOL = SQRT ( EPIOL**2 - AMPIO**2 )
         ETOTX  = EPIOL + AMTAR
         AAFACT = ECHCK + ETOTX
         BBFACT = PPIOL - PZCHCK
         DDENOM = ETOTX * AAFACT - PPIOL * BBFACT
         GAM1 = ( ECHCK * AAFACT + PPIOL * BBFACT ) / DDENOM
         BGZ1 = - BBFACT * AAFACT / DDENOM
         BGX1 = PXCHCK * ( GAM1 + 1.D+00 ) / AAFACT
         BGY1 = PYCHCK * ( GAM1 + 1.D+00 ) / AAFACT
         CALL HADEVV ( NHAD, KPIO, KTARG, PPIOL, EPIOL, AMCH )
*   Restore the original b3bamj parameter
         B3BAMJ = B3SAVE
C     PRINT 888,PX1,PY1,PZ1,EH2EX
C        PRINT 888,PXX1,PYY1,PZZ1,EXXX
C 888    FORMAT(4F12.4)
*  | The following to go back to the original (lab) frame
*  |  +----------------------------------------------------------------*
*  |  |                      Looping over the produced particles
         DO 800 I=1,NHAD
            CALL ALTRA ( GAM1, BGX1, BGY1, BGZ1, PXH(I), PYH(I), PZH(I),
     &                   HEPH(I), PLR, PLRX, PLRY, PLRZ, ELR )
            PXH(I) = PLRX
            PYH(I) = PLRY
            PZH(I) = PLRZ
            HEPH(I) = ELR
*  |  |  Updating conservation counters
 800     CONTINUE
*  |
*  +-------------------------------------------------------------------*
*  |  Add the original projectile to the final particles, transforming
*  |  it back to the lab system
         NHAD = NHAD+1
         PXH  (NHAD) = PX1
         PYH  (NHAD) = PY1
         PZH  (NHAD) = GAMCM * PZ1 + BGCM * EH1S
         HEPH (NHAD) = GAMCM * EH1S + BGCM * PZ1
         AMH  (NHAD) = AMPROJ
         ICHH (NHAD) = ICH  (KPROJ)
         IBARH(NHAD) = IBAR (KPROJ)
         ANH  (NHAD) = ANAME(KPROJ)
         NREH (NHAD) = KPROJ
      GO TO 600
*  |  end of excited target treatment !!
*  +-------------------------------------------------------------------*
 
*  +-------------------------------------------------------------------*
*  |  Excited projectile !!!!!!
 260  CONTINUE
         LPRDIF = .FALSE.
C=======================================================================
C   THE TARGET PARTICLE GETS XHAD , THE PROJECTILE BECOMES EXCITED
C   WE GO TO THE PROJECTILE REST FRAME
C=======================================================================
         IIDIF = 2
         AMCH  = SQRT (XPIO) * UMO
         MK = 0
         DO 270 IQ = 1, 3
            MK = MK + ABS ( IQSCHR ( MQUARK ( IQ, KPTOIP(KPROJ) ) ) )
 270     CONTINUE
         BITBIT = 0.5D+00 + 0.2D+00 * MK
*  |  +----------------------------------------------------------------*
*  |  | The following condition roughly guarantees 1 GeV for the total
*  |  | energy of the pseudo-pion if the projectile has no strangeness
*  |  | a bit more if it has
         IF ( AMCH .LE. AMPROJ + BITBIT ) THEN
            AMCH = AMPROJ + BITBIT
            XPIO = ( AMCH / UMO )**2
            XHAD = 1.D+00 - XPIO
         END IF
*  |  |
*  |  +----------------------------------------------------------------*
*  | The following instructions make the division of the invariant
*  | mass of the system between the two particles, the two resulting
*  | energies being Eh1s and Eh2ex and the momentum Ph1s: the two
*  | particles are the target nucleon and the excited projectile,
*  | ("mass" = amch)
         EH1S = (UMO**2 + AMTAR**2 -  AMCH**2) / (2.D+00*UMO)
         IF (EH1S .LE. AMTAR) GO TO 211
         EH2EX = UMO-EH1S
         PH1S  = SQRT (EH1S**2 - AMTAR**2)
C*** AND  INT. CHAIN TRANSVERSE MOMENTA
         B3SAVE = B3BAMJ
         CALL GRNDM(RNDM,2)
         ES  = -2.D0/(B3BAMJ**2)*LOG(RNDM(1)*RNDM(2))
         HPS = SQRT(ES*ES+2.D0*ES*AMPROJ)
         CALL SFECFE(SFE,CFE)
         PTXCH1 = HPS * CFE
         PTYCH1 = HPS * SFE
6181  CONTINUE
         PX1 = PTXCH1
         PY1 = PTYCH1
         ACH = PH1S**2 - PX1**2 - PY1**2
         IF (ACH .LE. 0.D+00) THEN
            PTXCH1 = 0.75D+00 * PTXCH1
            PTYCH1 = 0.75D+00 * PTYCH1
            GO TO 6181
         END IF
         PZ1  = SQRT (ACH)
*  | Now transform back the excited projectile to the lab system
         ECHCK  = GAMCM * EH2EX + BGCM * PZ1
         PXCHCK = PX1
         PYCHCK = PY1
         PZCHCK = GAMCM * PZ1   + BGCM * EH2EX
*  | Now ..chck are the kinematical variables of the excited projectile
*  | in the lab frame, the invariant mass is always Amch
         CALL GRNDM(RNDM,1)
         IF (RNDM (1) .LE. 0.5D+00 ) THEN
            KPIO = 26
         ELSE
            KPIO = 23
         END IF
         AMPIO = AM (KPIO)
         EPIOL = 0.5D+00 * ( AMCH**2 - AMPIO**2 - AMPROJ**2 ) / AMPROJ
         PPIOL = SQRT ( EPIOL**2 - AMPIO**2 )
         ETOTX  = EPIOL + AMPROJ
         AAFACT = ECHCK + ETOTX
         BBFACT = PPIOL - PZCHCK
         DDENOM = ETOTX * AAFACT - PPIOL * BBFACT
         GAM1 = ( ECHCK * AAFACT + PPIOL * BBFACT ) / DDENOM
         BGZ1 = - BBFACT * AAFACT / DDENOM
         BGX1 = PXCHCK * ( GAM1 + 1.D+00 ) / AAFACT
         BGY1 = PYCHCK * ( GAM1 + 1.D+00 ) / AAFACT
         CALL HADEVV ( NHAD, KPIO, KPROJ, PPIOL, EPIOL, AMCH )
*   Restore the original b3bamj parameter
         B3BAMJ = B3SAVE
C     PRINT 888,PX1,PY1,PZ1,EH2EX
C        PRINT 888,PXX1,PYY1,PZZ1,EXXX
*  | The following to go back to the original (lab) frame
*  |  +----------------------------------------------------------------*
*  |  |                      Looping over the produced particles
         DO 900 I=1,NHAD
            CALL ALTRA ( GAM1, BGX1, BGY1, BGZ1, PXH(I), PYH(I), PZH(I),
     &                   HEPH(I), PLR, PLRX, PLRY, PLRZ, ELR )
            PXH(I) = PLRX
            PYH(I) = PLRY
            PZH(I) = PLRZ
            HEPH(I) = ELR
*  |  |  Updating conservation counters
 900     CONTINUE
*  |
*  +-------------------------------------------------------------------*
*  |  Add the target nucleon to the final particles, transforming
*  |  it back to the lab system
         NHAD = NHAD + 1
         PXH  (NHAD) = - PX1
         PYH  (NHAD) = - PY1
         PZH  (NHAD) = - GAMCM * PZ1  + BGCM * EH1S
         HEPH (NHAD) =   GAMCM * EH1S - BGCM * PZ1
         AMH  (NHAD) = AMTAR
         ICHH (NHAD) = ICH  (KTARG)
         IBARH(NHAD) = IBAR (KTARG)
         ANH  (NHAD) = ANAME(KTARG)
         NREH (NHAD) = KTARG
      GO TO 600
*  |  end of excited projectile !!!
*  +-------------------------------------------------------------------*
  600 CONTINUE
C
C********************************************************************
C
C*** PRINT AND TEST ENERGY CONSERVATION
C
C********************************************************************
C
      PUZZ = 0.D+00
      EUZZ = 0.D+00
      PUXX = 0.D+00
      PUYY = 0.D+00
      ICUU = 0
      IBUU = 0
      DO 82 I=1,NHAD
         PUXX = PUXX + PXH(I)
         PUYY = PUYY + PYH(I)
         PUZZ = PUZZ + PZH(I)
         EUZZ = EUZZ + HEPH(I)
         ICUU = ICUU + ICHH(I)
         IBUU = IBUU + IBARH(I)
  82  CONTINUE
      ICHTOT=ICH(KPROJ)+ICH(KTARG)
      IBTOT =IBAR(KPROJ)+IBAR(KTARG)
      PCHMIN = 1.D-10 * PPROJ
      IF ((ABS(PUXX) .GE. PCHMIN) .OR. (ABS(PUYY) .GE. PCHMIN) .OR.
     &    (ABS(PUZZ-PPROJ) .GE. PCHMIN) .OR. (ABS(EPROJ+AMTAR-EUZZ) .GE.
     &     1.D-10*EUZZ) .OR. (ICHTOT .NE. ICUU) .OR. (IBTOT .NE. IBUU))
     &   THEN
         WRITE(LUNERR,*)
     &               ' Difevt: failure!!!: NHAD, KPROJ, KTARG, IIDIF: ',
     &                 NHAD, KPROJ, KTARG, IIDIF
         WRITE(LUNERR,*)'                     ',
     &                  'ICHTOT, ICUU, IBTOT, IBUU: ',
     &                   ICHTOT, ICUU, IBTOT, IBUU
         WRITE(LUNERR,*)'                     ',
     &                  'PPROJ, PUXX, PUYY, PUZZ: ',
     &                   PPROJ, PUXX, PUYY, PUZZ
         WRITE(LUNERR,*)'                     EPROJ, EUZZ: ',
     &                   EPROJ, EUZZ
      END IF
      IF (IPRI .NE. 1) GO TO 90
      DO 84 I=1,NHAD
         WRITE(LUNERR,85)I,NREH(I),ICHH(I),IBARH(I),ANH(I),
     &                   PXH(I),PYH(I),PZH(I),HEPH(I),AMH(I)
  85     FORMAT (4I5,A8,5F12.6)
  84  CONTINUE
  90  CONTINUE
      RETURN
      END
Commit	Line	Data
fe4da5cc	1	*
	2	* $Id$
	3	*
	4	* $Log$
	5	* Revision 1.1.1.1 1995/10/24 10:19:54 cernlib
	6	* Geant
	7	*
	8	*
	9	#include "geant321/pilot.h"
	10	*CMZ : 3.21/02 29/03/94 15.41.42 by S.Giani
	11	*-- Author :
	12	*$ CREATE DIFEVV.FOR
	13	*COPY DIFEVV
	14	*
	15	=== difevv ===========================================================
	16	*
	17	SUBROUTINE DIFEVV ( NHAD, KPROJ, KTARG, PPROJ, EPROJ, UMO )
	18
	19	#include "geant321/dblprc.inc"
	20	#include "geant321/dimpar.inc"
	21	#include "geant321/iounit.inc"
	22	----------------------------------------------------------------------
	23	C
	24	C DIFFRACTIVE HADRON -HADRON COLLISIONS
	25	C GENERATE HADRON PRODUCTION EVENT IN KPROJ - KTARG COLLISION
	26	C WITH LAB PROJECTILE MOMENTUM PPROJ
	27	C
	28	C********************************************************************
	29	C
	30	#include "geant321/auxpar.inc"
	31	#include "geant321/balanc.inc"
	32	#include "geant321/cmsres.inc"
	33	#include "geant321/finpar.inc"
	34	#include "geant321/hadpar.inc"
	35	#include "geant321/inpdat2.inc"
	36	#include "geant321/part.inc"
	37	#include "geant321/qquark.inc"
	38	COMMON /FKPRIN/ IPRI, INIT
	39	REAL RNDM(2)
	40	C
	41	C*******************************************************************"
	42	C
	43	C KINEMATICS
	44	C
	45	C********************************************************************
	46	C
	47	IPRI = 0
	48	AMPROJ = AM(KPROJ)
	49	AMTAR = AM(KTARG)
	50	* The following are the Lorentz parameters to come from the system
	51	* (projectile + target) rest frame to the starting one, which is the
	52	* one where the target is at rest and the projectile is moving
	53	* along the +z direction with Pproj: from now down to 600 continue
	54	* we are working in the system rest frame !!!
	55	GAMCM = (EPROJ+AMTAR)/UMO
	56	BGCM = PPROJ/UMO
	57	C
	58	*or IF(IPRI.EQ.1) WRITE(LUNOUT,101)KPROJ,KTARG,PPROJ,AMPROJ,AMTAR,
	59	or EPROJ,UMO,GAMCM,BGCM
	60	*or101 FORMAT(2I5,10F11.5)
	61	C
	62	C
	63	IBPROJ = IBAR(KPROJ)
	64	IBTARG = IBAR(KTARG)
65	C
66	C
67	C=====================================================================
68	C
69	C SAMPLE X-VALUES OF QUARK-ANTIQUARK PAIRS
70	C
71	C======================================================================
72	IF ( KPROJ .GT. 2 ) THEN
73	UNOSEA = 5.D+00
74	ELSE
75	UNOSEA = 3.D+00
76	END IF
77	* Come here if we need to resample xsea and xasea!!
78	211 CONTINUE
79	TMP005 = 0.05D+00
80	XSEA = BETARN(TMP005,UNOSEA)
81	XASEA = BETARN(TMP005,UNOSEA)
82	XPIO = XSEA+XASEA
83	IF (XPIO .GE. 1.D+00) GO TO 211
84	XHAD = 1.D+00 - XPIO
85	CALL GRNDM(RNDM,1)
86	ISAM = 2.D+00 * RNDM(1) + 1.D+00
87	*or IF (IPRI.EQ.1) WRITE(LUNOUT,371)XSEA,XASEA,XPIO,XHAD,ISAM
88	*or 371 FORMAT (' XSEA,XASEA,XPIO,XHAD,ISAM',4F10.5,I10)
89	C=====================================================================
90	C
91	C CALL EQUIVALENT PIO HADRON COLLISIONS
92	C
93	C=====================================================================
94	GO TO (250,260),ISAM
95	* +-------------------------------------------------------------------*
96	* \| Target excited !!!!
97	250 CONTINUE
98	LPRDIF = .TRUE.
99	C=======================================================================
100	C PROJECTILE MOVING WITH XHAD, TARGET EXCITED
101	C=======================================================================
102	IIDIF = 1
103	AMCH = SQRT (XPIO) * UMO
104	BITBIT = 0.5D+00
105	* \| +----------------------------------------------------------------*
106	* \| \| The following condition roughly guarantees 1 GeV for the total
107	* \| \| energy of the pseudo-pion!
108	IF ( AMCH .LE. AMTAR + BITBIT ) THEN
109	AMCH = AMTAR + BITBIT
110	XPIO = ( AMCH / UMO )**2
111	XHAD = 1.D+00 - XPIO
112	END IF
113	* \| \|
114	* \| +----------------------------------------------------------------*
115	* \| The following instructions make the division of the invariant
116	* \| mass of the system between the two particles, the two resulting
117	* \| energies being Eh1s and Eh2ex and the momentum Ph1s: the two
118	* \| particles are the original projectile and the excited target,
119	* \| ("mass" = amch)
120	EH1S = (UMO2 + AMPROJ2 - AMCH*2) / (2.D0UMO)
121	IF (EH1S .LE. AMPROJ) GO TO 211
122	EH2EX = UMO - EH1S
123	PH1S = SQRT (EH1S2 - AMPROJ2)
124	C*** AND INT. CHAIN TRANSVERSE MOMENTA
125	B3SAVE = B3BAMJ
126	CALL GRNDM(RNDM,2)
127	ES = -2.D0/(B3BAMJ*2)LOG(RNDM(1)*RNDM(2))
128	HPS = SQRT(ESES+2.D0ES*AMTAR)
129	CALL SFECFE(SFE,CFE)
130	PTXCH1 = HPS * CFE
131	PTYCH1 = HPS * SFE
132	6171 CONTINUE
133	PX1 = PTXCH1
134	PY1 = PTYCH1
135	ACH = PH1S2 - PX12 - PY1**2
136	IF (ACH .LE. 0.D+00) THEN
137	PTXCH1 = 0.75D+00 * PTXCH1
138	PTYCH1 = 0.75D+00 * PTYCH1
139	GO TO 6171
140	END IF
141	PZ1 = SQRT (ACH)
142	* \| Now transform back the excited target to the lab system
143	ECHCK = GAMCM * EH2EX - BGCM * PZ1
144	PXCHCK = - PX1
145	PYCHCK = - PY1
146	PZCHCK = - GAMCM * PZ1 + BGCM * EH2EX
147	* \| Now ..chck are the kinematical variables of the excited target
148	* \| in the lab frame, the invariant mass is always Amch
149	CALL GRNDM(RNDM,1)
150	IF (RNDM (1) .LE. 0.5D+00 ) THEN
151	KPIO = 26
152	ELSE
153	KPIO = 23
154	END IF
155	AMPIO = AM (KPIO)
156	EPIOL = 0.5D+00 * ( AMCH2 - AMPIO2 - AMTAR**2 ) / AMTAR
157	PPIOL = SQRT ( EPIOL2 - AMPIO2 )
158	ETOTX = EPIOL + AMTAR
159	AAFACT = ECHCK + ETOTX
160	BBFACT = PPIOL - PZCHCK
161	DDENOM = ETOTX * AAFACT - PPIOL * BBFACT
162	GAM1 = ( ECHCK * AAFACT + PPIOL * BBFACT ) / DDENOM
163	BGZ1 = - BBFACT * AAFACT / DDENOM
164	BGX1 = PXCHCK * ( GAM1 + 1.D+00 ) / AAFACT
165	BGY1 = PYCHCK * ( GAM1 + 1.D+00 ) / AAFACT
166	CALL HADEVV ( NHAD, KPIO, KTARG, PPIOL, EPIOL, AMCH )
167	* Restore the original b3bamj parameter
168	B3BAMJ = B3SAVE
169	C PRINT 888,PX1,PY1,PZ1,EH2EX
170	C PRINT 888,PXX1,PYY1,PZZ1,EXXX
171	C 888 FORMAT(4F12.4)
172	* \| The following to go back to the original (lab) frame
173	* \| +----------------------------------------------------------------*
174	* \| \| Looping over the produced particles
175	DO 800 I=1,NHAD
176	CALL ALTRA ( GAM1, BGX1, BGY1, BGZ1, PXH(I), PYH(I), PZH(I),
177	& HEPH(I), PLR, PLRX, PLRY, PLRZ, ELR )
178	PXH(I) = PLRX
179	PYH(I) = PLRY
180	PZH(I) = PLRZ
181	HEPH(I) = ELR
182	* \| \| Updating conservation counters
183	800 CONTINUE
184	* \|
185	* +-------------------------------------------------------------------*
186	* \| Add the original projectile to the final particles, transforming
187	* \| it back to the lab system
188	NHAD = NHAD+1
189	PXH (NHAD) = PX1
190	PYH (NHAD) = PY1
191	PZH (NHAD) = GAMCM * PZ1 + BGCM * EH1S
192	HEPH (NHAD) = GAMCM * EH1S + BGCM * PZ1
193	AMH (NHAD) = AMPROJ
194	ICHH (NHAD) = ICH (KPROJ)
195	IBARH(NHAD) = IBAR (KPROJ)
196	ANH (NHAD) = ANAME(KPROJ)
197	NREH (NHAD) = KPROJ
198	GO TO 600
199	* \| end of excited target treatment !!
200	* +-------------------------------------------------------------------*
201
202	* +-------------------------------------------------------------------*
203	* \| Excited projectile !!!!!!
204	260 CONTINUE
205	LPRDIF = .FALSE.
206	C=======================================================================
207	C THE TARGET PARTICLE GETS XHAD , THE PROJECTILE BECOMES EXCITED
208	C WE GO TO THE PROJECTILE REST FRAME
209	C=======================================================================
210	IIDIF = 2
211	AMCH = SQRT (XPIO) * UMO
212	MK = 0
213	DO 270 IQ = 1, 3
214	MK = MK + ABS ( IQSCHR ( MQUARK ( IQ, KPTOIP(KPROJ) ) ) )
215	270 CONTINUE
216	BITBIT = 0.5D+00 + 0.2D+00 * MK
217	* \| +----------------------------------------------------------------*
218	* \| \| The following condition roughly guarantees 1 GeV for the total
219	* \| \| energy of the pseudo-pion if the projectile has no strangeness
220	* \| \| a bit more if it has
221	IF ( AMCH .LE. AMPROJ + BITBIT ) THEN
222	AMCH = AMPROJ + BITBIT
223	XPIO = ( AMCH / UMO )**2
224	XHAD = 1.D+00 - XPIO
225	END IF
226	* \| \|
227	* \| +----------------------------------------------------------------*
228	* \| The following instructions make the division of the invariant
229	* \| mass of the system between the two particles, the two resulting
230	* \| energies being Eh1s and Eh2ex and the momentum Ph1s: the two
231	* \| particles are the target nucleon and the excited projectile,
232	* \| ("mass" = amch)
233	EH1S = (UMO2 + AMTAR2 - AMCH*2) / (2.D+00UMO)
234	IF (EH1S .LE. AMTAR) GO TO 211
235	EH2EX = UMO-EH1S
236	PH1S = SQRT (EH1S2 - AMTAR2)
237	C*** AND INT. CHAIN TRANSVERSE MOMENTA
238	B3SAVE = B3BAMJ
239	CALL GRNDM(RNDM,2)
240	ES = -2.D0/(B3BAMJ*2)LOG(RNDM(1)*RNDM(2))
241	HPS = SQRT(ESES+2.D0ES*AMPROJ)
242	CALL SFECFE(SFE,CFE)
243	PTXCH1 = HPS * CFE
244	PTYCH1 = HPS * SFE
245	6181 CONTINUE
246	PX1 = PTXCH1
247	PY1 = PTYCH1
248	ACH = PH1S2 - PX12 - PY1**2
249	IF (ACH .LE. 0.D+00) THEN
250	PTXCH1 = 0.75D+00 * PTXCH1
251	PTYCH1 = 0.75D+00 * PTYCH1
252	GO TO 6181
253	END IF
254	PZ1 = SQRT (ACH)
255	* \| Now transform back the excited projectile to the lab system
256	ECHCK = GAMCM * EH2EX + BGCM * PZ1
257	PXCHCK = PX1
258	PYCHCK = PY1
259	PZCHCK = GAMCM * PZ1 + BGCM * EH2EX
260	* \| Now ..chck are the kinematical variables of the excited projectile
261	* \| in the lab frame, the invariant mass is always Amch
262	CALL GRNDM(RNDM,1)
263	IF (RNDM (1) .LE. 0.5D+00 ) THEN
264	KPIO = 26
265	ELSE
266	KPIO = 23
267	END IF
268	AMPIO = AM (KPIO)
269	EPIOL = 0.5D+00 * ( AMCH2 - AMPIO2 - AMPROJ**2 ) / AMPROJ
270	PPIOL = SQRT ( EPIOL2 - AMPIO2 )
271	ETOTX = EPIOL + AMPROJ
272	AAFACT = ECHCK + ETOTX
273	BBFACT = PPIOL - PZCHCK
274	DDENOM = ETOTX * AAFACT - PPIOL * BBFACT
275	GAM1 = ( ECHCK * AAFACT + PPIOL * BBFACT ) / DDENOM
276	BGZ1 = - BBFACT * AAFACT / DDENOM
277	BGX1 = PXCHCK * ( GAM1 + 1.D+00 ) / AAFACT
278	BGY1 = PYCHCK * ( GAM1 + 1.D+00 ) / AAFACT
279	CALL HADEVV ( NHAD, KPIO, KPROJ, PPIOL, EPIOL, AMCH )
280	* Restore the original b3bamj parameter
281	B3BAMJ = B3SAVE
282	C PRINT 888,PX1,PY1,PZ1,EH2EX
283	C PRINT 888,PXX1,PYY1,PZZ1,EXXX
284	* \| The following to go back to the original (lab) frame
285	* \| +----------------------------------------------------------------*
286	* \| \| Looping over the produced particles
287	DO 900 I=1,NHAD
288	CALL ALTRA ( GAM1, BGX1, BGY1, BGZ1, PXH(I), PYH(I), PZH(I),
289	& HEPH(I), PLR, PLRX, PLRY, PLRZ, ELR )
290	PXH(I) = PLRX
291	PYH(I) = PLRY
292	PZH(I) = PLRZ
293	HEPH(I) = ELR
294	* \| \| Updating conservation counters
295	900 CONTINUE
296	* \|
297	* +-------------------------------------------------------------------*
298	* \| Add the target nucleon to the final particles, transforming
299	* \| it back to the lab system
300	NHAD = NHAD + 1
301	PXH (NHAD) = - PX1
302	PYH (NHAD) = - PY1
303	PZH (NHAD) = - GAMCM * PZ1 + BGCM * EH1S
304	HEPH (NHAD) = GAMCM * EH1S - BGCM * PZ1
305	AMH (NHAD) = AMTAR
306	ICHH (NHAD) = ICH (KTARG)
307	IBARH(NHAD) = IBAR (KTARG)
308	ANH (NHAD) = ANAME(KTARG)
309	NREH (NHAD) = KTARG
310	GO TO 600
311	* \| end of excited projectile !!!
312	* +-------------------------------------------------------------------*
313	600 CONTINUE
314	C
315	C********************************************************************
316	C
317	C*** PRINT AND TEST ENERGY CONSERVATION
318	C
319	C********************************************************************
320	C
321	PUZZ = 0.D+00
322	EUZZ = 0.D+00
323	PUXX = 0.D+00
324	PUYY = 0.D+00
325	ICUU = 0
326	IBUU = 0
327	DO 82 I=1,NHAD
328	PUXX = PUXX + PXH(I)
329	PUYY = PUYY + PYH(I)
330	PUZZ = PUZZ + PZH(I)
331	EUZZ = EUZZ + HEPH(I)
332	ICUU = ICUU + ICHH(I)
333	IBUU = IBUU + IBARH(I)
334	82 CONTINUE
335	ICHTOT=ICH(KPROJ)+ICH(KTARG)
336	IBTOT =IBAR(KPROJ)+IBAR(KTARG)
337	PCHMIN = 1.D-10 * PPROJ
338	IF ((ABS(PUXX) .GE. PCHMIN) .OR. (ABS(PUYY) .GE. PCHMIN) .OR.
339	& (ABS(PUZZ-PPROJ) .GE. PCHMIN) .OR. (ABS(EPROJ+AMTAR-EUZZ) .GE.
340	& 1.D-10*EUZZ) .OR. (ICHTOT .NE. ICUU) .OR. (IBTOT .NE. IBUU))
341	& THEN
342	WRITE(LUNERR,*)
343	& ' Difevt: failure!!!: NHAD, KPROJ, KTARG, IIDIF: ',
344	& NHAD, KPROJ, KTARG, IIDIF
345	WRITE(LUNERR,*)' ',
346	& 'ICHTOT, ICUU, IBTOT, IBUU: ',
347	& ICHTOT, ICUU, IBTOT, IBUU
348	WRITE(LUNERR,*)' ',
349	& 'PPROJ, PUXX, PUYY, PUZZ: ',
350	& PPROJ, PUXX, PUYY, PUZZ
351	WRITE(LUNERR,*)' EPROJ, EUZZ: ',
352	& EPROJ, EUZZ
353	END IF
354	IF (IPRI .NE. 1) GO TO 90
355	DO 84 I=1,NHAD
356	WRITE(LUNERR,85)I,NREH(I),ICHH(I),IBARH(I),ANH(I),
357	& PXH(I),PYH(I),PZH(I),HEPH(I),AMH(I)
358	85 FORMAT (4I5,A8,5F12.6)
359	84 CONTINUE
360	90 CONTINUE
361	RETURN
362	END