JHUGen MELA  v2.4.1
Matrix element calculations as used in JHUGen. MELA is an important tool that was used for the Higgs boson discovery and for precise measurements of its structure and interactions. Please see the website https://spin.pha.jhu.edu/ and papers cited there for more details, and kindly cite those papers when using this code.
Functions/Subroutines
Cteq61Pdf.f File Reference

Go to the source code of this file.

Functions/Subroutines

function ctq6pdf (Iparton, X, Q)
 
subroutine setctq6 (Iset)
 
subroutine readtbl (Nu)
 
function nextun ()
 
subroutine polint (XA, YA, N, X, Y, DY)
 
function partonx6 (IPRTN, XX, QQ)
 

Function/Subroutine Documentation

◆ ctq6pdf()

function ctq6pdf (   Iparton,
  X,
  Q 
)

Definition at line 93 of file Cteq61Pdf.f.

93  Implicit Double Precision (a-h,o-z)
94  Logical Warn
95  Common
96  > / ctqpar2 / nx, nt, nfmx
97  > / qcdtable / alambda, nfl, iorder
98 
99  Data warn /.true./
100  save warn
101 
102  If (x .lt. 0d0 .or. x .gt. 1d0) Then
103  print *, 'X out of range in Ctq6Pdf: ', x
104  stop
105  Endif
106  If (q .lt. alambda) Then
107  print *, 'Q out of range in Ctq6Pdf: ', q
108  stop
109  Endif
110  If ((iparton .lt. -nfmx .or. iparton .gt. nfmx)) Then
111  If (warn) Then
112 C put a warning for calling extra flavor.
113  print *, 'Warning: Iparton out of range in Ctq6Pdf! '
114  print *, 'Iparton, MxFlvN0: ', iparton, nfmx
115  Endif
116  ctq6pdf = 0d0
117  Return
118  Endif
119 
120  ctq6pdf = partonx6(iparton, x, q)
121  if(ctq6pdf.lt.0.d0) ctq6pdf = 0.d0
122  if(isnan(ctq6pdf)) then
123  if (warn) then
124  print *, 'Warning: Ctq6Pdf is NaN!'
125  print *, 'Iparton, MxFlvN0, X, Q: ', iparton, nfmx, x, q
126  endif
127  ctq6pdf = 0.d0
128  endif
129 
130  if (warn) warn = .false.
131 
132  Return
133 
134 C ********************

◆ nextun()

function nextun

Definition at line 224 of file Cteq61Pdf.f.

224 C Returns an unallocated FORTRAN i/o unit.
225  Logical EX
226 C
227  Do 10 n = 10, 300
228  INQUIRE (unit=n, opened=ex)
229  If (.NOT. ex) then
230  nextun = n
231  Return
232  Endif
233  10 Continue
234  stop ' There is no available I/O unit. '
235 C *************************

◆ partonx6()

function partonx6 (   IPRTN,
  XX,
  QQ 
)

Definition at line 281 of file Cteq61Pdf.f.

