LHAPDF veraion 5.9.1
[u/mrichter/AliRoot.git] / 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