Geant4_10
dpm25nuc4.f
Go to the documentation of this file.
1 *CMZ : 1.00/00 27/11/91 17.09.38 by H.-J. M¶hring
2 *-- Author :
3 C-------------------------------------------------------------------
4 C
5 C FILE DMNUC4.FOR
6 C
7 C-------------------------------------------------------------------
8 C
9 C Sample according to Poisson distribution
10  FUNCTION npoiss(AVN)
11  IMPLICIT DOUBLE PRECISION (a-h,o-z)
12  expavn=exp(-avn)
13  k=1
14  a=1.
15  10 CONTINUE
16  u=rndm(v)
17  a=u*a
18  IF(a.LT.expavn)THEN
19  npoiss=k-1
20  RETURN
21  ELSE
22  k=k+1
23  go to 10
24  ENDIF
25  END
26 C
27 C--------------------------------------------------------------------
28 C------------------------------------------------------------------
29 C
30 C FILE DPMSEWEW FORTRAN
31 C
32 C------------------------------------------------------------------
33  SUBROUTINE sewew(IOP,NHKKH1)
34  IMPLICIT DOUBLE PRECISION (a-h,o-z)
35  SAVE
36 C
37 C*** 1=P, 2=N, 3=PI+, 4=PI-, 5=PIO, 6=GAM+HYP, 7=K, 8=ANUC, 9=CHARGED
38 C*** 10=TOT, 11=TOTHAD
39 C
40 *KEEP,HKKEVT.
41 c INCLUDE (HKKEVT)
42  parameter(tiny= 1.d-5)
43  parameter(nmxhkk= 89998)
44 c PARAMETER (NMXHKK=25000)
45  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
46  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
47  +(4,nmxhkk)
48  COMMON /extevt/ idres(nmxhkk),idxres(nmxhkk),nobam(nmxhkk),
49  & idbam(nmxhkk),idch(nmxhkk),npoint(10)
50  COMMON /wewevt/ iywew(nmxhkk),irwew(nmxhkk),iiywew(2600),
51  & iirwew(2600),wewew(2600),idwew1(2600),
52  & idwew2(2600),idwew(2600),iblock(nmxhkk),
53  & idweww(2600)
54 C
55 C
56 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
57 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
58 C THE POSITIONS OF THE PROJECTILE NUCLEONS
59 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
60 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
61 C COMPLETELY CONSISTENT. THE TIMES IN THE
62 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
63 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
64 C
65 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
66 C
67 C NMXHKK: maximum numbers of entries (partons/particles) that can be
68 C stored in the commonblock.
69 C
70 C NHKK: the actual number of entries stored in current event. These are
71 C found in the first NHKK positions of the respective arrays below.
72 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
73 C entry.
74 C
75 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
76 C = 0 : null entry.
77 C = 1 : an existing entry, which has not decayed or fragmented.
78 C This is the main class of entries which represents the
79 C "final state" given by the generator.
80 C = 2 : an entry which has decayed or fragmented and therefore
81 C is not appearing in the final state, but is retained for
82 C event history information.
83 C = 3 : a documentation line, defined separately from the event
84 C history. (incoming reacting
85 C particles, etc.)
86 C = 4 - 10 : undefined, but reserved for future standards.
87 C = 11 - 20 : at the disposal of each model builder for constructs
88 C specific to his program, but equivalent to a null line in the
89 C context of any other program. One example is the cone defining
90 C vector of HERWIG, another cluster or event axes of the JETSET
91 C analysis routines.
92 C = 21 - : at the disposal of users, in particular for event tracking
93 C in the detector.
94 C
95 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
96 C standard.
97 C
98 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
99 C The value is 0 for initial entries.
100 C
101 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
102 C one mother exist, in which case the value 0 is used. In cluster
103 C fragmentation models, the two mothers would correspond to the q
104 C and qbar which join to form a cluster. In string fragmentation,
105 C the two mothers of a particle produced in the fragmentation would
106 C be the two endpoints of the string (with the range in between
107 C implied).
108 C
109 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
110 C entry has not decayed, this is 0.
111 C
112 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
113 C entry has not decayed, this is 0. It is assumed that the daughters
114 C of a particle (or cluster or string) are stored sequentially, so
115 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
116 C daughters. Even in cases where only one daughter is defined (e.g.
117 C K0 -> K0S) both values should be defined, to make for a uniform
118 C approach in terms of loop constructions.
119 C
120 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
121 C
122 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
123 C
124 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
125 C
126 C PHKK(4,IHKK) : energy, in GeV.
127 C
128 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
129 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
130 C
131 C VHKK(1,IHKK) : production vertex x position, in mm.
132 C
133 C VHKK(2,IHKK) : production vertex y position, in mm.
134 C
135 C VHKK(3,IHKK) : production vertex z position, in mm.
136 C
137 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
138 C********************************************************************
139 *KEEP,DPAR.
140 C /DPAR/ CONTAINS PARTICLE PROPERTIES
141 C ANAME = LITERAL NAME OF THE PARTICLE
142 C AAM = PARTICLE MASS IN GEV
143 C GA = DECAY WIDTH
144 C TAU = LIFE TIME OF INSTABLE PARTICLES
145 C IICH = ELECTRIC CHARGE OF THE PARTICLE
146 C IIBAR = BARYON NUMBER
147 C K1,K1 = BIGIN AND END OF DECAY CHANNELS OF PARTICLE
148 C
149  CHARACTER*8 aname
150  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
151  +iibar(210),k1(210),k2(210)
152 C------------------
153 *KEEP,NUCC.
154  COMMON /nucc/ it,itz,ip,ipz,ijproj,ibproj,ijtarg,ibtarg
155 *KEEP,TNUCFI.
156  COMMON /tnucfi/ amtfin,amfino,eexcta,ptfinx,ptfiny,ptfinz,ptfine,
157  + itfin,itzfin,nntar,nptar,nhtar,nqtar,nparf
158 *KEEP,DPRIN.
159  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
160 *KEND.
161 
162 C modified DPMJET
163  COMMON /bufues/ bnnvv,bnnss,bnnsv,bnnvs,bnncc,
164  * bnndv,bnnvd,bnnds,bnnsd,
165  * bnnhh,bnnzz,
166  * bptvv,bptss,bptsv,bptvs,bptcc,bptdv,
167  * bptvd,bptds,bptsd,
168  * bpthh,bptzz,
169  * beevv,beess,beesv,beevs,beecc,beedv,
170  * beevd,beeds,beesd,
171  * beehh,beezz
172  * ,bnndi,bptdi,beedi
173  * ,bnnzd,bnndz,bptzd,bptdz,beezd,beedz
174  COMMON /ncoucs/ bcouvv,bcouss,bcousv,bcouvs,
175  * bcouzz,bcouhh,bcouds,bcousd,
176  * bcoudz,bcouzd,bcoudi,
177  * bcoudv,bcouvd,bcoucc
178  COMMON /bufueh/ annvv,annss,annsv,annvs,anncc,
179  * anndv,annvd,annds,annsd,
180  * annhh,annzz,
181  * ptvv,ptss,ptsv,ptvs,ptcc,ptdv,ptvd,ptds,ptsd,
182  * pthh,ptzz,
183  * eevv,eess,eesv,eevs,eecc,eedv,eevd,eeds,eesd,
184  * eehh,eezz
185  * ,anndi,ptdi,eedi
186  * ,annzd,anndz,ptzd,ptdz,eezd,eedz
187  COMMON /ncouch/ acouvv,acouss,acousv,acouvs,
188  * acouzz,acouhh,acouds,acousd,
189  * acoudz,acouzd,acoudi,
190  * acoudv,acouvd,acoucc
191 C---------------------
192  COMMON /eventa/idumtp
193  COMMON /hboo/ihbook
194 C 1 2 Transverse rapidity densities (rap, p,pi-, trans. dist. (fermi))
195 C 3 4 Transverse rapidity densities (rap,ap,pi+, trans. dist. (fermi))
196 C 5 6 Transverse rapidity densities (rap,La,aLa, trans. dist. (fermi))
197 C 7 8 Transverse rapidity densities (rap,n ,an , trans. dist. (fermi))
198 C 9 Transverse rapidity densities (rap,pi0 , trans. dist. (fermi))
199  dimension yyltra(40,9,41),transa(40),dyylam(40),dyalam(40),
200  *xyltra(40,9,41)
201 C
202  COMMON /sigla/siglau
203  dimension indx(28)
204 
205  DATA ipriop/0/
206  DATA indx/1,8,10,10,10,10,7,2,7,10,10,7,3,4,5,6,
207  * 7,7,7,7,7,7,7,7,7,7,7,7/
208 C----------------------------------------------------------------------
209  DATA kpl /1/
210  DATA ievl /0/
211 C----------------------------------------------------------------------
212  zero=0
213  oneone=1
214  twotwo=2
215  nip=ip
216  aip=ip
217 C Transverse Rapidity density
218 C DRTRA in Fermi
219  drtra=0.5d0
220  dytra=2.65d0
221  rdytra=rndm(v)*dytra
222  DO 1711 j=1,40
223  dyylam(j)=1.d-18
224  dyalam(j)=1.d-18
225  DO 1711 jj=1,9
226  DO 1711 jjj=1,41
227  yyltra(j,jj,jjj)=1.d-18
228  xyltra(j,jj,jjj)=-12.00+rdytra+j*dytra
229  1711 CONTINUE
230  DO 1712 j=1,40
231  aj=j
232  transa(j)=3.14d0*(2.d0*aj-1.d0)*drtra**2
233  1712 CONTINUE
234 C COMMON /WEWEVT/ IYWEW(NMXHKK),IRWEW(NMXHKK),IIYWEW(500),
235 C & IIRWEW(500),WEWEW(500),IDWEW1(500),IDWEW2(500),
236 C & IDWEW(500)
237  DO 1812 i=1,2600
238  iiywew(i)=0
239  iirwew(i)=0
240  wewew(i)=0.d0
241  idwew1(i)=0
242  idwew2(i)=0
243  idwew(i)=0
244  idweww(i)=0
245  1812 CONTINUE
246  DO 1813 i=nhkkh1,nhkk
247  iywew(i)=0
248  irwew(i)=0
249  iblock(i)=0
250  1813 CONTINUE
251 C RETURN
252 C
253 C-------------------------------------------------------------------
254 C
255  2 CONTINUE
256 C
257  DO 21 i=nhkkh1,nhkk
258  IF (isthkk(i).EQ.1)THEN
259 C j.r. temp for s-s
260  ptt=phkk(1,i)**2+phkk(2,i)**2+0.00001
261  pt=sqrt(ptt)+0.001
262  nhad=nhad+1
263  nrhkk=mcihad(idhkk(i))
264  IF (nrhkk.LE.0.OR.nrhkk.GT.210)THEN
265  nrhkk=1
266  ENDIF
267  nre=nrhkk
268  ichhkk=iich(nrhkk)
269  IF (nre.GT.160) nre=28
270  IF (nre.LT. 1) nre=28
271  nrex=nre
272  IF(nrex.GE.28)nrex=28
273  ni=indx(nrex)
274  nix=ni
275  IF (nrhkk.EQ.9)nix=12
276  IF (nrhkk.EQ.17.OR.nrhkk.EQ.22)nix=13
277  IF (nrhkk.LE.22.AND.nrhkk.GE.20)nix=14
278  IF (nrhkk.EQ.18.OR.nrhkk.EQ.100)nix=15
279  IF (nrhkk.LE.101.AND.nrhkk.GE.99)nix=16
280  IF (nrhkk.EQ.98)nix=17
281  IF (nrhkk.EQ.103)nix=18
282  IF (nrhkk.EQ.12.OR.nrhkk.EQ.19)nix=19
283  IF (nrhkk.EQ.24.OR.nrhkk.EQ.25)nix=19
284  IF(nre.EQ.28) ni=7
285  pt=sqrt(ptt)+0.001
286  amt=sqrt(ptt+phkk(5,i)**2)
287  yl=log((abs(phkk(3,i)+sqrt(phkk(3,i)**2+amt**2)))/amt+1.e-18)
288  yllps=log((abs(phkk(3,i)+sqrt(phkk(3,i)**2+ptt)))/pt
289  * +1.e-18)
290 C Transverse rapidity Densities
291  rtra=sqrt(vhkk(1,i)**2+vhkk(2,i)**2)*1.d12
292  irtra=rtra/drtra+1.d0
293  IF(irtra.LT.1)irtra=1
294  IF(irtra.GT.40)irtra=40
295  iytra=(yl+12.00-rdytra)/dytra+1.d0
296  IF(iytra.LT.1)iytra=1
297  IF(iytra.GT.40)iytra=40
298  iywew(i)=iytra
299  irwew(i)=irtra
300  IF(nix.EQ.4)THEN
301  yyltra(iytra,2,irtra)=yyltra(iytra,2,irtra)+
302  * 1.d0/transa(irtra)
303  ELSEIF(nix.EQ.3)THEN
304  yyltra(iytra,4,irtra)=yyltra(iytra,4,irtra)+
305  * 1.d0/transa(irtra)
306  ELSEIF(nix.EQ.1)THEN
307  yyltra(iytra,1,irtra)=yyltra(iytra,1,irtra)+
308  * 1.d0/transa(irtra)
309  ELSEIF(nix.EQ.8)THEN
310  yyltra(iytra,3,irtra)=yyltra(iytra,3,irtra)+
311  * 1.d0/transa(irtra)
312  ELSEIF(nix.EQ.13)THEN
313  yyltra(iytra,5,irtra)=yyltra(iytra,5,irtra)+
314  * 1.d0/transa(irtra)
315  ELSEIF(nix.EQ.12)THEN
316  yyltra(iytra,8,irtra)=yyltra(iytra,8,irtra)+
317  * 1.d0/transa(irtra)
318  ELSEIF(nix.EQ.2)THEN
319  yyltra(iytra,7,irtra)=yyltra(iytra,7,irtra)+
320  * 1.d0/transa(irtra)
321  ELSEIF(nix.EQ.15)THEN
322  yyltra(iytra,6,irtra)=yyltra(iytra,6,irtra)+
323  * 1.d0/transa(irtra)
324  ELSEIF(nix.EQ.23)THEN
325  yyltra(iytra,9,irtra)=yyltra(iytra,9,irtra)+
326  * 1.d0/transa(irtra)
327  ENDIF
328  ENDIF
329  21 CONTINUE
330 C RETURN
331 C
332 C--------------------------------------------------------------------
333 C
334  3 CONTINUE
335 C Normalization and Printing
336 C------------------------------------------------------------------
337 C------------------------------------------------------------------
338 C Transverse Rapidity Densities
339  DO 1133 i=1,40
340  DO 1133 ii=1,9
341  DO 1133 iii=1,40
342  yyltra(i,ii,iii)=yyltra(i,ii,iii)/(dytra)
343  yyltra(i,ii,41) =yyltra(i,ii,41) +
344  * yyltra(i,ii,iii)*transa(iii)
345  1133 CONTINUE
346 C transverse Rapidity Densities
347  IF(ipev.GE.2)THEN
348  WRITE(6,'(A)')' Transverse Rapidity Density of proton '
349  WRITE(6,37)xyltra(1,1,1),(transa(i),i=1,10)
350  37 FORMAT(f10.2,10e11.3)
351  DO 1136 j=1,40
352  WRITE(6,37)xyltra(j,1,1),(yyltra(j,1,i),i=1,9),yyltra(j,1,41)
353  1136 CONTINUE
354  WRITE(6,'(A)')' Transverse Rapidity Density of pi- '
355  WRITE(6,37)xyltra(1,1,1),(transa(i),i=1,10)
356  DO 1137 j=1,40
357  WRITE(6,37)xyltra(j,1,1),(yyltra(j,2,i),i=1,9),yyltra(j,2,41)
358  1137 CONTINUE
359  WRITE(6,'(A)')' Transverse Rapidity Density of aproton '
360  WRITE(6,37)xyltra(1,1,1),(transa(i),i=1,10)
361  DO 2236 j=1,40
362  WRITE(6,37)xyltra(j,1,1),(yyltra(j,3,i),i=1,9),yyltra(j,3,41)
363  2236 CONTINUE
364  WRITE(6,'(A)')' Transverse Rapidity Density of pi+ '
365  WRITE(6,37)xyltra(1,1,1),(transa(i),i=1,10)
366  DO 2237 j=1,40
367  WRITE(6,37)xyltra(j,1,1),(yyltra(j,4,i),i=1,9),yyltra(j,4,41)
368  2237 CONTINUE
369  WRITE(6,'(A)')' Transverse Rapidity Density of pi0 '
370  WRITE(6,37)xyltra(1,1,1),(transa(i),i=1,10)
371  DO 2238 j=1,40
372  WRITE(6,37)xyltra(j,1,1),(yyltra(j,9,i),i=1,9),yyltra(j,9,41)
373  2238 CONTINUE
374  ENDIF
375  iwew=0
376  DO 1134 iii=1,40
377  DO 1134 i=1,40
378 C pi- + p --> Lam + K0
379  dellla=
380  * (yyltra(i,1,iii)-yyltra(i,5,iii))*yyltra(i,2,iii)*
381  * transa(iii)*0.53d0
382  dellla=dellla*dytra
383 C pi+ + n --> Lam + K+
384  delllb=
385  * (yyltra(i,7,iii)-yyltra(i,5,iii))*yyltra(i,4,iii)*
386  * transa(iii)*0.53d0
387  delllb=delllb*dytra
388 C pi0 + n --> Lam + K0
389  delllc=
390  * (yyltra(i,7,iii)-yyltra(i,5,iii))*yyltra(i,9,iii)*
391  * transa(iii)*0.53d0
392  delllc=delllc*dytra
393 C pi0 + p --> Lam + K+
394  dellld=
395  * (yyltra(i,1,iii)-yyltra(i,5,iii))*yyltra(i,9,iii)*
396  * transa(iii)*0.53d0
397  dellld=dellld*dytra
398 C 3.11.95 Add 0.15* to reduce Lambdas
399 C * 0.15*
400  delllo=
401  * (yyltra(i,1,iii)+yyltra(i,7,iii)-1.6*yyltra(i,5,iii))*
402  * transa(iii)/3.2
403  delllo=delllo*dytra
404  IF(dellla.GT.delllo)dellla=delllo
405  IF(delllb.GT.delllo)delllb=delllo
406  IF(delllc.GT.delllo)delllc=delllo
407  IF(dellld.GT.delllo)dellld=delllo
408  IF(delllo.LE.0.d0)dellla=0.d0
409  IF(delllo.LE.0.d0)delllb=0.d0
410  IF(delllo.LE.0.d0)delllc=0.d0
411  IF(delllo.LE.0.d0)dellld=0.d0
412  IF(dellla.LE.0.d0)dellla=0.d0
413  IF(delllb.LE.0.d0)delllb=0.d0
414  IF(delllc.LE.0.d0)delllc=0.d0
415  IF(dellld.LE.0.d0)dellld=0.d0
416 C pi- + p --> Lam + K0
417  IF(dellla.GT.tiny)THEN
418  idelll=dellla
419  tellla=dellla-idelll
420  IF(rndm(v).LE.tellla)idelll=idelll+1
421  dellla=idelll
422  IF (dellla.LT.tiny)go to 2233
423  DO 556 kkk=1,idelll
424  iwew=iwew+1
425  wewew(iwew)=dellla
426  idwew(iwew)=17
427  idweww(iwew)=24
428  idwew1(iwew)=1
429  idwew2(iwew)=14
430  IF(ipev.GE.2)THEN
431  WRITE (6,'(A,3F10.3,7I5)')' LAMBDA ',
432  * dellla,yyltra(i,1,iii),yyltra(i,2,iii),i,iii,
433  * idwew(iwew),idweww(iwew),idwew1(iwew),idwew2(iwew),
434  * iwew
435  ENDIF
436  dyylam(i)=dyylam(i)+ dellla/dytra
437  iiywew(iwew)=i
438  iirwew(iwew)=iii
439  556 CONTINUE
440 C IYWEW(IWEW)=I
441 C IRWEW(IWEW)=III
442  2233 CONTINUE
443  ENDIF
444 C pi+ + n --> Lam + K+
445  IF(delllb.GT.tiny)THEN
446  idelll=delllb
447  telllb=delllb-idelll
448  IF(rndm(v).LE.telllb)idelll=idelll+1
449  delllb=idelll
450  IF (delllb.LT.tiny)go to 2243
451  DO 557 kkk=1,idelll
452  iwew=iwew+1
453  wewew(iwew)=delllb
454  idwew(iwew)=17
455  idweww(iwew)=15
456  idwew1(iwew)=8
457  idwew2(iwew)=13
458  IF(ipev.GE.2)THEN
459  WRITE (6,'(A,3F10.3,7I5)')' LAMBDA ',
460  * delllb,yyltra(i,1,iii),yyltra(i,2,iii),i,iii,
461  * idwew(iwew),idweww(iwew),idwew1(iwew),idwew2(iwew),
462  * iwew
463  ENDIF
464  dyylam(i)=dyylam(i)+ delllb/dytra
465  iiywew(iwew)=i
466  iirwew(iwew)=iii
467  557 CONTINUE
468 C IYWEW(IWEW)=I
469 C IRWEW(IWEW)=III
470  2243 CONTINUE
471  ENDIF
472 C pi0 + n --> Lam + K0
473  IF(delllc.GT.tiny)THEN
474  idelll=delllc
475  telllc=delllc-idelll
476  IF(rndm(v).LE.telllc)idelll=idelll+1
477  delllc=idelll
478  IF (delllc.LT.tiny)go to 2244
479  DO 558 kkk=1,idelll
480  iwew=iwew+1
481  wewew(iwew)=delllc
482  idwew(iwew)=17
483  idweww(iwew)=24
484  idwew1(iwew)=8
485  idwew2(iwew)=23
486  IF(ipev.GE.2)THEN
487  WRITE (6,'(A,3F10.3,7I5)')' LAMBDA ',
488  * delllc,yyltra(i,1,iii),yyltra(i,2,iii),i,iii,
489  * idwew(iwew),idweww(iwew),idwew1(iwew),idwew2(iwew),
490  * iwew
491  ENDIF
492  dyylam(i)=dyylam(i)+ delllc/dytra
493  iiywew(iwew)=i
494  iirwew(iwew)=iii
495  558 CONTINUE
496 C IYWEW(IWEW)=I
497 C IRWEW(IWEW)=III
498  2244 CONTINUE
499  ENDIF
500 C pi0 + p --> Lam + K+
501  IF(dellld.GT.tiny)THEN
502  idelll=dellld
503  tellld=dellld-idelll
504  IF(rndm(v).LE.tellld)idelll=idelll+1
505  dellld=idelll
506  IF (dellld.LT.tiny)go to 2245
507  DO 559 kkk=1,idelll
508  iwew=iwew+1
509  wewew(iwew)=dellld
510  idwew(iwew)=17
511  idweww(iwew)=15
512  idwew1(iwew)=1
513  idwew2(iwew)=23
514  IF(ipev.GE.2)THEN
515  WRITE (6,'(A,3F10.3,7I5)')' LAMBDA ',
516  * dellld,yyltra(i,1,iii),yyltra(i,2,iii),i,iii,
517  * idwew(iwew),idweww(iwew),idwew1(iwew),idwew2(iwew),
518  * iwew
519  ENDIF
520  dyylam(i)=dyylam(i)+ dellld/dytra
521  iiywew(iwew)=i
522  iirwew(iwew)=iii
523  559 CONTINUE
524 C IYWEW(IWEW)=I
525 C IRWEW(IWEW)=III
526  2245 CONTINUE
527  ENDIF
528 C COMMON /WEWEVT/ IYWEW(NMXHKK),IRWEW(NMXHKK),IIYWEW(500),
529 C & IIRWEW(500),WEWEW(500),IDWEW1(500),IDWEW2(500),
530 C & IDWEW(500)
531 C pi+ + pb --> aLam + aK0
532  delala=
533  * (yyltra(i,3,iii)-yyltra(i,6,iii))*yyltra(i,4,iii)*
534  * transa(iii)*0.53d0
535  delala=delala*dytra
536 C pi- + nb --> aLam + aK-
537  delalb=
538  * (yyltra(i,8,iii)-yyltra(i,6,iii))*yyltra(i,2,iii)*
539  * transa(iii)*0.53d0
540  delalb=delalb*dytra
541 C 3.11.95 Add 0.15* to reduce Lambdas
542 C * 0.15*
543  delalo=
544  * (yyltra(i,3,iii)+yyltra(i,8,iii)-1.6*yyltra(i,6,iii))*
545  * transa(iii)/3.2
546  delalo=delalo*dytra
547  IF(delala.GT.delalo)delala=delalo
548  IF(delalb.GT.delalo)delalb=delalo
549  IF(delalo.LE.0.d0)delala=0.d0
550  IF(delalo.LE.0.d0)delalb=0.d0
551  IF(delala.LE.0.d0)delala=0.d0
552  IF(delalb.LE.0.d0)delalb=0.d0
553 C pi+ + pb --> aLam + aK0
554  IF(delala.GT.tiny)THEN
555  idelal=delala
556  telala=delala-idelal
557  IF(rndm(v).LE.telala)idelal=idelal+1
558  delala=idelal
559  IF (delala.LT.tiny)go to 2234
560  DO 554 kkk=1,idelal
561  iwew=iwew+1
562  wewew(iwew)=delala
563  idwew(iwew)=18
564  idweww(iwew)=25
565  idwew1(iwew)=2
566  idwew2(iwew)=13
567  IF(ipev.GE.2)THEN
568  WRITE (6,'(A,3F10.3,7I5)')' ALAMBDA ',
569  * delala,yyltra(i,3,iii),yyltra(i,6,iii),i,iii,
570  * idwew(iwew),idweww(iwew),idwew1(iwew),idwew2(iwew),
571  * iwew
572  ENDIF
573  dyalam(i)=dyalam(i)+ delala/dytra
574  iiywew(iwew)=i
575  iirwew(iwew)=iii
576  554 CONTINUE
577 C IYWEW(IWEW)=I
578 C IRWEW(IWEW)=III
579  2234 CONTINUE
580  ENDIF
581 C pi- + nb --> aLam + aK-
582  IF(delalb.GT.tiny)THEN
583  idelal=delalb
584  telalb=delalb-idelal
585  IF(rndm(v).LE.telalb)idelal=idelal+1
586  delalb=idelal
587  IF (delalb.LT.tiny)go to 2247
588  DO 555 kkk=1,idelal
589  iwew=iwew+1
590  wewew(iwew)=delalb
591  idwew(iwew)=18
592  idweww(iwew)=16
593  idwew1(iwew)=9
594  idwew2(iwew)=14
595  IF(ipev.GE.2)THEN
596  WRITE (6,'(A,3F10.3,7I5)')' ALAMBDA ',
597  * delalb,yyltra(i,3,iii),yyltra(i,6,iii),i,iii,
598  * idwew(iwew),idweww(iwew),idwew1(iwew),idwew2(iwew),
599  * iwew
600  ENDIF
601  dyalam(i)=dyalam(i)+ delalb/dytra
602  iiywew(iwew)=i
603  iirwew(iwew)=iii
604  555 CONTINUE
605 C IYWEW(IWEW)=I
606 C IRWEW(IWEW)=III
607  2247 CONTINUE
608  ENDIF
609  1134 CONTINUE
610  iwewma=iwew
611  ddlamn=0.d0
612  dalamn=0.d0
613  DO 1139 j=1,40
614  ddlamn=ddlamn+dyylam(j)*dytra
615  dalamn=dalamn+dyalam(j)*dytra
616  1139 CONTINUE
617  IF(ipev.GE.2)THEN
618  WRITE(6,'(A,I10,F10.3)')' IWEWMA,Extra Rap. Den. of Lam ,DLAM= '
619  * ,iwewma ,ddlamn
620  WRITE(6,'(A,I10,F10.3)')' IWEWMA,Extra Rap. Den. ofALam ,DLAM= '
621  * ,iwewma ,dalamn
622  WRITE(6,37)xyltra(1,1,1),(transa(i),i=1,10)
623  DO 1138 j=1,40
624  WRITE(6,37)xyltra(j,1,1),dyylam(j),dyalam(j)
625  1138 CONTINUE
626  ENDIF
627  i=0
628 C Reduce IWEWMA forPb-Pb runs
629  iwewma=iwewma
630  iloma=0
631  iloml=0
632  itoma=0
633  itoml=0
634  DO 20 ii=1,iwewma
635  id1=mpdgha(idwew1(ii))
636  id2=mpdgha(idwew2(ii))
637  id3=mpdgha(idwew(ii))
638  id4=mpdgha(idweww(ii))
639 C WRITE(6,'(5I5)')II,ID1,ID2,ID3,ID4
640  ihkk1=0
641  ihkk2=0
642  kk1=0
643  kk2=0
644  DO 211 jj=nhkkh1,nhkk
645  IF(iblock(jj).EQ.0)THEN
646  IF((idhkk(jj).EQ.id1).AND.(iywew(jj).EQ.iiywew(ii))
647  * .AND.(irwew(jj).EQ.iirwew(ii)).AND.(kk1.EQ.0))THEN
648  ihkk1=jj
649  jj1=jj
650  kk1=1
651  ENDIF
652  IF((idhkk(jj).EQ.id2).AND.(iywew(jj).EQ.iiywew(ii))
653  * .AND.(irwew(jj).EQ.iirwew(ii)).AND.(kk2.EQ.0))THEN
654  ihkk2=jj
655  jj2=jj
656  kk2=1
657  ENDIF
658  IF(ihkk1.EQ.0)go to 211
659  IF(ihkk2.EQ.0)go to 211
660  IF(jmohkk(1,ihkk1).EQ.jmohkk(1,ihkk2))THEN
661  ihkk2=0
662  jj2=0
663  kk2=0
664  go to 211
665  ENDIF
666  CALL pinkla(idwew1(ii),idwew2(ii),idwew(ii),idweww(ii),
667  & ihkk1,ihkk2,iloml,iloma,itoml,itoma,irej)
668  IF(ipev.GE.2)
669  * WRITE(6,'(A,5I10)')' IWEWMA,ILOML,ILOMA,ITOML,ITOMA '
670  * ,iwewma ,iloml,iloma,itoml,itoma
671  IF(irej.EQ.0)THEN
672  IF(ipev.GE.2)THEN
673  WRITE(6,'(A,3I7)')' Extra Positions ',ihkk1,ihkk2,ii
674  ENDIF
675  iblock(jj1)=1
676  iblock(jj2)=1
677  kk1=0
678  kk2=0
679  jj1=0
680  jj2=0
681  ihkk1=0
682  ihkk2=0
683  go to 20
684  ENDIF
685  IF(jj1.GE.1)iblock(jj1)=1
686  ihkk1=0
687  jj1=0
688  kk1=0
689  IF(jj2.GE.1)iblock(jj2)=1
690  ihkk2=0
691  jj2=0
692  kk2=0
693  ENDIF
694  211 CONTINUE
695  20 CONTINUE
696  IF(ipev.GE.2)
697  *WRITE(6,'(A,5I10)')' IWEWMA,ILOML,ILOMA,ITOML,ITOMA '
698  * ,iwewma ,iloml,iloma,itoml,itoma
699 
700  RETURN
701  END
702  SUBROUTINE pinkla(I1,I2,I3,I4,IHKK1,IHKK2,ILOML,ILOMA,
703  *itoml,itoma,irej)
704 C
705  IMPLICIT DOUBLE PRECISION (a-h,o-z)
706  SAVE
707 C--------------------------------------------------------------
708 C
709 C I1 + I2 ---> I3 + I4 final state interaction
710 C
711 C--------------------------------------------------------------
712 *KEEP,HKKEVT.
713 c INCLUDE (HKKEVT)
714  parameter(tiny= 1.d-5)
715  parameter(nmxhkk= 89998)
716 c PARAMETER (NMXHKK=25000)
717  COMMON /hkkevt/ nhkk,nevhkk,isthkk(nmxhkk),idhkk(nmxhkk), jmohkk
718  +(2,nmxhkk),jdahkk(2,nmxhkk), phkk(5,nmxhkk),vhkk(4,nmxhkk),whkk
719  +(4,nmxhkk)
720  COMMON /extevt/ idres(nmxhkk),idxres(nmxhkk),nobam(nmxhkk),
721  & idbam(nmxhkk),idch(nmxhkk),npoint(10)
722  COMMON /wewevt/ iywew(nmxhkk),irwew(nmxhkk),iiywew(2600),
723  & iirwew(2600),wewew(2600),idwew1(2600),
724  & idwew2(2600),idwew(2600),iblock(nmxhkk),
725  & idweww(2600)
726 C
727 C
728 C WHKK(4,NMXHKK) GIVES POSITIONS AND TIMES IN
729 C PROJECTILE FRAME, THE CHAINS ARE CREATED ON
730 C THE POSITIONS OF THE PROJECTILE NUCLEONS
731 C IN THE PROJECTILE FRAME (TARGET NUCLEONS IN
732 C TARGET FRAME) BOTH POSITIONS ARE THREFORE NOT
733 C COMPLETELY CONSISTENT. THE TIMES IN THE
734 C PROJECTILE FRAME HOWEVER ARE OBTAINED BY
735 C LORENTZ TRANSFORMING FROM THE LAB SYSTEM.
736 C
737 C Based on the proposed standard COMMON block (Sjostrand Memo 17.3,89)
738 C
739 C NMXHKK: maximum numbers of entries (partons/particles) that can be
740 C stored in the commonblock.
741 C
742 C NHKK: the actual number of entries stored in current event. These are
743 C found in the first NHKK positions of the respective arrays below.
744 C Index IHKK, 1 <= IHKK <= NHKK, is below used to denote a given
745 C entry.
746 C
747 C ISTHKK(IHKK): status code for entry IHKK, with following meanings:
748 C = 0 : null entry.
749 C = 1 : an existing entry, which has not decayed or fragmented.
750 C This is the main class of entries which represents the
751 C "final state" given by the generator.
752 C = 2 : an entry which has decayed or fragmented and therefore
753 C is not appearing in the final state, but is retained for
754 C event history information.
755 C = 3 : a documentation line, defined separately from the event
756 C history. (incoming reacting
757 C particles, etc.)
758 C = 4 - 10 : undefined, but reserved for future standards.
759 C = 11 - 20 : at the disposal of each model builder for constructs
760 C specific to his program, but equivalent to a null line in the
761 C context of any other program. One example is the cone defining
762 C vector of HERWIG, another cluster or event axes of the JETSET
763 C analysis routines.
764 C = 21 - : at the disposal of users, in particular for event tracking
765 C in the detector.
766 C
767 C IDHKK(IHKK) : particle identity, according to the Particle Data Group
768 C standard.
769 C
770 C JMOHKK(1,IHKK) : pointer to the position where the mother is stored.
771 C The value is 0 for initial entries.
772 C
773 C JMOHKK(2,IHKK) : pointer to position of second mother. Normally only
774 C one mother exist, in which case the value 0 is used. In cluster
775 C fragmentation models, the two mothers would correspond to the q
776 C and qbar which join to form a cluster. In string fragmentation,
777 C the two mothers of a particle produced in the fragmentation would
778 C be the two endpoints of the string (with the range in between
779 C implied).
780 C
781 C JDAHKK(1,IHKK) : pointer to the position of the first daughter. If an
782 C entry has not decayed, this is 0.
783 C
784 C JDAHKK(2,IHKK) : pointer to the position of the last daughter. If an
785 C entry has not decayed, this is 0. It is assumed that the daughters
786 C of a particle (or cluster or string) are stored sequentially, so
787 C that the whole range JDAHKK(1,IHKK) - JDAHKK(2,IHKK) contains
788 C daughters. Even in cases where only one daughter is defined (e.g.
789 C K0 -> K0S) both values should be defined, to make for a uniform
790 C approach in terms of loop constructions.
791 C
792 C PHKK(1,IHKK) : momentum in the x direction, in GeV/c.
793 C
794 C PHKK(2,IHKK) : momentum in the y direction, in GeV/c.
795 C
796 C PHKK(3,IHKK) : momentum in the z direction, in GeV/c.
797 C
798 C PHKK(4,IHKK) : energy, in GeV.
799 C
800 C PHKK(5,IHKK) : mass, in GeV/c**2. For spacelike partons, it is allowed
801 C to use a negative mass, according to PHKK(5,IHKK) = -sqrt(-m**2).
802 C
803 C VHKK(1,IHKK) : production vertex x position, in mm.
804 C
805 C VHKK(2,IHKK) : production vertex y position, in mm.
806 C
807 C VHKK(3,IHKK) : production vertex z position, in mm.
808 C
809 C VHKK(4,IHKK) : production time, in mm/c (= 3.33*10**(-12) s).
810 C********************************************************************
811 *KEEP,DPAR.
812 C /DPAR/ CONTAINS PARTICLE PROPERTIES
813 C ANAME = LITERAL NAME OF THE PARTICLE
814 C AAM = PARTICLE MASS IN GEV
815 C GA = DECAY WIDTH
816 C TAU = LIFE TIME OF INSTABLE PARTICLES
817 C IICH = ELECTRIC CHARGE OF THE PARTICLE
818 C IIBAR = BARYON NUMBER
819 C K1,K1 = BIGIN AND END OF DECAY CHANNELS OF PARTICLE
820 C
821  CHARACTER*8 aname
822  COMMON /dpar/ aname(210),aam(210),ga(210),tau(210),iich(210),
823  +iibar(210),k1(210),k2(210)
824  COMMON /dprin/ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
825 C------------------
826  IF(i3.EQ.17)itoml=itoml+1
827  IF(i3.EQ.18)itoma=itoma+1
828 C Check Inv Mass of two particle system
829  ama12 = sqrt((phkk(4,ihkk1)+phkk(4,ihkk2))**2
830  + -(phkk(1,ihkk1)+phkk(1,ihkk2))**2
831  + -(phkk(2,ihkk1)+phkk(2,ihkk2))**2
832  + -(phkk(3,ihkk1)+phkk(3,ihkk2))**2)
833  amamin = aam(i3)+aam(i4)
834  irej=1
835  IF(ama12.GE.amamin+0.05d0)THEN
836  IF(ipev.GE.2)THEN
837  WRITE(6,'(A,2F10.3)')' AMAMIN,AMA12 ',amamin,ama12
838  WRITE(6,'(A,I5,A,I5,A,I5,A,I5)')' Reaction ',
839  + i1,' + ',i2,' ---> ',i3,' + ',i4
840  ENDIF
841 C Lorentz parameters of IHKK1,IHKK2 system
842  e12 = phkk(4,ihkk1)+phkk(4,ihkk2)
843  px12 = phkk(1,ihkk1)+phkk(1,ihkk2)
844  py12 = phkk(2,ihkk1)+phkk(2,ihkk2)
845  pz12 = phkk(3,ihkk1)+phkk(3,ihkk2)
846  g12 = e12/ama12
847  bgx12 = px12/ama12
848  bgy12 = py12/ama12
849  bgz12 = pz12/ama12
850  IF(ipev.GE.2)
851  +WRITE(6,'(A,4E12.3)')' G12,BGX12,BGY12,BGZ12 ',
852  +g12,bgx12,bgy12,bgz12
853 C Decay AMA12 into I3 + I4
854  e3 = (ama12**2-aam(i4)**2+aam(i3)**2)/(2.*ama12)
855  e4 = (ama12**2-aam(i3)**2+aam(i4)**2)/(2.*ama12)
856  p3 = sqrt(e3**2-aam(i3)**2)
857  p4 = sqrt(e4**2-aam(i4)**2)
858  IF(ipev.GE.2)
859  +WRITE(6,'(A,4E12.3)')' E3,E4,P3,P4 ',
860  +e3,e4,p3,p4
861 C Uniform rotation
862  CALL dpoli(polc,pols)
863  CALL dsfecf(sfe,cfe)
864  IF(phkk(3,ihkk1).GE.0.d0)THEN
865  pols=0.d0
866  polc=1.d0
867  ELSEIF(phkk(3,ihkk1).LE.0.d0)THEN
868  pols=0.d0
869  polc=-1.d0
870  ENDIF
871 C WRITE(6,*)' POLS,POLC ',POLS,POLC
872  cx=pols*cfe
873  cy=pols*sfe
874  cz=polc
875  px3=cx*p3
876  py3=cy*p3
877  pz3=cz*p3
878  px4=-cx*p4
879  py4=-cy*p4
880  pz4=-cz*p4
881  IF(ipev.GE.2)
882 C +WRITE(6,'(A,8E12.3)')' 4 mom cms 3 4 ',
883 C +PX3,PY3,PZ3,E3,PX4,PY4,PZ4,E4
884  +WRITE(6,'(A,4E12.3)')' 4 mom cms 3 4 ',
885  +pz3,e3,pz4,e4
886 C Lorentz transformation back
887  CALL daltra(g12,bgx12,bgy12,bgz12,px3,py3,pz3,e3,pc3,pcx3,
888  + pcy3,pcz3,ec3)
889  CALL daltra(g12,bgx12,bgy12,bgz12,px4,py4,pz4,e4,pc4,pcx4,
890  + pcy4,pcz4,ec4)
891  IF(ipev.GE.2)
892 C +WRITE(6,'(A,8E12.3)')' 4 mom 1 2 ',
893 C +PHKK(1,IHKK1),PHKK(2,IHKK1),PHKK(3,IHKK1),PHKK(4,IHKK1),
894 C +PHKK(1,IHKK2),PHKK(2,IHKK2),PHKK(3,IHKK2),PHKK(4,IHKK2)
895  +WRITE(6,'(A,4E12.3)')' 4 mom 1 2 ',
896  +phkk(3,ihkk1),phkk(4,ihkk1),
897  +phkk(3,ihkk2),phkk(4,ihkk2)
898  IF(ipev.GE.2)
899 C +WRITE(6,'(A,8E12.3)')' 4 mom 3 4 ',
900 C +PCX3,PCY3,PCZ3,EC3,PCX4,PCY4,PCZ4,EC4
901  +WRITE(6,'(A,4E12.3)')' 4 mom 3 4 ',
902  +pcz3,ec3,pcz4,ec4
903 C---------------------------------------------------
904  irej=0
905 C WRITE(6,*)' I1+I2-->I3+I4 ',I1,I2,I3,I4
906 C Put particles I3 and I4 into COMMON Block
907  isthkk(ihkk1)=512
908  isthkk(ihkk2)=512
909  nhkk=nhkk+1
910  isthkk(nhkk)=1
911  IF((i3.EQ.24).OR.(i3.EQ.25))THEN
912  i3=12
913  IF(rndm(vv).LE.0.5d0)i3=19
914  ENDIF
915  idhkk(nhkk)=mpdgha(i3)
916  jmohkk(1,nhkk)=ihkk1
917  jmohkk(2,nhkk)=ihkk2
918  jdahkk(1,ihkk1)=nhkk
919  jdahkk(1,ihkk2)=nhkk
920  phkk(1,nhkk)=pcx3
921  phkk(2,nhkk)=pcy3
922  phkk(3,nhkk)=pcz3
923  phkk(4,nhkk)=ec3
924  phkk(5,nhkk)=aam(i3)
925  DO 11 iiii=1,4
926  vhkk(iiii,nhkk)=vhkk(iiii,ihkk1)
927  whkk(iiii,nhkk)=whkk(iiii,ihkk1)
928  11 CONTINUE
929  nhkk=nhkk+1
930  isthkk(nhkk)=1
931  IF((i4.EQ.24).OR.(i4.EQ.25))THEN
932  i4=12
933  IF(rndm(vv).LE.0.5d0)i4=19
934  ENDIF
935  idhkk(nhkk)=mpdgha(i4)
936  jmohkk(1,nhkk)=ihkk1
937  jmohkk(2,nhkk)=ihkk2
938  jdahkk(2,ihkk1)=nhkk
939  jdahkk(2,ihkk2)=nhkk
940  phkk(1,nhkk)=pcx4
941  phkk(2,nhkk)=pcy4
942  phkk(3,nhkk)=pcz4
943  phkk(4,nhkk)=ec4
944  phkk(5,nhkk)=aam(i4)
945  DO 12 iiii=1,4
946  vhkk(iiii,nhkk)=vhkk(iiii,ihkk2)
947  whkk(iiii,nhkk)=whkk(iiii,ihkk2)
948  12 CONTINUE
949  ELSE
950  IF(i3.EQ.17)iloml=iloml+1
951  IF(i3.EQ.18)iloma=iloma+1
952  ENDIF
953  RETURN
954  END
955 
956 C++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
957 C
958  SUBROUTINE hadjck(NHAD,AMCH,PPR,PTA,GAM,BGX,BGY,BGZ, IFB1,IFB2,
959  +ifb3,ifb4,i1,i2,nobam,nnch,norig,irej)
960  IMPLICIT DOUBLE PRECISION (a-h,o-z)
961  SAVE
962 C
963 C HADJCK Fragmentation of vv,vs and sv chains CK method
964 C j.r.6/96
965 C
966 C HADJET DOES ALL THE NECESSARY ROTATIONS AND LORENTZ TRANSFORMS AND
967 C CALL CALBAM (BAMJET)
968 C
969 C NHAD = NUMBER OF HADRONS CREATED (OUTPUT) = IHAD (CALBAM)
970 C******** PRODUCED PARTICLES IN COMMON /CMSRES/
971 C NOTE: NOW IN /FINPAR/ HJM 21/06/90
972 C AMCH = INVARIANT MASS OF JET (INPUT)
973 C PPR = 4-MOMENTUM OF FORWARD GOING PARTON (PROJECTILE)(INPUT)
974 C PTA = 4-MOMENTUM OF BACKWARD GOING PARTON (TARGET)(INPUT)
975 C GAM,BGX,BGY,BGZ = LORENTZ PARAMETERS OF JET CMS RELATIVE TO
976 C COLLISION CMS (INPUT)
977 C
978 C--------------------------------------------------------------------
979 C CALBAM(NNCH,I1,I2,IFB1,IFB2,IFB3,IFB4,AMCH,NOBAM,IHAD)
980 C SAMPLING OF Q-AQ, Q-QQ, QQ-AQQ CHAINS
981 C USING BAMJET(IHAD,IFB1,IFB2,IFB3,IFB4,AMCH,NOBAM)-----FOR NNCH=0
982 C OR PARJET(IHAD,ICH1=I1 OR I2)------FOR NNCH=-1 OR +1
983 C-------------------------------------------------------------------
984 C IHAD : NUMBER OF PRODUCED PARTICLES
985 C NNCH : CALL BAMJET FOR NNCH=0
986 C CALL PARJET FOR NNCH=+1 ICH1=I1
987 C FOR NNCH=-1 ICH1=I2
988 C jet not existing for NNCH=+/-99, i.e. IHAD=0
989 C PRODUCED PARTICLES IN CHAIN REST FRAME ARE IN COMMON /FINPAR/
990 C AMCH : INVARIANT MASS OF CHAIN (BAMJET)
991 C
992 C NOBAM : = 3 QUARK-ANTIQUARK JET QUARK FLAVORS : IFB1,IFB2
993 C OR ANTIQUARK-QUARK JET IN ANY ORDER
994 C
995 C = 4 QUARK-DIQUARK JET, FLAVORS : QU : IFB1, DIQU :IFB2,IFB
996 C OR ANTIQUARK-ANTIDIQUARK JET
997 C
998 C
999 C = 5 DIQUARK-ANTIDIQUARK JET
1000 C OR ANTIDIQUARK-DIQUARK JET
1001 C FLAVORS : DIQU : IFB1,IFB2, ANTIDIQU : IFB3,IFB4
1002 C IN ANY ORDER
1003 C
1004 C = 6 DIQUARK-QUARK JET, FLAVORS : DIQU : IFB1,IFB2 QU: IFB
1005 C OR ANTIDIQUARK-ANTIQUARK JET
1006 C
1007 C IFBI : FLAVORS : 1,2,3,4 = U,D,S,C 7,8,9,10 = AU,AD,AS,AC
1008 C
1009 C I1,I2 : NUMBER LABEL OF A HADRON CREATED BY PARJET
1010 C
1011 C NORMALLY IN BAMJET THE QUARKS MOVE FORWARD (POSITIVE Z-DIRECTION)
1012 C THE QUARK FLAVORS ARE FIRST GIVEN
1013 C CALBAM ALLOWS EITHER THE QUARK OR ANTIQUARK (DIQU) TO MOVE FORWARD
1014 C THE FORWARD GOING FLAVORS ARE GIVEN FIRST
1015 C
1016 C =====================================================================
1017 *KEEP,DFINPA.
1018  CHARACTER*8 anf,anff
1019  parameter(nfimax=249)
1020  COMMON /dfinpa/ anf(nfimax),pxf(nfimax),pyf(nfimax),pzf(nfimax),
1021  +hef(nfimax),amf(nfimax), ichf(nfimax),ibarf(nfimax),nref(nfimax)
1022  COMMON /dfinpz/iormo(nfimax),idaug1(nfimax),idaug2(nfimax),
1023  * istath(nfimax)
1024  dimension anff(nfimax),pxff(nfimax),pyff(nfimax),pzff(nfimax),
1025  +heff(nfimax),amff(nfimax),
1026  *ichff(nfimax),ibarff(nfimax),nreff(nfimax)
1027  CHARACTER*80 title
1028  CHARACTER*8 projty,targty
1029 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
1030 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
1031  COMMON /user1/title,projty,targty
1032  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
1033 *KEEP,DPRIN.
1034  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
1035 *KEND.
1036 *--------------------------------------------------------------------
1037 *-------------------------------------------------------------------
1038  COMMON /jspart/pxp(1000),pyp(1000),pzp(1000),hepp(1000),nnnp
1039  COMMON /jspar/pxj(1000),pyj(1000),pzj(1000),hej(1000),nnnpj
1040 ************************************************************************
1041 ************************************************************************
1042  common/ifragm/ifrag
1043  dimension ppr(4),pta(4)
1044 *KEEP,XSEADI.
1045  COMMON /xseadi/ xseacu,unon,unom,unosea,cvq,cdq,csea,ssmima,
1046  +ssmimq,vvmthr
1047  dimension ppp1(4),ppp2(4),ppp3(4),ppp4(4)
1048 C----------------------------------------------------------------------
1049 C
1050 C WRITE(6,'(A,3I5)')' Jet0 IFB1,IFB2,IFB3 ',IFB1,IFB2,IFB3
1051  irej=0
1052  IF(ipco.GE.3) THEN
1053  WRITE(6,'(4E12.4)') gam,bgx,bgy,bgz,ppr
1054  WRITE(6,1000) nhad,amch,pta, ifb1,ifb2,ifb3,ifb4,i1,i2,nobam,
1055  + nnch,norig
1056  1000 FORMAT(10x,i10,5f10.3/10x,9i10)
1057  ENDIF
1058  IF(abs(nnch).EQ.99) THEN
1059  nhad=0
1060 C IPCO=0
1061  RETURN
1062  ENDIF
1063 C=================================================================
1064 C=================================================================
1065 C Different options:
1066  IF(nobam.EQ.6)THEN
1067 C NOBAM=6: Diquark--Quark Jet
1068 C
1069 C Work out approximate x-Values of diquark and quark
1070  xdiq=ppr(4)*2.d0/cmener
1071  xqua=pta(4)*2.d0/cmener
1072 C WRITE(6,'(A,2E12.4)')' HADJCK:XDIQ,XQUA ',XDIQ,XQUA
1073 C Select valence quark x for one of the diquark quarks
1074 C But only if XDIQ > 1.5*XQUA drop this condition!
1075 C IF(XDIQ.LE.1.5D0*XQUA)THEN
1076 C IREJ=1
1077 C RETURN
1078 C ENDIF
1079 C Select x values of sea quark pair
1080  icou=0
1081  1234 CONTINUE
1082  icou=icou+1
1083  IF(icou.GE.100)THEN
1084  irej=1
1085  RETURN
1086  ENDIF
1087  CALL xseapa(cmener,xqua/2.d0,isq,isaq,xsq,xsaq,irej)
1088  IF(xsaq.GE.xqua/2.d0)go to 1234
1089 C WRITE(6,'(A,2E12.4)')' HADJCK:XSQ,XSAQ ',XSQ,XSAQ
1090 * projectile valence quarks xpvqi
1091  xvthro=cvq/cmener
1092  ivthr=0
1093  3465 CONTINUE
1094  IF(ivthr.EQ.50)THEN
1095  irej=1
1096  WRITE(6,*)' HADJSE 3465 reject IVTHR 50'
1097  RETURN
1098  ENDIF
1099  ivthr=ivthr+1
1100  xvthr=xvthro/(51-ivthr)
1101  unoprv=unon
1102  80 CONTINUE
1103  IF(xvthr.GT.0.05)THEN
1104  IF(xvthr.GT.0.66d0*xdiq)THEN
1105  irej=1
1106  RETURN
1107  ENDIF
1108  xpvqi=betrej(0.5,unoprv,xvthr,0.66d0*xdiq)
1109  ELSE
1110  n90=0
1111  90 CONTINUE
1112  n90=n90+1
1113  IF(n90.GE.20)THEN
1114  irej=1
1115  RETURN
1116  ENDIF
1117  xpvqi=dbetar(0.5,unoprv)
1118  IF ((xpvqi.LT.xvthr).OR.(xpvqi.GT.0.66d0*xdiq))
1119  * goto 90
1120  ENDIF
1121  xdiqq=xdiq-xpvqi
1122 C WRITE(6,'(A,2E12.4)')' HADJCK:XDIQQ,XPVQI ',XDIQQ,XPVQI
1123 C
1124 C We form two strings(1-2) q-aq and (3-4) qq-q
1125 C 1 q with XDIQ-XPVQI-XSQ
1126 C 1 q with XDIQ-XPVQI (j.r.5.5.96)
1127 C 2 aq with XSAQ
1128 C 3 qq with XPVQI+XSQ
1129 C 3 qq with XPVQI (j.r.5.5.96)
1130 C 4 q with XQUA-XSAQ
1131 C with 4-momenta PPP1,PPP2,PPP3,PPP4
1132  DO 2345 i=1,4
1133 C PPP1(I)=PPR(I)*(XDIQ-XPVQI-XSQ)/XDIQ
1134  ppp1(i)=ppr(i)*(xdiq-xpvqi)/xdiq
1135 C PPP3(I)=PPR(I)*(XPVQI+XSQ)/XDIQ
1136  ppp3(i)=ppr(i)*(xpvqi)/xdiq
1137  ppp2(i)=pta(i)*xsaq/xqua
1138  ppp4(i)=pta(i)*(xqua-xsaq)/xqua
1139  2345 CONTINUE
1140  2346 FORMAT(a,4e12.4)
1141 C WRITE(6,2346)' PPR ',PPR
1142 C WRITE(6,2346)' PPP1 ',PPP1
1143 C WRITE(6,2346)' PPP3 ',PPP3
1144 C WRITE(6,2346)' PTA ',PTA
1145 C WRITE(6,2346)' PPP2 ',PPP2
1146 C WRITE(6,2346)' PPP4 ',PPP4
1147 C Invariant Masses of chains
1148 C ORIGINAL CHAIN
1149  amchor=sqrt((ppr(4)+pta(4))**2-(ppr(3)+pta(3))**2
1150  * -(ppr(2)+pta(2))**2-(ppr(1)+pta(1))**2)
1151  amchn1=sqrt((ppp1(4)+ppp2(4))**2-(ppp1(3)+ppp2(3))**2
1152  * -(ppp1(2)+ppp2(2))**2-(ppp1(1)+ppp2(1))**2)
1153 C Chain 1 is q-aq restrict mass >1.5 GeV
1154  chamal=1.d0
1155  IF(ifb1.GE.3.OR.isaq.GE.9)chamal=1.5d0
1156  IF(amchn1.LE.chamal)THEN
1157 C IREJ=1
1158 C RETURN
1159  go to 3465
1160  ENDIF
1161  amchn2=sqrt((ppp3(4)+ppp4(4))**2-(ppp3(3)+ppp4(3))**2
1162  * -(ppp3(2)+ppp4(2))**2-(ppp3(1)+ppp4(1))**2)
1163 C Chain 2 is qq-q restrict mass >2.5 GeV
1164  chamal=1.8d0
1165  IF(ifb2.GE.3.OR.ifb3.GE.3.OR.isq.GE.3)chamal=2.5d0
1166  IF(amchn2.LE.chamal)THEN
1167 C IREJ=1
1168 C RETURN
1169  go to 3465
1170  ENDIF
1171  pxchk=ppr(1)+pta(1)-ppp1(1)-ppp2(1)-ppp3(1)-ppp4(1)
1172  pychk=ppr(2)+pta(2)-ppp1(2)-ppp2(2)-ppp3(2)-ppp4(2)
1173  pzchk=ppr(3)+pta(3)-ppp1(3)-ppp2(3)-ppp3(3)-ppp4(3)
1174  pechk=ppr(4)+pta(4)-ppp1(4)-ppp2(4)-ppp3(4)-ppp4(4)
1175  IF(ipco.GE.1)THEN
1176  WRITE(6,'(A/8E12.4,I5)')
1177  * ' Chain masses AMCH,AMCHOR,AMCHN1,AMCHN2,PZCHK,PECHK,
1178  * PXCHK,PYCHK,NOBAM ',
1179  * amch,amchor,amchn1,amchn2,pzchk,pechk,pxchk,pychk,nobam
1180  ENDIF
1181 C Lorentz parameters of chains
1182  gamor=(ppr(4)+pta(4))/amch
1183  bgxor=(ppr(1)+pta(1))/amch
1184  bgyor=(ppr(2)+pta(2))/amch
1185  bgzor=(ppr(3)+pta(3))/amch
1186  gamch1=(ppp1(4)+ppp2(4))/amchn1
1187  bgxch1=(ppp1(1)+ppp2(1))/amchn1
1188  bgych1=(ppp1(2)+ppp2(2))/amchn1
1189  bgzch1=(ppp1(3)+ppp2(3))/amchn1
1190  gamch2=(ppp3(4)+ppp4(4))/amchn2
1191  bgxch2=(ppp3(1)+ppp4(1))/amchn2
1192  bgych2=(ppp3(2)+ppp4(2))/amchn2
1193  bgzch2=(ppp3(3)+ppp4(3))/amchn2
1194  IF(ipco.GE.1)THEN
1195  WRITE(6,2346)' L.Parm in ',gam,bgx,bgy,bgz
1196  WRITE(6,2346)' L.Parm OR ',gamor,bgxor,bgyor,bgzor
1197  WRITE(6,2346)' L.Parm C1 ',gamch1,bgxch1,bgych1,bgzch1
1198  WRITE(6,2346)' L.Parm C2 ',gamch2,bgxch2,bgych2,bgzch2
1199  ENDIF
1200  peor=amch*gamor
1201  pzor=amch*bgzor
1202  pech1=amchn1*gamch1
1203  pzch1=amchn1*bgzch1
1204  pech2=amchn2*gamch2
1205  pzch2=amchn2*bgzch2
1206  IF(ipco.GE.1)THEN
1207  WRITE(6,'(A,6E12.4)')' PEOR,PECH1,PECH2,PZOR,PZCH1,PZCH2',
1208  * peor,pech1,pech2,pzor,pzch1,pzch2
1209  ENDIF
1210  noba1=3
1211  noba2=6
1212 C WRITE(6,'(A,2I5)')' Jet1 IFB1,ISAQ ',IFB1,ISAQ
1213  CALL hadjet(nhad1,amchn1,ppp1,ppp2,gamch1,bgxch1,bgych1,
1214  * bgzch1, ifb1,isaq,ifb3,ifb4,i1,i2,noba1,nnch,norig)
1215 C DO 42 I=1,NHAD1
1216 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
1217 C + IBARF(I),NREF(I),ANF(I)
1218 C 42 CONTINUE
1219 C DIMENSION ANFF(NFIMAX),PXFF(NFIMAX),PYFF(NFIMAX),PZFF(NFIMAX),
1220 C +HEFF(NFIMAX),AMFF(NFIMAX),
1221 C *ICHFF(NFIMAX),IBARFF(NFIMAX),NREFF(NFIMAX)
1222 C Intermediate store of hadrons from jet 1
1223  DO 2348 i=1,nhad1
1224  anff(i)=anf(i)
1225  pxff(i)=pxf(i)
1226  pyff(i)=pyf(i)
1227  pzff(i)=pzf(i)
1228  heff(i)=hef(i)
1229  amff(i)=amf(i)
1230  ichff(i)=ichf(i)
1231  ibarff(i)=ibarf(i)
1232  nreff(i)=nref(i)
1233  2348 CONTINUE
1234 C WRITE(6,'(A,3I5)')' Jet2 IFB2,ISQ,IFB3 ',IFB2,ISQ,IFB3
1235  CALL hadjet(nhad2,amchn2,ppp3,ppp4,gamch2,bgxch2,bgych2,
1236  * bgzch2, ifb2,isq,ifb3,ifb4,i1,i2,noba2,nnch,norig)
1237 C DO 41 I=1,NHAD2
1238 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
1239 C + IBARF(I),NREF(I),ANF(I)
1240 C 41 CONTINUE
1241  nhad=nhad1+nhad2
1242  DO 2349 i=nhad2+1,nhad
1243  ii=i-nhad2
1244  anf(i)=anff(ii)
1245  pxf(i)=pxff(ii)
1246  pyf(i)=pyff(ii)
1247  pzf(i)=pzff(ii)
1248  hef(i)=heff(ii)
1249  amf(i)=amff(ii)
1250  ichf(i)=ichff(ii)
1251  ibarf(i)=ibarff(ii)
1252  nref(i)=nreff(ii)
1253  2349 CONTINUE
1254  ENDIF
1255 C==================================================================
1256 C==================================================================
1257 C Different options:
1258  IF(nobam.EQ.4)THEN
1259 C NOBAM=4: Quark-DIQUARK Jet
1260 C
1261 C Work out approximate x-Values of quark and diquark
1262  xdiq=pta(4)*2.d0/cmener
1263  xqua=ppr(4)*2.d0/cmener
1264 C WRITE(6,'(A,2E12.4)')' HADJCK:XDIQ,XQUA ',XDIQ,XQUA
1265 C Select valence quark x for one of the diquark quarks
1266 C But only if XDIQ > 1.5*XQUA
1267 C IF(XDIQ.LE.1.5D0*XQUA)THEN
1268 C IREJ=1
1269 C RETURN
1270 C ENDIF
1271 C Select x values of sea quark pair
1272  icou=0
1273  2234 CONTINUE
1274  icou=icou+1
1275  IF(icou.GE.100)THEN
1276  irej=1
1277  RETURN
1278  ENDIF
1279  CALL xseapa(cmener,xqua/2.d0,isq,isaq,xsq,xsaq,irej)
1280  IF(xsaq.GE.xqua/2.d0)go to 2234
1281  IF(ipco.GE.1)THEN
1282  WRITE(6,'(A,2E12.4)')' HADJCK:XSQ,XSAQ ',xsq,xsaq
1283  ENDIF
1284 * target valence quarks xtvqi
1285  xvthro=cvq/cmener
1286  ivthr=0
1287  3466 CONTINUE
1288  IF(ivthr.EQ.50)THEN
1289  irej=1
1290  WRITE(6,*)' HADJSE 3466 reject IVTHR 50'
1291  RETURN
1292  ENDIF
1293  ivthr=ivthr+1
1294  xvthr=xvthro/(51-ivthr)
1295  unoprv=unon
1296  380 CONTINUE
1297  IF(xvthr.GT.0.05)THEN
1298  IF(xvthr.GT.0.66d0*xdiq)THEN
1299  irej=1
1300  RETURN
1301  ENDIF
1302  xtvqi=betrej(0.5,unoprv,xvthr,0.66d0*xdiq)
1303  ELSE
1304  n90=0
1305  390 CONTINUE
1306  n90=n90+1
1307  IF(n90.GE.20)THEN
1308  irej=1
1309  RETURN
1310  ENDIF
1311  xtvqi=dbetar(0.5,unoprv)
1312  IF ((xtvqi.LT.xvthr).OR.(xtvqi.GT.0.66d0*xdiq))
1313  * goto 390
1314  ENDIF
1315  xdiqq=xdiq-xtvqi
1316  IF(ipco.GE.1)THEN
1317  WRITE(6,'(A,2E12.4)')' HADJCK:XDIQQ,XTVQI ',xdiqq,xtvqi
1318  ENDIF
1319 C We form two strings(2-1) aq-q and (4-3) q-qq
1320 C 1 q with XDIQ-XTVQI-XSQ
1321 C 1 q with XDIQ-XTVQI (j.r.5.5.96)
1322 C 2 aq with XSAQ
1323 C 3 qq with XTVQI+XSQ
1324 C 3 qq with XTVQI (j.r.5.5.96)
1325 C 4 q with XQUA-XSAQ
1326 C with 4-momenta PPP1,PPP2,PPP3,PPP4
1327  DO 3345 i=1,4
1328 C PPP1(I)=PTA(I)*(XDIQ-XTVQI-XSQ)/XDIQ
1329  ppp1(i)=pta(i)*(xdiq-xtvqi)/xdiq
1330 C PPP3(I)=PTA(I)*(XTVQI+XSQ)/XDIQ
1331  ppp3(i)=pta(i)*(xtvqi)/xdiq
1332  ppp2(i)=ppr(i)*xsaq/xqua
1333  ppp4(i)=ppr(i)*(xqua-xsaq)/xqua
1334  3345 CONTINUE
1335  3346 FORMAT(a,4e12.4)
1336  IF(ipco.GE.1)THEN
1337  WRITE(6,3346)' PPR ',ppr
1338  WRITE(6,3346)' PPP1 ',ppp1
1339  WRITE(6,3346)' PPP3 ',ppp3
1340  WRITE(6,3346)' PTA ',pta
1341  WRITE(6,3346)' PPP2 ',ppp2
1342  WRITE(6,3346)' PPP4 ',ppp4
1343  ENDIF
1344 C Invariant Masses of chains
1345 C ORIGINAL CHAIN
1346  amchor=sqrt((ppr(4)+pta(4))**2-(ppr(3)+pta(3))**2
1347  * -(ppr(2)+pta(2))**2-(ppr(1)+pta(1))**2)
1348  amchn1=sqrt((ppp1(4)+ppp2(4))**2-(ppp1(3)+ppp2(3))**2
1349  * -(ppp1(2)+ppp2(2))**2-(ppp1(1)+ppp2(1))**2)
1350 C Chain 1 is q-aq restrict mass >1.5 GeV
1351  chamal=1.d0
1352  IF(ifb3.GE.3.OR.isaq.GE.9)chamal=1.5d0
1353  IF(amchn1.LE.chamal)THEN
1354 C IREJ=1
1355 C RETURN
1356  go to 3466
1357  ENDIF
1358  amchn2=sqrt((ppp3(4)+ppp4(4))**2-(ppp3(3)+ppp4(3))**2
1359  * -(ppp3(2)+ppp4(2))**2-(ppp3(1)+ppp4(1))**2)
1360 C Chain 2 is qq-q restrict mass >2.5 GeV
1361  chamal=1.8d0
1362  IF(ifb2.GE.3.OR.ifb1.GE.3.OR.isq.GE.3)chamal=2.5d0
1363  IF(amchn2.LE.chamal)THEN
1364 C IREJ=1
1365 C RETURN
1366  go to 3466
1367  ENDIF
1368  pxchk=ppr(1)+pta(1)-ppp1(1)-ppp2(1)-ppp3(1)-ppp4(1)
1369  pychk=ppr(2)+pta(2)-ppp1(2)-ppp2(2)-ppp3(2)-ppp4(2)
1370  pzchk=ppr(3)+pta(3)-ppp1(3)-ppp2(3)-ppp3(3)-ppp4(3)
1371  pechk=ppr(4)+pta(4)-ppp1(4)-ppp2(4)-ppp3(4)-ppp4(4)
1372  IF(ipco.GE.1)THEN
1373  WRITE(6,'(A/8E12.4,I5)')
1374  * ' Chain masses AMCH,AMCHOR,AMCHN1,AMCHN2,PZCHK,PECHK,
1375  * PXCHK,PYCHK,NOBAM ',
1376  * amch,amchor,amchn1,amchn2,pzchk,pechk,pxchk,pychk,nobam
1377  ENDIF
1378 C Lorentz parameters of chains
1379  gamor=(ppr(4)+pta(4))/amch
1380  bgxor=(ppr(1)+pta(1))/amch
1381  bgyor=(ppr(2)+pta(2))/amch
1382  bgzor=(ppr(3)+pta(3))/amch
1383  gamch1=(ppp1(4)+ppp2(4))/amchn1
1384  bgxch1=(ppp1(1)+ppp2(1))/amchn1
1385  bgych1=(ppp1(2)+ppp2(2))/amchn1
1386  bgzch1=(ppp1(3)+ppp2(3))/amchn1
1387  gamch2=(ppp3(4)+ppp4(4))/amchn2
1388  bgxch2=(ppp3(1)+ppp4(1))/amchn2
1389  bgych2=(ppp3(2)+ppp4(2))/amchn2
1390  bgzch2=(ppp3(3)+ppp4(3))/amchn2
1391  IF(ipco.GE.1)THEN
1392  WRITE(6,3346)' L.Parm in ',gam,bgx,bgy,bgz
1393  WRITE(6,3346)' L.Parm OR ',gamor,bgxor,bgyor,bgzor
1394  WRITE(6,3346)' L.Parm C1 ',gamch1,bgxch1,bgych1,bgzch1
1395  WRITE(6,3346)' L.Parm C2 ',gamch2,bgxch2,bgych2,bgzch2
1396  ENDIF
1397  peor=amch*gamor
1398  pzor=amch*bgzor
1399  pech1=amchn1*gamch1
1400  pzch1=amchn1*bgzch1
1401  pech2=amchn2*gamch2
1402  pzch2=amchn2*bgzch2
1403  IF(ipco.GE.1)THEN
1404  WRITE(6,'(A,6E12.4)')' PEOR,PECH1,PECH2,PZOR,PZCH1,PZCH2',
1405  * peor,pech1,pech2,pzor,pzch1,pzch2
1406  ENDIF
1407  noba1=3
1408  noba2=4
1409 C WRITE(6,'(A,2I5)')' Jet1 ISAQ,IFB3 ',ISAQ,IFB3
1410  CALL hadjet(nhad1,amchn1,ppp2,ppp1,gamch1,bgxch1,bgych1,
1411  * bgzch1, isaq,ifb3,ifb1,ifb4,i1,i2,noba1,nnch,norig)
1412 C DO 342 I=1,NHAD1
1413 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
1414 C + IBARF(I),NREF(I),ANF(I)
1415 C 342 CONTINUE
1416 C DIMENSION ANFF(NFIMAX),PXFF(NFIMAX),PYFF(NFIMAX),PZFF(NFIMAX),
1417 C +HEFF(NFIMAX),AMFF(NFIMAX),
1418 C *ICHFF(NFIMAX),IBARFF(NFIMAX),NREFF(NFIMAX)
1419 C Intermediate store of hadrons from jet 1
1420  DO 3348 i=1,nhad1
1421  anff(i)=anf(i)
1422  pxff(i)=pxf(i)
1423  pyff(i)=pyf(i)
1424  pzff(i)=pzf(i)
1425  heff(i)=hef(i)
1426  amff(i)=amf(i)
1427  ichff(i)=ichf(i)
1428  ibarff(i)=ibarf(i)
1429  nreff(i)=nref(i)
1430  3348 CONTINUE
1431 C WRITE(6,'(A,3I5)')' Jet2 IFB1,ISQ,IFB2 ',IFB1,ISQ,IFB2
1432  CALL hadjet(nhad2,amchn2,ppp4,ppp3,gamch2,bgxch2,bgych2,
1433  * bgzch2, ifb1,isq,ifb2,ifb4,i1,i2,noba2,nnch,norig)
1434 C DO 341 I=1,NHAD2
1435 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
1436 C + IBARF(I),NREF(I),ANF(I)
1437 C 341 CONTINUE
1438  nhad=nhad1+nhad2
1439  DO 3349 i=nhad2+1,nhad
1440  ii=i-nhad2
1441  anf(i)=anff(ii)
1442  pxf(i)=pxff(ii)
1443  pyf(i)=pyff(ii)
1444  pzf(i)=pzff(ii)
1445  hef(i)=heff(ii)
1446  amf(i)=amff(ii)
1447  ichf(i)=ichff(ii)
1448  ibarf(i)=ibarff(ii)
1449  nref(i)=nreff(ii)
1450  3349 CONTINUE
1451  ENDIF
1452 C==================================================================
1453  heft=0.d0
1454  IF(ipco.GE.1)THEN
1455  DO 40 i=1,nhad
1456  IF(ibarf(i).EQ.500)go to 4040
1457  heft=heft+hef(i)
1458  WRITE(6,4050)i,pxf(i),pyf(i),pzf(i),hef(i),amf(i), ichf(i),
1459  + ibarf(i),nref(i),anf(i)
1460  4050 FORMAT(' JET ',i5,5f12.4,3i5,a10)
1461  4040 CONTINUE
1462  40 CONTINUE
1463  ENDIF
1464  hefff=amch*gam
1465  IF(ipco.GE.1)THEN
1466  WRITE(6,'(A,I5,2E12.4)')'1IREJ,HEFT,HEFFF ',
1467  * irej,heft,hefff
1468  ENDIF
1469 C IREJ=1
1470  IF(irej.EQ.1)RETURN
1471 C
1472  RETURN
1473  END
1474 
1475 ************************************************************************
1476 ************************************************************************
1477 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
1478 C
1479  SUBROUTINE hadjse(NHAD,AMCH,PPR,PTA,GAM,BGX,BGY,BGZ, IFB1,IFB2,
1480  +ifb3,ifb4,i1,i2,nobam,nnch,norig,irej,iissqq)
1481  IMPLICIT DOUBLE PRECISION (a-h,o-z)
1482  SAVE
1483 C
1484 C HADJSE Fragmentation of vv,vs and sv chains CK method
1485 C using Glauber sea quarks
1486 C j.r.7/96
1487 C
1488 C HADJET DOES ALL THE NECESSARY ROTATIONS AND LORENTZ TRANSFORMS AND
1489 C CALL CALBAM (BAMJET)
1490 C
1491 C NHAD = NUMBER OF HADRONS CREATED (OUTPUT) = IHAD (CALBAM)
1492 C******** PRODUCED PARTICLES IN COMMON /CMSRES/
1493 C NOTE: NOW IN /FINPAR/ HJM 21/06/90
1494 C AMCH = INVARIANT MASS OF JET (INPUT)
1495 C PPR = 4-MOMENTUM OF FORWARD GOING PARTON (PROJECTILE)(INPUT)
1496 C PTA = 4-MOMENTUM OF BACKWARD GOING PARTON (TARGET)(INPUT)
1497 C GAM,BGX,BGY,BGZ = LORENTZ PARAMETERS OF JET CMS RELATIVE TO
1498 C COLLISION CMS (INPUT)
1499 C
1500 C--------------------------------------------------------------------
1501 C CALBAM(NNCH,I1,I2,IFB1,IFB2,IFB3,IFB4,AMCH,NOBAM,IHAD)
1502 C SAMPLING OF Q-AQ, Q-QQ, QQ-AQQ CHAINS
1503 C USING BAMJET(IHAD,IFB1,IFB2,IFB3,IFB4,AMCH,NOBAM)-----FOR NNCH=0
1504 C OR PARJET(IHAD,ICH1=I1 OR I2)------FOR NNCH=-1 OR +1
1505 C-------------------------------------------------------------------
1506 C IHAD : NUMBER OF PRODUCED PARTICLES
1507 C NNCH : CALL BAMJET FOR NNCH=0
1508 C CALL PARJET FOR NNCH=+1 ICH1=I1
1509 C FOR NNCH=-1 ICH1=I2
1510 C jet not existing for NNCH=+/-99, i.e. IHAD=0
1511 C PRODUCED PARTICLES IN CHAIN REST FRAME ARE IN COMMON /FINPAR/
1512 C AMCH : INVARIANT MASS OF CHAIN (BAMJET)
1513 C
1514 C NOBAM : = 3 QUARK-ANTIQUARK JET QUARK FLAVORS : IFB1,IFB2
1515 C OR ANTIQUARK-QUARK JET IN ANY ORDER
1516 C
1517 C = 4 QUARK-DIQUARK JET, FLAVORS : QU : IFB1, DIQU :IFB2,IFB
1518 C OR ANTIQUARK-ANTIDIQUARK JET
1519 C
1520 C
1521 C = 5 DIQUARK-ANTIDIQUARK JET
1522 C OR ANTIDIQUARK-DIQUARK JET
1523 C FLAVORS : DIQU : IFB1,IFB2, ANTIDIQU : IFB3,IFB4
1524 C IN ANY ORDER
1525 C
1526 C = 6 DIQUARK-QUARK JET, FLAVORS : DIQU : IFB1,IFB2 QU: IFB
1527 C OR ANTIDIQUARK-ANTIQUARK JET
1528 C
1529 C IFBI : FLAVORS : 1,2,3,4 = U,D,S,C 7,8,9,10 = AU,AD,AS,AC
1530 C
1531 C I1,I2 : NUMBER LABEL OF A HADRON CREATED BY PARJET
1532 C
1533 C NORMALLY IN BAMJET THE QUARKS MOVE FORWARD (POSITIVE Z-DIRECTION)
1534 C THE QUARK FLAVORS ARE FIRST GIVEN
1535 C CALBAM ALLOWS EITHER THE QUARK OR ANTIQUARK (DIQU) TO MOVE FORWARD
1536 C THE FORWARD GOING FLAVORS ARE GIVEN FIRST
1537 C
1538 C =====================================================================
1539 *KEEP,DFINPA.
1540  CHARACTER*8 anf,anff
1541  parameter(nfimax=249)
1542  COMMON /dfinpa/ anf(nfimax),pxf(nfimax),pyf(nfimax),pzf(nfimax),
1543  +hef(nfimax),amf(nfimax), ichf(nfimax),ibarf(nfimax),nref(nfimax)
1544  COMMON /dfinpz/iormo(nfimax),idaug1(nfimax),idaug2(nfimax),
1545  * istath(nfimax)
1546  dimension anff(nfimax),pxff(nfimax),pyff(nfimax),pzff(nfimax),
1547  +heff(nfimax),amff(nfimax),
1548  *ichff(nfimax),ibarff(nfimax),nreff(nfimax),iormoo(nfimax)
1549  CHARACTER*80 title
1550  CHARACTER*8 projty,targty
1551 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
1552 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
1553  COMMON /user1/title,projty,targty
1554  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
1555 *KEEP,DPRIN.
1556  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
1557 *KEND.
1558 *--------------------------------------------------------------------
1559 *-------------------------------------------------------------------
1560  COMMON /jspart/pxp(1000),pyp(1000),pzp(1000),hepp(1000),nnnp
1561  COMMON /jspar/pxj(1000),pyj(1000),pzj(1000),hej(1000),nnnpj
1562 ************************************************************************
1563 ************************************************************************
1564  common/ifragm/ifrag
1565  dimension ppr(4),pta(4)
1566 *KEEP,XSEADI.
1567  COMMON /xseadi/ xseacu,unon,unom,unosea,cvq,cdq,csea,ssmima,
1568  +ssmimq,vvmthr
1569  COMMON /hdjase/nhse1,nhse2,nhse3,nhase1,nhase2,nhase3
1570  dimension ppp1(4),ppp2(4),ppp3(4),ppp4(4)
1571  DATA tiny/1.0d-3/
1572 C----------------------------------------------------------------------
1573 C
1574 C WRITE(6,'(A,3I5)')' Jet0 IFB1,IFB2,IFB3 ',IFB1,IFB2,IFB3
1575  IF(ipco.GE.0)WRITE(6,*)
1576  *' HADJSE(NHAD,AMCH,PPR,PTA,GAM,BGX,BGY,BGZ, IFB1,IFB2,',
1577  +'IFB3,IFB4,I1,I2,NOBAM,NNCH,NORIG,IREJ),IPCO',
1578  *nhad,amch,ppr,pta,gam,bgx,bgy,bgz, ifb1,ifb2,
1579  +ifb3,ifb4,i1,i2,nobam,nnch,norig,irej,ipco
1580  irej=0
1581  IF(ipco.GE.0) THEN
1582  WRITE(6,'(A,3I5)')' HADJSE Jet0 IFB1,IFB2,IFB3',ifb1,ifb2,ifb3
1583  WRITE(6,'(4E12.4)') gam,bgx,bgy,bgz,ppr
1584  WRITE(6,1000) nhad,amch,pta, ifb1,ifb2,ifb3,ifb4,i1,i2,nobam,
1585  + nnch,norig
1586  1000 FORMAT(10x,i10,5f10.3/10x,9i10)
1587  ENDIF
1588  IF(abs(nnch).EQ.99) THEN
1589  nhad=0
1590 C IPCO=0
1591  RETURN
1592  ENDIF
1593 C=================================================================
1594 C=================================================================
1595 C DIQUARK-QUARK NOBAM=6
1596 C=================================================================
1597 C=================================================================
1598 C Different options:
1599  IF(nobam.EQ.6)THEN
1600  IF(ipco.GE.0)WRITE(6,*)' DIQUARK-QUARK NOBAM=6'
1601 C NOBAM=6: Diquark--Quark Jet
1602 C
1603 C Work out approximate x-Values of diquark and quark
1604  IF(cmener.LT.1.d-3)THEN
1605  WRITE(6,*)' CMENER=0. HADJSE',cmener
1606  stop
1607  ENDIF
1608  xdiq=ppr(4)*2.d0/cmener
1609  xqua=pta(4)*2.d0/cmener
1610  IF(ipco.GE.0) THEN
1611  WRITE(6,'(A,2E12.4)')' HADJSE:XDIQ,XQUA ',xdiq,xqua
1612  ENDIF
1613 C Select sea quark x for one of the diquark quarks
1614 C Select x values of sea quark pair
1615  icou=0
1616  1234 CONTINUE
1617  icou=icou+1
1618  IF(icou.GE.200)THEN
1619  irej=1
1620  IF(ipco.GE.0)WRITE(6,*)' HADJSE reject icou 100'
1621  RETURN
1622  ENDIF
1623  CALL xseapa(cmener,xqua/2.d0,isq,isaq,xsq,xsaq,irej)
1624  iissqq=isq
1625  IF(irej.GE.1)THEN
1626  IF(ipco.GE.0)WRITE(6,*)' HADJSE reject XSEAPA'
1627  RETURN
1628  ENDIF
1629  IF(xsaq.GE.2.d0*xqua/3.d0)go to 1234
1630  IF(ipco.GE.0) THEN
1631  WRITE(6,*)' HADJSE,XSEAPA:XSQ,XSAQ,ISQ,ISAQ ',
1632  * xsq,xsaq,isq,isaq
1633  ENDIF
1634 * projectile valence quarks xpvqi
1635 C here selected as sea quark! from 1/x-sea
1636  xvthro=cvq/cmener
1637  ivthr=0
1638  3465 CONTINUE
1639  IF(ivthr.EQ.100)THEN
1640  irej=1
1641  IF(isq.EQ.3)irej=3
1642  IF(ipco.GE.0)WRITE(6,*)' HADJSE 3465 reject IVTHR 50'
1643  RETURN
1644  ENDIF
1645  ivthr=ivthr+1
1646  xvthr=xvthro/(101-ivthr)
1647  unoprv=unon
1648  80 CONTINUE
1649  IF(xvthr.GT.0.66d0*xdiq)THEN
1650  IF(ipco.GE.0)WRITE(6,*)' HADJSE 80 reject XVTHR too great',
1651  * xvthr
1652  irej=1
1653  IF(isq.EQ.3)irej=3
1654  RETURN
1655  ENDIF
1656  xpvqi=sampey(xvthr,0.66d0*xdiq)
1657  xdiqq=xdiq-xpvqi
1658  IF(ipco.GE.0) THEN
1659  WRITE(6,'(A,2E12.4)')' HADJSE:XDIQQ,XPVQI ',xdiqq,xpvqi
1660  ENDIF
1661 C
1662 C We form two strings(1-2) q-aq and (3-4) qq-q
1663 C 1 q with XDIQ-XPVQI-XSQ
1664 C 1 q with XDIQ-XPVQI (j.r.5.5.96)
1665 C 2 aq with XSAQ
1666 C 3 qq with XPVQI+XSQ
1667 C 3 qq with XPVQI (j.r.5.5.96)
1668 C 4 q with XQUA-XSAQ
1669 C with 4-momenta PPP1,PPP2,PPP3,PPP4
1670  DO 2345 i=1,4
1671  IF(xdiq.LE.1.d-15.OR.xqua.LT.1.d-15)THEN
1672  irej=1
1673  IF(isq.EQ.3)irej=3
1674  IF(ipco.GE.0)WRITE(6,*)' HADJSE reject 2345 XDIQ,',
1675  * 'XQUA too small ',
1676  * xdiq,xqua
1677  RETURN
1678  ENDIF
1679 C PPP1(I)=PPR(I)*(XDIQ-XPVQI-XSQ)/XDIQ
1680  ppp1(i)=ppr(i)*(xdiq-xpvqi)/xdiq
1681 C PPP3(I)=PPR(I)*(XPVQI+XSQ)/XDIQ
1682  ppp3(i)=ppr(i)*(xpvqi)/xdiq
1683  ppp2(i)=pta(i)*xsaq/xqua
1684  ppp4(i)=pta(i)*(xqua-xsaq)/xqua
1685  2345 CONTINUE
1686  2346 FORMAT(a,5e12.4)
1687  IF(ipco.GE.0) THEN
1688  WRITE(6,2346)' PPR ',ppr
1689  WRITE(6,2346)' PPP1 ',ppp1
1690  WRITE(6,2346)' PPP3 ',ppp3
1691  WRITE(6,2346)' PTA ',pta
1692  WRITE(6,2346)' PPP2 XSAQ',ppp2,xsaq
1693  WRITE(6,2346)' PPP4 ',ppp4
1694  ENDIF
1695 C Invariant Masses of chains
1696 C ORIGINAL CHAIN
1697  amchor=sqrt((ppr(4)+pta(4))**2-(ppr(3)+pta(3))**2
1698  * -(ppr(2)+pta(2))**2-(ppr(1)+pta(1))**2)
1699  amchn1=sqrt((ppp1(4)+ppp2(4))**2-(ppp1(3)+ppp2(3))**2
1700  * -(ppp1(2)+ppp2(2))**2-(ppp1(1)+ppp2(1))**2)
1701  IF(amchor.LE.tiny)THEN
1702  WRITE(6,2346)' PPR ',ppr
1703  WRITE(6,2346)' PPP1 ',ppp1
1704  WRITE(6,2346)' PPP3 ',ppp3
1705  WRITE(6,2346)' PTA ',pta
1706  WRITE(6,2346)' PPP2 ',ppp2
1707  WRITE(6,2346)' PPP4 ',ppp4
1708  ENDIF
1709  IF(amchn1.LE.tiny)THEN
1710  WRITE(6,2346)' PPR ',ppr
1711  WRITE(6,2346)' PPP1 ',ppp1
1712  WRITE(6,2346)' PPP3 ',ppp3
1713  WRITE(6,2346)' PTA ',pta
1714  WRITE(6,2346)' PPP2 ',ppp2
1715  WRITE(6,2346)' PPP4 ',ppp4
1716  ENDIF
1717 C Chain 1 is q-aq restrict mass >1.5 GeV
1718  chamal=0.8d0
1719  IF(ifb1.GE.3.OR.isaq.GE.9)chamal=1.2d0
1720  IF(amchn1.LE.chamal)THEN
1721  IF(ipco.GE.0)WRITE (6,*).LE.'HADJSE jump1AMCHN1CHAMAL AMCHOR',
1722  * amchn1,chamal,amchor
1723  go to 3465
1724  ENDIF
1725  amchn2=sqrt((ppp3(4)+ppp4(4))**2-(ppp3(3)+ppp4(3))**2
1726  * -(ppp3(2)+ppp4(2))**2-(ppp3(1)+ppp4(1))**2)
1727 C Chain 2 is qq-q restrict mass >2.5 GeV
1728  chamal=1.5d0
1729  IF(ifb2.GE.3.OR.ifb3.GE.3.OR.isq.GE.3)chamal=2.0d0
1730  IF(amchn2.LE.chamal)THEN
1731 C IREJ=1
1732 C RETURN
1733  IF(ipco.GE.0)WRITE (6,*).LE.'HADJSE jump AMCHN2CHAMAL ',
1734  * amchn2,chamal
1735  go to 3465
1736  ENDIF
1737  pxchk=ppr(1)+pta(1)-ppp1(1)-ppp2(1)-ppp3(1)-ppp4(1)
1738  pychk=ppr(2)+pta(2)-ppp1(2)-ppp2(2)-ppp3(2)-ppp4(2)
1739  pzchk=ppr(3)+pta(3)-ppp1(3)-ppp2(3)-ppp3(3)-ppp4(3)
1740  pechk=ppr(4)+pta(4)-ppp1(4)-ppp2(4)-ppp3(4)-ppp4(4)
1741  IF(ipco.GE.0)THEN
1742  WRITE(6,'(A/8E12.4,I5)')
1743  * ' Chain masses AMCH,AMCHOR,AMCHN1,AMCHN2,PZCHK,PECHK,
1744  * PXCHK,PYCHK,NOBAM ',
1745  * amch,amchor,amchn1,amchn2,pzchk,pechk,pxchk,pychk,nobam
1746  ENDIF
1747 C Lorentz parameters of chains
1748  IF(amch.LE.1.d-15)THEN
1749  irej=1
1750  IF(isq.EQ.3)irej=3
1751  IF(ipco.GE.0)WRITE(6,*)' HADJSE rejection AMCH too small ',
1752  * amch
1753  RETURN
1754  ENDIF
1755  gamor=(ppr(4)+pta(4))/amch
1756  bgxor=(ppr(1)+pta(1))/amch
1757  bgyor=(ppr(2)+pta(2))/amch
1758  bgzor=(ppr(3)+pta(3))/amch
1759  IF(amchn1.LE.1.d-15)THEN
1760  IF(ipco.GE.0)WRITE(6,*)' HADJSE rejection AMCHN1 too small ',
1761  * amchn1
1762  irej=1
1763  IF(isq.EQ.3)irej=3
1764  RETURN
1765  ENDIF
1766  gamch1=(ppp1(4)+ppp2(4))/amchn1
1767  bgxch1=(ppp1(1)+ppp2(1))/amchn1
1768  bgych1=(ppp1(2)+ppp2(2))/amchn1
1769  bgzch1=(ppp1(3)+ppp2(3))/amchn1
1770  IF(amchn2.LE.1.d-15)THEN
1771  IF(ipco.GE.0)WRITE(6,*)' HADJSE rejection AMCHN2 too small ',
1772  * amchn2
1773  irej=1
1774  IF(isq.EQ.3)irej=3
1775  RETURN
1776  ENDIF
1777  gamch2=(ppp3(4)+ppp4(4))/amchn2
1778  bgxch2=(ppp3(1)+ppp4(1))/amchn2
1779  bgych2=(ppp3(2)+ppp4(2))/amchn2
1780  bgzch2=(ppp3(3)+ppp4(3))/amchn2
1781  IF(ipco.GE.0)THEN
1782  WRITE(6,2346)' L.Parm in ',gam,bgx,bgy,bgz
1783  WRITE(6,2346)' L.Parm OR ',gamor,bgxor,bgyor,bgzor
1784  WRITE(6,2346)' L.Parm C1 ',gamch1,bgxch1,bgych1,bgzch1
1785  WRITE(6,2346)' L.Parm C2 ',gamch2,bgxch2,bgych2,bgzch2
1786  ENDIF
1787  peor=amch*gamor
1788  pzor=amch*bgzor
1789  pech1=amchn1*gamch1
1790  pzch1=amchn1*bgzch1
1791  pech2=amchn2*gamch2
1792  pzch2=amchn2*bgzch2
1793  IF(ipco.GE.0)THEN
1794  WRITE(6,'(A,6E12.4)')' PEOR,PECH1,PECH2,PZOR,PZCH1,PZCH2',
1795  * peor,pech1,pech2,pzor,pzch1,pzch2
1796  ENDIF
1797  noba1=3
1798  noba2=6
1799  IF(ipco.GE.-1)THEN
1800 C WRITE(6,'(A,2I5)')' Jet1 IFB1,ISAQ ',IFB1,ISAQ
1801  ENDIF
1802  CALL hadjet(nhad1,amchn1,ppp1,ppp2,gamch1,bgxch1,bgych1,
1803  * bgzch1, ifb1,isaq,ifb3,ifb4,i1,i2,noba1,nnch,norig)
1804 C DO 42 I=1,NHAD1
1805 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
1806 C + IBARF(I),NREF(I),ANF(I)
1807 C 42 CONTINUE
1808  IF(ipco.GE.-1)THEN
1809 C WRITE(6,'(A,2I5)')' Jet1 NHAD1 ',NHAD1
1810  ENDIF
1811 C DIMENSION ANFF(NFIMAX),PXFF(NFIMAX),PYFF(NFIMAX),PZFF(NFIMAX),
1812 C + HEFF(NFIMAX),AMFF(NFIMAX),
1813 C * ICHFF(NFIMAX),IBARFF(NFIMAX),NREFF(NFIMAX)
1814 C Intermediate store of hadrons from jet 1
1815  DO 2348 i=1,nhad1
1816  anff(i)=anf(i)
1817  pxff(i)=pxf(i)
1818  pyff(i)=pyf(i)
1819  pzff(i)=pzf(i)
1820  heff(i)=hef(i)
1821  amff(i)=amf(i)
1822  ichff(i)=ichf(i)
1823  ibarff(i)=ibarf(i)
1824  nreff(i)=nref(i)
1825  iormoo(i)=iormo(i)
1826  2348 CONTINUE
1827  IF(ipco.GE.0)THEN
1828  WRITE(6,'(A,3I5)')' Jet2 IFB2,ISQ,IFB3 ',ifb2,isq,ifb3
1829  ENDIF
1830  CALL hadjet(nhad2,amchn2,ppp3,ppp4,gamch2,bgxch2,bgych2,
1831  * bgzch2, ifb2,isq,ifb3,ifb4,i1,i2,noba2,nnch,norig)
1832  IF(ipco.GE.-1)THEN
1833 C WRITE(6,'(A,2I5)')' Jet2 NHAD2 ',NHAD2
1834  ENDIF
1835 C DO 41 I=1,NHAD2
1836 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
1837 C + IBARF(I),NREF(I),ANF(I)
1838 C IF(IORMO(I).NE.999)IORMO(I+NHAD1)=IORMO(I)+NHAD1
1839 C 41 CONTINUE
1840 C Intermediate store of hadrons from jet 2
1841  DO 2448 i=nhad1+1,nhad1+nhad2
1842  ii=i-nhad1
1843  anff(i)=anf(ii)
1844  pxff(i)=pxf(ii)
1845  pyff(i)=pyf(ii)
1846  pzff(i)=pzf(ii)
1847  heff(i)=hef(ii)
1848  amff(i)=amf(ii)
1849  ichff(i)=ichf(ii)
1850  ibarff(i)=ibarf(ii)
1851  nreff(i)=nref(ii)
1852  iormoo(i)=iormo(ii)
1853  IF(iormoo(i).NE.999)iormoo(i)=iormoo(i)+nhad1
1854  2448 CONTINUE
1855  nhad=nhad1+nhad2
1856  DO 2349 i=1,nhad
1857  ii=i
1858  anf(i)=anff(ii)
1859  pxf(i)=pxff(ii)
1860  pyf(i)=pyff(ii)
1861  pzf(i)=pzff(ii)
1862  hef(i)=heff(ii)
1863  amf(i)=amff(ii)
1864  ichf(i)=ichff(ii)
1865  ibarf(i)=ibarff(ii)
1866  nref(i)=nreff(ii)
1867  iormo(i)=iormoo(ii)
1868  2349 CONTINUE
1869  ENDIF
1870 C=================================================================
1871 C=================================================================
1872 C QUARK-DIQUARK NOBAM=4
1873 C=================================================================
1874 C==================================================================
1875 C Different options:
1876  IF(nobam.EQ.4)THEN
1877  IF(ipco.GE.0)WRITE(6,*)' QUARK-DIQUARK NOBAM=4'
1878 C NOBAM=4: Quark-DIQUARK Jet
1879 C
1880 C Work out approximate x-Values of quark and diquark
1881  IF(cmener.LT.1.d-3)THEN
1882  IF(ipco.GE.0)WRITE(6,*)' CMENER=0. HADJSE',cmener
1883  stop
1884  ENDIF
1885  xdiq=pta(4)*2.d0/cmener
1886  xqua=ppr(4)*2.d0/cmener
1887  IF(ipco.GE.0)THEN
1888  WRITE(6,'(A,2E12.4)')' HADJSE:XDIQ,XQUA ',xdiq,xqua
1889  ENDIF
1890 C Select valence quark x for one of the diquark quarks
1891 C Select x values of sea quark pair
1892  icou=0
1893  2234 CONTINUE
1894  icou=icou+1
1895  IF(icou.GE.200)THEN
1896  irej=1
1897  IF(isq.EQ.3)irej=3
1898  IF(ipco.GE.0)WRITE(6,*)' HADJSE Rejection 2234 ICOU. GT.100'
1899  RETURN
1900  ENDIF
1901  IF(ipco.GE.0)WRITE(6,*)' XSEAPA: CMENER,XQUA ',cmener,xqua
1902  CALL xseapa(cmener,xqua/2.d0,isq,isaq,xsq,xsaq,irej)
1903  iissqq=isq
1904  IF(irej.GE.1)THEN
1905  IF(ipco.GE.0)WRITE(6,*)' HADJSE reject XSEAPA'
1906  RETURN
1907  ENDIF
1908  IF(xsaq.GE.2.d0*xqua/3.d0)go to 2234
1909  IF(ipco.GE.0)THEN
1910  WRITE(6,'(A,2E12.4)')' HADJSE:XSQ,XSAQ ',xsq,xsaq
1911  ENDIF
1912 * target valence quarks xtvqi
1913 C here selected as sea quark! from 1/x-sea
1914  xvthro=cvq/cmener
1915  ivthr=0
1916  3466 CONTINUE
1917  IF(ivthr.EQ.100)THEN
1918  irej=1
1919  IF(isq.EQ.3)irej=3
1920  IF(ipco.GE.0)WRITE(6,*)' HADJSE 3466 reject IVTHR 50'
1921  RETURN
1922  ENDIF
1923  ivthr=ivthr+1
1924  xvthr=xvthro/(101-ivthr)
1925  unoprv=unon
1926  380 CONTINUE
1927  IF(xvthr.GT.0.66d0*xdiq)THEN
1928  irej=1
1929  IF(isq.EQ.3)irej=3
1930  IF(ipco.GE.0)WRITE(6,*)' HADJSE Rejection 380 XVTHR large ',
1931  * xvthr
1932  RETURN
1933  ENDIF
1934  xtvqi=sampey(xvthr,0.66d0*xdiq)
1935  xdiqq=xdiq-xtvqi
1936  IF(ipco.GE.0)THEN
1937  WRITE(6,'(A,2E12.4)')' HADJCK:XDIQQ,XTVQI ',xdiqq,xtvqi
1938  ENDIF
1939 C We form two strings(2-1) aq-q and (4-3) q-qq
1940 C 1 q with XDIQ-XTVQI-XSQ
1941 C 1 q with XDIQ-XTVQI (j.r.5.5.96)
1942 C 2 aq with XSAQ
1943 C 3 qq with XTVQI+XSQ
1944 C 3 qq with XTVQI (j.r.5.5.96)
1945 C 4 q with XQUA-XSAQ
1946 C with 4-momenta PPP1,PPP2,PPP3,PPP4
1947  DO 3345 i=1,4
1948  IF(xdiq.LE.1.d-15.OR.xqua.LT.1.d-15)THEN
1949  irej=1
1950  IF(isq.EQ.3)irej=3
1951  IF(ipco.GE.0)WRITE (6,*)' HSDJSE Rejection 3345 XDIQ,XQUA ',
1952  * xdiq,xqua
1953  RETURN
1954  ENDIF
1955 C PPP1(I)=PTA(I)*(XDIQ-XTVQI-XSQ)/XDIQ
1956  ppp1(i)=pta(i)*(xdiq-xtvqi)/xdiq
1957 C PPP3(I)=PTA(I)*(XTVQI+XSQ)/XDIQ
1958  ppp3(i)=pta(i)*(xtvqi)/xdiq
1959  ppp2(i)=ppr(i)*xsaq/xqua
1960  ppp4(i)=ppr(i)*(xqua-xsaq)/xqua
1961  3345 CONTINUE
1962  3346 FORMAT(a,5e12.4)
1963  IF(ipco.GE.0)THEN
1964  WRITE(6,3346)' PPR ',ppr
1965  WRITE(6,3346)' PPP1 ',ppp1
1966  WRITE(6,3346)' PPP3 ',ppp3
1967  WRITE(6,3346)' PTA ',pta
1968  WRITE(6,3346)' PPP2 XSAQ',ppp2,xsaq
1969  WRITE(6,3346)' PPP4 ',ppp4
1970  ENDIF
1971 C Invariant Masses of chains
1972 C ORIGINAL CHAIN
1973  amchor=sqrt((ppr(4)+pta(4))**2-(ppr(3)+pta(3))**2
1974  * -(ppr(2)+pta(2))**2-(ppr(1)+pta(1))**2)
1975  IF(ipco.GE.0)WRITE(6,2346)'AMCHOR ',amchor
1976  amchn1=sqrt((ppp1(4)+ppp2(4))**2-(ppp1(3)+ppp2(3))**2
1977  * -(ppp1(2)+ppp2(2))**2-(ppp1(1)+ppp2(1))**2)
1978  IF(ipco.GE.0)WRITE(6,2346)'AMCHN1 ',amchn1
1979 C Chain 1 is q-aq restrict mass >1.5 GeV
1980  IF(amchor.LE.tiny)THEN
1981  IF(ipco.GE.0)WRITE(6,2346)' PPR ',ppr
1982  WRITE(6,2346)' PPP1 ',ppp1
1983  WRITE(6,2346)' PPP3 ',ppp3
1984  WRITE(6,2346)' PTA ',pta
1985  WRITE(6,2346)' PPP2 ',ppp2
1986  WRITE(6,2346)' PPP4 ',ppp4
1987  WRITE(6,2346)'AMCHOR ',amchor
1988  ENDIF
1989  IF(amchn1.LE.tiny)THEN
1990  WRITE(6,2346)' PPR ',ppr
1991  WRITE(6,2346)' PPP1 ',ppp1
1992  WRITE(6,2346)' PPP3 ',ppp3
1993  WRITE(6,2346)' PTA ',pta
1994  WRITE(6,2346)' PPP2 ',ppp2
1995  WRITE(6,2346)' PPP4 ',ppp4
1996  WRITE(6,2346)'AMCHOR ',amchor
1997  ENDIF
1998  chamal=0.8d0
1999  IF(ifb3.GE.3.OR.isaq.GE.9)chamal=1.2d0
2000  IF(amchn1.LE.chamal)THEN
2001  IF(ipco.GE.0)WRITE(6 ,*).LE.'HADJSE jump2AMCHN1CHAMAL AMCHOR',
2002  * amchn1,chamal,amchor
2003  go to 3466
2004  ENDIF
2005  amchn2=sqrt((ppp3(4)+ppp4(4))**2-(ppp3(3)+ppp4(3))**2
2006  * -(ppp3(2)+ppp4(2))**2-(ppp3(1)+ppp4(1))**2)
2007 C Chain 2 is qq-q restrict mass >2.5 GeV
2008  chamal=1.5d0
2009  IF(ifb2.GE.3.OR.ifb1.GE.3.OR.isq.GE.3)chamal=2.0d0
2010  IF(amchn2.LE.chamal)THEN
2011 C IREJ=1
2012 C RETURN
2013  IF(ipco.GE.0)WRITE(6 ,*).LE.' HADJSE jump AMCHN2CHAMAL',
2014  * amchn2,chamal
2015  go to 3466
2016  ENDIF
2017  pxchk=ppr(1)+pta(1)-ppp1(1)-ppp2(1)-ppp3(1)-ppp4(1)
2018  pychk=ppr(2)+pta(2)-ppp1(2)-ppp2(2)-ppp3(2)-ppp4(2)
2019  pzchk=ppr(3)+pta(3)-ppp1(3)-ppp2(3)-ppp3(3)-ppp4(3)
2020  pechk=ppr(4)+pta(4)-ppp1(4)-ppp2(4)-ppp3(4)-ppp4(4)
2021  IF(ipco.GE.0)THEN
2022  WRITE(6,'(A/8E12.4,I5)')
2023  * ' Chain masses AMCH,AMCHOR,AMCHN1,AMCHN2,PZCHK,PECHK,
2024  * PXCHK,PYCHK,NOBAM ',
2025  * amch,amchor,amchn1,amchn2,pzchk,pechk,pxchk,pychk,nobam
2026  ENDIF
2027 C Lorentz parameters of chains
2028  IF(amch.LE.1.d-15)THEN
2029  irej=1
2030  IF(isq.EQ.3)irej=3
2031  IF(ipco.GE.0)WRITE(6,*)' HSDJSE Rejection AMCH too small',amch
2032  RETURN
2033  ENDIF
2034  gamor=(ppr(4)+pta(4))/amch
2035  bgxor=(ppr(1)+pta(1))/amch
2036  bgyor=(ppr(2)+pta(2))/amch
2037  bgzor=(ppr(3)+pta(3))/amch
2038  IF(amchn1.LE.1.d-15)THEN
2039  IF(ipco.GE.0)WRITE(6,*)' HSDJSE Rejection AMCHN1 too small',
2040  * amchn1
2041  irej=1
2042  IF(isq.EQ.3)irej=3
2043  RETURN
2044  ENDIF
2045  gamch1=(ppp1(4)+ppp2(4))/amchn1
2046  bgxch1=(ppp1(1)+ppp2(1))/amchn1
2047  bgych1=(ppp1(2)+ppp2(2))/amchn1
2048  bgzch1=(ppp1(3)+ppp2(3))/amchn1
2049  IF(amchn2.LE.1.d-15)THEN
2050  irej=1
2051  IF(isq.EQ.3)irej=3
2052  IF(ipco.GE.0)WRITE(6,*)' HSDJSE Rejection AMCHN2 too small',
2053  * amchn2
2054  RETURN
2055  ENDIF
2056  gamch2=(ppp3(4)+ppp4(4))/amchn2
2057  bgxch2=(ppp3(1)+ppp4(1))/amchn2
2058  bgych2=(ppp3(2)+ppp4(2))/amchn2
2059  bgzch2=(ppp3(3)+ppp4(3))/amchn2
2060  IF(ipco.GE.0)THEN
2061  WRITE(6,3346)' L.Parm in ',gam,bgx,bgy,bgz
2062  WRITE(6,3346)' L.Parm OR ',gamor,bgxor,bgyor,bgzor
2063  WRITE(6,3346)' L.Parm C1 ',gamch1,bgxch1,bgych1,bgzch1
2064  WRITE(6,3346)' L.Parm C2 ',gamch2,bgxch2,bgych2,bgzch2
2065  ENDIF
2066  peor=amch*gamor
2067  pzor=amch*bgzor
2068  pech1=amchn1*gamch1
2069  pzch1=amchn1*bgzch1
2070  pech2=amchn2*gamch2
2071  pzch2=amchn2*bgzch2
2072  IF(ipco.GE.0)THEN
2073  WRITE(6,'(A,6E12.4)')' PEOR,PECH1,PECH2,PZOR,PZCH1,PZCH2',
2074  * peor,pech1,pech2,pzor,pzch1,pzch2
2075  ENDIF
2076  noba1=3
2077  noba2=4
2078  IF(ipco.GE.0)THEN
2079  WRITE(6,'(A,2I5)')' Jet1 ISAQ,IFB3 ',isaq,ifb3
2080  ENDIF
2081  CALL hadjet(nhad1,amchn1,ppp2,ppp1,gamch1,bgxch1,bgych1,
2082  * bgzch1, isaq,ifb3,ifb1,ifb4,i1,i2,noba1,nnch,norig)
2083 C DO 342 I=1,NHAD1
2084 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
2085 C + IBARF(I),NREF(I),ANF(I)
2086 C 342 CONTINUE
2087 C DIMENSION ANFF(NFIMAX),PXFF(NFIMAX),PYFF(NFIMAX),PZFF(NFIMAX),
2088 C + HEFF(NFIMAX),AMFF(NFIMAX),
2089 C * ICHFF(NFIMAX),IBARFF(NFIMAX),NREFF(NFIMAX)
2090 C Intermediate store of hadrons from jet 1
2091  DO 3348 i=1,nhad1
2092  anff(i)=anf(i)
2093  pxff(i)=pxf(i)
2094  pyff(i)=pyf(i)
2095  pzff(i)=pzf(i)
2096  heff(i)=hef(i)
2097  amff(i)=amf(i)
2098  ichff(i)=ichf(i)
2099  ibarff(i)=ibarf(i)
2100  nreff(i)=nref(i)
2101  iormoo(i)=iormo(i)
2102  3348 CONTINUE
2103  IF(ipco.GE.0)THEN
2104  WRITE(6,'(A,3I5)')' Jet2 IFB1,ISQ,IFB2 ',ifb1,isq,ifb2
2105  ENDIF
2106  CALL hadjet(nhad2,amchn2,ppp4,ppp3,gamch2,bgxch2,bgych2,
2107  * bgzch2, ifb1,isq,ifb2,ifb4,i1,i2,noba2,nnch,norig)
2108 C DO 341 I=1,NHAD2
2109 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
2110 C + IBARF(I),NREF(I),ANF(I)
2111 C IF(IORMO(I).NE.999)IORMO(I+NHAD1)=IORMO(I)+NHAD1
2112 C 341 CONTINUE
2113 C Intermediate store of hadrons from jet 2
2114  DO 3448 i=nhad1+1,nhad1+nhad2
2115  ii=i-nhad1
2116  anff(i)=anf(ii)
2117  pxff(i)=pxf(ii)
2118  pyff(i)=pyf(ii)
2119  pzff(i)=pzf(ii)
2120  heff(i)=hef(ii)
2121  amff(i)=amf(ii)
2122  ichff(i)=ichf(ii)
2123  ibarff(i)=ibarf(ii)
2124  nreff(i)=nref(ii)
2125  iormoo(i)=iormo(ii)
2126  IF(iormoo(i).NE.999)iormoo(i)=iormoo(i)+nhad1
2127  3448 CONTINUE
2128  nhad=nhad1+nhad2
2129  DO 3349 i=1,nhad
2130  ii=i
2131  anf(i)=anff(ii)
2132  pxf(i)=pxff(ii)
2133  pyf(i)=pyff(ii)
2134  pzf(i)=pzff(ii)
2135  hef(i)=heff(ii)
2136  amf(i)=amff(ii)
2137  ichf(i)=ichff(ii)
2138  ibarf(i)=ibarff(ii)
2139  nref(i)=nreff(ii)
2140  iormo(i)=iormoo(ii)
2141  3349 CONTINUE
2142  ENDIF
2143 C==================================================================
2144  heft=0.d0
2145  IF(ipco.GE.0)THEN
2146  DO 40 i=1,nhad
2147  IF(ibarf(i).EQ.500)go to 4040
2148  heft=heft+hef(i)
2149  WRITE(6,4050)i,pxf(i),pyf(i),pzf(i),hef(i),amf(i), ichf(i),
2150  + ibarf(i),nref(i),anf(i)
2151  4050 FORMAT(' JET ',i5,5f12.4,3i5,a10)
2152  4040 CONTINUE
2153  40 CONTINUE
2154  ENDIF
2155  hefff=amch*gam
2156  IF(ipco.GE.0)THEN
2157  WRITE(6,'(A,I5,2E12.4)')' HADJSE 2IREJ,HEFT,HEFFF ',
2158  * irej,heft,hefff
2159  ENDIF
2160 C IREJ=1
2161  IF(irej.GE.1)RETURN
2162 C
2163 C WRITE(6,*)' HADJSE: IREJ,ISQ ',IREJ,ISQ
2164  IF(isq.EQ.1)nhse1=nhse1+1
2165  IF(isq.EQ.2)nhse2=nhse2+1
2166  IF(isq.EQ.3)nhse3=nhse3+1
2167  RETURN
2168  END
2169 
2170 ************************************************************************
2171 ************************************************************************
2172 ************************************************************************
2173 ************************************************************************
2174 C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2175 C
2176  SUBROUTINE hadjase(NHAD,AMCH,PPR,PTA,GAM,BGX,BGY,BGZ, IFB1,IFB2,
2177  +ifb3,ifb4,i1,i2,nobam,nnch,norig,irej,iissqq)
2178  IMPLICIT DOUBLE PRECISION (a-h,o-z)
2179  SAVE
2180 C
2181 C HADJASE Fragmentation of vv,vs and sv chains CK method
2182 C using Glauber sea quarks
2183 C This times for aqaq--aq and aq--aqaq chains
2184 C j.r.3/99
2185 C
2186 C HADJET DOES ALL THE NECESSARY ROTATIONS AND LORENTZ TRANSFORMS AND
2187 C CALL CALBAM (BAMJET)
2188 C
2189 C NHAD = NUMBER OF HADRONS CREATED (OUTPUT) = IHAD (CALBAM)
2190 C******** PRODUCED PARTICLES IN COMMON /CMSRES/
2191 C NOTE: NOW IN /FINPAR/ HJM 21/06/90
2192 C AMCH = INVARIANT MASS OF JET (INPUT)
2193 C PPR = 4-MOMENTUM OF FORWARD GOING PARTON (PROJECTILE)(INPUT)
2194 C PTA = 4-MOMENTUM OF BACKWARD GOING PARTON (TARGET)(INPUT)
2195 C GAM,BGX,BGY,BGZ = LORENTZ PARAMETERS OF JET CMS RELATIVE TO
2196 C COLLISION CMS (INPUT)
2197 C
2198 C--------------------------------------------------------------------
2199 C CALBAM(NNCH,I1,I2,IFB1,IFB2,IFB3,IFB4,AMCH,NOBAM,IHAD)
2200 C SAMPLING OF Q-AQ, Q-QQ, QQ-AQQ CHAINS
2201 C USING BAMJET(IHAD,IFB1,IFB2,IFB3,IFB4,AMCH,NOBAM)-----FOR NNCH=0
2202 C OR PARJET(IHAD,ICH1=I1 OR I2)------FOR NNCH=-1 OR +1
2203 C-------------------------------------------------------------------
2204 C IHAD : NUMBER OF PRODUCED PARTICLES
2205 C NNCH : CALL BAMJET FOR NNCH=0
2206 C CALL PARJET FOR NNCH=+1 ICH1=I1
2207 C FOR NNCH=-1 ICH1=I2
2208 C jet not existing for NNCH=+/-99, i.e. IHAD=0
2209 C PRODUCED PARTICLES IN CHAIN REST FRAME ARE IN COMMON /FINPAR/
2210 C AMCH : INVARIANT MASS OF CHAIN (BAMJET)
2211 C
2212 C NOBAM : = 3 QUARK-ANTIQUARK JET QUARK FLAVORS : IFB1,IFB2
2213 C OR ANTIQUARK-QUARK JET IN ANY ORDER
2214 C
2215 C = 4 QUARK-DIQUARK JET, FLAVORS : aQU : IFB1, aDIQU :IFB2,IFB3
2216 C OR ANTIQUARK-ANTIDIQUARK JET
2217 C
2218 C
2219 C = 5 DIQUARK-ANTIDIQUARK JET
2220 C OR ANTIDIQUARK-DIQUARK JET
2221 C FLAVORS : DIQU : IFB1,IFB2, ANTIDIQU : IFB3,IFB4
2222 C IN ANY ORDER
2223 C
2224 C = 6 DIQUARK-QUARK JET, FLAVORS : DIQU : IFB1,IFB2 QU: IFB3
2225 C OR ANTIDIQUARK-ANTIQUARK JET
2226 C
2227 C IFBI : FLAVORS : 1,2,3,4 = U,D,S,C 7,8,9,10 = AU,AD,AS,AC
2228 C
2229 C I1,I2 : NUMBER LABEL OF A HADRON CREATED BY PARJET
2230 C
2231 C NORMALLY IN BAMJET THE QUARKS MOVE FORWARD (POSITIVE Z-DIRECTION)
2232 C THE QUARK FLAVORS ARE FIRST GIVEN
2233 C CALBAM ALLOWS EITHER THE QUARK OR ANTIQUARK (DIQU) TO MOVE FORWARD
2234 C THE FORWARD GOING FLAVORS ARE GIVEN FIRST
2235 C
2236 C =====================================================================
2237 *KEEP,DFINPA.
2238  CHARACTER*8 anf,anff
2239  parameter(nfimax=249)
2240  COMMON /dfinpa/ anf(nfimax),pxf(nfimax),pyf(nfimax),pzf(nfimax),
2241  +hef(nfimax),amf(nfimax), ichf(nfimax),ibarf(nfimax),nref(nfimax)
2242  COMMON /dfinpz/iormo(nfimax),idaug1(nfimax),idaug2(nfimax),
2243  * istath(nfimax)
2244  dimension anff(nfimax),pxff(nfimax),pyff(nfimax),pzff(nfimax),
2245  +heff(nfimax),amff(nfimax),
2246  *ichff(nfimax),ibarff(nfimax),nreff(nfimax),iormoo(nfimax)
2247  CHARACTER*80 title
2248  CHARACTER*8 projty,targty
2249 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
2250 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
2251  COMMON /user1/title,projty,targty
2252  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
2253 *KEEP,DPRIN.
2254  COMMON /dprin/ ipri,ipev,ippa,ipco,init,iphkk,itopd,ipaupr
2255 *KEND.
2256 *--------------------------------------------------------------------
2257 *-------------------------------------------------------------------
2258  COMMON /jspart/pxp(1000),pyp(1000),pzp(1000),hepp(1000),nnnp
2259  COMMON /jspar/pxj(1000),pyj(1000),pzj(1000),hej(1000),nnnpj
2260 ************************************************************************
2261 ************************************************************************
2262  common/ifragm/ifrag
2263  dimension ppr(4),pta(4)
2264 *KEEP,XSEADI.
2265  COMMON /xseadi/ xseacu,unon,unom,unosea,cvq,cdq,csea,ssmima,
2266  +ssmimq,vvmthr
2267  COMMON /hdjase/nhse1,nhse2,nhse3,nhase1,nhase2,nhase3
2268  dimension ppp1(4),ppp2(4),ppp3(4),ppp4(4)
2269  DATA tiny/1.0d-3/
2270 C----------------------------------------------------------------------
2271 C
2272 C WRITE(6,'(A,3I5)')' Jet0 IFB1,IFB2,IFB3 ',IFB1,IFB2,IFB3
2273  IF(ipco.GE.0)WRITE(6,*)
2274  *' HADJASE(NHAD,AMCH,PPR,PTA,GAM,BGX,BGY,BGZ, IFB1,IFB2,',
2275  +'IFB3,IFB4,I1,I2,NOBAM,NNCH,NORIG,IREJ),IPCO',
2276  *nhad,amch,ppr,pta,gam,bgx,bgy,bgz, ifb1,ifb2,
2277  +ifb3,ifb4,i1,i2,nobam,nnch,norig,irej,ipco
2278  irej=0
2279  IF(ipco.GE.0) THEN
2280  WRITE(6,'(A,3I5)')' HADJASE Jet0 IFB1,IFB2,IFB3',ifb1,ifb2,ifb3
2281  WRITE(6,'(4E12.4)') gam,bgx,bgy,bgz,ppr
2282  WRITE(6,1000) nhad,amch,pta, ifb1,ifb2,ifb3,ifb4,i1,i2,nobam,
2283  + nnch,norig
2284  1000 FORMAT(10x,i10,5f10.3/10x,9i10)
2285  ENDIF
2286  IF(abs(nnch).EQ.99) THEN
2287  nhad=0
2288 C IPCO=0
2289  RETURN
2290  ENDIF
2291 C=================================================================
2292 C=================================================================
2293 C DIQUARK-QUARK NOBAM=6
2294 C=================================================================
2295 C=================================================================
2296 C Different options:
2297  IF(nobam.EQ.6)THEN
2298  IF(ipco.GE.0)WRITE(6,*)' DIQUARK-QUARK NOBAM=6'
2299 C NOBAM=6: Diquark--Quark Jet
2300 C
2301 C Work out approximate x-Values of diquark and quark
2302  IF(cmener.LT.1.d-3)THEN
2303  WRITE(6,*)' CMENER=0. HADJASE',cmener
2304  stop
2305  ENDIF
2306  xdiq=ppr(4)*2.d0/cmener
2307  xqua=pta(4)*2.d0/cmener
2308  IF(ipco.GE.0) THEN
2309  WRITE(6,'(A,2E12.4)')' HADJASE:XDIQ,XQUA ',xdiq,xqua
2310  ENDIF
2311 C Select sea quark x for one of the diquark quarks
2312 C Select x values of sea quark pair
2313  icou=0
2314  1234 CONTINUE
2315  icou=icou+1
2316  IF(icou.GE.500)THEN
2317  irej=1
2318  IF(ipco.GE.0)WRITE(6,*)' HADJASE reject icou 100'
2319  RETURN
2320  ENDIF
2321  CALL xseapa(cmener,xqua/2.d0,isq,isaq,xsq,xsaq,irej)
2322  iissqq=isq
2323  IF(irej.GE.1)THEN
2324  IF(ipco.GE.0)WRITE(6,*)' HADJASE reject XSEAPA'
2325  RETURN
2326  ENDIF
2327  IF(xsaq.GE.2.d0*xqua/3.d0)go to 1234
2328  IF(ipco.GE.0) THEN
2329  WRITE(6,*)' HADJASE,XSEAPA:XSQ,XSAQ,ISQ,ISAQ ',
2330  * xsq,xsaq,isq,isaq
2331  ENDIF
2332 * projectile valence quarks xpvqi
2333 C here selected as sea quark! from 1/x-sea
2334  xvthro=cvq/cmener
2335  ivthr=0
2336  3465 CONTINUE
2337  IF(ivthr.EQ.200)THEN
2338  irej=1
2339  IF(isq.EQ.3)irej=3
2340  IF(ipco.GE.0)WRITE(6,*)' HADJASE 3465 reject IVTHR 50'
2341  RETURN
2342  ENDIF
2343  ivthr=ivthr+1
2344  xvthr=xvthro/(201-ivthr)
2345  unoprv=unon
2346  80 CONTINUE
2347  IF(xvthr.GT.0.66d0*xdiq)THEN
2348  IF(ipco.GE.0)WRITE(6,*)' HADJASE 80 reject XVTHR too great',
2349  * xvthr
2350  irej=1
2351  IF(isq.EQ.3)irej=3
2352  RETURN
2353  ENDIF
2354  xpvqi=sampey(xvthr,0.66d0*xdiq)
2355  xdiqq=xdiq-xpvqi
2356  IF(ipco.GE.0) THEN
2357  WRITE(6,'(A,2E12.4)')' HADJASE:XDIQQ,XPVQI ',xdiqq,xpvqi
2358  ENDIF
2359 C
2360 C We form two strings(1-2) aq-q and (3-4) aqaq-aq
2361 C 1 aq with XDIQ-XPVQI (j.r.5.5.96)
2362 C 2 q with XSAQ
2363 C 3 aqaq with XPVQI (j.r.5.5.96)
2364 C 4 aq with XQUA-XSAQ
2365 C with 4-momenta PPP1,PPP2,PPP3,PPP4
2366  DO 2345 i=1,4
2367  IF(xdiq.LE.1.d-15.OR.xqua.LT.1.d-15)THEN
2368  irej=1
2369  IF(isq.EQ.3)irej=3
2370  IF(ipco.GE.0)WRITE(6,*)' HADJASE reject 2345 XDIQ,',
2371  * 'XQUA too small ',
2372  * xdiq,xqua
2373  RETURN
2374  ENDIF
2375  ppp1(i)=ppr(i)*(xdiq-xpvqi)/xdiq
2376  ppp3(i)=ppr(i)*(xpvqi)/xdiq
2377  ppp2(i)=pta(i)*xsaq/xqua
2378  ppp4(i)=pta(i)*(xqua-xsaq)/xqua
2379  2345 CONTINUE
2380  2346 FORMAT(a,5e12.4)
2381  IF(ipco.GE.0) THEN
2382  WRITE(6,2346)' PPR ',ppr
2383  WRITE(6,2346)' PPP1 ',ppp1
2384  WRITE(6,2346)' PPP3 ',ppp3
2385  WRITE(6,2346)' PTA ',pta
2386  WRITE(6,2346)' PPP2 XSAQ',ppp2,xsaq
2387  WRITE(6,2346)' PPP4 ',ppp4
2388  ENDIF
2389 C Invariant Masses of chains
2390 C ORIGINAL CHAIN
2391  amchor=sqrt((ppr(4)+pta(4))**2-(ppr(3)+pta(3))**2
2392  * -(ppr(2)+pta(2))**2-(ppr(1)+pta(1))**2)
2393  amchn1=sqrt((ppp1(4)+ppp2(4))**2-(ppp1(3)+ppp2(3))**2
2394  * -(ppp1(2)+ppp2(2))**2-(ppp1(1)+ppp2(1))**2)
2395  IF(amchor.LE.tiny)THEN
2396  WRITE(6,2346)' PPR ',ppr
2397  WRITE(6,2346)' PPP1 ',ppp1
2398  WRITE(6,2346)' PPP3 ',ppp3
2399  WRITE(6,2346)' PTA ',pta
2400  WRITE(6,2346)' PPP2 ',ppp2
2401  WRITE(6,2346)' PPP4 ',ppp4
2402  ENDIF
2403  IF(amchn1.LE.tiny)THEN
2404  WRITE(6,2346)' PPR ',ppr
2405  WRITE(6,2346)' PPP1 ',ppp1
2406  WRITE(6,2346)' PPP3 ',ppp3
2407  WRITE(6,2346)' PTA ',pta
2408  WRITE(6,2346)' PPP2 ',ppp2
2409  WRITE(6,2346)' PPP4 ',ppp4
2410  ENDIF
2411 C Chain 1 is q-aq restrict mass >1.5 GeV
2412  chamal=0.8d0
2413  IF(ifb1.GE.9.OR.isq.GE.3)chamal=1.1d0
2414  IF(amchn1.LE.chamal)THEN
2415  IF(ipco.GE.0)WRITE (6,*).LE.'HADJASE jump1AMCHN1CHAMAL AMCHOR'
2416  * ,amchn1,chamal,amchor
2417  go to 3465
2418  ENDIF
2419  amchn2=sqrt((ppp3(4)+ppp4(4))**2-(ppp3(3)+ppp4(3))**2
2420  * -(ppp3(2)+ppp4(2))**2-(ppp3(1)+ppp4(1))**2)
2421 C Chain 2 is qq-q restrict mass >2.5 GeV
2422  chamal=1.3d0
2423  IF(ifb2.GE.9.OR.ifb3.GE.9.OR.isaq.GE.9)chamal=1.80d0
2424  IF(amchn2.LE.chamal)THEN
2425 C IREJ=1
2426 C RETURN
2427  IF(ipco.GE.0)WRITE (6,*).LE.'HADJASE jump AMCHN2CHAMAL ',
2428  * amchn2,chamal
2429  go to 3465
2430  ENDIF
2431  pxchk=ppr(1)+pta(1)-ppp1(1)-ppp2(1)-ppp3(1)-ppp4(1)
2432  pychk=ppr(2)+pta(2)-ppp1(2)-ppp2(2)-ppp3(2)-ppp4(2)
2433  pzchk=ppr(3)+pta(3)-ppp1(3)-ppp2(3)-ppp3(3)-ppp4(3)
2434  pechk=ppr(4)+pta(4)-ppp1(4)-ppp2(4)-ppp3(4)-ppp4(4)
2435  IF(ipco.GE.0)THEN
2436  WRITE(6,'(A/8E12.4,I5)')
2437  * ' Chain masses AMCH,AMCHOR,AMCHN1,AMCHN2,PZCHK,PECHK,
2438  * PXCHK,PYCHK,NOBAM ',
2439  * amch,amchor,amchn1,amchn2,pzchk,pechk,pxchk,pychk,nobam
2440  ENDIF
2441 C Lorentz parameters of chains
2442  IF(amch.LE.1.d-15)THEN
2443  irej=1
2444  IF(isq.EQ.3)irej=3
2445  IF(ipco.GE.0)WRITE(6,*)' HADJASE rejection AMCH too small ',
2446  * amch
2447  RETURN
2448  ENDIF
2449  gamor=(ppr(4)+pta(4))/amch
2450  bgxor=(ppr(1)+pta(1))/amch
2451  bgyor=(ppr(2)+pta(2))/amch
2452  bgzor=(ppr(3)+pta(3))/amch
2453  IF(amchn1.LE.1.d-15)THEN
2454  IF(ipco.GE.0)WRITE(6,*)' HADJASE rejection AMCHN1 too small ',
2455  * amchn1
2456  irej=1
2457  IF(isq.EQ.3)irej=3
2458  RETURN
2459  ENDIF
2460  gamch1=(ppp1(4)+ppp2(4))/amchn1
2461  bgxch1=(ppp1(1)+ppp2(1))/amchn1
2462  bgych1=(ppp1(2)+ppp2(2))/amchn1
2463  bgzch1=(ppp1(3)+ppp2(3))/amchn1
2464  IF(amchn2.LE.1.d-15)THEN
2465  IF(ipco.GE.0)WRITE(6,*)' HADJASE rejection AMCHN2 too small ',
2466  * amchn2
2467  irej=1
2468  IF(isq.EQ.3)irej=3
2469  RETURN
2470  ENDIF
2471  gamch2=(ppp3(4)+ppp4(4))/amchn2
2472  bgxch2=(ppp3(1)+ppp4(1))/amchn2
2473  bgych2=(ppp3(2)+ppp4(2))/amchn2
2474  bgzch2=(ppp3(3)+ppp4(3))/amchn2
2475  IF(ipco.GE.0)THEN
2476  WRITE(6,2346)' L.Parm in ',gam,bgx,bgy,bgz
2477  WRITE(6,2346)' L.Parm OR ',gamor,bgxor,bgyor,bgzor
2478  WRITE(6,2346)' L.Parm C1 ',gamch1,bgxch1,bgych1,bgzch1
2479  WRITE(6,2346)' L.Parm C2 ',gamch2,bgxch2,bgych2,bgzch2
2480  ENDIF
2481  peor=amch*gamor
2482  pzor=amch*bgzor
2483  pech1=amchn1*gamch1
2484  pzch1=amchn1*bgzch1
2485  pech2=amchn2*gamch2
2486  pzch2=amchn2*bgzch2
2487  IF(ipco.GE.0)THEN
2488  WRITE(6,'(A,6E12.4)')' PEOR,PECH1,PECH2,PZOR,PZCH1,PZCH2',
2489  * peor,pech1,pech2,pzor,pzch1,pzch2
2490  ENDIF
2491  noba1=3
2492  noba2=6
2493  IF(ipco.GE.-1)THEN
2494 C WRITE(6,'(A,2I5)')' Jet1 IFB1,ISQ ',IFB1,ISQ
2495  ENDIF
2496  CALL hadjet(nhad1,amchn1,ppp1,ppp2,gamch1,bgxch1,bgych1,
2497  * bgzch1, ifb1,isq,ifb3,ifb4,i1,i2,noba1,nnch,norig)
2498 C DO 42 I=1,NHAD1
2499 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
2500 C + IBARF(I),NREF(I),ANF(I)
2501 C 42 CONTINUE
2502  IF(ipco.GE.-1)THEN
2503 C WRITE(6,'(A,2I5)')' Jet1 NHAD1 ',NHAD1
2504  ENDIF
2505 C DIMENSION ANFF(NFIMAX),PXFF(NFIMAX),PYFF(NFIMAX),PZFF(NFIMAX),
2506 C + HEFF(NFIMAX),AMFF(NFIMAX),
2507 C * ICHFF(NFIMAX),IBARFF(NFIMAX),NREFF(NFIMAX)
2508 C Intermediate store of hadrons from jet 1
2509  DO 2348 i=1,nhad1
2510  anff(i)=anf(i)
2511  pxff(i)=pxf(i)
2512  pyff(i)=pyf(i)
2513  pzff(i)=pzf(i)
2514  heff(i)=hef(i)
2515  amff(i)=amf(i)
2516  ichff(i)=ichf(i)
2517  ibarff(i)=ibarf(i)
2518  nreff(i)=nref(i)
2519  iormoo(i)=iormo(i)
2520  2348 CONTINUE
2521  IF(ipco.GE.0)THEN
2522  WRITE(6,'(A,3I5)')' Jet2 IFB2,ISAQ,IFB3 ',ifb2,isaq,ifb3
2523  ENDIF
2524  CALL hadjet(nhad2,amchn2,ppp3,ppp4,gamch2,bgxch2,bgych2,
2525  * bgzch2, ifb2,isaq,ifb3,ifb4,i1,i2,noba2,nnch,norig)
2526  IF(ipco.GE.-1)THEN
2527 C WRITE(6,'(A,2I5)')' Jet2 NHAD2 ',NHAD2
2528  ENDIF
2529 C DO 41 I=1,NHAD2
2530 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
2531 C + IBARF(I),NREF(I),ANF(I)
2532 C IF(IORMO(I).NE.999)IORMO(I+NHAD1)=IORMO(I)+NHAD1
2533 C 41 CONTINUE
2534 C Intermediate store of hadrons from jet 2
2535  DO 2448 i=nhad1+1,nhad1+nhad2
2536  ii=i-nhad1
2537  anff(i)=anf(ii)
2538  pxff(i)=pxf(ii)
2539  pyff(i)=pyf(ii)
2540  pzff(i)=pzf(ii)
2541  heff(i)=hef(ii)
2542  amff(i)=amf(ii)
2543  ichff(i)=ichf(ii)
2544  ibarff(i)=ibarf(ii)
2545  nreff(i)=nref(ii)
2546  iormoo(i)=iormo(ii)
2547  IF(iormoo(i).NE.999)iormoo(i)=iormoo(i)+nhad1
2548  2448 CONTINUE
2549  nhad=nhad1+nhad2
2550  DO 2349 i=1,nhad
2551  ii=i
2552  anf(i)=anff(ii)
2553  pxf(i)=pxff(ii)
2554  pyf(i)=pyff(ii)
2555  pzf(i)=pzff(ii)
2556  hef(i)=heff(ii)
2557  amf(i)=amff(ii)
2558  ichf(i)=ichff(ii)
2559  ibarf(i)=ibarff(ii)
2560  nref(i)=nreff(ii)
2561  iormo(i)=iormoo(ii)
2562  2349 CONTINUE
2563  ENDIF
2564 C=================================================================
2565 C=================================================================
2566 C QUARK-DIQUARK NOBAM=4
2567 C=================================================================
2568 C==================================================================
2569 C Different options:
2570  IF(nobam.EQ.4)THEN
2571  IF(ipco.GE.0)WRITE(6,*)' AQUARK-ADIQUARK NOBAM=4'
2572 C NOBAM=4: AQuark-ADIQUARK Jet
2573 C
2574 C Work out approximate x-Values of quark and diquark
2575  IF(cmener.LT.1.d-3)THEN
2576  IF(ipco.GE.0)WRITE(6,*)' CMENER=0. HADJASE',cmener
2577  stop
2578  ENDIF
2579  xdiq=pta(4)*2.d0/cmener
2580  xqua=ppr(4)*2.d0/cmener
2581  IF(ipco.GE.0)THEN
2582  WRITE(6,'(A,2E12.4)')' HADJASE:XDIQ,XQUA ',xdiq,xqua
2583  ENDIF
2584 C Select valence quark x for one of the diquark quarks
2585 C Select x values of sea quark pair
2586  icou=0
2587  2234 CONTINUE
2588  icou=icou+1
2589  IF(icou.GE.500)THEN
2590  irej=1
2591  IF(isq.EQ.3)irej=3
2592  IF(ipco.GE.0)WRITE(6,*)' HADJASE Rejection 2234 ICOU. GT.100'
2593  RETURN
2594  ENDIF
2595  IF(ipco.GE.0)WRITE(6,*)' XSEAPA: CMENER,XQUA ',cmener,xqua
2596  CALL xseapa(cmener,xqua/2.d0,isq,isaq,xsq,xsaq,irej)
2597  iissqq=isq
2598  IF(irej.GE.1)THEN
2599  IF(ipco.GE.0)WRITE(6,*)' HADJASE reject XSEAPA'
2600  RETURN
2601  ENDIF
2602  IF(xsaq.GE.2.d0*xqua/3.d0)go to 2234
2603  IF(ipco.GE.0)THEN
2604  WRITE(6,'(A,2E12.4)')' HADJASE:XSQ,XSAQ ',xsq,xsaq
2605  ENDIF
2606 * target valence quarks xtvqi
2607 C here selected as sea quark! from 1/x-sea
2608  xvthro=cvq/cmener
2609  ivthr=0
2610  3466 CONTINUE
2611  IF(ivthr.EQ.200)THEN
2612  irej=1
2613  IF(isq.EQ.3)irej=3
2614  IF(ipco.GE.0)WRITE(6,*)' HADJASE 3466 reject IVTHR 50'
2615  RETURN
2616  ENDIF
2617  ivthr=ivthr+1
2618  xvthr=xvthro/(201-ivthr)
2619  unoprv=unon
2620  380 CONTINUE
2621  IF(xvthr.GT.0.66d0*xdiq)THEN
2622  irej=1
2623  IF(isq.EQ.3)irej=3
2624  IF(ipco.GE.0)WRITE(6,*)' HADJASE Rejection 380 XVTHR large ',
2625  * xvthr
2626  RETURN
2627  ENDIF
2628  xtvqi=sampey(xvthr,0.66d0*xdiq)
2629  xdiqq=xdiq-xtvqi
2630  IF(ipco.GE.0)THEN
2631  WRITE(6,'(A,2E12.4)')' HADJCK:XDIQQ,XTVQI ',xdiqq,xtvqi
2632  ENDIF
2633 C We form two strings(2-1) q-aq and (4-3) aq-aqaq
2634 C 1 aq with XDIQ-XTVQI (j.r.5.5.96)
2635 C 2 q with XSAQ
2636 C 3 aqaq with XTVQI (j.r.5.5.96)
2637 C 4 aq with XQUA-XSAQ
2638 C with 4-momenta PPP1,PPP2,PPP3,PPP4
2639  DO 3345 i=1,4
2640  IF(xdiq.LE.1.d-15.OR.xqua.LT.1.d-15)THEN
2641  irej=1
2642  IF(isq.EQ.3)irej=3
2643  IF(ipco.GE.0)WRITE (6,*)' HSDJSE Rejection 3345 XDIQ,XQUA ',
2644  * xdiq,xqua
2645  RETURN
2646  ENDIF
2647  ppp1(i)=pta(i)*(xdiq-xtvqi)/xdiq
2648  ppp3(i)=pta(i)*(xtvqi)/xdiq
2649  ppp2(i)=ppr(i)*xsaq/xqua
2650  ppp4(i)=ppr(i)*(xqua-xsaq)/xqua
2651  3345 CONTINUE
2652  3346 FORMAT(a,5e12.4)
2653  IF(ipco.GE.0)THEN
2654  WRITE(6,3346)' PPR ',ppr
2655  WRITE(6,3346)' PPP1 ',ppp1
2656  WRITE(6,3346)' PPP3 ',ppp3
2657  WRITE(6,3346)' PTA ',pta
2658  WRITE(6,3346)' PPP2 XSAQ',ppp2,xsaq
2659  WRITE(6,3346)' PPP4 ',ppp4
2660  ENDIF
2661 C Invariant Masses of chains
2662 C ORIGINAL CHAIN
2663  amchor=sqrt((ppr(4)+pta(4))**2-(ppr(3)+pta(3))**2
2664  * -(ppr(2)+pta(2))**2-(ppr(1)+pta(1))**2)
2665  IF(ipco.GE.0)WRITE(6,2346)'AMCHOR ',amchor
2666  amchn1=sqrt((ppp1(4)+ppp2(4))**2-(ppp1(3)+ppp2(3))**2
2667  * -(ppp1(2)+ppp2(2))**2-(ppp1(1)+ppp2(1))**2)
2668  IF(ipco.GE.0)WRITE(6,2346)'AMCHN1 ',amchn1
2669 C Chain 1 is q-aq restrict mass >1.5 GeV
2670  IF(amchor.LE.tiny)THEN
2671  IF(ipco.GE.0)WRITE(6,2346)' PPR ',ppr
2672  WRITE(6,2346)' PPP1 ',ppp1
2673  WRITE(6,2346)' PPP3 ',ppp3
2674  WRITE(6,2346)' PTA ',pta
2675  WRITE(6,2346)' PPP2 ',ppp2
2676  WRITE(6,2346)' PPP4 ',ppp4
2677  WRITE(6,2346)'AMCHOR ',amchor
2678  ENDIF
2679  IF(amchn1.LE.tiny)THEN
2680  WRITE(6,2346)' PPR ',ppr
2681  WRITE(6,2346)' PPP1 ',ppp1
2682  WRITE(6,2346)' PPP3 ',ppp3
2683  WRITE(6,2346)' PTA ',pta
2684  WRITE(6,2346)' PPP2 ',ppp2
2685  WRITE(6,2346)' PPP4 ',ppp4
2686  WRITE(6,2346)'AMCHOR ',amchor
2687  ENDIF
2688  chamal=0.8d0
2689  IF(ifb3.GE.9.OR.isq.GE.3)chamal=1.2d0
2690  IF(amchn1.LE.chamal)THEN
2691  IF(ipco.GE.0)WRITE(6 ,*).LE.'HADJASE jump2AMCHN1CHAMAL AMCHOR'
2692  * ,amchn1,chamal,amchor
2693  go to 3466
2694  ENDIF
2695  amchn2=sqrt((ppp3(4)+ppp4(4))**2-(ppp3(3)+ppp4(3))**2
2696  * -(ppp3(2)+ppp4(2))**2-(ppp3(1)+ppp4(1))**2)
2697 C Chain 2 is qq-q restrict mass >2.5 GeV
2698  chamal=1.4d0
2699  IF(ifb2.GE.9.OR.ifb1.GE.9.OR.isaq.GE.9)chamal=1.8d0
2700  IF(amchn2.LE.chamal)THEN
2701 C IREJ=1
2702 C RETURN
2703  IF(ipco.GE.0)WRITE(6 ,*).LE.' HADJASE jump AMCHN2CHAMAL',
2704  * amchn2,chamal
2705  go to 3466
2706  ENDIF
2707  pxchk=ppr(1)+pta(1)-ppp1(1)-ppp2(1)-ppp3(1)-ppp4(1)
2708  pychk=ppr(2)+pta(2)-ppp1(2)-ppp2(2)-ppp3(2)-ppp4(2)
2709  pzchk=ppr(3)+pta(3)-ppp1(3)-ppp2(3)-ppp3(3)-ppp4(3)
2710  pechk=ppr(4)+pta(4)-ppp1(4)-ppp2(4)-ppp3(4)-ppp4(4)
2711  IF(ipco.GE.0)THEN
2712  WRITE(6,'(A/8E12.4,I5)')
2713  * ' Chain masses AMCH,AMCHOR,AMCHN1,AMCHN2,PZCHK,PECHK,
2714  * PXCHK,PYCHK,NOBAM ',
2715  * amch,amchor,amchn1,amchn2,pzchk,pechk,pxchk,pychk,nobam
2716  ENDIF
2717 C Lorentz parameters of chains
2718  IF(amch.LE.1.d-15)THEN
2719  irej=1
2720  IF(isq.EQ.3)irej=3
2721  IF(ipco.GE.0)WRITE(6,*)' HSDJSE Rejection AMCH too small',amch
2722  RETURN
2723  ENDIF
2724  gamor=(ppr(4)+pta(4))/amch
2725  bgxor=(ppr(1)+pta(1))/amch
2726  bgyor=(ppr(2)+pta(2))/amch
2727  bgzor=(ppr(3)+pta(3))/amch
2728  IF(amchn1.LE.1.d-15)THEN
2729  IF(ipco.GE.0)WRITE(6,*)' HSDJSE Rejection AMCHN1 too small',
2730  * amchn1
2731  irej=1
2732  IF(isq.EQ.3)irej=3
2733  RETURN
2734  ENDIF
2735  gamch1=(ppp1(4)+ppp2(4))/amchn1
2736  bgxch1=(ppp1(1)+ppp2(1))/amchn1
2737  bgych1=(ppp1(2)+ppp2(2))/amchn1
2738  bgzch1=(ppp1(3)+ppp2(3))/amchn1
2739  IF(amchn2.LE.1.d-15)THEN
2740  irej=1
2741  IF(isq.EQ.3)irej=3
2742  IF(ipco.GE.0)WRITE(6,*)' HSDJSE Rejection AMCHN2 too small',
2743  * amchn2
2744  RETURN
2745  ENDIF
2746  gamch2=(ppp3(4)+ppp4(4))/amchn2
2747  bgxch2=(ppp3(1)+ppp4(1))/amchn2
2748  bgych2=(ppp3(2)+ppp4(2))/amchn2
2749  bgzch2=(ppp3(3)+ppp4(3))/amchn2
2750  IF(ipco.GE.0)THEN
2751  WRITE(6,3346)' L.Parm in ',gam,bgx,bgy,bgz
2752  WRITE(6,3346)' L.Parm OR ',gamor,bgxor,bgyor,bgzor
2753  WRITE(6,3346)' L.Parm C1 ',gamch1,bgxch1,bgych1,bgzch1
2754  WRITE(6,3346)' L.Parm C2 ',gamch2,bgxch2,bgych2,bgzch2
2755  ENDIF
2756  peor=amch*gamor
2757  pzor=amch*bgzor
2758  pech1=amchn1*gamch1
2759  pzch1=amchn1*bgzch1
2760  pech2=amchn2*gamch2
2761  pzch2=amchn2*bgzch2
2762  IF(ipco.GE.0)THEN
2763  WRITE(6,'(A,6E12.4)')' PEOR,PECH1,PECH2,PZOR,PZCH1,PZCH2',
2764  * peor,pech1,pech2,pzor,pzch1,pzch2
2765  ENDIF
2766  noba1=3
2767  noba2=4
2768  IF(ipco.GE.0)THEN
2769  WRITE(6,'(A,2I5)')' Jet1 ISQ,IFB3 ',isq,ifb3
2770  ENDIF
2771  CALL hadjet(nhad1,amchn1,ppp2,ppp1,gamch1,bgxch1,bgych1,
2772  * bgzch1, isq,ifb3,ifb1,ifb4,i1,i2,noba1,nnch,norig)
2773 C DO 342 I=1,NHAD1
2774 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
2775 C + IBARF(I),NREF(I),ANF(I)
2776 C 342 CONTINUE
2777 C DIMENSION ANFF(NFIMAX),PXFF(NFIMAX),PYFF(NFIMAX),PZFF(NFIMAX),
2778 C + HEFF(NFIMAX),AMFF(NFIMAX),
2779 C * ICHFF(NFIMAX),IBARFF(NFIMAX),NREFF(NFIMAX)
2780 C Intermediate store of hadrons from jet 1
2781  DO 3348 i=1,nhad1
2782  anff(i)=anf(i)
2783  pxff(i)=pxf(i)
2784  pyff(i)=pyf(i)
2785  pzff(i)=pzf(i)
2786  heff(i)=hef(i)
2787  amff(i)=amf(i)
2788  ichff(i)=ichf(i)
2789  ibarff(i)=ibarf(i)
2790  nreff(i)=nref(i)
2791  iormoo(i)=iormo(i)
2792  3348 CONTINUE
2793  IF(ipco.GE.0)THEN
2794  WRITE(6,'(A,3I5)')' Jet2 IFB1,ISAQ,IFB2 ',ifb1,isaq,ifb2
2795  ENDIF
2796  CALL hadjet(nhad2,amchn2,ppp4,ppp3,gamch2,bgxch2,bgych2,
2797  * bgzch2, ifb1,isaq,ifb2,ifb4,i1,i2,noba2,nnch,norig)
2798 C DO 341 I=1,NHAD2
2799 C WRITE(6,1050)I,PXF(I),PYF(I),PZF(I),HEF(I),AMF(I), ICHF(I),
2800 C + IBARF(I),NREF(I),ANF(I)
2801 C IF(IORMO(I).NE.999)IORMO(I+NHAD1)=IORMO(I)+NHAD1
2802 C 341 CONTINUE
2803 C Intermediate store of hadrons from jet 2
2804  DO 3448 i=nhad1+1,nhad1+nhad2
2805  ii=i-nhad1
2806  anff(i)=anf(ii)
2807  pxff(i)=pxf(ii)
2808  pyff(i)=pyf(ii)
2809  pzff(i)=pzf(ii)
2810  heff(i)=hef(ii)
2811  amff(i)=amf(ii)
2812  ichff(i)=ichf(ii)
2813  ibarff(i)=ibarf(ii)
2814  nreff(i)=nref(ii)
2815  iormoo(i)=iormo(ii)
2816  IF(iormoo(i).NE.999)iormoo(i)=iormoo(i)+nhad1
2817  3448 CONTINUE
2818  nhad=nhad1+nhad2
2819  DO 3349 i=1,nhad
2820  ii=i
2821  anf(i)=anff(ii)
2822  pxf(i)=pxff(ii)
2823  pyf(i)=pyff(ii)
2824  pzf(i)=pzff(ii)
2825  hef(i)=heff(ii)
2826  amf(i)=amff(ii)
2827  ichf(i)=ichff(ii)
2828  ibarf(i)=ibarff(ii)
2829  nref(i)=nreff(ii)
2830  iormo(i)=iormoo(ii)
2831  3349 CONTINUE
2832  ENDIF
2833 C==================================================================
2834  heft=0.d0
2835  IF(ipco.GE.0)THEN
2836  DO 40 i=1,nhad
2837  IF(ibarf(i).EQ.500)go to 4040
2838  heft=heft+hef(i)
2839  WRITE(6,4050)i,pxf(i),pyf(i),pzf(i),hef(i),amf(i), ichf(i),
2840  + ibarf(i),nref(i),anf(i)
2841  4050 FORMAT(' JET ',i5,5f12.4,3i5,a10)
2842  4040 CONTINUE
2843  40 CONTINUE
2844  ENDIF
2845  hefff=amch*gam
2846  IF(ipco.GE.0)THEN
2847  WRITE(6,'(A,I5,2E12.4)')'3IREJ,HEFT,HEFFF ',
2848  * irej,heft,hefff
2849  ENDIF
2850 C IREJ=1
2851  IF(irej.GE.1)RETURN
2852 C
2853 C WRITE(6,*)' HADJASE: IREJ,ISQ ',IREJ,ISQ
2854  IF(isq.EQ.1)nhase1=nhase1+1
2855  IF(isq.EQ.2)nhase2=nhase2+1
2856  IF(isq.EQ.3)nhase3=nhase3+1
2857  RETURN
2858  END
2859 
2860 ************************************************************************
2861 ************************************************************************
mpdgha
function mpdgha(MCIND)
Definition: dpm25nulib.f:386
c2_transformation::zero
static float_type zero(float_type)
utility function f(x)=0 useful in axis transforms
Definition: c2_function.hh:1203
z
Double_t z
Definition: plot.C:279
hadjse
subroutine hadjse(NHAD, AMCH, PPR, PTA, GAM, BGX, BGY, BGZ, IFB1, IFB2, IFB3, IFB4, I1, I2, NOBAM, NNCH, NORIG, IREJ, IISSQQ)
Definition: dpm25nuc4.f:1479
mcihad
function mcihad(MCIND)
Definition: dpm25nulib.f:364
d
Float_t d
Definition: plot.C:237
rndm
DOUBLE PRECISION function rndm(RDUMMY)
Definition: dpm25nulib.f:1460
pyp
function pyp(I, J)
Definition: pythia61.f:38097
pinkla
subroutine pinkla(I1, I2, I3, I4, IHKK1, IHKK2, ILOML, ILOMA, ITOML, ITOMA, IREJ)
Definition: dpm25nuc4.f:702
a
G4double a
Definition: TRTMaterials.hh:39
dpoli
subroutine dpoli(CS, SI)
Definition: dpm25nulib.f:597
nmxhkk
const int nmxhkk
Definition: G4DPMJET2_5Interface.hh:78
dbetar
DOUBLE PRECISION function dbetar(GAM, ETA)
Definition: dpm25nulib.f:289
G4StokesVector::p3
G4double p3() const
Definition: G4StokesVector.hh:76
sewew
subroutine sewew(IOP, NHKKH1)
Definition: dpm25nuc4.f:33
testem0.App.init
def init
Definition: testem0.py:56
daltra
subroutine daltra(GA, BGX, BGY, BGZ, PCX, PCY, PCZ, EC, P, PX, PY, PZ, E)
Definition: dpm25nulib.f:542
hadjck
subroutine hadjck(NHAD, AMCH, PPR, PTA, GAM, BGX, BGY, BGZ, IFB1, IFB2, IFB3, IFB4, I1, I2, NOBAM, NNCH, NORIG, IREJ)
Definition: dpm25nuc4.f:958
x
Double_t x
Definition: plot.C:279
pt
TMarker * pt
Definition: egs.C:22
sampey
DOUBLE PRECISION function sampey(X1, X2)
Definition: dpm25nulib.f:1107
dsfecf
subroutine dsfecf(SFE, CFE)
Definition: dpm25nuc7.f:3354
title
subroutine title(NA, NB, NCA, NCB)
Definition: dpm25nuc7.f:1744
c2_factory::log
static c2_log_p< float_type > & log()
make a *new object
Definition: c2_factory.hh:138
npoiss
function npoiss(AVN)
Definition: dpm25nuc4.f:10
G4Abla::tau
G4double tau(G4double bet, G4double homega, G4double ef, G4double t)
Definition: G4Abla.cc:2586
c2_factory::sqrt
static c2_sqrt_p< float_type > & sqrt()
make a *new object
Definition: c2_factory.hh:142
xseapa
subroutine xseapa(ECM, XXXX, IPSQ1, IPSAQ1, XPSQ1, XPSAQ1, IREJ)
Definition: dpm25nuc7.f:1323
betrej
DOUBLE PRECISION function betrej(GAM, ETA, XMIN, XMAX)
Definition: dpm25nulib.f:344
hadjase
subroutine hadjase(NHAD, AMCH, PPR, PTA, GAM, BGX, BGY, BGZ, IFB1, IFB2, IFB3, IFB4, I1, I2, NOBAM, NNCH, NORIG, IREJ, IISSQQ)
Definition: dpm25nuc4.f:2176
hadjet
subroutine hadjet(NHAD, AMCH, PPR, PTA, GAM, BGX, BGY, BGZ, IFB1, IFB2, IFB3, IFB4, I1, I2, NOBAM, NNCH, NORIG)
Definition: dpm25nuc3.f:5651
c2_factory::exp
static c2_exp_p< float_type > & exp()
make a *new object
Definition: c2_factory.hh:140
Generated on Sat Dec 14 2013 14:33:20 for Geant4_10 by   doxygen 1.8.5