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