← Back to branch summary

~alifson/chiralityflow/trunk

« back to all changes in this revision

Viewing changes to Template/NLO/MCatNLO/srcCommon/mcatnlo_mlmtopdf.f

  • Committer: andrew.lifson at lu
  • Date: 2021-09-01 15:34:39 UTC
  • Revision ID: andrew.lifson@thep.lu.se-20210901153439-7fasjhav4cp4m88r
testing a new repository of a madgraph folder

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/NLO/MCatNLO/srcCommon/mcatnlo_mlmtopdf.f
 
1
c Interface to pdflib routines, version 4.00 and beyond (PDF set
 
2
c is identified by three parameters, NPTYPE, NGROUP, NSET).
 
3
c Includes also electron distribution functions of the jet package
 
4
c Modified on 8/11/01, since pdflib version 8.04 has in one case
 
5
c 100 pdf sets for one author group. Bug relevant to special treatment
 
6
c of AFG set fixed on 24/10/05 (603->6003)
 
7
c
 
8
      subroutine mlmpdf(mode,ih,q2,x,fx,nf)
 
9
      implicit real * 8 (a-h,o-z)
 
10
      real * 4 q2,x,fx(-nf:nf)
 
11
      common/trans/nptype,ngroup,nset
 
12
c set print flag. Set iflprt=2 in order to print the content
 
13
c of the common blocks w50511, w50512, and w50513
 
14
      common/w50510/iflprt
 
15
c x*pdf(x) in PDG format
 
16
      dimension fxp(-6:6)
 
17
c used by pdfset
 
18
      dimension val(20)
 
19
c used by pdfset
 
20
      character * 20 parm(20)
 
21
      logical ini
 
22
      data ini/.true./,imode/0/
 
23
c
 
24
      if(x.le.0.or.x.ge.1) then
 
25
        do j=-nf,nf
 
26
          fx(j) = 0
 
27
        enddo
 
28
        return
 
29
      endif
 
30
      if(ih.ne.5)then
 
31
c incoming particle is not an electron: use PDFLIB
 
32
        if(ini.or.imode.ne.mode) then
 
33
          ini = .false.
 
34
          imode = mode
 
35
c pass from our conventions to PDFLIB conventions for parameter settings. 
 
36
c See subroutine newmode for comments on the conventions adopted
 
37
          call newmode(ih,mode)
 
38
          parm(1) = 'NPTYPE'
 
39
          val (1) =  nptype
 
40
          parm(2) = 'NGROUP'
 
41
          val (2) =  ngroup
 
42
          parm(3) = 'NSET'
 
43
          val (3) =  nset
 
44
          call pdfset(parm,val)
 
45
        endif
 
46
        xd = dble(x)
 
47
        qd = dble(sqrt(q2))
 
48
        call pftopdg(xd,qd,fxp)
 
49
c in our conventions 1=up, 2=down; PDFLIB 1=down, 2=up. With
 
50
c f(1)<-->f(2) we mean also f(-1)<-->f(-2)
 
51
c in the following lines, deals with particles only (no antiparticles)
 
52
c proton(ih=1) ==> f(1)<-->f(2)
 
