VBF: libstructf_12.f

File libstructf_12.f, 23.4 KB (added by trac, 7 years ago)
Line 
1c--------couplings----second argument 1->W-, 2->w+, 3->Z, 4->gamma,
2c-----                                  5->gammaZ
3      subroutine SetCouplings()
4      implicit none
5      integer nf,i,j
6      double precision couplings(-6:6,5), couplings3(-6:6,5),sthw
7      common/coupl/couplings
8      common/coupl3/couplings3
9      common/nflav/nf
10      sthw=0.2315d0
11      do i=-6,6
12       do j=1,5
13       couplings(i,j)=0d0
14       couplings3(i,j)=0d0
15       end do
16      end do
17c---------- W- ------------------     
18      do i=1,3
19      couplings(2*i,1)=2d0
20      couplings3(2*i,1)=2d0
21      couplings(-2*i+1,1)=2d0
22      couplings3(-2*i+1,1)=2d0
23      end do
24c--------- W+ ------------------     
25      do i=-6,6
26      couplings(i,2)=couplings(-i,1)
27      couplings3(i,2)=couplings3(-i,1)
28      end do
29c---------gamma-----------------
30      do i=1,3
31      end do
32c--------Z0--------------------
33      do i=1,3
34      couplings(2*i,3)=2d0*(1d0/4d0+(1d0/2d0-4d0/3d0*sthw)**2)
35      couplings(-2*i,3)=2d0*(1d0/4d0+(1d0/2d0-4d0/3d0*sthw)**2)
36      couplings(2*i-1,3)=2d0*(1d0/4d0+(1d0/2d0-2d0/3d0*sthw)**2)
37      couplings(-2*i+1,3)=2d0*(1d0/4d0+(1d0/2d0-2d0/3d0*sthw)**2)
38
39      couplings3(2*i,3)=(1d0-8d0/3d0*sthw)
40      couplings3(-2*i,3)=(1d0-8d0/3d0*sthw)
41      couplings3(2*i-1,3)=(1d0-4d0/3d0*sthw)
42      couplings3(-2*i+1,3)=(1d0-4d0/3d0*sthw)
43      end do
44      do i=1,nf
45       do j=1,5
46       couplings(0,j)=couplings(0,j)+couplings(-i,j)+couplings(i,j)
47       couplings3(0,j)=couplings3(0,j)+couplings3(-i,j)+couplings3(i,j)
48       end do
49      end do
50      return
51      end
52
53ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
54c--------------------------------------------------------------------------
55c  PDF gives the value of the parton density, i.e. q(x) instead of x q(x), with the coupling included 
56
57      subroutine PDF(x,q,f)
58      implicit none
59      integer i,v,f3c
60      double precision x,q,f(-6:6),couplings(-6:6,5),couplings3(-6:6,5)
61      common/coupl/couplings
62      common/coupl3/couplings3
63      common/vect/v
64      common/f3call/f3c
65      call Evolvepdf(x,q,f)
66      if (f3c.eq.0) then
67      do i=-6,6
68       f(i)=f(i)*couplings(i,v)/x
69      end do
70      else if(f3c.eq.1) then
71      do i=-6,6
72       f(i)=f(i)*couplings3(i,v)/x
73      end do
74      endif
75
76      return
77      end
78
79      double precision function pdfg(z)
80      implicit none
81      double precision z,x,q,f(-6:6)
82      common/pdfpar/x,q
83      call pdf(x/z,q,f)
84      pdfg=f(0)/z
85      return
86      end
87
88      double precision function Sing(z)
89      implicit none
90      integer i,nf
91      common/nflav/nf
92      double precision z,x,Q, quark, antiq, f(-6:6)
93      common/pdfpar/x,q
94      call PDF(x/z,q,f)
95      quark=0d0
96      antiq=0d0
97      do i=1,nf
98      quark=quark+f(i)
99      antiq=antiq+f(-i)
100      end do
101      sing=(quark+antiq)/z
102      return
103      end
104
105      double precision function NSing(z)
106      implicit none
107      integer i,nf
108      common/nflav/nf
109      double precision z,x,Q, quark, antiq, f(-6:6)
110      common/pdfpar/x,q
111      call PDF(x/z,q,f)
112      quark=0d0
113      antiq=0d0
114      do i=1,nf
115      quark=quark+f(i)
116      antiq=antiq+f(-i)
117      end do
118      nsing=(quark-antiq)/z
119      return
120      end
121
122      double precision function SingReg(z)
123      implicit none
124      integer i,nf
125      common/nflav/nf
126      double precision z,x,Q, quark, antiq, f(-6:6)
127      common/pdfpar/x,q
128      quark=0d0
129      antiq=0d0
130      call PDF(x/z,q,f)
131       do i=1,nf
132      quark=quark+f(i)
133      antiq=antiq+f(-i)
134      end do
135      singreg=(quark+antiq)/z
136
137      quark=0d0
138      antiq=0d0
139      call PDF(x,q,f)
140       do i=1,nf
141      quark=quark+f(i)
142      antiq=antiq+f(-i)
143      end do
144      singreg=singreg-(quark+antiq)
145      return
146      end
147
148      double precision function NSingReg(z)
149      implicit none
150      integer i,nf
151      common/nflav/nf
152      double precision z,x,Q, quark, antiq, f(-6:6)
153      common/pdfpar/x,q
154      quark=0d0
155      antiq=0d0
156      call PDF(x/z,q,f)
157       do i=1,nf
158      quark=quark+f(i)
159      antiq=antiq+f(-i)
160      end do
161      nsingreg=(quark-antiq)/z
162
163      quark=0d0
164      antiq=0d0
165      call PDF(x,q,f)
166       do i=1,nf
167      quark=quark+f(i)
168      antiq=antiq+f(-i)
169      end do
170      nsingreg=nsingreg-(quark-antiq)
171      return
172      end
173
174
175ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
176c--------------------------------------------------------------------------
177c------------------------   FL   ------------------------------------------
178
179c-----NLO----------------------------------------------------------
180 
181      double precision function CLNLOa(z)
182      implicit none
183      double precision z
184      clnloa=8d0/3d0 *z
185      return
186      end
187
188      double precision function cLNLOga(z)
189      implicit none
190      double precision z
191      clnloga=2d0*z*(1d0-z)
192      return
193      end
194     
195c-----NNLO---------------------------------------------------------
196c------------------------------------------------------------------
197c----------SINGLET-------------------------------------------------
198
199      DOUBLE PRECISION FUNCTION CLNNLOSA(Y)
200      IMPLICIT DOUBLE PRECISION (A-Z)
201      INTEGER NF
202      COMMON/NFLAV/NF
203      DL  = LOG (Y)
204      DL1 = LOG (1.D0-Y)
205      CLNNLOSA = NF * ( (15.94D0 - 5.212D0 * Y) * (1.D0-Y)**2 * DL1
206     1    + (0.421D0 + 1.520D0 * Y) * DL**2 + 28.09D0 * (1.D0-Y) * DL
207     2         - (2.370D0/Y - 19.27D0) * (1.D0-Y)**3 )
208      RETURN
209      END
210
211      DOUBLE PRECISION FUNCTION CLNNLOGA (Y)
212      IMPLICIT DOUBLE PRECISION (A-Z)
213      DL  = LOG (Y)
214      DL1 = LOG (1.D0-Y)
215      CLNNLOGA = ( (94.74D0 - 49.20D0 * Y) * (1.D0-Y) * DL1**2
216     1 + 864.8D0 * (1.D0-Y) * DL1 + 1161.D0* Y * DL * DL1
217     2         + 60.060 * Y * DL**2 + 39.66D0 * (1.D0-Y) * DL
218     3         - 5.333D0 * (1.D0/Y - 1.D0) )
219      RETURN
220      END
221
222c-------------------------------------------------------------------
223c---------NON SINGLET-----------------------------------------------
224C-----CHARGED CURRENT--------------     
225      DOUBLE PRECISION FUNCTION CLNNLONSA(Y)
226      IMPLICIT DOUBLE PRECISION (A-Z)
227      INTEGER NF
228      COMMON/NFLAV/NF
229      DL  = LOG (Y)
230      DL1 = LOG (1.D0-Y)
231      CLNNLONSA =
232     1   - 52.27D0 + 100.8D0 * Y
233     2   + (23.29D0 * Y - 0.043D0) * DL**2 - 22.21D0 * DL
234     3   + 13.30D0 * DL1**2 - 59.12D0 * DL1 - 141.7D0 * DL * DL1
235     4   + NF * 16.D0/27.D0 *
236     5   ( 6.D0* Y*DL1 - 12.D0* Y*DL - 25.D0* Y + 6.D0)
237      RETURN
238      END
239
240
241c--------------neutral current---------------
242      DOUBLE PRECISION  FUNCTION CLNNLONSA_nc (Y)
243      IMPLICIT DOUBLE PRECISION (A-Z)
244      INTEGER NF
245      COMMON/NFLAV/NF
246      DL  = LOG (Y)
247      DL1 = LOG (1.D0-Y)
248      CLNNLONSA_NC =
249     1          - 40.41D0 + 97.48D0
250     2          + (26.56D0 * Y - 0.031D0) * DL**2 - 14.85D0 * DL
251     3          + 13.62D0 * DL1**2 - 55.79D0 * DL1 - 150.5D0 * DL * DL1
252     4 + NF * 16.D0/27.D0 * ( 6.D0* Y*DL1 - 12.D0* Y*DL - 25.* Y + 6.D0)
253      RETURN
254      END
255
256
257      DOUBLE PRECISION FUNCTION CLNNLONSC (Y)
258      IMPLICIT DOUBLE PRECISION (A-Z)
259      integer v
260      common/vect/v
261      if (v.le.2) then
262      CLNNLONSC = -0.150D0
263      else
264      CLNNLONSC=-0.164d0
265      end if
266      RETURN
267      END
268
269      double precision function intLnloA(z)
270      implicit none
271      double precision z,cLnloa,sing
272      intLnloa=cLnloA(z)*sing(z)
273      return
274      end
275
276      double precision function intLnlogA(z)
277      implicit none
278      double precision z,cLnloga, pdfg
279      intLnloga=cLnloga(z)*pdfg(z)
280      return
281      end
282
283      double precision function IntLNNLOsa(z)
284      implicit none
285      double precision z, clnnlosa, sing
286      intlnnlosa=clnnlosa(z)*sing(z)
287      return
288      end
289
290      double precision function IntLNNLOga(z)
291      implicit none
292      double precision z, clnnloga, pdfg
293      intlnnloGa=clnnloGa(z)*pdfg(z)
294      return
295      end
296     
297      double precision function IntLNNLOnsa(z)
298      implicit none
299      double precision z, clnnlonsa,clnnlonsa_nc, nsing
300      integer v
301      common/vect/v
302      if(v.le.2) then
303      intlnnlonsa=clnnlonsa(z)*nsing(z)
304      else
305      intlnnlonsa=clnnlonsa_nc(z)*nsing(z)
306      end if
307      return
308      end
309
310
311      double precision function FL(x,q,ord,v,zz)
312      implicit none
313      double precision x,q,zz, f(-6:6), quark, antiq,alphaspdf,pi,eps
314      integer ord, nf,v,i,vv,f2,f1
315      double precision intlnloa, intlnloga
316      double precision intlnnlosa, intlnnloga, intlnnlonsa
317      double precision clnnlonsc
318      common/prec/eps
319      double precision xx,qq
320      common/pdfpar/xx,qq
321      common/nflav/nf
322      common/vect/vv
323      integer f3c
324      common/f3call/f3c
325      double precision z,z0
326       z=(1d0-eps-x)*zz +x
327       z0=zz*x
328
329      f3c=0
330
331      xx=x
332      qq=q
333      vv=v
334      if(1d0-x.lt.eps) then
335      fl=0d0
336      else
337
338      quark=0d0
339      antiq=0d0
340      call pdf(x,q,f)
341      do i=1,nf
342      quark=quark+f(i)
343      antiq=antiq+f(-i)
344      end do
345
346      pi=3.1415926535
347
348      if (ord.eq.1) then
349      FL=0d0
350
351      else if (ord.eq.2) then
352      fL=x*alphaspdf(q)/(4d0*pi)*2d0*(1d0-eps-x)*
353     1 (intlnloa(z)+intlnloga(z))
354     
355      else if (ord.eq.3) then
356      fl=x*(alphaspdf(q)/(4d0*pi))**2* (
357     1    (1d0-x-eps)*(intlnnlosa(z)+intlnnloga(z)
358     2                +intlnnlonsa(z))
359     3   +(quark-antiq)*cLnnlonsc(x)
360     4  )
361
362      endif
363      endif
364     
365      return
366      end
367     
368
369
370ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
371c--------------------------------------------------------------------------
372c------------------------   F1   ------------------------------------------
373cccc need to correct for NC
374
375
376      double precision function F1(x,q,ord,v,zz)
377      implicit none
378      double precision x,q,zz, f(-6:6), quark, antiq,alphaspdf,pi,eps
379      integer ord, nf,v,i,vv
380      common/prec/eps
381      double precision xx,qq
382      common/pdfpar/xx,qq
383      common/nflav/nf
384      common/vect/vv
385      double precision f2,fl
386      integer f3c
387      common/f3call/f3c
388      double precision z, z0
389      z=(1d0-eps-x)*zz +x
390      z0=zz*x
391      f3c=0
392
393      xx=x
394      qq=q
395      vv=v
396      if(1d0-x.lt.eps) then
397      f1=0d0
398      else
399
400      quark=0d0
401      antiq=0d0
402      call pdf(x,q,f)
403      do i=1,nf
404      quark=quark+f(i)
405      antiq=antiq+f(-i)
406      end do
407
408      pi=3.1415926535
409
410      if (ord.eq.1) then
411      F1=(quark+antiq)/2d0
412
413      else if (ord.eq.2) then
414      f1=(f2(x,q,ord,v,zz)-fl(x,q,ord,v,zz))/ (2d0*x)
415
416      else if (ord.eq.3) then
417      f1=(f2(x,q,ord,v,zz)-fl(x,q,ord,v,zz))/ (2d0*x)
418      endif
419
420      endif
421      return
422      end
423 
424
425ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
426c--------------------------------------------------------------------------
427c------------------------   F2   ------------------------------------------
428
429      double precision function C2NLOa(z)
430      implicit none
431      double precision z
432      C2NLOa= 4d0/3d0*( 3d0+2d0*z-(1d0+z)*log(1-z)
433     1   -(1d0+z**2)/(1d0-z)*log(z))
434      return
435      end
436
437      double precision function C2NLOb(z)
438      implicit none
439      double precision z
440      C2NLOb= 4d0/3d0*(2d0*log(1d0-z)/(1d0-z)-3d0/2d0 /(1d0-z) )
441      return
442      end
443
444      double precision function C2NLOc(z)
445      implicit none
446      double precision z,pi
447      pi=3.1415926535
448      C2NLOc= -4d0/3d0*(pi**2/3d0+9d0/2d0 )
449      return
450      end
451
452      double precision function C2NLOg(z)
453      implicit none
454      double precision z
455c     C2NLOg=1d0/2d0*z*( (z**2+(1-z)**2)*log((1-z)/z)-1d0+8*z*(1-z) )
456      C2NLOg=1d0/2d0*( (z**2+(1-z)**2)*log((1-z)/z)-1d0+8*z*(1-z) )
457      return
458      end
459
460
461      double precision FUNCTION C2NNLOSA(Y)
462      IMPLICIT double precision (A-Z)
463      INTEGER NF
464      COMMON/NFLAV/NF
465      DL  = LOG (Y)
466      DL1 = LOG (1.d0-Y)
467      C2NNLOSA =   NF * ( 5.290d0 * (1.d0/Y-1.d0) + 4.310d0 * DL**3   
468     1   - 2.086d0 * DL**2 + 39.78d0 * DL - 0.101d0 * (1.d0-Y) * DL1**3
469     2   - (24.75d0 - 13.80d0 * Y) * DL**2 * DL1 + 30.23d0 * DL * DL1 )
470      RETURN
471      END
472
473
474      double precision FUNCTION C2NNLOGA (Y)
475      IMPLICIT double precision (A-Z)
476      DL  = LOG (Y)
477      DL1 = LOG (1.D0-Y)
478      C2NNLOGA =
479     1  ( 1.d0/Y * (11.90d0 + 1494.d0* DL1) + 5.319d0 * DL**3 
480     1       - 59.48d0 * DL**2 - 284.8d0 * DL + 392.4d0 - 1483.d0* DL1
481     2     + (6.445d0 + 209.4d0 * (1.d0-Y)) * DL1**3 - 24.00d0 * DL1**2
482     3        - 724.1d0 * DL**2 * DL1 - 871.8d0 * DL * DL1**2 )
483      RETURN
484      END
485
486      double precision FUNCTION C2NNLOGC (Y)
487      IMPLICIT double precision (A-Z)
488      C2NNLOGC=-0.28d0
489      RETURN
490      END
491c------------------NON SINGLET----------------------
492c---------------charged current-----------------------
493      double precision function C2NNLONSA(Y)
494      IMPLICIT double precision (A-Z)
495      INTEGER NF
496      COMMON/NFLAV/NF
497      DL  = LOG (Y)
498      DL1 = LOG (1.D0-Y)
499      C2NNLONSA = - 84.18d0 - 1010.d0* Y
500     2 -3.748d0 * DL**3 - 19.56d0 * DL**2 - 1.235d0 * DL
501     3 - 17.19d0 * DL1**3 + 71.08d0 * DL1**2 - 663.0d0 * DL1
502     4 - 192.4d0 * DL * DL1**2 + 80.41d0  * DL**2 * DL1
503     5 + NF * ( - 5.691d0 - 37.91d0 *Y
504     6 + 2.244d0 * DL**2 + 5.770d0 * DL
505     7 - 1.707d0* DL1**2  + 22.95d0 * DL1
506     8 + 3.036d0 * DL**2 * DL1 + 17.97d0 * DL * DL1 )     
507      RETURN
508      END
509
510      DOUBLE PRECISION FUNCTION C2NNLONSB (Y)
511      IMPLICIT DOUBLE PRECISION (A-Z)
512      INTEGER NF
513      COMMON/NFLAV/NF
514      DL1 = LOG (1.D0-Y)
515      DM  = 1./(1.D0-Y)
516      C2NNLONSB =
517     1  + 14.2222d0 * DL1**3 - 61.3333d0 * DL1**2- 31.105d0 * DL1
518     2  + 188.64d0
519     3  + NF * ( 1.77778d0 * DL1**2 - 8.5926d0 *  DL1 + 6.3489 )
520      C2NNLONSB = DM * C2NNLONSB
521      RETURN
522      END
523
524      DOUBLE PRECISION FUNCTION C2NNLONSC (Y)
525      IMPLICIT DOUBLE PRECISION (A-Z)
526      INTEGER NF
527      COMMON/NFLAV/NF
528      DL1 = LOG (1.D0-Y)
529      C2NNLONSC =
530     1 + 3.55555D0 * DL1**4 - 20.4444D0 * DL1**3 - 15.5525D0 * DL1**2
531     2 + 188.64D0 * DL1 - 338.531D0 + 0.537D0
532     3 + NF * (0.592593D0 * DL1**3 - 4.2963D0 * DL1**2
533     4 + 6.3489D0 * DL1 + 46.844D0 - 0.0035D0)
534      RETURN
535      END
536
537c------------------------ neutral current
538      DOUBLE PRECISION  FUNCTION C2NNLONSA_NC (Y)
539      IMPLICIT DOUBLE PRECISION (A-Z)
540      INTEGER NF
541      COMMON/NFLAV/NF
542      DL  = LOG (Y)
543      DL1 = LOG (1.D0-Y)
544      C2NNLONSA_NC =
545     1          - 69.59D0 - 1008.D0* Y
546     2          - 2.835D0 * DL**3 - 17.08D0 * DL**2 + 5.986D0 * DL
547     3          - 17.19D0 * DL1**3 + 71.08D0 * DL1**2 - 660.7D0 * DL1
548     4          - 174.8D0 * DL * DL1**2 + 95.09D0 * DL**2 * DL1
549     5        + NF * ( - 5.691D0 - 37.91D0 * Y
550     6          + 2.244D0 * DL**2 + 5.770D0 * DL
551     7          - 1.707D0 * DL1**2  + 22.95D0 * DL1
552     8          + 3.036D0 * DL**2 * DL1 + 17.97D0 * DL * DL1 )     
553      RETURN
554      END
555
556
557      DOUBLE PRECISION FUNCTION C2NNLONSC_NC (Y)
558      IMPLICIT REAL*8 (A-Z)
559      INTEGER NF
560      COMMON/NFLAV/NF
561      DL1 = LOG (1.D0-Y)
562      C2NNLONSC_NC =
563     1 + 3.55555D0 * DL1**4 - 20.4444D0 * DL1**3 - 15.5525D0 * DL1**2
564     2 + 188.64D0 * DL1 - 338.531D0 + 0.485D0
565     3 + NF * (0.592593D0 * DL1**3 - 4.2963D0 * DL1**2
566     4 + 6.3489D0 * DL1 + 46.844D0 - 0.0035D0)
567      RETURN
568      END
569
570
571c--------------------------------
572c--------------------------------
573c--------------------------------
574      double precision function Int2NLOa(z)
575      implicit none
576      doubLe precision z,c2nloa,sing
577      int2nloa=c2nloa(z)*sing(z)
578      return
579      end
580
581      double precision function Int2NLOb(z)
582      implicit none
583      doubLe precision z,c2nlob,singreg
584      int2nlob=c2nlob(z)*singreg(z)
585      return
586      end
587
588      double precision function Int2NLOg(z)
589      implicit none
590      double precision z,c2nlog,pdfg
591      int2nlog=c2nlog(z)*pdfg(z)
592      return
593      end
594
595      double precision function Int2NNLOsA(z)
596      implicit none
597      double precision z,sing, c2nnlosa
598      int2nnlosa=sing(z)*c2nnlosa(z)
599      return
600      end
601
602      double precision function Int2NNLOgA(z)
603      implicit none
604      double precision z,pdfg, c2nnloga
605      int2nnloga=pdfg(z)*c2nnloga(z)
606      return
607      end
608
609      double precision function Int2NNLOnsA(z)
610      implicit none
611      double precision z,nsing, c2nnlonsa, c2nnlonsa_nc
612      integer v
613      if (v.le.2) then
614      int2nnlonsa=nsing(z)*c2nnlonsa(z)
615      else
616      int2nnlonsa=nsing(z)*c2nnlonsa_nc(z)
617      end if
618      return
619      end
620
621      double precision function Int2NNLOnsB(z)
622      implicit none
623      double precision z,nsingreg, c2nnlonsB
624      int2nnlonsb=nsingreg(z)*c2nnlonsB(z)
625      return
626      end
627
628
629cccc need to correct for NC
630      double precision function F2(x,q,ord,v,zz)
631      implicit none
632      double precision x,q,zz, f(-6:6), quark, antiq,alphaspdf,pi,eps,
633     1  c2nnlogc,c2nnlonsc,c2nnlonsc_nc,c2nnlonscVAL
634      integer ord, nf,v,i,vv
635      double precision int2nloa,int2nlob, c2nlob,c2nloc, int2nlog
636      double precision int2nnlosa,int2nnloga,
637     1 int2nnlonsa, int2nnlonsb
638      common/prec/eps
639      double precision xx,qq
640      common/pdfpar/xx,qq
641      common/nflav/nf
642      common/vect/vv
643      integer f3c
644      common/f3call/f3c
645      double precision z,z0
646      z=(1d0-eps-x)*zz +x
647      z0=zz*x
648      f3c=0
649
650
651      xx=x
652      qq=q
653      vv=v
654      if(1d0-x.lt.eps) then
655      f2=0d0
656      else
657
658      quark=0d0
659      antiq=0d0
660      call pdf(x,q,f)
661      do i=1,nf
662      quark=quark+f(i)
663      antiq=antiq+f(-i)
664      end do
665
666      pi=3.1415926535
667
668      if (ord.eq.1) then
669      F2=x*(quark+antiq)
670
671      else if (ord.eq.2) then
672      f2=alphaspdf(q)/(4d0*pi)*2d0*x*((1d0-eps-x)*
673     1 (int2nloa(z)+int2nlob(z) +int2nlog(z))
674     2 -x*(quark+antiq)*c2nlob(z0)
675     3 +(quark+antiq)*c2nloc(z)
676     1 )
677     
678      else if (ord.eq.3) then
679       if (v.le.2) then
680        c2nnlonscVAL=c2nnlonsc(x)
681       else
682        c2nnlonscVAL=c2nnlonsc_nc(x)
683       end if
684      f2=(alphaspdf(q)/(4d0*pi))**2 *x*(
685     1         (1d0-eps-x)*( int2nnlosa(z)+int2nnloga(z)
686     2                      + int2nnlonsa(z)+int2nnlonsb(z))
687     3      + f(0)*c2nnlogc(x)
688     4      + (quark-antiq)*c2nnlonscVAL
689     5      )
690       endif
691
692      endif
693      return
694      end
695   
696 
697
698cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc--------------------------------------------------------------------------
699c------------------------   F3   ------------------------------------------
700
701      double precision function C3NLOb(z)
702      implicit none
703      double precision z,c2nlob
704      c3nlob=c2nlob(z)
705      return
706      end
707
708      double precision function C3NLOa(z)
709      implicit none
710      double precision z, c2nloa
711      c3nloa=c2nloa(z)-4d0/3d0 *(1d0+z)
712      return
713      end
714
715      DOUBLE PRECISION  FUNCTION C3NNLOSUMA (Y)
716      IMPLICIT DOUBLE PRECISION (A-Z)
717      INTEGER NF
718      COMMON/NFLAV/NF
719      DL  = LOG (Y)
720      DL1 = LOG (1.D0-Y)
721      C3NNLOSUMA =
722     1          - 206.1D0 - 576.8D0 * Y
723     2          - 3.922D0 * DL**3 - 33.31D0 * DL**2 - 67.60D0 * DL
724     3          - 15.20D0 * DL1**3 + 94.61D0 * DL1**2 - 409.6D0 * DL1
725     4          - 147.9D0 * DL * DL1**2
726     5          + NF * ( - 6.337D0 - 14.97D0 * Y
727     6          + 2.207D0 * DL**2 + 8.683D0 * DL
728     7          + 0.042D0 * DL1**3 - 0.808D0 * DL1**2 + 25.00D0 * DL1
729     8          + 9.684D0 * DL * DL1 )     
730      RETURN
731      END
732
733
734      DOUBLE PRECISION FUNCTION C3NNLODIFA (Y)
735      IMPLICIT DOUBLE PRECISION (A-Z)
736      INTEGER NF
737      COMMON/NFLAV/NF
738      DL  = LOG (Y)
739      DL1 = LOG (1.-Y)
740      C3NNLODIFA =
741     1          - 242.9D0 - 467.2D0 * Y
742     2          - 3.049D0 * DL**3 - 30.14D0 * DL**2 - 79.14D0 * DL
743     3          - 15.20D0 * DL1**3 + 94.61D0 * DL1**2 - 396.1D0 * DL1
744     4          - 92.43D0 * DL * DL1**2
745     5          + NF * ( - 6.337D0 - 14.97D0 * Y
746     6          + 2.207D0 * DL**2 + 8.683D0 * DL
747     7          + 0.042D0 * DL1**3 - 0.808D0 * DL1**2  + 25.00D0 * DL1
748     8          + 9.684D0 * DL * DL1 )     
749      RETURN
750      END
751
752      DOUBLE PRECISION  FUNCTION C3NNLOSUMC (Y)
753      IMPLICIT DOUBLE PRECISION (A-Z)
754      INTEGER NF
755      COMMON/NFLAV/NF
756      DL1 = LOG (1.D0-Y)
757      C3NNLOSUMC =
758     1    + 3.55555D0 * DL1**4 - 20.4444D0 * DL1**3 - 15.5525D0 * DL1**2
759     2        + 188.64D0 * DL1 - 338.531D0 - 0.104D0
760     3        + NF * (0.592593D0 * DL1**3 - 4.2963D0 * DL1**2
761     4        + 6.3489D0 * DL1 + 46.844D0 + 0.013D0)
762      RETURN
763      END
764
765
766      DOUBLE PRECISION FUNCTION C3NNLODIFC (Y)
767      IMPLICIT DOUBLE PRECISION (A-Z)
768      INTEGER NF
769      COMMON/NFLAV/NF
770      DL1  = LOG (1D0-Y)
771      C3NNLODIFC =
772     1  + 3.55555D0 * DL1**4 - 20.4444D0 * DL1**3 - 15.5525D0 * DL1**2
773     2        + 188.64D0 * DL1 - 338.531D0 - 0.104D0
774     3        + NF * (0.592593D0 * DL1**3 - 4.2963D0 * DL1**2
775     4        + 6.3489D0 * DL1 + 46.844D0 + 0.013D0)
776      RETURN
777      END
778
779      DOUBLE PRECISION FUNCTION C3NNLOB(Y)
780      IMPLICIT DOUBLE PRECISION (A-Z)
781      INTEGER NF
782      COMMON/NFLAV/NF
783      DL1 = LOG (1.-Y)
784      DM  = 1./(1.-Y)
785      C3NNLOB =
786     1   + 14.2222D0 * DL1**3 - 61.3333D0 * DL1**2 - 31.105D0 * DL1
787     2          + 188.64D0
788     3        + NF * ( 1.77778D0 * DL1**2 - 8.5926D0 * DL1 + 6.3489D0 )
789      C3NNLOB = DM * C3NNLOB     
790      RETURN
791      END
792
793      double precision function int3nloB(z)
794      implicit none
795      double precision z,c3nlob,nsingreg
796      int3nlob=c3nlob(z)*nsingreg(z)
797      return
798      end
799
800      double precision function int3nloA(z)
801      implicit none
802      double precision z,c3nloa,nsing
803      int3nloa=c3nloA(z)*nsing(z)
804      return
805      end
806
807      double precision function int3nnloB(z)
808      implicit none
809      double precision z,c3nnlob,nsingreg
810      int3nnlob=c3nnlob(z)*nsingreg(z)
811      return
812      end
813
814      double precision function int3nnloA(z)
815      implicit none
816      integer v
817      double precision z,c3nnloa,c3nnlosuma, c3nnlodifa,nsing
818      common/vect/v
819      if ((v.eq.1).or.(v.ge.3)) then
820      c3nnloa=(c3nnlosuma(z)+c3nnlodifa(z))/2d0
821      else if(v.eq.2) then
822      c3nnloa=(c3nnlosuma(z)-c3nnlodifa(z))/2d0
823      end if
824      int3nnloa=c3nnloA*nsing(z)
825      return
826      end
827
828
829
830cccc need to correct for NC
831      double precision function F3(x,q,ord,v,zz)
832      implicit none
833      double precision x,q,zz,f(-6:6),quark,antiq,alphaspdf,pi,eps,
834     1 c3nnloc, c3nnlodifc, c3nnlosumc
835      integer ord, nf,v,i,vv
836      double precision int3nloa, int3nlob, c3nlob,c2nloc
837      double precision int3nnloa, int3nnlob
838      common/prec/eps
839      double precision xx,qq
840      common/pdfpar/xx,qq
841      common/nflav/nf
842      common/vect/vv
843      integer f3c
844      common/f3call/f3c
845      double precision z,z0
846      z=(1d0-eps-x)*zz +x
847      z0=zz*x
848      f3c=1
849
850
851      xx=x
852      qq=q
853      vv=v
854      if(1d0-x.lt.eps) then
855      f3=0d0
856      else
857
858      if (v.eq.4) then
859      f3=0d0
860      else
861
862      quark=0d0
863      antiq=0d0
864      call pdf(x,q,f)
865      do i=1,nf
866      quark=quark+f(i)
867      antiq=antiq+f(-i)
868      end do
869
870      pi=3.1415926535
871
872      if (ord.eq.1) then
873      F3=(quark-antiq)
874
875      else if (ord.eq.2) then
876       f3=alphaspdf(q)/(4d0*pi)*2d0*(
877     1  (1d0-eps-x)*(int3nloa(z)+int3nlob(z))
878     2  -x*(quark-antiq)*c3nlob(z0)
879     3  +(quark-antiq)*c2nloc(z)
880     4    )
881
882      else if (ord.eq.3) then
883       if ((v.eq.1).or.(v.ge.3)) then
884        c3nnloc=(c3nnlosumc(x)+c3nnlodifc(x))/2d0
885       else if(v.eq.2) then
886        c3nnloc=(c3nnlosumc(x)-c3nnlodifc(x))/2d0
887       end if
888      f3=(alphaspdf(q)/(4d0*pi))**2 *(
889     1  (1d0-eps-x)*(int3nnloa(z)  +int3nnlob(z))
890     2 + (quark-antiq)*c3nnloc
891     3  )
892
893
894c     1 daind(x,1d0, int3nnloa,eps,key,maxcnt,cnt, err )
895c     2 +daind(x,1d0, int3nnlob,eps,key,maxcnt,cnt, err )
896c     3 + (quark-antiq)*c3nnloc
897c     4  )
898      endif
899     
900      endif
901      endif
902      return
903      end