Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
dtu25lap.f
Go to the documentation of this file.
1 C------- name of the file ----------------------------------------------
2 C DTULAP.FOR
3 C______________________________________________________________________
4 C
5 C originally: JTDTU.FOR program
6 C connection between DTU and JT ( hard scattering )
7 *
8 C first parameters are taken from DTU and set to the own
9 C parameter-commons of JT, then JT will be initialized
10 C
11 C______________________________________________________________________
12 * revision 3.92: adjust COMMONS,
13 * caraful: DTU90 was based on a older version of DTULAP
14 * ********************************************************************
15  SUBROUTINE jtdtu(IOPT)
16 *
17  IMPLICIT DOUBLE PRECISION(a-h,o-z)
18  SAVE
19  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
20  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
21  COMMON /haenvi/ nindep
22  COMMON /haoutl/ noutl,nouter,noutco
23  COMMON /hapadi/ npdm
24  COMMON /haqqap/ aqqal,aqqpd,nqqal,nqqpd
25  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
26 *
27  CHARACTER*80 title
28  CHARACTER*8 projty,targty
29 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
30 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
31  COMMON /user1/title,projty,targty
32  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
33 *
34  common/collis/s,ijproj,ijtar,ptthr,ptthr2,iophrd,ijprlu,ijtalu
35 C repl. COMMON /COLLIS/ECMDTU,S,IJPROJ,IJTAR,PTTHR,IOPHRD,IJPRLU,IJTALU,PTTHR2
36  COMMON /strufu/istrum,istrut
37  COMMON /ptlarg/xsmax
38  COMMON /haxsum/xshmx
39 C
40 C===read default values
41 C
42  CALL hastrt
43 C
44 C===read parameters from DTU
45 C
46 C max. # of flavors
47  nf = 4
48 C partondistributions
49  npd = istruf
50  npdm = istrum
51 C correct scale for CTEQ PDFs
52  IF((istruf.GE.16).OR.(istruf.LE.20)) THEN
53  aqqal = 1.d0
54  aqqpd = 1.d0
55  ENDIF
56 C hadron a
57  nha = ijproj
58  IF ( ijproj.EQ.2 ) nha =-1
59 C hadron b
60  nhb = ijtar
61  IF ( ijtar .EQ.2 ) nhb =-1
62 C output level
63  noutl = ioutpa
64 C cms-energy ( GeV )
65 C repl ECM=ECMDTU
66  ecm = cmener
67 C pt-cut ( GeV )
68  ptini(1) = ptthr
69  ptini(2) = ptthr2
70  ptini(3) = 0.0
71  ptini(4) = 0.0
72 C maximum sum of hard x
73  xshmx = xsmax
74 C is program called from DTU ( NINDEP=0 ) or independent ( NINDEP=1 )
75  nindep = iopt
76 C
77 C===initialize JT
78 C
79  CALL hisini
80  IF ( iopt.EQ.0 ) CALL harini
81  RETURN
82  END
83 C******************************************************************************
84  SUBROUTINE selhrd(MHARD,IJPVAL,IJTVAL,PTTHRE)
85 C
86 C select the initial parton x-fractions and flavors and the final flavors
87 C for an event with mhard hard or semihard scatterings
88 C
89 C IJPVAL,IJTVAL =0 valence quarks of projectile or target not involved
90 C in hard scattering
91 C IJPVAL,IJTVAL =1 valence quarks of projectile or target involved
92 C in hard scattering
93 C
94 C the results are in COMMON /ABRHRD/
95 C XH1(I),XH2(I): x-values of initial partons
96 C IJHI1(I),IJHI2(I): flavor of initial parton
97 C 0 gluon
98 C 1,2 valence u,d quarks
99 C 11,12,13,14 sea udsc-quarks
100 C negative anti s or v quarks
101 C IJHF1(I),IJHF2(I): flavor of final state partons
102 C PHARD1(I,J),PHARD2(I,J): final part. momentum and energy
103 C J=1 PX
104 C =2 PY
105 C =3 PZ
106 C =4 ENERGY (massless partons)
107 C-----------------------------------------------------------------------
108  IMPLICIT DOUBLE PRECISION(a-h,o-z)
109  SAVE
110  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
111  CHARACTER*80 title
112  CHARACTER*8 projty,targty
113 C COMMON /USER/TITLE,PROJTY,TARGTY,CMENER,ISTRUF
114 C & ,ISINGD,IDUBLD,SDFRAC,PTLAR
115  COMMON /user1/title,projty,targty
116  COMMON /user2/cmener,sdfrac,ptlar,istruf,isingd,idubld
117  common/collis/s,ijproj,ijtar,ptthr,ptthr2,iophrd,ijprlu,ijtalu
118  COMMON /abrhrd/xh1(mscahd),xh2(mscahd),ijhi1(mscahd),
119  *ijhi2(mscahd),ijhf1(mscahd),ijhf2(mscahd),phard1(mscahd,4),
120  *phard2(mscahd,4)
121  COMMON /outlev/ioutpo,ioutpa,iouxev,ioucol
122  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
123  COMMON /haoutl/ noutl,nouter,noutco
124  COMMON /haevtr/ line,lin,lrec1(mline),lrec2(mline),prec(0:3,mline)
125  COMMON /harslt/ lscahd,lsc1hd,
126  & etahd(mscahd,2) ,pthd(mscahd),
127  & xhd(mscahd,2) ,vhd(mscahd) ,x0hd(mscahd,2),
128  & ninhd(mscahd,2) ,nouthd(mscahd,2),
129  & n0inhd(mscahd,2),nbrahd(mscahd,2),nprohd(mscahd)
130  DATA x1su/0./ , x2su/0./
131 C
132  ijpval =0
133  ijtval =0
134  IF (ioutpa.GE.3) WRITE(6,221)
135  * mhard,ijpval,ijtval
136  221 FORMAT (' SELHRD ',3i10)
137 C call of hard scattering routines
138  CALL harevt(mhard,ptthr2)
139 C select information from event-record
140 C number of hard scatterings reached
141  mhard = lscahd
142 C initial partons
143  DO 10 n=1,lscahd
144 C X-values
145  xh1(n) = xhd(n,1)
146  xh2(n) = xhd(n,2)
147  x1su = x1su + xh1(n)
148  x2su = x2su + xh2(n)
149  IF( ioutpa.GT. 6 )WRITE(6,*)n,x1su,x2su,xh1(n),xh2(n)
150 C flavors
151  iii = ninhd(n,1)
152  iiia = abs(iii)
153  IF ( iiia.GT. 0 .AND. iiia.LT.10 ) iii = sign(iiia+10,iii)
154  IF ( iiia.GE.10 ) iii = sign(iiia-10,iii)
155  IF ( iiia.GE.10 ) ijpval = 1
156  ijhi1(n) = iii
157  iii = ninhd(n,2)
158  iiia = abs(iii)
159  IF ( iiia.GT. 0 .AND. iiia.LT.10 ) iii = sign(iiia+10,iii)
160  IF ( iiia.GE.10 ) iii = sign(iiia-10,iii)
161  IF ( iiia.GE.10 ) ijtval = 1
162  ijhi2(n) = iii
163 10 CONTINUE
164 C final partons
165  DO 30 n=1,lscahd
166  i3 = 4*n-1
167  i4 = 4*n
168 C flavors
169  ijhf1(n) = nouthd(n,1)
170  ijhf2(n) = nouthd(n,2)
171 C four momentum
172  DO 20 j=1,3
173  phard1(n,j) = prec(j,i3)
174 20 phard2(n,j) = prec(j,i4)
175  phard1(n,4) = prec(0,i3)
176  phard2(n,4) = prec(0,i4)
177 30 CONTINUE
178 C
179 C output ( optional )
180 C
181  IF (ioutpa.GE.3)WRITE (6,101)
182  101 FORMAT(' SELHRD OUTPUT FOR INITIAL STATE SCATTERED PARTONS')
183  DO 102 i=1,lscahd
184  IF (ioutpa.GE.3)
185  * WRITE (6,103)i,ijpval,ijtval,ijhi1(i),ijhi2(i),xh1(i),xh2(i)
186  103 FORMAT (' I,IJPVAL,IJTVAL,IJHI1,IJHI2,XH1,XH2= ',5i5,2f12.6)
187  102 CONTINUE
188  IF (ioutpa.GE.3)WRITE (6,301)
189  301 FORMAT(' SELHRD OUTPUT FOR FINAL STATE SCATTERED PARTONS')
190  DO 302 i=1,lscahd
191  IF (ioutpa.GE.3)
192  * WRITE (6,303)i,ijhf1(i),ijhf2(i),(phard1(i,iii),iii=1,4)
193  IF (ioutpa.GE.3)
194  * WRITE (6,303)i,ijhf1(i),ijhf2(i),(phard2(i,iii),iii=1,4)
195  303 FORMAT (' I,IJHI1,IJHI2,PHARD1 OR PHARD2 ',3i5,4f16.6)
196  302 CONTINUE
197  RETURN
198  END
199 *
200 C______________________________________________________________________
201 C
202 C originally: JTWORK.FOR
203 C procedures to simulate a single event
204 C
205 C
206 C ( this procedures work independent to procedures in other blocks
207 C if initialization was done )
208 C
209 C______________________________________________________________________
210 *
211 * ********************************************************************
212  SUBROUTINE harevt(MHARD,PT1IN)
213 *
214  IMPLICIT DOUBLE PRECISION(a-h,o-z)
215  SAVE
216  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
217  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
218  COMMON /haenvi/ nindep
219  COMMON /haevnt/ pt1,pt2,nhard,ntry,ihard,itry,irejev
220 
221  pt1 = max(pt1in,ptini(1))
222  pt2 = ptini(1)
223  nhard = mhard
224  ihard = 0
225 C NTRY = 5
226  ntry = 200
227 C NTRY = 1000
228  itry = 0
229  irejev = 0
230  CALL hamult
231  CALL haoutp
232  IF ( nindep.EQ.1 ) CALL hisfil
233  RETURN
234  END
235 C_______________________________________________________________________
236 C===========================================================================
237 C
238 C THE FOLLOWING 5 SUBROUTINES ARE REWRITTEN BY
239 C BY I.KAWRAKOW IN ORDER TO BE ABLE
240 C TO PRODUCE GREAT NUMBER OF HARD POMERONS (>100) IN A
241 C SHORT TIME
242 C VERSION IK.1 - 01.93
243 C--------------------------------------------------------------------
244  SUBROUTINE hamult
245 C--------------------------------------------------------------------
246  IMPLICIT DOUBLE PRECISION(a-h,o-z)
247  SAVE
248  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
249  parameter( tiny= 1.d-30, one=1.d0, zsmall=1.d-3 )
250  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
251  COMMON /hapdco/ npdcor
252  COMMON /haoutl/ noutl,nouter,noutco
253  COMMON /haevnt/ pt1,pt2,nhard,ntry,ihard,itry,irejev
254  COMMON /hasca / ptwant,a,aln,z1max,z1dif,z2max,z2dif,
255  & pt,etac,etad,x1,x2,v,u,w,w1,axx,weight,mspr,irejsc
256  COMMON /haxik / xrest,yrest,zmax,axxmax,wemax
257  COMMON /haevtr/ line,lin,lrec1(mline),lrec2(mline),prec(0:3,mline)
258  COMMON /haxsum/xshmx
259  INTEGER mxsect
260  COMMON /haxsec/ xsecta(2,-1:maxpro,4),xsect(5,-1:maxpro),
261  & mxsect(0:2,-1:maxpro)
262  COMMON /harslt/ lscahd,lsc1hd,
263  & etahd(mscahd,2) ,pthd(mscahd),
264  & xhd(mscahd,2) ,vhd(mscahd) ,x0hd(mscahd,2),
265  & ninhd(mscahd,2) ,nouthd(mscahd,2),
266  & n0inhd(mscahd,2),nbrahd(mscahd,2),nprohd(mscahd)
267  itype(l) = mod(lrec1(l),100)-50
268 
269  line = 0
270  lscahd = 0
271  aa = (2.*pt2/ecm)**2
272  sa = sqrt(aa)
273 C
274 C loop until event is accepted or too many attempts ( more then NTRY )
275 C
276 5 itry = 0
277 20 itry = itry+1
278  IF(itry.GT.ntry) goto 301
279  line = 0
280  xrest = xshmx-nhard*sa
281  yrest = xshmx-nhard*sa
282  IF(xrest*yrest.LT.aa) THEN
283  WRITE(6,*) ' ****************** HAMULT ****************** '
284  WRITE(6,*) ' IT IS NOT POSSIBLE TO PRODUCE ',nhard,' POMERONS '
285  nhard=0
286  RETURN
287 C STOP
288  ENDIF
289  zmax=xrest*yrest
290  axxmax=aa/zmax
291  wemax =sqrt(1-axxmax)
292  x1s = 0.0
293  x2s = 0.0
294  ihard = 0
295  ptwant = pt1
296 10 CONTINUE
297  a = (2.*ptwant/ecm)**2
298  sa = sqrt(a)
299  i = 5
300 50 i = i-1
301  IF ( pt1.LT.ptini(i) .AND. i.GT.1 ) goto 50
302  DO 60 m=-1,maxpro
303  xsect(1,m) = xsecta(1,m,i)
304  xsect(2,m) = xsecta(2,m,i)
305 60 CONTINUE
306  CALL harsca
307  x1s = x1s+x1
308  x2s = x2s+x2
309  xrest=xrest-x1+sa
310  yrest=yrest-x2+sa
311  zmax=xrest*yrest
312  ihard=ihard+1
313  lscahd = ihard
314  xhd(ihard,1) = x1
315  xhd(ihard,2) = x2
316  vhd(ihard) = v
317  etahd(ihard,1) = etac
318  etahd(ihard,2) = etad
319  pthd(ihard) = pt
320  nprohd(ihard) = mspr
321  if(zmax/a-one.lt.zsmall) THEN
322  CALL xcheck(x1s,x2s,linmax)
323  goto 10
324  ENDIF
325  axxmax=a/zmax
326  wemax=sqrt(1.-axxmax)
327  ptwant = pt2
328  IF(ihard.LT.nhard) goto 10
329 C-------------------------------------------------- NOW THE REQUIRED NUMBER
330 C OF POMERTONS IS CREATED
331  IF ( npdcor.EQ.1 .AND.
332  & ihard .GT.1 .AND.
333  & (1.-x1s)*(1.-x2s).LT.rndm(ai)*(1.-aa*ihard)**2 ) goto 5
334  301 CONTINUE
335 C
336 C end of loop
337 C
338 C check choice of valence quarks
339  DO 120 k=1,2
340  ival = 0
341  DO 110 i=1,ihard
342  ind = 4*(i-1)+k
343  it = itype(ind)
344  IF ( abs(it).GT.10 .AND. ival.EQ.0 ) THEN
345  ival = 1
346  ELSEIF ( abs(it).GT.10 .AND. ival.EQ.1 ) THEN
347  it = sign(abs(it)-10,it)
348  lrec1(ind) = (lrec1(ind)/100)*100+50+it
349  ENDIF
350 C fill COMMON HARSLT
351  ninhd(i,k) = it
352  nouthd(i,k) = itype(ind+2)
353 110 CONTINUE
354 120 CONTINUE
355 C
356 C information if HAMULT is not able to produce the required # of scatt.
357 C
358  IF ( ihard.NE.nhard .AND. nouter.EQ.1 ) THEN
359  WRITE(6,1010) nhard,ihard
360 1010 FORMAT(' ###### HAMULT : CANNOT PRODUCE',i3,' HARD SCATT.',
361  & '; ONLY',i3,' ARE PRODUCED !!!')
362  ENDIF
363  RETURN
364  END
365 
366 C______________________________________________________________________
367  SUBROUTINE recchk( LINMAX,X,IOPT )
368  IMPLICIT DOUBLE PRECISION(a-h,o-z)
369  SAVE
370  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
371  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
372  COMMON /hapdco/ npdcor
373  COMMON /haoutl/ noutl,nouter,noutco
374  COMMON /haevnt/ pt1,pt2,nhard,ntry,ihard,itry,irejev
375  COMMON /hasca / ptwant,a,aln,z1max,z1dif,z2max,z2dif,
376  & pt,etac,etad,x1,x2,v,u,w,w1,axx,weight,mspr,irejsc
377  COMMON /haevtr/ line,lin,lrec1(mline),lrec2(mline),prec(0:3,mline)
378  INTEGER mxsect
379  COMMON /haxsec/ xsecta(2,-1:maxpro,4),xsect(5,-1:maxpro),
380  & mxsect(0:2,-1:maxpro)
381  COMMON /harslt/ lscahd,lsc1hd,
382  & etahd(mscahd,2) ,pthd(mscahd),
383  & xhd(mscahd,2) ,vhd(mscahd) ,x0hd(mscahd,2),
384  & ninhd(mscahd,2) ,nouthd(mscahd,2),
385  & n0inhd(mscahd,2),nbrahd(mscahd,2),nprohd(mscahd)
386 C
387  IF( iopt.EQ.0 ) THEN
388  lstart = linmax + 1
389  DO 1 l = lstart,line
390  lp = l - 4
391  prec(1,lp) = prec(1,l)
392  prec(2,lp) = prec(2,l)
393  prec(3,lp) = prec(3,l)
394  prec(0,lp) = prec(0,l)
395  lrec1( lp) = lrec1( l)
396  lrec2( lp) = lrec2( l)
397  1 CONTINUE
398  line = line - 4
399  RETURN
400  ELSEIF( iopt.EQ.1 ) THEN
401  qtest = 0.5*ecm*x
402  DO 2 l=1,line
403  ptest = prec(0,l)
404  IF( ptest.EQ.qtest ) THEN
405  linmax = l
406  RETURN
407  ENDIF
408  2 CONTINUE
409  WRITE(6,*)' RECCHK: NO NEW LINMAX FOUND - LINMAX=',linmax
410  RETURN
411  ENDIF
412  WRITE(6,*)' RECCHK: IOPT OUT OF RANGE - 0 OR 1 - IOPT=',iopt
413  RETURN
414  END
415 C______________________________________________________________________
416  SUBROUTINE xcheck( X1S, X2S, LINMAX )
417  IMPLICIT DOUBLE PRECISION(a-h,o-z)
418  SAVE
419  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
420  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
421  COMMON /hapdco/ npdcor
422  COMMON /haoutl/ noutl,nouter,noutco
423  COMMON /haevnt/ pt1,pt2,nhard,ntry,ihard,itry,irejev
424  COMMON /hasca / ptwant,a,aln,z1max,z1dif,z2max,z2dif,
425  & pt,etac,etad,x1,x2,v,u,w,w1,axx,weight,mspr,irejsc
426  COMMON /haxik / xrest,yrest,zmax,axxmax,wemax
427  COMMON /haevtr/ line,lin,lrec1(mline),lrec2(mline),prec(0:3,mline)
428  COMMON /haxsum/xshmx
429  INTEGER mxsect
430  COMMON /haxsec/ xsecta(2,-1:maxpro,4),xsect(5,-1:maxpro),
431  & mxsect(0:2,-1:maxpro)
432  COMMON /harslt/ lscahd,lsc1hd,
433  & etahd(mscahd,2) ,pthd(mscahd),
434  & xhd(mscahd,2) ,vhd(mscahd) ,x0hd(mscahd,2),
435  & ninhd(mscahd,2) ,nouthd(mscahd,2),
436  & n0inhd(mscahd,2),nbrahd(mscahd,2),nprohd(mscahd)
437  parameter(one=1d0, zsmall=1d-3)
438 C
439 50 CONTINUE
440  IF(ihard.LT.1) THEN
441  WRITE(6,*) ' ERROR IN XCHECK : IHARD < 1 ',ihard
442  stop
443  ENDIF
444 C-------------------------------------- FIND PROCESS WITH THE MAX. X
445  imax=0
446  xmax=0.
447  DO 10 i=1,ihard
448  IF(xhd(i,1).GT.xmax) THEN
449  imax=i
450  xmax=xhd(i,1)
451  ENDIF
452  IF(xhd(i,2).GT.xmax) THEN
453  imax=i
454  xmax=xhd(i,2)
455  ENDIF
456 10 CONTINUE
457 C--------------------------------------- REJECT THIS PROCESS
458  x1s=x1s-xhd(imax,1)
459  x2s=x2s-xhd(imax,2)
460  xrest=xrest+xhd(imax,1)-sqrt(a)
461  yrest=yrest+xhd(imax,2)-sqrt(a)
462  zmax=xrest*yrest
463  axxmax=a/zmax
464  wemax=sqrt(1.-axxmax)
465  mh=0
466  DO 20 i=1,ihard
467  IF(i.NE.imax) THEN
468  mh=mh+1
469  xhd(mh,1) = xhd(i,1)
470  xhd(mh,2) = xhd(i,2)
471  vhd(mh) = vhd(i)
472  etahd(mh,1) = etahd(i,1)
473  etahd(mh,2) = etahd(i,2)
474  pthd(mh) = pthd(i)
475  nprohd(mh) = nprohd(i)
476  ENDIF
477 20 CONTINUE
478  CALL recchk( 4*imax,xhd1,0)
479  ihard=ihard-1
480  lscahd=ihard
481  IF(zmax/a-one.LT.zsmall) goto 50
482  RETURN
483  END
484 C_______________________________________________________
485  SUBROUTINE hax1x2
486  IMPLICIT DOUBLE PRECISION(a-h,o-z)
487  SAVE
488  COMMON /haevnt/ pt1,pt2,nhard,ntry,ihard,itry,irejev
489  COMMON /hasca / ptwant,a,aln,z1max,z1dif,z2max,z2dif,
490  & pt,etac,etad,x1,x2,v,u,w,w1,axx,weight,mspr,irejsc
491 C COMMON /HAXIK / XREST,YREST,ZMAX
492  COMMON /haxik / xrest,yrest,zmax,axxmax,wemax
493  parameter( tiny= 1.d-30, one=1.d0 ,tiny6=1.d-06)
494 
495  sa=sqrt(a)
496 12 continue
497 c--------------------------------- sample z=x*y
498  z=a*exp(rndm(1.)*log(zmax/a))
499  xm=xrest
500  ym=yrest
501  if(xm.lt.yrest) then
502  xm=yrest
503  ym=xrest
504  endif
505  ww=log(xm**2/z)/log(xm**2/a)
506  if(rndm(1.1).gt.ww) goto 12
507 c--------------------------------- sample u=x+y
508  umin=sqrt(4.*z)
509  umax=xm+z/xm
510  cc=umax**2-4.*z
511  if(cc.lt.0.) cc=0.
512 13 continue
513  c=exp(rndm(2.)*log((umax+sqrt(cc))/umin))
514  uu=umin*(c**2+1.)/2./c
515  if(uu.gt.2.*ym.and.uu.lt.ym+z/ym) goto 13
516 c------------------------------------- x,y from u,z
517  c=uu**2-4.*z
518  if(c.lt.0.) c=0.
519  c=sqrt(c)
520  xtemp=(uu+c)/2.
521  ytemp=(uu-c)/2.
522  if(xrest.ge.yrest) then
523  x=xtemp
524  y=ytemp
525  if(xrest.eq.yrest) then
526  if(rndm(3.).gt.0.5) then
527  x=ytemp
528  y=xtemp
529  endif
530  endif
531  else
532  x=ytemp
533  y=xtemp
534  endif
535  x1=x
536  x2=y
537  axx = a/(x1*x2)
538  w = sqrt(max(tiny,one-axx))
539  w1 = axx/(1.+w)
540  RETURN
541  END
542 C______________________________________________________________________
543  SUBROUTINE harkin
544  IMPLICIT DOUBLE PRECISION(a-h,o-z)
545  SAVE
546  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
547  parameter( tiny= 1.d-30, one=1.d0 ,tiny6=1.d-06)
548  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
549  COMMON /haevnt/ pt1,pt2,nhard,ntry,ihard,itry,irejev
550  COMMON /hasca / ptwant,a,aln,z1max,z1dif,z2max,z2dif,
551  & pt,etac,etad,x1,x2,v,u,w,w1,axx,weight,mspr,irejsc
552  dimension rm(-1:maxpro)
553  COMMON /haxik / xrest,yrest,zmax,axxmax,wemax
554  DATA rm / 3.31, 0.0,
555  & 3.80, 0.65, 2.00, 0.65, 0.89, 0.45, 0.445, 0.89 /
556  m = mspr
557  IF ( m.EQ.1 ) THEN
558 10 CALL hax1x2
559  v =-0.5*w1/(w1+rndm(ai)*w)
560  u =-1.-v
561  r = (1.+w)*2.25*(v*v*(3.-u*v-v/(u*u))-u)
562  rmax=rm(1)*wemax*(1.+wemax)
563  wik=r*w/rmax
564  IF(wik.GT.1.d0) WRITE(6,*) ' HARKIN : WIK > 1 : ',m,r
565 C IF ( R*W.LT.RM(1)*RNDM(AI) ) GOTO 10
566  IF(wik.LT.rndm(ai)) goto 10
567  IF ( rndm(aj).LE.0.5d0 ) v = u
568  ELSEIF ( m.EQ.2 .OR. m.EQ.4 ) THEN
569 20 CALL hax1x2
570  wl = log(w1)
571  v =-exp(-0.6931472+rndm(ai)*wl)
572  u =-1.-v
573  r = (u*u+v*v)*((16./27.)/u-(4./3.)*v)*(wl/w)*axx
574  IF ( r*w.LT.rm(m)*rndm(ai) ) goto 20
575  IF ( rndm(aj).LE.0.5d0 ) v = u
576  ELSEIF ( m.EQ.3 ) THEN
577 30 CALL hax1x2
578  v =-0.5*w1/(w1+rndm(ai)*w)
579  u =-1.-v
580  r = (1.+w)*(1.+u*u)*(1.-(4./9.)*v*v/u)
581  rmax=rm(3)*wemax*(1.+wemax)
582  wik=r*w/rmax
583  IF(wik.GT.1.d0) WRITE(6,*) ' HARKIN : WIK > 1 : ',m,r
584 C IF ( R*W.LT.RM(3)*RNDM(AI) ) GOTO 30
585  IF(wik.LT.rndm(ai)) goto 30
586  ELSEIF ( m.EQ.5 ) THEN
587 50 CALL hax1x2
588  v =-0.5*axx/(w1+2.*rndm(ai)*w)
589  u =-1.-v
590  r = (4./9.)*(1.+u*u+v*v*(u*u+v*v))-(8./27.)*u*u*v
591  rmax=rm(5)*wemax
592  wik=r*w/rmax
593  IF(wik.GT.1.d0) WRITE(6,*) ' HARKIN : WIK > 1 : ',m,r
594 C IF ( R*W.LT.RM(5)*RNDM(AI) ) GOTO 50
595  IF(wik.LT.rndm(ai)) goto 50
596  ELSEIF ( m.EQ.6 ) THEN
597 60 CALL hax1x2
598  v =-0.5*(1.+w)+rndm(ai)*w
599  u =-1.-v
600  r = (4./9.)*(u*u+v*v)*axx
601  IF ( r*w.LT.rm(6)*rndm(ai) ) goto 60
602  ELSEIF ( m.EQ.7 ) THEN
603 70 CALL hax1x2
604  v =-0.5*w1/(w1+rndm(ai)*w)
605  u =-1.-v
606  r = (1.+w)*((2./9.)*(1.+u*u+(1.+v*v)*v*v/(u*u))-(4./27.)*v/u)
607  rmax=rm(7)*wemax*(1.+wemax)
608  wik=r*w/rmax
609  IF(wik.GT.1.d0) WRITE(6,*) ' HARKIN : WIK > 1 : ',m,r
610 C IF ( R*W.LT.RM(7)*RNDM(AI) ) GOTO 70
611  IF(wik.LT.rndm(ai)) goto 70
612  IF ( rndm(aj).LE.0.5d0 ) v = u
613  ELSEIF ( m.EQ.8 ) THEN
614 80 CALL hax1x2
615  v =-0.5*axx/(w1+2.*rndm(ai)*w)
616  u =-1.-v
617  r = (4./9.)*(1.+u*u)
618  rmax=rm(8)*wemax
619  wik=r*w/rmax
620  IF(wik.GT.1.d0) WRITE(6,*) ' HARKIN : WIK > 1 : ',m,r
621 C IF ( R*W.LT.RM(8)*RNDM(AI) ) GOTO 80
622  IF(wik.LT.rndm(ai)) goto 80
623  ELSEIF ( m.EQ.-1 ) THEN
624 90 CALL hax1x2
625  wl = log(w1)
626  v =-exp(-0.6931472+rndm(ai)*wl)
627  u =-1.-v
628  r = (1.+v*v)*(v/(u*u)-(4./9.))*(wl/w)*axx
629  IF ( r*w.LT.rm(-1)*rndm(ai) ) goto 90
630  ENDIF
631 C PARAMETER ( TINY= 1.D-30, ONE=1.D0 ,TINY6 =1.D-06)
632  v = max(min( v,-tiny6 ),-1.+tiny6 )
633  u = max(min(-1.e0-v,-tiny6 ),-1.+tiny6 )
634  pt = sqrt(u*v*x1*x2)*ecm
635  etac = 0.5*log((u*x1)/(v*x2))
636  etad = 0.5*log((v*x1)/(u*x2))
637  RETURN
638  END
639 C-------------------------------------------- END OF CHANGES BY IK 01.93
640 C===========================================================================
641 C_______________________________________________________________________
642 
643  SUBROUTINE hachek(IOPT)
644  IMPLICIT DOUBLE PRECISION(a-h,o-z)
645  SAVE
646  COMMON /hacuts/ ptl,ptu,etacl,etacu,etadl,etadu
647  COMMON /hasca / ptwant,a,aln,z1max,z1dif,z2max,z2dif,
648  & pt,etac,etad,x1,x2,v,u,w,w1,axx,weight,mspr,irejsc
649 
650  iopt = 1
651  IF ( pt .LT.ptl .OR. pt .GT.ptu
652  & .OR. etac.LT.etacl .OR. etac.GT.etacu
653  & .OR. etad.LT.etadl .OR. etad.GT.etadu ) iopt = 0
654  RETURN
655  END
656 C______________________________________________________________________
657  SUBROUTINE hafdis(PDS,PDA,PDB,FDISTR)
658  IMPLICIT DOUBLE PRECISION(a-h,o-z)
659  SAVE
660  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
661  parameter( tiny= 1.d-30, one=1.d0 ,tiny6=1.d-06)
662  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
663  COMMON /haqqap/ aqqal,aqqpd,nqqal,nqqpd
664  COMMON /haevnt/ pt1,pt2,nhard,ntry,ihard,itry,irejev
665  COMMON /hasca / ptwant,a,aln,z1max,z1dif,z2max,z2dif,
666  & pt,etac,etad,x1,x2,v,u,w,w1,axx,weight,mspr,irejsc
667  dimension pda(-6:6),pdb(-6:6)
668  INTEGER mxsect
669  COMMON /haxsec/ xsecta(2,-1:maxpro,4),xsect(5,-1:maxpro),
670  & mxsect(0:2,-1:maxpro)
671 
672  fdistr = 0.0
673 C set hard scale QQ for alpha and partondistr.
674  IF ( nqqal.EQ.1 ) THEN
675  qqal = aqqal*pt*pt
676  ELSEIF ( nqqal.EQ.2 ) THEN
677  qqal = aqqal*x1*x2*ecm*ecm
678  ELSEIF ( nqqal.EQ.3 ) THEN
679  qqal = aqqal*x1*x2*ecm*ecm*(u*v)**(1./3.)
680  ELSEIF ( nqqal.EQ.4 ) THEN
681  qqal = aqqal*x1*x2*ecm*ecm*u*v/(1.+v*v+u*u)
682  ENDIF
683  IF ( nqqpd.EQ.1 ) THEN
684  qqpd = aqqpd*pt*pt
685  ELSEIF ( nqqpd.EQ.2 ) THEN
686  qqpd = aqqpd*x1*x2*ecm*ecm
687  ELSEIF ( nqqpd.EQ.3 ) THEN
688  qqpd = aqqpd*x1*x2*ecm*ecm*(u*v)**(1./3.)
689  ELSEIF ( nqqpd.EQ.4 ) THEN
690  qqpd = aqqpd*x1*x2*ecm*ecm*u*v/(1.+v*v+u*u)
691  ENDIF
692  alpha = bqcd/log(max(qqal/alasqr,1.1*one))
693  f = xsect(1,mspr)*alpha**2
694 C calculate partondistributions
695  CALL jtpdis(x1,qqpd,nha,mspr,pda)
696  CALL jtpdis(x2,qqpd,nhb,mspr,pdb)
697 C calculate full distribution FDISTR
698  IF ( mspr.EQ.1 .OR. mspr.EQ.4 ) THEN
699  pds = pda(0)*pdb(0)
700  ELSE
701  s2 = 0.0
702  s3 = 0.0
703  s4 = 0.0
704  s5 = 0.0
705  DO 10 i=1,nf
706  s2 = s2+pda(i)*pdb(-i)+pda(-i)*pdb( i)
707  s3 = s3+pda(i)*pdb( i)+pda(-i)*pdb(-i)
708  s4 = s4+pda(i)+pda(-i)
709  s5 = s5+pdb(i)+pdb(-i)
710 10 CONTINUE
711  IF ( mspr.EQ.2 .OR. mspr.EQ.5 .OR. mspr.EQ.6 ) THEN
712  pds = s2
713  ELSEIF ( mspr.EQ.3 .OR. mspr.EQ.-1 ) THEN
714  pds = pda(0)*s5+pdb(0)*s4
715  ELSEIF ( mspr.EQ.7 ) THEN
716  pds = s3
717  ELSEIF ( mspr.EQ.8 ) THEN
718  pds = s4*s5-(s2+s3)
719  ENDIF
720  ENDIF
721  fdistr = f*pds
722  RETURN
723  END
724 C______________________________________________________________________
725  SUBROUTINE harsca
726 C HARSCA determines the type of hard subprocess, the partons taking
727 C part in subprocess and the kinematic variables
728  IMPLICIT DOUBLE PRECISION(a-h,o-z)
729  SAVE
730  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
731  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
732  COMMON /haoutl/ noutl,nouter,noutco
733  COMMON /haevnt/ pt1,pt2,nhard,ntry,ihard,itry,irejev
734  COMMON /hasca / ptwant,a,aln,z1max,z1dif,z2max,z2dif,
735  & pt,etac,etad,x1,x2,v,u,w,w1,axx,weight,mspr,irejsc
736  dimension pda(-6:6),pdb(-6:6)
737  COMMON /haevtr/ line,lin,lrec1(mline),lrec2(mline),prec(0:3,mline)
738  INTEGER mxsect
739  COMMON /haxsec/ xsecta(2,-1:maxpro,4),xsect(5,-1:maxpro),
740  & mxsect(0:2,-1:maxpro)
741 
742  mxsect(0,0) = 0
743  xsect(2,0) = 0.0
744  DO 15 m=-1,maxpro
745  IF ( mxsect(0,m).EQ.1 ) xsect(2,0) = xsect(2,0)+xsect(2,m)
746 15 CONTINUE
747 C
748 C -------------------------------------------I
749 C begin of iteration loop I
750 C I
751  irejsc = 0
752 10 CONTINUE
753  irejsc = irejsc+1
754  irejev = irejev+1
755 C find subprocess
756  b = rndm(ai)*xsect(2,0)
757  mspr =-2
758  sum = 0.0
759 20 mspr = mspr+1
760  IF ( mxsect(0,mspr).EQ.1 ) sum = sum+xsect(2,mspr)
761  IF ( sum.LT.b .AND. mspr.LT.maxpro ) goto 20
762 C find kin. variables X1,X2 and V
763  CALL harkin
764 C check kin. cuts eventually given by user
765  CALL hachek(iopt)
766  IF ( iopt.EQ.0 ) goto 10
767 C calculate remaining distribution
768  CALL hafdis(pds,pda,pdb,f)
769 C actualize counter for cross-section calculation
770  IF( f .LE. 1.d-15 ) f=0.
771  xsect(3,mspr) = xsect(3,mspr)+f
772  xsect(4,mspr) = xsect(4,mspr)+f*f
773  mxsect(1,mspr) = mxsect(1,mspr)+1
774 C check F against FMAX
775 C
776  weight = f/xsect(2,mspr)
777  IF ( weight.LT.rndm(ai) ) goto 10
778 C-------------------------------------------------------------------
779 C IF(WEIGHT.GT.1.D0) WRITE(6,1234)F,XSECT(2,MSPR),WEIGHT
780 C1234 FORMAT(' HARSCA: MONTE-CARLO WEIGHT FUNCTION H/HMAX GT 1 !',/
781 C * ' H = SUM OVER A,B FOR PROCESS M OF:',/
782 C * ' E(M)*ALPHAS**2*XA*FA(XA,Q**2)*XB*FB(XB,Q**2)',/
783 C * ' F(=H),XSECT(2,MSPR)(=HMAX), WEIGHT = H/HMAX',3E12.5)
784 C-------------------------------------------------------------------
785 C I
786 C end of iteration loop I
787 C -------------------------------------------I
788 C
789 C the event is accepted now
790 C
791 C actualize counter for accepted events
792  mxsect(2,mspr) = mxsect(2,mspr)+1
793  IF ( mspr.EQ.-1 ) mspr = 3
794 C find initial partons
795  sum = 0.0
796  scheck = rndm(ai)*pds
797  IF ( mspr.EQ.1 .OR. mspr.EQ.4 ) THEN
798  ia = 0
799  ib = 0
800  ELSEIF ( mspr.EQ.2 .OR. mspr.EQ.5 .OR. mspr.EQ.6 ) THEN
801  DO 610 ia=-nf,nf
802  IF ( ia.EQ.0 ) goto 610
803  sum = sum+pda(ia)*pdb(-ia)
804  IF ( sum.GE.scheck ) goto 620
805 610 CONTINUE
806 620 ib =-ia
807  ELSEIF ( mspr.EQ.3 ) THEN
808  ib = 0
809  DO 630 ia=-nf,nf
810  IF ( ia.EQ.0 ) goto 630
811  sum = sum+pda(0)*pdb(ia)
812  IF ( sum.GE.scheck ) goto 640
813  sum = sum+pda(ia)*pdb(0)
814  IF ( sum.GE.scheck ) goto 650
815 630 CONTINUE
816 640 ib = ia
817  ia = 0
818 650 CONTINUE
819  ELSEIF ( mspr.EQ.7 ) THEN
820  DO 660 ia=-nf,nf
821  IF ( ia.EQ.0 ) goto 660
822  sum = sum+pda(ia)*pdb(ia)
823  IF ( sum.GE.scheck ) goto 670
824 660 CONTINUE
825 670 ib = ia
826  ELSEIF ( mspr.EQ.8 ) THEN
827  DO 690 ia=-nf,nf
828  IF ( ia.EQ.0 ) goto 690
829  DO 680 ib=-nf,nf
830  IF ( abs(ib).EQ.abs(ia) .OR. ib.EQ.0 ) goto 680
831  sum = sum+pda(ia)*pdb(ib)
832  IF ( sum.GE.scheck ) goto 700
833 680 CONTINUE
834 690 CONTINUE
835 700 CONTINUE
836  ENDIF
837 C find final partons
838  ic = ia
839  id = ib
840  IF ( mspr.EQ.2 ) THEN
841  ic = 0
842  id = 0
843  ELSEIF ( mspr.EQ.4 ) THEN
844  ic = int(float(nf+nf)*rndm(ai))+1
845  IF ( ic.GT.nf ) ic = nf-ic
846  id =-ic
847  ELSEIF ( mspr.EQ.6 ) THEN
848  ic = int(float(nf+nf-2)*rndm(ai))+1
849  IF ( ic.GT.nf-1 ) ic = nf-1-ic
850  IF ( abs(ic).EQ.abs(ia) ) ic = sign(nf,ic)
851  id =-ic
852  ENDIF
853 C
854 30 a1 = rndm(ai)
855  a2 = rndm(ai)
856  IF ( ((a1*a1)+(a2*a2)).GT.1.0d0 ) goto 30
857  cosphi = ((a1*a1)-(a2*a2))/((a1*a1)+(a2*a2))
858  sinphi = sign(((a1*a2)+(a1*a2))/((a1*a1)+(a2*a2)),rndm(ai)-0.5)
859 C
860  IF ( rndm(ai)*pda(ia).GT.pda(-ia) ) ia = sign(abs(ia)+10,ia)
861  IF ( rndm(aj)*pdb(ib).GT.pdb(-ib) ) ib = sign(abs(ib)+10,ib)
862 C fill event record
863  line = line+1
864  prec(1,line) = 0.0
865  prec(2,line) = 0.0
866  prec(3,line) = 0.5*ecm*x1
867  prec(0,line) = prec(3,line)
868  lrec1(line) = ia+50+100*mspr
869  lrec2(line) = 01000
870  line = line+1
871  prec(1,line) = 0.0
872  prec(2,line) = 0.0
873  prec(3,line) =-0.5*ecm*x2
874  prec(0,line) =-prec(3,line)
875  lrec1(line) = ib+50
876  lrec2(line) = 01000
877  line = line+1
878  prec(1,line) = pt*cosphi
879  prec(2,line) = pt*sinphi
880  prec(3,line) =-0.5*ecm*(u*x1-v*x2)
881  prec(0,line) =-0.5*ecm*(u*x1+v*x2)
882  lrec1(line) = ic+50
883  lrec2(line) = 11000
884  line = line+1
885  prec(1,line) =-pt*cosphi
886  prec(2,line) =-pt*sinphi
887  prec(3,line) =-0.5*ecm*(v*x1-u*x2)
888  prec(0,line) =-0.5*ecm*(v*x1+u*x2)
889  lrec1(line) = id+50
890  lrec2(line) = 11000
891  RETURN
892  END
893 C_______________________________________________________________________
894  SUBROUTINE haoutp
895  IMPLICIT DOUBLE PRECISION(a-h,o-z)
896  SAVE
897  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
898  COMMON /haoutl/ noutl,nouter,noutco
899  COMMON /haevnt/ pt1,pt2,nhard,ntry,ihard,itry,irejev
900  COMMON /haevtr/ line,lin,lrec1(mline),lrec2(mline),prec(0:3,mline)
901  COMMON /harslt/ lscahd,lsc1hd,
902  & etahd(mscahd,2) ,pthd(mscahd),
903  & xhd(mscahd,2) ,vhd(mscahd) ,x0hd(mscahd,2),
904  & ninhd(mscahd,2) ,nouthd(mscahd,2),
905  & n0inhd(mscahd,2),nbrahd(mscahd,2),nprohd(mscahd)
906 
907 C output of data for hard scattering
908  IF ( noutl.GE.4 ) THEN
909  WRITE(6,1010) nhard,ihard,irejev
910 1010 FORMAT(' ===HARD EVENT=== NHARD,NTRUE,REJECTIONS ',3i5,/
911  &' IA IB IC ID XA XB PT YC YD',
912  &' PHI')
913  DO 10 n=1,lscahd
914  phi = atan2(prec(1,4*n-1),prec(2,4*n-1))
915  WRITE(6,1020) ninhd(n,1),ninhd(n,2),nouthd(n,1),nouthd(n,2),
916  & xhd(n,1),xhd(n,2),pthd(n),etahd(n,1),etahd(n,2),phi
917 1020 FORMAT(1x,4i3,2f11.7,4f9.3)
918 10 CONTINUE
919  ENDIF
920  IF ( noutl.GE.6 ) THEN
921 C output of eventrecord
922  WRITE(6,1030)
923 1030 FORMAT(' EVENTRECORD')
924  DO 20 l=1,line
925  WRITE(6,1040) lrec1(l),lrec2(l),(prec(i,l),i=0,3)
926 20 CONTINUE
927 1040 FORMAT(2i12,4(1pe12.4))
928  WRITE(6,1050)
929 1050 FORMAT(/)
930  ENDIF
931  RETURN
932  END
933 C_____________________________________________________________________
934 C
935 C original title: JTPADI.FOR
936 C_______________________________________________________________________
937 *
938 * ********************************************************************
939  SUBROUTINE jtpdis(X,QQ,IHATYP,MSPR,PD)
940 *
941 * Parton distributios
942 * NPD=ISTRUF (elsewhere)
943 *
944  IMPLICIT DOUBLE PRECISION(a-h,o-z)
945  SAVE
946  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
947  dimension pd(-6:6)
948  DATA iset / 0 /
949  DO 10 i=-6,6
950 10 pd(i) = 0.0
951 C SET MAXFL - THE MAX. NUMBER OF FLAVORS TO CALCULATE
952  maxfl = nf
953  IF ( mspr.EQ.1 .OR. mspr.EQ.4 ) maxfl = 0
954 C------------------------------------------------------
955 C CALL ROUTINES CALCULATING THE DISTRIBUTIONS
956 C EHLQ 1/2, MRS 1/2/3, GRV, HMRS 1/2, KMRS 1/2/3/4,
957 C WHERE HMRS 1/2, KMRS 1/2/3/4 CORRESP. TO HKMRS 1-6
958 C------------------------------------------------------
959  IF ( npd.EQ.1 .OR. npd.EQ.2 ) THEN
960 C CALL PDEHLQ(X,QQ,MAXFL,PD)
961  WRITE(6,*) ' unsupported PDF number: ',npd
962  ELSEIF ( npd.GE.3 .AND. npd.LE.5 ) THEN
963 C CALL PDMRS(X,QQ,MAXFL,PD)
964  WRITE(6,*) ' unsupported PDF number: ',npd
965  ELSEIF(npd.EQ.6)THEN
966 C CALL PDGRV(X,QQ,PD)
967  WRITE(6,*) ' unsupported PDF number: ',npd
968  ELSEIF(npd.EQ.7)THEN
969 C CALL PHKMRS(X,QQ,PD,1)
970  WRITE(6,*) ' unsupported PDF number: ',npd
971  ELSEIF(npd.EQ.8)THEN
972 C CALL PHKMRS(X,QQ,PD,2)
973  WRITE(6,*) ' unsupported PDF number: ',npd
974  ELSEIF(npd.EQ.9)THEN
975 C CALL PHKMRS(X,QQ,PD,3)
976  WRITE(6,*) ' unsupported PDF number: ',npd
977  ELSEIF(npd.EQ.10)THEN
978 C CALL PHKMRS(X,QQ,PD,4)
979  WRITE(6,*) ' unsupported PDF number: ',npd
980  ELSEIF(npd.EQ.11)THEN
981 C CALL PHKMRS(X,QQ,PD,5)
982  WRITE(6,*) ' unsupported PDF number: ',npd
983  ELSEIF(npd.EQ.12)THEN
984 C CALL PHKMRS(X,QQ,PD,6)
985  WRITE(6,*) ' unsupported PDF number: ',npd
986 * updates of April 92, April 93
987  ELSEIF((npd.GE.13).AND.(npd.LE.20)) THEN
988 C CALL PHKMRS(X,QQ,PD,NPD-6)
989  WRITE(6,*) ' unsupported PDF number: ',npd
990  ELSEIF((npd.GE.21).AND.(npd.LE.23)) THEN
991  CALL phkmrs(x,qq,pd,npd-6)
992  ELSE
993  WRITE(6,*) ' unsupported PDF number: ',npd
994  stop
995  ENDIF
996  DO 20 i=-maxfl,maxfl
997  IF ( pd(i).LT.1.d-15 ) pd(i) = 0.0
998 20 CONTINUE
999 C IF ANTIPROTON CHANGE QUARK <---> ANTIQUARK
1000  IF ( ihatyp.EQ.-1 ) THEN
1001  DO 50 i=1,6
1002  tttt = pd(-i)
1003  pd(-i) = pd(i)
1004 50 pd( i) = tttt
1005  ENDIF
1006  RETURN
1007  END
1008 
1009 C_______________________________________________________________________
1010  SUBROUTINE phkmrs(XQ,QQ,PD,MODE)
1011 C***************************************************************C
1012 C C
1013 C ORIGINAL NAME: MRSEB( ... ) C
1014 C C
1015 C ----- VARIABLE QUARKS AND GLUONS AT SMALL X ---- C
1016 C C
1017 C NEW VERSIONS !!!! JANUARY 1990 (AS DESCRIBED IN C
1018 C "PARTON DISTRIBUTIONS ... " P.N. HARRIMAN, A.D. MARTIN, C
1019 C R.G. ROBERTS AND W.J. STIRLING PREPRINT DTP-90-04 ) C
1020 C C
1021 C NEW VERSIONS !!!! JULY 1990 C
1022 C "........................ " J. KWIECINSKI, A.D. MARTIN, C
1023 C R.G. ROBERTS AND W.J. STIRLING PREPRINT DTP-90-46 ) C
1024 C
1025 C****************************************************************
1026 C All modes 1-14 dropped in dpmjet-II.4.2
1027 C***************************************************************
1028 C MODE 1 CORRESPONDS TO HARRIMAN, C
1029 C MARTIN, ROBERTS, STIRLING (EMC FIT) WITH LAMBDA= 100 MEV C
1030 C ORIGINAL: STRC27 NOW: PHMRS1 C
1031 C C
1032 C MODE 2 CORRESPONDS TO HARRIMAN, C
1033 C MARTIN, ROBERTS, STIRLING (BCDMS FIT) WITH LAMBDA= 190 MEV C
1034 C WITH SMALL X BEHAVIOUR DETERMINED FROM FIT "HB FIT" C
1035 C ORIGINAL: STRC28 NOW: PHMRS2 C
1036 C C
1037 C MODE 3 CORRESPONDS TO KWIECINSKI, C
1038 C MARTIN, ROBERTS, STIRLING (BCDMS FIT) WITH LAMBDA= 190 MEV C
1039 C AND XG,XQ --> CONSTANT AS X--> 0 AT Q0**2 "B0 FIT" C
1040 C ORIGINAL: STRC38 NOW: PKMRS1 C
1041 C C
1042 C MODE 4 CORRESPONDS TO KWIECINSKI, C
1043 C MARTIN, ROBERTS, STIRLING (BCDMS FIT) WITH LAMBDA= 190 MEV C
1044 C AND XG,XQ --> X**-1/2 AS X--> 0 AT Q0**2 "B- FIT" C
1045 C ORIGINAL: STRC48 NOW: PKMRS2 C
1046 C C
1047 C MODE 5 CORRESPONDS TO KWIECINSKI, C
1048 C MARTIN, ROBERTS, STIRLING (BCDMS FIT) WITH LAMBDA= 190 MEV C
1049 C AND XG,XQ --> X**-1/2 AS X--> 0 AT Q0**2 "B-(5) FIT" C
1050 C I.E. WITH WEAK (R=5 GEV-1) SHADOWING INCLUDED C
1051 C ORIGINAL: STRC58 NOW: PKMRS3 C
1052 C C
1053 C MODE 6 CORRESPONDS TO KWIECINSKI, C
1054 C MARTIN, ROBERTS, STIRLING (BCDMS FIT) WITH LAMBDA= 190 MEV C
1055 C AND XG,XQ --> X**-1/2 AS X--> 0 AT Q0**2 "B-(2) FIT" C
1056 C I.E. WITH STRONG (R=2 GEV-1) SHADOWING INCLUDED C
1057 C ORIGINAL: STRC68 NOW: PKMRS4 C
1058 C C
1059 C****************************************************************
1060 C All modes 1-14 dropped in dpmjet-II.4.2
1061 C***************************************************************
1062 C C
1063 C -*- C
1064 C C
1065 C (NOTE THAT X TIMES THE PARTON DISTRIBUTION FUNCTION C
1066 C IS RETURNED I.E. G(X) = GLU/X ETC, AND THAT "SEA" C
1067 C IS THE LIGHT QUARK SEA I.E. UBAR(X)=DBAR(X) C
1068 C = SEA/X FOR A PROTON. IF IN DOUBT, CHECK THE C
1069 C MOMENTUM SUM RULE! NOTE ALSO THAT SCALE=Q**2 IN GEV**2) C
1070 C C
1071 C -*- C
1072 C C
1073 C***************************************************************C
1074  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1075  SAVE
1076 C REAL XQ,QQ,PD
1077 C to get the same as outside
1078 C in the D IMPLICIT DOUBLE PRECISION(A-H,O-Z) world
1079  real
1080  & *8
1081  & xq,qq,pd
1082  dimension pd(-6:6)
1083  dimension pdff(-6:2)
1084  x=dble( xq )
1085  scale=dble( qq )
1086 C-------------------------------------------------------------------
1087 C****************************************************************
1088 C All modes 1-14 dropped in dpmjet-II.4.2
1089 C***************************************************************
1090 C****************************************************************
1091 C All modes 1-14 dropped in dpmjet-II.4.2
1092 C***************************************************************
1093 C--------------------------------------------------------------------
1094 C updates Oct 98 replace DOR94LO by PO_GRV98LO (R.Engel)
1095  IF((mode.EQ.15)) THEN
1096  scale2=scale
1097 C CALL DOR94LO(X,SCALE2,UV, DV, DEL, UDB, SB, GL)
1098  CALL po_grv98lo(iset,x,scale2,uv,dv,us,ds,ss,gl)
1099  cb = 0.d0
1100  bb = 0.d0
1101  pd(-5) = bb
1102  pd(-4) = cb
1103  pd(-3) = ss
1104  pd(-2) = us
1105  pd(-1) = ds
1106  pd(0) = gl
1107  pd(1) = dv+ds
1108  pd(2) = uv+us
1109  pd(3) = ss
1110  pd(4) = pd(-4)
1111  pd(5) = pd(-5)
1112  xq= x
1113  qq= scale
1114 C WRITE(6,*)' PD ',PD
1115  RETURN
1116  ENDIF
1117 C updates Oct 98 replace DOR94LO by PO_GRV98LO (R.Engel)
1118  IF((mode.EQ.16)) THEN
1119  scale2=scale
1120 C CALL DOR94LO(X,SCALE2,UV, DV, DEL, UDB, SB, GL)
1121  CALL po_grv98lo(iset,x,scale2,uv,dv,us,ds,ss,gl)
1122  cb = 0.d0
1123  bb = 0.d0
1124  pd(-5) = bb
1125  pd(-4) = cb
1126  pd(-3) = ss
1127  pd(-2) = us
1128  pd(-1) = ds
1129  pd(0) = gl
1130  pd(1) = dv+ds
1131  pd(2) = uv+us
1132  pd(3) = ss
1133  pd(4) = pd(-4)
1134  pd(5) = pd(-5)
1135  xq= x
1136  qq= scale
1137 C WRITE(6,*)' PD ',PD
1138  RETURN
1139  ENDIF
1140 C--------------------------------------------------------------------
1141 C updates Feb. 97
1142  IF((mode.EQ.17)) THEN
1143  CALL structm(x,scale,upv,dnv,usea,dsea,str,chm,bot,top,glu)
1144  pd(0)= glu
1145  pd(1)= usea+upv
1146  pd(2)= dsea+dnv
1147  pd(3)= str
1148  pd(4)= chm
1149  pd(5)= bot
1150  pd(-5)= bot
1151  pd(-4)= chm
1152  pd(-3)= str
1153  pd(-2)= dsea
1154  pd(-1)= usea
1155  xq= x
1156  qq= scale
1157  RETURN
1158  ENDIF
1159 C--------------------------------------------------------------------
1160 C--------------------------------------------------------------------
1161  pd(0)= glu
1162  pd(1)= sea+upv
1163  pd(2)= sea+dnv
1164  pd(3)= str
1165  pd(4)= chm
1166  pd(5)= bot
1167  pd(-5)= bot
1168  pd(-4)= chm
1169  pd(-3)= str
1170  pd(-2)= sea
1171  pd(-1)= sea
1172  xq= x
1173  qq= scale
1174 C--------------------------------------------------------------------
1175  RETURN
1176  END
1177 C
1178 C*********************************************************************
1179 C-----original seperate file with the name------------------------------------------------------------------
1180 C DTULPTPE.FOR
1181 C
1182 C------------------------------------------------------------------------
1183 C_______________________________________________________________________
1184 C
1185 C PROGRAM FOR SIMULATION OF HARD SCATTERING OF HADRONIC PARTICLES
1186 C
1187 C AUTHOR K.HAHN LEIPZIG GDR
1188 C_______________________________________________________________________
1189  SUBROUTINE laptab
1190  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1191  SAVE
1192  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
1193  COMMON /hacons/ pi,pi2,pi4,gevtmb
1194  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
1195  COMMON /hapadi/ npdm
1196  COMMON /hapdco/ npdcor
1197  COMMON /haqqap/ aqqal,aqqpd,nqqal,nqqpd
1198  COMMON /haenvi/ nindep
1199  COMMON /haoutl/ noutl,nouter,noutco
1200  COMMON /hacuts/ ptl,ptu,etacl,etacu,etadl,etadu
1201  COMMON /hagaup/ ngaup1,ngaup2,ngauet,ngauin
1202  COMMON /haevtr/ line,lin,lrec1(mline),lrec2(mline),prec(0:3,mline)
1203  COMMON /harslt/ lscahd,lsc1hd,
1204  & etahd(mscahd,2) ,pthd(mscahd),
1205  & xhd(mscahd,2) ,vhd(mscahd) ,x0hd(mscahd,2),
1206  & ninhd(mscahd,2) ,nouthd(mscahd,2),
1207  & n0inhd(mscahd,2),nbrahd(mscahd,2),nprohd(mscahd)
1208  INTEGER mxsect
1209  COMMON /haxsec/ xsecta(2,-1:maxpro,4),xsect(5,-1:maxpro),
1210  & mxsect(0:2,-1:maxpro)
1211  CHARACTER*8 cw
1212  CHARACTER*79 commnt
1213  CHARACTER*70 inp
1214  dimension what(6)
1215 C======================================================================
1216 C THE AVAILABLE CODE WORDS ARE
1217 C
1218 C END , COMMENT , ENERGYPT, PARDISTR, CUTS
1219 C INTPOINT, FLAVOR , PARTICLE, OUTPUT , INIT ,
1220 C TESTINCL, TESTMC , SUBPRON , SUBPROFF, HISOUT ,
1221 C HISINI , HARDSCAL, USER , PARDISCO
1222 C_______________________________________________________________________
1223  WRITE(6,2000)
1224 2000 FORMAT('1***************************************************'
1225  & ,/, ' MONTE-CARLO GENERATION OF HARD HADRONIC SCATTERINGS'
1226  & ,/, ' ***************************************************',/)
1227  CALL jtdtu(1)
1228 C read a control card
1229 10 CONTINUE
1230  READ (5,1010) inp
1231  IF ( inp(1:1).EQ.'-' ) goto 10
1232  WRITE(6,1011) inp
1233  READ(inp,1012,err=99) cw,what
1234  goto 15
1235 99 WRITE(6,1013)
1236  goto 10
1237 1010 FORMAT(a70)
1238 1011 FORMAT(' *********.* CONTROL.CARD*****.',4(9x,'.'),/,1x,a70,/)
1239 1012 FORMAT(a8,2x,6e10.0)
1240 1013 FORMAT(' CARD IS INCORRECT, IGNORE AND TRY NEXT CARD',/)
1241 15 CONTINUE
1242 C======================================================================
1243 C======================================================================
1244  IF ( cw.EQ.'END ' ) THEN
1245 C
1246 C STOPS THE PROGRAM
1247 C_______________________________________________________________________
1248  WRITE(6,1030)
1249 1030 FORMAT(' ******** END OF PROGRAM EXECUTION ********')
1250  RETURN
1251 C======================================================================
1252  ELSEIF ( cw.EQ.'COMMENT ' ) THEN
1253 C
1254 C TO INCLUDE COMMENTS IN PROGRAM OUTPUT
1255 C WHAT(1) - # OF CARDS FOLLOWING THE CONTROLCARD AND
1256 C CONTAINING THE COMMENT
1257 C_______________________________________________________________________
1258  n = max(1,int(what(1)))
1259  DO 20 i=1,n
1260  READ(5,1040) commnt
1261 20 WRITE(6,1050) commnt
1262 1040 FORMAT(a79)
1263 1050 FORMAT(1x,a79)
1264 C======================================================================
1265  ELSEIF ( cw.EQ.'ENERGYPT' ) THEN
1266 C
1267 C READ CMS-ENERGY AND MIN. TRANSVERSE MOMENTUM FOR JETS
1268 C WHAT(1) - CMS-ENERGY ( IN GEV ) DEFAULT 540.
1269 C WHAT(2) - MIN. PT ( IN GEV ) DEFAULT 2.
1270 C WHAT(3) - MIN. PT ( IN GEV ) DEFAULT 2.
1271 C WHAT(4) - MIN. PT ( IN GEV ) DEFAULT 2.
1272 C WHAT(5) - MIN. PT ( IN GEV ) DEFAULT 2.
1273 C
1274 C_______________________________________________________________________
1275  IF ( what(1).GT.0.0d0 ) ecm = what(1)
1276  DO 22 i=1,4
1277  ptini(i) = what(i+1)
1278 22 CONTINUE
1279 C======================================================================
1280  ELSEIF ( cw.EQ.'PARDISTR' ) THEN
1281 C
1282 C DEFINE THE PARTON DISTRIBUTION SET
1283 C WHAT(1) - NPD DEFAULT 3.
1284 C NPD = 1,2 : EICHTEN,HINCHLIFFE,LANE,QUIGG
1285 C = 3,4,5 : MARTIN,ROBERTS,STIRLING
1286 C = 6 : GLUECK,REYA,VOGT
1287 C = 7,8 : HARRIMAN,MARTIN,ROBERTS,STIRLING
1288 C = 9 - 12 : KWIECINSKI,MARTIN,ROBERTS,STIRLING
1289 C
1290 C WHAT(2) - NPDM DEFAULT 0.
1291 C NPDM = 0,1 : CORRECTS THE PARTON DISTRIBUTION IN CASE
1292 C OF NPD = 10 TO THE 1/SQRT(X) BEHAVIOUR
1293 C FOR X < XMIN=1.E-05 IF NPDM = 1 AND NOT
1294 C IF NPDM = 0 ( M STANDS FOR MODIFICATION )
1295 C_______________________________________________________________________
1296  ipd = int(what(1))
1297  ipdm = int(what(2))
1298  npd = 3
1299  npdm = 0
1300 * IF ( IPD.GE.1 .AND. IPD.LE.12 ) NPD = IPD
1301  IF ( ipd.GE.1 .AND. ipd.LE.15 ) npd = ipd
1302  IF ( ipdm.EQ.1 ) npdm = ipdm
1303 C======================================================================
1304  ELSEIF ( cw.EQ.'CUTS ' ) THEN
1305 C
1306 C set kinematic cuts on partonlevel
1307 C WHAT(1) - PTL min. pt for parton DEFAULT 0.0
1308 C WHAT(2) - PTU max. pt for parton DEFAULT 1.E+30
1309 C WHAT(3) - ETACL min. rapidity for parton c DEFAULT -1.E+30
1310 C WHAT(4) - ETACU max. rapidity for parton c DEFAULT 1.E+30
1311 C WHAT(5) - ETADL min. rapidity for parton d DEFAULT -1.E+30
1312 C WHAT(6) - ETADU max. rapidity for parton d DEFAULT 1.E+30
1313 C
1314 C this cuts are additional cuts which not affect the pt-cut given
1315 C in the ENERGYPT card;
1316 C this cuts are checked during event generation and events violating
1317 C the cuts are rejected;
1318 C the defaults are such that there is really no cutting;
1319 C_______________________________________________________________________
1320  ptl = what(1)
1321  ptu = what(2)
1322  etacl = what(3)
1323  etacu = what(4)
1324  etadl = what(5)
1325  etadu = what(6)
1326  IF ( ptu .LE.ptl ) ptu = ptl +1.0
1327  IF ( etacu.LE.etacl ) etacu = etacl+1.0
1328  IF ( etadu.LE.etadl ) etadu = etadl+1.0
1329 C======================================================================
1330  ELSEIF ( cw.EQ.'INTPOINT' ) THEN
1331 C
1332 C define number of integration points
1333 C WHAT(1) - NGAUP1 DEFAULT 8.
1334 C WHAT(2) - NGAUP2 DEFAULT 8.
1335 C WHAT(3) - NGAUET DEFAULT 8.
1336 C WHAT(4) - NGAUIN DEFAULT 8.
1337 C
1338 C NGAUP1,NGAUP2,NGAUET : used in inclusive x-section calculation;
1339 C first, second pt-integration and
1340 C eta-integration respectivly
1341 C NGAUIN : used in initialization ( SUBROUTINE HABINT )
1342 C
1343 C max. value allowed for integration is 32 ! - For NGAUP1,P2,ET
1344 C max. value allowed for integration is 1000 ! - For NGAUIN.
1345 C (4.6.91)
1346 C_______________________________________________________________________
1347  IF ( what(1).GE.1.d0.AND.what(1).LE.32.d0) ngaup1 = int(what(1))
1348  IF ( what(2).GE.1.d0.AND. what(2).LE.32.d0) ngaup2 = int(what(2))
1349  IF ( what(3).GE.1.d0.AND. what(3).LE.32.d0) ngauet = int(what(3))
1350  IF(what(4).GE.1.d0.AND. what(3).LE.1000.d0) ngauin = int(what(4))
1351 C======================================================================
1352  ELSEIF ( cw.EQ.'FLAVOR ' ) THEN
1353 C DEFINES ACTIVE FLAVORS
1354 C WHAT(1) - # OF FLAVORS DEFAULT 4
1355 C_______________________________________________________________________
1356  nff = int(what(1))
1357  IF ( nff.GE.0 .AND. nff .LE.6 ) nf = nff
1358 C======================================================================
1359  ELSEIF ( cw.EQ.'PARTICLE' ) THEN
1360 C
1361 C TARGET AND BEAM PARTICLES
1362 C WHAT(1) - BEAM ( POSITIVE Z-DIRECTION ) DEFAULT 1
1363 C WHAT(2) - TARGET DEFAULT 1
1364 C 1 : PROTON
1365 C -1 : ANTIPROTON
1366 C_______________________________________________________________________
1367  iha = int(what(1))
1368  IF ( abs(iha).EQ.1 ) nha = iha
1369  ihb = int(what(2))
1370  IF ( abs(ihb).EQ.1 ) nhb = ihb
1371 C======================================================================
1372  ELSEIF ( cw.EQ.'OUTPUT ' ) THEN
1373 C
1374 C set output level
1375 C WHAT(1) - NOUTL output level 0.... DEFAULT 1.
1376 C WHAT(2) - NOUTER error messages ( 0 - no ; 1 - yes ) DEFAULT 1.
1377 C_______________________________________________________________________
1378  IF ( what(1).GE.0.d0 ) noutl = int(what(1))
1379  IF ( what(2).EQ.0.d0.OR.what(2).EQ.1.d0)nouter = int(what(2))
1380  IF ( what(3).GE.0.d0 ) noutco = int(what(3))
1381 C======================================================================
1382  ELSEIF ( cw.EQ.'INIT ' ) THEN
1383 C
1384 C DEMANDS A NEW INITIALIZATION
1385 C_______________________________________________________________________
1386  CALL harini
1387 C======================================================================
1388  ELSEIF ( cw.EQ.'TESTINCL' ) THEN
1389 C
1390 C TEST OF INCLUSIVE JET PRODUCTION
1391 C PLOT OF PARTONDISTRIBUTIONS USED ( GLUONS )
1392 C PLOT OF DIFFERENTIAL JET CROSS SECTIONS
1393 C
1394 C_______________________________________________________________________
1395  DO 35 i=1,6
1396  j = int(what(i))
1397  IF ( j.GE.1 .AND. j.LE.4 ) CALL hatest(j)
1398 35 CONTINUE
1399 C======================================================================
1400  ELSEIF ( cw.EQ.'TESTMC ' ) THEN
1401 C
1402 C TEST OF MONTE-CARLO JET PRODUCTION
1403 C PLOT OF DIFFERENTIAL JET CROSS SECTIONS
1404 C WHAT(1) - # OF EVENTS TO PRODUCE IN TEST DEFAULT 100
1405 C WHAT(2) - # OF HARD SCATTERINGS AT HARD EVENT
1406 C WHAT(3) - MIN. PT FOR FIRST HARD SCATTERING
1407 C
1408 C_______________________________________________________________________
1409  nevt = int(what(1))
1410  IF ( nevt.LE.0 ) nevt = 100
1411  nhard = max(1,int(what(2)))
1412  pt1 = what(3)
1413  CALL timdat
1414  DO 36 i=1,nevt
1415  mhard = nhard
1416  CALL harevt(mhard,pt1)
1417 36 CONTINUE
1418  CALL timdat
1419 C======================================================================
1420  ELSEIF ( cw.EQ.'SUBPRON ' ) THEN
1421 C
1422 C SWITCH SUBPROCESSES ON
1423 C
1424 C______________________________________________________________________
1425  DO 40 i=1,6
1426  m = int(what(i))
1427  IF ( m.GE.1 .AND. m.LE.maxpro ) mxsect(0,m) = 1
1428 40 CONTINUE
1429  mxsect(0,-1) = mxsect(0,3)
1430 C======================================================================
1431  ELSEIF ( cw.EQ.'SUBPROFF' ) THEN
1432 C
1433 C SWITCH SUBPROCESSES OFF
1434 C
1435 C_______________________________________________________________________
1436  DO 50 i=1,6
1437  m = int(what(i))
1438  IF ( m.GE.1 .AND. m.LE.maxpro ) mxsect(0,m) = 0
1439 50 CONTINUE
1440  mxsect(0,-1) = mxsect(0,3)
1441 C======================================================================
1442  ELSEIF ( cw.EQ.'HISOUT ' ) THEN
1443 C
1444 C OUTPUT OF MC RESULTS
1445 C WHAT()- 1 : TOTAL CROSS SECTION
1446 C WHAT()- 2 : PT-DISTR.
1447 C WHAT()- 3 : PT-DISTR.
1448 C WHAT()- 4 : PT-DISTR.
1449 C WHAT()- 5 : RAPIDITY-DISTR.
1450 C WHAT()- 6 : RAPIDITY-DISTR.
1451 C
1452 C_______________________________________________________________________
1453  DO 60 i=1,6
1454  j = int(what(i))
1455  IF ( j.GE.1 .AND. j.LE.6 ) CALL hisout(j)
1456 60 CONTINUE
1457 C======================================================================
1458  ELSEIF ( cw.EQ.'HISINI ' ) THEN
1459 C
1460 C INITIALIZE THE HISTOGRAMS FOR MC TESTING
1461 C_______________________________________________________________________
1462  CALL hisini
1463 C======================================================================
1464  ELSEIF ( cw.EQ.'HARDSCAL' ) THEN
1465 C
1466 C define hard scale
1467 C WHAT(1) - NQQAL : definition of hard scale in coupling constant
1468 C WHAT(2) - AQQAL : factor multiplied to hard scale in coupling
1469 C WHAT(3) - NQQPD : definition of hard scale in parton distr.
1470 C WHAT(4) - AQQPD : factor multiplied to hard scale in parton distr.
1471 C the possible NQQAL(PD) are
1472 C 1 : QQ = AQQAL(PD) * PT**2
1473 C 2 : QQ = AQQAL(PD) * SP
1474 C 3 : QQ = AQQAL(PD) * (SP*TP*UP)**(1./3.)
1475 C 4 : QQ = AQQAL(PD) * (SP*TP*UP)/(SP**2+TP**2+UP**2)
1476 C
1477 C DEFAULT is NQQAL = NQQPD = 1 and AQQAL = AQQPD = 1./4.
1478 C_______________________________________________________________________
1479  IF ( what(1).GE.1.d0.AND.what(1).LE.4.d0)nqqal = int(what(1))
1480  IF ( what(2).GT.0.d0 ) aqqal = what(2)
1481  IF ( what(3).GE.1.d0.AND.what(3).LE.4.d0)nqqpd = int(what(3))
1482  IF ( what(4).GT.0.d0 ) aqqpd = what(4)
1483 C C======================================================================
1484 C ELSEIF ( CW.EQ.'USER ' ) THEN
1485 C C
1486 C C_______________________________________________________________________
1487 C WRITE(6,9050)
1488 C 9050 FORMAT(' -----> GIVE CONTROL TO USER ROUTINE')
1489 C CALL HAUSER(WHAT)
1490 C WRITE(6,9060)
1491 C 9060 FORMAT(' -----> CONTROL COMES BACK FROM USER ROUTINE')
1492 C======================================================================
1493  ELSEIF ( cw.EQ.'PARDISCO' ) THEN
1494 C
1495 C define partoncorrelationfunction
1496 C WHAT(1) - NPDCOR : =0 no correlations, =1 correlations
1497 C
1498 C DEFAULT is NPDCOR = 0
1499 C_______________________________________________________________________
1500  IF ( what(1).EQ.0.d0.OR.what(1).EQ.1.d0)npdcor = int(what(1))
1501 C======================================================================
1502  ELSE
1503  WRITE(6,9999)
1504 9999 FORMAT(' ##### UNKNOWN CODEWORD; CARD IS IGNORED ###',/)
1505  ENDIF
1506  goto 10
1507  END
1508 C
1509 C-----originally seperate file with the name ---------------------------
1510 C JTINCL.FOR
1511 C
1512 C______________________________________________________________________
1513 C
1514 C PROCEDURES FOR CALCULATION OF INCLUSIVE CROSS SECTIONS
1515 C
1516 C______________________________________________________________________
1517 C______________________________________________________________________
1518  SUBROUTINE csj2m(PT,ETAC,ETAD,DSIGMM)
1519 C CALCULATION OF DIFFERENTIAL CROSS SECTION DSIG/(DETAC*DETAD*DPT)
1520 C FOR DIFFERENT SUBPROCESSES
1521  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1522  SAVE
1523  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
1524  parameter( tiny= 1.d-30, onep1=1.1d0 ,tiny6=1.d-06)
1525  COMMON /hacons/ pi,pi2,pi4,gevtmb
1526  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
1527  COMMON /haqqap/ aqqal,aqqpd,nqqal,nqqpd
1528  DOUBLE PRECISION ec,ed,xa,xb,sp,tp,up,tt,uu,
1529  * factor,dsigm(0:maxpro)
1530  dimension dsigmm(0:maxpro),pda(-6:6),pdb(-6:6)
1531 
1532  DO 10 i=0,maxpro
1533  dsigm(i) = 0.0
1534 10 CONTINUE
1535  ec = exp(etac)
1536  ed = exp(etad)
1537 C KINETICS
1538  xa = pt*(ec+ed)/ecm
1539  xb = xa/(ec*ed)
1540  IF ( xa.GE.1.d0 .OR. xb.GE.1.d0 ) RETURN
1541  sp = xa*xb*ecm*ecm
1542  up =-ecm*pt*ec*xb
1543  up = up/sp
1544  tp =-(1.+up)
1545  uu = up*up
1546  tt = tp*tp
1547 C set hard scale QQ for alpha and partondistr.
1548  IF ( nqqal.EQ.1 ) THEN
1549  qqal = aqqal*pt*pt
1550  ELSEIF ( nqqal.EQ.2 ) THEN
1551  qqal = aqqal*sp
1552  ELSEIF ( nqqal.EQ.3 ) THEN
1553  qqal = aqqal*sp*(up*tp)**(1./3.)
1554  ELSEIF ( nqqal.EQ.4 ) THEN
1555  qqal = aqqal*sp*up*tp/(1.+tt+uu)
1556  ENDIF
1557  IF ( nqqpd.EQ.1 ) THEN
1558  qqpd = aqqpd*pt*pt
1559  ELSEIF ( nqqpd.EQ.2 ) THEN
1560  qqpd = aqqpd*sp
1561  ELSEIF ( nqqpd.EQ.3 ) THEN
1562  qqpd = aqqpd*sp*(up*tp)**(1./3.)
1563  ELSEIF ( nqqpd.EQ.4 ) THEN
1564  qqpd = aqqpd*sp*up*tp/(1.+tt+uu)
1565  ENDIF
1566 C PARAMETER ( TINY= 1.D-30, ONEP1=1.1D0 ,TINY6=1.D-06)
1567  alpha = bqcd/log(max(qqal/alasqr,onep1))
1568  factor = pi2*gevtmb*pt*(alpha/sp)**2
1569 C PARTONDISTRIBUTIONS ( MULTIPLIED BY X )
1570  x1 = xa
1571  x2 = xb
1572  CALL jtpdis(x1,qqpd,nha,0,pda)
1573  CALL jtpdis(x2,qqpd,nhb,0,pdb)
1574  s1 = pda(0)*pdb(0)
1575  s2 = 0.0
1576  s3 = 0.0
1577  s4 = 0.0
1578  s5 = 0.0
1579  DO 20 i=1,nf
1580  s2 = s2+pda(i)*pdb(-i)+pda(-i)*pdb( i)
1581  s3 = s3+pda(i)*pdb( i)+pda(-i)*pdb(-i)
1582  s4 = s4+pda(i)+pda(-i)
1583  s5 = s5+pdb(i)+pdb(-i)
1584 20 CONTINUE
1585 C CROSS-SECTIONS ( INCLUDING MATRIX ELEMENTS, SYMMETRY-FACTORS AND
1586 C FACTORS FOR FINALSTATE-SUMMATION )
1587  dsigm(1) = 2.25*(3.-((up*tp)+up/tt+tp/uu))
1588  dsigm(6) = (4./9.)*(uu+tt)
1589  dsigm(8) = (4./9.)*(1.+uu)/tt
1590  dsigm(2) = (16./27.)*(uu+tt)/(up*tp)-3.*dsigm(6)
1591  dsigm(3) = ((1.+uu)/tt)-(4./9.)*(1.+uu)/up
1592  dsigm(4) = (9./32.)*dsigm(2)
1593  dsigm(5) = dsigm(6)+dsigm(8)-(8./27.)*uu/tp
1594  dsigm(7) = 0.5*(dsigm(8)+(4./9.)*(1.+tt)/uu-(8./27.)/(up*tp))
1595 
1596  dsigm(1) = factor*dsigm(1)*s1
1597  dsigm(2) = factor*dsigm(2)*s2
1598  dsigm(3) = factor*dsigm(3)*(pda(0)*s5+pdb(0)*s4)
1599  dsigm(4) = factor*dsigm(4)*s1*nf
1600  dsigm(5) = factor*dsigm(5)*s2
1601  dsigm(6) = factor*dsigm(6)*s2*max(0,(nf-1))
1602  dsigm(7) = factor*dsigm(7)*s3
1603  dsigm(8) = factor*dsigm(8)*(s4*s5-(s2+s3))
1604 C sum over processes
1605  DO 40 m=1,maxpro
1606  dsigm(0) = dsigm(0)+dsigm(m)
1607 40 CONTINUE
1608  DO 50 m=0,maxpro
1609  dsigmm(m) = dsigm(m)
1610 50 CONTINUE
1611  RETURN
1612  END
1613 C______________________________________________________________________
1614  SUBROUTINE csj1m(PT,ETAC,DSIGM)
1615 C ONE JET CROSS SECTION
1616 C ( CSJ2M INTEGRATED OVER ETAD )
1617  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1618  SAVE
1619  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
1620 C PARAMETER ( TINY= 1.D-30 )
1621 C CHANGED 14.10.92 AFTER COMPARISON WITH THE ORIGINAL
1622 C VERSION OF DTULAP FROM 1.D-30 BACK TO THE ORIGINAL 1.D-20
1623  parameter( tiny= 1.d-20 )
1624  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
1625  COMMON /hagaup/ ngaup1,ngaup2,ngauet,ngauin
1626  dimension dsigm(0:maxpro),dsig1(0:maxpro)
1627  dimension absz(32),weig(32)
1628 
1629  DO 10 m=0,maxpro
1630  dsigm(m) = 0.0
1631 10 CONTINUE
1632  ec = exp(etac)
1633  arg = ecm/pt
1634  IF ( arg.LE.ec .OR. arg.LE.1./ec ) RETURN
1635  edu = log(arg-ec)
1636  edl =-log(arg-1./ec)
1637  npoint = ngauet
1638  CALL gset(edl,edu,npoint,absz,weig)
1639  DO 30 i=1,npoint
1640  CALL csj2m(pt,etac,absz(i),dsig1)
1641  DO 20 m=0,maxpro
1642 C PCTRL= DSIG1(M)/1.E-20
1643  pctrl= dsig1(m)/tiny
1644  pctrl= abs( pctrl )
1645  IF( pctrl.GE.1.d0 ) THEN
1646  dsigm(m) = dsigm(m)+weig(i)*dsig1(m)
1647  ENDIF
1648 20 CONTINUE
1649 30 CONTINUE
1650  RETURN
1651  END
1652 C______________________________________________________________________
1653  SUBROUTINE csj1mi(PT,DSIGM)
1654 C ONE JET CROSS SECTION
1655 C ( CSJ2M INTEGRATED OVER ETAD AND ETAC )
1656  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1657  SAVE
1658  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
1659  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
1660  COMMON /hagaup/ ngaup1,ngaup2,ngauet,ngauin
1661  dimension dsigm(0:maxpro),dsig1(0:maxpro)
1662  dimension absz(32),weig(32)
1663 
1664  DO 10 m=0,maxpro
1665  dsigm(m) = 0.0
1666 10 CONTINUE
1667  amt = 2.*pt/ecm
1668  IF ( amt.GE.1.d0 ) RETURN
1669  ecu = log((sqrt(1.-amt*amt)+1.)/amt)
1670  ecl =-ecu
1671  npoint = ngauet
1672  CALL gset(ecl,ecu,npoint,absz,weig)
1673  DO 30 i=1,npoint
1674  CALL csj1m(pt,absz(i),dsig1)
1675  DO 20 m=0,maxpro
1676  dsigm(m) = dsigm(m)+weig(i)*dsig1(m)
1677 20 CONTINUE
1678 30 CONTINUE
1679  RETURN
1680  END
1681 C______________________________________________________________________
1682  SUBROUTINE csharm(DSIGM)
1683 C
1684 C TOTAL HARD CROSS SECTION
1685 C ( CSJ2M INTEGRATED OVER PT, ETAD AND ETAC )
1686 C
1687  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1688  SAVE
1689  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
1690  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
1691  COMMON /hagaup/ ngaup1,ngaup2,ngauet,ngauin
1692  COMMON /xsecpt/ ptcut,sigs,dsigh
1693  dimension dsigm(0:maxpro),dsig1(0:maxpro)
1694  dimension absz(32),weig(32)
1695  DATA fac / 3.0 /
1696 
1697  DO 10 m=0,maxpro
1698  dsigm(m) = 0.0
1699 10 CONTINUE
1700  IF ( ptini(1).GE.ecm/2.d0 ) RETURN
1701  ptmin = ptini(1)
1702  ptcut = ptmin
1703  ptmax = min(fac*ptmin,ecm/2.)
1704  npoint = ngaup1
1705  CALL csj1mi(ptmin,dsig1)
1706  sig1 = dsig1(0)
1707 C WRITE(6,1000) SIG1
1708  1000 FORMAT(1x,' d sigma/ p_t d p_t ',e12.5)
1709  dsigh = sig1
1710  ptmxx = 0.95*ptmax
1711  CALL csj1mi(ptmxx,dsig1)
1712  ex = log(sig1/(dsig1(0)+1.e-30))/log(fac)
1713  ex1 = 1.0-ex
1714  DO 50 k=1,2
1715  IF ( ptmin.GE.ptmax ) goto 40
1716  rl = ptmin**ex1
1717  ru = ptmax**ex1
1718  CALL gset(rl,ru,npoint,absz,weig)
1719  DO 30 i=1,npoint
1720  r = absz(i)
1721  pt = r**(1.0/ex1)
1722  CALL csj1mi(pt,dsig1)
1723  f = weig(i)*pt/(r*ex1)
1724  DO 20 m=0,maxpro
1725  dsigm(m) = dsigm(m)+f*dsig1(m)
1726 20 CONTINUE
1727 30 CONTINUE
1728 40 ptmin = ptmax
1729  ptmax = ecm/2.0
1730  npoint = ngaup2
1731 50 CONTINUE
1732  RETURN
1733  END
1734 C-------originally seperate file ----------------------------------------
1735 C JTHIST.FOR
1736 C
1737 C________________________________________________________________________
1738 
1739 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1740 C
1741 C PROCEDURES TO PRINT OUT HISTOGRAMS AND INCLUSICE TESTCALCULATIONS
1742 C ( AND TO INITIALIZE AND FILL HISTOGRAMS ALSO )
1743 C
1744 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1745 C______________________________________________________________________
1746  SUBROUTINE hisout(IOUT)
1747  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1748  SAVE
1749  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
1750  parameter( tiny= 1.d-30, one=1.d0 ,tiny6=1.d-06,zero=0.)
1751  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
1752  INTEGER mxsect
1753  COMMON /haxsec/ xsecta(2,-1:maxpro,4),xsect(5,-1:maxpro),
1754  & mxsect(0:2,-1:maxpro)
1755  CHARACTER*18 proc
1756  CHARACTER*11 pdset,partic
1757  COMMON /peproc/ proc(0:maxpro),pdset(23),partic(-1:1)
1758 C LENGTH OF HISTO : 15000 REAL*4
1759  COMMON /histo / pt10,dpt1,eta10,deta1,pt20,dpt2,eta20,deta2,
1760  & x(50,-5:5),ab(50,-5:5),hpe(50,-5:5),hep(50,5),
1761  & hpm(50,8),hem(50,8),hp(50),he(50),
1762  & sig(maxpro),stdev(maxpro),
1763  & fill(12176)
1764 C
1765  xsmin =-6.
1766  xsstep= 0.1
1767  DO 5 j=-5,5
1768  DO 5 i=1,50
1769 5 x(i,j) =-100.
1770  xsect(2,0) = xsect(2,-1)
1771  mxsect(1,0) = mxsect(1,-1)
1772  mxsect(2,0) = mxsect(2,-1)
1773  DO 7 m=1,maxpro
1774  mxsect(1,0) = mxsect(1,0)+mxsect(1,m)
1775  mxsect(2,0) = mxsect(2,0)+mxsect(2,m)
1776 7 xsect(2,0) = xsect(2,0)+xsect(2,m)
1777  WRITE(6,1010) iout
1778 1010 FORMAT(1x,20('=='),' HISTO-OUTPUT ',i2,1x,10('=='),/)
1779  IF ( iout.EQ.1 ) THEN
1780  WRITE(6,1040)
1781 1040 FORMAT(' PROCESS',15x,'EVENTS',22x,'HARD CROSS SECTION',/,
1782  & 25x,'TOTAL ACCEPT.',10x,'MONTE-CARLO',11x,'INCLUSIVE')
1783  sigsum = 0.0
1784  stdevs = 0.0
1785  DO 20 m=1,maxpro
1786  sig(m) = 0.0
1787  stdev(m) = 0.0
1788  IF ( mxsect(1,m).GT.0 ) THEN
1789  sig(m) = xsect(3,m)/mxsect(1,m)
1790  stdev(m) = sqrt(max(zero,xsect(4,m)-xsect(3,m)*sig(m)))/
1791  & mxsect(1,m)
1792  ENDIF
1793  IF ( m.EQ.3 .AND. mxsect(1,-1).GT.0 ) THEN
1794  sigg = xsect(3,-1)/mxsect(1,-1)
1795  sig(3) = sig(3)+sigg
1796  stdev(3) = stdev(3)
1797  & +sqrt(max(zero,xsect(4,-1)-xsect(3,-1)*sigg))/mxsect(1,-1)
1798  ENDIF
1799  sigsum = sigsum+sig(m)
1800  stdevs = stdevs+stdev(m)
1801 20 CONTINUE
1802  mxsect(1,3) = mxsect(1,3)+mxsect(1,-1)
1803  mxsect(2,3) = mxsect(2,3)+mxsect(2,-1)
1804  WRITE(6,1050) proc(0),(mxsect(l,0),l=0,2),
1805  & sigsum,stdevs,xsect(5,0)
1806  DO 25 m=1,maxpro
1807  IF ( mxsect(0,m).EQ.1 ) WRITE(6,1050) proc(m),
1808  & (mxsect(l,m),l=0,2),sig(m),stdev(m),xsect(5,m)
1809 25 CONTINUE
1810 1050 FORMAT(a19,i3,2i8,e14.4,' +- ',e10.4,e14.4)
1811  mxsect(1,3) = mxsect(1,3)-mxsect(1,-1)
1812  mxsect(2,3) = mxsect(2,3)-mxsect(2,-1)
1813  ELSEIF ( iout.EQ.2 ) THEN
1814  fac = xsect(2,0)/(dpt1*mxsect(1,0))
1815  DO 30 i=1,50
1816  ab(i,1) = pt10+(i-1)*dpt1
1817  IF ( hp(i).GT.1.d-35 ) x(i,1) = log10(fac*hp(i))
1818 30 CONTINUE
1819  WRITE(6,1060)
1820 1060 FORMAT(' JET CROSS SECTION PT-DISTRIBUTION',/)
1821  CALL plot(ab(1,1),x(1,1),50,1,50,pt10,dpt1,xsmin,xsstep)
1822  ELSEIF ( iout.EQ.3 ) THEN
1823  fac = xsect(2,0)/(dpt1*mxsect(1,0))
1824  DO 50 i=1,50
1825  pt = pt10+(i-1)*dpt1
1826  DO 40 j=1,8
1827  ab(i,j-6) = pt
1828  IF ( hpm(i,j).GT.1.d-35 ) x(i,j-6) = log10(fac*hpm(i,j))
1829 40 CONTINUE
1830 50 CONTINUE
1831  WRITE(6,1070)
1832 1070 FORMAT(' JET CROSS SECTION PT-DISTRIBUTION',/,
1833  & ' FOR THE DIFF. SUBPROCESSES',/)
1834  CALL plot(ab,x,400,8,50,pt10,dpt1,xsmin,xsstep)
1835  ELSEIF ( iout.EQ.4 ) THEN
1836  fac = xsect(2,0)/(dpt1*deta1*mxsect(1,0))
1837  DO 70 i=1,50
1838  pt = pt10+(i-1)*dpt1
1839  DO 60 j=-5,5
1840  ab(i,j) = pt
1841  IF ( hpe(i,j).GT.1.d-35 ) x(i,j) = log10(fac*hpe(i,j))
1842 60 CONTINUE
1843 70 CONTINUE
1844  WRITE(6,1080) eta10,-eta10
1845 1080 FORMAT(' JET CROSS SECTION PT-DISTRIBUTION',/,
1846  & ' RAP.=',f5.2,'...',f4.2,/)
1847  CALL plot(ab,x,550,11,50,pt10,dpt1,xsmin,xsstep)
1848  ELSEIF ( iout.EQ.5 ) THEN
1849  fac = xsect(2,0)/(deta2*dpt2*mxsect(1,0))
1850  DO 80 i=1,50
1851  eta = eta20+(i-1)*deta2
1852  DO 75 j=1,5
1853  ab(i,j) = eta
1854  IF ( hep(i,j).GT.1.d-35 ) x(i,j) = log10(fac*hep(i,j))
1855 75 CONTINUE
1856 80 CONTINUE
1857  WRITE(6,1090) pt20,pt20+4.*dpt2
1858 1090 FORMAT(' JET CROSS SECTION RAP.-DISTRIBUTION',/,
1859  & ' PT=',f6.2,'...',f6.2,/)
1860  CALL plot(ab(1,1),x(1,1),250,5,50,eta20,deta2,xsmin,xsstep)
1861  ELSEIF ( iout.EQ.6 ) THEN
1862  fac = xsect(2,0)/(deta2*mxsect(1,0))
1863  DO 100 i=1,50
1864  eta = eta20+(i-1)*deta2
1865  DO 90 j=1,8
1866  ab(i,j-6) = eta
1867  IF ( hem(i,j).GT.1.d-35 ) x(i,j-6) = log10(fac*hem(i,j))
1868 90 CONTINUE
1869 100 CONTINUE
1870  WRITE(6,1100)
1871 1100 FORMAT(' JET CROSS SECTION RAP.-DISTRIBUTION',/,
1872  & ' FOR THE DIFF. SUBPROCESSES',/)
1873  CALL plot(ab,x,400,8,50,eta20,deta2,xsmin,xsstep)
1874  ENDIF
1875  WRITE(6,1110)
1876 1110 FORMAT(/)
1877  RETURN
1878  END
1879 
1880 C______________________________________________________________________
1881  SUBROUTINE hisini
1882  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1883  SAVE
1884  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
1885 C LENGTH OF HISTO : 15000 REAL*4
1886  COMMON /histo / pt10,dpt1,eta10,deta1,pt20,dpt2,eta20,deta2,
1887  & x(50,-5:5),ab(50,-5:5),hpe(50,-5:5),hep(50,5),
1888  & hpm(50,8),hem(50,8),hp(50),he(50),
1889  & fill(12192)
1890 C
1891 C INITIALIZE HISTOGRAMS
1892  dpt1 = 1.
1893  deta1 = 1.
1894  dpt2 = 2.
1895  deta2 = 0.2
1896  pt10 = 1.
1897  eta10 =-5.*deta1
1898  pt20 = 2.
1899  eta20 =-25.*deta2
1900  DO 40 i=1,50
1901  hp(i) = 0.0
1902  he(i) = 0.0
1903  DO 10 j=-5,5
1904 10 hpe(i,j) = 0.0
1905  DO 20 j=1,5
1906 20 hep(i,j) = 0.0
1907  DO 30 j=1,8
1908  hpm(i,j) = 0.0
1909 30 hem(i,j) = 0.0
1910 40 CONTINUE
1911  RETURN
1912  END
1913 C______________________________________________________________________
1914  SUBROUTINE hisfil
1915  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1916  SAVE
1917  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
1918  COMMON /haevtr/ line,lin,lrec1(mline),lrec2(mline),prec(0:3,mline)
1919  COMMON /harslt/ lscahd,lsc1hd,
1920  & etahd(mscahd,2) ,pthd(mscahd),
1921  & xhd(mscahd,2) ,vhd(mscahd) ,x0hd(mscahd,2),
1922  & ninhd(mscahd,2) ,nouthd(mscahd,2),
1923  & n0inhd(mscahd,2),nbrahd(mscahd,2),nprohd(mscahd)
1924 C LENGTH OF HISTO : 15000 REAL*4
1925  COMMON /histo / pt10,dpt1,eta10,deta1,pt20,dpt2,eta20,deta2,
1926  & x(50,-5:5),ab(50,-5:5),hpe(50,-5:5),hep(50,5),
1927  & hpm(50,8),hem(50,8),hp(50),he(50),
1928  & fill(12192)
1929 C
1930 C fill histogram
1931  DO 20 n=1,lscahd
1932  mspr = nprohd(n)
1933  DO 10 k=1,2
1934  ipt1 = int((pthd(n)-pt10)/dpt1)+1
1935  ieta1 = int((etahd(n,k)-eta10)/deta1+0.5)-5
1936  ipt2 = int((pthd(n)-pt20)/dpt2)+1
1937  ieta2 = int((etahd(n,k)-eta20)/deta2+0.5)
1938  IF ( ipt1.GE. 1 .AND. ipt1.LE.50 ) THEN
1939  hpm(ipt1,mspr) = hpm(ipt1,mspr)+1.
1940  hp(ipt1) = hp(ipt1)+1.
1941  IF ( abs(ieta1).LE.5 ) hpe(ipt1,ieta1) = hpe(ipt1,ieta1)+1.
1942  ENDIF
1943  IF ( ieta2.GE. 1 .AND. ieta2.LE.50 ) THEN
1944  hem(ieta2,mspr) = hem(ieta2,mspr)+1.
1945  he(ieta2) = he(ieta2)+1.
1946  IF ( ipt2.GE.1 .AND. ipt2.LE.5 ) hep(ieta2,ipt2) =
1947  & hep(ieta2,ipt2)+1.
1948  ENDIF
1949 10 CONTINUE
1950 20 CONTINUE
1951  RETURN
1952  END
1953 C_______________________________________________________________________
1954  SUBROUTINE hatest(IOUT)
1955  IMPLICIT DOUBLE PRECISION(a-h,o-z)
1956  SAVE
1957  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
1958  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
1959  CHARACTER*18 proc
1960  CHARACTER*11 pdset,partic
1961  COMMON /peproc/ proc(0:maxpro),pdset(23),partic(-1:1)
1962 C LENGTH OF HISTO : 15000 REAL*4
1963  COMMON /histo / vvv(50),xs(50,6),ab(50,6),dsig(0:maxpro),pd(-6:6),
1964  & fill(14328)
1965  IF ( iout.EQ.1 ) THEN
1966  CALL csharm(dsig)
1967  WRITE(6,1010) ecm,ptini(1),(proc(m),dsig(m),m=0,maxpro)
1968 1010 FORMAT(' HARD CROSS SECTIONS FOR SINGLE PROCESSES',/,
1969  & ' AT CM-ENERGY=',e8.1,' AND PTMIN=',f5.1,/,9(a25,e14.6,/))
1970  ELSEIF ( iout.EQ.2 ) THEN
1971 C PLOT PARTON DISTRIBUTIONS
1972  pdmin = 0.0
1973  pdstep = 0.04
1974  ymax = 5.0
1975  dy = ymax/50.
1976  qq = 10.0
1977  DO 10 j=1,50
1978  y = float(j-1)*dy
1979  vvv(j) = 10.0**(-y)
1980  DO 10 i=1,5
1981  xs(j,i) =-1.e+30
1982  ab(j,i) = y
1983 10 CONTINUE
1984  qq = 1.0
1985  DO 20 i=1,5
1986  qq = qq*10.
1987  DO 20 j=1,50
1988  CALL jtpdis(vvv(j),qq,1,1,pd)
1989  IF ( pd(0).GT.1.d-30 ) xs(j,i) = log10(pd(0))
1990 20 CONTINUE
1991  WRITE(6,1020)
1992 1020 FORMAT(' GLUONDISTRIBUTION OVER LOG10(X) ( Q**2=10**I;',
1993  & ' I=1...5 )')
1994  CALL plot(ab,xs,250,5,50,ymax,-dy,pdmin,pdstep)
1995  ELSEIF ( iout.EQ.3 ) THEN
1996  qqmin = 1.0
1997  qqstep = 0.1
1998  DO 30 i=1,50
1999  b = float(i-1)*qqstep+qqmin
2000  vvv(i) = 10.0**b
2001  DO 30 j=1,5
2002  xs(i,j) = -1.e+30
2003  ab(i,j) = b
2004 30 CONTINUE
2005  x = 1.0
2006  DO 40 i=1,50
2007  x = 1.0
2008  DO 40 j=1,4
2009  x = x*0.1
2010  CALL jtpdis(x,vvv(i),1,1,pd)
2011  IF ( pd(0).GT.1.d-30 ) xs(i,j) = log10(pd(0))
2012 40 CONTINUE
2013  WRITE(6,1030)
2014 1030 FORMAT(' GLUONDISTRIBUTION OVER LOG10(Q**2) ( X=10**(-I)'
2015  & ,'; I=1...4')
2016  CALL plot(ab,xs,200,4,50,qqmin,qqstep,pdmin,pdstep)
2017  ELSEIF ( iout.EQ.4 ) THEN
2018 C PLOT DIFFERENTIAL ONE JET CROSS SECTION
2019  xsmin =-6.
2020  xsstep = 0.1
2021  ptmin = 1.0
2022  ptstep = 1.0
2023  DO 50 i=1,50
2024  pt = (i-1)*ptstep+ptmin
2025  xs(i,1) =-35.0
2026  DO 50 j=1,6
2027 50 ab(i,j) = pt
2028 C DO 60 J=1,6
2029 C ETAC = (J-1)*0.5
2030  etac = 0.0
2031  DO 60 i=1,50
2032  pt = ab(i,1)
2033  CALL csj1m(pt,etac,dsig)
2034  IF ( dsig(0).GT.1.d-30 ) xs(i,1) = log10(dsig(0))
2035 60 CONTINUE
2036  WRITE(6,1040)
2037 1040 FORMAT(' DIFFERENTIAL HARD CROSS SECTION OVER PT , RAP.=0.')
2038  CALL plot(ab,xs,50,1,50,ptmin,ptstep,xsmin,xsstep)
2039  ENDIF
2040  RETURN
2041  END
2042 C-----------------------------------------------------------------------
2043 C
2044 C JTINIT.FOR
2045 C
2046 C---------------------------------------------------------------------
2047 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2048 C
2049 C PROCEDURES TO INITIALIZE PROGRAM
2050 C ( NONE OF THIS PROCEDURES IS USED DURING EVENT SIMULATION,
2051 C BEFORE DEMANDING EVENT SIMULATION SUBROUTINES
2052 C HASTRT, HARINI AND HISINI MUST BE CALLED )
2053 C
2054 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2055 C_______________________________________________________________________
2056  SUBROUTINE harini
2057  IMPLICIT DOUBLE PRECISION(a-h,o-z)
2058  SAVE
2059  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
2060  COMMON /hacons/ pi,pi2,pi4,gevtmb
2061  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
2062  COMMON /hapdco/ npdcor
2063  COMMON /haqqap/ aqqal,aqqpd,nqqal,nqqpd
2064  COMMON /haoutl/ noutl,nouter,noutco
2065  INTEGER mxsect
2066  COMMON /haxsec/ xsecta(2,-1:maxpro,4),xsect(5,-1:maxpro),
2067  & mxsect(0:2,-1:maxpro)
2068  CHARACTER*18 proc
2069  CHARACTER*11 pdset,partic
2070  COMMON /peproc/ proc(0:maxpro),pdset(23),partic(-1:1)
2071  dimension dsig(0:maxpro),alam(23),q0s(23)
2072  DATA alam / 0.20, 0.29, 0.107, 0.250, 0.178, 0.25,
2073  * 0.10, 0.19, 0.190, 0.190, 0.190, 0.19,
2074  * 0.215,0.215,0.215,
2075  * 0.231,0.231,0.322, 0.247, 0.168,0.2,0.2,0.202 /
2076  DATA q0s / 5.0 , 5.0 , 5.0 , 5.0 , 5.0 , 0.2,
2077  * 5.0 , 5.0 , 5.0 , 5.0 , 5.0 , 5.0,
2078  * 5.0 , 5.0 , 5.0 , 4.0 , 4.0 , 4.0,
2079  * 4.0 , 4.0 , 0.4 ,0.4 ,1.60 /
2080 C
2081  IF ( noutl.GE.1 )CALL timdat
2082  alasqr = alam(npd)**2
2083  q0sqr = q0s(npd)
2084  bqcd = pi4/(11.-(2./3.)*nf)
2085 C check and sort the min. pt values in PTINI(1..4)
2086  ini = 0
2087  DO 30 i=1,4
2088  IF ( ptini(i).LE..5d0.OR.ptini(i).GE.ecm*.5d0)ptini(i)=1.d+30
2089  IF ( ptini(i).NE.1.d+30 ) ini = ini+1
2090 30 CONTINUE
2091  DO 50 i=1,3
2092  DO 40 j=i+1,4
2093  IF ( ptini(j).LT.ptini(i) ) THEN
2094  ttt = ptini(j)
2095  ptini(j) = ptini(i)
2096  ptini(i) = ttt
2097  ENDIF
2098 40 CONTINUE
2099 50 CONTINUE
2100 C calculate constant weights
2101  DO 10 m=-1,maxpro
2102  xsect(3,m) = 0.0
2103  xsect(4,m) = 0.0
2104  mxsect(1,m) = 0
2105  mxsect(2,m) = 0
2106 10 CONTINUE
2107  DO 20 i = 1,4
2108  DO 20 m =-1,maxpro
2109  DO 20 j = 1,2
2110  xsecta(j,m,i) = 0.0
2111 20 CONTINUE
2112  DO 70 i=ini,1,-1
2113  CALL habint(i)
2114  IF ( noutl.GE.10 ) WRITE(6,1060) ptini(i)
2115 1060 FORMAT(' NORMALIZATION FOR PTMIN=',f10.4,' CALCULATED')
2116  CALL hamaxi(i)
2117  IF ( noutl.GE.10 ) WRITE(6,1070) ptini(i)
2118 1070 FORMAT(' MAXIMA FOR PTMIN=',f10.4,' CALCULATED')
2119  xsecta(1,0,i) = ptini(i)
2120  DO 60 m=-1,maxpro
2121  xsecta(1,m,i) = xsect(1,m)
2122  xsecta(2,m,i) = xsect(2,m)
2123 60 CONTINUE
2124 70 CONTINUE
2125 C calculate inclusive cross-sections
2126  CALL csharm(dsig)
2127  DO 80 m=0,maxpro
2128  xsect(5,m) = dsig(m)
2129 80 CONTINUE
2130 C
2131 C print results of initialization
2132 C
2133  IF ( noutl.GE.10 ) WRITE(6,'(/,1X,70(1H*))')
2134  WRITE(6,1057) ptini(1),pdset(npd),sqrt(alasqr),q0sqr
2135 1057 FORMAT(/,
2136  & ' --- parameters of the hard scattering program ---',/,
2137  & ' MIN. PT :',f15.1,/,
2138  & ' PARTON-DISTR. :',a15,/,
2139  & ' LAMBDA :',f15.3,/,
2140  & ' Q0**2 :',f15.3,/)
2141  IF ( noutl.GE.1 ) THEN
2142  WRITE(6,1050) partic(nha),partic(nhb),ecm,ptini(1),pdset(npd),
2143  & sqrt(alasqr),q0sqr,npdcor,nqqal,aqqal,nqqpd,aqqpd
2144 1050 FORMAT(/,1x,70('*'),/,
2145  & ' HARD SCATTERING PROGRAM IS INITIALIZED FOR',/,
2146  & ' PROJECTILE :',a15,/,
2147  & ' TARGET :',a15,/,
2148  & ' CM-ENERGY :',f15.1,/,
2149  & ' MIN. PT :',f15.1,/,
2150  & ' PARTON-DISTR. :',a15,/,
2151  & ' LAMBDA :',f15.3,/,
2152  & ' Q0**2 :',f15.3,/,
2153  & ' NPDCOR :',i15,/,
2154  & ' NQQAL :',i15,/,
2155  & ' AQQAL :',f15.3,/,
2156  & ' NQQPD :',i15,/,
2157  & ' AQQPD :',f15.3,/)
2158  CALL timdat
2159  ENDIF
2160  RETURN
2161  END
2162 C_______________________________________________________________________
2163  SUBROUTINE habint(IND)
2164  IMPLICIT DOUBLE PRECISION(a-h,o-z)
2165  SAVE
2166  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
2167  parameter( mxabwt = 1000 )
2168  parameter( zero=0.d0, one=1.d0)
2169  COMMON /hacons/ pi,pi2,pi4,gevtmb
2170  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
2171  COMMON /hagaup/ ngaup1,ngaup2,ngauet,ngauin
2172  INTEGER mxsect
2173  COMMON /haxsec/ xsecta(2,-1:maxpro,4),xsect(5,-1:maxpro),
2174  & mxsect(0:2,-1:maxpro)
2175  dimension absz(mxabwt),weig(mxabwt)
2176  dimension s(-1:maxpro),s1(-1:maxpro),s2(-1:maxpro),f124(-1:maxpro)
2177  DATA f124 / 1.,0.,4.,2.,2.,2.,4.,1.,4.,4. /
2178 
2179  a = (2.*ptini(ind)/ecm)**2
2180  aln = log(a)
2181  hln = log(0.5)
2182  npoint = ngauin
2183  CALL gset(zero,one,npoint,absz,weig)
2184  DO 10 m=-1,maxpro
2185  s1(m) = 0.0
2186 10 CONTINUE
2187  DO 80 i1=1,npoint
2188  z1 = absz(i1)
2189  x1 = exp(aln*z1)
2190  DO 20 m=-1,maxpro
2191  s2(m) = 0.0
2192 20 CONTINUE
2193  DO 60 i2=1,npoint
2194  z2 = (1.-z1)*absz(i2)
2195  x2 = exp(aln*z2)
2196  faxx = a/(x1*x2)
2197  w = sqrt(1.-faxx)
2198  w1 = faxx/(1.+w)
2199  wlog = log(w1)
2200  fww = faxx*wlog/w
2201  DO 30 m=-1,maxpro
2202  s(m) = 0.0
2203 30 CONTINUE
2204  DO 40 i=1,npoint
2205  z = absz(i)
2206  va =-0.5*w1/(w1+z*w)
2207  ua =-1.-va
2208  vb =-0.5*faxx/(w1+2.*w*z)
2209  ub =-1.-vb
2210  vc =-exp(hln+z*wlog)
2211  uc =-1.-vc
2212  ve =-0.5*(1.+w)+z*w
2213  ue =-1.-ve
2214  s(1) = s(1)+(1.+w)*2.25*(va*va*(3.-ua*va-va/(ua*ua))-ua)*
2215  & weig(i)
2216  s(2) = s(2)+(vc*vc+uc*uc)*((16./27.)/uc-(4./3.)*vc)*fww*
2217  & weig(i)
2218  s(3) = s(3)+(1.+w)*(1.+ua*ua)*(1.-(4./9.)*va*va/ua)*weig(i)
2219  s(5) = s(5)+((4./9.)*(1.+ub*ub+(ub*ub+vb*vb)*vb*vb)-
2220  & (8./27.)*ua*ua*va)*weig(i)
2221  s(6) = s(6)+(4./9.)*(ue*ue+ve*ve)*faxx*weig(i)
2222  s(7) = s(7)+(1.+w)*((2./9.)*(1.+ua*ua+(1.+va*va)*va*va/
2223  & (ua*ua))-(4./27.)*va/ua)*weig(i)
2224  s(8) = s(8)+(4./9.)*(1.+ub*ub)*weig(i)
2225  s(-1) = s(-1)+(1.+vc*vc)*(vc/(uc*uc)-(4./9.))*fww*weig(i)
2226 40 CONTINUE
2227  s(4) = s(2)*(9./32.)
2228  DO 50 m=-1,maxpro
2229  s2(m) = s2(m)+s(m)*weig(i2)*w
2230 50 CONTINUE
2231 60 CONTINUE
2232  DO 70 m=-1,maxpro
2233  s1(m) = s1(m)+s2(m)*(1.-z1)*weig(i1)
2234 70 CONTINUE
2235 80 CONTINUE
2236  fff = pi*gevtmb*aln*aln/(a*ecm*ecm)
2237  DO 90 m=-1,maxpro
2238  xsect(1,m) = fff*f124(m)*s1(m)
2239 90 CONTINUE
2240  xsect(1,4) = xsect(1,4)*nf
2241  xsect(1,6) = xsect(1,6)*max(0,nf-1)
2242  RETURN
2243  END
2244 C______________________________________________________________________
2245  SUBROUTINE hamaxi(IND)
2246  IMPLICIT DOUBLE PRECISION(a-h,o-z)
2247  SAVE
2248  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
2249  parameter( nkm = 5 )
2250  parameter( tiny= 1.d-30 )
2251  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
2252  INTEGER mxsect
2253  COMMON /haxsec/ xsecta(2,-1:maxpro,4),xsect(5,-1:maxpro),
2254  & mxsect(0:2,-1:maxpro)
2255  dimension z(3),d(3),ff(nkm)
2256 
2257  DO 40 nkon=1,nkm
2258  z(1) = 0.99
2259  z(2) = 0.5
2260  z(3) = 0.0
2261  d(1) =-1.0
2262  d(2) = 2.0
2263  d(3) = 2.5
2264  it = 0
2265  CALL hafdi1(nkon,z,f2,ind)
2266 10 it = it+1
2267  fold = f2
2268  DO 30 i=1,3
2269  d(i) = d(i)/5.
2270  z(i) = z(i)+d(i)
2271  CALL hafdi1(nkon,z,f3,ind)
2272  IF ( f2.GT.f3 ) z(i) = z(i)-d(i)
2273  IF ( f2.GT.f3 ) d(i) =-d(i)
2274 20 f1 = min(f2,f3)
2275  f2 = max(f2,f3)
2276  z(i) = z(i)+d(i)
2277  CALL hafdi1(nkon,z,f3,ind)
2278  IF ( f3.GT.f2 ) goto 20
2279  zz = z(i)-d(i)
2280  z(i) = zz+0.5*d(i)*(f3-f1)/max(tiny,f2+f2-f1-f3)
2281  IF ( abs(zz-z(i)).GT.d(i)*0.1d0)CALL hafdi1(nkon,z,f1,ind)
2282  IF ( f1.LE.f2 ) z(i) = zz
2283  f2 = max(f1,f2)
2284 30 CONTINUE
2285  IF ( abs(fold-f2)/f2.GT.0.002d0.OR. it.LT.3 ) goto 10
2286  ff(nkon) = f2
2287 40 CONTINUE
2288  xsect(2,1) = ff(1)*xsect(1,1)
2289  xsect(2,2) = ff(2)*xsect(1,2)
2290  xsect(2,3) = ff(4)*xsect(1,3)
2291  xsect(2,4) = ff(1)*xsect(1,4)
2292  xsect(2,5) = ff(2)*xsect(1,5)
2293  xsect(2,6) = ff(2)*xsect(1,6)
2294  xsect(2,7) = ff(3)*xsect(1,7)
2295  xsect(2,8) = ff(5)*xsect(1,8)
2296  xsect(2,-1)= ff(4)*xsect(1,-1)
2297  RETURN
2298  END
2299 C______________________________________________________________________
2300  SUBROUTINE hafdi1(NKON,Z,FDIS,IND)
2301  IMPLICIT DOUBLE PRECISION(a-h,o-z)
2302  SAVE
2303  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
2304  parameter( nkm = 5 )
2305  parameter( tiny= 1.d-30, one=1.d0 ,tiny6=1.d-06,zero=0.d0)
2306  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
2307  COMMON /haqqap/ aqqal,aqqpd,nqqal,nqqpd
2308  dimension f(nkm),pda(-6:6),pdb(-6:6),z(3)
2309 
2310  fdis = 0.0
2311 C check input values
2312  IF ( z(1).LE.0.0d0 .OR. z(1).GE.1.0d0 ) RETURN
2313  IF ( z(2).LE.0.0d0 .OR. z(2).GE.1.0d0 ) RETURN
2314  IF ( z(3).LT.0.0d0 .OR. z(3).GT.1.0d0 ) RETURN
2315  a = (2.*ptini(ind)/ecm)**2
2316  aln = log(a)
2317  y1 = exp(aln*z(1))
2318  y2 =-(1.-y1)+2.*(1.-y1)*z(2)
2319  x1 = 0.5*(y2+sqrt(y2*y2+4.*y1))
2320  x2 = x1-y2
2321  w = sqrt(max(tiny,1.-a/y1))
2322  v =-0.5+w*(z(3)-0.5)
2323  u =-(1.+v)
2324  pt = max(ptini(ind),sqrt(u*v*y1*ecm*ecm))
2325 C set hard scale QQ for alpha and partondistr.
2326  IF ( nqqal.EQ.1 ) THEN
2327  qqal = aqqal*pt*pt
2328  ELSEIF ( nqqal.EQ.2 ) THEN
2329  qqal = aqqal*y1*ecm*ecm
2330  ELSEIF ( nqqal.EQ.3 ) THEN
2331  qqal = aqqal*y1*ecm*ecm*(u*v)**(1./3.)
2332  ELSEIF ( nqqal.EQ.4 ) THEN
2333  qqal = aqqal*y1*ecm*ecm*u*v/(1.+v*v+u*u)
2334  ENDIF
2335  IF ( nqqpd.EQ.1 ) THEN
2336  qqpd = aqqpd*pt*pt
2337  ELSEIF ( nqqpd.EQ.2 ) THEN
2338  qqpd = aqqpd*y1*ecm*ecm
2339  ELSEIF ( nqqpd.EQ.3 ) THEN
2340  qqpd = aqqpd*y1*ecm*ecm*(u*v)**(1./3.)
2341  ELSEIF ( nqqpd.EQ.4 ) THEN
2342  qqpd = aqqpd*y1*ecm*ecm*u*v/(1.+v*v+u*u)
2343  ENDIF
2344  factor = (bqcd/log(max(qqal/alasqr,1.1*one)))**2
2345 C calculate partondistributions
2346  CALL jtpdis(x1,qqpd,nha,0,pda)
2347  CALL jtpdis(x2,qqpd,nhb,0,pdb)
2348 C calculate full distribution FDIS
2349  DO 10 n=1,nkm
2350  f(n) = 0.0
2351 10 CONTINUE
2352  DO 20 i=1,nf
2353  f(2) = f(2)+pda(i)*pdb(-i)+pda(-i)*pdb( i)
2354  f(3) = f(3)+pda(i)*pdb( i)+pda(-i)*pdb(-i)
2355  f(4) = f(4)+pda(i)+pda(-i)
2356  f(5) = f(5)+pdb(i)+pdb(-i)
2357 20 CONTINUE
2358  f(1) = pda(0)*pdb(0)
2359  t = pda(0)*f(5)+pdb(0)*f(4)
2360  f(5) = f(4)*f(5)-(f(2)+f(3))
2361  f(4) = t
2362  fdis = max(zero,f(nkon)*factor)
2363  RETURN
2364  END
2365 C______________________________________________________________________
2366  SUBROUTINE hastrt
2367  IMPLICIT DOUBLE PRECISION(a-h,o-z)
2368  SAVE
2369  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
2370  COMMON /hacons/ pi,pi2,pi4,gevtmb
2371  COMMON /hapara/ ecm,ptini(4),q0sqr,alasqr,bqcd,npd,nf,nha,nhb
2372  COMMON /hapadi/ npdm
2373  COMMON /hapdco/ npdcor
2374  COMMON /haqqap/ aqqal,aqqpd,nqqal,nqqpd
2375  COMMON /haoutl/ noutl,nouter,noutco
2376  COMMON /hacuts/ ptl,ptu,etacl,etacu,etadl,etadu
2377  COMMON /hagaup/ ngaup1,ngaup2,ngauet,ngauin
2378  COMMON /haevtr/ line,lin,lrec1(mline),lrec2(mline),prec(0:3,mline)
2379  INTEGER mxsect
2380  COMMON /haxsec/ xsecta(2,-1:maxpro,4),xsect(5,-1:maxpro),
2381  & mxsect(0:2,-1:maxpro)
2382  COMMON /haxsum/xshmx
2383 C initialize COMMON /HAXSUM/
2384  xshmx = 1.
2385 C initialize COMMON /HACONS/
2386  pi = 3.1415927
2387  pi2 = 2.*pi
2388  pi4 = 4.*pi
2389  gevtmb = 0.389365
2390 C initialize COMMON /HAPARA/
2391  ecm = 1800.
2392  ptini(1) = 2.0
2393  ptini(2) = 0.0
2394  ptini(3) = 0.0
2395  ptini(4) = 0.0
2396  q0sqr = 5.0
2397  alasqr = 0.107**2
2398  nf = 4
2399  bqcd = pi4/(11.0-(2./3.)*nf)
2400  npd = 3
2401  nha = 1
2402  nhb =-1
2403 C initialize COMMON /HAPADI/
2404  npdm = 0
2405 C initialize COMMON /HAPDCO/
2406  npdcor = 0
2407 C initialize COMMON /HAQQAP/
2408  nqqal = 1
2409  aqqal = 0.25
2410  nqqpd = 1
2411  aqqpd = 0.25
2412 C initialize COMMON /HAOUTL/
2413  noutl = 1
2414  nouter = 1
2415  noutco = 0
2416 C initialize COMMON /HAGAUP/
2417  ngaup1 = 8
2418  ngaup2 = 8
2419  ngauet = 8
2420  ngauin = 8
2421 C initialize COMMON /HACUTS/
2422  ptl = 0.e+00
2423  ptu = 1.e+30
2424  etacl =-1.e+30
2425  etacu = 1.e+30
2426  etadl =-1.e+30
2427  etadu = 1.e+30
2428 C initialize COMMON /HAEVTR/
2429  line = 0
2430  lin = 0
2431  DO 20 l = 1,mline
2432  lrec1(l) = 0
2433  lrec2(l) = 0
2434  DO 10 i=0,3
2435  prec(i,l) = 0.0
2436 10 CONTINUE
2437 20 CONTINUE
2438 C initialize COMMON /HAXSEC/
2439  DO 40 m=-1,maxpro
2440  DO 30 i=1,5
2441  xsect(i,m) = 0.0
2442 30 CONTINUE
2443  mxsect(1,m) = 0
2444  mxsect(2,m) = 0
2445  mxsect(0,m) = 1
2446 40 CONTINUE
2447  RETURN
2448  END
2449 C----------------------------------------------------------------------
2450  BLOCK DATA jtdata
2451  IMPLICIT DOUBLE PRECISION(a-h,o-z)
2452  SAVE
2453  parameter( maxpro = 8 , mline = 1000 , mscahd = 250 )
2454  CHARACTER*18 proc
2455  CHARACTER*11 pdset,partic
2456  COMMON /peproc/ proc(0:maxpro),pdset(23),partic(-1:1)
2457 
2458  DATA proc / 'SUM OVER PROCESSES', 'G +G --> G +G ',
2459  & 'Q +QB --> G +G ', 'G +Q --> G +Q ',
2460  & 'G +G --> Q +QB ', 'Q +QB --> Q +QB ',
2461  & 'Q +QB --> QS +QBS', 'Q +Q --> Q +Q ',
2462  & 'Q +QS --> Q +QS ' /
2463  DATA pdset / ' EHLQ SET 1',' EHLQ SET 2',' MRS SET 1',
2464  & ' MRS SET 2',' MRS SET 3',' GRV LO ',
2465  & ' HMRS SET 1',' HMRS SET 2',' KMRS SET 1',
2466  & ' KMRS SET 2',' KMRS SET 3',' KMRS SET 4',
2467  & ' MRS(S0) ',' MRS(D0) ',' MRS(D-) ',
2468  & ' CTEQ 1M ',' CTEQ 1MS ',' CTEQ 1ML ',
2469  & ' CTEQ 1D ',' CTEQ 1L ',' GRV94LO1 ' ,
2470  & ' GRV98LO ',' CTEQ96 '/
2471  DATA partic / ' ANTIPROTON',' ',' PROTON' /
2472  END
2473 C***********************************************************************
2474 C
2475 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
2476 * *
2477 * G R V - P R O T O N - P A R A M E T R I Z A T I O N S *
2478 * *
2479 * 1994 UPDATE *
2480 * *
2481 * FOR A DETAILED EXPLANATION SEE *
2482 * M. GLUECK, E.REYA, A.VOGT : *
2483 * DO-TH 94/24 = DESY 94-206 *
2484 * (TO APPEAR IN Z. PHYS. C) *
2485 * *
2486 * THE PARAMETRIZATIONS ARE FITTED TO THE EVOLVED PARTONS FOR *
2487 * Q**2 / GEV**2 BETWEEN 0.4 AND 1.E6 *
2488 * X BETWEEN 1.E-5 AND 1. *
2489 * LARGE-X REGIONS, WHERE THE DISTRIBUTION UNDER CONSIDERATION *
2490 * IS NEGLIGIBLY SMALL, WERE EXCLUDED FROM THE FIT. *
2491 * *
2492 * HEAVY QUARK THRESHOLDS Q(H) = M(H) IN THE BETA FUNCTION : *
2493 * M(C) = 1.5, M(B) = 4.5 *
2494 * CORRESPONDING LAMBDA(F) VALUES IN GEV FOR Q**2 > M(H)**2 : *
2495 * LO : LAMBDA(3) = 0.232, LAMBDA(4) = 0.200, *
2496 * LAMBDA(5) = 0.153, *
2497 * NLO : LAMBDA(3) = 0.248, LAMBDA(4) = 0.200, *
2498 * LAMBDA(5) = 0.131. *
2499 * THE NUMBER OF ACTIVE QUARK FLAVOURS IS NF = 3 EVERYWHERE *
2500 * EXCEPT IN THE BETA FUNCTION, I.E. THE HEAVY QUARKS C,B,... *
2501 * ARE NOT PRESENT AS PARTONS IN THE Q2-EVOLUTION. *
2502 * IF NEEDED, HEAVY QUARK DENSITIES CAN BE TAKEN FROM THE 1991 *
2503 * GRV PARAMETRIZATION. *
2504 * *
2505 * NLO DISTRIBUTIONS ARE GIVEN IN MS-BAR FACTORIZATION SCHEME *
2506 * (SUBROUTINE GRV94HO) AS WELL AS IN THE DIS SCHEME (GRV94DI), *
2507 * THE LEADING ORDER PARAMETRIZATION IS PROVIDED BY "GRV94LO". *
2508 * *
2509 * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
2510 *
2511 *...INPUT PARAMETERS :
2512 *
2513 * X = MOMENTUM FRACTION
2514 * Q2 = SCALE Q**2 IN GEV**2
2515 *
2516 *...OUTPUT (ALWAYS X TIMES THE DISTRIBUTION) :
2517 *
2518 * UV = U(VAL) = U - U(BAR)
2519 * DV = D(VAL) = D - D(BAR)
2520 * DEL = D(BAR) - U(BAR)
2521 * UDB = U(BAR) + D(BAR)
2522 * SB = S = S(BAR)
2523 * GL = GLUON
2524 *
2525 *...LO PARAMETRIZATION :
2526 *
2527  SUBROUTINE dor94lo (X, Q2, UV, DV, DEL, UDB, SB, GL)
2528  IMPLICIT DOUBLE PRECISION (a - z)
2529  SAVE
2530  mu2 = 0.23
2531  lam2 = 0.2322 * 0.2322
2532  s = log(log(q2/lam2) / log(mu2/lam2))
2533  ds = sqrt(s)
2534  s2 = s * s
2535  s3 = s2 * s
2536 *...UV :
2537  nu = 2.284 + 0.802 * s + 0.055 * s2
2538  aku = 0.590 - 0.024 * s
2539  bku = 0.131 + 0.063 * s
2540  au = -0.449 - 0.138 * s - 0.076 * s2
2541  bu = 0.213 + 2.669 * s - 0.728 * s2
2542  cu = 8.854 - 9.135 * s + 1.979 * s2
2543  du = 2.997 + 0.753 * s - 0.076 * s2
2544  uv = dor94fv(x, nu, aku, bku, au, bu, cu, du)
2545 *...DV :
2546  nd = 0.371 + 0.083 * s + 0.039 * s2
2547  akd = 0.376
2548  bkd = 0.486 + 0.062 * s
2549  ad = -0.509 + 3.310 * s - 1.248 * s2
2550  bd = 12.41 - 10.52 * s + 2.267 * s2
2551  cd = 6.373 - 6.208 * s + 1.418 * s2
2552  dd = 3.691 + 0.799 * s - 0.071 * s2
2553  dv = dor94fv(x, nd, akd, bkd, ad, bd, cd, dd)
2554 *...DEL :
2555  ne = 0.082 + 0.014 * s + 0.008 * s2
2556  ake = 0.409 - 0.005 * s
2557  bke = 0.799 + 0.071 * s
2558  ae = -38.07 + 36.13 * s - 0.656 * s2
2559  be = 90.31 - 74.15 * s + 7.645 * s2
2560  ce = 0.0
2561  de = 7.486 + 1.217 * s - 0.159 * s2
2562  del = dor94fv(x, ne, ake, bke, ae, be, ce, de)
2563 *...UDB :
2564  alx = 1.451
2565  bex = 0.271
2566  akx = 0.410 - 0.232 * s
2567  bkx = 0.534 - 0.457 * s
2568  agx = 0.890 - 0.140 * s
2569  bgx = -0.981
2570  cx = 0.320 + 0.683 * s
2571  dx = 4.752 + 1.164 * s + 0.286 * s2
2572  ex = 4.119 + 1.713 * s
2573  esx = 0.682 + 2.978 * s
2574  udb=dor94fw(x, s, alx, bex, akx, bkx, agx, bgx, cx, dx, ex, esx)
2575 *...SB :
2576  als = 0.914
2577  bes = 0.577
2578  aks = 1.798 - 0.596 * s
2579  as = -5.548 + 3.669 * ds - 0.616 * s
2580  bs = 18.92 - 16.73 * ds + 5.168 * s
2581  dst = 6.379 - 0.350 * s + 0.142 * s2
2582  est = 3.981 + 1.638 * s
2583  ess = 6.402
2584  sb = dor94fs(x, s, als, bes, aks, as, bs, dst, est, ess)
2585 *...GL :
2586  alg = 0.524
2587  beg = 1.088
2588  akg = 1.742 - 0.930 * s
2589  bkg = - 0.399 * s2
2590  ag = 7.486 - 2.185 * s
2591  bg = 16.69 - 22.74 * s + 5.779 * s2
2592  cg = -25.59 + 29.71 * s - 7.296 * s2
2593  dg = 2.792 + 2.215 * s + 0.422 * s2 - 0.104 * s3
2594  eg = 0.807 + 2.005 * s
2595  esg = 3.841 + 0.316 * s
2596  gl =dor94fw(x, s, alg, beg, akg, bkg, ag, bg, cg, dg, eg, esg)
2597  RETURN
2598  END
2599 *
2600 *...NLO PARAMETRIZATION (MS(BAR)) :
2601 *
2602  SUBROUTINE dor94ho (X, Q2, UV, DV, DEL, UDB, SB, GL)
2603  IMPLICIT DOUBLE PRECISION (a - z)
2604  SAVE
2605  mu2 = 0.34
2606  lam2 = 0.248 * 0.248
2607  s = log(log(q2/lam2) / log(mu2/lam2))
2608  ds = sqrt(s)
2609  s2 = s * s
2610  s3 = s2 * s
2611 *...UV :
2612  nu = 1.304 + 0.863 * s
2613  aku = 0.558 - 0.020 * s
2614  bku = 0.183 * s
2615  au = -0.113 + 0.283 * s - 0.321 * s2
2616  bu = 6.843 - 5.089 * s + 2.647 * s2 - 0.527 * s3
2617  cu = 7.771 - 10.09 * s + 2.630 * s2
2618  du = 3.315 + 1.145 * s - 0.583 * s2 + 0.154 * s3
2619  uv = dor94fv(x, nu, aku, bku, au, bu, cu, du)
2620 *...DV :
2621  nd = 0.102 - 0.017 * s + 0.005 * s2
2622  akd = 0.270 - 0.019 * s
2623  bkd = 0.260
2624  ad = 2.393 + 6.228 * s - 0.881 * s2
2625  bd = 46.06 + 4.673 * s - 14.98 * s2 + 1.331 * s3
2626  cd = 17.83 - 53.47 * s + 21.24 * s2
2627  dd = 4.081 + 0.976 * s - 0.485 * s2 + 0.152 * s3
2628  dv = dor94fv(x, nd, akd, bkd, ad, bd, cd, dd)
2629 *...DEL :
2630  ne = 0.070 + 0.042 * s - 0.011 * s2 + 0.004 * s3
2631  ake = 0.409 - 0.007 * s
2632  bke = 0.782 + 0.082 * s
2633  ae = -29.65 + 26.49 * s + 5.429 * s2
2634  be = 90.20 - 74.97 * s + 4.526 * s2
2635  ce = 0.0
2636  de = 8.122 + 2.120 * s - 1.088 * s2 + 0.231 * s3
2637  del = dor94fv(x, ne, ake, bke, ae, be, ce, de)
2638 *...UDB :
2639  alx = 0.877
2640  bex = 0.561
2641  akx = 0.275
2642  bkx = 0.0
2643  agx = 0.997
2644  bgx = 3.210 - 1.866 * s
2645  cx = 7.300
2646  dx = 9.010 + 0.896 * ds + 0.222 * s2
2647  ex = 3.077 + 1.446 * s
2648  esx = 3.173 - 2.445 * ds + 2.207 * s
2649  udb=dor94fw(x, s, alx, bex, akx, bkx, agx, bgx, cx, dx, ex, esx)
2650 *...SB :
2651  als = 0.756
2652  bes = 0.216
2653  aks = 1.690 + 0.650 * ds - 0.922 * s
2654  as = -4.329 + 1.131 * s
2655  bs = 9.568 - 1.744 * s
2656  dst = 9.377 + 1.088 * ds - 1.320 * s + 0.130 * s2
2657  est = 3.031 + 1.639 * s
2658  ess = 5.837 + 0.815 * s
2659  sb = dor94fs(x, s, als, bes, aks, as, bs, dst, est, ess)
2660 *...GL :
2661  alg = 1.014
2662  beg = 1.738
2663  akg = 1.724 + 0.157 * s
2664  bkg = 0.800 + 1.016 * s
2665  ag = 7.517 - 2.547 * s
2666  bg = 34.09 - 52.21 * ds + 17.47 * s
2667  cg = 4.039 + 1.491 * s
2668  dg = 3.404 + 0.830 * s
2669  eg = -1.112 + 3.438 * s - 0.302 * s2
2670  esg = 3.256 - 0.436 * s
2671  gl =dor94fw(x, s, alg, beg, akg, bkg, ag, bg, cg, dg, eg, esg)
2672  RETURN
2673  END
2674 *
2675 *...NLO PARAMETRIZATION (DIS) :
2676 *
2677  SUBROUTINE dor94di (X, Q2, UV, DV, DEL, UDB, SB, GL)
2678  IMPLICIT DOUBLE PRECISION (a - z)
2679  SAVE
2680  mu2 = 0.34
2681  lam2 = 0.248 * 0.248
2682  s = log(log(q2/lam2) / log(mu2/lam2))
2683  ds = sqrt(s)
2684  s2 = s * s
2685  s3 = s2 * s
2686 *...UV :
2687  nu = 2.484 + 0.116 * s + 0.093 * s2
2688  aku = 0.563 - 0.025 * s
2689  bku = 0.054 + 0.154 * s
2690  au = -0.326 - 0.058 * s - 0.135 * s2
2691  bu = -3.322 + 8.259 * s - 3.119 * s2 + 0.291 * s3
2692  cu = 11.52 - 12.99 * s + 3.161 * s2
2693  du = 2.808 + 1.400 * s - 0.557 * s2 + 0.119 * s3
2694  uv = dor94fv(x, nu, aku, bku, au, bu, cu, du)
2695 *...DV :
2696  nd = 0.156 - 0.017 * s
2697  akd = 0.299 - 0.022 * s
2698  bkd = 0.259 - 0.015 * s
2699  ad = 3.445 + 1.278 * s + 0.326 * s2
2700  bd = -6.934 + 37.45 * s - 18.95 * s2 + 1.463 * s3
2701  cd = 55.45 - 69.92 * s + 20.78 * s2
2702  dd = 3.577 + 1.441 * s - 0.683 * s2 + 0.179 * s3
2703  dv = dor94fv(x, nd, akd, bkd, ad, bd, cd, dd)
2704 *...DEL :
2705  ne = 0.099 + 0.019 * s + 0.002 * s2
2706  ake = 0.419 - 0.013 * s
2707  bke = 1.064 - 0.038 * s
2708  ae = -44.00 + 98.70 * s - 14.79 * s2
2709  be = 28.59 - 40.94 * s - 13.66 * s2 + 2.523 * s3
2710  ce = 84.57 - 108.8 * s + 31.52 * s2
2711  de = 7.469 + 2.480 * s - 0.866 * s2
2712  del = dor94fv(x, ne, ake, bke, ae, be, ce, de)
2713 *...UDB :
2714  alx = 1.215
2715  bex = 0.466
2716  akx = 0.326 + 0.150 * s
2717  bkx = 0.956 + 0.405 * s
2718  agx = 0.272
2719  bgx = 3.794 - 2.359 * ds
2720  cx = 2.014
2721  dx = 7.941 + 0.534 * ds - 0.940 * s + 0.410 * s2
2722  ex = 3.049 + 1.597 * s
2723  esx = 4.396 - 4.594 * ds + 3.268 * s
2724  udb=dor94fw(x, s, alx, bex, akx, bkx, agx, bgx, cx, dx, ex, esx)
2725 *...SB :
2726  als = 0.175
2727  bes = 0.344
2728  aks = 1.415 - 0.641 * ds
2729  as = 0.580 - 9.763 * ds + 6.795 * s - 0.558 * s2
2730  bs = 5.617 + 5.709 * ds - 3.972 * s
2731  dst = 13.78 - 9.581 * s + 5.370 * s2 - 0.996 * s3
2732  est = 4.546 + 0.372 * s2
2733  ess = 5.053 - 1.070 * s + 0.805 * s2
2734  sb = dor94fs(x, s, als, bes, aks, as, bs, dst, est, ess)
2735 *...GL :
2736  alg = 1.258
2737  beg = 1.846
2738  akg = 2.423
2739  bkg = 2.427 + 1.311 * s - 0.153 * s2
2740  ag = 25.09 - 7.935 * s
2741  bg = -14.84 - 124.3 * ds + 72.18 * s
2742  cg = 590.3 - 173.8 * s
2743  dg = 5.196 + 1.857 * s
2744  eg = -1.648 + 3.988 * s - 0.432 * s2
2745  esg = 3.232 - 0.542 * s
2746  gl = dor94fw(x, s, alg, beg, akg, bkg, ag, bg, cg, dg, eg, esg)
2747  RETURN
2748  END
2749 *
2750 *...FUNCTIONAL FORMS OF THE PARAMETRIZATIONS :
2751 *
2752  FUNCTION dor94fv (X, N, AK, BK, A, B, C, D)
2753  IMPLICIT DOUBLE PRECISION (a - z)
2754  SAVE
2755  dx = sqrt(x)
2756  dor94fv=n * x**ak * (1.+ a*x**bk + x * (b + c*dx)) * (1.- x)**d
2757  RETURN
2758  END
2759 *
2760  FUNCTION dor94fw (X, S, AL, BE, AK, BK, A, B, C, D, E, ES)
2761  IMPLICIT DOUBLE PRECISION (a - z)
2762  SAVE
2763  lx = log(1./x)
2764  dor94fw = (x**ak * (a + x * (b + x*c)) * lx**bk + s**al
2765  1 * dexp(-e + sqrt(es * s**be * lx))) * (1.- x)**d
2766  RETURN
2767  END
2768 *
2769  FUNCTION dor94fs (X, S, AL, BE, AK, AG, B, D, E, ES)
2770  IMPLICIT DOUBLE PRECISION (a - z)
2771  SAVE
2772  dx = sqrt(x)
2773  lx = log(1./x)
2774  dor94fs = s**al / lx**ak * (1.+ ag*dx + b*x) * (1.- x)**d
2775  1 * dexp(-e + sqrt(es * s**be * lx))
2776  RETURN
2777  END
2778 *
2779 
2780 C---------------- end of file -----------------------------------------
Generated on Sat May 25 2013 14:32:39 for Geant4 by   doxygen 1.8.4