From d1a6f1ef7aa347d7b398caf5d31bd3077342c101 Mon Sep 17 00:00:00 2001 From: Valerio Bertone Date: Thu, 2 Feb 2017 16:59:43 +0100 Subject: [PATCH] Adjusting code for the computation of the derivative of PDFs --- src/Evolution/LHAPDFgridDerivative.f | 498 ++++++++++----------------- 1 file changed, 190 insertions(+), 308 deletions(-) diff --git a/src/Evolution/LHAPDFgridDerivative.f b/src/Evolution/LHAPDFgridDerivative.f index f5662501..1257dd05 100644 --- a/src/Evolution/LHAPDFgridDerivative.f +++ b/src/Evolution/LHAPDFgridDerivative.f @@ -21,7 +21,6 @@ subroutine LHAPDFgridDerivative(Nrep,fname) 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" @@ -29,34 +28,28 @@ subroutine LHAPDFgridDerivative(Nrep,fname) * Input Variables * integer Nrep - logical islhapdf6 character*50 fname ** * Intenal Variables * integer ln - integer i,ix,iq2,iq2c,ipdf,ilep,krep,iq2in,iq2fi + integer ix,iq2,iq2c,ipdf,krep,iq2in,iq2fi integer nfin,nffi integer isg,nQ(3:7) integer nfmax - double precision qref,mth(4:6) double precision lambda3,lambda4,lambda5,lambdaNF - double precision xbLHA(nxmax),q2LHA(nq2max) + 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 dxPDF,dxgamma double precision AlphaQCD - character*30 str + character*30 str + character*3 ids(-6:6) parameter(eps=1d-4) - parameter(eps2=1d-12) + 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 * @@ -64,10 +57,11 @@ subroutine LHAPDFgridDerivative(Nrep,fname) 1 1d0,50,1d0,1d10) call SetQLimits(dsqrt(q2minLHA)-offset,dsqrt(q2maxLHA)+offset) call SetNumberOfGrids(3) - call SetGridParameters(1,150,3,xminLHA) - call SetGridParameters(2,100,5,xmLHA) - call SetGridParameters(3,100,5,0.65d0) + call SetGridParameters(1,85,3,xminLHA) + call SetGridParameters(2,75,5,xmLHA) + call SetGridParameters(3,55,5,0.65d0) call LockGrids(.true.) + call SetPDFEvolution("exactmu") * * Initialize APFEL * @@ -112,6 +106,12 @@ subroutine LHAPDFgridDerivative(Nrep,fname) * 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 @@ -135,12 +135,6 @@ subroutine LHAPDFgridDerivative(Nrep,fname) nffi = 3 endif if(nffi.gt.nfmax) nffi = nfmax -* -* Number of points per subgrid -* - do isg=3,7 - nQ(isg) = 0 - enddo * isg = nfin do iq2=1,nq2LHA @@ -156,45 +150,54 @@ subroutine LHAPDFgridDerivative(Nrep,fname) endif enddo * -* If the LHAPDF6 format is to be used redefine the grid with the subgrids +* Make sure that all subgrids have at least two points. * - if(islhapdf6())then - lnQmin = dlog( ( q2minLHA ) / Lambda2 ) - if(nfin.eq.nffi)then - lnQmax = dlog( q2maxLHA / Lambda2 ) - else - lnQmax = dlog( ( m2th(nfin+1) ) / Lambda2 ) + 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 * - 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 +* 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 * - 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) + 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 - nfin = Nf_FF - nffi = Nf_FF +* do iq2=1,nq2LHA q2LHA(iq2) = Lambda2 * dexp( lnQmin 1 * dexp( dble( iq2 - 1 ) / dble( nq2LHA - 1 ) @@ -202,97 +205,140 @@ subroutine LHAPDFgridDerivative(Nrep,fname) enddo endif * -* LHAPDF5 output -* - if(.not.islhapdf6())then -* -* Open the output file -* - ln = index(fname,char(0)) - 1 - if(ln.eq.-1) ln = index(fname,char(32)) - 1 - open(unit=13,status="unknown",file=fname(1:ln)//".LHgrid") -* -* Write header -* - write(13,*) "'Version'"," '5.8.9'" - write(13,*) "'Description: set generated with APFEL'" - write(13,*) "'",fname(1:ln),".LHgrid'" - write(13,*) "'Reference: arXiv:1310.1394'" - write(13,*) "'APFEL: A PDF Evolution Library'" - write(13,*) "'V. Bertone, S. Carrazza, J. Rojo'" - write(13,*)"' ***********************'" - write(13,*)"' *** PDF DERIVATIVES ***'" - write(13,*)"' ***********************'" - write(13,*) "'This set has ",Nrep," member PDF'" - write(13,*) "'mem=0 --> average on replicas. '" - write(13,*) "'Evolution: pol. interpolation on the LH grid.'" -* - write(13,*) "'Alphas:'" - if(ipt.eq.0)then - write(13,*) "'Variable',","'lo',","'EvolCode'" - elseif(ipt.eq.1)then - write(13,*) "'Variable',","'nlo',","'EvolCode'" - elseif(ipt.eq.2)then - write(13,*) "'Variable',","'nnlo',","'EvolCode'" - endif +* Compute alphas on the grid being careful with the grids * - qref = dsqrt(q2_ref_qcd) - do i=4,6 - mth(i) = dsqrt(m2th(i)) + 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 * - write(13,*) 1,",",qref,",",mth(4),",",mth(5),",",mth(6) +* 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" * - write(13,*) "'MinMax:'" - write(13,*) Nrep,",",1 - write(13,*) xminLHA,",",xmaxLHA,",",q2minLHA,",",q2maxLHA +* LHAPDF6 output * - write(13,*) "'QCDparams:'" - write(13,*) Nrep,",",1 - write(13,*) lambda4,",",lambda5 + 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(13,*) "'Parameterlist:'" - write(13,*) "'list',",Nrep,",",1 +* Write header * - do i=0,Nrep - write(13,*) alpha_ref_qcd - enddo + 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 + write(13,*) "Flavors: [", + 1 (ids(ipdf),",",ipdf=-nfMaxPDFs,nfMaxPDFs), + 2 " 22]" + 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)) + 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) * - write(13,* ) "'Evolution:'" - if(ipt.eq.0)then - write(13,*) "'lo',",q2minLHA,",",kren - elseif(ipt.eq.1)then - write(13,*) "'nlo',",q2minLHA,",",kren - elseif(ipt.eq.2)then - write(13,*) "'nnlo',",q2minLHA,",",kren +* 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(Th.eq."QCD")then - write(13,*) "'NNPDF20int'" + if (krep.eq.0) then + write(13,"(a)") "PdfType: central" else - write(13,*) "'NNPDF20intqed'" + write(13,"(a)") "PdfType: replica" endif - write(13,*) Nrep,",",1 + write(13,"(a)") "Format: lhagrid1" + write(13,"(a)") "---" * -* Compute and write the grids in x and Q2 -* - write(13,*) nxLHA - do ix=1,nxLHA - write(13,*) xbLHA(ix) - enddo + write(6,*) "Evaluating replica",krep," ..." + call SetReplica(krep) * - write(13,*) nq2LHA - write(13,*) q2minLHA - do iq2=1,nq2LHA - write(13,*) q2LHA(iq2) - enddo + iq2in = 1 + iq2fi = nQ(nfin) + do isg=nfin,nffi + write(13,*) (xbLHA(ix),ix=1,nxLHA) + write(13,*) (dsqrt(q2LHA(iq2)),iq2=iq2in,iq2fi) * -* Compute derivarive and write it on file + if(Th.eq."QCD") then + write(13,*) (ids(ipdf),ipdf=-nfMaxPDFs,nfMaxPDFs) + else + write(13,*) (ids(ipdf),ipdf=-nfMaxPDFs,nfMaxPDFs), + 1 " 22" + endif * - write(13,*) Nrep - do krep=0,Nrep - write(6,*) "Evaluating replica",krep," ..." - call SetReplica(krep) - do iq2=1,nq2LHA + do iq2=iq2in,iq2fi call DeriveAPFEL(dsqrt(q2LHA(iq2))) do ix=1,nxLHA do ipdf=-6,6 @@ -304,197 +350,33 @@ subroutine LHAPDFgridDerivative(Nrep,fname) * if(Th.eq."QCD")then do ix=1,nxLHA - do iq2=1,nq2LHA - write(13,44) (xpdfLHA(ipdf,ix,iq2),ipdf=-6,6) + do iq2=iq2in,iq2fi + write(13,40) (xpdfLHA(ipdf,ix,iq2), + 1 ipdf=-nfMaxPDFs,nfMaxPDFs) enddo enddo else do ix=1,nxLHA - do iq2=1,nq2LHA - write(13,45) (xpdfLHA(ipdf,ix,iq2),ipdf=-6,6), - 1 xgammalha(ix,iq2) + do iq2=iq2in,iq2fi + write(13,50)(xpdfLHA(ipdf,ix,iq2), + 1 ipdf=-nfMaxPDFs,nfMaxPDFs), + 2 xgammaLHA(ix,iq2) enddo enddo endif + write(13,"(a)") "---" + iq2in = iq2in + nQ(isg) + iq2fi = iq2fi + nQ(isg+1) enddo -* - write(13,* ) "'End:'" - close(13) -* - write(6,*) "File ",fname(1:ln),".LHgrid produced!" - write(6,*) " " - else -* -* 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 - ",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: [-6, -5, -4, -3, -2, -1,", - 1 "21, 1, 2, 3, 4, 5, 6]" - else - if(LeptEvol)then - write(13,*) "Flavors: [-6, -5, -4, -3, -2, -1, 21,", - 1 " 1, 2, 3, 4, 5, 6, 22, ", - 2 "11, -11, 13, -13, 15, -15]" - else - write(13,*) "Flavors: [-6, -5, -4, -3, -2, -1, 21,", - 1 " 1, 2, 3, 4, 5, 6, 22]" - endif - endif - write(13,*) "OrderQCD:",ipt - write(13,*) "FlavorScheme: variable" - write(13,*) "NumFlavors: 6" - 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(m2th(4)) - write(13,*) "MBottom:",dsqrt(m2th(5)) - write(13,*) "MTop:",dsqrt(m2th(6)) - write(13,*) "AlphaS_MZ:",alpha_ref_qcd - write(13,*) "AlphaS_OrderQCD:",ipt - write(13,*) "AlphaS_Type: ipol" - write(13,*)"AlphaS_Qs: [",(dsqrt(q2LHA(iq2)),",",iq2=1,nq2LHA), - 1 "]" - write(13,*) "AlphaS_Vals: [",(AlphaQCD(dsqrt(q2LHA(iq2))), - 1 ",",iq2=1,nq2LHA), "]" - write(13,*) "AlphaS_Lambda4:", lambda4 - write(13,*) "AlphaS_Lambda5:", lambda5 close(13) + enddo * -* 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,*) "PdfType: central" - else - write(13,*) "PdfType: replica" - endif - write(13,*) "Format: lhagrid1" - write(13,*) "---" -* - write(6,*) "Evaluating replica",krep," ..." - call SetReplica(krep) -* - 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,"(a)") "-6 -5 -4 -3 -2 -1 21 1 2 3 4 5 6" - else - if(LeptEvol)then - write(13,"(a,a,a)") "-6 -5 -4 -3 -2 -1 21", - 1 " 1 2 3 4 5 6 22 ", - 2 "11 -11 13 -13 15 -15" - else - write(13,"(a,a)") "-6 -5 -4 -3 -2 -1 21", - 1 " 1 2 3 4 5 6 22" - endif - endif -* - do iq2=iq2in,iq2fi - if(iq2.eq.iq2in.and.isg.ne.nfin)then - call DeriveAPFEL(dsqrt(q2LHA(iq2)+eps2)) - elseif(iq2.eq.iq2fi.and.isg.ne.nffi)then - call DeriveAPFEL(dsqrt(q2LHA(iq2)-eps2)) - else - call DeriveAPFEL(dsqrt(q2LHA(iq2))) - endif - do ix=1,nxLHA - do ipdf=-6,6 - xpdfLHA(ipdf,ix,iq2) = dxPDF(ipdf,xbLHA(ix)) - enddo - do ilep=-3,3 - xlepLHA(ilep,ix,iq2) = 0d0 - enddo - xgammaLHA(ix,iq2) = dxgamma(xbLHA(ix)) - enddo - enddo -* - if(Th.eq."QCD")then - do ix=1,nxLHA - do iq2=iq2in,iq2fi - write(13,44) (xpdfLHA(ipdf,ix,iq2),ipdf=-6,6) - enddo - enddo - else - if(LeptEvol)then - do ix=1,nxLHA - do iq2=iq2in,iq2fi - write(13,46)(xpdfLHA(ipdf,ix,iq2),ipdf=-6,6), - 1 xgammaLHA(ix,iq2), - 2 xlepLHA(1,ix,iq2), - 3 xlepLHA(-1,ix,iq2), - 4 xlepLHA(2,ix,iq2), - 5 xlepLHA(-2,ix,iq2), - 6 xlepLHA(3,ix,iq2), - 7 xlepLHA(-3,ix,iq2) - enddo - enddo - else - do ix=1,nxLHA - do iq2=iq2in,iq2fi - write(13,45)(xpdfLHA(ipdf,ix,iq2),ipdf=-6,6), - 1 xgammaLHA(ix,iq2) - enddo - enddo - endif - endif - write(13,*) "---" - iq2in = iq2in + nQ(isg) - iq2fi = iq2fi + nQ(isg+1) - enddo - close(13) - enddo -* - write(6,*) "File ",fname(1:ln)," grid produced!" - write(6,*) " " - endif + write(6,*) achar(27)//"[1;32m" + write(6,*) "File ",fname(1:ln)," grid produced" + write(6,*) achar(27)//"[0m" * - 44 format(13(1x,es14.7)) - 45 format(14(1x,es14.7)) - 46 format(20(1x,es14.7)) + 40 format(13(es14.7,1x)) + 50 format(14(es14.7,1x)) * return end