************************************************************************ * * LHAPDFgrid.f: * * If LHAPDF5 is found, this subrotine produces the *.LHgrid file * that can be used with LHAPDF v5.9. If instead LHAPDF6 is found * a folder contained the PDF set in a suitable format is produced. * Note that if "Qin" < 0 then the internal evolution of the input * PDF set will be used. * ************************************************************************ subroutine LHAPDFgrid(Nrep,Qin,fname) * implicit none * include "../commons/grid.h" include "../commons/ipt.h" include "../commons/alpha_ref_QCD.h" include "../commons/m2th.h" include "../commons/kren.h" include "../commons/Th.h" include "../commons/LHAgrid.h" include "../commons/Evs.h" include "../commons/MaxFlavourPDFs.h" include "../commons/MaxFlavourAlpha.h" include "../commons/LeptEvol.h" include "../commons/Nf_FF.h" include "../commons/pdfset.h" include "../commons/f0ph.h" ** * Input Variables * integer Nrep double precision Qin character*50 fname ** * Intenal Variables * integer ln integer ix,iq2,iq2c,ipdf,ilep,krep,iq2in,iq2fi integer nfin,nffi integer isg,nQ(3:7) integer nfmax integer alpha double precision lambda3,lambda4,lambda5,lambdaNF double precision xbLHA(nxmax),q2LHA(nq2max),as(nxmax) double precision lnQmin,lnQmax double precision eps,eps2,offset double precision xpdfLHA(-6:6,nxmax,nq2max) double precision xgammaLHA(nxmax,nq2max) double precision xlepLHA(-6:6,nxmax,nq2max) double precision xPDFj,xgammaj,xLeptonj double precision AlphaQCD character*30 str character*100 pdfsetbkp character*3 ids(-6:6) parameter(eps=1d-4) parameter(eps2=1d-8) parameter(offset=1d-2) * double precision fqpre1(ngrid_max,-6:6,0:nint_max) double precision flpre1(ngrid_max,-3:3,0:nint_max) common / pretab1APFEL / fqpre1,flpre1 * * Specify initialization for the LHgrid * if(InLHgrid.ne."done") call SetLHgridParameters(100,50,1d-9,1d-1, 1 1d0,50,1d0,1d10) call SetQLimits(dsqrt(q2minLHA)-offset,dsqrt(q2maxLHA)+offset) * * Adjust internal grid only if no external grid has been provided. * if(.not.ThereAreExtGrids)then call SetNumberOfGrids(3) call SetGridParameters(1,85,3,xminLHA) call SetGridParameters(2,75,5,xmLHA) call SetGridParameters(3,55,5,0.65d0) call LockGrids(.true.) endif * * Initialize APFEL * call initializeAPFEL * * Report grid parameters * write(6,*) "Report of the LHAPDF grid parameters:" write(6,"(a,i3,a,f11.9,a,f6.4,a)") " - x-space grid: ",nxLHA, 1 " points in [",xminLHA," :",xmaxLHA,"]" write(6,"(a,f6.4)") " transition from log to lin in x = ",xmLHA write(6,"(a,i3,a,f5.2,a,f13.1,a)") " - Q2-space grid: ",nq2LHA, 1 " points in [",q2minLHA,":",q2maxLHA,"] GeV^2" write(6,*) " " * * Compute the values of lambdaQCD * lambda5 = lambdaNF(5,alpha_ref_qcd,q2_ref_qcd) lambda4 = lambdaNF(4,alpha_ref_qcd,q2_ref_qcd) lambda3 = lambdaNF(3,alpha_ref_qcd,q2_ref_qcd) * * Define maximun number of flavours * nfmax = max(nfMaxPDFs,nfMaxAlpha) * * Compute x-space grid * do ix=1,nxLHA if(ix.le.nxmLHA)then xbLHA(ix) = xminLHA * ( xmLHA / xminLHA ) 1 **( 2d0 * dble( ix-1 ) / dble( nxLHA - 1 ) ) else xbLHA(ix) = xmLHA + ( xmaxLHA - xmLHA ) 1 * ( dble( ix - nxmLHA - 1 ) 2 / dble( nxLHA - nxmLHA - 1 ) ) endif enddo * * Compute Q2 grid. * Use a distribution of the Q2 nodes uniform in ln(ln(Q2/Lambda2)) * and that has nodes on the heavy quark thresholds. * lnQmin = dlog( q2minLHA / Lambda2 ) lnQmax = dlog( q2maxLHA / Lambda2 ) * * Initialize number of points per subgrid * do isg=3,7 nQ(isg) = 0 enddo * if(Evs.eq."VF")then if(q2minLHA.gt.m2th(6))then nfin = 6 elseif(q2minLHA.gt.m2th(5))then nfin = 5 elseif(q2minLHA.gt.m2th(4))then nfin = 4 else nfin = 3 endif if(nfin.gt.nfmax) nfin = nfmax * if(q2maxLHA.gt.m2th(6))then nffi = 6 elseif(q2maxLHA.gt.m2th(5))then nffi = 5 elseif(q2maxLHA.gt.m2th(4))then nffi = 4 else nffi = 3 endif if(nffi.gt.nfmax) nffi = nfmax * isg = nfin do iq2=1,nq2LHA q2LHA(iq2) = Lambda2 * dexp( lnQmin 1 * dexp( dble( iq2 - 1 ) / dble( nq2LHA - 1 ) 2 * dlog( lnQmax / lnQmin ) ) ) if(q2LHA(iq2).lt.m2th(isg+1)+eps)then nQ(isg) = nQ(isg) + 1 else isg = isg + 1 if(isg.gt.nfmax) isg = nfmax nQ(isg) = nQ(isg) + 1 endif enddo * * Make sure that all subgrids have at least two points. * do isg=nfin,nffi-1 if(nQ(isg).lt.2)then nQ(isg) = nQ(isg) + 1 nQ(isg+1) = nQ(isg+1) - 1 endif enddo * * Redefine the grid with the subgrids * lnQmin = dlog( ( q2minLHA ) / Lambda2 ) if(nfin.eq.nffi)then lnQmax = dlog( q2maxLHA / Lambda2 ) else lnQmax = dlog( ( m2th(nfin+1) ) / Lambda2 ) endif * iq2c = 0 do isg=nfin,nffi do iq2=1,nQ(isg) iq2c = iq2c + 1 q2LHA(iq2c) = Lambda2 * dexp( lnQmin 1 * dexp( dble( iq2 - 1 ) 2 / dble( nQ(isg) - 1 ) 3 * dlog( lnQmax / lnQmin ) ) ) enddo lnQmin = dlog( ( q2LHA(iq2c) ) / Lambda2 ) if(isg.eq.nffi-1)then lnQmax = dlog( q2maxLHA / Lambda2 ) else lnQmax = dlog( ( m2th(isg+2) ) / Lambda2 ) endif enddo * if(iq2c.ne.nq2LHA)then write(6,*) "In LHAPDFgrid.f:" write(6,*) "Mismatch in the Number of Q2 nodes" write(6,*) "- Expected = ",nq2LHA write(6,*) "- Found = ",iq2c call exit(-10) endif else nfin = Nf_FF nffi = Nf_FF * nQ(Nf_FF) = nq2LHA * do iq2=1,nq2LHA q2LHA(iq2) = Lambda2 * dexp( lnQmin 1 * dexp( dble( iq2 - 1 ) / dble( nq2LHA - 1 ) 2 * dlog( lnQmax / lnQmin ) ) ) enddo endif * * Compute alphas on the grid being careful with the grids * iq2in = 1 iq2fi = nQ(nfin) do isg=nfin,nffi do iq2=iq2in,iq2fi if(iq2.eq.iq2in.and.isg.ne.nfin)then as(iq2) = AlphaQCD(dsqrt(q2LHA(iq2)+eps2)) elseif(iq2.eq.iq2fi.and.isg.ne.nffi)then as(iq2) = AlphaQCD(dsqrt(q2LHA(iq2)-eps2)) else as(iq2) = AlphaQCD(dsqrt(q2LHA(iq2))) endif enddo iq2in = iq2in + nQ(isg) iq2fi = iq2fi + nQ(isg+1) enddo * * Define quark IDs * ids(-6) = " -6" ids(-5) = " -5" ids(-4) = " -4" ids(-3) = " -3" ids(-2) = " -2" ids(-1) = " -1" ids(0) = " 21" ids(1) = " 1" ids(2) = " 2" ids(3) = " 3" ids(4) = " 4" ids(5) = " 5" ids(6) = " 6" * * LHAPDF6 output * ln = index(fname,char(0)) - 1 if(ln.eq.-1) ln = index(fname,char(32)) - 1 * creating main folder call mkdir(fname(1:ln)) * creating info file open(unit=13,status="unknown",file=fname(1:ln)//"/" 1 //fname(1:ln)//".info") * * Write header * write(13,*) 'SetDesc: "set generated with APFEL - ', 1 fname(1:ln),'"' write(13,*) "Authors: V. Bertone, S. Carrazza, J. Rojo" write(13,*) "Reference: ArXiv:1310.1394" write(13,*) "Format: lhagrid1" write(13,*) "DataVersion: 1" write(13,*) "NumMembers:",Nrep+1 write(13,*) "Particle: 2212" if(Th.eq."QCD") then write(13,*) "Flavors: [", 1 (ids(ipdf),",",ipdf=-nfMaxPDFs,nfMaxPDFs-1), 2 ids(nfMaxPDFs),"]" else if(LeptEvol)then write(13,*) "Flavors: [", 1 (ids(ipdf),",",ipdf=-nfMaxPDFs,nfMaxPDFs), 2 " 22, 11, -11, 13, -13, 15, -15]" else write(13,*) "Flavors: [", 1 (ids(ipdf),",",ipdf=-nfMaxPDFs,nfMaxPDFs), 2 " 22]" endif endif write(13,*) "OrderQCD:",ipt write(13,*) "FlavorScheme: variable" write(13,*) "NumFlavors: ",nfMaxPDFs write(13,*) "ErrorType: replicas" write(13,*) "XMin:",xminLHA write(13,*) "XMax:",xmaxLHA write(13,*) "QMin:",dsqrt(q2minLHA) write(13,*) "QMax:",dsqrt(q2maxLHA) write(13,*) "MZ:",dsqrt(q2_ref_qcd) write(13,*) "MUp: 0" write(13,*) "MDown: 0" write(13,*) "MStrange: 0" write(13,*) "MCharm:",dsqrt(m2ph(4)) write(13,*) "MBottom:",dsqrt(m2ph(5)) write(13,*) "MTop:",dsqrt(m2ph(6)) if(nfmax.ge.4) then write(13,*) "ThresholdCharm:",dsqrt(m2th(4)) else write(13,*) "ThresholdCharm:",dsqrt(q2maxLHA)+1d0 endif if(nfmax.ge.5) then write(13,*) "ThresholdBottom:",dsqrt(m2th(5)) else write(13,*) "ThresholdBottom:",dsqrt(q2maxLHA)+2d0 endif if(nfmax.ge.6) then write(13,*) "ThresholdTop:",dsqrt(m2th(6)) else write(13,*) "ThresholdTop:",dsqrt(q2maxLHA)+3d0 endif write(13,*) "AlphaS_MZ:",alpha_ref_qcd write(13,*) "AlphaS_OrderQCD:",ipt write(13,*) "AlphaS_Type: ipol" write(13,*) "AlphaS_Qs: [",(dsqrt(q2LHA(iq2)),",", 1 iq2=1,nq2LHA-1),dsqrt(q2LHA(nq2LHA)),"]" write(13,*) "AlphaS_Vals: [",(as(iq2),",", 1 iq2=1,nq2LHA-1),as(nq2LHA),"]" write(13,*) "AlphaS_Lambda4:", lambda4 write(13,*) "AlphaS_Lambda5:", lambda5 close(13) * * Now Loop over all replicas and print to file * do krep=0,Nrep if (krep.lt.10) then write(str,'(i1)' ) krep open(unit=13,status="unknown",file=fname(1:ln)//"/" 1 //fname(1:ln)//"_000"//str(1:1)//".dat") elseif (krep.lt.100) then write(str,'(i2)' ) krep open(unit=13,status="unknown",file=fname(1:ln)//"/" 1 //fname(1:ln)//"_00"//str(1:2)//".dat") elseif (krep.lt.1000) then write(str,'(i3)' ) krep open(unit=13,status="unknown",file=fname(1:ln)//"/" 1 //fname(1:ln)//"_0"//str(1:3)//".dat") else write(str,'(i4)' ) krep open(unit=13,status="unknown",file=fname(1:ln)//"/" 1 //fname(1:ln)//"_"//str(1:4)//".dat") endif if (krep.eq.0) then write(13,"(a)") "PdfType: central" else write(13,"(a)") "PdfType: replica" endif write(13,"(a)") "Format: lhagrid1" write(13,"(a)") "---" * write(6,*) "Evaluating replica",krep," ..." call SetReplica(krep) * * Tabulate PDFs at the initial scale on the grid * if(Qin.gt.0d0)then do igrid=1,ngrid call initPDFs(Qin**2) do alpha=0,nin(igrid) do ipdf=-6,6 fqpre1(igrid,ipdf,alpha) = f0ph(ipdf,alpha) enddo do ipdf=-3,3 flpre1(igrid,ipdf,alpha) = f0lep(ipdf,alpha) enddo enddo enddo * Back up PDF name and set pretabulated PDFs as input pdfsetbkp = pdfset call SetPDFSet("pretabulated1") endif * iq2in = 1 iq2fi = nQ(nfin) do isg=nfin,nffi write(13,*) (xbLHA(ix),ix=1,nxLHA) write(13,*) (dsqrt(q2LHA(iq2)),iq2=iq2in,iq2fi) * if(Th.eq."QCD") then write(13,*) (ids(ipdf),ipdf=-nfMaxPDFs,nfMaxPDFs) else if(LeptEvol)then write(13,*) (ids(ipdf),ipdf=-nfMaxPDFs,nfMaxPDFs), 2 " 22 11 -11 13 -13 15 -15" else write(13,*) (ids(ipdf),ipdf=-nfMaxPDFs,nfMaxPDFs), 1 " 22" endif endif * do iq2=iq2in,iq2fi if(Qin.gt.0d0)then if(iq2.eq.iq2in.and.isg.ne.nfin)then call EvolveAPFEL(Qin,dsqrt(q2LHA(iq2)+eps2)) elseif(iq2.eq.iq2fi.and.isg.ne.nffi)then call EvolveAPFEL(Qin,dsqrt(q2LHA(iq2)-eps2)) else call EvolveAPFEL(Qin,dsqrt(q2LHA(iq2))) endif else if(iq2.eq.iq2in.and.isg.ne.nfin)then call EvolveAPFEL(dsqrt(q2LHA(iq2)+eps2), 1 dsqrt(q2LHA(iq2)+eps2)) elseif(iq2.eq.iq2fi.and.isg.ne.nffi)then call EvolveAPFEL(dsqrt(q2LHA(iq2)-eps2), 1 dsqrt(q2LHA(iq2)-eps2)) else call EvolveAPFEL(dsqrt(q2LHA(iq2)), 1 dsqrt(q2LHA(iq2))) endif endif do ix=1,nxLHA do ipdf=-6,6 xpdfLHA(ipdf,ix,iq2) = xPDFj(ipdf,xbLHA(ix)) enddo do ilep=-3,3 xlepLHA(ilep,ix,iq2) = xLeptonj(ilep,xbLHA(ix)) enddo xgammaLHA(ix,iq2) = xgammaj(xbLHA(ix)) enddo enddo * if(Th.eq."QCD")then do ix=1,nxLHA do iq2=iq2in,iq2fi write(13,40) (xpdfLHA(ipdf,ix,iq2), 1 ipdf=-nfMaxPDFs,nfMaxPDFs) enddo enddo else if(LeptEvol)then do ix=1,nxLHA do iq2=iq2in,iq2fi write(13,60)(xpdfLHA(ipdf,ix,iq2), 1 ipdf=-nfMaxPDFs,nfMaxPDFs), 2 xgammaLHA(ix,iq2), 3 xlepLHA(1,ix,iq2), 4 xlepLHA(-1,ix,iq2), 5 xlepLHA(2,ix,iq2), 6 xlepLHA(-2,ix,iq2), 7 xlepLHA(3,ix,iq2), 8 xlepLHA(-3,ix,iq2) enddo enddo else do ix=1,nxLHA do iq2=iq2in,iq2fi write(13,50)(xpdfLHA(ipdf,ix,iq2), 1 ipdf=-nfMaxPDFs,nfMaxPDFs), 2 xgammaLHA(ix,iq2) enddo enddo endif endif write(13,"(a)") "---" iq2in = iq2in + nQ(isg) iq2fi = iq2fi + nQ(isg+1) enddo close(13) * Restore PDF name if(Qin.gt.0d0) call SetPDFSet(pdfsetbkp) enddo * write(6,*) achar(27)//"[1;32m" write(6,*) "File ",fname(1:ln)," grid produced" write(6,*) achar(27)//"[0m" * 40 format(13(es14.7,1x)) 50 format(14(es14.7,1x)) 60 format(20(es14.7,1x)) * return end * ************************************************************************ subroutine pretabulatedPDFs1(jgrid,alpha,fq,fl) * implicit none * include "../commons/grid.h" ** * Input Variables * integer jgrid,alpha ** * Internal Variables * integer ifl,ilept double precision fqpre1(ngrid_max,-6:6,0:nint_max) double precision flpre1(ngrid_max,-3:3,0:nint_max) common / pretab1APFEL / fqpre1,flpre1 ** * Output Variables * double precision fq(-6:6) double precision fl(-3:3) * do ifl=-6,6 fq(ifl) = fqpre1(jgrid,ifl,alpha) enddo do ilept=-3,3 fl(ilept) = flpre1(jgrid,ilept,alpha) enddo * return end