Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
dpm25diff.f
Go to the documentation of this file.
1 *
2 *===sdiff================================================================*
3 *
4  SUBROUTINE sdiff(EPROJ,PPROJ,KPROJ,NHKKH1,IQQDD)
5 
6 **************************************************************************
7 * Version November 1993 by Stefan Roesler *
8 * University of Leipzig *
9 * This subroutine calls one single-diffractive event depending on the *
10 * single-diffractive cross section for a given hadron-hadron interaction.*
11 **************************************************************************
12 
13  IMPLICIT DOUBLE PRECISION (a-h,o-z)
14  SAVE
15  parameter(lout=6,llook=9)
16  parameter(nmxhkk= 89998,intmx=2488)
17  CHARACTER*8 aname,anc
18  COMMON /diqi/ ipvq(248), ippv1(248), ippv2(248), itvq(248),
19  & ittv1(248), ittv2(248), ipsq(intmx),ipsq2(intmx),
20  & ipsaq(intmx),ipsaq2(intmx),itsq(intmx),itsq2(intmx),
21  & itsaq(intmx),itsaq2(intmx),kkproj(248),kktarg(248)
22  COMMON /hkkevt/ nhkk,nevhkk, isthkk(nmxhkk), idhkk(nmxhkk),
23  & jmohkk(2,nmxhkk),jdahkk(2,nmxhkk),phkk(5,nmxhkk),
24  & vhkk(4,nmxhkk), whkk(4,nmxhkk)
25  COMMON /diffra/ isingd,idiftp,ioudif,iflagd
26  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
27  COMMON /cmhico/ cmhis
28  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
29  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
30  & iibar(210),k1(210),k2(210)
31  COMMON /difpar/ pxc(902), pyc(902),pzc(902),
32  & hec(902), amc(902),ichc(902),
33  & ibarc(902),anc(902),nrc(902)
34 *----------------- S. Roesler 07/07/93
35 C COMMON /COUNTEV/NCDIFF,NCSDIF
36 *---------- S. Roesler 5-11-93
37  COMMON /ifragm/ifrag
38  common/xxlmdd/ijlmdd,kdlmdd
39  COMMON /nuccms/gamv,bglv,dummy(7)
40 
41 *
42 
43  DATA ncirej /0/
44 C DATA NCDIFF /0/
45  DATA sigdif,sigdih /6.788790702d0,3.631283998d0/
46 *
47  ijlmdd= 0
48  irej = 0
49  ktarg = 0
50  pprox = 0.0d0
51  pproy = 0.0d0
52  pprol = pproj
53  epro = eproj
54  ampro = aam(kproj)
55  IF(ipri.GE.1)WRITE(6,'(A,2E20.8)')' SDIFF:EPROJ,PPROJ ',
56  * eproj,pproj
57 C WRITE(6,'(A,2E20.8)')' SDIFF:EPROJ,PPROJ ',
58 C * EPROJ,PPROJ
59 *
60  DO 10 i=1,it
61  ihkk = i+ip
62  IF (isthkk(ihkk).EQ.12) THEN
63  ktarg = kktarg(i)
64  ptarx = phkk(1,ihkk)
65  ptary = phkk(2,ihkk)
66  ptarl = phkk(3,ihkk)
67  etar = phkk(4,ihkk)
68  amtar = phkk(5,ihkk)
69  ptotx = pprox + ptarx
70  ptoty = pproy + ptary
71  ptotl = pprol + ptarl
72  etot = epro + etar
73  ptot = sqrt(ptotx**2+ptoty**2+ptotl**2)
74  amtot = sqrt(abs(etot-ptot)*(etot+ptot))
75  gam = etot /(amtot)
76  bgx = ptotx/(amtot)
77  bgy = ptoty/(amtot)
78  bgl = ptotl/(amtot)
79  IF (ipev.GE.6) THEN
80  WRITE(lout,1000) pprox,pproy,pprol,epro,ampro,kproj
81  WRITE(lout,1001) ptarx,ptary,ptarl,etar,amtar,ktarg
82  WRITE(lout,1011) ptotx,ptoty,ptotl,etot,amtot,ktarg
83  WRITE(lout,1002) amtot,gam,bgx,bgy,bgl
84  WRITE(lout,1702) gamv,bglv
85  1000 FORMAT('SDIFF: PPROX,PPROY,PPROL,EPRO,AMPRO,KPROJ',
86  & 2e15.5,2e15.5,e15.6,i2)
87  1001 FORMAT('SDIFF: PTARX,PTARY,PTARL,ETAR,AMTAR,KTARG',
88  & 5e15.6,i2)
89  1011 FORMAT('SDIFF: PTOTX,PTOTY,PTOTL,ETOT,AMTOT,KTARG',
90  & 5e15.6,i2)
91  1002 FORMAT('SDIFF: AMTOT,GAM,BGX,BGY,BGL',5e15.6)
92  1702 FORMAT('SDIFF: GAMV,BGLV',2e15.6)
93  ENDIF
94 C This is the hadron-hadron cms mot the overall cms
95 C (due to Fermi momenta)
96  CALL daltra(gam,-bgx,-bgy,-bgl,pprox,pproy,pprol,epro,
97  & ppcm,ppcmx,ppcmy,ppcml,epcm)
98  CALL daltra(gam,-bgx,-bgy,-bgl,ptarx,ptary,ptarl,etar,
99  & ptcm,ptcmx,ptcmy,ptcml,etcm)
100 C Back to lab frame
101  CALL daltra(gam,bgx,bgy,bgl,ppcmx,ppcmy,ppcml,epcm,
102  & ppla,pplax,pplay,pplal,epla)
103  CALL daltra(gam,bgx,bgy,bgl,ptcmx,ptcmy,ptcml,etcm,
104  & ptla,ptlax,ptlay,ptlal,etla)
105 C And now into the real cms frame
106  epcms=gamv*epla-bglv*pplal
107  pplcms=gamv*pplal-bglv*epla
108  etcms=gamv*etla-bglv*ptlal
109  ptlcms=gamv*ptlal-bglv*etla
110  IF (ipev.GE.6) THEN
111  WRITE(lout,1003) ppcm,ppcmx,ppcmy,ppcml,epcm
112  WRITE(lout,1004) ptcm,ptcmx,ptcmy,ptcml,etcm
113  1003 FORMAT('SDIFF: PPCM,PPCMX,PPCMY,PPCML,EPCM',5e15.5)
114  1004 FORMAT('SDIFF: PTCM,PTCMX,PTCMY,PTCML,ETCM',5e15.5)
115  WRITE(lout,1703) ppla,pplax,pplay,pplal,epla
116  WRITE(lout,1704) ptla,ptlax,ptlay,ptlal,etla
117  1703 FORMAT('SDIFF: PPLA,PPLAX,PPLAY,PPLAL,EPLA',5e15.5)
118  1704 FORMAT('SDIFF: PTLA,PTLAX,PTLAY,PTLAL,ETLA',5e15.5)
119  WRITE(lout,1803) ppcm,ppcmx,ppcmy,pplcms,epcms
120  WRITE(lout,1804) ptcm,ptcmx,ptcmy,ptlcms,etcms
121  1803 FORMAT('SDIFF: PPCM,PPCMX,PPCMY,PPLCMS,EPCMS',5e15.5)
122  1804 FORMAT('SDIFF: PTCM,PTCMX,PTCMY,PTLCMS,ETCMS',5e15.5)
123  ENDIF
124  cod = ppcml/ppcm
125  cod2 = cod**2
126  IF (cod2.GT.0.999999999999d0) cod2 = 0.999999999999d0
127  sid = sqrt((1.0d0-cod)*(1.0d0+cod))
128  cof = 1.0d0
129  sif = 0.0d0
130  IF (ppcm*sid.GT.1.0d-9) THEN
131  cof = ppcmx/(sid*ppcm)
132  sif = ppcmy/(sid*ppcm)
133  anorf = sqrt(cof**2+sif**2)
134  cof = cof/anorf
135  sif = sif/anorf
136  ENDIF
137  IF (ipev.GE.6) THEN
138  WRITE(lout,1005) cod,sid,cof,sif
139  1005 FORMAT('SDIFF: COD,SID,COF,SIF',4e15.5)
140  ENDIF
141 C This is the hadron-hadron cms mot the overall cms
142 C (due to Fermi momenta)
143  ecm = epcm + etcm
144  IF (ipev.GE.6) THEN
145  WRITE(lout,1705)ecm,epcm,etcm
146  1705 FORMAT('SDIFF:ECM,EPCM,ETCM',3e15.5)
147  ENDIF
148 C------------------------------------------------------------------
149 C
150 C Test J.R.2/94
151 C
152 C------------------------------------------------------------------
153  CALL sihndi(ecm,kproj,ktarg,sigdif,sigdih)
154 C------------------------------------------------------------------
155 C
156 C Test J.R.2/94
157 C
158 C------------------------------------------------------------------
159  fakk=1.9
160 C ------------------------------------------------------------------
161 C
162 C further modification j.r.6.1.95
163 C
164 C--------------------------------------------------------------------
165  IF(ecm.LE.10.d0)THEN
166  fakk=1.d0
167  ELSEIF(ecm.GE.30.d0)THEN
168  fakk=1.9d0
169  ELSE
170  fak=(ecm-10.d0)/20.d0
171  fakk=1.d0+fak*0.9d0
172  ENDIF
173  aitt=it
174  IF((isingd.LE.2).AND.(kproj.EQ.1.OR.kproj.EQ.8))THEN
175  aadiff=fakk*aitt**0.17d0
176  sigdif=aadiff*sigdif
177  ELSEIF((isingd.LE.2).AND.
178  * ((kproj.EQ.13).OR.(kproj.EQ.14).OR.(kproj.EQ.23)))THEN
179  aadiff=fakk*aitt**0.15d0
180  sigdif=aadiff*sigdif
181  ELSEIF((isingd.LE.2).AND.
182  * ((kproj.EQ.15).OR.(kproj.EQ.16).OR.
183  * (kproj.EQ.24).OR.(kproj.EQ.25)))THEN
184  aadiff=fakk*aitt**0.13d0
185  sigdif=aadiff*sigdif
186  ENDIF
187 C------------------------------------------------------------------
188 C
189  sigin = dshnto(kproj,ktarg,ecm)-dshnel(kproj,ktarg,ecm)
190  IF (ipev.GE.6) THEN
191  WRITE(lout,1060)kproj,ktarg,ecm,sigdif/sigin
192  WRITE(lout,1006)sigdif,sigdif-sigdih,sigdih,sigin
193  1006 FORMAT('SDIFF: SIGDIF,SIGDIL,SIGDIH,SIGIN',4f10.5)
194  1060 FORMAT('SDIFF: KPROJ,KTARG,ECM,SIGDIF/SIGIN',2i3,2f10.5)
195  ENDIF
196  iflagd = 0
197  r = rndm(v)
198  IF ((r.LE.(sigdif/sigin)).OR.(isingd.GE.2)) THEN
199 C NCDIFF = NCDIFF+1
200  2000 CONTINUE
201 C IFLAGD = 1
202 C j.r.11/98 drop the following two lines
203 C GAMV=GAM
204 C BGLV=BGL
205  r = rndm(v)
206 *---------------------- S.Roesler 5/26/93
207  IF (((r.LT.(sigdih/sigdif)).OR.(isingd.EQ.5).OR.
208  & (isingd.EQ.6)).AND.(isingd.LE.6)) THEN
209 *
210 *--------------------- call high mass single diffractive event
211 *
212  IF(isingd.GE.1)THEN
213  CALL vahmsd(ihkk,ecm,kproj,ktarg,irej)
214  iflagd=1
215  ENDIF
216 *
217  ELSE IF((r.GE.(sigdih/sigdif)).OR.(isingd.EQ.7).OR.
218  & (isingd.EQ.8)) THEN
219 *
220 *
221 *--------------------- call low mass single diffractive event
222 *
223 C----------------------------------------------------------------
224 C
225 C j.r. test 2/94
226 C
227 C----------------------------------------------------------------
228  rrrr=rndm(vv)
229  IF(isingd.GE.3.AND.isingd.LE.6)rrrr=1.d0
230  IF((rrrr.LE.0.33d0).AND.
231  * ((kproj.EQ.15).OR.(kproj.EQ.16).OR.
232  * (kproj.EQ.1).OR.(kproj.EQ.8).OR.
233  * (kproj.EQ.13).OR.(kproj.EQ.14).OR.(kproj.EQ.23).OR.
234  * (kproj.EQ.24).OR.(kproj.EQ.25)))THEN
235  CALL valmdd(ihkk,ecm,kproj,ktarg,irej)
236  iflagd=1
237  ELSE
238 C----------------------------------------------------------------
239  IF(isingd.GE.1)THEN
240  CALL valmsd(ihkk,ecm,kproj,ktarg,irej)
241  iflagd=1
242  ENDIF
243  ENDIF
244 *
245  ENDIF
246  IF (irej.GT.0) THEN
247  ncirej = ncirej + 1
248  IF (mod(ncirej,1000).EQ.0) THEN
249  WRITE(lout,1007) ncirej
250  1007 FORMAT('SDIFF: REJECTION, NCIREJ = ',i8)
251  ENDIF
252  goto 2000
253  ENDIF
254 *
255  IF(iflagd.EQ.1)THEN
256  CALL hadrdi(naux,kproj,ktarg,nhkkh1)
257 *
258  iihkk = nhkk -naux
259  DO 11 j=1,naux
260  CALL dtrans(pxc(j),pyc(j),pzc(j),
261  & cod,sid,cof,sif,pxx,pyy,pll)
262  iihkk = iihkk + 1
263  IF (ipev.GE.6) THEN
264  WRITE(lout,1008) iihkk,pxx,pyy,pll
265  1008 FORMAT('SDIFF: NHKK,PXX,PYY,PLL',i4,3f10.5)
266  ENDIF
267  phkk(1,iihkk) = pxx
268  phkk(2,iihkk) = pyy
269  phkk(3,iihkk) = pll
270 C IF (CMHIS.EQ.0.0D0) THEN
271 C CALL DALTRA(GAM,BGX,BGY,BGL,PXX,PYY,PLL,HEC(J),
272 C & PLAB,PHKK(1,IIHKK),PHKK(2,IIHKK),
273 C & PHKK(3,IIHKK),PHKK(4,IIHKK))
274 C Test j.r. 11/98
275 C---------------------------------------------------------------
276 C Back to lab frame
277 C DO 1277 III=NHKKH1,NHKK
278  iii=iihkk
279  IF(isthkk(iii).EQ.1)THEN
280  CALL daltra(gam,bgx,bgy,bgl,phkk(1,iii),phkk(2,iii),
281  & phkk(3,iii),phkk(4,iii),
282  & ppla,pplax,pplay,pplal,epla)
283 C And now into the real cms frame
284  epcms=gamv*epla-bglv*pplal
285  pplcms=gamv*pplal-bglv*epla
286  phkk(1,iii)=pplax
287  phkk(2,iii)=pplay
288  phkk(3,iii)=pplcms
289  phkk(4,iii)=epcms
290  IF (ipev.GE.6) THEN
291  WRITE(lout,1903) pplcms,epcms
292  1903 FORMAT('SDIFF TEST cms: PPLCMS,EPCMS',2e15.5)
293  ENDIF
294  ENDIF
295 C1277 CONTINUE
296 C---------------------------------------------------------------
297  IF (ipev.GE.6) THEN
298  WRITE(lout,1009) iihkk,
299  & phkk(1,iihkk),phkk(2,iihkk),phkk(3,iihkk),
300  & phkk(4,iihkk)
301  1009 FORMAT('SDIFF: NHKK,PHKK(1..4)',i4,4e15.5)
302  ENDIF
303 C ENDIF
304  11 CONTINUE
305  ENDIF
306  ENDIF
307  ENDIF
308  10 CONTINUE
309  IF (ktarg.EQ.0) WRITE(lout,*)'SDIFF: NO INTERACTION'
310  9999 CONTINUE
311  RETURN
312  END
313 *
314 *===vahmsd===============================================================*
315 *
316  SUBROUTINE vahmsd(ITAPOI,ECM,KPROJ,KTARG,IREJ)
317 
318 **************************************************************************
319 * Version November 1993 by Stefan Roesler *
320 * University of Leipzig *
321 * This subroutine selects x-values, flavors and 4-momenta of partons *
322 * in high-mass single diffractive chains. *
323 **************************************************************************
324 
325  IMPLICIT DOUBLE PRECISION (a-h,o-z)
326  SAVE
327  parameter(lout=6,llook=9)
328  parameter(nmxhkk= 89998)
329  CHARACTER*8 aname
330  COMMON /hkkevt/ nhkk,nevhkk, isthkk(nmxhkk), idhkk(nmxhkk),
331  & jmohkk(2,nmxhkk),jdahkk(2,nmxhkk),phkk(5,nmxhkk),
332  & vhkk(4,nmxhkk), whkk(4,nmxhkk)
333  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
334  COMMON /diffra/ isingd,idiftp,ioudif,iflagd
335  COMMON /abrdif/ xdq1,xdq2,xddq1,xddq2,
336  & ikvq1,ikvq2,ikd1q1,ikd2q1,ikd1q2,ikd2q2,
337  & idiffp,idifap,
338  & amdch1,amdch2,amdch3,gamdc1,gamdc2,gamdc3,
339  & pgxvc1,pgyvc1,pgzvc1,pgxvc2,pgyvc2,pgzvc2,
340  & pgxvc3,pgyvc3,pgzvc3,ndch1,ndch2,ndch3,
341  & ikdch1,ikdch2,ikdch3,
342  & pdq1(4),pdq2(4),pdd1(4),pdd2(4),pdfq1(4)
343  COMMON /dpar/ aname(210),am(210),ga(210),tau(210),ich(210),
344  & ibar(210),k1(210),k2(210)
345  COMMON /trafop/ gamp,bgamp,betp
346  COMMON /enerin/ eproj,etarg
347  COMMON /sdflag/ isd
348  COMMON /xdidid/xdidi
349 *
350  dimension mquark(3,30),ihkkq(-6:6),ihkkqq(-3:3,-3:3),
351  & idx(-4:4)
352  DATA idx /-4,-3,-1,-2,0,2,1,3,4/
353  DATA ihkkq /-6,-5,-4,-3,-1,-2,0,2,1,3,4,5,6/
354  DATA ihkkqq/-3301,-3103,-3203,0, 0,0,0,
355  & -3103,-1103,-2103,0, 0, 0, 0,
356  & -3203,-2103,-2203,0, 0, 0, 0,
357  & 0, 0, 0,0, 0, 0, 0,
358  & 0, 0, 0,0,2203,2103,3202,
359  & 0, 0, 0,0,2103,1103,3103,
360  & 0, 0, 0,0,3203,3103,3303/
361 *
362 *----------------------------------- quark content of hadrons:
363 * 1, 2, 3, 4 - u, d, s, c
364 * -1,-2,-3,-4 - au,ad,as,ac
365 *
366  DATA mquark/
367  & 1,1,2, -1,-1,-2, 0,0,0, 0,0,0, 0,0,0,
368  & 0,0,0, 0,0,0, 1,2,2, -1,-2,-2, 0,0,0,
369  & 0,0,0, 0,0,0, 1,-2,0, 2,-1,0, 1,-3,0,
370  & 3,-1,0, 1,2,3, -1,-2,-3, 0,0,0, 2,2,3,
371  & 1,1,3, 1,2,3, 1,-1,0, 2,-3,0, 3,-2,0,
372  & 2,-2,0, 3,-3,0, 0,0,0, 0,0,0, 0,0,0/
373  DATA unon/2.0/
374  DATA ncrej, ncxdi, ncxp, ncxt /0, 0, 0, 0/
375 *
376  isd = 1
377  irej = 0
378  iirej = 0
379  eproj = (am(kproj)**2-am(ktarg)**2+ecm**2)/(2.0d0*ecm)
380  etarg = (am(ktarg)**2-am(kproj)**2+ecm**2)/(2.0d0*ecm)
381  IF(ipev.GE.2) WRITE(lout,1001) eproj,etarg
382  1001 FORMAT('VAHMSD: EPROJ,ETARG ',2f10.5)
383  ibproj = ibar(kproj)
384  ibtarg = ibar(ktarg)
385  IF (ibtarg.LE.0) THEN
386  WRITE(lout,1002) ibtarg
387  1002 FORMAT('VAHMSD: NO HMSD FOR TARGET WITH BARYON-CHARGE',i4)
388  iirej=1
389  goto 9999
390  ENDIF
391  iqp1 = mquark(1,kproj)
392  iqp2 = mquark(2,kproj)
393  iqp3 = mquark(3,kproj)
394  iqt1 = mquark(1,ktarg)
395  iqt2 = mquark(2,ktarg)
396  iqt3 = mquark(3,ktarg)
397  IF(ipev.GE.2) WRITE(lout,1003)
398  & ibproj,ibtarg,iqp1,iqp2,iqp3,iqt1,iqt2,iqt3
399  1003 FORMAT('VAHMSD: IBPROJ,IBTARG,IQP1,IQP2,IQP3,IQT1,IQT2,IQT3 ',8i4)
400 *
401  IF (ibproj.NE.0) THEN
402 *
403 *-------------------- q-qq (aq-aqaq) - flavors of projectile (baryon)
404 *
405  isam = 1.0d0+2.999d0*rndm(v)
406  goto(10,11,12) isam
407  10 CONTINUE
408  iqp = iqp1
409  idiqp1 = iqp2
410  idiqp2 = iqp3
411  goto 13
412  11 CONTINUE
413  iqp = iqp2
414  idiqp1 = iqp1
415  idiqp2 = iqp3
416  goto 13
417  12 CONTINUE
418  iqp = iqp3
419  idiqp1 = iqp1
420  idiqp2 = iqp2
421  13 CONTINUE
422 *
423  ELSE IF (ibproj.EQ.0) THEN
424 *
425 *-------------------- q-aq - flavors of projectile (meson)
426 *
427  isam = 1.0d0+1.999d0*rndm(v)
428  goto(14,15) isam
429  14 CONTINUE
430  iqp = iqp1
431  idiqp1 = iqp2
432  goto 16
433  15 CONTINUE
434  iqp = iqp2
435  idiqp1 = iqp1
436  16 CONTINUE
437  idiqp2 = 0
438 *
439  ENDIF
440 *
441 *-------------------- q-qq - flavors of target (baryon)
442 *
443  isam = 1.0d0+2.999d0*rndm(v)
444  goto(17,18,19) isam
445  17 CONTINUE
446  iqt = iqt1
447  idiqt1 = iqt2
448  idiqt2 = iqt3
449  goto 20
450  18 CONTINUE
451  iqt = iqt2
452  idiqt1 = iqt1
453  idiqt2 = iqt3
454  goto 20
455  19 CONTINUE
456  iqt = iqt3
457  idiqt1 = iqt1
458  idiqt2 = iqt2
459  20 CONTINUE
460 *
461  ikvq1 = iqp
462  ikd1q1 = idiqp1
463  ikd2q1 = idiqp2
464 *
465  ikvq2 = iqt
466  ikd1q2 = idiqt1
467  ikd2q2 = idiqt2
468 *
469  IF (ipev.GE.2) WRITE(lout,1004)
470  & ikvq1,ikd1q1,ikd2q1,ikvq2,ikd1q2,ikd2q2
471  1004 FORMAT('VAHMSD: IKVQ1,IKD1Q1,IKD2Q1,IKVQ2,IKD1Q2,IKD2Q2 ',6i4)
472 *
473 *-------------------- q-aq - flavors of diffractive parton
474 *
475  idiffp = 1.0d0+2.3d0*rndm(v)
476  idifap = -idiffp
477  IF (ipev.GE.2) WRITE(lout,1005) idiffp,idifap
478  1005 FORMAT('VAHMSD: IDIFFP,IDIFAP ',2i4)
479 *
480 *-------------------- IDIFTP = 1 target (backward) hadron excited
481 * IDIFTP = 2 projectile (forward) hadron excited
482 *
483  idiftp = 1.0d0+1.999d0*rndm(v)
484 *----------------------- S.Roesler 5/26/93
485  IF ((isingd.EQ.3).OR.(isingd.EQ.5)) idiftp = 1
486  IF ((isingd.EQ.4).OR.(isingd.EQ.6)) idiftp = 2
487 *
488  IF (ipev.GE.2) WRITE(lout,1006) idiftp
489  1006 FORMAT('VAHMSD: IDIFTP ',i4)
490  IF ((idiftp.NE.1).AND.(idiftp.NE.2)) THEN
491  IF (ipev.GE.2) WRITE(lout,'(A19)') 'VAHMSD-ERROR: IDIFTP'
492  goto 9999
493  ENDIF
494 *
495 *-------------------- momentum fractions of quarks and diquarks
496 *
497  xxmax = 0.8d0
498  30 CONTINUE
499  xp = dbetar(0.5d0,unon)
500  IF (xp.GE.xxmax) THEN
501  ncxp = ncxp+1
502  IF(mod(ncxp,500).EQ.0) WRITE(lout,1007) ncxp
503  1007 FORMAT('VAHMSD: INEFFICIENT XP-SELECTION, NCXP=',i8)
504  goto 30
505  ENDIF
506  xxp = 1.0d0-xp
507  31 CONTINUE
508  xt = dbetar(0.5d0,unon)
509  IF (xt.GE.xxmax) THEN
510  ncxt = ncxt+1
511  IF(mod(ncxt,2500).EQ.0) WRITE(lout,1008) ncxt
512  1008 FORMAT('VAHMSD: INEFFICIENT XT-SELECTION, NCXT=',i8)
513  goto 31
514  ENDIF
515  xxt = 1.0d0-xt
516  IF (ipev.GE.2) WRITE(lout,1010) xp,xxp,xt,xxt
517  1010 FORMAT('VAHMSD: XP,XXP,XT,XXT ',4d10.5)
518 *-------------------- x - values of diffractive partons
519 *
520  ncxdi = 0
521  32 CONTINUE
522  r = rndm(v)
523  IF (idiftp.EQ.1) THEN
524  xdimin = (3.0d0+400.0d0*(r**2))/(4.0d0*(etarg**2)*xxt)
525  IF (ecm.LE.300.0d0) THEN
526  rr = (1.0d0-exp(-((ecm/140.0d0)**4)))
527  xdimin = (3.0d0+400.0d0*(r**2)*rr)/(4.0d0*(etarg**2)*xxt)
528  ENDIF
529  ELSE IF (idiftp.EQ.2) THEN
530  xdimin = (3.0d0+400.0d0*(r**2))/(4.0d0*(eproj**2)*xxp)
531  IF (ecm.LE.300.0d0) THEN
532  rr = (1.0d0-exp(-((ecm/140.0d0)**4)))
533  xdimin = (3.0d0+400.0d0*(r**2)*rr)/(4.0d0*(eproj**2)*xxp)
534  ENDIF
535  ENDIF
536 C-----------------------------------------------------------------
537 C original version
538 C-----------------------------------------------------------------
539 C XDIMAX = 0.05D0
540 C IF (ECM.LE.1000.0D0) THEN
541 C XDIMAX = 0.05D0*(1.0D0+EXP(-((ECM/420.0D0)**2)))
542 C IF (IBPROJ.EQ.0) XDIMAX = 0.05D0*
543 CC---------------- change mass-cuts
544 CC & (1.0D0+4.0D0*EXP(-((ECM/420.0D0)**2)))
545 C & (1.0D0+2.0D0*EXP(-((ECM/420.0D0)**2)))
546 C ENDIF
547 C-----------------------------------------------------------------
548 C-----------------------------------------------------------------
549 C version j.r.28.1.94
550 C extend diffraction beyond limits of single diffractive experiments
551 C-----------------------------------------------------------------
552  xdimaa = 0.15d0
553  xdimax = xdimaa
554  IF (ecm.LE.10000.0d0) THEN
555  xdimax = xdimaa*(1.0d0+exp(-((ecm/420.0d0)**2)))
556  IF (ibproj.EQ.0) xdimax = xdimaa*
557 C----------------- change mass-cuts
558 C & (1.0D0+4.0D0*EXP(-((ECM/420.0D0)**2)))
559  & (1.0d0+2.0d0*exp(-((ecm/420.0d0)**2)))
560  ENDIF
561 C-----------------------------------------------------------------
562  40 CONTINUE
563  IF (xdimin.GE.xdimax) THEN
564  ncxdi = ncxdi+1
565  IF (ncxdi.EQ.200) goto 9999
566  goto 32
567  ENDIF
568  xditot = sampey(xdimin,xdimax)
569  r = rndm(v)**6
570  dx = xditot-xdimin
571  xdi = dx*r
572 C IF (XDI.LT.XMINQ1) THEN
573 C NCXDI = NCXDI+1
574 C IF (MOD(NCXDI,2000).EQ.0) WRITE(LOUT,1011) NCXDI
575 C1011 FORMAT('VAHMSD: INEFFICIENT XDI-SELECTION, NCXDI=',I8)
576 C GOTO 40
577 C ENDIF
578  xxdi = xdimin+(1.0d0-r)*dx
579  IF (idiftp.EQ.1) THEN
580  xxp = xxp-xdi-xxdi
581  IF (xxp.LT.8.0d-2) goto 9999
582  ELSE IF (idiftp.EQ.2) THEN
583  xxt = xxt-xdi-xxdi
584  IF (xxt.LT.8.0d-2) goto 9999
585  ENDIF
586  IF (ipev.GE.2) WRITE(lout,1012) xp,xxp,xt,xxt,xdi,xxdi
587  1012 FORMAT('VAHMSD: XP,XXP,XT,XXT,XDI,XXDI ',6f10.5)
588  xdidi=xdi+xxdi
589  amdidi=sqrt(xdidi*ecm**2)
590  IF(ipev.GE.2)WRITE(lout,*)'HM AMDIDI,XDIDI ',amdidi,xdidi
591 *
592 *-------------------- kinematical parameters of three chains in CMS
593 *
594  IF ((ibproj.EQ.-1).AND.(idiftp.EQ.1)) THEN
595 *
596 *-------------------- target: baryon, projectile: antibaryon
597 * (excited)
598 *
599  xdifap = xdi
600  xdiffp = xxdi
601  CALL diffch(xdifap, idifap, xt, ikvq2, 99,
602  & xdiffp, idiffp, xxt, xxp, etarg,
603  & amdch1, ech1, pch1, gamdc1, pgvc1,
604  & ndch1, ikdch1, eproj, nuno, iirej, 1)
605  IF (iirej.EQ.1) goto 9999
606  CALL diffch(xdiffp, idiffp, xxt, ikd1q2, ikd2q2,
607  & dum, idum, dum, xxp, etarg,
608  & amdch2, ech2, pch2, gamdc2, pgvc2,
609  & ndch2, ikdch2, eproj, nuno, iirej, 2)
610  IF (iirej.EQ.1) goto 9999
611  CALL diffch( xp, ikvq1, xxp, ikd1q1, ikd2q1,
612  & dum, idum, dum, dum, eproj,
613  & amdch3, ech3, pch3, gamdc3, pgvc3,
614  & ndch3, idum, dum, nuno, iirej, 3)
615  IF (iirej.EQ.1) goto 9999
616 *
617 *
618  ELSE IF ((ibproj.EQ.-1).AND.(idiftp.EQ.2)) THEN
619 *
620 *-------------------- target: baryon, projectile: antibaryon
621 * (excited)
622 *
623  xdifap = xxdi
624  xdiffp = xdi
625  CALL diffch(xdiffp, idiffp, xp, ikvq1, 99,
626  & xdifap, idifap, xxp, xxt, eproj,
627  & amdch1, ech1, pch1, gamdc1, pgvc1,
628  & ndch1, ikdch1, etarg, nuno, iirej, 1)
629  IF (iirej.EQ.1) goto 9999
630  CALL diffch(xdifap, idifap, xxp, ikd1q1, ikd2q1,
631  & dum, idum, dum, xxt, eproj,
632  & amdch2, ech2, pch2, gamdc2, pgvc2,
633  & ndch2, ikdch2, etarg, nuno, iirej, 2)
634  IF (iirej.EQ.1) goto 9999
635  CALL diffch( xt, ikvq2, xxt, ikd1q2, ikd2q2,
636  & dum, idum, dum, dum, etarg,
637  & amdch3, ech3, pch3, gamdc3, pgvc3,
638  & ndch3, idum, dum, nuno, iirej, 3)
639  IF (iirej.EQ.1) goto 9999
640 *
641 *
642  ELSE IF ((ibproj.EQ.0).AND.(idiftp.EQ.1)) THEN
643 *
644 *-------------------- target: baryon, projectile: meson
645 * (excited)
646 *
647  xdifap = xdi
648  xdiffp = xxdi
649  CALL diffch(xdifap, idifap, xt, ikvq2, 99,
650  & xdiffp, idiffp, xxt, xxp, etarg,
651  & amdch1, ech1, pch1, gamdc1, pgvc1,
652  & ndch1, ikdch1, eproj, nuno, iirej, 1)
653  IF (iirej.EQ.1) goto 9999
654  CALL diffch(xdiffp, idiffp, xxt, ikd1q2, ikd2q2,
655  & dum, idum, dum, xxp, etarg,
656  & amdch2, ech2, pch2, gamdc2, pgvc2,
657  & ndch2, ikdch2, eproj, nuno, iirej, 2)
658  IF (iirej.EQ.1) goto 9999
659  CALL diffch( xp, ikvq1, xxp, ikd1q1, 99,
660  & dum, idum, dum, dum, eproj,
661  & amdch3, ech3, pch3, gamdc3, pgvc3,
662  & ndch3, idum, dum, nuno, iirej, 3)
663  IF (iirej.EQ.1) goto 9999
664 *
665 *
666  ELSEIF ((ibproj.EQ.0).AND.(idiftp.EQ.2)) THEN
667 *
668 *-------------------- target: baryon, projectile: meson
669 * (excited)
670 *
671  IF (ikd1q1.LT.0) THEN
672  xdifap = xdi
673  xdiffp = xxdi
674  CALL diffch(xdifap, idifap, xp, ikvq1, 99,
675  & xdiffp, idiffp, xxp, xxt, eproj,
676  & amdch1, ech1, pch1, gamdc1, pgvc1,
677  & ndch1, ikdch1, etarg, nuno, iirej, 1)
678  IF (iirej.EQ.1) goto 9999
679  CALL diffch(xdiffp, idiffp, xxp, ikd1q1, 99,
680  & xdifap, idifap, xp, xxt, eproj,
681  & amdch2, ech2, pch2, gamdc2, pgvc2,
682  & ndch2, ikdch2, etarg, nuno, iirej, 1)
683  IF (iirej.EQ.1) goto 9999
684  ELSE
685  xdifap = xxdi
686  xdiffp = xdi
687  CALL diffch(xdiffp, idiffp, xp, ikvq1, 99,
688  & xdifap, idifap, xxp, xxt, eproj,
689  & amdch1, ech1, pch1, gamdc1, pgvc1,
690  & ndch1, ikdch1, etarg, nuno, iirej, 1)
691  IF (iirej.EQ.1) goto 9999
692  CALL diffch(xdifap, idifap, xxp, ikd1q1, 99,
693  & xdiffp, idiffp, xp, xxt, eproj,
694  & amdch2, ech2, pch2, gamdc2, pgvc2,
695  & ndch2, ikdch2, etarg, nuno, iirej, 1)
696  IF (iirej.EQ.1) goto 9999
697  ENDIF
698  CALL diffch( xt, ikvq2, xxt, ikd1q2, ikd2q2,
699  & dum, idum, dum, dum, etarg,
700  & amdch3, ech3, pch3, gamdc3, pgvc3,
701  & ndch3, idum, dum, nuno, iirej, 3)
702  IF (iirej.EQ.1) goto 9999
703 *
704 *
705  ELSEIF ((ibproj.EQ.1).AND.(idiftp.EQ.1)) THEN
706 *
707 *-------------------- target: baryon, projectile: baryon
708 * (excited)
709 *
710  xdifap = xdi
711  xdiffp = xxdi
712  CALL diffch(xdifap, idifap, xt, ikvq2, 99,
713  & xdiffp, idiffp, xxt, xxp, etarg,
714  & amdch1, ech1, pch1, gamdc1, pgvc1,
715  & ndch1, ikdch1, eproj, nuno, iirej, 1)
716  IF (iirej.EQ.1) goto 9999
717  CALL diffch(xdiffp, idiffp, xxt, ikd1q2, ikd2q2,
718  & dum, idum, dum, xxp, etarg,
719  & amdch2, ech2, pch2, gamdc2, pgvc2,
720  & ndch2, ikdch2, eproj, nuno, iirej, 2)
721  IF (iirej.EQ.1) goto 9999
722  CALL diffch( xp, ikvq1, xxp, ikd1q1, ikd2q1,
723  & dum, idum, dum, dum, eproj,
724  & amdch3, ech3, pch3, gamdc3, pgvc3,
725  & ndch3, idum, dum, nuno, iirej, 3)
726  IF (iirej.EQ.1) goto 9999
727 *
728 *
729  ELSE IF ((ibproj.EQ.1).AND.(idiftp.EQ.2)) THEN
730 *
731 *-------------------- target: baryon, projectile: baryon
732 * (excited)
733 *
734  xdifap = xdi
735  xdiffp = xxdi
736  CALL diffch(xdifap, idifap, xp, ikvq1, 99,
737  & xdiffp, idiffp, xxp, xxt, eproj,
738  & amdch1, ech1, pch1, gamdc1, pgvc1,
739  & ndch1, ikdch1, etarg, nuno, iirej, 1)
740  IF (iirej.EQ.1) goto 9999
741  CALL diffch(xdiffp, idiffp, xxp, ikd1q1, ikd2q1,
742  & dum, idum, dum, xxt, eproj,
743  & amdch2, ech2, pch2, gamdc2, pgvc2,
744  & ndch2, ikdch2, etarg, nuno, iirej, 2)
745  IF (iirej.EQ.1) goto 9999
746  CALL diffch( xt, ikvq2, xxt, ikd1q2, ikd2q2,
747  & dum, idum, dum, dum, etarg,
748  & amdch3, ech3, pch3, gamdc3, pgvc3,
749  & ndch3, idum, dum, nuno, iirej, 3)
750  IF (iirej.EQ.1) goto 9999
751 *
752  ENDIF
753 *
754 *-------------------- store results in common block /ABRDIF/
755 *
756  xdq1 = xp
757  xdq2 = xt
758  xddq1 = xxp
759  xddq2 = xxt
760 *
761  IF (ipev.GE.2) THEN
762  WRITE(lout,1013) amdch1,ech1,pch1,gamdc1,pgvc1,ikdch1,ndch1
763  1013 FORMAT('VAHMSD: AMDCH1,ECH1,PCH1,GAMDC1,PGVC1,IKDCH1,NDCH1 ',
764  & 5f10.5,2i4)
765  WRITE(lout,1014) amdch2,ech2,pch2,gamdc2,pgvc2,ikdch2,ndch2
766  1014 FORMAT('VAHMSD: AMDCH2,ECH2,PCH2,GAMDC2,PGVC2,IKDCH2,NDCH2 ',
767  & 5f10.5,2i4)
768  WRITE(lout,1015) amdch3,ech3,pch3,gamdc3,pgvc3,ndch3
769  1015 FORMAT('VAHMSD: AMDCH3,ECH3,PCH3,GAMDC3,PGVC3,NDCH3 ',
770  & 5f10.5,i4)
771  ENDIF
772 *
773 *-------------------- select transverse momenta
774 *
775  IF (idiftp.EQ.1) THEN
776  CALL diffpt( ech1, pch1, xdifap, xt,
777  & ech2, pch2, xdiffp, xxt,
778  & ech3, pch3, eproj, kproj,
779  & etarg, eproj, iirej)
780  ELSE IF (idiftp.EQ.2) THEN
781  IF ((ibproj.LT.0).OR.(ikvq1.LT.0)) THEN
782  CALL diffpt( ech1, pch1, xdiffp, xp,
783  & ech2, pch2, xdifap, xxp,
784  & ech3, pch3, etarg, ktarg,
785  & eproj, etarg, iirej)
786  ELSE
787  CALL diffpt( ech1, pch1, xdifap, xp,
788  & ech2, pch2, xdiffp, xxp,
789  & ech3, pch3, etarg, ktarg,
790  & eproj, etarg, iirej)
791  ENDIF
792  ENDIF
793  IF (iirej.EQ.1) goto 9999
794 *
795 *-------------------- store results in common-block /HKKEVT/
796 *
797 * partonen of projectile
798 *
799  nhkk = nhkk+1
800  isthkk(nhkk) = 21
801  idhkk(nhkk) = ihkkq(idx(ikvq1))
802  jmohkk(1,nhkk) = 1
803  jmohkk(2,nhkk) = 0
804  jdahkk(1,nhkk) = 0
805  jdahkk(2,nhkk) = 0
806  phkk(1,nhkk) = 0.0d0
807  phkk(2,nhkk) = 0.0d0
808  phkk(3,nhkk) = xdq1
809  phkk(4,nhkk) = xdq1
810  phkk(5,nhkk) = 0.0d0
811  vhkk(1,nhkk) = vhkk(1,1)
812  vhkk(2,nhkk) = vhkk(2,1)
813  vhkk(3,nhkk) = vhkk(3,1)
814  vhkk(4,nhkk) = vhkk(4,1)
815 *
816  nhkk = nhkk+1
817  isthkk(nhkk) = 21
818  idhkk(nhkk) = ihkkqq(idx(ikd1q1),idx(ikd2q1))
819  jmohkk(1,nhkk) = 1
820  jmohkk(2,nhkk) = 0
821  jdahkk(1,nhkk) = 0
822  jdahkk(2,nhkk) = 0
823  phkk(1,nhkk) = 0.0d0
824  phkk(2,nhkk) = 0.0d0
825  phkk(3,nhkk) = xddq1
826  phkk(4,nhkk) = xddq1
827  phkk(5,nhkk) = 0.0d0
828  vhkk(1,nhkk) = vhkk(1,1)
829  vhkk(2,nhkk) = vhkk(2,1)
830  vhkk(3,nhkk) = vhkk(3,1)
831  vhkk(4,nhkk) = vhkk(4,1)
832 *
833 * partonen of target
834 *
835  nhkk = nhkk+1
836  isthkk(nhkk) = 22
837  idhkk(nhkk) = ihkkq(idx(ikvq2))
838  jmohkk(1,nhkk) = itapoi
839  jmohkk(2,nhkk) = 0
840  jdahkk(1,nhkk) = 0
841  jdahkk(2,nhkk) = 0
842  phkk(1,nhkk) = 0.0d0
843  phkk(2,nhkk) = 0.0d0
844  phkk(3,nhkk) = xdq2
845  phkk(4,nhkk) = xdq2
846  phkk(5,nhkk) = 0.0d0
847  vhkk(1,nhkk) = vhkk(1,itapoi)
848  vhkk(2,nhkk) = vhkk(2,itapoi)
849  vhkk(3,nhkk) = vhkk(3,itapoi)
850  vhkk(4,nhkk) = vhkk(4,itapoi)
851 *
852  nhkk = nhkk+1
853  isthkk(nhkk) = 22
854  idhkk(nhkk) = ihkkqq(idx(ikd1q2),idx(ikd2q2))
855  jmohkk(1,nhkk) = itapoi
856  jmohkk(2,nhkk) = 0
857  jdahkk(1,nhkk) = 0
858  jdahkk(2,nhkk) = 0
859  phkk(1,nhkk) = 0.0d0
860  phkk(2,nhkk) = 0.0d0
861  phkk(3,nhkk) = xddq2
862  phkk(4,nhkk) = xddq2
863  phkk(5,nhkk) = 0.0d0
864  vhkk(1,nhkk) = vhkk(1,itapoi)
865  vhkk(2,nhkk) = vhkk(2,itapoi)
866  vhkk(3,nhkk) = vhkk(3,itapoi)
867  vhkk(4,nhkk) = vhkk(4,itapoi)
868 *
869 * sea q-aq pair
870 *
871  nhkk = nhkk+1
872  isthkk(nhkk) = 30+idiftp
873  idhkk(nhkk) = ihkkq(idx(idiffp))
874  IF (idiftp.EQ.1) jmohkk(1,nhkk) = 1
875  IF (idiftp.EQ.2) jmohkk(1,nhkk) = itapoi
876  jmohkk(2,nhkk) = 0
877  jdahkk(1,nhkk) = 0
878  jdahkk(2,nhkk) = 0
879  phkk(1,nhkk) = 0.0d0
880  phkk(2,nhkk) = 0.0d0
881  phkk(3,nhkk) = xxdi
882  phkk(4,nhkk) = xxdi
883  phkk(5,nhkk) = 0.0d0
884  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))
885  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))
886  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
887  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
888 *
889  nhkk = nhkk+1
890  isthkk(nhkk) = 30+idiftp
891  idhkk(nhkk) = ihkkq(idx(idifap))
892  IF (idiftp.EQ.1) jmohkk(1,nhkk) = 1
893  IF (idiftp.EQ.2) jmohkk(1,nhkk) = itapoi
894  jmohkk(2,nhkk) = 0
895  jdahkk(1,nhkk) = 0
896  jdahkk(2,nhkk) = 0
897  phkk(1,nhkk) = 0.0d0
898  phkk(2,nhkk) = 0.0d0
899  phkk(3,nhkk) = xdi
900  phkk(4,nhkk) = xdi
901  phkk(5,nhkk) = 0.0d0
902  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))
903  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))
904  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
905  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
906 *
907 * ends of chain 1
908 *
909  nhkk = nhkk+1
910  isthkk(nhkk) = 123-idiftp
911  ind = nhkk-2-2*idiftp
912  idhkk(nhkk) = idhkk(ind)
913  jmohkk(1,nhkk) = ind
914  jmohkk(2,nhkk) = jmohkk(1,ind)
915  jdahkk(1,nhkk) = nhkk+2
916  jdahkk(2,nhkk) = nhkk+2
917  phkk(1,nhkk) = pdq1(1)
918  phkk(2,nhkk) = pdq1(2)
919  phkk(3,nhkk) = pdq1(3)
920  phkk(4,nhkk) = pdq1(4)
921  phkk(5,nhkk) = 0.0d0
922 C Add position of parton in hadron
923  CALL qinnuc(xxpp,yypp)
924  vhkk(1,nhkk) = vhkk(1,ind)+xxpp
925  vhkk(2,nhkk) = vhkk(2,ind)+yypp
926  vhkk(3,nhkk) = vhkk(3,ind)
927  vhkk(4,nhkk) = vhkk(4,ind)
928 *
929  nhkk = nhkk+1
930  isthkk(nhkk) = 130+idiftp
931  IF (idhkk(nhkk-1).GT.0) jmohkk(1,nhkk) = nhkk-2
932  IF (idhkk(nhkk-1).LE.0) jmohkk(1,nhkk) = nhkk-3
933  jmohkk(2,nhkk) = jmohkk(1,nhkk-3)
934  jdahkk(1,nhkk) = nhkk+1
935  jdahkk(2,nhkk) = nhkk+1
936  idhkk(nhkk) = idhkk(jmohkk(1,nhkk))
937  phkk(1,nhkk) = pdd2(1)
938  phkk(2,nhkk) = pdd2(2)
939  phkk(3,nhkk) = pdd2(3)
940  phkk(4,nhkk) = pdd2(4)
941  phkk(5,nhkk) = 0.0d0
942 C Add position of parton in hadron
943  CALL qinnuc(xxpp,yypp)
944  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))+xxpp
945  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))+yypp
946  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
947  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
948 *
949 * chain 1
950 *
951  nhkk = nhkk+1
952  isthkk(nhkk) = 199
953  idhkk(nhkk) = 88888
954  jmohkk(1,nhkk) = nhkk-2
955  jmohkk(2,nhkk) = nhkk-1
956  jdahkk(1,nhkk) = 0
957  jdahkk(2,nhkk) = 0
958  phkk(5,nhkk) = amdch1
959  vhkk(1,nhkk) = vhkk(1,nhkk-1)
960  vhkk(2,nhkk) = vhkk(2,nhkk-1)
961  vhkk(3,nhkk) = vhkk(3,nhkk-1)
962  IF ((betp.NE.0.0d0).AND.(bgamp.NE.0.0d0))
963  &vhkk(4,nhkk) = vhkk(3,nhkk)/betp-vhkk(3,nhkk-2)/bgamp
964 *
965 * ends of chain 2
966 *
967  nhkk = nhkk+1
968  isthkk(nhkk) = 123-idiftp
969  ind = nhkk-4-2*idiftp
970  idhkk(nhkk) = idhkk(ind)
971  jmohkk(1,nhkk) = ind
972  jmohkk(2,nhkk) = jmohkk(1,ind)
973  jdahkk(1,nhkk) = nhkk+2
974  jdahkk(2,nhkk) = nhkk+2
975  phkk(1,nhkk) = pdd1(1)
976  phkk(2,nhkk) = pdd1(2)
977  phkk(3,nhkk) = pdd1(3)
978  phkk(4,nhkk) = pdd1(4)
979  phkk(5,nhkk) = 0.0d0
980 C Add position of parton in hadron
981  CALL qinnuc(xxpp,yypp)
982  vhkk(1,nhkk) = vhkk(1,ind)+xxpp
983  vhkk(2,nhkk) = vhkk(2,ind)+yypp
984  vhkk(3,nhkk) = vhkk(3,ind)
985  vhkk(4,nhkk) = vhkk(4,ind)
986 *
987  nhkk = nhkk+1
988  isthkk(nhkk) = 130+idiftp
989  jmohkk(1,nhkk) = 2*nhkk-11-jmohkk(1,nhkk-3)
990  jmohkk(2,nhkk) = jmohkk(1,nhkk-6)
991  jdahkk(1,nhkk) = nhkk+1
992  jdahkk(2,nhkk) = nhkk+1
993  idhkk(nhkk) = idhkk(jmohkk(1,nhkk))
994  phkk(1,nhkk) = pdq2(1)
995  phkk(2,nhkk) = pdq2(2)
996  phkk(3,nhkk) = pdq2(3)
997  phkk(4,nhkk) = pdq2(4)
998  phkk(5,nhkk) = 0.0d0
999 C Add position of parton in hadron
1000  CALL qinnuc(xxpp,yypp)
1001  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk-6))+xxpp
1002  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk-6))+yypp
1003  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk-6))
1004  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk-6))
1005 *
1006 * chain 2
1007 *
1008  nhkk = nhkk+1
1009  isthkk(nhkk) = 199
1010  idhkk(nhkk) = 88888
1011  jmohkk(1,nhkk) = nhkk-2
1012  jmohkk(2,nhkk) = nhkk-1
1013  jdahkk(1,nhkk) = 0
1014  jdahkk(2,nhkk) = 0
1015  phkk(5,nhkk) = amdch2
1016  vhkk(1,nhkk) = vhkk(1,nhkk-1)
1017  vhkk(2,nhkk) = vhkk(2,nhkk-1)
1018  vhkk(3,nhkk) = vhkk(3,nhkk-1)
1019  IF ((betp.NE.0.0d0).AND.(bgamp.NE.0.0d0))
1020  &vhkk(4,nhkk) = vhkk(3,nhkk)/betp-vhkk(3,nhkk-2)/bgamp
1021 *
1022 * diffractive nucleon
1023 *
1024  nhkk = nhkk+1
1025  isthkk(nhkk) = 199
1026  idhkk(nhkk) = 88888
1027  jmohkk(1,nhkk) = nhkk-14+2*idiftp
1028  jmohkk(2,nhkk) = nhkk-13+2*idiftp
1029  jdahkk(1,nhkk) = 0
1030  jdahkk(2,nhkk) = 0
1031  phkk(1,nhkk) = pdfq1(1)
1032  phkk(2,nhkk) = pdfq1(2)
1033  phkk(3,nhkk) = pdfq1(3)
1034  phkk(4,nhkk) = pdfq1(4)
1035  phkk(5,nhkk) = amdch3
1036  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))
1037  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))
1038  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
1039  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
1040 C WRITE(6,*)' Diffr. Nucleon'
1041 C *,PHKK(3,NHKK),PHKK(4,NHKK),PHKK(5,NHKK)
1042 *
1043 *---------- S. Roesler 21-10-93
1044 * check energy-momentum conservation
1045  IF (ipev.GE.2) THEN
1046  px = pdq1(1)+pdq2(1)+pdd1(1)+pdd2(1)+pdfq1(1)
1047  py = pdq1(2)+pdq2(2)+pdd1(2)+pdd2(2)+pdfq1(2)
1048  pz = pdq1(3)+pdq2(3)+pdd1(3)+pdd2(3)+pdfq1(3)
1049  ee = ecm-(pdq1(4)+pdq2(4)+pdd1(4)+pdd2(4)+pdfq1(4))
1050  WRITE(lout,*)'VAHMSD: ENERGY-MOMENTUM-CHECK (PX,PY,PZ,E)'
1051  WRITE(lout,'(5F12.6)')px,py,pz,ee,ecm
1052  ENDIF
1053 *
1054  RETURN
1055 *
1056  9999 CONTINUE
1057  ncrej = ncrej+1
1058  IF (mod(ncrej,2500).EQ.0) WRITE(lout,9900) ncrej
1059  9900 FORMAT('REJECTION IN VAHMSD ',i10)
1060  iirej = 0
1061  irej = 1
1062  RETURN
1063  END
1064 *
1065 *===diffpt===============================================================*
1066 *
1067  SUBROUTINE diffpt( ECH1, PCH1, XDICH1, XCH1,
1068  & ech2, pch2, xdich2, xch2,
1069  & ech3, pch3, ech3i, kch3i,
1070  & ee, ees, iirej)
1071 
1072 **************************************************************************
1073 * Version November 1993 by Stefan Roesler *
1074 * University of Leipzig *
1075 * This subroutine is called by VAHMSD (HMSD), selects the transverse *
1076 * momenta of the chains and the diffractive nucleon. *
1077 * The results are stored in the common-block ABRDIF. *
1078 **************************************************************************
1079 
1080  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1081  SAVE
1082  parameter(lout=6,llook=9)
1083  parameter(nmxhkk= 89998)
1084  CHARACTER*8 aname
1085  COMMON /hkkevt/ nhkk,nevhkk, isthkk(nmxhkk), idhkk(nmxhkk),
1086  & jmohkk(2,nmxhkk),jdahkk(2,nmxhkk),phkk(5,nmxhkk),
1087  & vhkk(4,nmxhkk), whkk(4,nmxhkk)
1088  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
1089  COMMON /diffra/ isingd,idiftp,ioudif,iflagd
1090  COMMON /abrdif/ xdq1,xdq2,xddq1,xddq2,
1091  & ikvq1,ikvq2,ikd1q1,ikd2q1,ikd1q2,ikd2q2,
1092  & idiffp,idifap,
1093  & amdch1,amdch2,amdch3,gamdc1,gamdc2,gamdc3,
1094  & pgxvc1,pgyvc1,pgzvc1,pgxvc2,pgyvc2,pgzvc2,
1095  & pgxvc3,pgyvc3,pgzvc3,ndch1,ndch2,ndch3,
1096  & ikdch1,ikdch2,ikdch3,
1097  & pdq1(4),pdq2(4),pdd1(4),pdd2(4),pdfq1(4)
1098  COMMON /dpar/ aname(210),am(210),ga(210),tau(210),ich(210),
1099  & ibar(210),k1(210),k2(210)
1100  DATA ncpt /0/
1101 *
1102  pch3i = sign(sqrt((ech3i-am(kch3i))*(ech3i+am(kch3i))),pch3)
1103 *
1104 *-------------------- select transverse momenta - diffractive particle
1105 *
1106  tttmin = (ech3-ech3i)**2-(pch3-pch3i)**2
1107  10 CONTINUE
1108 *
1109  xdiff = xdich1+xdich2
1110  btp0 = 3.7d0
1111  alph = 0.24d0
1112  slope = btp0-2.0d0*alph*log(xdiff)
1113  y = rndm(v)
1114  ttt = -log(1.0d0-y)/slope
1115  IF (ipev.GE.2) THEN
1116  WRITE(lout,1009) slope,ttt
1117  1009 FORMAT('VAHMSD/DIFFPT: SLOPE,TTT ',2e10.5)
1118  ENDIF
1119 *
1120  IF (ttt.LE.abs(tttmin)) goto 10
1121  ttt = -abs(ttt)
1122  IF (ipev.GE.2) WRITE(lout,1000) pch3i,tttmin,ttt
1123  1000 FORMAT('VAHMSD/DIFFPT: PCH3I,TTTMIN,TTT ',3f10.5)
1124  pch3f = abs(ttt-(ech3-ech3i)**2+pch3**2+pch3i**2)/(2*pch3i)
1125  IF ((pch3**2).LE.(pch3f**2)) THEN
1126  ncpt = ncpt+1
1127  IF(mod(ncpt,5000).EQ.0) WRITE(lout,1001) ncpt
1128  1001 FORMAT('VAHMSD/DIFFPT: INEFFICIENT PT-SELECTION 1, NCPT=',i8)
1129  goto 10
1130  ENDIF
1131  p3xyf = sqrt(pch3**2-pch3f**2)
1132  CALL dsfecf(sfe,cfe)
1133  pxch3f = p3xyf*cfe
1134  pych3f = p3xyf*sfe
1135 *
1136 *-------------------- select transverse momenta for partons of chain 1
1137 *
1138  IF (ndch1.EQ.-99) THEN
1139  eaq2 = 0.0d0
1140  ptxva2 = 0.0d0
1141  ptyva2 = 0.0d0
1142  plaq2 = 0.0d0
1143  eq1 = 0.0d0
1144  ptxdq1 = 0.0d0
1145  ptydq1 = 0.0d0
1146  plq1 = 0.0d0
1147  pxch1f = 0.0d0
1148  pych1f = 0.0d0
1149  plch1f = 0.0d0
1150  goto 11
1151  ENDIF
1152  pych1f = -abs(pch1/pch3)*pych3f
1153  pxch1f = -abs(pych1f/pych3f)*pxch3f
1154  temp = pxch1f**2+pych1f**2
1155  IF ((pch1**2).LT.temp) THEN
1156  ncpt = ncpt+1
1157  IF(mod(ncpt,5000).EQ.0) WRITE(lout,1002) ncpt
1158  1002 FORMAT('VAHMSD/DIFFPT: INEFFICIENT PT-SELECTION 2, NCPT=',i8)
1159  goto 10
1160  ENDIF
1161  plch1f = sign(sqrt(pch1**2-temp),pch1)
1162  eaq2 = ees*xdich1
1163  temp = abs(eaq2/pch1)
1164  ptxva2 = -pxch1f*temp
1165  ptyva2 = -pych1f*temp
1166  plaq2 = -plch1f*temp
1167  eq1 = ee*xch1
1168  temp = abs(eq1/pch1)
1169  ptxdq1 = pxch1f*temp
1170  ptydq1 = pych1f*temp
1171  plq1 = plch1f*temp
1172 *
1173 *-------------------- select transverse momenta for partons of chain 2
1174 *
1175  11 CONTINUE
1176  IF (ndch2.EQ.-99) THEN
1177  eq2 = 0.0d0
1178  ptxdq2 = 0.0d0
1179  ptydq2 = 0.0d0
1180  plq2 = 0.0d0
1181  eaq1 = 0.0d0
1182  ptxva1 = 0.0d0
1183  ptyva1 = 0.0d0
1184  plaq1 = 0.0d0
1185  pxch2f = 0.0d0
1186  pych2f = 0.0d0
1187  plch2f = 0.0d0
1188  goto 12
1189  ENDIF
1190  pych2f = -abs(pch2/pch3)*pych3f
1191  pxch2f = -abs(pych2f/pych3f)*pxch3f
1192  temp = pxch2f**2+pych2f**2
1193  IF ((pch2**2).LT.temp) THEN
1194  ncpt = ncpt+1
1195  IF(mod(ncpt,5000).EQ.0) WRITE(lout,1003) ncpt
1196  1003 FORMAT('VAHMSD/DIFFPT: INEFFICIENT PT-SELECTION 3, NCPT=',i8)
1197  goto 10
1198  ENDIF
1199  plch2f = sign(sqrt(pch2**2-temp),pch2)
1200  eq2 = ees*xdich2
1201  temp = abs(eq2/pch2)
1202  ptxdq2 = -pxch2f*temp
1203  ptydq2 = -pych2f*temp
1204  plq2 = -plch2f*temp
1205  eaq1 = ee*xch2
1206  temp = abs(eaq1/pch2)
1207  ptxva1 = pxch2f*temp
1208  ptyva1 = pych2f*temp
1209  plaq1 = plch2f*temp
1210 *
1211  12 CONTINUE
1212  ptxch3 = pxch3f
1213  ptych3 = pych3f
1214  pch3 = pch3f
1215 *---------- S. Roesler 11/4/93
1216 * calculate diffractive x-value
1217  IF (ipev.GE.2) THEN
1218  tempx = ptxva2+ptxdq1+ptxdq2+ptxva1
1219  tempy = ptyva2+ptydq1+ptydq2+ptyva1
1220  tempz = plaq2+plq1+plq2+plaq1
1221  tempe = eaq1+eaq2+eq1+eq2
1222  tempm = sqrt(tempe**2-tempx**2-tempy**2-tempz**2)
1223  tempp = tempm**2/((ees+ee)**2)
1224  WRITE(*,*)'diffractive x-value before energy/momentum corr.',
1225  & tempp
1226  ENDIF
1227 *
1228 *---------- S. Roesler 10/26/93
1229 *
1230 * introduce off-shell partons in order to ensure
1231 * energy-momentum conservation
1232 *
1233  diffx = pxch3f+ptxva2+ptxdq1+ptxdq2+ptxva1
1234  diffy = pych3f+ptyva2+ptydq1+ptydq2+ptyva1
1235  diffz = pch3f+plaq2+plq1+plq2+plaq1
1236  IF (ipev.GE.2) THEN
1237  WRITE(lout,1011) diffx,diffy,diffz
1238  1011 FORMAT('VAHMSD/DIFFPT: DIFFX,DIFFY,DIFFZ ',3e15.5)
1239  ENDIF
1240 *
1241  IF ((ndch1.EQ.-99).OR.(ndch1.EQ.1).OR.(ndch1.EQ.-1)) THEN
1242 *
1243  ptxdq2 = ptxdq2-diffx/2.0d0
1244  ptxva1 = ptxva1-diffx/2.0d0
1245 *
1246  ptydq2 = ptydq2-diffy/2.0d0
1247  ptyva1 = ptyva1-diffy/2.0d0
1248 *
1249  plq2 = plq2 -diffz/2.0d0
1250  plaq1 = plaq1 -diffz/2.0d0
1251 *
1252  ELSEIF ((ndch2.EQ.-99).OR.(ndch2.EQ.1).OR.(ndch2.EQ.-1)) THEN
1253 *
1254  ptxva2 = ptxva2-diffx/2.0d0
1255  ptxdq1 = ptxdq1-diffx/2.0d0
1256 *
1257  ptyva2 = ptyva2-diffy/2.0d0
1258  ptydq1 = ptydq1-diffy/2.0d0
1259 *
1260  plaq2 = plaq2 -diffz/2.0d0
1261  plq1 = plq1 -diffz/2.0d0
1262 *
1263  ELSE
1264 *
1265  ptxva2 = ptxva2-diffx/4.0d0
1266  ptxdq1 = ptxdq1-diffx/4.0d0
1267  ptxdq2 = ptxdq2-diffx/4.0d0
1268  ptxva1 = ptxva1-diffx/4.0d0
1269 *
1270  ptyva2 = ptyva2-diffy/4.0d0
1271  ptydq1 = ptydq1-diffy/4.0d0
1272  ptydq2 = ptydq2-diffy/4.0d0
1273  ptyva1 = ptyva1-diffy/4.0d0
1274 *
1275  plaq2 = plaq2 -diffz/4.0d0
1276  plq1 = plq1 -diffz/4.0d0
1277  plq2 = plq2 -diffz/4.0d0
1278  plaq1 = plaq1 -diffz/4.0d0
1279 *
1280  ENDIF
1281 *---------- S. Roesler 11/4/93
1282 * calculate diffractive x-value
1283  IF (ipev.GE.2) THEN
1284  tempx = ptxva2+ptxdq1+ptxdq2+ptxva1
1285  tempy = ptyva2+ptydq1+ptydq2+ptyva1
1286  tempz = plaq2+plq1+plq2+plaq1
1287  tempe = eaq1+eaq2+eq1+eq2
1288  tempm = sqrt(tempe**2-tempx**2-tempy**2-tempz**2)
1289  tempp = tempm**2/((ees+ee)**2)
1290  WRITE(*,*)'diffractive x-value after energy/momentum corr.',
1291  & tempp
1292  ENDIF
1293 *
1294 * recalculate chain masses...
1295 *
1296  amdch1 = sqrt(abs((eaq2+eq1)**2-(ptxva2+ptxdq1)**2
1297  & -(ptyva2+ptydq1)**2-(plaq2+plq1)**2)+1.0d-8)
1298  amdch2 = sqrt(abs((eq2+eaq1)**2-(ptxdq2+ptxva1)**2
1299  & -(ptydq2+ptyva1)**2-(plq2+plaq1)**2)+1.0d-8)
1300 *
1301  IF (ipev.GE.2) THEN
1302  WRITE(lout,1012) amdch1,amdch2
1303  1012 FORMAT('VAHMSD/DIFFPT: AMDCH1,AMDCH2 ',2e10.5)
1304  ENDIF
1305 *
1306 * ...and the Lorentz-parameter gamma
1307 *
1308  gamdc1 = (eaq2+eq1)/amdch1
1309  gamdc2 = (eq2+eaq1)/amdch2
1310  IF (ipev.GE.2) THEN
1311  WRITE(lout,1013) gamdc1,gamdc2
1312  1013 FORMAT('VAHMSD/DIFFPT: GAMDC1,GAMDC2 ',2e10.5)
1313  ENDIF
1314 *
1315 *-------------------- store results in common block /ABRDIF/
1316 *
1317  pdq1(1) = ptxdq1
1318  pdq1(2) = ptydq1
1319  pdq1(3) = plq1
1320  pdq1(4) = eq1
1321  pdd1(1) = ptxva1
1322  pdd1(2) = ptyva1
1323  pdd1(3) = plaq1
1324  pdd1(4) = eaq1
1325  pdfq1(1) = ptxch3
1326  pdfq1(2) = ptych3
1327  pdfq1(3) = pch3
1328  pdfq1(4) = ech3
1329  pdq2(1) = ptxdq2
1330  pdq2(2) = ptydq2
1331  pdq2(3) = plq2
1332  pdq2(4) = eq2
1333  pdd2(1) = ptxva2
1334  pdd2(2) = ptyva2
1335  pdd2(3) = plaq2
1336  pdd2(4) = eaq2
1337  IF (ipev.GE.2) THEN
1338  WRITE(lout,1004)
1339  1004 FORMAT('VAHMSD/DIFFPT: PDQ1,PDD1,PDQ2,PDD2,PDFQ1 ')
1340  WRITE(lout,1005)
1341  & (pdq1(j),pdd1(j),pdq2(j),pdd2(j),pdfq1(j),j=1,4)
1342  1005 FORMAT(12x,5f10.5)
1343  ENDIF
1344  pgxvc1 = (ptxdq1+ptxva2)/amdch1
1345  pgyvc1 = (ptydq1+ptyva2)/amdch1
1346  pgzvc1 = (plq1+plaq2)/amdch1
1347  pgxvc2 = (ptxva1+ptxdq2)/amdch2
1348  pgyvc2 = (ptyva1+ptydq2)/amdch2
1349  pgzvc2 = (plaq1+plq2)/amdch2
1350  pgxvc3 = ptxch3/amdch3
1351  pgyvc3 = ptych3/amdch3
1352  pgzvc3 = pch3/amdch3
1353 *
1354  pxch1f = ptxdq1+ptxva2
1355  pych1f = ptydq1+ptyva2
1356  plch1f = plq1+plaq2
1357  pxch2f = ptxva1+ptxdq2
1358  pych2f = ptyva1+ptydq2
1359  plch2f = plaq1+plq2
1360 *
1361  IF (ipev.GE.2) THEN
1362  WRITE(lout,1006) pgxvc1,pgyvc1,pgzvc1
1363  1006 FORMAT('VAHMSD/DIFFPT: PGXVC1,PGYVC1,PGZVC1 ',3f10.5)
1364  WRITE(lout,1007) pgxvc2,pgyvc2,pgzvc2
1365  1007 FORMAT('VAHMSD/DIFFPT: PGXVC2,PGYVC2,PGZVC2 ',3f10.5)
1366  WRITE(lout,1008) pgxvc3,pgyvc3,pgzvc3
1367  1008 FORMAT('VAHMSD/DIFFPT: PGXVC3,PGYVC3,PGZVC3 ',3f10.5)
1368  ENDIF
1369  IF ((ndch1.EQ.-99).AND.(ndch2.EQ.-99)) THEN
1370  WRITE(lout,*) 'REJECT IN DIFFPT: NO CHAINS CREATED'
1371  iirej = 1
1372  RETURN
1373  ENDIF
1374 *
1375 *-------------------- store results in common block /HKKEVT/
1376 *
1377  phkk(1,nhkk+9) = pxch1f
1378  phkk(2,nhkk+9) = pych1f
1379  phkk(3,nhkk+9) = plch1f
1380  phkk(4,nhkk+9) = ech1
1381 *
1382  phkk(1,nhkk+12) = pxch2f
1383  phkk(2,nhkk+12) = pych2f
1384  phkk(3,nhkk+12) = plch2f
1385  phkk(4,nhkk+12) = ech2
1386 *
1387  RETURN
1388  END
1389 *
1390 *===diffch================================================================
1391 *
1392  SUBROUTINE diffch ( XSEA, IFSEA, XPAR, IFPARA, IFPARB,
1393  & xsea2, ifsea2, xpar12, xpar22, ee,
1394  & rm, ech, pch, gamma, beta,
1395  & ndch, ich, ees, nuno, iirej, iopt)
1396 
1397 **************************************************************************
1398 * Version November 1993 by Stefan Roesler *
1399 * University of Leipzig *
1400 * This subroutine is called by VAHMSD (HMSD) and calculates the kine- *
1401 * matical parameters of chain IOPT from sampled x-values and flavors *
1402 * in CMS. *
1403 **************************************************************************
1404 
1405  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1406  SAVE
1407  parameter(lout=6,llook=9)
1408  CHARACTER*8 aname
1409  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
1410  COMMON /diffra/ isingd,idiftp,ioudif,iflagd
1411  COMMON /dpar/ aname(210),am(210),ga(210),tau(210),iich(210),
1412  & ibar(210),k1(210),k2(210)
1413  common/dinpda/imps(6,6),imve(6,6),ib08(6,21),ib10(6,21),ia08(6,21)
1414  &,ia10(6,21),a1,b1,b2,b3,lt,le,bet,as,b8,ame,diq,isu
1415  DATA ncnoth,ncwd /0,0/
1416 *
1417  goto(1,2,3) iopt
1418 *
1419 *-------------------- kinematical parameters of a q-aq (in any order) chain
1420 * XSEA - XPAR, IFSEA - IFPAR
1421 *
1422  1 CONTINUE
1423  iirej = 0
1424  ifq = ifpara
1425  ifaq = iabs(ifsea)
1426  IF (ifq.LT.0) THEN
1427  ifq = ifsea
1428  ifaq = iabs(ifpara)
1429  ENDIF
1430  ips = imps(ifaq,ifq)
1431  iv = imve(ifaq,ifq)
1432  rmps = am(ips)
1433  rmv = am(iv)
1434  rmbb = rmv+3.0d-1
1435  ndch = 0
1436 *
1437  IF ((xsea.LE.0.0d0).OR.(xpar.LE.0.0d0)) THEN
1438  WRITE(lout,*) 'REJECTION IN DIFFCH: 1, XSEA,XPAR ',xsea,xpar
1439  iirej=1
1440  RETURN
1441  ENDIF
1442  rm = 2.0d0*sqrt(ees*xsea*ee*xpar)
1443  IF (ipev.GE.2) THEN
1444  WRITE(lout,1000) 'XSEA, XSEA2, XPAR, XPAR12, RM, EE, EES'
1445  WRITE(lout,1001) xsea, xsea2, xpar, xpar12, rm, ee, ees
1446  ENDIF
1447 *
1448  IF (rm.LT.rmps) THEN
1449 * produce nothing
1450  ndch = -99
1451  xsea2 = xsea2 +xsea
1452  xpar12 = xpar12+xpar
1453  ifsea2 = ifpara
1454  xsea = 0.0d0
1455  xpar = 0.0d0
1456  rm = 1.0d-4
1457  ech = 0.0d0
1458  pch = 0.0d0
1459  gamma = 1.0d0
1460  beta = 0.0d0
1461  ncnoth = ncnoth+1
1462  IF (mod(ncnoth,5000).EQ.0) WRITE(lout,1002) ncnoth
1463  1002 FORMAT('VAHMSD/DIFFCH: PRODUCE NOTHING 1, NCNOTH=',i8)
1464  RETURN
1465 *
1466  ELSE IF (rm.LT.rmv) THEN
1467 * produce RMPS
1468  rm = rmps
1469  ndch = -1
1470  ich = ips
1471  xsq = rm/(2.0d0*sqrt(ee*ees))
1472  xseaol = xsea
1473  xsea = xsq**2/xpar
1474  xpar22 = xpar22+xseaol-xsea
1475 *
1476  ELSE IF (rm.LT.rmbb) THEN
1477 * produce RMV
1478  rm = rmv
1479  ndch = 1
1480  ich = iv
1481  xsq = rm/(2.0d0*sqrt(ee*ees))
1482  xseaol = xsea
1483  xsea = xsq**2/xpar
1484  xpar22 = xpar22+xseaol-xsea
1485  ENDIF
1486 *
1487  IF (idiftp.EQ.1) THEN
1488  pch = ees*xsea-ee*xpar
1489  IF (pch.GE.0.0d0) THEN
1490  ncwd = ncwd+1
1491  IF (mod(ncwd,500).EQ.0) WRITE(lout,1003) ncwd
1492  1003 FORMAT('VAHMSD/DIFFCH: WRONG SIGN OF PCH (1), NCWD=',i8)
1493  iirej = 1
1494  ENDIF
1495  ELSE IF (idiftp.EQ.2) THEN
1496  pch = ee*xpar-ees*xsea
1497  IF (pch.LE.0.0d0) THEN
1498  ncwd = ncwd+1
1499  IF (mod(ncwd,500).EQ.0) WRITE(lout,1004) ncwd
1500  1004 FORMAT('VAHMSD/DIFFCH: WRONG SIGN OF PCH (2), NCWD=',i8)
1501  iirej = 1
1502  ENDIF
1503  ELSE
1504  WRITE(lout,*) 'ERROR IN VAHMSD/DIFFCH, IDIFTP=', idiftp
1505  ENDIF
1506  IF (ipev.GE.2) THEN
1507  WRITE(lout,1000) 'EES,XSEA,EE,XPAR'
1508  WRITE(lout,1001) ees,xsea,ee,xpar
1509  ENDIF
1510  ech = ees*xsea+ee*xpar
1511  gamma = ech/rm
1512  beta = pch/rm
1513  IF (ipev.GE.2) THEN
1514  WRITE(lout,1000) 'RM,ECH,PCH,NDCH'
1515  WRITE(lout,1001) rm,ech,pch
1516  WRITE(lout,'(I4)') ndch
1517  ENDIF
1518  1000 FORMAT('VAHMSD/DIFFCH-CHAIN Q-AQ',a34)
1519  1001 FORMAT(10f10.5)
1520  RETURN
1521 *
1522 *-------------------- kinematical parameters of a q-qq or aq-aqaq (in any
1523 * order) chain
1524 * XSEA - XPAR, IFSEA - IFPARA/IFPARB
1525 *
1526  2 CONTINUE
1527  iirej = 0
1528  CALL dbklas(ifsea, ifpara, ifparb, i8, i10)
1529  rm8 = am( i8)
1530  rm10 = am(i10)
1531  rmbb = rm10+3.0d-1
1532  ndch = 0
1533 *
1534  IF ((xsea.LE.0.0d0).OR.(xpar.LE.0.0d0)) THEN
1535  iirej=1
1536  WRITE(lout,*) 'REJECTION IN DIFFCH: 3, XSEA,XPAR ',xsea,xpar
1537  RETURN
1538  ENDIF
1539  rm = 2.0d0*sqrt(ees*xsea*ee*xpar)
1540  IF (ipev.GE.2) THEN
1541  WRITE(lout,2000) 'XSEA, XSEA2, XPAR, XPAR12, RM, EE, EES'
1542  WRITE(lout,2001) xsea, xsea2, xpar, xpar12, rm, ee, ees
1543  ENDIF
1544 *
1545  IF (rm.LT.rm10) THEN
1546 * reject
1547  iirej = 1
1548  ncnoth = ncnoth+1
1549  IF (mod(ncnoth,20).EQ.0) WRITE(lout,2002) ncnoth
1550  2002 FORMAT('VAHMSD/DIFFCH: PRODUCE NOTHING 2, NCNOTH=',i8)
1551  RETURN
1552  ELSE IF (rm.LT.rmbb) THEN
1553 * produce RM10
1554  ndch = 1
1555  rm = rm10
1556  ich = i10
1557  xsq = rm/(2.0d0*sqrt(ee*ees))
1558  xseaol = xsea
1559  xsea = xsq**2/xpar
1560  xpar22 = xpar22+xseaol-xsea
1561  ENDIF
1562 *
1563  IF (idiftp.EQ.1) THEN
1564  pch = ees*xsea-ee*xpar
1565  IF (pch.GE.0.0d0) THEN
1566  ncwd = ncwd+1
1567  IF (mod(ncwd,500).EQ.0) WRITE(lout,2003) ncwd
1568  2003 FORMAT('VAHMSD/DIFFCH: WRONG SIGN OF PCH (3), NCWD=',i8)
1569  iirej = 1
1570  ENDIF
1571  ELSE IF (idiftp.EQ.2) THEN
1572  pch = ee*xpar-ees*xsea
1573  IF (pch.LE.0.0d0) THEN
1574  ncwd = ncwd+1
1575  IF (mod(ncwd,500).EQ.0) WRITE(lout,2004) ncwd
1576  2004 FORMAT('VAHMSD/DIFFCH: WRONG SIGN OF PCH (4), NCWD=',i8)
1577  iirej = 1
1578  ENDIF
1579  ELSE
1580  WRITE(lout,*) 'ERROR IN VAHMSD/DIFFCH, IDIFTP=', idiftp
1581  ENDIF
1582  ech = ees*xsea+ee*xpar
1583  gamma = ech/rm
1584  beta = pch/rm
1585  IF (ipev.GE.2) THEN
1586  WRITE(lout,2000) 'RM,ECH,PCH,NDCH'
1587  WRITE(lout,2001) rm,ech,pch
1588  WRITE(lout,'(I4)') ndch
1589  ENDIF
1590  2000 FORMAT('VAHMSD/DIFFCH-CHAIN Q-QQ',a34)
1591  2001 FORMAT(10f10.5)
1592  RETURN
1593 *
1594 *-------------------- kinematical parameters of a baryon/meson
1595 *
1596  3 CONTINUE
1597  iirej = 0
1598  IF (ifparb.EQ.99) THEN
1599  IF (ifsea.LT.0) THEN
1600  ifaq = iabs(ifsea)
1601  ifq = ifpara
1602  ELSE
1603  ifaq = iabs(ifpara)
1604  ifq = ifsea
1605  ENDIF
1606  im = imps(ifaq,ifq)
1607  ELSE
1608  CALL dbklas(ifsea, ifpara, ifparb, im, idum)
1609  ENDIF
1610  rm = am(im)
1611  ndch = -1
1612  IF ((xsea.LE.0.0d0).OR.(xpar.LE.0.0d0)) THEN
1613  iirej=1
1614  WRITE(lout,*) 'REJECTION IN DIFFCH: 5, XSEA,XPAR ',xsea,xpar
1615  RETURN
1616  ENDIF
1617  IF (ipev.GE.2) THEN
1618  WRITE(lout,3000) 'XSEA, XPAR, RM'
1619  WRITE(lout,3001) xsea, xpar, rm
1620  ENDIF
1621  ech = ee*(xsea+xpar)
1622  IF (idiftp.EQ.1) THEN
1623  pch = sqrt((ech-rm)*(ech+rm))
1624  ELSE IF (idiftp.EQ.2) THEN
1625  pch = -sqrt((ech-rm)*(ech+rm))
1626  ELSE
1627  WRITE(lout,*) 'ERROR IN VAHMSD/DIFFCH, IDIFTP=', idiftp
1628  ENDIF
1629  gamma = ech/rm
1630  beta = pch/rm
1631  IF (ipev.GE.2) THEN
1632  WRITE(lout,3000) 'ECH, PCH, GAMMA, BETA'
1633  WRITE(lout,3001) ech, pch, gamma, beta
1634  ENDIF
1635  3000 FORMAT('VAHMSD/DIFFCH-BARYON/MESON ',a27)
1636  3001 FORMAT(10f10.5)
1637  RETURN
1638  END
1639 *
1640 *===valmsd===============================================================*
1641 *
1642  SUBROUTINE valmsd(ITAPOI,ECM,KPROJ,KTARG,IREJ)
1643 
1644 **************************************************************************
1645 * Version November 1993 by Stefan Roesler *
1646 * University of Leipzig *
1647 * This subroutine selects flavors and 4-momenta of partons in low-mass *
1648 * single diffractive chains. *
1649 **************************************************************************
1650 
1651  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1652  SAVE
1653  parameter(lout=6,llook=9)
1654  parameter(nmxhkk= 89998)
1655  CHARACTER*8 aname
1656  COMMON /hkkevt/ nhkk,nevhkk, isthkk(nmxhkk), idhkk(nmxhkk),
1657  & jmohkk(2,nmxhkk),jdahkk(2,nmxhkk),phkk(5,nmxhkk),
1658  & vhkk(4,nmxhkk), whkk(4,nmxhkk)
1659  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
1660  COMMON /diffra/ isingd,idiftp,ioudif,iflagd
1661  COMMON /abrdif/ xdq1,xdq2,xddq1,xddq2,
1662  & ikvq1,ikvq2,ikd1q1,ikd2q1,ikd1q2,ikd2q2,
1663  & idiffp,idifap,
1664  & amdch1,amdch2,amdch3,gamdc1,gamdc2,gamdc3,
1665  & pgxvc1,pgyvc1,pgzvc1,pgxvc2,pgyvc2,pgzvc2,
1666  & pgxvc3,pgyvc3,pgzvc3,ndch1,ndch2,ndch3,
1667  & ikdch1,ikdch2,ikdch3,
1668  & pdq1(4),pdq2(4),pdd1(4),pdd2(4),pdfq1(4)
1669  COMMON /dpar/ aname(210),am(210),ga(210),tau(210),ich(210),
1670  & ibar(210),k1(210),k2(210)
1671  common/dinpda/imps(6,6),imve(6,6),ib08(6,21),ib10(6,21),ia08(6,21)
1672  &,ia10(6,21),a1,b1,b2,b3,lt,le,bet,as,b8,ame,diq,isu
1673  COMMON /trafop/ gamp,bgamp,betp
1674  COMMON /enerin/ eproj,etarg
1675  COMMON /sdflag/ isd
1676  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
1677  COMMON /xdidid/xdidi
1678 *
1679  dimension mquark(3,30),ihkkq(-6:6),ihkkqq(-3:3,-3:3),
1680  & idx(-4:4)
1681  DATA idx /-4,-3,-1,-2,0,2,1,3,4/
1682  DATA ihkkq /-6,-5,-4,-3,-1,-2,0,2,1,3,4,5,6/
1683  DATA ihkkqq/-3301,-3103,-3203,0, 0,0,0,
1684  & -3103,-1103,-2103,0, 0, 0, 0,
1685  & -3203,-2103,-2203,0, 0, 0, 0,
1686  & 0, 0, 0,0, 0, 0, 0,
1687  & 0, 0, 0,0,2203,2103,3202,
1688  & 0, 0, 0,0,2103,1103,3103,
1689  & 0, 0, 0,0,3203,3103,3303/
1690 *
1691 *----------------------------------- quark content of hadrons:
1692 * 1, 2, 3, 4 - u, d, s, c
1693 * -1,-2,-3,-4 - au,ad,as,ac
1694 *
1695  DATA mquark/
1696  & 1,1,2, -1,-1,-2, 0,0,0, 0,0,0, 0,0,0,
1697  & 0,0,0, 0,0,0, 1,2,2, -1,-2,-2, 0,0,0,
1698  & 0,0,0, 0,0,0, 1,-2,0, 2,-1,0, 1,-3,0,
1699  & 3,-1,0, 1,2,3, -1,-2,-3, 0,0,0, 2,2,3,
1700  & 1,1,3, 1,2,3, 1,-1,0, 2,-3,0, 3,-2,0,
1701  & 2,-2,0, 3,-3,0, 0,0,0, 0,0,0, 0,0,0/
1702  DATA unon/2.0/
1703  DATA ncrej, ncpt /0, 0/
1704 *
1705  isd = 2
1706  irej = 0
1707  iirej = 0
1708  ibproj = ibar(kproj)
1709  ibtarg = ibar(ktarg)
1710  eproj = (am(kproj)**2-am(ktarg)**2+ecm**2)/(2.0d0*ecm)
1711  etarg = (am(ktarg)**2-am(kproj)**2+ecm**2)/(2.0d0*ecm)
1712  IF(ipev.GE.2) WRITE(lout,1014)eproj,etarg
1713  1014 FORMAT('VALMSD: EPROJ,ETARG',2f10.5)
1714  IF (ibtarg.LE.0) THEN
1715  WRITE(lout,1001) ibtarg
1716  1001 FORMAT('VALMSD: NO LMSD FOR TARGET WITH BARYON-CHARGE',i4)
1717  iirej = 1
1718  goto 9999
1719  ENDIF
1720  iqp1 = mquark(1,kproj)
1721  iqp2 = mquark(2,kproj)
1722  iqp3 = mquark(3,kproj)
1723  iqt1 = mquark(1,ktarg)
1724  iqt2 = mquark(2,ktarg)
1725  iqt3 = mquark(3,ktarg)
1726  IF(ipev.GE.2) WRITE(lout,1002)
1727  & ibproj,ibtarg,iqp1,iqp2,iqp3,iqt1,iqt2,iqt3
1728  1002 FORMAT('VALMSD: IBPROJ,IBTARG,IQP1,IQP2,IQP3,IQT1,IQT2,IQT3 ',8i4)
1729 *
1730  IF (ibproj.NE.0) THEN
1731 *
1732 *-------------------- q-qq (aq-aqaq) - flavors of projectile (baryon)
1733 *
1734  isam = 1.0d0+2.999d0*rndm(v)
1735  goto(10,11,12) isam
1736  10 CONTINUE
1737  iqp = iqp1
1738  idiqp1 = iqp2
1739  idiqp2 = iqp3
1740  goto 13
1741  11 CONTINUE
1742  iqp = iqp2
1743  idiqp1 = iqp1
1744  idiqp2 = iqp3
1745  goto 13
1746  12 CONTINUE
1747  iqp = iqp3
1748  idiqp1 = iqp1
1749  idiqp2 = iqp2
1750  13 CONTINUE
1751 *
1752  ELSE IF (ibproj.EQ.0) THEN
1753 *
1754 *-------------------- q-aq - flavors of projectile (meson)
1755 *
1756  isam = 1.0d0+1.999d0*rndm(v)
1757  goto(14,15) isam
1758  14 CONTINUE
1759  iqp = iqp1
1760  idiqp1 = iqp2
1761  goto 16
1762  15 CONTINUE
1763  iqp = iqp2
1764  idiqp1 = iqp1
1765  16 CONTINUE
1766  idiqp2 = 0
1767 *
1768  ENDIF
1769 *
1770 *-------------------- q-qq - flavors of target (baryon)
1771 *
1772  isam = 1.0d0+2.999d0*rndm(v)
1773  goto(17,18,19) isam
1774  17 CONTINUE
1775  iqt = iqt1
1776  idiqt1 = iqt2
1777  idiqt2 = iqt3
1778  goto 20
1779  18 CONTINUE
1780  iqt = iqt2
1781  idiqt1 = iqt1
1782  idiqt2 = iqt3
1783  goto 20
1784  19 CONTINUE
1785  iqt = iqt3
1786  idiqt1 = iqt1
1787  idiqt2 = iqt2
1788  20 CONTINUE
1789 *
1790  ikvq1 = iqp
1791  ikd1q1 = idiqp1
1792  ikd2q1 = idiqp2
1793 *
1794  ikvq2 = iqt
1795  ikd1q2 = idiqt1
1796  ikd2q2 = idiqt2
1797 *
1798  IF (ipev.GE.2) WRITE(lout,1003)
1799  & ikvq1,ikd1q1,ikd2q1,ikvq2,ikd1q2,ikd2q2
1800  1003 FORMAT('VALMSD: IKVQ1,IKD1Q1,IKD2Q1,IKVQ2,IKD1Q2,IKD2Q2 ',6i4)
1801 *
1802 *-------------------- IDIFTP = 1 target (backward) hadron excited
1803 * IDIFTP = 2 projectile (forward) hadron excited
1804 *
1805  idiftp = 1.0d0+1.999d0*rndm(v)
1806 *-------------------- S.Roesler 5/26/93
1807  IF ((isingd.EQ.3).OR.(isingd.EQ.7)) idiftp = 1
1808  IF ((isingd.EQ.4).OR.(isingd.EQ.8)) idiftp = 2
1809 *
1810  IF (ipev.GE.2) WRITE(lout,1004) idiftp
1811  1004 FORMAT('VALMSD: IDIFTP ',i4)
1812  IF ((idiftp.NE.1).AND.(idiftp.NE.2)) THEN
1813  IF (ipev.GE.2) WRITE(lout,'(A19)') 'VALMSD-ERROR: IDIFTP'
1814  goto 9999
1815  ENDIF
1816 *
1817 *-------------------- diffractive mass
1818 *
1819 31 CONTINUE
1820  amo = 1.5d0
1821  r = rndm(v)
1822  IF ((ibproj.EQ.0).AND.(idiftp.EQ.1))
1823  & amo = 1.5d0*r+2.83d0*(1.0d0-r)
1824  IF ((ibproj.EQ.0).AND.(idiftp.EQ.2)) amo = 1.0d0
1825  sam = 1.0d0
1826  IF (ecm.LE.300.0d0) sam = 1.0d0-exp(-((ecm/200.0d0)**4))
1827  r = rndm(v)*sam
1828  amax= (1.0d0-sam)*sqrt(0.1d0*ecm**2)+sam*sqrt(400.0d0)
1829  amu = r*sqrt(100.0d0)+(1.0d0-r)*amax
1830  IF (ibproj.EQ.0) THEN
1831 C----------------- change mass-cuts
1832 C AMAX= (1.0D0-SAM)*SQRT(0.25D0*ECM**2)+SAM*SQRT(400.0D0)
1833  amax= (1.0d0-sam)*sqrt(0.15d0*ecm**2)+sam*sqrt(400.0d0)
1834  amu = r*sqrt(100.0d0)+(1.0d0-r)*amax
1835  ENDIF
1836 C---------------------------------------------------------------
1837 C
1838 C j.r. test 2/94
1839 C
1840 C---------------------------------------------------------------
1841  amu=2.d0*amu
1842 C---------------------------------------------------------------
1843  r = rndm(v)
1844  IF (ecm.LE.50.0d0) THEN
1845  amdiff = amo*(amu/amo)**r
1846  ELSE
1847  a = 0.7d0
1848  IF (ecm.LE.300.0d0) a = 0.7d0*(1.0d0-exp(-((ecm/100.0d0)**2)))
1849  amdiff = 1.0d0/((r/(amu**a)+(1.0d0-r)/(amo**a))
1850  & **(1.0d0/a))
1851  ENDIF
1852  IF(amdiff.GT.0.5d0*ecm)go to 31
1853  IF (iouxev.GE.2) WRITE(lout,1005) amdiff
1854  1005 FORMAT('VALMSD: AMDIFF',e10.5)
1855  xdidi=amdiff**2/ecm**2
1856  IF (iouxev.GE.2) WRITE(lout,*)'LM AMDIFF,XDIDI ',amdiff,xdidi
1857 *
1858 *-------------------- kinematical parameters of
1859 * diffractive projectile (IDIFTP = 1)
1860 * or diffractive target (IDIFTP = 2)
1861 *
1862  IF (idiftp.EQ.1) THEN
1863  amdch3 = am(kproj)
1864  ELSE IF (idiftp.EQ.2) THEN
1865  amdch3 = am(ktarg)
1866  ENDIF
1867  ech3 = (ecm**2+amdch3**2-amdiff**2)/(2.0d0*ecm)
1868  IF (ech3.LE.amdch3) THEN
1869  ncmerr = ncmerr+1
1870  IF (mod(ncmerr,2200).EQ.0)
1871  & WRITE(lout,*) 'LMSD: INEFFICIENT SELECTION OF AMDIFF(1),
1872  & NCMERR = ',ncmerr
1873  goto 31
1874  ENDIF
1875 *
1876  pch3 = sqrt(abs(ech3-amdch3))*sqrt(ech3+amdch3)
1877  IF (idiftp.EQ.2) pch3 = -pch3
1878  ndch3 = -1
1879  gamdc3 = ech3/amdch3
1880  pgvc3 = pch3/amdch3
1881  IF (iouxev.GE.2) WRITE(lout,1006) amdch3,ech3,pch3,gamdc3,pgvc3
1882  1006 FORMAT('VALMSD: AMDCH3,ECH3,PCH3,GAMDC3,PGVC3',5e15.5)
1883 *
1884  30 CONTINUE
1885  b33 = 8.0d0
1886  es = -2.0d0/(b33**2)*log(rndm(v)*rndm(v))
1887  hps = sqrt(es*es+2.0d0*es*0.94d0)
1888  CALL dsfecf(sfe,cfe)
1889  pxch3 = hps*cfe
1890  pych3 = hps*sfe
1891  ptch3 = sqrt(pxch3**2+pych3**2)
1892  IF (ptch3.GT.abs(pch3)) THEN
1893  ncpt = ncpt+1
1894  IF (mod(ncpt,500).EQ.0) WRITE(lout,1007) ncpt
1895  1007 FORMAT('VALMSD: INEFFICIENT PT-SELECTION 1, NCPT=',i8)
1896  goto 30
1897  ENDIF
1898  plch3 = sign(sqrt(abs(pch3-ptch3))*sqrt(abs(pch3+ptch3)),pch3)
1899  IF (iouxev.GE.2) WRITE(lout,1008) es,hps,pxch3,pych3,plch3
1900  1008 FORMAT('VALMSD: ES,HPS,PXCH3,PYCH3,PLCH3',5e15.5)
1901 *
1902 *-------------------- no chain 1
1903 *
1904  ndch1 = -99
1905  amdch1 = 0.0d0
1906  ech1 = 0.0d0
1907  pch1 = 0.0d0
1908  gamdc1 = 1.0d0
1909  pgvc1 = 0.0d0
1910  pgxvc1 = 0.0d0
1911  pgyvc1 = 0.0d0
1912  pgzvc1 = 0.0d0
1913  DO 40 i=1,4
1914  pdd2(i) = 0.0d0
1915  pdq1(i) = 0.0d0
1916  40 CONTINUE
1917 *
1918 *-------------------- kinematical parameters of
1919 * excited target (IDIFTP = 1)
1920 * or excited projectile (IDIFTP = 2)
1921 *
1922  amdch2 = amdiff
1923  ech2 = (ecm**2+amdch2**2-amdch3**2)/(2.0d0*ecm)
1924  pch2 = -sign(sqrt(abs(ech2-amdch2))*sqrt(ech2+amdch2),pch3)
1925  ndch2 = 0
1926  gamdc2 = ech2/amdch2
1927  pgvc2 = pch2/amdch2
1928 *
1929 *-------------------- (IDIFFP,IDIFAP in analogy to HMSD in VAHMSD)
1930 * IDIFTP = 1 : IKD1Q2/IKD2Q2 - IDIFFP = IKVQ2
1931 * IDIFTP = 2 : projectile - baryon
1932 * IKD1Q1/IKD2Q1 - IDIFFP = IKVQ1
1933 * projectile - antibaryon
1934 * IKD1Q1/IKD2Q1 - IDIFAP = IKVQ1
1935 * projectile - meson
1936 * IKD1Q1 > 0 - IDIFAP = IKVQ1
1937 * IKD1Q1 < 0 - IDIFFP = IKVQ1
1938 *
1939  IF (idiftp.EQ.1) idiffp = ikvq2
1940  IF (idiftp.EQ.2) THEN
1941  IF (ibproj.GT.0) idiffp = ikvq1
1942  IF (ibproj.LT.0) idifap = ikvq1
1943  IF ((ibproj.EQ.0).AND.(ikd1q1.GT.0)) idifap = ikvq1
1944  IF ((ibproj.EQ.0).AND.(ikd1q1.LT.0)) idiffp = ikvq1
1945  ENDIF
1946  IF (iouxev.GE.2) WRITE(lout,1009) amdch2,ech2,pch2,gamdc2,pgvc2,
1947  & idiffp,idifap
1948  1009 FORMAT('VALMSD: AMDCH2,ECH2,PCH2,GAMDC2,PGVC2,IDIFFP,IDIFAP',
1949  & 5e15.5,2i4)
1950  pxch2 = -pxch3
1951  pych2 = -pych3
1952  ptch2 = sqrt(pxch2**2+pych2**2)
1953  IF (ptch2.GT.abs(pch2)) THEN
1954  ncpt = ncpt+1
1955  IF (mod(ncpt,500).EQ.0) WRITE(lout,1010) ncpt
1956  1010 FORMAT('VALMSD: INEFFICIENT PT-SELECTION 2, NCPT=',i8)
1957  goto 30
1958  ENDIF
1959  plch2 = sign(sqrt(abs(pch2-ptch2))*sqrt(abs(pch2+ptch2)),pch2)
1960  IF (iouxev.GE.2) WRITE(lout,1011) pxch2,pych2,plch2
1961  1011 FORMAT('VALMSD: PXCH2,PYCH2,PLCH2',3e15.5)
1962  cc = amdch2/(2.0d0*ech2)
1963 *--------------------------- S.R. 11/12/92
1964  IF (cc.GE.0.5d0) THEN
1965  ncmerr = ncmerr+1
1966  IF (mod(ncmerr,200).EQ.0)
1967  & WRITE(lout,*) 'LMSD: INEFFICIENT SELECTION OF AMDIFF(2),
1968  & NCMERR = ',ncmerr
1969  goto 31
1970  ENDIF
1971 *
1972  xdiq = 0.5d0+sqrt(abs((0.5d0-cc)*(0.5d0+cc)))
1973  xq = 1.0d0-xdiq
1974  ediq = ech2*xdiq
1975  eq = ech2*xq
1976  temp = abs(ediq/pch2)
1977  pxdiq = pxch2*temp
1978  pydiq = pych2*temp
1979  pldiq = plch2*temp
1980  temp = abs(eq/pch2)
1981  pxq = -pxch2*temp
1982  pyq = -pych2*temp
1983  plq = -plch2*temp
1984 *
1985 *-------------------- store results in COMMON-block /ABRDIF/
1986 *
1987  pdq2(1) = pxq
1988  pdq2(2) = pyq
1989  pdq2(3) = plq
1990  pdq2(4) = eq
1991  pdd1(1) = pxdiq
1992  pdd1(2) = pydiq
1993  pdd1(3) = pldiq
1994  pdd1(4) = ediq
1995  pdfq1(1) = pxch3
1996  pdfq1(2) = pych3
1997  pdfq1(3) = plch3
1998  pdfq1(4) = ech3
1999  pgxvc2 = pxch2/amdch2
2000  pgyvc2 = pych2/amdch2
2001  pgzvc2 = plch2/amdch2
2002  pgxvc3 = pxch3/amdch3
2003  pgyvc3 = pych3/amdch3
2004  pgzvc3 = plch3/amdch3
2005 *
2006  IF (ipev.GE.2) THEN
2007  WRITE(lout,*) 'VALMSD: PDQ1, PDD2, PDD1, PDQ2, PDFQ1'
2008  WRITE(lout,1012)
2009  & (pdq1(i),pdd2(i),pdd1(i),pdq2(i),pdfq1(i),i=1,4)
2010  1012 FORMAT(5e15.5)
2011  WRITE(lout,1013) 'PGXVC1,PGYVC1,PGZVC1,GAMDC1',
2012  & pgxvc1,pgyvc1,pgzvc1,gamdc1
2013  WRITE(lout,1013) 'PGXVC2,PGYVC2,PGZVC2,GAMDC2',
2014  & pgxvc2,pgyvc2,pgzvc2,gamdc2
2015  WRITE(lout,1013) 'PGXVC3,PGYVC3,PGZVC3,GAMDC3',
2016  & pgxvc3,pgyvc3,pgzvc3,gamdc3
2017  1013 FORMAT(a27,4e15.5)
2018  ENDIF
2019 *
2020 *------------------- store results in common block /HKKEVT/
2021 *
2022 * partonen of projectile/target
2023 *
2024  nhkk = nhkk + 1
2025  isthkk(nhkk) = 23-idiftp
2026  IF (idiftp.EQ.1) idhkk(nhkk) = ihkkq(idx(ikvq2))
2027  IF (idiftp.EQ.2) idhkk(nhkk) = ihkkq(idx(ikvq1))
2028  IF (idiftp.EQ.1) jmohkk(1,nhkk) = itapoi
2029  IF (idiftp.EQ.2) jmohkk(1,nhkk) = 1
2030  jmohkk(2,nhkk) = 0
2031  jdahkk(1,nhkk) = 0
2032  jdahkk(2,nhkk) = 0
2033  phkk(1,nhkk) = 0.0d0
2034  phkk(2,nhkk) = 0.0d0
2035  phkk(3,nhkk) = xq
2036  phkk(4,nhkk) = xq
2037  phkk(5,nhkk) = 0.0d0
2038  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))
2039  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))
2040  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
2041  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
2042 *
2043  nhkk = nhkk + 1
2044  isthkk(nhkk) = 23-idiftp
2045  IF (idiftp.EQ.1) idhkk(nhkk) = ihkkqq(idx(ikd1q2),idx(ikd2q2))
2046  IF (idiftp.EQ.2) idhkk(nhkk) = ihkkqq(idx(ikd1q1),idx(ikd2q1))
2047  IF (idiftp.EQ.1) jmohkk(1,nhkk) = itapoi
2048  IF (idiftp.EQ.2) jmohkk(1,nhkk) = 1
2049  jmohkk(2,nhkk) = 0
2050  jdahkk(1,nhkk) = 0
2051  jdahkk(2,nhkk) = 0
2052  phkk(1,nhkk) = 0.0d0
2053  phkk(2,nhkk) = 0.0d0
2054  phkk(3,nhkk) = xdiq
2055  phkk(4,nhkk) = xdiq
2056  phkk(5,nhkk) = 0.0d0
2057  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))
2058  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))
2059  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
2060  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
2061 *
2062 * ends of chain 2
2063 *
2064  nhkk = nhkk + 1
2065  isthkk(nhkk) = 100+isthkk(nhkk-2)
2066  idhkk(nhkk) = idhkk(nhkk-2)
2067  jmohkk(1,nhkk) = nhkk-2
2068  jmohkk(2,nhkk) = jmohkk(1,nhkk-2)
2069  jdahkk(1,nhkk) = nhkk+2
2070  jdahkk(2,nhkk) = nhkk+2
2071  phkk(1,nhkk) = pdq2(1)
2072  phkk(2,nhkk) = pdq2(2)
2073  phkk(3,nhkk) = pdq2(3)
2074  phkk(4,nhkk) = pdq2(4)
2075  phkk(5,nhkk) = 0.0d0
2076 C Add position of parton in hadron
2077  CALL qinnuc(xxpp,yypp)
2078  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))+xxpp
2079  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))+yypp
2080  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
2081  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
2082 *
2083  nhkk = nhkk + 1
2084  isthkk(nhkk) = 100+isthkk(nhkk-2)
2085  idhkk(nhkk) = idhkk(nhkk-2)
2086  jmohkk(1,nhkk) = nhkk-2
2087  jmohkk(2,nhkk) = jmohkk(1,nhkk-2)
2088  jdahkk(1,nhkk) = nhkk+1
2089  jdahkk(2,nhkk) = nhkk+1
2090  phkk(1,nhkk) = pdd1(1)
2091  phkk(2,nhkk) = pdd1(2)
2092  phkk(3,nhkk) = pdd1(3)
2093  phkk(4,nhkk) = pdd1(4)
2094  phkk(5,nhkk) = 0.0d0
2095 C Add position of parton in hadron
2096  CALL qinnuc(xxpp,yypp)
2097  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))+xxpp
2098  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))+yypp
2099  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
2100  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
2101 *
2102 * chain 2
2103 *
2104  nhkk = nhkk + 1
2105  isthkk(nhkk) = 177
2106  idhkk(nhkk) = 88888
2107  jmohkk(1,nhkk) = nhkk-2
2108  jmohkk(2,nhkk) = nhkk-1
2109  jdahkk(1,nhkk) = 0
2110  jdahkk(2,nhkk) = 0
2111  phkk(1,nhkk) = pxch2
2112  phkk(2,nhkk) = pych2
2113  phkk(3,nhkk) = plch2
2114  phkk(4,nhkk) = ech2
2115  phkk(5,nhkk) = amdch2
2116  vhkk(1,nhkk) = vhkk(1,nhkk-1)
2117  vhkk(2,nhkk) = vhkk(2,nhkk-1)
2118  vhkk(3,nhkk) = vhkk(3,nhkk-1)
2119  IF ((betp.NE.0.0d0).AND.(bgamp.NE.0.0d0))
2120  &vhkk(4,nhkk) = vhkk(3,nhkk)/betp-vhkk(3,nhkk-2)/bgamp
2121 *
2122 * diffractive nucleon
2123 *
2124  nhkk = nhkk + 1
2125  isthkk(nhkk) = 177
2126  idhkk(nhkk) = 88888
2127  IF (idiftp.EQ.2) jmohkk(1,nhkk) = itapoi
2128  IF (idiftp.EQ.1) jmohkk(1,nhkk) = 1
2129  jmohkk(2,nhkk) = 0
2130  jdahkk(1,nhkk) = 0
2131  jdahkk(2,nhkk) = 0
2132  phkk(1,nhkk) = pdfq1(1)
2133  phkk(2,nhkk) = pdfq1(2)
2134  phkk(3,nhkk) = pdfq1(3)
2135  phkk(4,nhkk) = pdfq1(4)
2136 C WRITE(6,*)' Diffr. Nucleon'
2137 C *,PHKK(3,NHKK),PHKK(4,NHKK),PHKK(5,NHKK)
2138  phkk(5,nhkk) = amdch3
2139  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))
2140  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))
2141  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
2142  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
2143 *
2144 *
2145 *---------- S. Roesler 21-10-93
2146 * check energy-momentum conservation
2147  IF (ipev.GE.2) THEN
2148  px = pdq1(1)+pdq2(1)+pdd1(1)+pdd2(1)+pdfq1(1)
2149  py = pdq1(2)+pdq2(2)+pdd1(2)+pdd2(2)+pdfq1(2)
2150  pz = pdq1(3)+pdq2(3)+pdd1(3)+pdd2(3)+pdfq1(3)
2151  ee = ecm-(pdq1(4)+pdq2(4)+pdd1(4)+pdd2(4)+pdfq1(4))
2152  WRITE(lout,*)'VALMSD: ENERGY-MOMENTUM-CHECK (PX,PY,PZ,E)'
2153  WRITE(lout,'(5F12.6)')px,py,pz,ee,ecm
2154  ENDIF
2155 *
2156  RETURN
2157  9999 CONTINUE
2158  ncrej = ncrej+1
2159  IF (mod(ncrej,500).EQ.0) WRITE(lout,9900) ncrej
2160  9900 FORMAT('REJECTION IN VALMSD ',i10)
2161  iirej = 0
2162  irej = 1
2163  RETURN
2164  END
2165 *
2166 *===valmdd===============================================================*
2167 *
2168  SUBROUTINE valmdd(ITAPOI,ECM,KPROJ,KTARG,IREJ)
2169 
2170 C Low mass double diffraction J.R. Feb/Mar 94
2171 
2172 **************************************************************************
2173 * Version November 1993 by Stefan Roesler *
2174 * University of Leipzig *
2175 * This subroutine selects flavors and 4-momenta of partons in low-mass *
2176 * single diffractive chains. *
2177 **************************************************************************
2178 
2179  IMPLICIT DOUBLE PRECISION (a-h,o-z)
2180  SAVE
2181  parameter(lout=6,llook=9)
2182  parameter(nmxhkk= 89998)
2183  CHARACTER*8 aname
2184  COMMON /hkkevt/ nhkk,nevhkk, isthkk(nmxhkk), idhkk(nmxhkk),
2185  & jmohkk(2,nmxhkk),jdahkk(2,nmxhkk),phkk(5,nmxhkk),
2186  & vhkk(4,nmxhkk), whkk(4,nmxhkk)
2187  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
2188  COMMON /diffra/ isingd,idiftp,ioudif,iflagd
2189  COMMON /abrdif/ xdq1,xdq2,xddq1,xddq2,
2190  & ikvq1,ikvq2,ikd1q1,ikd2q1,ikd1q2,ikd2q2,
2191  & idiffp,idifap,
2192  & amdch1,amdch2,amdch3,gamdc1,gamdc2,gamdc3,
2193  & pgxvc1,pgyvc1,pgzvc1,pgxvc2,pgyvc2,pgzvc2,
2194  & pgxvc3,pgyvc3,pgzvc3,ndch1,ndch2,ndch3,
2195  & ikdch1,ikdch2,ikdch3,
2196  & pdq1(4),pdq2(4),pdd1(4),pdd2(4),pdfq1(4)
2197  COMMON /dpar/ aname(210),am(210),ga(210),tau(210),ich(210),
2198  & ibar(210),k1(210),k2(210)
2199  common/dinpda/imps(6,6),imve(6,6),ib08(6,21),ib10(6,21),ia08(6,21)
2200  &,ia10(6,21),a1,b1,b2,b3,lt,le,bet,as,b8,ame,diq,isu
2201  COMMON /trafop/ gamp,bgamp,betp
2202  COMMON /enerin/ eproj,etarg
2203  COMMON /sdflag/ isd
2204  common/xxlmdd/ijlmdd,kdlmdd
2205  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
2206 *
2207  dimension mquark(3,30),ihkkq(-6:6),ihkkqq(-3:3,-3:3),
2208  & idx(-4:4)
2209  DATA idx /-4,-3,-1,-2,0,2,1,3,4/
2210  DATA ihkkq /-6,-5,-4,-3,-1,-2,0,2,1,3,4,5,6/
2211  DATA ihkkqq/-3301,-3103,-3203,0, 0,0,0,
2212  & -3103,-1103,-2103,0, 0, 0, 0,
2213  & -3203,-2103,-2203,0, 0, 0, 0,
2214  & 0, 0, 0,0, 0, 0, 0,
2215  & 0, 0, 0,0,2203,2103,3202,
2216  & 0, 0, 0,0,2103,1103,3103,
2217  & 0, 0, 0,0,3203,3103,3303/
2218 *
2219 *----------------------------------- quark content of hadrons:
2220 * 1, 2, 3, 4 - u, d, s, c
2221 * -1,-2,-3,-4 - au,ad,as,ac
2222 *
2223  DATA mquark/
2224  & 1,1,2, -1,-1,-2, 0,0,0, 0,0,0, 0,0,0,
2225  & 0,0,0, 0,0,0, 1,2,2, -1,-2,-2, 0,0,0,
2226  & 0,0,0, 0,0,0, 1,-2,0, 2,-1,0, 1,-3,0,
2227  & 3,-1,0, 1,2,3, -1,-2,-3, 0,0,0, 2,2,3,
2228  & 1,1,3, 1,2,3, 1,-1,0, 2,-3,0, 3,-2,0,
2229  & 2,-2,0, 3,-3,0, 0,0,0, 0,0,0, 0,0,0/
2230  DATA unon/2.0/
2231  DATA ncrej, ncpt /0, 0/
2232 *
2233  ijlmdd= 1
2234  isd = 2
2235  irej = 0
2236  iirej = 0
2237  ibproj = ibar(kproj)
2238  ibtarg = ibar(ktarg)
2239  eproj = (am(kproj)**2-am(ktarg)**2+ecm**2)/(2.0d0*ecm)
2240  etarg = (am(ktarg)**2-am(kproj)**2+ecm**2)/(2.0d0*ecm)
2241  IF(ipev.GE.2) WRITE(lout,1014)eproj,etarg
2242  1014 FORMAT('VALMSD: EPROJ,ETARG',2f10.5)
2243  IF (ibtarg.LE.0) THEN
2244  WRITE(lout,1001) ibtarg
2245  1001 FORMAT('VALMSD: NO HMSD FOR TARGET WITH BARYON-CHARGE',i4)
2246  iirej = 1
2247  goto 9999
2248  ENDIF
2249  iqp1 = mquark(1,kproj)
2250  iqp2 = mquark(2,kproj)
2251  iqp3 = mquark(3,kproj)
2252  iqt1 = mquark(1,ktarg)
2253  iqt2 = mquark(2,ktarg)
2254  iqt3 = mquark(3,ktarg)
2255  IF(ipev.GE.2) WRITE(lout,1002)
2256  & ibproj,ibtarg,iqp1,iqp2,iqp3,iqt1,iqt2,iqt3
2257  1002 FORMAT('VALMSD: IBPROJ,IBTARG,IQP1,IQP2,IQP3,IQT1,IQT2,IQT3 ',8i4)
2258 *
2259  IF (ibproj.NE.0) THEN
2260 *
2261 *-------------------- q-qq (aq-aqaq) - flavors of projectile (baryon)
2262 *
2263  isam = 1.0d0+2.999d0*rndm(v)
2264  goto(10,11,12) isam
2265  10 CONTINUE
2266  iqp = iqp1
2267  idiqp1 = iqp2
2268  idiqp2 = iqp3
2269  goto 13
2270  11 CONTINUE
2271  iqp = iqp2
2272  idiqp1 = iqp1
2273  idiqp2 = iqp3
2274  goto 13
2275  12 CONTINUE
2276  iqp = iqp3
2277  idiqp1 = iqp1
2278  idiqp2 = iqp2
2279  13 CONTINUE
2280 *
2281  ELSE IF (ibproj.EQ.0) THEN
2282 *
2283 *-------------------- q-aq - flavors of projectile (meson)
2284 *
2285  isam = 1.0d0+1.999d0*rndm(v)
2286  goto(14,15) isam
2287  14 CONTINUE
2288  iqp = iqp1
2289  idiqp1 = iqp2
2290  goto 16
2291  15 CONTINUE
2292  iqp = iqp2
2293  idiqp1 = iqp1
2294  16 CONTINUE
2295  idiqp2 = 0
2296 *
2297  ENDIF
2298 *
2299 *-------------------- q-qq - flavors of target (baryon)
2300 *
2301  isam = 1.0d0+2.999d0*rndm(v)
2302  goto(17,18,19) isam
2303  17 CONTINUE
2304  iqt = iqt1
2305  idiqt1 = iqt2
2306  idiqt2 = iqt3
2307  goto 20
2308  18 CONTINUE
2309  iqt = iqt2
2310  idiqt1 = iqt1
2311  idiqt2 = iqt3
2312  goto 20
2313  19 CONTINUE
2314  iqt = iqt3
2315  idiqt1 = iqt1
2316  idiqt2 = iqt2
2317  20 CONTINUE
2318 *
2319  ikvq1 = iqp
2320  ikd1q1 = idiqp1
2321  ikd2q1 = idiqp2
2322 *
2323  ikvq2 = iqt
2324  ikd1q2 = idiqt1
2325  ikd2q2 = idiqt2
2326 *
2327  IF (ipev.GE.2) WRITE(lout,1003)
2328  & ikvq1,ikd1q1,ikd2q1,ikvq2,ikd1q2,ikd2q2
2329  1003 FORMAT('VALMSD: IKVQ1,IKD1Q1,IKD2Q1,IKVQ2,IKD1Q2,IKD2Q2 ',6i4)
2330 *
2331 *-------------------- IDIFTP = 1 target (backward) hadron excited
2332 * IDIFTP = 2 projectile (forward) hadron excited
2333 *
2334  idiftp = 1.0d0+1.999d0*rndm(v)
2335 *-------------------- S.Roesler 5/26/93
2336  IF ((isingd.EQ.3).OR.(isingd.EQ.7)) idiftp = 1
2337  IF ((isingd.EQ.4).OR.(isingd.EQ.8)) idiftp = 2
2338 *
2339  IF (ipev.GE.2) WRITE(lout,1004) idiftp
2340  1004 FORMAT('VALMSD: IDIFTP ',i4)
2341  IF ((idiftp.NE.1).AND.(idiftp.NE.2)) THEN
2342  IF (ipev.GE.2) WRITE(lout,'(A19)') 'VALMSD-ERROR: IDIFTP'
2343  goto 9999
2344  ENDIF
2345 *
2346 *-------------------- diffractive mass
2347 *
2348 31 CONTINUE
2349  amo = 1.5d0
2350 C------------------------------------------------------------
2351 C
2352 C first simple test j.r. 2/94
2353 C
2354 C-----------------------------------------------------------
2355 C AMO = 4.0D0
2356 C-----------------------------------------------------------
2357  r = rndm(v)
2358  IF ((ibproj.EQ.0).AND.(idiftp.EQ.1))
2359  & amo = 1.5d0*r+2.83d0*(1.0d0-r)
2360  IF ((ibproj.EQ.0).AND.(idiftp.EQ.2)) amo = 1.0d0
2361  sam = 1.0d0
2362  IF (ecm.LE.300.0d0) sam = 1.0d0-exp(-((ecm/200.0d0)**4))
2363  r = rndm(v)*sam
2364  amax= (1.0d0-sam)*sqrt(0.1d0*ecm**2)+sam*sqrt(400.0d0)
2365  amu = r*sqrt(100.0d0)+(1.0d0-r)*amax
2366  IF (ibproj.EQ.0) THEN
2367 C----------------- change mass-cuts
2368 C AMAX= (1.0D0-SAM)*SQRT(0.25D0*ECM**2)+SAM*SQRT(400.0D0)
2369  amax= (1.0d0-sam)*sqrt(0.15d0*ecm**2)+sam*sqrt(400.0d0)
2370  amu = r*sqrt(100.0d0)+(1.0d0-r)*amax
2371  ENDIF
2372 C---------------------------------------------------------------
2373 C
2374 C j.r. test 2/94
2375 C
2376 C---------------------------------------------------------------
2377  amu=2.d0*amu
2378 C---------------------------------------------------------------
2379  r = rndm(v)
2380  IF (ecm.LE.50.0d0) THEN
2381  amdiff = amo*(amu/amo)**r
2382  ELSE
2383  a = 0.7d0
2384  IF (ecm.LE.300.0d0) a = 0.7d0*(1.0d0-exp(-((ecm/100.0d0)**2)))
2385  amdiff = 1.0d0/((r/(amu**a)+(1.0d0-r)/(amo**a))
2386  & **(1.0d0/a))
2387  ENDIF
2388  IF(amdiff.GT.0.5d0*ecm)go to 31
2389  IF (iouxev.GE.2) WRITE(lout,1005) amdiff
2390  1005 FORMAT('VALMSD: AMDIFF',e10.5)
2391 *
2392 *-------------------- kinematical parameters of
2393 * diffractive projectile (IDIFTP = 1)
2394 * or diffractive target (IDIFTP = 2)
2395 *
2396  IF (idiftp.EQ.1) THEN
2397  amdch3 = am(kproj)
2398  IF(kproj.EQ.1)THEN
2399  kdlmdd = 61
2400  amdch3 = am(61)
2401 C 19.4.99
2402  kdlmdd = 54
2403  amdch3 = am(54)
2404  ELSEIF(kproj.EQ.2)THEN
2405 C 19.4.99
2406  kdlmdd = 68
2407  amdch3 = am(68)
2408  ELSEIF(kproj.EQ.8)THEN
2409  amdch3 = am(62)
2410  kdlmdd = 62
2411 C 19.4.99
2412  amdch3 = am(55)
2413  kdlmdd = 55
2414  ELSEIF(kproj.EQ.13)THEN
2415  amdch3 = am(186)
2416  kdlmdd = 186
2417  ELSEIF(kproj.EQ.14)THEN
2418  amdch3 = am(188)
2419  kdlmdd = 188
2420  ELSEIF(kproj.EQ.15)THEN
2421  amdch3 = am(190)
2422  kdlmdd = 190
2423  ELSEIF(kproj.EQ.16)THEN
2424  amdch3 = am(191)
2425  kdlmdd = 191
2426  ELSEIF(kproj.EQ.23)THEN
2427  amdch3 = am(187)
2428  kdlmdd = 187
2429  ELSEIF(kproj.EQ.24)THEN
2430  amdch3 = am(192)
2431  kdlmdd = 192
2432  ELSEIF(kproj.EQ.25)THEN
2433  amdch3 = am(193)
2434  kdlmdd = 193
2435  ENDIF
2436  ELSE IF (idiftp.EQ.2) THEN
2437  amdch3 = am(ktarg)
2438  IF(ktarg.EQ.1)THEN
2439  kdlmdd = 61
2440  amdch3 = am(61)
2441 C 19.4.99
2442  kdlmdd = 54
2443  amdch3 = am(54)
2444  ELSEIF(ktarg.EQ.2)THEN
2445 C 19.4.99
2446  kdlmdd = 68
2447  amdch3 = am(68)
2448  ELSEIF(ktarg.EQ.8)THEN
2449  amdch3 = am(62)
2450  kdlmdd = 62
2451 C 19.4.99
2452  amdch3 = am(55)
2453  kdlmdd = 55
2454  ELSEIF(ktarg.EQ.13)THEN
2455  amdch3 = am(186)
2456  kdlmdd = 186
2457  ELSEIF(ktarg.EQ.14)THEN
2458  amdch3 = am(188)
2459  kdlmdd = 188
2460  ELSEIF(ktarg.EQ.15)THEN
2461  amdch3 = am(190)
2462  kdlmdd = 190
2463  ELSEIF(ktarg.EQ.16)THEN
2464  amdch3 = am(191)
2465  kdlmdd = 191
2466  ELSEIF(ktarg.EQ.23)THEN
2467  amdch3 = am(187)
2468  kdlmdd = 187
2469  ELSEIF(ktarg.EQ.24)THEN
2470  amdch3 = am(192)
2471  kdlmdd = 192
2472  ELSEIF(ktarg.EQ.25)THEN
2473  amdch3 = am(193)
2474  kdlmdd = 193
2475  ENDIF
2476  ENDIF
2477  ech3 = (ecm**2+amdch3**2-amdiff**2)/(2.0d0*ecm)
2478  IF (ech3.LE.amdch3) THEN
2479  ncmerr = ncmerr+1
2480  IF (mod(ncmerr,2200).EQ.0)
2481  & WRITE(lout,*) 'LMSD: INEFFICIENT SELECTION OF AMDIFF(1),
2482  & NCMERR = ',ncmerr
2483  goto 31
2484  ENDIF
2485 *
2486  pch3 = sqrt(abs(ech3-amdch3))*sqrt(ech3+amdch3)
2487  IF (idiftp.EQ.2) pch3 = -pch3
2488  ndch3 = -1
2489  gamdc3 = ech3/amdch3
2490  pgvc3 = pch3/amdch3
2491  IF (iouxev.GE.2) WRITE(lout,1006) amdch3,ech3,pch3,gamdc3,pgvc3
2492  1006 FORMAT('VALMSD: AMDCH3,ECH3,PCH3,GAMDC3,PGVC3',5e10.5)
2493 *
2494  30 CONTINUE
2495  b33 = 8.0d0
2496  es = -2.0d0/(b33**2)*log(rndm(v)*rndm(v))
2497  hps = sqrt(es*es+2.0d0*es*0.94d0)
2498  CALL dsfecf(sfe,cfe)
2499  pxch3 = hps*cfe
2500  pych3 = hps*sfe
2501  ptch3 = sqrt(pxch3**2+pych3**2)
2502  IF (ptch3.GT.abs(pch3)) THEN
2503  ncpt = ncpt+1
2504  IF (mod(ncpt,500).EQ.0) WRITE(lout,1007) ncpt
2505  1007 FORMAT('VALMSD: INEFFICIENT PT-SELECTION 1, NCPT=',i8)
2506  goto 30
2507  ENDIF
2508  plch3 = sign(sqrt(abs(pch3-ptch3))*sqrt(abs(pch3+ptch3)),pch3)
2509  IF (iouxev.GE.2) WRITE(lout,1008) es,hps,pxch3,pych3,plch3
2510  1008 FORMAT('VALMSD: ES,HPS,PXCH3,PYCH3,PLCH3',5e10.5)
2511 *
2512 *-------------------- no chain 1
2513 *
2514  ndch1 = -99
2515  amdch1 = 0.0d0
2516  ech1 = 0.0d0
2517  pch1 = 0.0d0
2518  gamdc1 = 1.0d0
2519  pgvc1 = 0.0d0
2520  pgxvc1 = 0.0d0
2521  pgyvc1 = 0.0d0
2522  pgzvc1 = 0.0d0
2523  DO 40 i=1,4
2524  pdd2(i) = 0.0d0
2525  pdq1(i) = 0.0d0
2526  40 CONTINUE
2527 *
2528 *-------------------- kinematical parameters of
2529 * excited target (IDIFTP = 1)
2530 * or excited projectile (IDIFTP = 2)
2531 *
2532  amdch2 = amdiff
2533  ech2 = (ecm**2+amdch2**2-amdch3**2)/(2.0d0*ecm)
2534  pch2 = -sign(sqrt(abs(ech2-amdch2))*sqrt(ech2+amdch2),pch3)
2535  ndch2 = 0
2536  gamdc2 = ech2/amdch2
2537  pgvc2 = pch2/amdch2
2538 *
2539 *-------------------- (IDIFFP,IDIFAP in analogy to HMSD in VAHMSD)
2540 * IDIFTP = 1 : IKD1Q2/IKD2Q2 - IDIFFP = IKVQ2
2541 * IDIFTP = 2 : projectile - baryon
2542 * IKD1Q1/IKD2Q1 - IDIFFP = IKVQ1
2543 * projectile - antibaryon
2544 * IKD1Q1/IKD2Q1 - IDIFAP = IKVQ1
2545 * projectile - meson
2546 * IKD1Q1 > 0 - IDIFAP = IKVQ1
2547 * IKD1Q1 < 0 - IDIFFP = IKVQ1
2548 *
2549  IF (idiftp.EQ.1) idiffp = ikvq2
2550  IF (idiftp.EQ.2) THEN
2551  IF (ibproj.GT.0) idiffp = ikvq1
2552  IF (ibproj.LT.0) idifap = ikvq1
2553  IF ((ibproj.EQ.0).AND.(ikd1q1.GT.0)) idifap = ikvq1
2554  IF ((ibproj.EQ.0).AND.(ikd1q1.LT.0)) idiffp = ikvq1
2555  ENDIF
2556  IF (iouxev.GE.2) WRITE(lout,1009) amdch2,ech2,pch2,gamdc2,pgvc2,
2557  & idiffp,idifap
2558  1009 FORMAT('VALMSD: AMDCH2,ECH2,PCH2,GAMDC2,PGVC2,IDIFFP,IDIFAP',
2559  & 5e10.5,2i4)
2560  pxch2 = -pxch3
2561  pych2 = -pych3
2562  ptch2 = sqrt(pxch2**2+pych2**2)
2563  IF (ptch2.GT.abs(pch2)) THEN
2564  ncpt = ncpt+1
2565  IF (mod(ncpt,500).EQ.0) WRITE(lout,1010) ncpt
2566  1010 FORMAT('VALMSD: INEFFICIENT PT-SELECTION 2, NCPT=',i8)
2567  goto 30
2568  ENDIF
2569  plch2 = sign(sqrt(abs(pch2-ptch2))*sqrt(abs(pch2+ptch2)),pch2)
2570  IF (iouxev.GE.2) WRITE(lout,1011) pxch2,pych2,plch2
2571  1011 FORMAT('VALMSD: PXCH2,PYCH2,PLCH2',3e10.5)
2572  cc = amdch2/(2.0d0*ech2)
2573 *--------------------------- S.R. 11/12/92
2574  IF (cc.GE.0.5d0) THEN
2575  ncmerr = ncmerr+1
2576  IF (mod(ncmerr,200).EQ.0)
2577  & WRITE(lout,*) 'LMSD: INEFFICIENT SELECTION OF AMDIFF(2),
2578  & NCMERR = ',ncmerr
2579  goto 31
2580  ENDIF
2581 *
2582  xdiq = 0.5d0+sqrt(abs((0.5d0-cc)*(0.5d0+cc)))
2583  xq = 1.0d0-xdiq
2584  ediq = ech2*xdiq
2585  eq = ech2*xq
2586  temp = abs(ediq/pch2)
2587  pxdiq = pxch2*temp
2588  pydiq = pych2*temp
2589  pldiq = plch2*temp
2590  temp = abs(eq/pch2)
2591  pxq = -pxch2*temp
2592  pyq = -pych2*temp
2593  plq = -plch2*temp
2594 *
2595 *-------------------- store results in COMMON-block /ABRDIF/
2596 *
2597  pdq2(1) = pxq
2598  pdq2(2) = pyq
2599  pdq2(3) = plq
2600  pdq2(4) = eq
2601  pdd1(1) = pxdiq
2602  pdd1(2) = pydiq
2603  pdd1(3) = pldiq
2604  pdd1(4) = ediq
2605  pdfq1(1) = pxch3
2606  pdfq1(2) = pych3
2607  pdfq1(3) = plch3
2608  pdfq1(4) = ech3
2609  pgxvc2 = pxch2/amdch2
2610  pgyvc2 = pych2/amdch2
2611  pgzvc2 = plch2/amdch2
2612  pgxvc3 = pxch3/amdch3
2613  pgyvc3 = pych3/amdch3
2614  pgzvc3 = plch3/amdch3
2615 *
2616  IF (ipev.GE.2) THEN
2617  WRITE(lout,*) 'VALMSD: PDQ1, PDD2, PDD1, PDQ2, PDFQ1'
2618  WRITE(lout,1012)
2619  & (pdq1(i),pdd2(i),pdd1(i),pdq2(i),pdfq1(i),i=1,4)
2620  1012 FORMAT(5f10.5)
2621  WRITE(lout,1013) 'PGXVC1,PGYVC1,PGZVC1,GAMDC1',
2622  & pgxvc1,pgyvc1,pgzvc1,gamdc1
2623  WRITE(lout,1013) 'PGXVC2,PGYVC2,PGZVC2,GAMDC2',
2624  & pgxvc2,pgyvc2,pgzvc2,gamdc2
2625  WRITE(lout,1013) 'PGXVC3,PGYVC3,PGZVC3,GAMDC3',
2626  & pgxvc3,pgyvc3,pgzvc3,gamdc3
2627  1013 FORMAT(a27,4f10.5)
2628  ENDIF
2629 *
2630 *------------------- store results in common block /HKKEVT/
2631 *
2632 * partonen of projectile/target
2633 *
2634  nhkk = nhkk + 1
2635  isthkk(nhkk) = 23-idiftp
2636  IF (idiftp.EQ.1) idhkk(nhkk) = ihkkq(idx(ikvq2))
2637  IF (idiftp.EQ.2) idhkk(nhkk) = ihkkq(idx(ikvq1))
2638  IF (idiftp.EQ.1) jmohkk(1,nhkk) = itapoi
2639  IF (idiftp.EQ.2) jmohkk(1,nhkk) = 1
2640  jmohkk(2,nhkk) = 0
2641  jdahkk(1,nhkk) = 0
2642  jdahkk(2,nhkk) = 0
2643  phkk(1,nhkk) = 0.0d0
2644  phkk(2,nhkk) = 0.0d0
2645  phkk(3,nhkk) = xq
2646  phkk(4,nhkk) = xq
2647  phkk(5,nhkk) = 0.0d0
2648  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))
2649  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))
2650  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
2651  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
2652 *
2653  nhkk = nhkk + 1
2654  isthkk(nhkk) = 23-idiftp
2655  IF (idiftp.EQ.1) idhkk(nhkk) = ihkkqq(idx(ikd1q2),idx(ikd2q2))
2656  IF (idiftp.EQ.2) idhkk(nhkk) = ihkkqq(idx(ikd1q1),idx(ikd2q1))
2657  IF (idiftp.EQ.1) jmohkk(1,nhkk) = itapoi
2658  IF (idiftp.EQ.2) jmohkk(1,nhkk) = 1
2659  jmohkk(2,nhkk) = 0
2660  jdahkk(1,nhkk) = 0
2661  jdahkk(2,nhkk) = 0
2662  phkk(1,nhkk) = 0.0d0
2663  phkk(2,nhkk) = 0.0d0
2664  phkk(3,nhkk) = xdiq
2665  phkk(4,nhkk) = xdiq
2666  phkk(5,nhkk) = 0.0d0
2667  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))
2668  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))
2669  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
2670  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
2671 *
2672 * ends of chain 2
2673 *
2674  nhkk = nhkk + 1
2675  isthkk(nhkk) = 100+isthkk(nhkk-2)
2676  idhkk(nhkk) = idhkk(nhkk-2)
2677  jmohkk(1,nhkk) = nhkk-2
2678  jmohkk(2,nhkk) = jmohkk(1,nhkk-2)
2679  jdahkk(1,nhkk) = nhkk+2
2680  jdahkk(2,nhkk) = nhkk+2
2681  phkk(1,nhkk) = pdq2(1)
2682  phkk(2,nhkk) = pdq2(2)
2683  phkk(3,nhkk) = pdq2(3)
2684  phkk(4,nhkk) = pdq2(4)
2685  phkk(5,nhkk) = 0.0d0
2686 C Add position of parton in hadron
2687  CALL qinnuc(xxpp,yypp)
2688  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))+xxpp
2689  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))+yypp
2690  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
2691  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
2692 *
2693  nhkk = nhkk + 1
2694  isthkk(nhkk) = 100+isthkk(nhkk-2)
2695  idhkk(nhkk) = idhkk(nhkk-2)
2696  jmohkk(1,nhkk) = nhkk-2
2697  jmohkk(2,nhkk) = jmohkk(1,nhkk-2)
2698  jdahkk(1,nhkk) = nhkk+1
2699  jdahkk(2,nhkk) = nhkk+1
2700  phkk(1,nhkk) = pdd1(1)
2701  phkk(2,nhkk) = pdd1(2)
2702  phkk(3,nhkk) = pdd1(3)
2703  phkk(4,nhkk) = pdd1(4)
2704  phkk(5,nhkk) = 0.0d0
2705 C Add position of parton in hadron
2706  CALL qinnuc(xxpp,yypp)
2707  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))+xxpp
2708  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))+yypp
2709  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
2710  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
2711 *
2712 * chain 2
2713 *
2714  nhkk = nhkk + 1
2715  isthkk(nhkk) = 188
2716  idhkk(nhkk) = 88888
2717  jmohkk(1,nhkk) = nhkk-2
2718  jmohkk(2,nhkk) = nhkk-1
2719  jdahkk(1,nhkk) = 0
2720  jdahkk(2,nhkk) = 0
2721  phkk(1,nhkk) = pxch2
2722  phkk(2,nhkk) = pych2
2723  phkk(3,nhkk) = plch2
2724  phkk(4,nhkk) = ech2
2725  phkk(5,nhkk) = amdch2
2726  vhkk(1,nhkk) = vhkk(1,nhkk-1)
2727  vhkk(2,nhkk) = vhkk(2,nhkk-1)
2728  vhkk(3,nhkk) = vhkk(3,nhkk-1)
2729  IF ((betp.NE.0.0d0).AND.(bgamp.NE.0.0d0))
2730  &vhkk(4,nhkk) = vhkk(3,nhkk)/betp-vhkk(3,nhkk-2)/bgamp
2731 *
2732 * diffractive nucleon
2733 *
2734  nhkk = nhkk + 1
2735  isthkk(nhkk) = 188
2736  idhkk(nhkk) = 88888
2737  IF (idiftp.EQ.2) jmohkk(1,nhkk) = itapoi
2738  IF (idiftp.EQ.1) jmohkk(1,nhkk) = 1
2739  jmohkk(2,nhkk) = 0
2740  jdahkk(1,nhkk) = 0
2741  jdahkk(2,nhkk) = 0
2742  phkk(1,nhkk) = pdfq1(1)
2743  phkk(2,nhkk) = pdfq1(2)
2744  phkk(3,nhkk) = pdfq1(3)
2745  phkk(4,nhkk) = pdfq1(4)
2746  phkk(5,nhkk) = amdch3
2747  vhkk(1,nhkk) = vhkk(1,jmohkk(1,nhkk))
2748  vhkk(2,nhkk) = vhkk(2,jmohkk(1,nhkk))
2749  vhkk(3,nhkk) = vhkk(3,jmohkk(1,nhkk))
2750  vhkk(4,nhkk) = vhkk(4,jmohkk(1,nhkk))
2751 *
2752 *
2753 *---------- S. Roesler 21-10-93
2754 * check energy-momentum conservation
2755  IF (ipev.GE.2) THEN
2756  px = pdq1(1)+pdq2(1)+pdd1(1)+pdd2(1)+pdfq1(1)
2757  py = pdq1(2)+pdq2(2)+pdd1(2)+pdd2(2)+pdfq1(2)
2758  pz = pdq1(3)+pdq2(3)+pdd1(3)+pdd2(3)+pdfq1(3)
2759  ee = ecm-(pdq1(4)+pdq2(4)+pdd1(4)+pdd2(4)+pdfq1(4))
2760  WRITE(lout,*)'VALMSD: ENERGY-MOMENTUM-CHECK (PX,PY,PZ,E)'
2761  WRITE(lout,'(5F12.6)')px,py,pz,ee,ecm
2762  ENDIF
2763 *
2764  RETURN
2765  9999 CONTINUE
2766  ncrej = ncrej+1
2767  IF (mod(ncrej,500).EQ.0) WRITE(lout,9900) ncrej
2768  9900 FORMAT('REJECTION IN VALMSD ',i10)
2769  iirej = 0
2770  irej = 1
2771  RETURN
2772  END
2773 *
2774 *===hadrdi===============================================================*
2775 *
2776  SUBROUTINE hadrdi(NAUX,IJPROJ,IJTAR,NHKKH1)
2777 
2778 **************************************************************************
2779 * Version November 1993 by Stefan Roesler *
2780 * University of Leipzig *
2781 * This subroutine prepares the hadronisation of diffractive chains *
2782 * (single-diffractive component), calls HADJET and stores the results *
2783 * in common block /HKKEVT/. *
2784 **************************************************************************
2785 
2786  IMPLICIT DOUBLE PRECISION (a-h,o-z)
2787  SAVE
2788  parameter(lout=6,llook=9)
2789  parameter(nmxhkk= 89998,intmx=2488,naumax=897,nfimax=249)
2790  CHARACTER*8 aname,anc,anf
2791  COMMON /hkkevt/ nhkk,nevhkk, isthkk(nmxhkk), idhkk(nmxhkk),
2792  & jmohkk(2,nmxhkk),jdahkk(2,nmxhkk),phkk(5,nmxhkk),
2793  & vhkk(4,nmxhkk), whkk(4,nmxhkk)
2794  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
2795  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
2796  & iibar(210),k1(210),k2(210)
2797  COMMON /diffra/ isingd,idiftp,ioudif,iflagd
2798  COMMON /difpar/ pxc(902), pyc(902),pzc(902),
2799  & hec(902), amc(902),ichc(902),
2800  & ibarc(902),anc(902),nrc(902)
2801  COMMON /dfinpa/ anf(nfimax),pxf(nfimax),pyf(nfimax),
2802  & pzf(nfimax),hef(nfimax),amf(nfimax),
2803  & ichf(nfimax),ibarf(nfimax),nref(nfimax)
2804  COMMON /dfinpz/iormo(nfimax),idaug1(nfimax),idaug2(nfimax),
2805  * istath(nfimax)
2806  COMMON /abrdif/ xdq1,xdq2,xddq1,xddq2,
2807  & ikvq1,ikvq2,ikd1q1,ikd2q1,ikd1q2,ikd2q2,
2808  & idiffp,idifap,
2809  & amdch1,amdch2,amdch3,gamdc1,gamdc2,gamdc3,
2810  & pgxvc1,pgyvc1,pgzvc1,pgxvc2,pgyvc2,pgzvc2,
2811  & pgxvc3,pgyvc3,pgzvc3,ndch1,ndch2,ndch3,
2812  & ikdch1,ikdch2,ikdch3,
2813  & pdq1(4),pdq2(4),pdd1(4),pdd2(4),pdfq1(4)
2814  COMMON /dinpda/ imps(6,6),imve(6,6),ib08(6,21),ib10(6,21),
2815  & ia08(6,21),ia10(6,21),
2816  & a1,b1,b2,b3,lt,le,bet,as,b8,ame,diq,isu
2817  COMMON /enerin/ eproj,etarg
2818  COMMON /difout/ amchdi,ttt,nnaux,kproj,ktarg
2819  COMMON /sdflag/ isd
2820  common/xxlmdd/ijlmdd,kdlmdd
2821 C modified DPMJET
2822  COMMON /bufueh/ annvv,annss,annsv,annvs,anncc,
2823  * anndv,annvd,annds,annsd,
2824  * annhh,annzz,
2825  * ptvv,ptss,ptsv,ptvs,ptcc,ptdv,ptvd,ptds,ptsd,
2826  * pthh,ptzz,
2827  * eevv,eess,eesv,eevs,eecc,eedv,eevd,eeds,eesd,
2828  * eehh,eezz
2829  * ,anndi,ptdi,eedi
2830  * ,annzd,anndz,ptzd,ptdz,eezd,eedz
2831 C---------------------
2832 *
2833  dimension poj1(4),pat1(4),poj2(4),pat2(4)
2834 *
2835  nhkkh1 = nhkk
2836  naux = 0
2837  nhad = 0
2838  pxdi = 0.0d0
2839  pydi = 0.0d0
2840  pldi = 0.0d0
2841  edi = 0.0d0
2842  imoh2 = nhkk-1
2843  imoh3 = nhkk
2844  IF (isd.EQ.1) imoh1 = nhkk-4
2845  ibproj = iibar(ijproj)
2846  ibtar = iibar(ijtar)
2847  kproj = ijproj
2848  ktarg = ijtar
2849 *
2850  IF (ipev.GE.2) WRITE(lout,1100) isd
2851  1100 FORMAT('HADRDI: ISD ',i3)
2852 *
2853  nunuc1 = 3
2854  nunuc2 = 3
2855  IF ((ibproj.NE.0).AND.(idiftp.EQ.2)) nunuc2 = 6
2856  IF ((ibtar .NE.0).AND.(idiftp.EQ.1)) nunuc2 = 4
2857  IF (idiftp.EQ.2) THEN
2858  DO 10 j=1,4
2859  poj1(j) = pdq1(j)
2860  pat1(j) = pdd2(j)
2861  poj2(j) = pdd1(j)
2862  pat2(j) = pdq2(j)
2863  10 CONTINUE
2864  ELSE
2865  DO 11 j=1,4
2866  poj1(j) = pdd2(j)
2867  pat1(j) = pdq1(j)
2868  poj2(j) = pdq2(j)
2869  pat2(j) = pdd1(j)
2870  11 CONTINUE
2871  ENDIF
2872  IF (idiftp.EQ.1) THEN
2873  ifb11 = idifap
2874  ifb12 = ikvq2
2875  ifb21 = idiffp
2876  ifb22 = ikd1q2
2877  ifb23 = ikd2q2
2878  ELSE IF ((idiftp.EQ.2).AND.(ikvq1.GT.0)) THEN
2879  ifb11 = ikvq1
2880  ifb12 = idifap
2881  ifb21 = ikd1q1
2882  ifb22 = ikd2q1
2883  ifb23 = idiffp
2884  ELSE IF ((idiftp.EQ.2).AND.(ikvq1.LE.0)) THEN
2885  ifb11 = ikvq1
2886  ifb12 = idiffp
2887  ifb21 = ikd1q1
2888  ifb22 = ikd2q1
2889  ifb23 = idifap
2890  ENDIF
2891  IF ((ifb22.EQ.0).AND.(ifb23.NE.0)) THEN
2892  ifb22 = ifb23
2893  ifb23 = 0
2894  ENDIF
2895  IF (ifb11.LT.0) ifb11 = iabs(ifb11)+6
2896  IF (ifb12.LT.0) ifb12 = iabs(ifb12)+6
2897  IF (ifb21.LT.0) ifb21 = iabs(ifb21)+6
2898  IF (ifb22.LT.0) ifb22 = iabs(ifb22)+6
2899  IF (ifb23.LT.0) ifb23 = iabs(ifb23)+6
2900 *
2901 *-------------------------- chain 1
2902 *
2903 C WRITE(LOUT,1001) (POJ1(I),I=1,4)
2904 C WRITE(LOUT,1002) (PAT1(I),I=1,4)
2905  IF (ipev.GE.2) THEN
2906  WRITE(lout,1000) nhad,ikdch1,nunuc1,ndch1,ifb11,ifb12,
2907  & ifb13,ifb14
2908  WRITE(lout,1001) (poj1(i),i=1,4)
2909  WRITE(lout,1002) (pat1(i),i=1,4)
2910  WRITE(lout,1003) gamdc1,pgxvc1,pgyvc1,pgzvc1,amdch1
2911  1000 FORMAT('HADRDI: NHAD,IKDCH1,NUNUC1,NDCH1,IFB11,IFB12',
2912  & ', IFB13,IFB14 ',8i4)
2913  1001 FORMAT('HADRDI: POJ1 ',4e15.5)
2914  1002 FORMAT('HADRDI: PAT1 ',4e15.5)
2915  1003 FORMAT('HADRDI: GAMDC1,PGXVC1,PGYVC1,PGZVC1,AMDCH1',5f10.5)
2916  ENDIF
2917  CALL hadjet(nhad,amdch1,poj1,pat1,gamdc1,pgxvc1,
2918  & pgyvc1,pgzvc1,ifb11,ifb12,ifb13,ifb14,
2919  & ikdch1,ikdch1,nunuc1,ndch1,7)
2920 *
2921  nhkkbe = nhkk+1
2922  IF (nhad.EQ.0) goto 13
2923  DO 12 j=1,nhad
2924  IF (nhkk.EQ.nmxhkk) THEN
2925  WRITE(lout,1004) nhkk
2926  1004 FORMAT(.EQ.'HADRDI: NHKKNMXHKK',i10)
2927  goto 9999
2928  ENDIF
2929  nhkk = nhkk+1
2930  naux = naux+1
2931  IF(ibarf(j).EQ.500)go to 776
2932  echeck = sqrt(pxf(j)**2+pyf(j)**2+pzf(j)**2+amf(j)**2)
2933  IF (abs(echeck-hef(j)).GT.0.5d-2) THEN
2934  WRITE(lout,1005) nhkk,echeck,hef(j),amf(j)
2935  1005 FORMAT('HADRDI: CHAIN 1 CORRECT INCONSISTENT ENERGY',
2936  & i10,3f10.5)
2937  hef(j) = echeck
2938  ENDIF
2939  776 CONTINUE
2940  anndi=anndi+1
2941  eedi=eedi+hef(j)
2942  ptdi=ptdi+sqrt(pxf(j)**2+pyf(j)**2)
2943  pxc(naux) = pxf(j)
2944  pyc(naux) = pyf(j)
2945  pzc(naux) = pzf(j)
2946  hec(naux) = hef(j)
2947  amc(naux) = amf(j)
2948  ichc(naux) = ichf(j)
2949  ibarc(naux) = ibarf(j)
2950  anc(naux) = anf(j)
2951  nrc(naux) = nref(j)
2952  pxdi = pxdi+pxf(j)
2953  pydi = pydi+pyf(j)
2954  pldi = pldi+pzf(j)
2955  edi = edi +hef(j)
2956  ediff = ediff+hef(j)
2957  ptdiff = ptdiff+sqrt(pxf(j)**2+pyf(j)**2)
2958  isthkk(nhkk) = 1
2959  IF(ibarf(j).EQ.500)isthkk(nhkk)=2
2960  idhkk(nhkk) = mpdgha(nref(j))
2961  IF(iormo(j).EQ.999)THEN
2962  jmohkk(1,nhkk) = imoh1
2963  ELSE
2964  jmohkk(1,nhkk)=nhkkbe+iormo(j)-1
2965  ENDIF
2966  jmohkk(2,nhkk) = 0
2967  jdahkk(1,nhkk) = 0
2968  jdahkk(2,nhkk) = 0
2969  phkk(1,nhkk) = pxf(j)
2970  phkk(2,nhkk) = pyf(j)
2971  phkk(3,nhkk) = pzf(j)
2972  phkk(4,nhkk) = hef(j)
2973  phkk(5,nhkk) = amf(j)
2974  imohkk = jmohkk(1,nhkk)
2975  vhkk(1,nhkk) = vhkk(1,imohkk)
2976  vhkk(2,nhkk) = vhkk(2,imohkk)
2977  vhkk(3,nhkk) = vhkk(3,imohkk)
2978  vhkk(4,nhkk) = vhkk(4,imohkk)
2979  IF (iphkk.GE.1) THEN
2980  WRITE(lout,1006) nhkk,isthkk(nhkk),idhkk(nhkk),
2981  & jmohkk(1,nhkk),jmohkk(2,nhkk),
2982  & jdahkk(1,nhkk),jdahkk(2,nhkk)
2983  WRITE(lout,1007) (phkk(k,nhkk),k=1,5)
2984  1006 FORMAT('HADRDI: NHKK,ISTHKK,IDHKK,JMOHKK,JDAHKK',7i6)
2985  1007 FORMAT('HADRDI: PHKK ',5f10.5)
2986  ENDIF
2987  12 CONTINUE
2988  13 CONTINUE
2989  IF ((nhad.GT.0).AND.(isd.EQ.1)) THEN
2990  jdahkk(1,imoh1) = nhkkbe
2991  jdahkk(2,imoh1) = nhkk
2992  ENDIF
2993 *
2994 *------------------------ chain 2
2995 *
2996  nhad = 0
2997 C WRITE(LOUT,1009) (POJ2(I),I=1,4)
2998 C WRITE(LOUT,1010) (PAT2(I),I=1,4)
2999  IF (ipev.GE.2) THEN
3000  WRITE(lout,1008) nhad,ikdch2,nunuc2,ndch2,ifb21,ifb22,
3001  & ifb23,ifb24
3002  WRITE(lout,1009) (poj2(i),i=1,4)
3003  WRITE(lout,1010) (pat2(i),i=1,4)
3004  WRITE(lout,1011) gamdc2,pgxvc2,pgyvc2,pgzvc2,amdch2
3005  1008 FORMAT('HADRDI: NHAD,IKDCH2,NUNUC2,NDCH2,IFB21,IFB22,
3006  & IFB23,IFB24 ',8i4)
3007  1009 FORMAT('HADRDI: POJ2 ',4e15.5)
3008  1010 FORMAT('HADRDI: PAT2 ',4e15.5)
3009  1011 FORMAT('HADRDI: GAMDC2,PGXVC2,PGYVC2,PGZVC2,AMDCH2',5f10.5)
3010  ENDIF
3011  CALL hadjet(nhad,amdch2,poj2,pat2,gamdc2,pgxvc2,
3012  & pgyvc2,pgzvc2,ifb21,ifb22,ifb23,ifb24,
3013  & ikdch2,ikdch2,nunuc2,ndch2,8)
3014 *
3015  nhkkbe = nhkk+1
3016  IF (nhad.EQ.0) goto 15
3017  DO 14 j=1,nhad
3018  IF (nhkk.EQ.nmxhkk) THEN
3019  WRITE(lout,1012) nhkk
3020  1012 FORMAT(.EQ.'HADRDI: NHKKNMXHKK',i10)
3021  goto 9999
3022  ENDIF
3023  nhkk = nhkk+1
3024  naux = naux+1
3025  IF(ibarf(j).EQ.500)go to 775
3026  echeck = sqrt(pxf(j)**2+pyf(j)**2+pzf(j)**2+amf(j)**2)
3027  IF (abs(echeck-hef(j)).GT.0.5d-2) THEN
3028  WRITE(lout,1013) nhkk,echeck,hef(j),amf(j)
3029  1013 FORMAT('HADRDI: CHAIN 2 CORRECT INCONSISTENT ENERGY',
3030  & i10,3f10.5)
3031  hef(j) = echeck
3032  ENDIF
3033  775 CONTINUE
3034  anndi=anndi+1
3035  eedi=eedi+hef(j)
3036  ptdi=ptdi+sqrt(pxf(j)**2+pyf(j)**2)
3037  pxc(naux) = pxf(j)
3038  pyc(naux) = pyf(j)
3039  pzc(naux) = pzf(j)
3040  hec(naux) = hef(j)
3041  amc(naux) = amf(j)
3042  ichc(naux) = ichf(j)
3043  ibarc(naux) = ibarf(j)
3044  anc(naux) = anf(j)
3045  nrc(naux) = nref(j)
3046  pxdi = pxdi+pxf(j)
3047  pydi = pydi+pyf(j)
3048  pldi = pldi+pzf(j)
3049  edi = edi +hef(j)
3050  ediff = ediff+hef(j)
3051  ptdiff = ptdiff+sqrt(pxf(j)**2+pyf(j)**2)
3052  isthkk(nhkk) = 1
3053  IF(ibarf(j).EQ.500)isthkk(nhkk)=2
3054  idhkk(nhkk) = mpdgha(nref(j))
3055  IF(iormo(j).EQ.999)THEN
3056  jmohkk(1,nhkk) = imoh2
3057  ELSE
3058  jmohkk(1,nhkk)=nhkkbe+iormo(j)-1
3059  ENDIF
3060  jmohkk(2,nhkk) = 0
3061  jdahkk(1,nhkk) = 0
3062  jdahkk(2,nhkk) = 0
3063  phkk(1,nhkk) = pxf(j)
3064  phkk(2,nhkk) = pyf(j)
3065  phkk(3,nhkk) = pzf(j)
3066  phkk(4,nhkk) = hef(j)
3067  phkk(5,nhkk) = amf(j)
3068  imohkk = jmohkk(1,nhkk)
3069  vhkk(1,nhkk) = vhkk(1,imohkk)
3070  vhkk(2,nhkk) = vhkk(2,imohkk)
3071  vhkk(3,nhkk) = vhkk(3,imohkk)
3072  vhkk(4,nhkk) = vhkk(4,imohkk)
3073  IF (iphkk.GE.2) THEN
3074  WRITE(lout,1014) nhkk,isthkk(nhkk),idhkk(nhkk),
3075  & jmohkk(1,nhkk),jmohkk(2,nhkk),
3076  & jdahkk(1,nhkk),jdahkk(2,nhkk)
3077  WRITE(lout,1015) (phkk(k,nhkk),k=1,5)
3078  1014 FORMAT('HADRDI: NHKK,ISTHKK,IDHKK,JMOHKK,JDAHKK',7i6)
3079  1015 FORMAT('HADRDI: PHKK ',5f10.5)
3080  ENDIF
3081  14 CONTINUE
3082  15 CONTINUE
3083  IF (nhad.GT.0) THEN
3084  jdahkk(1,imoh2) = nhkkbe
3085  jdahkk(2,imoh2) = nhkk
3086  ENDIF
3087 *
3088 *----------------------------- diffractive nucleon/antinucleon
3089 *
3090  naux = naux+1
3091  IF (idiftp.EQ.1) THEN
3092  IF(ijlmdd.EQ.1)THEN
3093  amc(naux) = aam(kdlmdd)
3094  ichc(naux) = iich(kdlmdd)
3095  ibarc(naux) = iibar(kdlmdd)
3096  anc(naux) = aname(kdlmdd)
3097  nrc(naux) =kdlmdd
3098  ELSEIF(ijlmdd.EQ.0)THEN
3099  amc(naux) = aam(ijproj)
3100  ichc(naux) = iich(ijproj)
3101  ibarc(naux) = iibar(ijproj)
3102  anc(naux) = aname(ijproj)
3103  nrc(naux) = ijproj
3104  ENDIF
3105  ELSEIF (idiftp.EQ.2) THEN
3106  IF(ijlmdd.EQ.1)THEN
3107  amc(naux) = aam(kdlmdd)
3108  ichc(naux) = iich(kdlmdd)
3109  ibarc(naux) = iibar(kdlmdd)
3110  anc(naux) = aname(kdlmdd)
3111  nrc(naux) =kdlmdd
3112  ELSEIF(ijlmdd.EQ.0)THEN
3113  amc(naux) = aam(ijtar)
3114  ichc(naux) = iich(ijtar)
3115  ibarc(naux) = iibar(ijtar)
3116  anc(naux) = aname(ijtar)
3117  nrc(naux) = ijtar
3118  ENDIF
3119  ENDIF
3120  pxc(naux) = pdfq1(1)
3121  pyc(naux) = pdfq1(2)
3122  pzc(naux) = pdfq1(3)
3123  hec(naux) = pdfq1(4)
3124  anndi=anndi+1
3125  eedi=eedi+hec(naux)
3126  ptdi=ptdi+sqrt(pxc(naux)**2+pyc(naux)**2)
3127  nhkk = nhkk+1
3128  isthkk(nhkk) = 1
3129 C IF(IBARF(J).EQ.500)ISTHKK(NHKK)=2
3130  idhkk(nhkk) = mpdgha(nrc(naux))
3131  IF(iormo(j).EQ.999)THEN
3132  jmohkk(1,nhkk) = imoh3
3133  ELSE
3134  jmohkk(1,nhkk)=nhkkbe+iormo(j)-1
3135  ENDIF
3136  jmohkk(2,nhkk) = 0
3137  jdahkk(1,nhkk) = 0
3138 *---------- S. Roesler 5-11-93
3139 * the following line is just to get the leading particle
3140 * detected
3141  jdahkk(2,nhkk) = 1000
3142 *
3143  phkk(1,nhkk) = pxc(naux)
3144  phkk(2,nhkk) = pyc(naux)
3145  phkk(3,nhkk) = pzc(naux)
3146  phkk(4,nhkk) = hec(naux)
3147  phkk(5,nhkk) = amc(naux)
3148  jdahkk(1,imoh3)= nhkk
3149  imohkk = jmohkk(1,nhkk)
3150  vhkk(1,nhkk) = vhkk(1,imohkk)
3151  vhkk(2,nhkk) = vhkk(2,imohkk)
3152  vhkk(3,nhkk) = vhkk(3,imohkk)
3153  vhkk(4,nhkk) = vhkk(4,imohkk)
3154  IF (iphkk.GE.2) THEN
3155  WRITE(lout,1016) nhkk,isthkk(nhkk),idhkk(nhkk),
3156  & jmohkk(1,nhkk),jmohkk(2,nhkk),
3157  & jdahkk(1,nhkk),jdahkk(2,nhkk)
3158  WRITE(lout,1017) (phkk(k,nhkk),k=1,5)
3159  1016 FORMAT('HADRDI: NHKK,ISTHKK,IDHKK,JMOHKK,JDAHKK',7i6)
3160  1017 FORMAT('HADRDI: PHKK ',5f10.5)
3161  ENDIF
3162  IF (ipev.GE.2) THEN
3163  DO 16 i=1,naux
3164  WRITE(lout,*) 'HADRDI: DIFFRACTIVE JET-HADRONS'
3165  WRITE(lout,1018)i,pxc(i),pyc(i),pzc(i),hec(i),amc(i),
3166  & ichc(i),ibarc(i),nrc(i),anc(i)
3167  1018 FORMAT(i5,5f12.4,3i5,a10)
3168  16 CONTINUE
3169  ENDIF
3170 *
3171  nnaux = naux
3172  amchdi = sqrt(edi**2-pxdi**2-pydi**2-pldi**2)
3173  IF (idiftp.EQ.1) ttt=2*aam(ijproj)**2-2.*(eproj*hec(naux)
3174  & -sqrt(eproj**2-aam(ijproj)**2)*pzc(naux))
3175  IF (idiftp.EQ.2) ttt=2*aam(ijtar)**2-2.*(etarg*hec(naux)
3176  & +sqrt(etarg**2-aam(ijtar)**2)*pzc(naux))
3177  IF (ipev.GE.2) THEN
3178  WRITE(lout,1019) amchdi,ttt
3179  1019 FORMAT('HADRDI: AMCHDI,TTT ',2f10.5)
3180  ENDIF
3181  9999 CONTINUE
3182  RETURN
3183  END
3184 *
3185 *===diadif===============================================================*
3186 *
3187  SUBROUTINE diadif(IOP,NHKKH1)
3188 
3189  IMPLICIT DOUBLE PRECISION (a-h,o-z)
3190  SAVE
3191  parameter(lout=6,llook=9)
3192  parameter(nmxhkk= 89998)
3193  CHARACTER*8 aname,anc
3194  COMMON /hkkevt/ nhkk,nevhkk, isthkk(nmxhkk), idhkk(nmxhkk),
3195  & jmohkk(2,nmxhkk),jdahkk(2,nmxhkk),phkk(5,nmxhkk),
3196  & vhkk(4,nmxhkk), whkk(4,nmxhkk)
3197  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
3198  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
3199  & iibar(210),k1(210),k2(210)
3200  COMMON /difpar/ pxc(902), pyc(902),pzc(902),
3201  & hec(902), amc(902),ichc(902),
3202  & ibarc(902),anc(902),nrc(902)
3203  COMMON /enerin/ eproj,etarg
3204  COMMON /nncms/ dgamcm,dbgcm,ecm,dpcm,deproj,dpproj
3205  COMMON /dhisto/ xyl(51,10),yyl(51,10),yylps(51,10),xxfl(51,10),
3206  & yxfl(51,10),tdtdm(40,24),dsdtdm(40,24),
3207  & amdm(24,40),dsdm(24,40),ave(30),avmult(30),
3208  & yxflch(51),yxflpi(51)
3209  COMMON /difout/ amch,tt,nhad,kproj,ktarg
3210  COMMON /diffra/ isingd,idiftp,ioudif,iflagd
3211  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
3212  COMMON /evflag/ numev
3213  dimension indx(28),px(902),py(902),pz(902),he(902),am(902),
3214  & ich(902),nr(902)
3215 C DATA INDX /1,8,10,10,10,10,7,2,7,10,10,7,3,4,5,6,
3216 C & 11,12,7,13,14,15,16,17,18,7,7,7/
3217  DATA indx/ 1,8,10,10,10, 10,7,2,7,10,
3218  & 10,7,3,4,5, 6,7,7,7,7,
3219  & 7,7,7,7,7, 7,7,7/
3220 *
3221  goto(1,2,3) iop
3222 *
3223  1 CONTINUE
3224  ncev = 0
3225  dxfl = 0.04d0
3226  dt = 0.1d0
3227  ddm = 10.0d0
3228  dy = 0.499999d0
3229  nevt = 0
3230  nch = 0
3231  npi = 0
3232  nhad = 0
3233 C KPROJ= 0
3234 C KTARG= 0
3235  amch = 0.0d0
3236  tt = 0.0d0
3237 *
3238  DO 10 j=1,51
3239  yxflch(j) = 0.0d0
3240  yxflpi(j) = 0.0d0
3241  DO 11 i=1,10
3242  xxfl(j,i) = j*dxfl-1.0d0
3243  yxfl(j,i) = 1.0d-18
3244  xyl(j,i) = (j-24)*dy-dy/2.0d0
3245  yyl(j,i) = 1.0d-18
3246  yylps(j,i)= 1.0d-18
3247  11 CONTINUE
3248  10 CONTINUE
3249  DO 12 i=1,40
3250  DO 13 j=1,24
3251  tdtdm(i,j) = i*dt
3252  dsdtdm(i,j)= 1.0d-8
3253  amdm(j,i) = (j*ddm-ddm/2.0d0)**2
3254  dsdm(j,i) = 1.0d-8
3255  amdm(j,i) = log10(amdm(j,i))
3256  13 CONTINUE
3257  12 CONTINUE
3258  DO 14 i=1,30
3259  ave(i) = 1.0d-18
3260  avmult(i) = 1.0d-18
3261  14 CONTINUE
3262  RETURN
3263 *
3264  2 CONTINUE
3265  ncev = ncev+1
3266  IF (mod(ncev,2000).EQ.0) WRITE(lout,*) ncev
3267  j = 0
3268 C IF (IFLAGD.EQ.1) RETURN
3269  DO 21 i=nhkkh1+1,nhkk
3270  IF ((isthkk(i).EQ.1).AND.(jmohkk(2,i).NE.100).AND.
3271  & (jdahkk(1,i).EQ.0)) THEN
3272  j = j+1
3273  px(j) = phkk(1,i)
3274  py(j) = phkk(2,i)
3275  pz(j) = phkk(3,i)
3276  he(j) = phkk(4,i)
3277  am(j) = phkk(5,i)
3278  nr(j) = mcihad(idhkk(i))
3279  ich(j)= iich(nr(j))
3280  ENDIF
3281  21 CONTINUE
3282  ihad = j
3283  IF (ipev.GE.2) THEN
3284  WRITE(lout,*) 'DIADIF: PX,PY,PZ,HE,AM,NR,ICH'
3285  DO 22 i=1,ihad
3286  WRITE(lout,*) px(i),py(i),pz(i),he(i),am(i),nr(i),ich(i)
3287  1000 FORMAT(5f12.5,2i4)
3288  22 CONTINUE
3289  ENDIF
3290 C IF (IDIFTP.EQ.1) THEN
3291 C P0 = SQRT(ETARG**2-AAM(KTARG)**2)
3292 C ELSE IF (IDIFTP.EQ.2) THEN
3293 C P0 = SQRT(EPROJ**2-AAM(KPROJ)**2)
3294 C ENDIF
3295  epcm = (aam(ijproj)**2-aam(1)**2+ecm**2)/(2.0d0*ecm)
3296  p0 = sqrt((epcm-aam(ijproj))*(epcm+aam(ijproj)))
3297  nevt = nevt+1
3298  avmult(30) = avmult(30)+ihad
3299  DO 20 i=1,ihad
3300  nre = nr(i)
3301  IF (nre.GT.25) nre = 28
3302  IF (nre.LT. 1) nre = 28
3303  ni = indx(nre)
3304  IF (nre.EQ.28) ni = 8
3305  ave(nre) = ave(nre)+he(i)
3306  ave(30) = ave(30) +he(i)
3307  IF (ni.NE.6) ave(29) = ave(29)+he(i)
3308  avmult(nre) = avmult(nre)+1.0d0
3309  IF (ni.NE.6) avmult(29) = avmult(29)+1.0d0
3310  IF (ich(i).NE.0) ave(27) = ave(27)+he(i)
3311  IF (ich(i).NE.0) avmult(27) = avmult(27)+1.0d0
3312 C XFL = PZ(I)/PO
3313  xfl = pz(i)/p0
3314  xfle = he(i)/p0
3315 C IF ((ICH(I).NE.0).AND.(XFL.GT.-0.84D0).AND.(XFL.LT.-0.44D0))
3316 C & WRITE(LOUT,*)'WRONG EVENT',NUMEV
3317  ixfl = xfl/dxfl+26
3318 C IF (XFL.LT.0.0D0) WRITE(LLOOK,'(2F15.5,I3)')XFL,PZ(I),IXFL
3319  IF (ixfl.LT.1 ) ixfl=1
3320  IF (ixfl.GT.50) ixfl=50
3321  xxxfl = abs(xfl)
3322  IF (nre.EQ.14) THEN
3323  npi = npi+1
3324  yxflpi(ixfl) = yxflpi(ixfl)+xfle
3325  ENDIF
3326  IF ((ich(i).GT.0).AND.(nre.NE.1)) THEN
3327  nch = nch+1
3328  yxflch(ixfl) = yxflch(ixfl)+xfle
3329  ENDIF
3330  IF (ich(i).NE.0) yxfl(ixfl,9) = yxfl(ixfl,9)+xxxfl
3331  yxfl(ixfl,ni) = yxfl(ixfl,ni)+xxxfl
3332  yxfl(ixfl,10) = yxfl(ixfl,10)+xxxfl
3333  ptt = px(i)**2+py(i)**2
3334  amt = sqrt(ptt+am(i)**2)
3335  yl = 0.5d0*log(abs((he(i)+pz(i)+1.d-10)
3336  & /(he(i)-pz(i)+1.d-10)))
3337  ylps = log(abs((pz(i)+sqrt(pz(i)**2+ptt))
3338  & /sqrt(ptt+1.d-6)+1.d-18))
3339  iylps = (ylps+25.0d0*dy)/dy
3340  IF (iylps.LT.1) iylps = 1
3341  IF (iylps.GT.51) iylps = 51
3342  yylps(iylps,ni) = yylps(iylps,ni)+1.0d0
3343  yylps(iylps,10) = yylps(iylps,10)+1.0d0
3344  IF (ich(i).NE.0) yylps(iylps,9) = yylps(iylps,9)+1.0d0
3345  iyl = (yl+25.0d0*dy)/dy
3346  IF (iyl.LT.1) iyl = 1
3347  IF (iyl.GT.51) iyl = 51
3348  IF (ich(i).NE.0) yyl(iyl,9) = yyl(iyl,9)+1.0d0
3349  yyl(iyl,ni) = yyl(iyl,ni)+1.0d0
3350  yyl(iyl,10) = yyl(iyl,10)+1.0d0
3351  20 CONTINUE
3352  kpl = (amch+ddm+ddm/2.0d0)/ddm
3353  IF (kpl.GT.24) kpl=24
3354  IF (kpl.LT.1) kpl=1
3355  itt = 10.0d0*abs(tt)+1
3356  IF (itt.GT.40) itt = 40
3357 C DSDTDM(ITT,KPL) = DSDTDM(ITT,KPL)+1.0D0/AMCH
3358  RETURN
3359 *
3360  3 CONTINUE
3361  DO 30 j=1,51
3362  IF(nch.NE.0) yxflch(j) = yxflch(j)/(nch*dxfl)
3363  IF(npi.NE.0) yxflpi(j) = yxflpi(j)/(npi*dxfl)
3364  DO 31 i=1,10
3365  yxfl(j,i) = log10(abs(yxfl(j,i) /(nevt*dxfl))+1.0d-8)
3366  yyl(j,i) = yyl(j,i) /(nevt*dy)
3367  yylps(j,i)= yylps(j,i)/(nevt*dy)
3368  31 CONTINUE
3369  30 CONTINUE
3370  DO 32 i=1,30
3371  avmult(i) = avmult(i)/nevt
3372  ave( i) = ave(i) /nevt
3373  32 CONTINUE
3374  DO 33 i=1,40
3375  DO 34 j=1,24
3376  dsdtdm(i,j) = dsdtdm(i,j)/nevt
3377  dsdtdm(i,j) = log10(dsdtdm(i,j))
3378  dsdm(j,i) = dsdtdm(i,j)
3379  34 CONTINUE
3380  33 CONTINUE
3381  WRITE(lout,*)'NAME,AVE,AVMULT'
3382  DO 105 i=1,30
3383  WRITE(lout,210) aname(i),ave(i),avmult(i)
3384  210 FORMAT(' ',a8,2f15.5)
3385  105 CONTINUE
3386  WRITE(lout,100)
3387  100 FORMAT('1 RAPIDITY DISTRIBUTION')
3388  DO 110 j=1,50
3389  WRITE(lout,200) xyl(j,1),(yyl(j,i),i=1,10)
3390 c WRITE(LOUT,200) XYL(J,1),(YYL(J,I),I=11,20)
3391  200 FORMAT (f10.2,10e11.3)
3392  110 CONTINUE
3393  CALL plot(xyl,yyl,510,10,51,-25.*dy,dy,0.,0.03)
3394  WRITE(lout,101)
3395  101 FORMAT('1 PSEUDORAPIDITY DISTRIBUTION')
3396  CALL plot(xyl,yylps,510,10,51,-25.*dy,dy,0.,0.03)
3397  WRITE(lout,102)
3398  102 FORMAT ('1 LONG MOMENTUM (SCALED) DISTRIBUTION (LOG)')
3399  CALL plot(xxfl,yxfl,510,10,51,-1.,dxfl,-3.5,0.05)
3400  WRITE(lout,103)
3401  103 FORMAT ('1DISTRIBUTION DS/DTDM AS FUNCTION OF T ')
3402  CALL plot(tdtdm,dsdtdm,960,24,40,0.,0.04,-5.,0.05)
3403  WRITE(lout,104)
3404  104 FORMAT ('1DISTRIBUTION DS/DTDM AS FUNCTION OF M**2 ')
3405  CALL plot(amdm,dsdm,960,40,24,2.,0.1,-5.,0.05)
3406  IF (ijproj.EQ.13) sig = 66.0d0
3407  IF (ijproj.EQ.14) sig = 66.0d0
3408  IF (ijproj.EQ.15) sig = 54.9d0
3409  IF (ijproj.EQ.1) sig = 95.0d0
3410  sigch = 84.8
3411  sigpim=0.23d0
3412 C OPEN(19,FILE='FEYNDI.OUT')
3413  WRITE(lout,*)'FEYNMAN-DISTRIBUTION FOR PION-'
3414  DO 35 i=1,51
3415  WRITE(lout,'(2F15.5)')
3416  & dxfl*(i-1)-1,sig*yxflpi(i)
3417  WRITE(lout,'(2F15.5)')
3418  & dxfl*i-1,sig*yxflpi(i)
3419  35 CONTINUE
3420  WRITE(lout,*)'FEYNMAN-DISTRIBUTION FOR CHARGED PARTICLES'
3421  DO 36 i=1,51
3422  WRITE(lout,'(2F15.5)')
3423  & dxfl*(i-1)-1,sigch*yxflch(i)
3424  WRITE(lout,'(2F15.5)')
3425  & dxfl*i-1,sigch*yxflch(i)
3426  36 CONTINUE
3427 C CLOSE(19)
3428  RETURN
3429  END
3430 *
3431 *===sihndi===============================================================*
3432 *
3433  SUBROUTINE sihndi(ECM,KPROJ,KTARG,SIGDIF,SIGDIH)
3434 
3435 **********************************************************************
3436 * Single diffractive hadron-nucleon cross sections *
3437 * S.Roesler 14/1/93 *
3438 * *
3439 * The cross sections are calculated from extrapolated single *
3440 * diffractive antiproton-proton cross sections (DTUJET92) using *
3441 * scaling relations between total and single diffractive cross *
3442 * sections. *
3443 **********************************************************************
3444 
3445  IMPLICIT DOUBLE PRECISION (a-h,o-z)
3446  SAVE
3447  CHARACTER*8 aname
3448  COMMON /dpar/ aname(210),am(210),ga(210),tau(210),iich(210),
3449  & iibar(210),k1(210),k2(210)
3450 *
3451  csd1 = 4.201483727
3452  csd4 = -0.4763103556e-02
3453  csd5 = 0.4324148297
3454 *
3455  chmsd1 = 0.8519297242
3456  chmsd4 = -0.1443076599e-01
3457  chmsd5 = 0.4014954567
3458 *
3459  epn = (ecm**2 - am(kproj)**2 - am(ktarg)**2)/(2.0d0*am(ktarg))
3460  ppn = sqrt((epn-am(kproj))*(epn+am(kproj)))
3461 *
3462  sdiapp = csd1+csd4*log(ppn)**2+csd5*log(ppn)
3463  shmsd = chmsd1+chmsd4*log(ppn)**2+chmsd5*log(ppn)
3464  frac = shmsd/sdiapp
3465 *
3466  goto( 10, 20,999,999,999,999,999, 10, 20,999,
3467  & 999, 20, 20, 20, 20, 20, 10, 20, 20, 10,
3468  & 10, 10, 20, 20, 20) kproj
3469 *
3470  10 CONTINUE
3471 *---------------------------- p - p , n - p , sigma0+- - p ,
3472 * Lambda - p
3473  csd1 = 6.004476070
3474  csd4 = -0.1257784606e-03
3475  csd5 = 0.2447335720
3476  sigdif = csd1+csd4*log(ppn)**2+csd5*log(ppn)
3477 C Replace SIGDIF with dpmjet result SIPPSD(ECM)
3478  sigdif=sippsd(ecm)
3479  sigdih = frac*sigdif
3480  RETURN
3481 *
3482  20 CONTINUE
3483 *
3484  kpscal = 2
3485  ktscal = 1
3486  f = sdiapp/dshnto(kpscal,ktscal,ecm)
3487 C F = SDIAPP/DSHNEL(KPSCAL,KTSCAL,ECM)
3488  kt = 1
3489  sigdif = dshnto(kproj,kt,ecm)*f
3490 C SIGDIF = DSHNEL(KPROJ,KT,ECM)*F
3491  sigdih = frac*sigdif
3492  RETURN
3493 *
3494  999 CONTINUE
3495 *-------------------------- leptons..
3496  sigdif = 1.e-10
3497  sigdih = 1.e-10
3498  RETURN
3499  END
3500 *
3501 *===dshnto===============================================================*
3502 *
3503  DOUBLE PRECISION FUNCTION dshnto(KPROJ,KTARG,UMO)
3504 
3505 **********************************************************************
3506 * Total hadron-nucleon cross sections S.Roesler 12/1/93 *
3507 * *
3508 * Fits from Rev. of Part. Prop.(1992) are used in the momentum *
3509 * ranges given there. Cross sections for momenta below them are *
3510 * calculated as in an earlier version of this function : *
3511 * SIG(el) + SIG(inel) = SIHNEL + SIHNIN . *
3512 * Total antiproton-proton cross sections calculated with DTUJET92 *
3513 * parametrize the cross sections at higher enrgies. *
3514 **********************************************************************
3515 
3516  IMPLICIT DOUBLE PRECISION (a-h,o-z)
3517  SAVE
3518  CHARACTER*8 aname
3519  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
3520  & iibar(210),k1(210),k2(210)
3521  COMMON /strufu/istrum,istrut
3522  dimension sqs(20),siiv(20),sqs22(40),siiv22(40)
3523  DATA sqs /20.,50.,100.,200.,500.,1000.,1500.,2000.,3000.,
3524  *4000.,6000.,8000.,10000.,15000.,20000.,30000.,40000.,
3525  *60000.,80000.,100000/
3526  DATA siiv /41.6,44.6,47.9,52.3,60.,67.2,71.4,75.,79.,82.4,
3527  *87.2,90.4,93.2,97.3,100.7,104.7,107.9,111.7,114.7,117.2/
3528  DATA sqs22 /
3529  *53.,69.,91.,119.,156.,205.,268.,351.,460.,603.,
3530  *790.,1036.,1357.,1778.,2329.,
3531  *3053.,3999.,5239.,6865.,8994.,
3532  *11785.,15441.,20232.,26509.,34733.,
3533  *45509.,59627.,78126.,102365.,134123.,
3534  *175734.,230255.,301690.,395288.,517925.,
3535  *678609.,889144.,1164997.,1526432.,2000000.
3536  */
3537  DATA siiv22 /
3538  *44.3,45.3,46.5,47.8,49.3,51.0,53.0,55.3,57.9,60.9,
3539  *64.3,68.0,72.0,76.3,80.8,85.4,90.0,94.7,99.3,103.9,
3540  *108.4,112.8,117.2,121.4,125.6,
3541  *129.8,133.9,138.0,142.0,146.1,
3542  *150.2,154.3,158.5,162.7,166.9,
3543  *171.2,175.6,180.0,184.6,189.1
3544  */
3545 *
3546  f1 = 1.0d0
3547  ca = 0.0d0
3548  cb = 0.0d0
3549  cc = 0.0d0
3550  cd = 0.0d0
3551  cn = 0.0d0
3552 *
3553  a1 = 0.0d0
3554  a2 = 0.0d0
3555  a3 = 1.0d0
3556  a4 = 0.0d0
3557  a5 = 0.0d0
3558  a6 = 0.0d0
3559 *
3560  param1 = 34.94235992
3561  param4 = 0.2104312854
3562  param5 = -0.4509592056e-01
3563 *
3564  ipio = 0
3565  spio1 = 0.0d0
3566  spio2 = 0.0d0
3567 *
3568  umo2 = umo**2
3569  epn = (umo**2-aam(kproj)**2-aam(ktarg)**2)/(2.0d0*aam(ktarg))
3570  po = sqrt((epn-aam(kproj))*(epn+aam(kproj)))
3571 *
3572  1 CONTINUE
3573 *
3574  IF (ktarg.EQ.8) THEN
3575  goto( 30, 40,999,999,999,999,999, 10, 20,999,
3576  & 999,140, 70, 60,150,160,100, 20,140, 10,
3577  & 10, 10,110,130,120) kproj
3578  ELSE
3579  goto( 10, 20,999,999,999,999,999, 30, 40,999,
3580  & 999, 50, 60, 70, 80, 90,100, 20, 50, 10,
3581  & 10, 10,110,120,130) kproj
3582  ENDIF
3583 *
3584  10 CONTINUE
3585 *---------------------------- p - p , sigma0+- - p
3586  IF (po.LE.3.0d0) THEN
3587  goto 500
3588  ELSE IF ((po.GT.3.0d0).AND.(umo.LE.100.0d0)) THEN
3589  ca = 48.0d0
3590  cc = 0.522d0
3591  cd = -4.51d0
3592  goto 600
3593  ELSE
3594  goto 700
3595  ENDIF
3596 *
3597  20 CONTINUE
3598 *---------------------------- pbar - p , Lambdabar - p
3599  IF (po.LE.5.0d0) THEN
3600  goto 500
3601  ELSE IF ((po.GT.5.0d0).AND.(umo.LE.200.0d0)) THEN
3602  ca = 38.4d0
3603  cb = 77.6
3604  cc = 0.26d0
3605  cd = -1.2d0
3606  cn = -0.64d0
3607  goto 600
3608  ELSE
3609  goto 700
3610  ENDIF
3611 *
3612  30 CONTINUE
3613 *---------------------------- n - p
3614  IF (po.LE.3.0d0) THEN
3615  goto 500
3616  ELSE IF ((po.GT.3.0d0).AND.(po.LE.370.0d0)) THEN
3617  ca = 47.3d0
3618  cc = 0.513d0
3619  cd = -4.27d0
3620  goto 600
3621  ELSE IF ((po.GT.370.0d0).AND.(umo.LE.110.0d0)) THEN
3622  a1 = 38.5d0
3623  a2 = 0.46d0
3624  a3 = 125.0d0
3625  a4 = 15.0d0
3626  goto 800
3627  ELSE
3628  goto 700
3629  ENDIF
3630 *
3631  40 CONTINUE
3632 *---------------------------- nbar - p
3633  IF (po.LE.1.1d0) THEN
3634  goto 500
3635  ELSE IF ((po.GT.1.1d0).AND.(po.LE.280.0d0)) THEN
3636  cb = 133.6d0
3637  cc = -1.22d0
3638  cd = 13.7d0
3639  cn = -0.7d0
3640  goto 600
3641  ELSE IF ((po.GT.280.0d0).AND.(umo.LE.110.0d0)) THEN
3642  a1 = 38.5d0
3643  a2 = 0.46d0
3644  a3 = 125.0d0
3645  a4 = 15.0d0
3646  a5 = 77.43d0
3647  a6 = -0.6d0
3648  goto 800
3649  ELSE
3650  goto 700
3651  ENDIF
3652 *
3653  50 CONTINUE
3654 *---------------------------- Klong - p , Kshort - p
3655  r = rndm(v)
3656  IF (r.LE.0.5d0) THEN
3657 * K+ - p
3658  goto 80
3659  ELSE
3660 * K- - p
3661  goto 90
3662  ENDIF
3663 *
3664  60 CONTINUE
3665 *---------------------------- pi+ - p
3666  IF (po.LE.4.0d0) THEN
3667  goto 500
3668  ELSE IF ((po.GT.4.0d0).AND.(po.LE.340.0d0)) THEN
3669  ca = 16.4d0
3670  cb = 19.3d0
3671  cc = 0.19d0
3672  cn = -0.42d0
3673  goto 600
3674  ELSE IF ((po.GT.340.0d0).AND.(umo.LE.47.0d0)) THEN
3675  a1 = 24.0d0
3676  a2 = 0.6d0
3677  a3 = 160.0d0
3678  a5 = -7.9d0
3679  a6 = -0.46d0
3680  goto 800
3681  ELSE
3682  f1 = 2.0d0/3.0d0
3683  goto 10
3684  ENDIF
3685 *
3686  70 CONTINUE
3687 *---------------------------- pi- - p
3688  IF (po.LE.2.5d0) THEN
3689  goto 500
3690  ELSE IF ((po.GT.2.5d0).AND.(po.LE.370.0d0)) THEN
3691  ca = 33.0d0
3692  cb = 14.0d0
3693  cc = 0.456d0
3694  cd = -4.03d0
3695  cn = -1.36d0
3696  goto 600
3697  ELSE IF ((po.GT.370.0d0).AND.(umo.LE.47.0d0)) THEN
3698  a1 = 24.0d0
3699  a2 = 0.6d0
3700  a3 = 160.0d0
3701  goto 800
3702  ELSE
3703  f1 = 2.0d0/3.0d0
3704  goto 10
3705  ENDIF
3706 *
3707  80 CONTINUE
3708 *---------------------------- K+ - p
3709  IF (po.LE.2.0d0) THEN
3710  goto 500
3711  ELSE IF ((po.GT.2.0d0).AND.(po.LE.310.0d0)) THEN
3712  ca = 18.1d0
3713  cc = 0.26d0
3714  cd = -1.0d0
3715  goto 600
3716  ELSE IF ((po.GT.310.0d0).AND.(umo.LE.110.0d0)) THEN
3717  a1 = 20.3d0
3718  a2 = 0.59d0
3719  a3 = 140.0d0
3720  a5 = -30.13d0
3721  a6 = -0.58d0
3722  goto 800
3723  ELSE
3724  f1 = 2.0d0/3.0d0
3725  goto 10
3726  ENDIF
3727 *
3728  90 CONTINUE
3729 *---------------------------- K- - p
3730  IF (po.LE.3.0d0) THEN
3731  goto 500
3732  ELSE IF ((po.GT.3.0d0).AND.(po.LE.310.0d0)) THEN
3733  ca = 32.1d0
3734  cc = 0.66d0
3735  cd = -5.6d0
3736  goto 600
3737  ELSE IF ((po.GT.310.0d0).AND.(umo.LE.110.0d0)) THEN
3738  a1 = 20.3d0
3739  a2 = 0.59d0
3740  a3 = 140.0d0
3741  goto 800
3742  ELSE
3743  f1 = 2.0d0/3.0d0
3744  goto 10
3745  ENDIF
3746 *
3747  100 CONTINUE
3748 *---------------------------- Lambda - p
3749  IF (po.LE.0.6d0) THEN
3750  goto 500
3751  ELSE IF ((po.GT.0.6d0).AND.(po.LE.21.0d0)) THEN
3752  ca = 30.4d0
3753  cd = 1.6d0
3754  goto 600
3755  ELSE
3756  goto 10
3757  ENDIF
3758 *
3759  110 CONTINUE
3760 *---------------------------- pi0 - p
3761 * 1/2(pi+p + pi-p)
3762  ipio = 1
3763  kproj = 13
3764  goto 1
3765 *
3766  120 CONTINUE
3767 *---------------------------- K0 - p
3768  IF (po.LE.2.0d0) THEN
3769  goto 500
3770  ELSE IF ((po.GT.2.0d0).AND.(po.LE.310.0d0)) THEN
3771 * K+ - n
3772  ca = 18.7d0
3773  cc = 0.21d0
3774  cd = -0.89d0
3775  goto 600
3776  ELSE
3777  goto 80
3778  ENDIF
3779 *
3780  130 CONTINUE
3781 *---------------------------- K0bar - p
3782  IF (po.LE.1.8d0) THEN
3783  goto 500
3784  ELSE IF ((po.GT.1.8d0).AND.(po.LE.310.0d0)) THEN
3785 * K- - n
3786  ca = 25.2d0
3787  cc = 0.38d0
3788  cd = -2.9d0
3789  goto 600
3790  ELSE
3791  goto 90
3792  ENDIF
3793 *
3794  140 CONTINUE
3795 *---------------------------- Klong - n , Kshort - n
3796  r = rndm(v)
3797  IF (r.LE.0.5d0) THEN
3798 * K+ - n
3799  goto 150
3800  ELSE
3801 * K- - n
3802  goto 160
3803  ENDIF
3804 *
3805  150 CONTINUE
3806 *---------------------------- K+ - n
3807  IF (po.LE.2.0d0) THEN
3808  goto 500
3809  ELSE IF ((po.GT.2.0d0).AND.(po.LE.310.0d0)) THEN
3810  ca = 18.7d0
3811  cc = 0.21d0
3812  cd = -0.89d0
3813  goto 600
3814  ELSE
3815  goto 90
3816  ENDIF
3817 *
3818  160 CONTINUE
3819 *---------------------------- K- - n
3820  IF (po.LE.1.8d0) THEN
3821  goto 500
3822  ELSE IF ((po.GT.1.8d0).AND.(po.LE.310.0d0)) THEN
3823  ca = 25.2d0
3824  cc = 0.38d0
3825  cd = -2.9d0
3826  goto 600
3827  ELSE
3828  goto 80
3829  ENDIF
3830 *
3831  500 CONTINUE
3832  CALL sihnel(kproj,ktarg,po,sel)
3833  CALL sihnin(kproj,ktarg,po,sin)
3834  stot = sel + sin
3835  goto 900
3836 *
3837  600 CONTINUE
3838  stot = f1*(ca+cb*po**cn+cc*(log(po))**2+cd*log(po))
3839  goto 900
3840 *
3841  700 CONTINUE
3842  IF(istrum.EQ.14.AND.istrut.EQ.2)THEN
3843  DO 701 i=1,20
3844  IF(umo.LE.sqs(i))go to 702
3845  701 CONTINUE
3846  i=20
3847  702 CONTINUE
3848  IF ((i.EQ.20).AND.(umo.GT.sqs(20)))THEN
3849  tepp=siiv(20)+(log(umo)-log(sqs(20)))*(siiv(20)-siiv(19))/
3850  * (log(sqs(20))-log(sqs(19)))
3851  stot=f1*tepp
3852  ELSE
3853  tepp=siiv(i-1)+(umo-sqs(i-1))*(siiv(i)-siiv(i-1))/
3854  * (sqs(i)-sqs(i-1))
3855  stot=f1*tepp
3856  ENDIF
3857  ELSEIF(istrum.EQ.22.AND.istrut.EQ.2)THEN
3858  DO 711 i=1,40
3859  IF(umo.LE.sqs22(i))go to 712
3860  711 CONTINUE
3861  i=40
3862  712 CONTINUE
3863  IF ((i.EQ.40).AND.(umo.GT.sqs22(40)))THEN
3864  tepp=siiv22(40)+(log(umo)-log(sqs22(40)))*
3865  * (siiv22(40)-siiv22(39))/
3866  * (log(sqs22(40))-log(sqs22(39)))
3867  stot=f1*tepp
3868  ELSE
3869  tepp=siiv22(i-1)+(umo-sqs22(i-1))*(siiv22(i)-siiv22(i-1))/
3870  * (sqs22(i)-sqs22(i-1))
3871  stot=f1*tepp
3872  ENDIF
3873  ELSE
3874  stot = f1*(param1+param4*log(po)**2+param5*log(po))
3875  ENDIF
3876  goto 900
3877 *
3878  800 CONTINUE
3879  stot = f1*(a1+a2*(log(umo2/a3))**2+a4/umo2+a5*umo2**a6)
3880 *
3881  900 CONTINUE
3882  IF ((ipio.EQ.1).AND.(kproj.EQ.13)) THEN
3883  spio1 = stot
3884  kproj = 14
3885  goto 1
3886  ENDIF
3887  IF ((ipio.EQ.1).AND.(kproj.EQ.14)) THEN
3888  spio2 = stot
3889  stot = 0.5d0*(spio1+spio2)
3890  ENDIF
3891  dshnto = stot
3892  RETURN
3893 *
3894  999 CONTINUE
3895 *-------------------------- leptons..
3896  dshnto = 1.e-10
3897  RETURN
3898  END
3899 *
3900 *===dshnel===============================================================*
3901 *
3902  DOUBLE PRECISION FUNCTION dshnel(KPROJ,KTARG,UMO)
3903 
3904 **********************************************************************
3905 * Elastic hadron-nucleon cross sections S.Roesler 13/1/93 *
3906 * *
3907 * SIHNEL is called for c.m. energies below 50 GeV. Otherwise *
3908 * a parametrisation of antiproton-proton elastic cross sections *
3909 * obtained from DTUJET92 is used. *
3910 **********************************************************************
3911 
3912  IMPLICIT DOUBLE PRECISION (a-h,o-z)
3913  SAVE
3914  CHARACTER*8 aname
3915  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
3916  & iibar(210),k1(210),k2(210)
3917 *
3918  f1 = 1.0d0
3919  ca = 0.0d0
3920  cb = 0.0d0
3921  cc = 0.0d0
3922  cd = 0.0d0
3923  cn = 0.0d0
3924 *
3925  ipio = 0
3926  spio1 = 0.0d0
3927  spio2 = 0.0d0
3928 *
3929  param1 = 7.789333344
3930  param4 = 0.7488331199e-01
3931  param5 = -0.6963931322
3932 *
3933  umo2 = umo**2
3934  epn = (umo**2-aam(kproj)**2-aam(ktarg)**2)/(2.0d0*aam(ktarg))
3935  po = sqrt((epn-aam(kproj))*(epn+aam(kproj)))
3936 *
3937  1 CONTINUE
3938 *
3939  goto( 10, 10,999,999,999,999,999, 10, 10,999,
3940  & 999, 30, 20, 20, 20, 20, 10, 10, 30, 10,
3941  & 10, 10, 40, 20, 20) kproj
3942 *
3943  10 CONTINUE
3944  IF (umo.LE.50.0d0) THEN
3945  CALL sihnel(kproj,ktarg,po,sel)
3946  goto 900
3947  ELSE
3948  f1 = 1.0d0
3949  goto 500
3950  ENDIF
3951 *
3952  20 CONTINUE
3953  IF (umo.LE.50.0d0) THEN
3954  CALL sihnel(kproj,ktarg,po,sel)
3955  goto 900
3956  ELSE
3957  f1 = 2.0d0/3.0d0
3958  goto 500
3959  ENDIF
3960 *
3961  30 CONTINUE
3962 *---------------------------- Klong - p , Kshort - p
3963  r = rndm(v)
3964  IF (r.LE.0.5d0) THEN
3965 * K+ - p
3966  kproj = 15
3967  ELSE
3968 * K- - p
3969  kproj = 16
3970  ENDIF
3971  goto 1
3972 *
3973  40 CONTINUE
3974 *---------------------------- pi0 - p
3975 * 1/2(pi+p + pi-p)
3976  ipio = 1
3977  kproj = 13
3978  goto 1
3979 *
3980  500 CONTINUE
3981  sel = f1*(param1+param4*log(po)**2+param5*log(po))
3982 *
3983  900 CONTINUE
3984  IF ((ipio.EQ.1).AND.(kproj.EQ.13)) THEN
3985  spio1 = sel
3986  kproj = 14
3987  goto 1
3988  ENDIF
3989  IF ((ipio.EQ.1).AND.(kproj.EQ.14)) THEN
3990  spio2 = sel
3991  sel = 0.5d0*(spio1+spio2)
3992  ENDIF
3993  dshnel = sel
3994  RETURN
3995 *
3996  999 CONTINUE
3997 *-------------------------- leptons..
3998  dshnel = 1.e-10
3999  RETURN
4000  END
Generated on Sat May 25 2013 14:32:27 for Geant4 by   doxygen 1.8.4