LHAPDF/lhapdf-5.9.1/src/wrapNNPDF20qedgrid.f

   1 !***********************************************
   2 !
   3 !
   4 !     wrapNNPDF20qedgrid.f
   5 !     Routine called by LHAPDF package for calculating
   6 !     the value of all xpdfs at x and Q from replica KREP
   7 !     in a (x,Q) point as called by NNPDF.LHgrid file.
   8 !
   9 !     In 'wrapevolve.f' the package calls:
  10 !   IF(NAME(NSET).EQ.'NNPDF20intqed') call NNPDFINT20qedevolve(x,Q,f)
  11 !   IF(NAME(NSET).EQ.'NNPDF20intqed') call NNPDFINT20qedread(nset)
  12 !   IF(NAME(NSET).EQ.'NNPDF20intqed') call NNPDFINT20qedalfa(alfas,Q)
  13 !   IF(NAME(NSET).EQ.'NNPDF20intqed') call NNPDFINT20qedinit(nset,Eorder,Q2fit)
  14 !
  15 !***********************************************
  16
  17       subroutine NNPDFINT20qedevolve(XIN,QIN,XPDF,XPHOTON)
  18       IMPLICIT none
  19 !
  20       include 'parmsetup.inc'
  21       character*16 name(nmxset)
  22       integer nmem(nmxset),ndef(nmxset),mmem
  23       common/NAME/name,nmem,ndef,mmem
  24       integer nset,iset,pdfmem
  25       real*8 parm(nopmax)
  26 !
  27       double precision gridx(nmxgridx),gridq(nmxgridq)
  28       integer ngridx,ngridq,jx,jq
  29 !
  30       INTEGER order
  31       REAL*8 alfas,alphaNNPDF
  32       REAL*8 Eorder,Q2fit,mass
  33 !
  34       INTEGER MXREP
  35       PARAMETER(MXREP=1e2)
  36       INTEGER NREP
  37       common/nnpdf20CNREP/NREP
  38 !
  39       INTEGER NX,NXMAX,NQ2,NPL
  40       PARAMETER(NX=100,NXMAX=100,NQ2=50)
  41       PARAMETER(NPL=5000)
  42       INTEGER NXX,NQQ2
  43       REAL*8 Q2MIN,Q2MAX,XPDFMIN,XPDFMAX
  44       REAL*8 XG(NXMAX,nmxset),Q2G(NQ2,nmxset),XPDFEV(NXMAX,NQ2,-6:7,0:MXREP)
  45       INTEGER IX,IQ2
  46 !     This common different from NNPDF1.X
  47 !      common/nnpdf20CPDFGR/XPDFEV,XG,Q2G,IX,IQ2
  48       common/nnpdfallCPDFGR/XPDFEV,XG,Q2G,IX,IQ2
  49 !
  50       INTEGER ipt,imodev,ivfn,itmc
  51 !     This common is the same as for NNPDF1.0
  52       COMMON/nnpdf10EVFLAGS/ipt,imodev,ivfn,itmc
  53       REAL*8 q0,alfas0
  54       REAL*8 q20,qth(4:6)
  55 !     This common is the same as for NNPDF1.0
  56       COMMON/nnpdf10EVSCALE/q20,q2
  57       REAL*8 q2th(4:6),asref,q2ref
  58 !     This common is the same as for NNPDF1.0
  59       COMMON/nnpdf10vfns/q2th,asref,q2ref
  60 !
  61       INTEGER I,J,K
  62       INTEGER IPDF,KREP
  63       INTEGER IDUM,JDUM
  64       REAL*8 X,XIN,Q,QIN,QAS,QQ,Q2,QQ2,XPDF(-6:6),XPHOTON
  65       REAL*8 XCH
  66       PARAMETER(XCH=1D-1)
  67 !
  68       integer m,n,nmax,mmax,minq,maxq,midq,maxx,minx,midx
  69 !     order of pol. interpolation
  70       parameter(m=4,n=2)
  71       parameter(nmax=1e3,mmax=1e3)
  72       double precision dy,x1,x2,y,x1a(mmax),x2a(nmax),ya(mmax,nmax)
  73       integer ix1a(m),ix2a(n)
  74       integer offset
  75
  76       double precision XPDFMIN_INTER
  77       parameter(XPDFMIN_INTER=1d-7)
  78 !
  79       KREP = mmem
  80       X=XIN
  81       Q=QIN
  82
  83 !     Set correct scale
  84       Q2=Q**2d0
  85
  86       call getnset(iset)
  87       offset = 101*(iset-1)
  88
  89 !     Check kinematic point is within allowed range
  90
  91       call GetXminM(iset,KREP,XPDFMIN)
  92       call GetXmaxM(iset,KREP,XPDFMAX)
  93       call GetQ2maxM(iset,KREP,Q2MAX)
  94       call GetQ2minM(iset,KREP,Q2MIN)
  95 !
  96       IF ( X.LT.XPDFMIN_INTER .OR. X.GT.XPDFMAX ) THEN
  97          WRITE(6,2000)
  98  2000    FORMAT (2X,'PARTON INTERPOLATION: X OUT OF RANGE -- STOP')
  99          IF ( X.LT.XPDFMIN_INTER ) X = XPDFMIN_INTER
 100          IF ( X.GT.XPDFMAX ) X = XPDFMAX
 101          write(6,*) "X= ",X," XMAX, XMIN = ",XPDFMAX,XPDFMIN
 102       ENDIF
 103 !
 104       IF ( Q2.LT.Q2MIN .OR. Q2.GT.Q2MAX ) THEN
 105          WRITE(6,2001)
 106  2001    FORMAT (2X,'PARTON INTERPOLATION: Q2 OUT OF RANGE -- STOP')
 107          write(6,*) "Q2 ,Q2MIN, Q2MAX = ",Q2,Q2MIN,Q2MAX
 108          if( Q2.LT.Q2MIN ) Q2 = Q2MIN
 109          if( Q2.GT.Q2MAX ) Q2 = Q2MAX
 110       ENDIF
 111
 112 !     FIND NEAREST POINTS IN THE GRID
 113       MINX = 1
 114       MAXX = NX+1
 115  10   CONTINUE
 116       MIDX = (MINX+MAXX)/2
 117       IF(X.LT.XG(MIDX,iset)) THEN
 118          MAXX=MIDX
 119       ELSE
 120          MINX=MIDX
 121       END IF
 122       IF((MAXX-MINX).GT.1) GO TO 10
 123       IX = MINX
 124
 125       MINQ = 1
 126       MAXQ = NQ2+1
 127  20   CONTINUE
 128       MIDQ = (MINQ+MAXQ)/2
 129       IF(Q2.LT.Q2G(MIDQ,iset)) THEN
 130          MAXQ=MIDQ
 131       ELSE
 132          MINQ=MIDQ
 133       END IF
 134       IF((MAXQ-MINQ).GT.1) GO TO 20
 135       IQ2 = MINQ
 136
 137 !
 138 !     POLYNOMIAL INTERPOLATION
 139 !
 140 !     uncomment for 3rd order interp.
 141 !     IF((IX+(M-1)/2).GT.NX) IX = NX - (M-1)/2
 142 !     IF((IX-(M-1)/2).LT.1) IX = (M+1)/2
 143 !     IDUM = 0
 144 !     DO I = -(M-1)/2,(M-1)/2,1
 145 !     IDUM = IDUM +1
 146 !     IX1A(IDUM) = IX + I
 147 !     ENDDO
 148 !     IF((IQ2+(N-1)/2).GT.NQ2) IQ2 = NQ2 - (N-1)/2
 149 !     IF((IQ2-(N-1)/2).LT.1) IQ2 = (N+1)/2
 150 !     JDUM = 0
 151 !     DO J = -(N-1)/2,(N-1)/2,1
 152 !     JDUM = JDUM +1
 153 !     IX2A(JDUM) = IQ2 + J
 154 !     ENDDO
 155
 156 !     Assign grid for interpolation. M, N -> order of polyN interpolation
 157       do I=1,M
 158          if(IX.ge.M/2.and.IX.le.(NX-M/2)) IX1A(I) = IX - M/2 + I
 159          if(IX.lt.M/2) IX1A(I) = I
 160          if(IX.gt.(NX-M/2)) IX1A(I) = (NX - M) + I
 161
 162 !     Check grids
 163          if(IX1A(I).le.0.or.IX1A(I).gt.NX) then
 164             write(6,*) "Error in grids! "
 165             write(6,*) "I, IXIA(I) = ",I, IX1A(I)
 166             call exit(-10)
 167          endif
 168       enddo
 169
 170       do J=1,N
 171          if(IQ2.ge.N/2.and.IQ2.le.(NQ2-N/2)) IX2A(J) = IQ2 - N/2 + J
 172          if(IQ2.lt.N/2) IX2A(J) = J
 173          if(IQ2.gt.(NQ2-N/2)) IX2A(J) = (NQ2 - N) + J
 174 !     Check grids
 175          if(IX2A(J).le.0.or.IX2A(J).gt.NQ2) then
 176             write(6,*) "Error in grids! "
 177             write(6,*) "J, IXIA(J) = ",J,IX2A(J)
 178             call exit(-10)
 179          endif
 180       enddo
 181
 182 !     Define points where to evaluate interpolation
 183 !     Choose between linear or logarithmic (x,Q2) interpolation
 184
 185       IF(X.LT.XCH)THEN
 186          X1=dlog(X)
 187       ELSE
 188          X1=X
 189       ENDIF
 190       X2=dlog(Q2)
 191
 192       DO IPDF = -6,7,1
 193
 194 !     Choose between linear or logarithmic (x,Q2) interpolation
 195
 196          DO I=1,M
 197             IF(X.LT.XCH)THEN
 198                X1A(I)= dlog(XG(IX1A(I),iset))
 199             ELSE
 200                X1A(I)= XG(IX1A(I),iset)
 201             ENDIF
 202             DO J=1,N
 203                X2A(J) = dlog(Q2G(IX2A(J),iset))
 204                YA(I,J) = XPDFEV(IX1A(I),IX2A(J),IPDF,KREP+offset)
 205             enddo
 206          enddo
 207
 208 !     2D polynomial interpolation
 209          call lh_polin2(x1a,x2a,ya,m,n,x1,x2,y,dy)
 210
 211          IF(IPDF.NE.7)THEN
 212             XPDF(IPDF) = y
 213          ELSE
 214             XPHOTON = y
 215          ENDIF
 216
 217       enddo
 218
 219       RETURN
 220
 221 !********************************************************
 222
 223 !********************************************************
 224       entry NNPDFINT20qedgetgrid(nset,ngridx,ngridq,gridx,gridq)
 225       do jx=1,nx
 226           gridx(jx)=xg(jx,nset)
 227       enddo
 228       do jq=1,nq2
 229           gridq(jq)=q2g(jq,nset)
 230       enddo
 231       ngridx=nx
 232       ngridq=nq2
 233       return
 234 !********************************************************
 235
 236
 237       ENTRY NNPDFINT20qedread(nset)
 238 !
 239       READ(1,*)nmem(nset),ndef(nset)
 240
 241 !     Set number of members
 242       call setnmem(nset,nmem)
 243
 244 !     Read the grid in x
 245       READ(1,*) nxx
 246       IF(NXX.NE.NX)WRITE(*,*)"WARNING CHANGE NX ACCORDING TO .LHgrid"
 247       DO ix = 1,nxx
 248          READ(1,*) xg(ix,nset)
 249       ENDDO
 250
 251 !     Read the grid in Q2
 252       READ(1,*) nqq2
 253       IF(NQQ2.NE.NQ2)WRITE(*,*)"WARNING CHANGE NQ2 ACCORDING TO .LHgrid"
 254       READ(1,*) q2min
 255       DO iq2 = 1,nqq2
 256          READ(1,*) q2g(iq2,nset)
 257       ENDDO
 258
 259 !     Read the number of replicas
 260       READ(1,*) NREP
 261
 262 !     Read the evolved xpdf grid for each replica
 263       offset = 101*(nset-1)
 264       DO K=0,NREP
 265          DO IX=1,NX
 266             DO IQ2=1,NQ2
 267                READ(1,*) ( xpdfev(ix,iq2,ipdf,k+offset), ipdf=-6,7,1 )
 268             ENDDO
 269          ENDDO
 270       ENDDO
 271 !
 272       RETURN
 273
 274 !********************************************************
 275
 276       ENTRY NNPDFINT20qedalfa(alfas,QAS)
 277 !
 278       QQ = QAS
 279       alfas =  alphaNNPDF(QQ)
 280 !
 281       RETURN
 282
 283 !********************************************************
 284
 285       ENTRY NNPDFINT20qedinit(nset,Eorder,Q2fit)
 286 !
 287       IMODEV = 0
 288       IVFN = 1
 289       ITMC = 1
 290 !
 291       CALL GetOrderPDFM(nset,order)
 292       IPT = order
 293 !
 294       CALL GetQ2fitM(nset,QQ2)
 295       Q2fit = QQ2
 296       Q20   = QQ2
 297 !
 298       call GetQmassM(nset,4,mass)
 299       QTH(4) = mass
 300       call GetQmassM(nset,5,mass)
 301       QTH(5) = mass
 302       call GetQmassM(nset,6,mass)
 303       QTH(6) = mass
 304 !
 305       DO i = 4,6
 306          q2th(i) = qth(i)**2d0
 307       ENDDO
 308 !
 309                                     ! added for filling Fparm->asref
 310       call initEVOLVEpdf(nset,mmem)
 311       CALL GetAlfas(nset,alfas0,Q0)
 312       asref = alfas0
 313       q2ref = q0**2d0
 314 !
 315       RETURN
 316
 317 !********************************************************
 318 !
 319       entry NNPDFINT20qedpdf(nset)
 320       pdfmem = nset
 321       IF (pdfmem.LT.0) THEN
 322          WRITE(*,*) 'NNPDF set:'
 323          WRITE(*,*) 'PDF member out of range:'
 324          WRITE(*,*) 'member = ',pdfmem
 325          STOP
 326       ENDIF
 327       RETURN
 328 !
 329       END
 330
 331
 332 !     THE ROUTINES FOR THE POLYNOMIAL INTERPOLATION ARE INSIDE
 333 !     THE wrapNNPDFgrid.f ROUTINE
 334