281 
282 c Given the parton distribution function in the array U in
283 c COMMON / PEVLDT / , this routine interpolates to find
284 c the parton distribution at an arbitray point in x and q.
285 c
286  Implicit Double Precision (a-h,o-z)
287 
288  parameter(mxx = 96, mxq = 20, mxf = 5)
289  parameter(mxqx= mxq * mxx, mxpqx = mxqx * (mxf+3))
290 
291  Common
292  > / ctqpar1 / al, xv(0:mxx), tv(0:mxq), upd(mxpqx)
293  > / ctqpar2 / nx, nt, nfmx
294  > / xqrange / qini, qmax, xmin
295 
296  dimension fvec(4), fij(4)
297  dimension xvpow(0:mxx)
298  Data onep / 1.00001 /
299  Data xpow / 0.3d0 / !**** choice of interpolation variable
300  Data nqvec / 4 /
301  Data ientry / 0 /
302  Save ientry,xvpow
303 
304 c store the powers used for interpolation on first call...
305  if(ientry .eq. 0) then
306  ientry = 1
307 
308  xvpow(0) = 0d0
309  do i = 1, nx
310  xvpow(i) = xv(i)**xpow
311  enddo
312  endif
313 
314  x = xx
315  q = qq
316  tt = log(log(q/al))
317 
318 c ------------- find lower end of interval containing x, i.e.,
319 c get jx such that xv(jx) .le. x .le. xv(jx+1)...
320  jlx = -1
321  ju = nx+1
322  11 If (ju-jlx .GT. 1) Then
323  jm = (ju+jlx) / 2
324  If (x .Ge. xv(jm)) Then
325  jlx = jm
326  Else
327  ju = jm
328  Endif
329  Goto 11
330  Endif
331 C Ix 0 1 2 Jx JLx Nx-2 Nx
332 C |---|---|---|...|---|-x-|---|...|---|---|
333 C x 0 Xmin x 1
334 C
335  If (jlx .LE. -1) Then
336  print '(A,1pE12.4)', 'Severe error: x <= 0 in PartonX6! x = ', x
337  stop
338  ElseIf (jlx .Eq. 0) Then
339  jx = 0
340  Elseif (jlx .LE. nx-2) Then
341 
342 C For interrior points, keep x in the middle, as shown above
343  jx = jlx - 1
344  Elseif (jlx.Eq.nx-1 .or. x.LT.onep) Then
345 
346 C We tolerate a slight over-shoot of one (OneP=1.00001),
347 C perhaps due to roundoff or whatever, but not more than that.
348 C Keep at least 4 points >= Jx
349  jx = jlx - 2
350  Else
351  print '(A,1pE12.4)', 'Severe error: x > 1 in PartonX6! x = ', x
352  stop
353  Endif
354 C ---------- Note: JLx uniquely identifies the x-bin; Jx does not.
355 
356 C This is the variable to be interpolated in
357  ss = x**xpow
358 
359  If (jlx.Ge.2 .and. jlx.Le.nx-2) Then
360 
361 c initiation work for "interior bins": store the lattice points in s...
362  svec1 = xvpow(jx)
363  svec2 = xvpow(jx+1)
364  svec3 = xvpow(jx+2)
365  svec4 = xvpow(jx+3)
366 
367  s12 = svec1 - svec2
368  s13 = svec1 - svec3
369  s23 = svec2 - svec3
370  s24 = svec2 - svec4
371  s34 = svec3 - svec4
372 
373  sy2 = ss - svec2
374  sy3 = ss - svec3
375 
376 c constants needed for interpolating in s at fixed t lattice points...
377  const1 = s13/s23
378  const2 = s12/s23
379  const3 = s34/s23
380  const4 = s24/s23
381  s1213 = s12 + s13
382  s2434 = s24 + s34
383  sdet = s12*s34 - s1213*s2434
384  tmp = sy2*sy3/sdet
385  const5 = (s34*sy2-s2434*sy3)*tmp/s12
386  const6 = (s1213*sy2-s12*sy3)*tmp/s34
387 
388  EndIf
389 
390 c --------------Now find lower end of interval containing Q, i.e.,
391 c get jq such that qv(jq) .le. q .le. qv(jq+1)...
392  jlq = -1
393  ju = nt+1
394  12 If (ju-jlq .GT. 1) Then
395  jm = (ju+jlq) / 2
396  If (tt .GE. tv(jm)) Then
397  jlq = jm
398  Else
399  ju = jm
400  Endif
401  Goto 12
402  Endif
403 
404  If (jlq .LE. 0) Then
405  jq = 0
406  Elseif (jlq .LE. nt-2) Then
407 C keep q in the middle, as shown above
408  jq = jlq - 1
409  Else
410 C JLq .GE. Nt-1 case: Keep at least 4 points >= Jq.
411  jq = nt - 3
412 
413  Endif
414 C This is the interpolation variable in Q
415 
416  If (jlq.GE.1 .and. jlq.LE.nt-2) Then
417 c store the lattice points in t...
418  tvec1 = tv(jq)
419  tvec2 = tv(jq+1)
420  tvec3 = tv(jq+2)
421  tvec4 = tv(jq+3)
422 
423  t12 = tvec1 - tvec2
424  t13 = tvec1 - tvec3
425  t23 = tvec2 - tvec3
426  t24 = tvec2 - tvec4
427  t34 = tvec3 - tvec4
428 
429  ty2 = tt - tvec2
430  ty3 = tt - tvec3
431 
432  tmp1 = t12 + t13
433  tmp2 = t24 + t34
434 
435  tdet = t12*t34 - tmp1*tmp2
436 
437  EndIf
438 
439 
440 c get the pdf function values at the lattice points...
441 
442  If (iprtn .GE. 3) Then
443  ip = - iprtn
444  Else
445  ip = iprtn
446  EndIf
447  jtmp = ((ip + nfmx)*(nt+1)+(jq-1))*(nx+1)+jx+1
448 
449  Do it = 1, nqvec
450 
451  j1 = jtmp + it*(nx+1)
452 
453  If (jx .Eq. 0) Then
454 C For the first 4 x points, interpolate x^2*f(x,Q)
455 C This applies to the two lowest bins JLx = 0, 1
456 C We can not put the JLx.eq.1 bin into the "interrior" section
457 C (as we do for q), since Upd(J1) is undefined.
458  fij(1) = 0
459  fij(2) = upd(j1+1) * xv(1)**2
460  fij(3) = upd(j1+2) * xv(2)**2
461  fij(4) = upd(j1+3) * xv(3)**2
462 C
463 C Use Polint which allows x to be anywhere w.r.t. the grid
464 
465  Call polint (xvpow(0), fij(1), 4, ss, fx, dfx)
466 
467  If (x .GT. 0d0) fvec(it) = fx / x**2
468 C Pdf is undefined for x.eq.0
469  ElseIf (jlx .Eq. nx-1) Then
470 C This is the highest x bin:
471 
472  Call polint (xvpow(nx-3), upd(j1), 4, ss, fx, dfx)
473 
474  fvec(it) = fx
475 
476  Else
477 C for all interior points, use Jon's in-line function
478 C This applied to (JLx.Ge.2 .and. JLx.Le.Nx-2)
479  sf2 = upd(j1+1)
480  sf3 = upd(j1+2)
481 
482  g1 = sf2*const1 - sf3*const2
483  g4 = -sf2*const3 + sf3*const4
484 
485  fvec(it) = (const5*(upd(j1)-g1)
486  & + const6*(upd(j1+3)-g4)
487  & + sf2*sy3 - sf3*sy2) / s23
488 
489  Endif
490 
491  enddo
492 C We now have the four values Fvec(1:4)
493 c interpolate in t...
494 
495  If (jlq .LE. 0) Then
496 C 1st Q-bin, as well as extrapolation to lower Q
497  Call polint (tv(0), fvec(1), 4, tt, ff, dfq)
498 
499  ElseIf (jlq .GE. nt-1) Then
500 C Last Q-bin, as well as extrapolation to higher Q
501  Call polint (tv(nt-3), fvec(1), 4, tt, ff, dfq)
502  Else
503 C Interrior bins : (JLq.GE.1 .and. JLq.LE.Nt-2)
504 C which include JLq.Eq.1 and JLq.Eq.Nt-2, since Upd is defined for
505 C the full range QV(0:Nt) (in contrast to XV)
506  tf2 = fvec(2)
507  tf3 = fvec(3)
508 
509  g1 = ( tf2*t13 - tf3*t12) / t23
510  g4 = (-tf2*t34 + tf3*t24) / t23
511 
512  h00 = ((t34*ty2-tmp2*ty3)*(fvec(1)-g1)/t12
513  & + (tmp1*ty2-t12*ty3)*(fvec(4)-g4)/t34)
514 
515  ff = (h00*ty2*ty3/tdet + tf2*ty3 - tf3*ty2) / t23
516  EndIf
517 
518  partonx6 = ff
519 
520  Return
521 C ********************

