(* :Title: NuclearConvolution *) (* :Summary: Generalized convolution integral for nuclear DIS structure functions. *) (* :Author: Sergey Kulagin *) (* :Package Version: Mar 12, 2006 *) (* :Sources: Packages Convolution2D.m Convolution3D.m For the discussion of the used approach see S.A. Kulagin and R. Petti, Nucl.Phys. A765 (2006) 126-187 [hep-ph/0412425] *) (* :Comments: Multiplication of the integrand by a large number caux and subsequent division of the result by caux was introduced to avoid a numerical problem with small numbers at large x calculations. *) BeginPackage["NuclearConvolution`"] prtcted = Unprotect[GenNConv2,GenNConv3,NConv2] GenNConv2::usage= "GenNConv2[n, MD, e0, m0, xBj, qsq, SF, SF2] is the convolution integral of a simplified version of nuclear spectral function, S(e,p)=2*Pi*MD(p)*delta(e-e0+p^2/(2m0)) with MD(p) the nuclear momentum distribution, and the nucleon structure function SF. The integral is calculated in terms of (p, cos(theta)) variables using the condition (p+q)^2 > M^2. The parameters xBj = the Bjorken x and qsq = Q^2 in GeV^2. The first argument n=1,2,3 defines the type of the structure function, i.e. xF1, F2, or xF3. SF is the bound nucleon structure function, which must be defined as a function of 3 variables: x, qsq, and the nucleon virtuality psq = (M+e)^2 - p^2. The last argument SF2 is optional and is only needed for n=1, in this case SF2 must be set to the nucleon structure function F2." NConv2::usage= "NConv2[xpdf, MD, e0, m0, xBj, qsq] is the convolution integral of the parton distribution xpdf with a simplified version of nuclear spectral function S(e,p)=2*Pi*MD(p)*delta(e-e0+p^2/(2m0)), where MD(p) is nuclear momentum distribution. The integral is calculated in terms of (p, cos(theta)) variables using the condition (p+q)^2 > M^2. The parameters xBj = Bjorken x and qsq = Q^2 in GeV^2. The xpdf is the parton distribution in a bound nucleon, which must be defined as a function of 3 variables: x, q2=qsq/hc2, and the nucleon virtuality psq = (M+e)^2 - p^2. Note that q2 and psq are in Fm units." GenNConv3::usage= "GenNConv3[n, S, e0, m0, xBj, qsq, SF, SF2] is the convolution integral of nuclear spectral function S(e,p,cos(theta)) and the nucleon structure function SF. The convolution is calculated in terms of (e,cos(theta),p) variables and taking into account the condition W^2 > M^2. The parameters xBj and qsq are Bjorken x and Q^2 in GeV^2. The parameters e0 and m0 determine the threshold value of nucleon removal energy as e0-p^2/(2 m0). The first argument (n=1,2,3) defines the type of the structure function, i.e. xF1, F2 or xF3. SF is the bound nucleon structure function, which must be defined as a function of 3 variables: x, qsq, and the nucleon virtuality psq = (M+e)^2 - p^2. The last argument SF2 is optional and is only needed for n=1, in this case SF2 must be set to the nucleon structure function F2." (* Messages: *) GenNConv2::toobad1 = "The nucleon structure function must be specified as the \ last argument." GenNConv2::toobad2 = "The structure functions xF1 together with F2 \ must be specified in the generalized convolution for xF1." GenNConv2::toobad3 = "Too many structure functions are given." GenNConv3::toobad1 = "The nucleon structure function must be specified." GenNConv3::toobad2 = "The nucleon F2 must also be specified in the convolution for xF1." Begin["`Private`"] (* ALL ENERGIES AND MOMENTA ARE IN FM UNITS *) (* hc = Mev fm conversion constant *) (* hc2 = Gev^2 Fm^2 conversion constant *) (* mN = Isoscalar nucleon mass in fm^-1 *) (* e0 & m0 are parameters of spectrum of residual nucleus *) (* pcut = cut in momentum integration *) hc = 197.327053; hc2=(hc/1000.)^2; mN = (938.27231 + 939.56563)/2. /hc; mpi= 139.57 /hc; (* delm2=(mN+mpi)^2 - mN^2;*) delm2=0; (* pcut = 1.7*mN;*) pcut = 1000. /hc; (* cut at 1 GeV the momentum integrals *) ecut = 1000. /hc; (* energy integral cut *) caux=10.^6; (* Options of numerical integration: *) (* intopts = {PrecisionGoal -> 4, Method -> MultiDimensional, MaxRecursion -> 30};*) intopts4 = {PrecisionGoal -> 4, MaxRecursion -> 33}; intopts3 = {PrecisionGoal -> 3, MaxRecursion -> 33}; (** GenNConv2 function **) GenNConv2[nsf_?IntegerQ, MD_, e0_?NumericQ, m0_?NumericQ, x_?NumericQ, qsq_?NumericQ, opts___] := ( (* opts are input nucleon structure functions *) If[Length[{opts}]==0, (Message[GenNConv2::toobad1]; Abort[ ])]; If[nsf==1 && Length[{opts}] < 2, (Message[GenNConv2::toobad2]; Abort[ ])]; If[Length[{opts}] > 2, (Message[GenNConv2::toobad3]; Abort[ ])]; (* Auxilary parameters: *) q2 = qsq/hc2; q0=q2/(2*mN*x); q=Sqrt[q2+q0^2]; gam = q/q0; gam2 = 1 + 2*mN^2*x/q2; mstar = m0*q/(m0+q0+mN); pstar = mN*(1 - x*(1+delm2/q2) + gam2*e0/mN)/gam; y0 = 2*pstar/mstar; cstar = Sqrt[Abs[y0]]; (* Below we calculate the integration limits for given x and q. Re was introduced here to assure that integration variable is real. The problem appears when one evaluates (Sqrt[Abs[y]])^2+y. It must equal 0 if y<0, but it is not due to finite precision of numerical calculations! As the result pminus and pplus become complex at the edge of the integration region. This causes trouble in the If operator etc. A poor man solution is to apply Re to the integration limits in y explicitly. *) pmax[y_] := (pplus = mstar*(y + Sqrt[y^2 + y0]) //Re; If[pplus < pcut, pplus, pcut]); pmin[y_] := (pminus= mstar*(y - Sqrt[y^2 + y0]) //Re; If[pminus < pcut, pminus, pcut]); (* Integrand in convolution formula. Note that q2 and psq are in Fm units. *) Clear[feff,x1,pz,t,psq]; feff := Which[nsf==2, (1+pz*gam/mN)*(1 + (4*psq+6*t)*x1^2/q2 )/gam^2 * {opts}[[1]][x1,q2,psq], nsf==1, (1+pz*gam/mN)*({opts}[[1]][x1,q2,psq] + x1^2*t/q2* {opts}[[2]][x1,q2,psq]), nsf==3, (1+pz/mN/gam)* {opts}[[1]][x1,q2,psq] ] * caux ; (* feff=1.; *) (* Now we integrate over y and p. We change the variable in momentum integral z=exp(-p) in order to avoid problems with large cutoffs. *) Which[ y0 > 0, NIntegrate[(p=Log[1/z];e = e0 - p^2/(2*m0); pz=p*y; t=p^2*(1-y^2); psq=(mN+e)^2 - p^2; x1=x/(1+(e+gam*pz)/mN); p^2*MD[p,y]* feff /z ), {y, -1,1},{z,Exp[-pmax[y]],1}, intopts3//Evaluate ]/2 , (y0 <= 0 && y0 > -1), NIntegrate[(p=Log[1/z];e = e0 - p^2/(2*m0); pz=p*y; t=p^2*(1-y^2); psq=(mN+e)^2 - p^2; x1=x/(1+(e+gam*pz)/mN); p^2*MD[p,y]* feff /z ), {y,cstar,1},{z,Exp[-pmax[y]],Exp[-pmin[y]]}, intopts3//Evaluate ]/2 , y0 <= -1, 0] /caux )(** End of GenNConv2 **) (*** NConv2 function Convolution integral for PDFs in kinematics corresponding to Bjorken limit. ***) NConv2[xpdf_, MD_, e0_?NumericQ, m0_?NumericQ, x_?NumericQ, qsq_?NumericQ] := ( q2 = qsq/hc2; q = q2/(2*mN*x); mstar = m0*q/(m0+q+mN); pstar = mN*(1 - x + e0/mN); y0 = 2*pstar/mstar; cstar = Sqrt[ Abs[y0] ]; (* Below we calculate the integration limits for given x and q. Re was introduced here to assure that integration variable is real. The problem appears when one evaluates (Sqrt[Abs[y]])^2+y. It must equal 0 if y<0, but it is not due to finite precision of numerical calculations! As the result pminus and pplus become complex at the edge of the integration region. This causes trouble in the If operator etc. A poor man solution is to apply Re to the integration limits in y explicitly. *) pmax[y_] := (pplus = mstar*(y + Sqrt[y^2 + y0]) //Re; If[pplus < pcut, pplus, pcut]); pmin[y_] := (pminus= mstar*(y - Sqrt[y^2 + y0]) //Re; If[pminus < pcut, pminus, pcut]); Clear[x1,pz,t,psq]; (* Now we integrate over y and p. We change the variable in momentum integral z=exp(-p) in order to avoid problems with large cutoffs. *) Which[ pstar > 0, NIntegrate[(p=Log[1/z];e = e0 - p^2/(2*m0); pz=p*y; t=p^2*(1-y^2); psq=(mN+e)^2 - p^2; x1=x/(1+(e+pz)/mN); p^2 * MD[p,y] * (1+pz/mN) * xpdf[x1,q2,psq] /z *caux ), {y, -1,1},{z,Exp[-pmax[y]],1}, intopts3//Evaluate ]/2 , (pstar <= 0 && pstar > -mstar/2), NIntegrate[(p=Log[1/z];e = e0 - p^2/(2*m0); pz=p*y; t=p^2*(1-y^2); psq=(mN+e)^2 - p^2; x1=x/(1+(e+pz)/mN); p^2 * MD[p,y] * (1+pz/mN) * xpdf[x1,q2,psq] /z *caux ), {y,cstar,1},{z,Exp[-pmax[y]],Exp[-pmin[y]]}, intopts3//Evaluate ]/2 , pstar <= -mstar/2, 0] /caux )(*** End of NConv2 ***) (*** GenNConv3 function ***) GenNConv3[nsf_?IntegerQ, Sdummy_, e0_?NumericQ, m0_?NumericQ, x_?NumericQ, q_?NumericQ, opts___] :=( If[Length[{opts}]==0, (Message[GenNConv3::toobad1]; Abort[])]; If[nsf==1 && Length[{opts}] < 2, (Message[GenNConv3::toobad2]; Abort[])]; (* Auxilary parameters: *) q2 = q/hc2; gam = Sqrt[1+ 4*mN^2*x^2/q2]; gam2 = 1 + 2*mN^2*x/q2; mstar = gam/(gam2/m0 + 2*mN*x/q2); mprim= mstar*gam2/gam; pstar = 2*mN*(1-x + e0*gam2/mN)/gam; ET=mstar*(mstar + pstar)/(2*mprim); Estar=pstar*mstar/(2*mprim); eprimmin1 = If[Estar > ecut, -ecut, -Estar]; eprimmin2 = If[ET > ecut, -ecut, -ET]; (* Below we calculate the integration limits for given x and q. Re was introduced here to assure that integration variable is real. The problem appears when one evaluates (Sqrt[Abs[y]])^2+y. It must equal 0 if y<0, but it is not (!) due to finite precision of calculations! As the result pminus and pplus become complex at the edge of the integration region. This causes trouble in the If operator etc. A poor man solution is to apply Re to the integration limits in y explicitly. *) Clear[eprim,y,p,feff,x1,pz,t,psq]; pstar1[eprim_] := pstar + 2*eprim*gam2/gam; pmx[y_,eprim_] := (pplus = mstar*(y+Sqrt[y^2 + pstar1[eprim]/mstar]) //Re; If[pplus < pcut, pplus, pcut]); pmn[y_,eprim_] := (pminus = mstar*(y-Sqrt[y^2 + pstar1[eprim]/mstar]) //Re; If[pminus < pcut, pminus, pcut]); (* The integrand *) feff := Which[nsf==2, (1+pz*gam/mN)*(1 + (4*psq+6*t)*x1^2/q2)/gam^2 * {opts}[[1]][x1,q2,psq], nsf==1, (1+pz*gam/mN)*({opts}[[1]][x1,q2,psq] + x1^2*t/q2* {opts}[[2]][x1,q2,psq]), nsf==3, (1+pz/mN/gam)* {opts}[[1]][x1,q2,psq] ]; (* feff=1; *) (* Now we integrate over y and p: *) Which[ pstar > 0, ( NIntegrate[( e = eprim + e0 - p^2/(2*m0); pz=p*y; t=p^2*(1-y^2); psq=mN^2 + 2*mN*e - p^2; x1=x/(1+(e+gam*pz)/mN); p^2*Sdummy[eprim,p,y]*feff), {eprim,eprimmin1, eprimmin1/2, 0}, {y,-1,0, 1}, {p, 0, pmx[y,eprim]/2, pmx[y,eprim]}, intopts3//Evaluate ] + If[Estar>ecut, 0,( NIntegrate[( e = eprim + e0 - p^2/(2*m0);pz=p*y; t=p^2*(1-y^2); psq=mN^2 + 2*mN*e - p^2; x1=x/(1+(e+gam*pz)/mN); p^2*Sdummy[eprim,p,y]*feff), {eprim, eprimmin2, (eprimmin2-Estar)/2, -Estar}, {y, Sqrt[Abs[pstar1[eprim]/mstar]] //Re, 1}, {p, pmn[y,eprim], (pmn[y,eprim]+pmx[y,eprim])/2, pmx[y,eprim]}, intopts3//Evaluate ]) ])/(4 Pi), (pstar <= 0 && pstar >= -mstar), ( NIntegrate[(e = eprim + e0 - p^2/(2*m0);pz=p*y; t=p^2*(1-y^2); psq=mN^2 + 2*mN*e - p^2; x1=x/(1+(e+gam*pz)/mN); p^2*Sdummy[eprim,p,y]*feff), {eprim, eprimmin2, 0}, {y, Sqrt[Abs[pstar1[eprim]/mstar]] //Re, 1}, {p, pmn[y,eprim], (pmn[y,eprim]+pmx[y,eprim])/2, pmx[y,eprim]}, intopts3//Evaluate ])/(4 Pi) ] )(*** End of GenNConv3 ***) End[] Protect[Evaluate[prtcted]] EndPackage[] Print["If you see no error messages the package NuclearConvolution \ has been successfully loaded. \ For more information type ?GenNConv2 GenNConv3 NConv2."]; (* End *)