This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / fluka / difevv.F
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 = (UMO**2 + AMPROJ**2 - AMCH**2) / (2.D0*UMO)
121          IF (EH1S .LE. AMPROJ) GO TO 211
122          EH2EX = UMO - EH1S
123          PH1S  = SQRT (EH1S**2 - AMPROJ**2)
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(ES*ES+2.D0*ES*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 = PH1S**2 - PX1**2 - 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 * ( AMCH**2 - AMPIO**2 - AMTAR**2 ) / AMTAR
157          PPIOL = SQRT ( EPIOL**2 - AMPIO**2 )
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 = (UMO**2 + AMTAR**2 -  AMCH**2) / (2.D+00*UMO)
234          IF (EH1S .LE. AMTAR) GO TO 211
235          EH2EX = UMO-EH1S
236          PH1S  = SQRT (EH1S**2 - AMTAR**2)
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(ES*ES+2.D0*ES*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 = PH1S**2 - PX1**2 - 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 * ( AMCH**2 - AMPIO**2 - AMPROJ**2 ) / AMPROJ
270          PPIOL = SQRT ( EPIOL**2 - AMPIO**2 )
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