◆ polint()

subroutine polint ( dimension(n)  XA,
dimension(n)  YA,
  N,
  X,
  Y,
  DY 
)

Definition at line 240 of file Cteq61Pdf.f.

240 
241  IMPLICIT DOUBLE PRECISION (a-h, o-z)
242 C Adapted from "Numerical Recipes"
243  parameter(nmax=10)
244  dimension xa(n),ya(n),c(nmax),d(nmax)
245  ns=1
246  dif=abs(x-xa(1))
247  DO 11 i=1,n
248  dift=abs(x-xa(i))
249  IF (dift.LT.dif) THEN
250  ns=i
251  dif=dift
252  ENDIF
253  c(i)=ya(i)
254  d(i)=ya(i)
255 11 CONTINUE
256  y=ya(ns)
257  ns=ns-1
258  DO 13 m=1,n-1
259  DO 12 i=1,n-m
260  ho=xa(i)-x
261  hp=xa(i+m)-x
262  w=c(i+1)-d(i)
263  den=ho-hp
264 ! IF(DEN.EQ.0.)PAUSE
265  den=w/den
266  d(i)=hp*den
267  c(i)=ho*den
268 12 CONTINUE
269  IF (2*ns.LT.n-m)THEN
270  dy=c(ns+1)
271  ELSE
272  dy=d(ns)
273  ns=ns-1
274  ENDIF
275  y=y+dy
276 13 CONTINUE
277  RETURN