53
c neutron(ih=2) ==> no action (f(1)<-->f(2) for PDFLIB convention and
 
54
c    f(1)<-->f(2) for isospin symmetry (u_proton=d_neutron....)
 
55
c pion+(ih=3) ==> f(2)<-->f(-2), since PDFLIB has d=u=q_v+q_sea, 
 
56
c    ubar=dbar=q_sea
 
57
c photon(ih=4) ==> f(-1)<-->f(-2) and f(i)=f(-i)/2, i=1,2 since PDFLIB 
 
58
c    has f(i)=2*f(-i), and f(1)<-->f(2)
 
59
c Notice that in the jet package pions and neutrons are not used. If
 
60
c selected, they are rejected by the routine pdfpar. This routine
 
61
c is therefore a completely general interface with PDFLIB
 
62
        if(abs(ih).eq.1) then
 
63
          tmp     = fxp(1)
 
64
          fxp(1)  = fxp(2)
 
65
          fxp(2)  = tmp
 
66
          tmp     = fxp(-1)
 
67
          fxp(-1) = fxp(-2)
 
68
          fxp(-2) = tmp
 
69
        elseif(abs(ih).eq.2) then
 
70
          continue
 
71
        elseif(abs(ih).eq.3) then
 
72
          tmp = fxp(-2)
 
73
          fxp(-2) = fxp(2)
 
74
          fxp(2)  = tmp
 
75
        elseif(abs(ih).eq.4) then
 
76
          tmp     = fxp(-1)
 
77
          fxp(-1) = fxp(-2)
 
78
          fxp(-2) = tmp
 
79
          fxp(1)=fxp(-1)
 
80
          fxp(2)=fxp(-2)
 
81
c this is (p+n)/2
 
82
        elseif(ih.eq.0) then
 
83
          va  = (fxp(1)+fxp(2))/2
 
84
          sea = (fxp(-1)+fxp(-2))/2
 
85
          fxp(1)  = va
 
86
          fxp(2)  = va
 
87
          fxp(-1) = sea
 
88
          fxp(-2) = sea
 
89
        else
 
90
          write(*,*) ' ih was', ih,' instead of 0, +-1, +-2, +-3, or 4'
 
91
          stop
 
92
        endif
 
93
c for particles, ich=1, for antiparticles, ich=-1
 
94
        if(ih.lt.0) then
 
95
          ich = -1
 
96
        else
 
97
          ich = 1
 
98
        endif
 
99
c divide by x and exchange q with qbar in the case of antiparticles
 
100
        do j=-nf,nf
 
101
          fx(j) = sngl(fxp(j*ich)/xd)
 
102
        enddo
 
103
c correct for a bug in PDFLIB specific to AFG densities in the photon
 
104
        if(ih.eq.4.and.mode.eq.6003)then
 
105
          do j=-nf,nf
 
106
            if(j.ne.0)fx(j)=fx(j)/2.d0
 
107
          enddo
 
108
        endif
 
109
      else
 
110
c incoming "hadron" is an electron
 
111
        if(ini.or.imode.ne.mode) then
 
112
          ini = .false.
 
113
          imode = mode
 
114
          iset=mode-50
 
115
        endif
 
116
        if(iset.eq.1) then
 
117
          call elpdf_lac1(q2,x,fx,nf)
 
118
        elseif(iset.eq.2) then
 
119
          call elpdf_grv(q2,x,fx,nf)
 
120
        elseif(iset.eq.3) then
 
121
          call elpdf_user(q2,x,fx,nf)
 
122
        else
 
123
          write(*,*)'Electron set non implemented'
 
124
          stop
 
125
        endif
 
126
      endif
 
127
      return
 
128
      end
 
129
 
 
130
 
 
131
      subroutine pdfpar(mode,ih,xlam,scheme,iret)
 
132
      implicit real * 8 (a-h,o-z)
 
133
      common/trans/nptype,ngroup,nset
 
134
      common/w50512/qcdl4,qcdl5
 
135
      character*10 pdfver(3)
 
136
      common/w505190/pdfver
 
137
      dimension val(20)
 
138
      character * 20 parm(20)
 
139
      character * 2 scheme
 
140
      logical ini
 
141
      data ini/.true./
 
142
      logical nopions
 
143
      data nopions/.false./
 
144
c iret#0 when problem occur
 
145
      iret = 0
 
146
      if(ini)then
 
147
        write(*,*)'This is an interface to PDFLIB'
 
148
        write(*,*)'PDFLIB version number:',pdfver(1)
 
149
        write(*,*)'Released:',pdfver(2)
 
150
        write(*,*)'On:',pdfver(3)
 
151
        ini = .false.
 
152
      endif
 
153
      if(nopions)then
 
154
        if(abs(ih).ne.1.and.ih.ne.4.and.ih.ne.5)then
 
155
          write(*,*) ' hadron tpye ',ih,' not implemented'
 
156
          iret=1
 
157
          return
 
158
        endif
 
159
      endif
 
160
c fake values. If kept, the main program crashes
 
161
      scheme='XX'
 
162
      xlam=0.0
 
163
c incoming particle is not an electron: use PDFLIB
 
164
      if(ih.ne.5)then
 
165
c the scheme of the PDFLIB set is not given in any common block.
 
166
c Set it by hand in the main program
 
167
        scheme = '**'
 
168
        call newmode(ih,mode)
 
169
        parm(1) = 'NPTYPE'
 
170
        val (1) =  nptype
 
171
        parm(2) = 'NGROUP'
 
172
        val (2) =  ngroup
 
173
        parm(3) = 'NSET'
 
174
        val (3) =  nset
 
175
c print the relevant parameters for the PDF set chosen
 
176
c        iflprt = 2
 
177
c set the parameters
 
178
        call pdfset(parm,val)
 
179
c Lambda_QCD_5, as given by PDFLIB
 
180
        xlam = qcdl5
 
181
      else
 
182
c incoming particle is an electron
 
183
        if(mode.eq.51)then
 
184
          scheme='MS'
 
185
          xlam=.130
 
186
        elseif(mode.eq.52)then
 
187
          scheme='DG'
 
188
          xlam=.130
 
189
        elseif(mode.eq.53)then
 
190
          scheme='**'
 
191
          xlam=.001
 
192
        else
 
193
          write(*,*)'Electron set not implemented'
 
194
          iret=1
 
195
        endif
 
196
      endif
 
197
      return
 
198
      end
 
199
 
 
200
 
 
201
      subroutine prntsf
 
202
c     prints details of the structure function sets
 
203
c
 
204
      write(*,100)                             
 
205
     #  '  For nucleons, pions and photons, enter the set number'
 
206
     # ,'  '
 
207
     # ,'              Set # = 1000*NGROUP+NSET'
 
208
     # ,'  '
 
209
     # ,'  where NGROUP and NSET are the parameters of the PDFLIB'
 
210
     # ,'  version you are using (please read the user manual of'
 
211
     # ,'  your version of PDFLIB to check the listing)'
 
212
     # ,'  '
 
213
      write(*,100)                             
 
214
     #  '  PDFLIB is not maintained by CERN any longer, and the last'
 
215
     # ,'  release (version 8.04) included up to MRST99 and CTEQ5 sets'
 
216
     # ,'  For more recent sets, we use the following labellings'
 
217
     # ,'  '
 
218
     # ,'Group #   Set #        Set'
 
219
     # ,'  '
 
220
      write(*,100)                             
 
221
     #  '  3        101      MRST 2001 NNLO average'
 
222
     # ,'  3        102      MRST 2001 NNLO fast'
 
223
     # ,'  3        103      MRST 2001 NNLO slow'
 
224
     # ,'  3        104      MRST 2001 NNLO jet'
 
225
     # ,'  3        105      MRST 2001 best fit'
 
226
     # ,'  3        106      MRST 2001 low as'
 
227
     # ,'  3        107      MRST 2001 high as'
 
228
     # ,'  3        108      MRST 2001 jet fit'
 
229
     # ,'  3        109      MRST 2001 leading order'
 
230
     # ,'  3        110      MRST 2002 best fit'
 
231
     # ,'  3        111      MRST 2002 NNLO'
 
232
     # ,'  3      112-142    MRST 2002 error sets'
 
233
     # ,'  '
 
234
      write(*,100)                             
 
235
     #  '  4         55      CTEQ5M1 parametrized version'
 
236
     # ,'  4         56      CTEQ6M'
 
237
     # ,'  4         57      CTEQ6D'
 
238
     # ,'  4         58      CTEQ6L'
 
239
     # ,'  4       59-98     CTEQ6 error sets'
 
240
     # ,'  '
 
241
      write(*,100)                             
 
242
     #  ' 10          1      Alekhin LO   nominal ffn'
 
243
     #, ' 10          2      Alekhin LO   nominal vfn'
 
244
     #, ' 10          3      Alekhin LO   mc=1.75 ffn'
 
245
     #, ' 10          4      Alekhin LO   mc=1.75 vfn'   
 
246
     #, ' 10          5      Alekhin LO   ss      ffn'   
 
247
     #, ' 10          6      Alekhin LO   ss      vfn'   
 
248
     #, ' 10          7      Alekhin NL   nominal ffn'
 
249
     #, ' 10          8      Alekhin NL   nominal vfn'
 
250
     #, ' 10          9      Alekhin NL   mc=1.75 ffn'
 
251
     #, ' 10         10      Alekhin NL   mc=1.75 vfn'   
 
252
     #, ' 10         11      Alekhin NL   ss      ffn'   
 
253
     #, ' 10         12      Alekhin NL   ss      vfn'   
 
254
     #, ' 10         13      Alekhin NNL  nominal ffn'
 
255
     #, ' 10         14      Alekhin NNL  nominal vfn'
 
256
     #, ' 10         15      Alekhin NNL  mc=1.75 ffn'
 
257
     #, ' 10         16      Alekhin NNL  mc=1.75 vfn'   
 
258
     #, ' 10         17      Alekhin NNL  ss      ffn'   
 
259
     #, ' 10         18      Alekhin NNL  ss      vfn'   
 
260
     #, ' 10         19      Alekhin NNL  slow ev ffn'   
 
261
     #, ' 10         20      Alekhin NNL  slow ev vfn'   
 
262
      write(*,100)                             
 
263
     #  '  For electrons, use'
 
264
     # ,'  '
 
265
     # ,'Set #        Set'
 
266
     # ,'  '
 
267
     # ,'  51         LAC1'
 
268
     # ,'  52         GRV-HO'
 
269
     # ,'  53         user defined'
 
270
 100  format(1x,a,100(/,1x,a))
 
271
      return
 
272
      end
 
273
 
 
274
 
 
275
      subroutine newmode(ih,mode)
 
276
c This subroutine converts our conventions for the identification
 
277
c of a set of PDFs into PDFLIB conventions. We use
 
278
c
 
279
c                     JET PACKAGE             PDFLIB
 
280
c 
 
281
c Particle type 
 
282
c
 
283
c  nucleons           -2,-1,0,1,2                1
 
284
c  pions                  -3,3                   2
 
285
c  photons                  4                    3
 
286
c 
 
287
c Given the particle type, our package uses a single number (MODE)
 
288
c to identify to PDF set, while PDFLIB uses two numbers (NGROUP and NSET).
 
289
c The value of MODE which corresponds to the values (NGROUP,NSET)
 
290
c is by definition given by
 
291
c
 
292
c                        MODE=1000*NGROUP+NSET
 
293
c
 
294
c This is working as long as every group has less the 1000 PDFs sets....
 
295
      implicit real * 8 (a-h,o-z)
 
296
      common/trans/nptype,ngroup,nset
 
297
c
 
298
      if(abs(ih).le.2)then
 
299
        nptype=1
 
300
      elseif(abs(ih).eq.3)then
 
301
        nptype=2
 
302
      elseif(ih.eq.4)then
 
303
        nptype=3
 
304
      else
 
305
        write(6,*)'hadron type not implemented in PDFLIB'
 
306
        stop
 
307
      endif
 
308
      ngroup=int(dfloat(mode)/1000.e0)
 
309
      nset=mode-1000*ngroup
 
310
      return
 
311
      end
 
312
 
 
313
 
 
314
      subroutine pdftomlm(ipdfih,ipdfgroup,ipdfndns,ihmlm,ndnsmlm)
 
315
c Performs the inverse operation of newmode
 
316
      implicit real * 8 (a-h,o-z)
 
317
c
 
318
      if(ipdfih.eq.1)then
 
319
        ihmlm=1
 
320
      elseif(ipdfih.eq.2)then
 
321
        ihmlm=3
 
322
      elseif(ipdfih.eq.3)then
 
323
        ihmlm=4
 
324
      else
 
325
        write(*,*)'Wrong hadron type in PDFLIB:',ipdfih
 
326
        stop
 
327
      endif
 
328
      ndnsmlm=1000*ipdfgroup+ipdfndns
 
329
      return
 
330
      end
 
331
 
 
332
 
 
333
      subroutine setlhacblk(strin)
 
334
c It should not be called when this file is linked
 
335
      write(*,*)'Call to setlhacblk must not occur'
 
336
      stop
 
337
      end
 
338
 
 
339
 
 
340
      subroutine getlamlha(xlam,xlamlha)
 
341
c It should not be called when this file is linked
 
342
      write(*,*)'Call to getlamlha must not occur'
 
343
      stop
 
344
      end
 
345
 
 
346
 
 
347
C------- ALPHA QCD -------------------------------------
 
348
c Program to calculate alfa strong with nf flavours,
 
349
c as a function of lambda with 5 flavors.
 
350
c The value of alfa is matched at the thresholds q = mq.
 
351
c When invoked with nf < 0 it chooses nf as the number of
 
352
c flavors with mass less then q.
 
353
c
 
354
      function alfas(q2,xlam,inf)
 
355
      implicit real * 8 (a-h,o-z)
 
356
      data olam/0.d0/,pi/3.14159d0/
 
357
      data xmb/5.d0/,xmc/1.5d0/
 
358
      if(xlam.ne.olam) then
 
359
        olam = xlam
 
360
        b5  = (33-2*5)/pi/12
 
361
        bp5 = (153 - 19*5) / pi / 2 / (33 - 2*5)
 
362
        b4  = (33-2*4)/pi/12
 
363
        bp4 = (153 - 19*4) / pi / 2 / (33 - 2*4)
 
364
        b3  = (33-2*3)/pi/12
 
365
        bp3 = (153 - 19*3) / pi / 2 / (33 - 2*3)
 
366
        xlc = 2 * log(xmc/xlam)
 
367
        xlb = 2 * log(xmb/xlam)
 
368
        xllc = log(xlc)
 
369
        xllb = log(xlb)
 
370
        c45  =  1/( 1/(b5 * xlb) - xllb*bp5/(b5 * xlb)**2 )
 
371
     #        - 1/( 1/(b4 * xlb) - xllb*bp4/(b4 * xlb)**2 )
 
372
        c35  =  1/( 1/(b4 * xlc) - xllc*bp4/(b4 * xlc)**2 )
 
373
     #        - 1/( 1/(b3 * xlc) - xllc*bp3/(b3 * xlc)**2 ) + c45
 
374
      endif
 
375
      q   = sqrt(q2)
 
376
      xlq = 2 * log( q/xlam )
 
377
      xllq = log( xlq )
 
378
      nf = inf
 
379
      if( nf .lt. 0) then
 
380
        if( q .gt. xmb ) then
 
381
          nf = 5
 
382
        elseif( q .gt. xmc ) then
 
383
          nf = 4
 
384
        else
 
385
          nf = 3
 
386
        endif
 
387
      endif
 
388
      if    ( nf .eq. 5 ) then
 
389
        alfas = 1/(b5 * xlq) -  bp5/(b5 * xlq)**2 * xllq
 
390
      elseif( nf .eq. 4 ) then
 
391
        alfas = 1/( 1/(1/(b4 * xlq) - bp4/(b4 * xlq)**2 * xllq) + c45 )
 
392
      elseif( nf .eq. 3 ) then
 
393
        alfas = 1/( 1/(1/(b3 * xlq) - bp3/(b3 * xlq)**2 * xllq) + c35 )
 
394
      else
 
395
        print *,'error in alfa: unimplemented # of light flavours',nf
 
396
        stop
 
397
      endif
 
398
      return
 
399
      end