Skip to content

Commit

Permalink
fixing bug in src/DIS/ComputeDISOperators.f in the computation of the…
Browse files Browse the repository at this point in the history 
… IC in the FONLL scheme for the NC structure functions
  • Loading branch information
Valerio Bertone committed Dec 14, 2015
1 parent 7e4e87c commit 1530ca8
Showing 1 changed file with 64 additions and 68 deletions.
132 changes: 64 additions & 68 deletions src/DIS/ComputeDISOperators.f
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,7 +38,7 @@ subroutine ComputeDISOperators(Q)
* Internal Variables
*
integer jgrid,ipdf,ihq
integer pt,ipt_FF,ipt_FF_min,iptbkp
integer pt,ipt_FF,iptbkp,ipt_max_IC
integer alpha,beta,gamma
integer nf
integer gbound
Expand Down Expand Up @@ -94,13 +94,7 @@ subroutine ComputeDISOperators(Q)
* (Remember that a_QCD takes as an argument the factorization scale
* and converts it internally into the renormalization scale).
*
muF = kfacQ*Q2
*
* Normally the NC coefficient functions in the FFNS start at O(alpha_s).
* In the presence of intrinsic charm insted they start at O(1).
*
ipt_FF_min = 1
if(IntrinsicCharm) ipt_FF_min = 0
muF = kfacQ * Q2
*
* Scale down the perturbative order of the alphas evolution if one
* of the FFNSs has been chosen
Expand Down Expand Up @@ -212,10 +206,13 @@ subroutine ComputeDISOperators(Q)
call ComputeChargesDIS(Q2,bq,dq,bqt)
*
* In case of intrinsic charm, compute the correction to apply
* to the non-singlet part of the massive F_L
* to the non-singlet part of the massive F_L and also limit the
* compuation of the IC component to NLO.
*
if(IntrinsicCharm)
1 cFLIC = ( 1d0 - ( 1d0 - bqt(4) / bq(4) ) / 2d0 )
if(IntrinsicCharm)then
cFLIC = ( 1d0 - ( 1d0 - bqt(4) / bq(4) ) / 2d0 )
ipt_max_IC = min(ipt,1)
endif
*
do jgrid=1,ngrid
*
Expand Down Expand Up @@ -341,7 +338,7 @@ subroutine ComputeDISOperators(Q)
if(Nf_FF.lt.6)then
do ihq=Nf_FF+1,6
if(Q2.ge.m2th(ihq))then
do pt=ipt_FF_min,ipt
do pt=1,ipt
C2g(ihq) = C2g(ihq) + as(pt)
1 * ( c0(ihq)
2 * SC2m0NC(jgrid,ixi(ihq),
Expand Down Expand Up @@ -379,7 +376,7 @@ subroutine ComputeDISOperators(Q)
*
if(IntrinsicCharm)then
if(Nf_FF.lt.4)then
do pt=ipt_FF_min,ipt
do pt=0,ipt_max_IC
CLnsp(4) = CLnsp(4) + as(pt)
1 * ( c0(4)
2 * SCLm0NC(jgrid,ixi(4),
Expand Down Expand Up @@ -459,7 +456,7 @@ subroutine ComputeDISOperators(Q)
if(Nf_FF.lt.6)then
do ihq=Nf_FF+1,6
if(W2.ge.4d0*m2th(ihq))then
do pt=ipt_FF_min,ipt
do pt=1,ipt
C2g(ihq) = C2g(ihq) + as(pt)
1 * ( c0(ihq)
2 * SC2mNC(jgrid,ixi(ihq),
Expand Down Expand Up @@ -497,7 +494,7 @@ subroutine ComputeDISOperators(Q)
*
if(IntrinsicCharm)then
if(Nf_FF.lt.4)then
do pt=ipt_FF_min,ipt
do pt=0,ipt_max_IC
CLnsp(4) = CLnsp(4) + as(pt)
1 * ( c0(4)
2 * SCLmNC(jgrid,ixi(4),
Expand Down Expand Up @@ -599,7 +596,7 @@ subroutine ComputeDISOperators(Q)
if(Nf_FF.lt.6)then
do ihq=Nf_FF+1,6
if(W2.ge.4d0*m2th(ihq))then
do pt=ipt_FF_min,ipt_FF
do pt=1,ipt_FF
C2g(ihq) = C2g(ihq) + as(pt)
1 * ( c0(ihq)
2 * SC2mNC(jgrid,ixi(ihq),
Expand Down Expand Up @@ -631,7 +628,7 @@ subroutine ComputeDISOperators(Q)
enddo
endif
if(Q2.ge.m2th(ihq))then
do pt=ipt_FF_min,ipt_FF
do pt=1,ipt_FF
C2g(ihq) = C2g(ihq) + as(pt)
1 * ( - damp(ihq) * ( c0(ihq)
2 * SC2m0NC(jgrid,ixi(ihq),
Expand Down Expand Up @@ -670,7 +667,7 @@ subroutine ComputeDISOperators(Q)
if(IntrinsicCharm)then
if(Nf_FF.lt.4)then
* FFNS
do pt=ipt_FF_min,ipt
do pt=0,ipt_max_IC
CLnsp(4) = CLnsp(4) + as(pt)
1 * ( c0(4)
2 * SCLmNC(jgrid,ixi(4),
Expand All @@ -688,16 +685,16 @@ subroutine ComputeDISOperators(Q)
enddo
CLnsp(4) = cFLIC * CLnsp(4)
* FFN0
do pt=ipt_FF_min,ipt
do pt=0,ipt_max_IC
c CLnsp(4) = CLnsp(4) - as(pt)
c 1 * damp(ihq) * ( c0(4)
c 1 * damp(4) * ( c0(4)
c 2 * SCLm0NC(jgrid,ixi(4),
c 3 3,pt,alpha,beta)
c 4 + c1(4)
c 5 * SCLm0NC(jgrid,ixi(4)+1,
c 6 3,pt,alpha,beta) )
C2nsp(4) = C2nsp(4) - as(pt)
1 * damp(ihq) * ( c0(4)
1 * damp(4) * ( c0(4)
2 * SC2m0NC(jgrid,ixi(4),
3 3,pt,alpha,beta)
4 + c1(4)
Expand Down Expand Up @@ -941,44 +938,44 @@ subroutine ComputeDISOperators(Q)
enddo
enddo
enddo
*
*
* Charged current structure functions
*
*
elseif(ProcessDIS.eq."CC")then
*
*
* Different projectiles
*
*
if(ProjectileDIS(1:8).eq."neutrino".or.
1 ProjectileDIS(1:8).eq."positron")then
ipr = 1
elseif(ProjectileDIS.eq."antineutrino".or.
1 ProjectileDIS(1:8).eq."electron")then
ipr = - 1
endif
*
*
* Define useful constants
*
*
Kl = V_ud2 + V_us2
Kc = V_cd2 + V_cs2
Kb = V_ub2 + V_cb2
Kt = V_td2 + V_ts2 + V_tb2
*
*
do jgrid=1,ngrid
*
*
* If an external grid is found compute the whole operator
*
*
gbound = 0
if(IsExt(jgrid)) gbound = nin(jgrid)
*
*
do alpha=0,gbound
do beta=alpha,nin(jgrid)
*
* Final state energy
*
W2 = Q2 * ( 1d0 - xg(jgrid,alpha) ) / xg(jgrid,alpha)
*
*
* Put together coefficient functions
*
*
do ihq=3,6
C2g(ihq) = 0d0
C2ps(ihq) = 0d0
Expand All @@ -992,9 +989,9 @@ subroutine ComputeDISOperators(Q)
C3nsp(ihq) = 0d0
C3nsm(ihq) = 0d0
enddo
*
*
* Construct coefficient functions according to the scheme chosen
*
*
if(MassScheme.eq."ZM-VFNS")then
do pt=0,ipt
C2g(3) = C2g(3)
Expand Down Expand Up @@ -1072,8 +1069,7 @@ subroutine ComputeDISOperators(Q)
*
* FFNS is NLO at most
*
ipt_FF = ipt
if(ipt.gt.1) ipt_FF = 1
ipt_FF = min(ipt,1)
*
* Heavy quark coefficient functions
*
Expand Down Expand Up @@ -1205,7 +1201,7 @@ subroutine ComputeDISOperators(Q)
enddo
endif
*
* FFNS is NLO at most
* FNNS is NLO at most
*
ipt_FF = ipt
if(ipt.gt.1) ipt_FF = 1
Expand Down Expand Up @@ -1545,9 +1541,9 @@ subroutine ComputeDISOperators(Q)
OpF2(jgrid,3,12,alpha,beta) = Kl * C2nsp(3) / 10d0
* T35
OpF2(jgrid,3,13,alpha,beta) = Kl * C2nsp(3) / 15d0
*
*
* Charm Component
*
*
* Singlet
OpF2(jgrid,4,1,alpha,beta) =
1 2d0 * Kc * ( C2ps(4) + C2nsp(4) / 6d0 )
Expand Down Expand Up @@ -1797,9 +1793,9 @@ subroutine ComputeDISOperators(Q)
1 OpF2(jgrid,4,13,alpha,beta) +
2 2d0 * Kc * c2ccIC / 60d0
endif
*
*
* Bottom Component
*
*
* Singlet
OpF2(jgrid,5,1,alpha,beta) =
1 2d0 * Kb * ( C2ps(5) + C2nsp(5) / 6d0 )
Expand Down Expand Up @@ -1993,9 +1989,9 @@ subroutine ComputeDISOperators(Q)
2 + Kb * C2qm(5) / 15d0
endif
endif
*
*
* Top Component
*
*
* Singlet
OpF2(jgrid,6,1,alpha,beta) =
1 2d0 * Kt * ( C2ps(6) + C2nsp(6) / 6d0 )
Expand Down Expand Up @@ -2215,21 +2211,21 @@ subroutine ComputeDISOperators(Q)
2 - 2d0 * Kt * C2qm(6) / 15d0
endif
endif
*
*
* Total
*
*
do ipdf=0,13
do ihq=3,6
OpF2(jgrid,7,ipdf,alpha,beta) =
1 OpF2(jgrid,7,ipdf,alpha,beta)
2 + OpF2(jgrid,ihq,ipdf,alpha,beta)
enddo
enddo
*
*
* FL
*
*
* Light Component
*
*
* Singlet
OpFL(jgrid,3,1,alpha,beta) =
1 2d0 * Kl * ( CLps(3) + CLnsp(3) / 6d0 )
Expand All @@ -2253,9 +2249,9 @@ subroutine ComputeDISOperators(Q)
OpFL(jgrid,3,12,alpha,beta) = Kl * CLnsp(3) / 10d0
* T35
OpFL(jgrid,3,13,alpha,beta) = Kl * CLnsp(3) / 15d0
*
*
* Charm Component
*
*
* Singlet
OpFL(jgrid,4,1,alpha,beta) =
1 2d0 * Kc * ( CLps(4) + CLnsp(4) / 6d0 )
Expand Down Expand Up @@ -2495,9 +2491,9 @@ subroutine ComputeDISOperators(Q)
1 OpFL(jgrid,4,13,alpha,beta) +
2 2d0 * Kc * cLccIC / 60d0
endif
*
*
* Bottom Component
*
*
* Singlet
OpFL(jgrid,5,1,alpha,beta) =
1 2d0 * Kb * ( CLps(5) + CLnsp(5) / 6d0 )
Expand Down Expand Up @@ -2691,9 +2687,9 @@ subroutine ComputeDISOperators(Q)
2 + Kb * CLqm(5) / 15d0
endif
endif
*
*
* Top Component
*
*
* Singlet
OpFL(jgrid,6,1,alpha,beta) =
1 2d0 * Kt * ( CLps(6) + CLnsp(6) / 6d0 )
Expand Down Expand Up @@ -2913,21 +2909,21 @@ subroutine ComputeDISOperators(Q)
2 - 2d0 * Kt * CLqm(6) / 15d0
endif
endif
*
*
* Total
*
*
do ipdf=0,13
do ihq=3,6
OpFL(jgrid,7,ipdf,alpha,beta) =
1 OpFL(jgrid,7,ipdf,alpha,beta)
2 + OpFL(jgrid,ihq,ipdf,alpha,beta)
enddo
enddo
*
*
* F3
*
*
* Light Component
*
*
* Gluon
OpF3(jgrid,3,2,alpha,beta) = ipr * 2d0 * Kl * C3g(3)
* Valence
Expand All @@ -2950,9 +2946,9 @@ subroutine ComputeDISOperators(Q)
* T8
OpF3(jgrid,3,10,alpha,beta) = - ipr
1 * V_us2 * C3nsp(3) / 2d0
*
*
* Charm Component
*
*
* Gluon
OpF3(jgrid,4,2,alpha,beta) = ipr * 2d0 * Kc * C3g(4)
* Valence
Expand Down Expand Up @@ -3191,9 +3187,9 @@ subroutine ComputeDISOperators(Q)
1 OpF3(jgrid,4,13,alpha,beta) -
2 ipr * 2d0 * Kc * c3ccIC / 60d0
endif
*
*
* Bottom Component
*
*
* Gluon
OpF3(jgrid,5,2,alpha,beta) = ipr * 2d0 * Kb * C3g(5)
* Valence
Expand Down Expand Up @@ -3386,9 +3382,9 @@ subroutine ComputeDISOperators(Q)
2 - ipr * Kb * C3qm(5) / 4d0
endif
endif
*
*
* Top Component
*
*
* Gluon
OpF3(jgrid,6,2,alpha,beta) = ipr * 2d0 * Kt * C3g(6)
* Valence
Expand Down Expand Up @@ -3608,9 +3604,9 @@ subroutine ComputeDISOperators(Q)
2 + ipr * Kt * C3qm(6) / 5d0
endif
endif
*
*
* Total
*
*
do ipdf=0,13
do ihq=3,6
OpF3(jgrid,7,ipdf,alpha,beta) =
Expand Down

0 comments on commit 1530ca8

Please sign in to comment.