◆ readtbl()

subroutine readtbl (   Nu)

Definition at line 179 of file Cteq61Pdf.f.

179  Implicit Double Precision (a-h,o-z)
180  Character Line*80
181  parameter(mxx = 96, mxq = 20, mxf = 5)
182  parameter(mxpqx = (mxf + 3) * mxq * mxx)
183  Common
184  > / ctqpar1 / al, xv(0:mxx), tv(0:mxq), upd(mxpqx)
185  > / ctqpar2 / nx, nt, nfmx
186  > / xqrange / qini, qmax, xmin
187  > / qcdtable / alambda, nfl, iorder
188  > / masstbl / amass(6)
189 
190  Read (nu, '(A)') line
191  Read (nu, '(A)') line
192  Read (nu, *) dr, fl, al, (amass(i),i=1,6)
193  iorder = nint(dr)
194  nfl = nint(fl)
195  alambda = al
196 
197  Read (nu, '(A)') line
198  Read (nu, *) nx, nt, nfmx
199 
200  Read (nu, '(A)') line
201  Read (nu, *) qini, qmax, (tv(i), i =0, nt)
202 
203  Read (nu, '(A)') line
204  Read (nu, *) xmin, (xv(i), i =0, nx)
205 
206  Do 11 iq = 0, nt
207  tv(iq) = log(log(tv(iq) /al))
208  11 Continue
209 C
210 C Since quark = anti-quark for nfl>2 at this stage,
211 C we Read out only the non-redundent data points
212 C No of flavors = NfMx (sea) + 1 (gluon) + 2 (valence)
213 
214  nblk = (nx+1) * (nt+1)
215  npts = nblk * (nfmx+3)
216  Read (nu, '(A)') line
217  Read (nu, *, iostat=iret) (upd(i), i=1,npts)
218 
219  Return
220 C ****************************

◆ setctq6()

subroutine setctq6 (   Iset)

Definition at line 138 of file Cteq61Pdf.f.

138  Implicit Double Precision (a-h,o-z)
139  parameter(isetmax0=5)
140  Character Flnm(Isetmax0)*6, nn*3, Tablefile*40
141  Data (flnm(i), i=1,isetmax0)
142  > / 'cteq6m', 'cteq6d', 'cteq6l', 'cteq6l','ctq61.'/
143  Data isetold, isetmin0, isetmin1, isetmax1 /-987,1,100,140/
144  Data isetmin2,isetmax2 /200,240/
145  save
146 
147 C If data file not initialized, do so.
148  If(iset.ne.isetold) then
149  iu= nextun()
150  If (iset.ge.isetmin0 .and. iset.le.3) Then
151  tablefile=flnm(iset)//'.tbl'
152  Elseif (iset.eq.4) Then
153  tablefile=flnm(iset)//'1.tbl'
154  Elseif (iset.ge.isetmin1 .and. iset.le.isetmax1) Then
155  write(nn,'(I3)') iset
156  tablefile=flnm(1)//nn//'.tbl'
157  Elseif (iset.ge.isetmin2 .and. iset.le.isetmax2) Then
158  write(nn,'(I3)') iset
159  tablefile=flnm(5)//nn(2:3)//'.tbl'
160  Else
161  print *, 'Invalid Iset number in SetCtq6 :', iset
162  stop
163  Endif
164  Open(iu, file='pdfs/'//tablefile
165  . ,status='OLD', err=100)
166  21 Call readtbl (iu)
167  Close (iu)
168  isetold=iset
169  Endif
170  Return
171 
172  100 print *, ' Data file ', tablefile, ' cannot be opened '
173  >//'in SetCtq6!!'
174  stop
175 C ********************
ctq6pdf
function ctq6pdf(Iparton, X, Q)
Definition: Cteq61Pdf.f:93
nextun
function nextun()
Definition: Cteq61Pdf.f:224
readtbl
subroutine readtbl(Nu)
Definition: Cteq61Pdf.f:179
partonx6
function partonx6(IPRTN, XX, QQ)
Definition: Cteq61Pdf.f:281
polint
subroutine polint(XA, YA, N, X, Y, DY)
Definition: Cteq61Pdf.f:240
Generated by   doxygen 1.8.17