API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Oct 10th, 2018
200
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
Fortran 56.34 KB | None | 0 0
raw download clone embed print report
  1. c
  2. c******************************************************************************
  3. c      
  4.       function dsigma(e1,costheta,omega)
  5. c
  6. c.....Returns the double differential neutrino-carbon cross section
  7. c     in units of 10^38 cm^2/GeV. The calculation is performed at
  8. c     fxed neutrino energy and lepton scattering angle, as a function
  9. c     of energy transfer.
  10. c     The spectral function is normalized to (A-Z) = Z = A/2
  11. c
  12. c     input variables
  13. c
  14. c     e1        : beam energy [GeV]
  15. c     costheta  : cos(theta), theta being the lepton scattering angle
  16. c     omega     : lepton energy loss [GeV]
  17. c
  18. c     Last revised: Roma, September 2016
  19. c
  20. c******************************************************************************
  21. c
  22.       implicit none
  23. c
  24.       integer nnn,nsamples,icount,iev,ig,ilept,iflag
  25. c      
  26.       parameter( nnn = 500000 )      
  27. c
  28.       real*8 e1,costheta,omega,W_MeV,e_lept,p_lept,xm_lept,zero      
  29.      $      ,xmp,xmn,xm,xmA,xmAm1,xmR_0,EA,ER,E_min,sigma_QE
  30.      $      ,Q2,q,q_MeV,dsigma ,p0,p0_MeV,xE,costheta_p,eps
  31.      $      ,xmRstar,pf,Ef,pf_MeV,Erec,arg2,diracdelta,e1_MeV
  32.      $      ,omega_thr,small,arg,stepf,k_F,sigma_RES,sigma_DIS
  33.      $      ,dsigma_QE,dsigma_RES,dsigma_DIS
  34. c
  35.       real*8 x(nnn),y(nnn),z(nnn)
  36. c      
  37.       common / config1 / nsamples
  38.       common / config2 / x,y,z      
  39. c
  40.       common / constants / zero,xmp,xmn,xm,xmA,xmAm1,xmR_0,xm_lept
  41.      $                    ,EA,ER,E_min,eps,k_F
  42. c
  43.       common / channels / dsigma_QE,dsigma_RES,dsigma_DIS
  44. c      
  45.       common / flags / ig,ilept
  46. c
  47.       small = 1.e-10      
  48. c
  49. c.....set lepton kinematics
  50. c      
  51.       e1_MeV = e1*1000.
  52.       w_MeV = omega*1000.
  53.       e_lept = e1 - omega                         ! energy of charged lepton [GeV]
  54.       arg = e_lept**2 - xm_lept**2
  55.         if(arg.ge.zero)then
  56.         p_lept = sqrt( e_lept**2 - xm_lept**2 )  ! momentum of charged lepton[GeV]
  57.         else
  58.         dsigma = zero
  59.         return      
  60.         end if  
  61.       Q2 = 2.*e1*e_lept*(1. - (p_lept/e_lept)*costheta )  - xm_lept**2 ! Q^2 [GeV^2]
  62.       omega_thr = Q2/2./xmA
  63.         if(omega.le.omega_thr)then
  64.         dsigma = zero
  65.         return
  66.         end if        
  67.       q = sqrt( Q2 + omega**2 )         ! three-momentum transfer [GeV]
  68.       q_MeV = 1000.*q
  69. c
  70. c.....start loop on samples
  71. c
  72.           dsigma_QE = zero
  73.           dsigma_RES = zero
  74.           dsigma_DIS = zero
  75. c          
  76.           do 1000 iev = 1,nsamples
  77. c
  78.           p0 = x(iev)                   ! nucleon momentum [GeV]
  79.           p0_MeV = 1000.*p0             ! change units to MeV
  80.           xE = y(iev)                   ! nucleon removal energy [GeV]
  81.           costheta_p = z(iev)           ! polar angle of the nucleon momentum
  82. c
  83.           xmRstar = xmR_0 + xE
  84.           Erec = sqrt(xmRstar**2+p0**2) ! energy of the recoiling system [GeV]
  85. c          
  86. c=========================================================================
  87. c
  88. c.....CCQE
  89. c
  90. c=========================================================================          
  91. c
  92. c.....set argument of energy conserving delta-function
  93. c          
  94.           pf = sqrt( p0**2 + q**2 + 2.*p0*q*costheta_p ) ! momentum of final state proton [GeV]
  95.           Ef = sqrt(pf**2 + xm**2)                       ! energy of final state proton [GeV]
  96.           pf_MeV = 1000.*pf                              ! change units to MeV
  97.           arg2 = omega + xmA - Ef - Erec                 ! argument of energy-conserving delta-function [GeV]      
  98. c
  99. c.....elementary CCQE cross section in units of millbarn/sr    
  100. c
  101.           sigma_QE = zero
  102.           call sig_nuN_CCQE(ilept,ig,q_MeV,w_MeV,p0_MeV,pf_MeV,e1_MeV
  103.      $        ,costheta,costheta_p,sigma_QE)
  104. c
  105. c         write(96,6161)q_MeV,w_MeV,sigma_QE
  106. c6161     format(3e15.5)
  107. c
  108. c.....change units mbarn/sr --> nbarn/sr --> 10^-38 cm^2/sr
  109. c          
  110.           sigma_QE = 1.e6*sigma_QE
  111.           sigma_QE = 1.e5*sigma_QE
  112. c        
  113.           stepf = 1.            ! step-function for Pauli blocking
  114.           if(pf.lt.k_F)then
  115.           stepf = 0.
  116.           end if
  117. c
  118. c......double-differential cross section in 10^-38 cm^2/sr/GeV
  119. c          
  120.           dsigma_QE = dsigma_QE + sigma_QE*diracdelta(eps,arg2)*stepf
  121. c
  122. c=========================================================================          
  123. c          
  124. c.....RES
  125. c
  126. c=========================================================================
  127. c
  128. c.....elementary resonance production cross section in units of millbarn/sr/GeV
  129. c
  130.           sigma_RES = zero
  131. c         call sig_nuN_RES(ilept,Erec,q_MeV,w_MeV,p0_MeV,pf_MeV,e1_MeV
  132. c    $        ,costheta,costheta_p,sigma_RES)
  133. c
  134. c.....change units mbarn/sr/GeV --> nbarn/sr/GeV --> 10^-38 cm^2/sr/GeV
  135. c          
  136.           sigma_RES = 1.e6*sigma_RES
  137.           sigma_RES = 1.e5*sigma_RES
  138. c
  139.           dsigma_RES = dsigma_RES + sigma_RES
  140. c
  141. c=========================================================================          
  142. c
  143. c.....DIS
  144. c
  145. c=========================================================================
  146. c          
  147. c.....elementary DIS cross section in units of millbarn/sr/GeV
  148. c
  149.           if(iev.eq.1)then
  150.              iflag = 0
  151.           else
  152.              iflag = 1
  153.           end if
  154. c
  155.           sigma_DIS = zero
  156. c         call sig_nuN_DIS(ilept,iflag,Erec,q_MeV,w_MeV,p0_MeV,pf_MeV
  157. c    $        ,e1_MeV,costheta,costheta_p,sigma_DIS)
  158. c
  159. c.....change units mbarn/sr/GeV --> nbarn/sr/GeV --> 10^-38 cm^2/sr/GeV
  160. c          
  161.           sigma_DIS = 1.e6*sigma_DIS
  162.           sigma_DIS = 1.e5*sigma_DIS
  163. c
  164.           dsigma_DIS = dsigma_DIS + sigma_DIS
  165. c
  166.  1000     continue        
  167. c
  168. c.....end loop on samples
  169. c
  170.       dsigma_QE  = 22.*dsigma_QE /nsamples
  171.       dsigma_RES = 22.*dsigma_RES/nsamples
  172.       dsigma_DIS = 40.*dsigma_DIS/nsamples
  173. c
  174.       dsigma = dsigma_QE + dsigma_RES + dsigma_DIS      
  175. c
  176.       if(dsigma.le.small)then
  177.       dsigma = zero
  178.       end if      
  179. c
  180.       return
  181. c
  182.       end
  183. c
  184. c******************************************************************************
  185. c      
  186.       subroutine sig_nuN_CCQE(ilept,ig,xq,w,xk,xp,e_nu,cos_theta
  187.      $          ,costheta_p,sig)
  188. c
  189. c.....returns the cross section of the process
  190. c     nu_ell + neutron ---> ell^- + proton in units of millibarn/sr
  191. c     all input energies and momenta are in MeV
  192. c
  193. c     ilept = 1 muon neutrino beam
  194. c     ilept = 2 electron neutrino beam
  195. c
  196.       implicit none
  197. c      
  198.       integer ilept,ig,i
  199. c      
  200.       real*8 A(5),WW(5)
  201. c
  202.       real*8 G_F,pi,c0,hc,xm_n,xm_p,xm_lept,e_lept,e_nu,xq,w,xk,xp
  203.      $      ,cos_theta,sig,xk_lept,c1,c2,c3,q2,Ek,Ep,c4
  204.      $      ,xkx2,xkz,xkz2,aux1,wt,q2t,summ,prod,diff,contr,ctot
  205.      $      ,xm_lept_MeV,zero,xm,xmA,xmAm1,xmR_0,EA,ER,E_min,eps
  206.      $      ,xmp,xmn,k_F,costheta_p,sinqtheta_p
  207. c
  208.       common / constants / zero,xmp,xmn,xm,xmA,xmAm1,xmR_0,xm_lept
  209.      $                    ,EA,ER,E_min,eps,k_F
  210. c
  211. c.....Fermi constant in fm^2 (Cabibbo's angle included)
  212. c
  213.       G_F = 1.14e-5*.1973**2
  214. c
  215.       pi = 4.*atan(1.)
  216.       c0 = 1./16./pi/pi
  217. c
  218.       hc = 197.3
  219.       xm_n = 1000.*xmn
  220.       xm_p = 1000.*xmp
  221.       xm_lept_MeV = 1000.*xm_lept
  222. c
  223.       e_lept = e_nu - w ! energy of charged lepton [MeV]
  224.       xk_lept = sqrt(e_lept**2 - xm_lept_MeV**2) ! momentum of charged lepton [MeV]
  225.       c1 = (xk_lept/e_lept)*cos_theta  
  226.       c2 = e_nu*e_lept      
  227.       c3 = xk_lept/e_nu      
  228. c
  229.       q2 = w**2 - xq**2    
  230.       Ek = sqrt(xk**2 + xm_n**2) ! neutron energy [MeV]
  231.       Ep = sqrt(xp**2 + xm_p**2) ! proton energy [MeV]
  232.       c4 = 1./4./Ek/Ep
  233.       sinqtheta_p = 1. - costheta_p**2
  234.       xkx2 = .5*xk*xk*sinqtheta_p
  235.       xkz = xk*costheta_p
  236.       xkz2= xkz**2            
  237.       aux1 = xm_lept_MeV**2 - q2  
  238.       wt = Ep - Ek
  239.       q2t = wt**2 - xq**2
  240. c
  241.       summ = (w/xq)*( 1. + c1 + xm_lept_MeV**2/w/e_lept )  ! (p_nu)z/E_nu + (p_mu)z/E_mu
  242. c      
  243.       prod = (1./xq/xq)*( w*w*c1 - c2*( 1. - c1)**2    
  244.      $     + xm_lept_MeV**2*( w/e_lept + 1. - c1 ) )       ! (p_nu)z/E_nu * (p_mu)z/E_mu
  245. c      
  246.       diff = ( (e_nu + e_lept)/xq )*(1. - c1)
  247.      $     - xm_lept_MeV**2/e_lept/xq                      ! (p_nu)z/E_nu - (p_mu)z/E_mu
  248. c    
  249.       A(1) = aux1/2.
  250.       A(2) = c2*( Ek**2 - Ek*xkz*summ  
  251.      $     + xkx2*(c2/xq/xq)*( (xk_lept/e_lept)**2 - c1**2 )
  252.      $     + xkz2*prod - (xm_n**2/2.)*( 1. - c1 ) )
  253.       A(3) = c2*( Ek*(xkz+xq) - Ep*xkz )*diff
  254.       A(4) = c2*(wt**2 - wt*xq*summ + xq*xq*prod - 0.5*q2t*( 1. - c1 ))
  255.       A(5) = c2*( 2.*Ek*wt + 2.*xkz*xq*prod - ( Ek*xq + wt*xkz )*summ
  256.      $     - ( Ek*wt - xkz*xq )*( 1. - c1 ) )
  257. c
  258.       call sf_CCQE(q2t,WW)
  259. c
  260. c     write(99,9999)q2t*1.e-6,WW(1),WW(2),WW(3),WW(4),WW(5)
  261. c9999 format(6e15.5)
  262. c
  263. c.....compute the contraction L_munu*W^munu in MeV^4 and
  264. c     change units to fm^-4
  265. c
  266.       contr = 0.
  267.       do 1000 i = 1,5
  268.  1000 contr = contr + 16.*A(i)*WW(i)/hc**4
  269. c
  270. c.....calculate coefficient in units of fm^6
  271. c
  272.       ctot = (G_F**2/2.)*c0*c3*(c4*hc*hc)
  273. c
  274. c.....get x-section in units of mbarn/sr
  275. c
  276.       sig = 10.*ctot*contr
  277. c
  278.       return
  279.       end
  280. c
  281. c******************************************************************************
  282. c      
  283.       subroutine sf_CCQE(q2,W)
  284. c
  285.       implicit real*8(a-h,o-z)
  286. c
  287.       common / flags / ig,ilept
  288. c
  289. c.....returns the structure functions W1, W2/m^2
  290. c     W3/m^2,  W4/m^2 and W5/m^2 as a function
  291. c     of the squared four-momentum transfer q2
  292. c
  293.       dimension W(5)
  294. c
  295.       xm = 939.57
  296. c
  297.       q2g = q2*1.e-6
  298.       call form_factors(ig,q2g,F1,F2,FA,FP)
  299. c     FP = 0.
  300. c
  301.       W(1) = 2.*( -(q2/2.)*( F1 + F2 )**2
  302.      $      + ( 2.*xm**2 - q2/2. )*FA**2 )
  303.       W(2) = 4.*( F1**2 -(q2/4./xm/xm)*F2**2 + FA**2 )
  304.       W(3) = -4.*( F1 + F2 )*FA
  305.       W(4) = -2.*( F1*F2 + ( 2.*xm**2 + q2/2. )*F2**2/4./xm/xm
  306.      $      + (q2/2.)*FP**2 - 2.*xm*FP*FA )
  307.       W(5) = W(2)/2.
  308. c
  309.       return
  310.       end
  311. c
  312. c******************************************************************************
  313. c      
  314.       subroutine form_factors(ig,q2,F_1,F_2,F_A,F_p)
  315. c
  316. c.....returns the form factors. The squared four momentum
  317. c     transfer q2 is given in units of GeV^2
  318.  
  319.       implicit real*8(a-h,o-z)
  320. c
  321.       integer ig      
  322. c
  323.       real*8 g_A,M_A,Q2fm,hc      
  324. c
  325.       hc = .1973
  326.       g_A = 1.267      
  327.       M_A = 1.03
  328. c
  329.       xm = .93957
  330.       tau = q2/4./xm/xm
  331.       Q2fm = -q2/hc/hc
  332.       cc = 1./(1. - tau)
  333. c
  334.       call nform(ig,Q2fm,G_ep,G_en,G_mp,G_mn)
  335. c
  336.       F_1_p = cc*( G_ep - tau*G_mp )
  337.       F_2_p = cc*(-G_ep +     G_mp )
  338.       F_1_n = cc*( G_en - tau*G_mn )
  339.       F_2_n = cc*(-G_en +     G_mn )
  340. c      
  341.       F_1 = F_1_p - F_1_n
  342.       F_2 = F_2_p - F_2_n
  343.       F_A = -g_A/(1. - q2/M_A/M_A)**2
  344.       F_p = 2.*F_A*xm / ( .1396**2 + q2)*0.001
  345. c     write(91,9191)q2,2.*F_A*xm,( .1396**2 + q2),F_p
  346. c9191 format(4e18.6)
  347. c
  348.       return
  349.       end
  350. c
  351. c*******************************************************************************
  352. c      
  353.       subroutine nform(ichoice,Q2,G_ep,G_en,G_mp,G_mn)
  354. c
  355. c.....nucleon form factors      
  356. c      
  357. c     revised January 2014, Virginia Tech
  358. c
  359. c     ichoice = 1  Dipole parametrization
  360. c     ichoice = 2  Kelly's parametrization [PRC 70, 068282 (2004)]
  361. c     ichoice = 3  BBBA parametrization [NPB Proc. Suppl., 159, 127 (2006)]
  362. c
  363. c     Q2 is given in units of fm^-2
  364. c
  365.       implicit none
  366. c
  367.       integer ichoice,i,j
  368. c
  369.       real*8 Q2,Q2G,G_ep,G_mp,G_en,G_mn,G_D,Lambdasq,hc,pm,zero
  370.      $      ,mu_p,mu_n,a_ep,a_mp,a_mn,tau,a_en,b_en,num_ep,den_ep
  371.      $      ,num_mp,den_mp,num_en,den_en,num_mn,den_mn,one,nm,xm
  372.       real*8 b_ep(3),b_mp(3),b_mn(3),aa_ep(4),aa_mp(4),aa_en(4)
  373.      $      ,aa_mn(4),bb_ep(4),bb_mp(4),bb_en(4),bb_mn(4)
  374. c      
  375.       data a_ep , a_mp ,a_mn  / -0.24 , 0.12 , 2.33 /
  376.       data b_ep / 10.89 , 12.82 , 21.97 /
  377.       data b_mp / 10.97 , 18.86 , 6.55 /
  378.       data b_mn / 14.72 , 24.20 , 84.1 /
  379.       data a_en , b_en / 1.70 , 3.30 /
  380. c
  381.       data aa_ep / 1.0 , -0.0578 , 0.0  , 0.0 /    
  382.       data aa_mp / 1.0 ,  0.0150 , 0.0  , 0.0 /
  383.       data aa_en / 0.0 ,  1.2500 , 1.30 , 0.0 /
  384.       data aa_mn / 1.0 ,  1.8100 , 0.0  , 0.0 /
  385. c
  386.       data bb_ep / 11.10 ,  13.60 ,   33.00 ,   0.0 /
  387.       data bb_mp / 11.10 ,  19.60 ,    7.54 ,   0.0 /
  388.       data bb_en / 9.86  , 305.00 , -758.00 , 802.0 /
  389.       data bb_mn / 14.10 ,  20.70 ,   68.70 ,   0.0 /
  390. c      
  391.       hc = .1973
  392.       pm = .93828         ! proton mass in GeV
  393.       nm = .93957         ! neutron mass in GeV
  394.       xm = .5*(pm + nm)   ! isoscalar nucleon mass in GeV
  395. c
  396.       zero = 0.
  397.       one = 1.
  398.       mu_p = 2.79278   ! proton magnetic moment
  399.       mu_n = -1.91315  ! neutron magnetic moment
  400.       Q2G = Q2*hc*hc
  401.       tau = Q2G/4./pm/pm
  402. c
  403.       Lambdasq = 0.71
  404.       G_D  = 1./(1. + Q2G/Lambdasq)**2
  405. c      
  406.       if(ichoice.eq.1)then
  407. c              
  408.         G_ep = G_D
  409.         G_en = zero
  410.         G_mp = mu_p*G_D
  411.         G_mn = mu_n*G_D
  412.       else if(ichoice.eq.2)then
  413.         G_ep = (1. + a_ep*tau)/
  414.      $         (1. + b_ep(1)*tau + b_ep(2)*tau**2 + b_ep(3)*tau**3)
  415.         G_en = G_D*a_en*tau/(1. + b_en*tau)
  416.         G_mp = mu_p*(1. + a_mp*tau)/
  417.      $              (1. + b_mp(1)*tau + b_mp(2)*tau**2 + b_mp(3)*tau**3)
  418.         G_mn = mu_n*(1. + a_mn*tau)/
  419.      $              (1. + b_mn(1)*tau + b_mn(2)*tau**2 + b_mn(3)*tau**3)
  420.       else if(ichoice.eq.3)then  
  421.         num_ep = zero
  422.         den_ep = one
  423.         num_mp = zero
  424.         den_mp = one
  425.         num_en = zero
  426.         den_en = one
  427.         num_mn = zero
  428.         den_mn = one
  429.         do i = 1,4
  430.         j = i-1
  431.         num_ep = num_ep + aa_ep(i)*tau**j
  432.         den_ep = den_ep + bb_ep(i)*tau**i
  433.         num_mp = num_mp + aa_mp(i)*tau**j
  434.         den_mp = den_mp + bb_mp(i)*tau**i
  435.         num_en = num_en + aa_en(i)*tau**j
  436.         den_en = den_en + bb_en(i)*tau**i
  437.         num_mn = num_mn + aa_mn(i)*tau**j
  438.         den_mn = den_mn + bb_mn(i)*tau**i
  439.         end do
  440.         G_ep = num_ep/den_ep
  441.         G_mp = mu_p*num_mp/den_mp
  442.         G_en = num_en/den_en
  443.         G_mn = mu_n*num_mn/den_mn
  444.       end if
  445. c
  446.       return
  447.       end      
  448. c
  449. c******************************************************************************
  450. c
  451. c      
  452.       subroutine sig_nuN_RES(ilept,Erec,xq,w,xk,xp,e_nu,cos_theta
  453.      $          ,costheta_p,sig)
  454. c
  455. c.....returns the cross section of the process
  456. c     nu_ell + neutron ---> ell^- + N^* in units of millibarn/sr/GeV
  457. c     all input energies and momenta are in MeV
  458. c
  459. c     The resonances included are the P33, D13, P11 and S11
  460. c
  461. c     ilept = 1 muon neutrino beam
  462. c     ilept = 2 electron neutrino beam
  463. c
  464. c******************************************************************************
  465. c
  466.       implicit none
  467. c      
  468.       integer ilept,i
  469. c      
  470.       real*8 A(5),WW(5)
  471. c
  472.       real*8 G_F,pi,c0,hc,xm_n,xm_p,xm_lept,e_lept,e_nu,xq,w,xk,xp
  473.      $      ,cos_theta,sig,xk_lept,c1,c2,c3,q2,Ek,Ep,c4,costheta_p
  474.      $      ,xkx2,xkz,xkz2,aux1,wt,q2t,summ,prod,diff,contr,ctot
  475.      $      ,xm_lept_MeV,zero,xm,xmA,xmAm1,xmR_0,EA,ER,E_min,eps
  476.      $      ,xmp,xmn,xmpi,WSQ_min,WSQ_max,WSQ,s,Ek_GeV,w_GeV,q2_GeV
  477.      $      ,xq_GeV,wt_GeV,Erec,sinqtheta_p,k_F
  478. c
  479. c     common / constants / zero,xmp,xmn,xm,xmA,xmAm1,xmR_0,xm_lept
  480. c    $                    ,EA,ER,E_min,eps      
  481.       common / constants / zero,xmp,xmn,xm,xmA,xmAm1,xmR_0,xm_lept
  482.      $                    ,EA,ER,E_min,eps,k_F
  483. c
  484.       common / RES / WSQ_min,WSQ_max,Ek_GeV      
  485. c
  486. c.....Fermi constant in fm^2 (Cabibbo's angle included)
  487. c
  488.       G_F = 1.14e-5*.1973**2
  489. c
  490.       pi = 4.*atan(1.)
  491.       c0 = 1./16./pi/pi
  492. c
  493.       hc = 197.3
  494.       xmpi = .13957             ! mass of charged pion [GeV]
  495.       xm_n = 1000.*xmn
  496.       xm_p = 1000.*xmp
  497.       xm_lept_MeV = 1000.*xm_lept
  498. c
  499.       e_lept = e_nu - w ! energy of charged lepton [MeV]
  500.       xk_lept = sqrt(e_lept**2 - xm_lept_MeV**2) ! momentum of charged lepton [MeV]
  501.       c1 = (xk_lept/e_lept)*cos_theta  
  502.       c2 = e_nu*e_lept      
  503.       c3 = xk_lept/e_nu      
  504. c
  505.       q2 = w**2 - xq**2    
  506.       w_GeV = w/1000.
  507.       xq_GeV = xq/1000.
  508.       q2_GeV = 1.e-6*q2
  509.       Ek = sqrt(xk**2 + xm_n**2) ! initial nucleon energy [MeV]
  510.       Ek_GeV = Ek/1000.
  511.       Ep = sqrt(xp**2 + xm_p**2) ! final nucleon energy [MeV]
  512.       c4 = xm_p/Ek
  513.       sinqtheta_p = 1. - costheta_p**2
  514.       xkx2 = .5*xk*xk*sinqtheta_p
  515.       xkz = xk*costheta_p
  516.       xkz2= xkz**2            
  517.       aux1 = xm_lept_MeV**2 - q2  
  518.       wt_GeV = w_GeV + xmA - Erec - Ek_GeV
  519.       wt = 1000.*wt_GeV
  520.       q2t = wt**2 - xq**2
  521.       s = xmA**2 + 2.*w_GeV*xmA + q2_GeV   ! squared CM energy [GeV^2]
  522. c      
  523.       WSQ_min = (xmp + xmpi)**2         ! squared pion production threshold [GeV^2]
  524.       WSQ = ( xm_p**2 + q2t + 2.*(wt*Ek - xq*xk*costheta_p) )/1.e6 ! squared mass of hadronic final state [GeV^2]
  525.       WSQ_max = ( sqrt(s) - xMA + xmp )**2
  526. c
  527.       summ = (w/xq)*( 1. + c1 + xm_lept_MeV**2/w/e_lept ) ! (p_nu)z/E_nu + (p_mu)z/E_mu
  528. c      
  529.       prod = (1./xq/xq)*( w*w*c1 - c2*( 1. - c1)**2    
  530.      $     + xm_lept_MeV**2*( w/e_lept + 1. - c1 ) )      ! (p_nu)z/E_nu * (p_mu)z/E_mu
  531. c      
  532.       diff = ( (e_nu + e_lept)/xq )*(1. - c1)
  533.      $     - xm_lept_MeV**2/e_lept/xq                     ! (p_nu)z/E_nu - (p_mu)z/E_mu
  534. c    
  535. c.....Alll A(i) are given in [MeV^2]. Note that A(2) A(4) nd A(5)are divided by m^2,
  536. c     A(3) is divided by 2*m^2
  537. c
  538.       A(1) = aux1/2.
  539.       A(2) = c2*( Ek**2 - Ek*xkz*summ  
  540.      $     + xkx2*(c2/xq/xq)*( (xk_lept/e_lept)**2 - c1**2 )
  541.      $     + xkz2*prod - (xm_n**2/2.)*( 1. - c1 ) )/xm_p**2
  542.       A(3) = c2*( Ek*(xkz+xq) - Ep*xkz )*diff/2./xm_p**2
  543.       A(4) = c2*(wt**2 - wt*xq*summ + xq*xq*prod
  544.      $     - 0.5*q2t*( 1. - c1) )/xm_p**2
  545.       A(5) = c2*( 2.*Ek*wt + 2.*xkz*xq*prod - ( Ek*xq + wt*xkz )*summ
  546.      $     - ( Ek*wt - xkz*xq )*( 1. - c1 ) )/xm_p**2
  547. c
  548.       call sf_RES(q2t,WSQ,WW)
  549. c
  550. c.....compute the contraction L_munu*W^munu in MeV^2/GeV and
  551. c     change units to fm^-2/GeV
  552. c
  553.       contr = zero
  554.       do 1000 i = 1,5
  555.  1000 contr = contr + 16.*A(i)*WW(i)/hc**2
  556. c
  557. c.....calculate coefficient in units of fm^4
  558. c
  559.       ctot = (4./3.)*(G_F**2/2.)*c0*c3*c4
  560. c
  561. c.....get x-section in units of mbarn/sr/GeV
  562. c
  563.       sig = 10.*ctot*contr
  564. c
  565.       return
  566.       end
  567. c      
  568. c******************************************************************************
  569. c      
  570.       subroutine sf_RES(q2,WSQ,W)
  571. c
  572.       implicit real*8(a-h,o-z)
  573. c
  574.       common / flags / ig,ilept
  575. c
  576.       common / RES / WSQ_min,WSQ_max,Ek_GeV
  577. c
  578. c.....returns the structure functions W1, W2, W3/(2m^2),
  579. c     W4 and W5 in units of GeV^-1, as a function
  580. c     of the squared four-momentum transfer q2
  581. c
  582.       dimension W(5),V1(4),V2(4),V3(4),V4(4),V5(4)
  583. c
  584.       xmP33 = 1.232           ! mass of P33 [GeV]                              
  585.       gammaP33 = 0.115        ! width of P33 [GeV]
  586.       xmD13 = 1.520           ! mass of D13 [GeV]      
  587.       gammaD13 = 0.125        ! width of D13 [GeV]  
  588.       xmP11 = 1.440           ! mass of P11 [GeV]  
  589.       gammaP11 = 0.350        ! width of P11 [GeV]
  590.       xmS11 = 1.535           ! mass of S11 [GeV]  
  591.       gammaS11 = 0.150        ! width ofSP11 [GeV]
  592. c
  593.       WW = sqrt(WSQ)
  594.       Q2g = -q2*1.e-6
  595. c
  596. c.....initialize
  597. c      
  598.       do i = 1,5      
  599.       W(i) = 0.
  600.       end do
  601. c
  602.       if(WSQ_min.ge.WSQ_max)then
  603.       return
  604.       end if      
  605. c
  606.       do i = 1,4
  607.       V1(i) = 0.
  608.       V2(i) = 0.
  609.       V3(i) = 0.
  610.       V4(i) = 0.
  611.       V5(i) = 0.
  612.       end do    
  613. c
  614. c.....P33
  615. c      
  616.       if(WSQ.le.WSQ_min.or.WSQ.ge.WSQ_max)then
  617.          return    
  618.       else    
  619.          call StructFuncP33(WW,Q2g,V1(1),V2(1),V3(1),V4(1),V5(1))
  620.          call bw(WW,xmP33,gammaP33,bw_P33)
  621. c
  622.          V1(1) = V1(1)*bw_P33/2./Ek_GeV
  623.          V2(1) = V2(1)*bw_P33/2./Ek_GeV
  624.          V3(1) = V3(1)*bw_P33/2./Ek_GeV
  625.          V4(1) = V4(1)*bw_P33/2./Ek_GeV
  626.          V5(1) = V5(1)*bw_P33/2./Ek_GeV        
  627. c
  628.          call StructFuncD13(WW,Q2g,V1(2),V2(2),V3(2),V4(2),V5(2))        
  629.          call bw(WW,xmD13,gammaD13,bw_D13)
  630. c
  631.          V1(2) = V1(2)*bw_D13/2./Ek_GeV
  632.          V2(2) = V2(2)*bw_D13/2./Ek_GeV
  633.          V3(2) = V3(2)*bw_D13/2./Ek_GeV
  634.          V4(2) = V4(2)*bw_D13/2./Ek_GeV
  635.          V5(2) = V5(2)*bw_D13/2./Ek_GeV
  636. c
  637.          call StructFuncP11(WW,Q2g,V1(3),V2(3),V3(3),V4(3),V5(3))
  638.          call bw2(WW,xmP11,gammaP11,bw_P11)
  639. c        
  640.          V1(3) = V1(3)*bw_P11/2./Ek_GeV
  641.          V2(3) = V2(3)*bw_P11/2./Ek_GeV
  642.          V3(3) = V3(3)*bw_P11/2./Ek_GeV
  643.          V4(3) = V4(3)*bw_P11/2./Ek_GeV
  644.          V5(3) = V5(3)*bw_P11/2./Ek_GeV
  645. c
  646.          call StructFuncS11(WW,Q2g,V1(4),V2(4),V3(4),V4(4),V5(4))
  647.          call bw2(WW,xmS11,gammaS11,bw_S11)
  648. c
  649.          V1(4) = V1(4)*bw_S11/2./Ek_GeV
  650.          V2(4) = V2(4)*bw_S11/2./Ek_GeV
  651.          V3(4) = V3(4)*bw_S11/2./Ek_GeV
  652.          V4(4) = V4(4)*bw_S11/2./Ek_GeV
  653.          V5(4) = V5(4)*bw_S11/2./Ek_GeV
  654. c
  655.          do i = 1,4
  656. c        do i = 2,4
  657.          W(1) = W(1) + V1(i)
  658.          W(2) = W(2) + V2(i)
  659.          W(3) = W(3) + V3(i)
  660.          W(4) = W(4) + V4(i)
  661.          W(5) = W(5) + V5(i)
  662.          end do
  663. c      
  664.       end if
  665. c      
  666.       return
  667.       end
  668. c
  669. c******************************************************************************
  670. c
  671.       subroutine StructFuncP33(xw,QQ,V1,V2,V3,V4,V5)
  672. c
  673. c.....invariant mass xw and Q^2 in GeV
  674. c.....the functions V1 and V2 are in GeV
  675. c
  676.       implicit none
  677. c      
  678.       real*8 xw,QQ,xm,xmr,pq,pqQ,pqQ2,V1,V2,V3,V4,V5,
  679.      $    P3v,P4v,P5v,P4a,P5a,P6a
  680. c
  681.       call coeffP33(QQ,P3v,P4v,P5v,P4a,P5a,P6a)
  682. c
  683.       xm = .5*(.93828 + .93957)             !GeV
  684.       xmr = 1.232                           !GeV
  685.       pq = .5*(xw**2 + QQ - xm**2)          !GeV
  686.       pqQ = QQ - pq
  687.       pqQ2 = 2.*pq - QQ
  688. c
  689.       V1 = 4.*P3v**2/xm**2/xmr**2*(pqQ**2*(pq+xm**2)+xmr**2
  690.      $     *(pq**2+QQ*xm**2+QQ*xm*xmr))
  691.      $     + 4.*(P4v*pqQ - P5v*pq)**2/xm**4*(pq + xm**2 - xm*xmr)
  692.      $     + 4.*(P3v*P4v*pqQ - P3v*P5v)/xm**3/xmr*(pqQ*
  693.      $     (pq + xm**2 - 2.*xm*xmr) - xmr**2*pq) + 4.*(pq + xm**2
  694.      $     + xm*xmr)*(P4a**2/xm**4*pqQ**2 + P5a**2
  695.      $     - 2.*P4a*P5a/xm**2*pqQ)
  696. c
  697.       V2 = 4.*P3v**2/xmr**2*QQ*(pq + xm**2 + xmr**2) + 4.*P4v**2/xm**2*
  698.      $     QQ*(pq + xm**2 - xm*xmr) + 4.*P3v*P4v/xm/xmr*QQ*(pq
  699.      $     + (xmr - xm)**2) + 4.*P3v*P5v*QQ/xm/xmr*(pq+QQ+(xmr-xm)**2) +
  700.      $     4.*P5v**2*QQ/xm**2/xmr**2*(xmr**2 + QQ)*(xm**2 - xm*xmr + pq)
  701.      $     + 8.*P4v*P5v*QQ/xm**2*(xm**2 - xm*xmr + pq)
  702.      $     + 4.*(P5a**2*xm**2/xmr**2 + P4a**2/xm**2
  703.      $     *QQ)*(pq + xm**2 + xm*xmr)
  704. c
  705.       V3 = 8.*((P4v*P4a/xm**2*pqQ - P4v*P5a)*pqQ + pq*
  706.      $     (P5v*P5a - P5v*P4a/xm**2*pqQ) + (P3v*P5a*xm/xmr -
  707.      $     P3v*P4a/xm/xmr*pqQ)*(2.*xmr**2+2.*xm*xmr+QQ-pq))
  708. c
  709.       V4 = 4.*P3v**2/xmr**2*(pqQ2*(pq + xm**2)-xmr**2*(xm**2 + xmr*xm))
  710.      $     + 4.*(P4v**2/xm**2*pqQ2 + P5v**2*pq**2/xm**2/xmr**2
  711.      $     + 2.*P4v*P5v/xm**2*pq)*(pq + xm**2 - xm*xmr)
  712.      $     + 4.*P3v*P4v/xm/xmr*(pqQ2*(pq + xm**2 - 2.*xm*xmr)+pq*xmr**2)
  713.      $     + 4.*P3v*P5v/xm/xmr*pq*(2.*pq + xm**2 - 2.*xm*xmr)
  714.      $     + 4.*(P5a**2*xm**2/xmr**2 + P4a**2/xm**2*pqQ2
  715.      $     + P6a**2/xm**2/xmr**2*(pqQ**2 + QQ*xmr**2) + 2.*P4a*P5a
  716.      $     - 2.*P4a*P6a/xm**2*pq - 2.*P5a*P6a/xmr**2*(xmr**2 + QQ - pq))
  717.      $     *(pq + xm**2 + xm*xmr)
  718. c
  719.       V5 = 4.*P3v**2/xmr**2*pq*(pq + xm**2 + xmr**2) + 4.*P3v*P5v/xm/xmr
  720.      $     *pq*(pq + (xmr - xm)**2 + QQ) + 4.*(P4v**2/xm**2
  721.      $     + P5v**2/xm**2/xmr**2*(QQ + xmr**2) + 2.*P4v*P5v/xm**2)*pq*
  722.      $     (pq + xm**2 - xm*xmr) + 4.*P3v*P4v/xm/xmr*pq*(pq
  723.      $     + (xmr-xm**2)) + 4.*(P4a**2/xm**2*pq + P5a**2*xm**2/xmr**2
  724.      $     + P4a*P5a - P4a*P6a/xm**2*QQ + P5a*P6a/xmr**2*(pq - QQ))*
  725.      $     (pq + xm**2 + xm*xmr)
  726. c
  727.       return
  728.       end
  729. c
  730. c******************************************************************************
  731. c
  732.       subroutine coeffP33(Q2,P3v,P4v,P5v,P4a,P5a,P6a)
  733. c
  734.       implicit none
  735.       real*8 Q2,xmv,xma,xm,xmpi,P3v,P4v,P5v,P4a,P5a,P6a
  736. c
  737.       xmv = .84                      !GeV
  738.       xma = 1.03                     !GeV
  739.       xm = .5*(.93828 + .93957)      !GeV
  740.       xmpi = .1396                   !GeV
  741. c
  742. c.....vector form factors
  743. c
  744.       P3v = 2.13/((1. + Q2/xmv**2)**2*(1. + Q2/4./xmv**2))
  745.       P4v = -1.51/(1. + Q2/xmv**2)**2/(1. + Q2/4./xmv**2)
  746.       P5v = 0.48/(1. + Q2/xmv**2)**2/(1. + Q2/0.776/xmv**2)
  747. c
  748. c.....axial form factors
  749. c
  750.       P5a = 1.2/(1. + Q2/xma**2)**2/(1. + Q2/3./xma**2)
  751.       P4a = -P5a/4.
  752.       P6a = P5a*xm**2/(Q2 + xmpi**2)
  753. c
  754.       return
  755.       end
  756. c
  757. c******************************************************************************
  758. c
  759.       subroutine StructFuncD13(xw,QQ,V1,V2,V3,V4,V5)
  760. c
  761. c...the invariant mass xw and Q^2 are in GeV
  762. c
  763.       implicit none
  764.       real*8 xw,QQ,xm,xmr,pq,pqQ,pqQ2,D3v,D4v,D5v,D4a,D5a,D6a
  765.      $     ,V1,V2,V3,V4,V5
  766. c
  767.       call coeffD13(QQ,D3v,D4v,D5v,D4a,D5a,D6a)
  768. c
  769.       xm = .5*(.93828 + .93957)         !GeV
  770.       xmr = 1.520                       !GeV
  771.       pq = .5*(xw**2+QQ-xm**2)          !GeV
  772.       pqQ = QQ - pq
  773.       pqQ2 = 2.*pq - QQ
  774. c
  775.       V1 = (4.*D3v**2/xm**2/xmr**2*(pqQ**2*(pq+xm**2)+xmr**2
  776.      $     *(pq**2+QQ*xm**2-QQ*xm*xmr))
  777.      $     + 4.*(D4v*pqQ - D5v*pq)**2/xm**4*(pq + xm**2 + xm*xmr)
  778.      $     + 4.*(D3v*D4v*pqQ - D3v*D5v)/xm**3/xmr*(pqQ*
  779.      $     (pq + xm**2 + 2.*xm*xmr) - xmr**2*pq) + 4.*(pq + xm**2
  780.      $     - xm*xmr)*(D4a**2/xm**4*pqQ**2 + D5a**2
  781.      $     - 2.*D4a*D5a/xm**2*pqQ))/3.
  782. c
  783.       V2 = (4.*D3v**2/xmr**2*QQ*(pq + xm**2 + xmr**2) + 4.*D4v**2/
  784.      $     xm**2*QQ*(pq + xm**2 + xm*xmr) + 4.*D3v*D4v/xm/xmr*QQ*(pq
  785.      $     +(xmr+xm)**2) + 4.*D3v*D5v*QQ/xm/xmr*(pq+QQ+(xmr+xm)**2)+
  786.      $     4.*D5v**2*QQ/xm**2/xmr**2*(xmr**2 + QQ)*(xm**2+xm*xmr + pq)
  787.      $     + 8.*D4v*D5v*QQ/xm**2*(xm**2 + xm*xmr + pq)
  788.      $     + 4.*(D5a**2*xm**2/xmr**2 + D4a**2/xm**2
  789.      $     *QQ)*(pq + xm**2 - xm*xmr))/3.
  790. c
  791.       V3 = 8.*((D4v*D4a/xm**2*pqQ - D4v*D5a)*pqQ + pq*
  792.      $     (D5v*D5a - D5v*D4a/xm**2*pqQ) + (D3v*D5a*xm/xmr -
  793.      $     D3v*D4a/xm/xmr*pqQ)*(2.*xmr**2+2.*xm*xmr+QQ-pq))/3.
  794. c
  795.       V4 = (4.*D3v**2/xmr**2*(pqQ2*(pq+xm**2)-xmr**2*(xm**2- xmr*xm))
  796.      $     + 4.*(D4v**2/xm**2*pqQ2 + D5v**2*pq**2/xm**2/xmr**2
  797.      $     + 2.*D4v*D5v/xm**2*pq)*(pq + xm**2 + xm*xmr)
  798.      $     + 4.*D3v*D4v/xm/xmr*(pqQ2*(pq+xm**2+2.*xm*xmr)+pq*xmr**2)
  799.      $     + 4.*D3v*D5v/xm/xmr*pq*(2.*pq + xm**2 + 2.*xm*xmr)
  800.      $     + 4.*(D5a**2*xm**2/xmr**2 + D4a**2/xm**2*pqQ2
  801.      $     + D6a**2/xm**2/xmr**2*(pqQ**2 + QQ*xmr**2) + 2.*D4a*D5a
  802.      $     - 2.*D4a*D6a/xm**2*pq - 2.*D5a*D6a/xmr**2*(xmr**2 + QQ - pq))
  803.      $     *(pq + xm**2 - xm*xmr))/3.
  804. c
  805.       V5 = (4.*D3v**2/xmr**2*pq*(pq + xm**2+xmr**2) + 4.*D3v*D5v
  806.      $     /xm/xmr*pq*(pq + (xmr + xm)**2 + QQ) + 4.*(D4v**2/xm**2
  807.      $     + D5v**2/xm**2/xmr**2*(QQ+xmr**2)+2.*D4v*D5v/xm**2)*pq*
  808.      $     (pq + xm**2 + xm*xmr) + 4.*D3v*D4v/xm/xmr*pq*(pq
  809.      $     + (xmr+xm**2)) + 4.*(D4a**2/xm**2*pq + D5a**2*xm**2/xmr**2
  810.      $     + D4a*D5a - D4a*D6a/xm**2*QQ + D5a*D6a/xmr**2*(pq - QQ))*
  811.      $     (pq + xm**2 - xm*xmr))/3.
  812. c
  813.       return
  814.       end
  815. c
  816. c******************************************************************************
  817. c  
  818.       subroutine StructFuncP11(xw,QQ,V1,V2,V3,V4,V5)
  819. c
  820.       implicit none
  821. c
  822.       real*8 xw,QQ,P1v,P2v,P1a,P3a,xm,xmr,pq,mu,V1,V2,V3,V4,V5
  823. c
  824.       call coeffP11(QQ,P1v,P2v,P1a,P3a)
  825. c    
  826.       xm = .5*(.93828 + .93957)         !GeV
  827.       xmr = 1.440                       !GeV
  828.       pq = .5*(xw**2 + QQ - xm**2)      !GeV
  829.       mu = xm + xmr                     !GeV
  830. c
  831.       V1 = 2.*(P1v**2/mu**4*QQ**2*(pq + xm**2 - xm*xmr) +
  832.      $     P2v**2/mu**2*(2.*pq**2 + QQ*(xm**2 + xm*xmr - pq)) +
  833.      $     P1v*P2v/mu**3*2.*QQ*(pq*(xmr - xm) + xm*QQ) +
  834.      $     P1a**2*(xm**2 + xm*xmr + pq))
  835. c
  836.       V2 = 2.*(2.*xm**2*(P1v**2/mu**4*QQ**2 + P2v**2/mu**2*QQ + P1a**2))
  837. c
  838.       V3 = 2.*(4.*xm**2*(P1v*P1a/mu**2*QQ + P2v*P1a/mu*(xmr + xm)))
  839. c
  840.       V4 = 2.*(xm**2*(P2v**2/mu**2*(pq - xm**2 - xm*xmr) +
  841.      $     P1v**2/mu**4*(2.*pq**2 - QQ*(pq + xm**2 - xm*xmr)) -
  842.      $     P1v*P2v/mu**3*(pq*(xmr - xm)+ xm*QQ) - 2.*P1a*P3a +
  843.      $     P3a**2/xm**2*(pq + xm**2 - xm*xmr)))
  844. c
  845.       V5 = 2.*(xm**2*(2.*P1v**2/mu**4*QQ*pq+2.*P2v**2/mu**2*pq + P1a**2+
  846.      $     P1a*P3a/xm*(xmr - xm)))
  847. c
  848.       return
  849.       end
  850. c
  851. c******************************************************************************
  852. c
  853.       subroutine StructFuncS11(xw,QQ,V1,V2,V3,V4,V5)
  854. c
  855.       implicit none
  856. c
  857.       real*8 xw,QQ,S1v,S2v,S1a,S3a,xm,xmr,pq,mu,V1,V2,V3,V4,V5
  858. c
  859.       call coeffS11(QQ,S1v,S2v,S1a,S3a)
  860. c    
  861.       xm = .5*(.93828 + .93957)         !GeV
  862.       xmr = 1.535                       !GeV
  863.       pq = .5*(xw**2 + QQ - xm**2)      !GeV
  864.       mu = xm + xmr                     !GeV
  865. c
  866.       V1 = 2.*(S1v**2/mu**4*QQ**2*(pq + xm**2 + xm*xmr) +
  867.      $     S2v**2/mu**2*(2.*pq**2 + QQ*(xm**2 - xm*xmr - pq)) +
  868.      $     S1v*S2v/mu**3*2.*QQ*(pq*(xmr + xm) - xm*QQ) +
  869.      $     S1a**2*(xm**2 - xmr*xm + pq))
  870. c
  871.       V2 = 2.*(2.*xm**2*(S1v**2/mu**4*QQ**2 + S2v**2/mu**2*QQ + S1a**2))
  872. c
  873.       V3 = 2.*(4.*xm**2*(S1v*S1a/mu**2*QQ + S2v*S1a/mu*(xmr - xm)))
  874. c
  875.       V4 = 2.*(xm**2*(S2v**2/mu**2*(pq - xm**2 + xm*xmr) +
  876.      $     S1v**2/mu**4*(2.*pq**2 - QQ*(pq + xm**2 + xm*xmr)) -
  877.      $     S1v*S2v/mu**3*(pq*(xmr + xm) - xm*QQ) + 2.*S1a*S3a +
  878.      $     S3a**2/xm**2*(pq + xm**2 + xm*xmr)))
  879. c
  880.       V5 = 2.*(xm**2*(2.*S1v**2/mu**4*QQ*pq + 2.*S2v**2/mu**2*pq +
  881.      $     S1a**2 + S1a*S3a/xm*(xmr + xm)))
  882. c
  883.       return
  884.       end
  885. c
  886. c******************************************************************************
  887. c
  888.       subroutine coeffD13(Q2,D3v,D4v,D5v,D4a,D5a,D6a)
  889. c
  890.       implicit none
  891.       real*8 Q2,xmv,xma,xm,xmpi,D3p_v,D4p_v,D5p_v,D3n_v,D4n_v,D5n_v
  892.      $     ,D4a,D5a,D6a,D3v,D4v,D5v
  893. c
  894.       xmv = 0.84                      !GeV
  895.       xma = 1.03                      !GeV
  896.       xm = .5*(.93828 + .93957)       !GeV
  897.       xmpi = .1396                    !GeV
  898. c
  899. c...vector form factors
  900. c
  901.       D3p_v = 2.95/(1.+ Q2/xmv**2)**2/(1. + Q2/(8.9*xmv**2))
  902.       D4p_v = -1.05/(1. + Q2/xmv**2)**2/(1. + Q2/(8.9*xmv**2))
  903.       D5p_v = -0.48/(1.+ Q2/xmv**2)**2
  904. c
  905.       D3n_v = -1.13/(1.+ Q2/xmv**2)**2/(1. + Q2/(8.9*xmv**2))
  906.       D4n_v = 0.46/(1. + Q2/xmv**2)**2/(1. + Q2/(8.9*xmv**2))
  907.       D5n_v = -0.17/(1.+ Q2/xmv**2)**2
  908. c
  909.       D3v = D3n_v - D3p_v
  910.       D4v = D4n_v - D4p_v
  911.       D5v = D5n_v - D5p_v
  912. c
  913. c...axial form factors
  914. c
  915.       D4a = 0.
  916.       D5a = -2.1/(1. + Q2/xma**2)**2/(1. + Q2/3./xma**2)
  917.       D6a = xm**2*D5a/(xmpi**2 + Q2)
  918. c
  919.       return
  920.       end
  921. c
  922. c******************************************************************************
  923. c
  924.       subroutine coeffP11(Q2,P1v,P2v,P1a,P3a)
  925. c
  926.       implicit none
  927.       real*8 Q2,xmv,xma,xm,xmr,xmpi,P1v,P2v,P1a,P3a
  928. c
  929.       xmv = 0.84                     !GeV
  930.       xma = 1.03                     !GeV
  931.       xm = .5*(.93828 + .93957)      !GeV
  932.       xmr = 1.440                    !GeV
  933.       xmpi = 0.1396                  !GeV
  934. c
  935. c...vector form factors
  936. c
  937.       P1v = -2.*2.3/(1. + Q2/xmv**2)**2/(1. + Q2/(4.3*xmv**2))
  938.       P2v = -2.*(-0.76/(1. + Q2/xmv**2)**2*(1. -  2.8*log(1. + Q2)))
  939. c
  940. c...axial form factors
  941. c
  942.       P1a = -0.51/(1. + Q2/xma**2)**2/(1. + Q2/3./xma**2)
  943.       P3a = P1a*xm*(xmr + xm)/(Q2 + xmpi**2)
  944. c
  945.       return
  946.       end
  947. c
  948. c******************************************************************************
  949. c
  950.       subroutine coeffS11(Q2,S1v,S2v,S1a,S3a)
  951. c
  952.       implicit none
  953.       real*8 Q2,xmv,xma,xm,xmr,xmpi,S1v,S2v,S1a,S3a
  954. c
  955.       xmv = 0.84                     !GeV
  956.       xma = 1.03                     !GeV
  957.       xm = .5*(.93828 + .93957)      !GeV
  958.       xmr = 1.535                    !GeV
  959.       xmpi = 0.1396                  !GeV
  960. c
  961. c...vector form factors
  962. c
  963.       S1v = -2.*(2./(1. + Q2/xmv**2)**2/(1. + Q2/(1.2*xmv**2))*
  964.      $     (1. + 7.2*log(1. + Q2)))
  965.       S2v = -2.*(0.84/(1. + Q2/xmv**2)**2*(1. + 0.11*log(1. + Q2)))
  966. c
  967. c...axial form factors
  968. c
  969.       S1a = -0.21/(1. + Q2/xma**2)**2/(1. + Q2/3./xma**2)
  970.       S3a = S1a*(xmr - xm)*xm/(Q2 + xmpi**2)
  971. c
  972.       return
  973.       end
  974. c
  975. c******************************************************************************
  976. c
  977.       subroutine bw(xw,xmr,gammar0,b_w)
  978. c
  979. c.....input values are in GeV
  980. c.....returns the Breit-Wigner function in units of 1/GeV^2
  981. c
  982.       implicit none
  983. c
  984.       real*8 pi,xw,gammar,xmr,b_w, xmpi,xmp,xmn,xm
  985.       real*8 gammar0,qPi,qPiMr,qPi0,r,fact,fact2,factor
  986. c
  987.       pi = 4.*atan(1.)
  988.       xmpi = .1396
  989.       xm = .93828
  990. c    
  991.       qPi0 = abs((xw**2-xm**2-xmpi**2)**2 - 4.*xm**2*xmpi**2)
  992.       qPi = .5/xw*sqrt(qPi0)
  993. c
  994.       qPiMr = 1./(2.*xmr)*sqrt((xmr**2-xm**2-xmpi**2)**2
  995.      $       -4.*xm**2*xmpi**2)
  996.       gammar = gammar0*(qPi/qPiMr)**3
  997.       b_w = gammar*xmr/
  998.      $     (pi*(xw**2-xmr**2)**2+pi*xmr**2*gammar**2)
  999. c
  1000.       return
  1001.       end
  1002. c      
  1003. c******************************************************************************
  1004. c
  1005.       subroutine bw2(xw,xmr,gammar0,b_w)
  1006. c
  1007. c...input values are in GeV
  1008. c...returs the Breit and Wigner in units of 1/GeV^2
  1009. c
  1010.       implicit none
  1011. c
  1012.       real*8 pi,xw,gammar,xmr,b_w, xmpi,xmp,xmn,xm
  1013.       real*8 gammar0,qPi,qPiMr,qPi0,factor
  1014. c
  1015.       pi = 4.*atan(1.)
  1016.       xmpi = .1396
  1017.       xm = .93828
  1018. c    
  1019.       qPi0 = abs((xw**2-xm**2-xmpi**2)**2 - 4.*xm**2*xmpi**2)
  1020.       qPi = .5/xw*sqrt(qPi0)
  1021. c
  1022.       qPiMr = 1./(2.*xmr)*sqrt((xmr**2-xm**2-xmpi**2)**2
  1023.      $       -4.*xm**2*xmpi**2)
  1024.       gammar = gammar0*(qPi/qPiMr)
  1025.       b_w = gammar*xmr/(pi*(xw**2-xmr**2)**2+pi*xmr**2*gammar**2)
  1026. c
  1027.       return
  1028.       end
  1029. c
  1030. c******************************************************************************
  1031. c
  1032.       function diracdelta(eps,x)
  1033. c
  1034.       implicit none      
  1035. c
  1036.       real*8 diracdelta,x,eps,pi,const
  1037. c
  1038.       pi = 4.*atan(1.)
  1039.       const =  1./2/sqrt(pi*eps)
  1040.       diracdelta = const*exp(-x*x/4./eps)
  1041. c
  1042.       return
  1043.       end
  1044. c
  1045. c******************************************************************************
  1046. c      
  1047.       subroutine sig_nuN_DIS(ilept,iflag,Erec,xq,w,xk,xp,e_nu,cos_theta
  1048.      $          ,costheta_p,sig)
  1049. c
  1050. c.....returns the cross section of the process
  1051. c     nu_ell + neutron ---> ell^- + X in the deep-inelastic sector
  1052. c     the cross section is given in units of millibarn/sr/GeV
  1053. c     all input energies and momenta are in MeV
  1054. c
  1055. c     ilept = 1 muon neutrino beam
  1056. c     ilept = 2 electron neutrino beam
  1057. c
  1058. c******************************************************************************
  1059. c
  1060.       implicit none
  1061. c      
  1062.       integer ilept,i,iflag
  1063. c      
  1064.       real*8 A(5),WW(5)
  1065. c
  1066.       real*8 G_F,pi,c0,hc,xm_n,xm_p,xm_lept,e_lept,e_nu,xq,w,xk,xp
  1067.      $      ,cos_theta,sig,xk_lept,c1,c2,c3,q2,Ek,Ep,c4,costheta_p
  1068.      $      ,xkx2,xkz,xkz2,aux1,wt,q2t,summ,prod,diff,contr,ctot
  1069.      $      ,xm_lept_MeV,zero,xm,xmA,xmAm1,xmR_0,EA,ER,E_min,eps
  1070.      $      ,xmp,xmn,xmpi,WSQ_min,WSQ_max,WSQ,s,Ek_GeV,w_GeV,q2_GeV
  1071.      $      ,xq_GeV,wt_GeV,Erec,sinqtheta_p,xm_i,k_F
  1072. c
  1073. c     common / constants / zero,xmp,xmn,xm,xmA,xmAm1,xmR_0,xm_lept
  1074. c    $                    ,EA,ER,E_min,eps      
  1075.       common / constants / zero,xmp,xmn,xm,xmA,xmAm1,xmR_0,xm_lept
  1076.      $                    ,EA,ER,E_min,eps,k_F
  1077. c
  1078. c     common / RES / WSQ_min,WSQ_max,Ek_GeV      
  1079. c
  1080. c.....Fermi constant in fm^2 (Cabibbo's angle included)
  1081. c
  1082.       G_F = 1.14e-5*.1973**2
  1083. c
  1084.       pi = 4.*atan(1.)
  1085.       c0 = 1./16./pi/pi
  1086. c
  1087.       hc = 197.3
  1088.       xmpi = .13957             ! mass of charged pion [GeV]
  1089.       xm_n = 1000.*xmn
  1090.       xm_p = 1000.*xmp
  1091.       xm_i = (xm_p + xm_n)/2.
  1092.       xm_lept_MeV = 1000.*xm_lept
  1093. c
  1094.       e_lept = e_nu - w ! energy of charged lepton [MeV]
  1095.       xk_lept = sqrt(e_lept**2 - xm_lept_MeV**2) ! momentum of charged lepton [MeV]
  1096.       c1 = (xk_lept/e_lept)*cos_theta  
  1097.       c2 = e_nu*e_lept      
  1098.       c3 = xk_lept/e_nu      
  1099. c
  1100.       q2 = w**2 - xq**2    
  1101.       w_GeV = w/1000.
  1102.       xq_GeV = xq/1000.
  1103.       q2_GeV = 1.e-6*q2
  1104.       Ek = sqrt(xk**2 + xm_n**2) ! initial nucleon energy [MeV]
  1105.       Ek_GeV = Ek/1000.
  1106.       Ep = sqrt(xp**2 + xm_p**2) ! final nucleon energy [MeV]
  1107.       c4 = xm_p/Ek
  1108.       sinqtheta_p = 1. - costheta_p**2
  1109.       xkx2 = .5*xk*xk*sinqtheta_p
  1110.       xkz = xk*costheta_p
  1111.       xkz2= xkz**2            
  1112.       aux1 = xm_lept_MeV**2 - q2  
  1113.       wt_GeV = w_GeV + xmA - Erec - Ek_GeV
  1114.       wt = 1000.*wt_GeV
  1115.       q2t = wt**2 - xq**2
  1116.       s = xmA**2 + 2.*w_GeV*xmA + q2_GeV   ! squared CM energy [GeV^2]
  1117. c      
  1118.       WSQ_min = (xmp + 2.*xmpi)**2         ! squared pion production threshold [GeV^2]
  1119.       WSQ = ( xm_p**2 + q2t + 2.*(wt*Ek - xq*xk*costheta_p) )/1.e6 ! squared mass of hadronic final state [GeV^2]
  1120.       WSQ_max = ( sqrt(s) - xMA + xmp )**2
  1121. c
  1122.       summ = (w/xq)*( 1. + c1 + xm_lept_MeV**2/w/e_lept ) ! (p_nu)z/E_nu + (p_mu)z/E_mu
  1123. c      
  1124.       prod = (1./xq/xq)*( w*w*c1 - c2*( 1. - c1)**2    
  1125.      $     + xm_lept_MeV**2*( w/e_lept + 1. - c1 ) )      ! (p_nu)z/E_nu * (p_mu)z/E_mu
  1126. c      
  1127.       diff = ( (e_nu + e_lept)/xq )*(1. - c1)
  1128.      $     - xm_lept_MeV**2/e_lept/xq                     ! (p_nu)z/E_nu - (p_mu)z/E_mu
  1129. c    
  1130. c.....compute A(1) and A(i>1)/m^2 [MeV^2]
  1131. c
  1132.       A(1) = aux1/2.
  1133.       A(2) = c2*( Ek**2 - Ek*xkz*summ  
  1134.      $     + xkx2*(c2/xq/xq)*( (xk_lept/e_lept)**2 - c1**2 )
  1135.      $     + xkz2*prod - (xm_n**2/2.)*( 1. - c1 ) )/xm_i**2
  1136.       A(3) = c2*( Ek*(xkz+xq) - Ep*xkz )*diff/xm_i**2
  1137.       A(4) = c2*(wt**2 - wt*xq*summ + xq*xq*prod
  1138.      $     - 0.5*q2t*( 1. - c1) )/xm_i**2
  1139.       A(5) = c2*( 2.*Ek*wt + 2.*xkz*xq*prod - ( Ek*xq + wt*xkz )*summ
  1140.      $     - ( Ek*wt - xkz*xq )*( 1. - c1 ) )/xm_i**2
  1141. c
  1142.       call sf_DIS(iflag,q2t,WSQ,WSQ_min,WW)
  1143. c
  1144. c.....compute the contraction L_munu*W^munu in MeV^2/GeV and
  1145. c     change units to fm^-2/GeV
  1146. c
  1147.       contr = zero
  1148.       do 1000 i = 1,5
  1149.  1000 contr = contr + 16.*A(i)*WW(i)/hc**2
  1150. c
  1151. c.....calculate coefficient in units of fm^4
  1152. c
  1153.       ctot = (G_F**2/2.)*c0*c3*c4
  1154. c
  1155. c.....get x-section in units of mbarn/sr/GeV
  1156. c
  1157.       sig = 10.*ctot*contr
  1158. c
  1159.       return
  1160.       end
  1161. c      
  1162. c******************************************************************************
  1163. c      
  1164.       subroutine sf_DIS(iflg,q2,WSQ,WSQ_min,W)
  1165. c
  1166.       implicit real*8 (a-h,o-z)
  1167.       integer iflg
  1168. c
  1169.       dimension W(5)      
  1170. c
  1171. c.....INTERPOLATION PROGRAM INITIALIZATION PARAMETER
  1172. c      
  1173.       COMMON / INTINIP / IINIP
  1174.       COMMON / INTINIF / IINIF
  1175. c
  1176.       Q2g_min = 0.8      
  1177. c
  1178. c.....initialization      
  1179. c
  1180.       do i = 1,5
  1181.       W(i) = 0.
  1182.       end do
  1183. c
  1184. c.....SELECTION OF THE SET OF PARTON DISTRIBUTIONS
  1185. c      
  1186.       ISET = 2
  1187. c
  1188.       if(iflg.eq.0)then      
  1189.       IINIP = 0
  1190.       IINIF = 0
  1191.       end if
  1192. c
  1193.       pi = 4.*atan(1.)
  1194.       amp = .93818
  1195.       amn = .93957
  1196.       am2 = ( (amp + amn)/2. )**2
  1197.       R = 0.18
  1198. c
  1199. c.....CALCULATE STRUCTURE FUNCTIONS FOR AN ISOSCALAR NUCLEON    
  1200. c
  1201.       w1_N = 0.
  1202.       w2_N = 0.
  1203.       w3_N = 0.
  1204. c      
  1205.       if (WSQ.LE.WSQ_min) return
  1206. c      
  1207.       Q2g = -q2*1.e-6
  1208.       xnum = Q2g
  1209.       xden = WSQ - am2 + Q2g
  1210.       x = xnum/xden
  1211. c     v = (WSQ - am2 + Q2g)/2./sqrt(am2)
  1212. c     vsq = v*v
  1213. c
  1214.       if(Q2g.le.Q2g_min)then
  1215.       Q2g_new = Q2g_min
  1216.       else
  1217.       Q2g_new = Q2g
  1218.       end if
  1219. c      
  1220.       CALL GRV98PA (ISET, X, Q2g_new, UV, DV, US, DS, SS, GL)
  1221. c
  1222.       F2_N = ( UV + 2.*US ) + ( DV + 2.*DS )
  1223.       F3_N = ( UV + DV )/X
  1224. c          
  1225. c     v = (WSQ - am2 + Q2g)/2./sqrt(am2)
  1226.       v = (WSQ - am2 + Q2g_new)/2./sqrt(am2)
  1227.       vsq = v*v
  1228.       w2_N = F2_N/v
  1229.       w1_N = w2_N*(1.0+vsq/Q2g)/(1.0+R)
  1230. c     w1_N = const*w2_N
  1231.       w3_N = F3_N/v
  1232. c                                            
  1233.       W(1) = w1_N
  1234.       W(2) = w2_N
  1235.       W(3) = w3_N
  1236.       W(4) = 0.
  1237.       W(5) = 0.
  1238. c      
  1239.       return
  1240.       end
  1241. c        
  1242. c******************************************************************************
  1243. * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  1244. *                                                                   *
  1245. *     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      *
  1246. *                                                                   *
  1247. *                          1998 UPDATE                              *
  1248. *                                                                   *
  1249. *                  For a detailed explanation see                   *
  1250. *                   M. Glueck, E. Reya, A. Vogt :                   *
  1251. *        hep-ph/9806404  =  DO-TH 98/07  =  WUE-ITP-98-019          *
  1252. *                  (To appear in Eur. Phys. J. C)                   *
  1253. *                                                                   *
  1254. *   This package contains subroutines returning the light-parton    *
  1255. *   distributions in NLO (for the MSbar and DIS schemes) and LO;    *
  1256. *   the respective light-parton, charm, and bottom contributions    *
  1257. *   to F2(electromagnetic); and the scale dependence of alpha_s.    *
  1258. *                                                                   *
  1259. *   The parton densities and F2 values are calculated from inter-   *
  1260. *   polation grids covering the regions                             *
  1261. *         Q^2/GeV^2  between   0.8   and  1.E6 ( 1.E4 for F2 )      *
  1262. *            x       between  1.E-9  and   1.                       *
  1263. *   Any call outside these regions stops the program execution.     *
  1264. *                                                                   *
  1265. *   At Q^2 = MZ^2, alpha_s reads  0.114 (0.125) in NLO (LO); the    *
  1266. *   heavy quark thresholds, QH^2 = mh^2, in the beta function are   *
  1267. *            mc = 1.4 GeV,  mb = 4.5 GeV,  mt = 175 GeV.            *
  1268. *   Note that the NLO alpha_s running is different from GRV(94).    *
  1269. *                                                                   *
  1270. *    Questions, comments etc to:  avogt@physik.uni-wuerzburg.de     *
  1271. *                                                                   *
  1272. * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
  1273. *
  1274. *
  1275. *
  1276. *
  1277.       SUBROUTINE GRV98PA (ISET, X, Q2, UV, DV, US, DS, SS, GL)
  1278. *********************************************************************
  1279. *                                                                   *
  1280. *   THE PARTON ROUTINE.                                             *
  1281. *                                     __                            *
  1282. *   INPUT:   ISET =  1 (LO),  2 (NLO, MS), or  3 (NLO, DIS)         *
  1283. *            X  =  Bjorken-x        (between  1.E-9 and 1.)         *
  1284. *            Q2 =  scale in GeV**2  (between  0.8 and 1.E6)         *
  1285. *                                                                   *
  1286. *   OUTPUT:  UV = u - u(bar),  DV = d - d(bar),  US = u(bar),       *
  1287. *            DS = d(bar),  SS = s = s(bar),  GL = gluon.            *
  1288. *            Always x times the distribution is returned.           *
  1289. *                                                                   *
  1290. *   COMMON:  The main program or the calling routine has to have    *
  1291. *            a common block  COMMON / INTINIP / IINIP , and the     *
  1292. *            integer variable  IINIP  has always to be zero when    *
  1293. *            GRV98PA is called for the first time or when  ISET     *
  1294. *            has been changed.                                      *
  1295. *                                                                   *
  1296. *   GRIDS:   1. grv98lo.grid, 2. grv98nlm.grid, 3. grv98nld.grid,   *
  1297. *            (1+1809 lines with 6 columns, 4 significant figures)   *
  1298. *                                                                   *
  1299. *******************************************************i*************
  1300. *
  1301.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1302.       PARAMETER (NPART=6, NX=68, NQ=27, NARG=2)
  1303. C     PARAMETER (NPART=6, NX=68, NQ=27, NARG=5)
  1304.       DIMENSION XUVF(NX,NQ), XDVF(NX,NQ), XDEF(NX,NQ), XUDF(NX,NQ),
  1305.      1          XSF(NX,NQ), XGF(NX,NQ), PARTON (NPART,NQ,NX-1),
  1306.      2          QS(NQ), XB(NX), XT(NARG), NA(NARG), ARRF(NX+NQ)
  1307.       CHARACTER*80 LINE
  1308.       COMMON / INTINIP / IINIP
  1309.       SAVE XUVF, XDVF, XDEF, XUDF, XSF, XGF, NA, ARRF
  1310. *
  1311. *...BJORKEN-X AND Q**2 VALUES OF THE GRID :
  1312.        DATA QS / 0.8E0,
  1313.      1           1.0E0, 1.3E0, 1.8E0, 2.7E0, 4.0E0, 6.4E0,
  1314.      2           1.0E1, 1.6E1, 2.5E1, 4.0E1, 6.4E1,
  1315.      3           1.0E2, 1.8E2, 3.2E2, 5.7E2,
  1316.      4           1.0E3, 1.8E3, 3.2E3, 5.7E3,
  1317.      5           1.0E4, 2.2E4, 4.6E4,
  1318.      6           1.0E5, 2.2E5, 4.6E5,
  1319.      7           1.E6 /
  1320.        DATA XB / 1.0E-9, 1.8E-9, 3.2E-9, 5.7E-9,
  1321.      A           1.0E-8, 1.8E-8, 3.2E-8, 5.7E-8,
  1322.      B           1.0E-7, 1.8E-7, 3.2E-7, 5.7E-7,
  1323.      C           1.0E-6, 1.4E-6, 2.0E-6, 3.0E-6, 4.5E-6, 6.7E-6,
  1324.      1           1.0E-5, 1.4E-5, 2.0E-5, 3.0E-5, 4.5E-5, 6.7E-5,
  1325.      2           1.0E-4, 1.4E-4, 2.0E-4, 3.0E-4, 4.5E-4, 6.7E-4,
  1326.      3           1.0E-3, 1.4E-3, 2.0E-3, 3.0E-3, 4.5E-3, 6.7E-3,
  1327.      4           1.0E-2, 1.4E-2, 2.0E-2, 3.0E-2, 4.5E-2, 0.06, 0.08,
  1328.      5           0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275,
  1329.      6           0.3, 0.325, 0.35, 0.375, 0.4,  0.45, 0.5, 0.55,
  1330.      7           0.6, 0.65,  0.7,  0.75,  0.8,  0.85, 0.9, 0.95, 1. /
  1331. *
  1332. *...CHECK OF X AND Q2 VALUES :
  1333.       IF ( (X.LT.0.99D-9) .OR. (X.GT.1.D0) ) THEN
  1334.          WRITE(6,91)
  1335.   91     FORMAT (2X,'PARTON INTERPOLATION: X OUT OF RANGE')
  1336.          STOP
  1337.       ENDIF
  1338.       IF ( (Q2.LT.0.799) .OR. (Q2.GT.1.01E6) ) THEN
  1339.          WRITE(6,92)
  1340.   92     FORMAT (2X,'PARTON INTERPOLATION: Q2 OUT OF RANGE')
  1341.          STOP
  1342.       ENDIF
  1343.       IF (IINIP .NE. 0) GOTO 16
  1344. *
  1345. *...INITIALIZATION, IF REQUIRED :
  1346. *
  1347. *    SELECTION AND READING OF THE GRID :
  1348. *    (COMMENT: FIRST NUMBER IN THE FIRST LINE OF THE GRID)
  1349.       IF (ISET .EQ. 1) THEN
  1350.         OPEN (11,FILE='grv98_grids/grv98lo.grid',STATUS='old')   ! 7.332E-05
  1351.       ELSE IF (ISET .EQ. 2) THEN
  1352.         OPEN (11,FILE='grv98_grids/grv98nlm.grid',STATUS='old')  ! 1.015E-04
  1353.       ELSE IF (ISET .EQ. 3) THEN
  1354.         OPEN (11,FILE='grv98_grids/grv98nld.grid',STATUS='old')  ! 1.238E-04
  1355.       ELSE
  1356.         WRITE(6,93)
  1357.   93    FORMAT (2X,'NO OR INVALID PARTON SET CHOICE')
  1358.         STOP
  1359.       END IF
  1360.       IINIP = 1
  1361.       READ(11,89) LINE
  1362.   89  FORMAT(A80)
  1363.       DO 15 M = 1, NX-1
  1364.       DO 15 N = 1, NQ
  1365.       READ(11,90) PARTON(1,N,M), PARTON(2,N,M), PARTON(3,N,M),
  1366.      1            PARTON(4,N,M), PARTON(5,N,M), PARTON(6,N,M)
  1367.   90  FORMAT (6(1PE10.3))
  1368.   15  CONTINUE
  1369.       CLOSE(11)
  1370. *
  1371. *....ARRAYS FOR THE INTERPOLATION SUBROUTINE :
  1372.       DO 10 IQ = 1, NQ
  1373.       DO 20 IX = 1, NX-1
  1374.         XB0V = XB(IX)**0.5
  1375.         XB0S = XB(IX)**(-0.2)
  1376.         XB1 = 1.-XB(IX)
  1377.         XUVF(IX,IQ) = PARTON(1,IQ,IX) / (XB1**3 * XB0V)
  1378.         XDVF(IX,IQ) = PARTON(2,IQ,IX) / (XB1**4 * XB0V)
  1379.         XDEF(IX,IQ) = PARTON(3,IQ,IX) / (XB1**7 * XB0V)
  1380.         XUDF(IX,IQ) = PARTON(4,IQ,IX) / (XB1**7 * XB0S)
  1381.         XSF(IX,IQ)  = PARTON(5,IQ,IX) / (XB1**7 * XB0S)
  1382.         XGF(IX,IQ)  = PARTON(6,IQ,IX) / (XB1**5 * XB0S)
  1383.   20  CONTINUE
  1384.         XUVF(NX,IQ) = 0.E0
  1385.         XDVF(NX,IQ) = 0.E0
  1386.         XDEF(NX,IQ) = 0.E0
  1387.         XUDF(NX,IQ) = 0.E0
  1388.         XSF(NX,IQ)  = 0.E0
  1389.         XGF(NX,IQ)  = 0.E0
  1390.   10  CONTINUE  
  1391.       NA(1) = NX
  1392.       NA(2) = NQ
  1393.       DO 30 IX = 1, NX
  1394.         ARRF(IX) = DLOG(XB(IX))
  1395.   30  CONTINUE
  1396.       DO 40 IQ = 1, NQ
  1397.         ARRF(NX+IQ) = DLOG(QS(IQ))
  1398.   40  CONTINUE
  1399. *
  1400. *...CONTINUATION, IF INITIALIZATION WAS DONE PREVIOUSLY.
  1401. *
  1402.   16  CONTINUE
  1403. *
  1404. *...INTERPOLATION :
  1405.       XT(1) = DLOG(X)
  1406.       XT(2) = DLOG(Q2)
  1407.       X1 = 1.- X
  1408.       XV = X**0.5
  1409.       XS = X**(-0.2)
  1410.       UV = FINT(NARG,XT,NA,ARRF,XUVF) * X1**3 * XV
  1411.       DV = FINT(NARG,XT,NA,ARRF,XDVF) * X1**4 * XV
  1412.       DE = FINT(NARG,XT,NA,ARRF,XDEF) * X1**7 * XV
  1413.       UD = FINT(NARG,XT,NA,ARRF,XUDF) * X1**7 * XS
  1414.       US = 0.5 * (UD - DE)
  1415.       DS = 0.5 * (UD + DE)
  1416.       SS = FINT(NARG,XT,NA,ARRF,XSF)  * X1**7 * XS
  1417.       GL = FINT(NARG,XT,NA,ARRF,XGF)  * X1**5 * XS
  1418. *
  1419.  60   RETURN
  1420.       END
  1421. *
  1422. *
  1423. *
  1424. *
  1425.       SUBROUTINE GRV98F2 (ISET, X, Q2, F2L, F2C, F2B, F2)
  1426. *********************************************************************
  1427. *                                                                   *
  1428. *   THE F2 ROUTINE.                                                 *
  1429. *                                                                   *
  1430. *   INPUT :  ISET = 1 (LO),  2 (NLO).                               *
  1431. *            X  =  Bjorken-x        (between  1.E-9 and 1.)         *
  1432. *            Q2 =  scale in GeV**2  (between  0.8 and 1.E4)         *
  1433. *                                                                   *
  1434. *   OUTPUT:  F2L = F2(light), F2C = F2(charm), F2B = F2(bottom,)    *
  1435. *            F2  = sum, all given for electromagnetic proton DIS.   *
  1436. *                                                                   *
  1437. *   COMMON:  The main program or the calling routine has to have    *
  1438. *            a common block  COMMON / INTINIF / IINIF , and the     *
  1439. *            integer variable  IINIF  has always to be zero when    *
  1440. *            GRV98F2 is called for the first time or when  ISET     *
  1441. *            has been changed.                                      *
  1442. *                                                                   *
  1443. *   GRIDS:   1. grv98lof.grid, 2. grv98nlf.grid.                    *
  1444. *            (1+1407 lines with 3 columns, 4 significant figures)   *
  1445. *                                                                   *
  1446. *******************************************************i*************
  1447. *
  1448.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1449.       PARAMETER (NSTRF=3, NX=68, NQ=21, NARG=2)
  1450. C     PARAMETER (NSTRF=3, NX=68, NQ=21, NARG=5)
  1451.       DIMENSION XF2LF(NX,NQ), XF2CF(NX,NQ), XF2BF(NX,NQ),
  1452.      1          STRFCT (NSTRF,NQ,NX-1), QS(NQ), XB(NX),
  1453.      3          XT(NARG), NA(NARG), ARRF(NX+NQ)
  1454.       CHARACTER*80 LINE
  1455.       COMMON / INTINIF / IINIF
  1456.       SAVE XF2LF, XF2CF, XF2BF, NA, ARRF
  1457. *
  1458. *...BJORKEN-X AND Q**2 VALUES OF THE GRID :
  1459.        DATA QS / 0.8E0,
  1460.      1           1.0E0, 1.3E0, 1.8E0, 2.7E0, 4.0E0, 6.4E0,
  1461.      2           1.0E1, 1.6E1, 2.5E1, 4.0E1, 6.4E1,
  1462.      3           1.0E2, 1.8E2, 3.2E2, 5.7E2,
  1463.      4           1.0E3, 1.8E3, 3.2E3, 5.7E3,
  1464.      5           1.0E4/
  1465.        DATA XB / 1.0E-9, 1.8E-9, 3.2E-9, 5.7E-9,
  1466.      A           1.0E-8, 1.8E-8, 3.2E-8, 5.7E-8,
  1467.      B           1.0E-7, 1.8E-7, 3.2E-7, 5.7E-7,
  1468.      C           1.0E-6, 1.4E-6, 2.0E-6, 3.0E-6, 4.5E-6, 6.7E-6,
  1469.      1           1.0E-5, 1.4E-5, 2.0E-5, 3.0E-5, 4.5E-5, 6.7E-5,
  1470.      2           1.0E-4, 1.4E-4, 2.0E-4, 3.0E-4, 4.5E-4, 6.7E-4,
  1471.      3           1.0E-3, 1.4E-3, 2.0E-3, 3.0E-3, 4.5E-3, 6.7E-3,
  1472.      4           1.0E-2, 1.4E-2, 2.0E-2, 3.0E-2, 4.5E-2, 0.06, 0.08,
  1473.      5           0.1, 0.125, 0.15, 0.175, 0.2, 0.225, 0.25, 0.275,
  1474.      6           0.3, 0.325, 0.35, 0.375, 0.4,  0.45, 0.5, 0.55,
  1475.      7           0.6, 0.65,  0.7,  0.75,  0.8,  0.85, 0.9, 0.95, 1. /
  1476. *
  1477. *...CHECK OF X AND Q2 VALUES :
  1478.       IF ( (X.LT.0.99D-9) .OR. (X.GT.1.D0) ) THEN
  1479.          WRITE(6,91)
  1480.   91     FORMAT (2X,'STR.FCT. INTERPOLATION: X OUT OF RANGE')
  1481.          STOP
  1482.       ENDIF
  1483. C     IF ( (Q2.LT.0.799) .OR. (Q2.GT.1.01E4) ) THEN
  1484. C        WRITE(6,92)
  1485. C 92     FORMAT (2X,'STR.FCT. INTERPOLATION: Q2 OUT OF RANGE')
  1486. C        STOP
  1487.       IF ( Q2.LT.0.799 ) THEN
  1488.       Q2 = 0.799        
  1489.       ENDIF
  1490.       IF (IINIF .NE. 0) GOTO 16
  1491. *
  1492. *...INITIALIZATION, IF REQUIRED :
  1493. *
  1494. *    SELECTION AND READING OF THE GRID :
  1495. *    (COMMENT: FIRST NUMBER IN THE FIRST LINE OF THE GRID)
  1496.       IF (ISET .EQ. 1) THEN
  1497.         OPEN (11,FILE='grv98_grids/grv98lof.grid',STATUS='old')  !  7.907E-01
  1498.       ELSE IF (ISET.EQ.2) THEN
  1499.         OPEN (11,FILE='grv98_grids/grv98nlf.grid',STATUS='old')  !  9.368E-01
  1500.       ELSE
  1501.         WRITE(6,93)
  1502.   93    FORMAT (2X,'NO OR INVALID STR.FCT. SET CHOICE')
  1503.         STOP
  1504.       END IF
  1505.       IINIF = 1
  1506.       READ(11,89) LINE
  1507.   89  FORMAT(A80)
  1508.       DO 15 M = 1, NX-1
  1509.       DO 15 N = 1, NQ
  1510.       READ(11,90) STRFCT(1,N,M), STRFCT(2,N,M), STRFCT(3,N,M)
  1511.   90  FORMAT (3(1PE10.3))
  1512.   15  CONTINUE
  1513.       CLOSE(11)
  1514. *
  1515. *....ARRAYS FOR THE INTERPOLATION SUBROUTINE :
  1516.       DO 10 IQ = 1, NQ
  1517.       DO 20 IX = 1, NX-1
  1518.         XBS = XB(IX)**0.2
  1519.         XB1 = 1.-XB(IX)
  1520.         XF2LF(IX,IQ) = STRFCT(1,IQ,IX) / XB1**3 * XBS
  1521.         XF2CF(IX,IQ) = STRFCT(2,IQ,IX) / XB1**7 * XBS
  1522.         XF2BF(IX,IQ) = STRFCT(3,IQ,IX) / XB1**7 * XBS
  1523.   20  CONTINUE
  1524.         XF2LF(NX,IQ) = 0.E0
  1525.         XF2CF(NX,IQ) = 0.E0
  1526.         XF2BF(NX,IQ) = 0.E0
  1527.   10  CONTINUE  
  1528.       NA(1) = NX
  1529.       NA(2) = NQ
  1530.       DO 30 IX = 1, NX
  1531.         ARRF(IX) = DLOG(XB(IX))
  1532.   30  CONTINUE
  1533.       DO 40 IQ = 1, NQ
  1534.         ARRF(NX+IQ) = DLOG(QS(IQ))
  1535.   40  CONTINUE
  1536. *
  1537. *...CONTINUATION, IF INITIALIZATION WAS DONE PREVIOUSLY.
  1538. *
  1539.   16  CONTINUE
  1540. *
  1541. *...INTERPOLATION :
  1542.       XT(1) = DLOG(X)
  1543.       XT(2) = DLOG(Q2)
  1544.       X1 = 1.- X
  1545.       XS = X**(-0.2)
  1546.       F2L = FINT(NARG,XT,NA,ARRF,XF2LF) * X1**3 * XS
  1547.       F2C = FINT(NARG,XT,NA,ARRF,XF2CF) * X1**7 * XS
  1548.       F2B = FINT(NARG,XT,NA,ARRF,XF2BF) * X1**7 * XS
  1549.       F2  = F2L + F2C + F2B
  1550. *
  1551.  60   RETURN
  1552.       END
  1553. *
  1554. *
  1555. *
  1556. *
  1557.       FUNCTION FINT(NARG,ARG,NENT,ENT,TABLE)
  1558. *********************************************************************
  1559. *                                                                   *
  1560. *   THE INTERPOLATION ROUTINE (CERN LIBRARY ROUTINE E104)           *
  1561. *                                                                   *
  1562. *********************************************************************
  1563.       IMPLICIT DOUBLE PRECISION (A-H, O-Z)
  1564. C     DIMENSION ARG(5),NENT(5),ENT(10),TABLE(10)
  1565.       DIMENSION ARG(2),NENT(2),ENT(10),TABLE(10)
  1566.       DIMENSION D(5),NCOMB(5),IENT(5)
  1567.       KD=1
  1568.       M=1
  1569.       JA=1
  1570.          DO 5 I=1,NARG
  1571.       NCOMB(I)=1
  1572.       JB=JA-1+NENT(I)
  1573.          DO 2 J=JA,JB
  1574.       IF (ARG(I).LE.ENT(J)) GO TO 3
  1575.     2 CONTINUE
  1576.       J=JB
  1577.     3 IF (J.NE.JA) GO TO 4
  1578.       J=J+1
  1579.     4 JR=J-1
  1580.       D(I)=(ENT(J)-ARG(I))/(ENT(J)-ENT(JR))
  1581.       IENT(I)=J-JA
  1582.       KD=KD+IENT(I)*M
  1583.       M=M*NENT(I)
  1584.     5 JA=JB+1
  1585.       FINT=0.
  1586.    10 FAC=1.
  1587.       IADR=KD
  1588.       IFADR=1
  1589.          DO 15 I=1,NARG
  1590.       IF (NCOMB(I).EQ.0) GO TO 12
  1591.       FAC=FAC*(1.-D(I))
  1592.       GO TO 15
  1593.    12 FAC=FAC*D(I)
  1594.       IADR=IADR-IFADR
  1595.    15 IFADR=IFADR*NENT(I)
  1596.       FINT=FINT+FAC*TABLE(IADR)
  1597.       IL=NARG
  1598.    40 IF (NCOMB(IL).EQ.0) GO TO 80
  1599.       NCOMB(IL)=0
  1600.       IF (IL.EQ.NARG) GO TO 10
  1601.       IL=IL+1
  1602.          DO 50  K=IL,NARG
  1603.    50 NCOMB(K)=1
  1604.       GO TO 10
  1605.    80 IL=IL-1
  1606.       IF(IL.NE.0) GO TO 40
  1607.       RETURN
  1608.       END
  1609. *
  1610. *
  1611. *
  1612. *
  1613.       FUNCTION ALPHAS (Q2, NAORD)
  1614. *********************************************************************
  1615. *                                                                   *
  1616. *   THE ALPHA_S ROUTINE.                                            *
  1617. *                                                                   *
  1618. *   INPUT :  Q2    =  scale in GeV**2  (not too low, of course);    *
  1619. *            NAORD =  1 (LO),  2 (NLO).                             *
  1620. *                                                                   *
  1621. *   OUTPUT:  alphas_s/(4 pi) for use with the GRV(98) partons.      *  
  1622. *                                                                   *
  1623. *******************************************************i*************
  1624. *
  1625.       IMPLICIT DOUBLE PRECISION (A - Z)
  1626.       INTEGER NF, K, I, NAORD
  1627.       DIMENSION LAMBDAL (3:6),  LAMBDAN (3:6), Q2THR (3)
  1628. *
  1629. *...HEAVY QUARK THRESHOLDS AND LAMBDA VALUES :
  1630.       DATA Q2THR   /  1.960,  20.25,  30625. /
  1631.       DATA LAMBDAL / 0.2041, 0.1750, 0.1320, 0.0665 /
  1632.       DATA LAMBDAN / 0.2994, 0.2460, 0.1677, 0.0678 /
  1633. *
  1634. *...DETERMINATION OF THE APPROPRIATE NUMBER OF FLAVOURS :
  1635.       NF = 3
  1636.       DO 10 K = 1, 3
  1637.       IF (Q2 .GT. Q2THR (K)) THEN
  1638.          NF = NF + 1
  1639.       ELSE
  1640.           GO TO 20
  1641.        END IF
  1642.   10   CONTINUE
  1643. *
  1644. *...LO ALPHA_S AND BETA FUNCTION FOR NLO CALCULATION :
  1645.   20   B0 = 11.- 2./3.* NF
  1646.        B1 = 102.- 38./3.* NF
  1647.        B10 = B1 / (B0*B0)
  1648.        IF (NAORD .EQ. 1) THEN
  1649.          LAM2 = LAMBDAL (NF) * LAMBDAL (NF)
  1650.          ALP  = 1./(B0 * DLOG (Q2/LAM2))
  1651.          GO TO 1
  1652.        ELSE IF (NAORD .EQ. 2) then
  1653.          LAM2 = LAMBDAN (NF) * LAMBDAN (NF)
  1654.          B1 = 102.- 38./3.* NF
  1655.          B10 = B1 / (B0*B0)
  1656.        ELSE
  1657.          WRITE (6,91)
  1658.   91     FORMAT ('INVALID CHOICE FOR ORDER IN ALPHA_S')
  1659.          STOP
  1660.        END IF
  1661. *
  1662. *...START VALUE FOR NLO ITERATION :
  1663.        LQ2 = DLOG (Q2 / LAM2)
  1664.        ALP = 1./(B0*LQ2) * (1.- B10*DLOG(LQ2)/LQ2)
  1665. *
  1666. *...EXACT NLO VALUE, FOUND VIA NEWTON PROCEDURE :
  1667.        DO 2 I = 1, 6
  1668.        XL  = DLOG (1./(B0*ALP) + B10)
  1669.        XLP = DLOG (1./(B0*ALP*1.01) + B10)
  1670.        XLM = DLOG (1./(B0*ALP*0.99) + B10)
  1671.        Y  = LQ2 - 1./ (B0*ALP) + B10 * XL
  1672.        Y1 = (- 1./ (B0*ALP*1.01) + B10 * XLP
  1673.      1       + 1./ (B0*ALP*0.99) - B10 * XLP) / (0.02D0*ALP)
  1674.        ALP = ALP - Y/Y1
  1675.   2    CONTINUE
  1676. *
  1677. *...OUTPUT :
  1678.   1    ALPHAS = ALP
  1679.        RETURN
  1680.        END
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!