← Back to branch summary

~maddevelopers/mg5amcnlo/2.5.4_run_py8_at_evtgen

« back to all changes in this revision

Viewing changes to Template/LO/Source/PDF/pdg2pdf_lhapdf6.f

  • Committer: olivier Mattelaer
  • Date: 2016-05-12 11:00:18 UTC
  • mfrom: (262.1.150 2.3.4)
  • Revision ID: olivier.mattelaer@uclouvain.be-20160512110018-sevb79f0wm4g8mpp
pass to 2.4.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/LO/Source/PDF/pdg2pdf_lhapdf6.f
 
1
      double precision function pdg2pdf(ih,ipdg,x,xmu)
 
2
c***************************************************************************
 
3
c     Based on pdf.f, wrapper for calling the pdf of MCFM
 
4
c***************************************************************************
 
5
      implicit none
 
6
c
 
7
c     Arguments
 
8
c
 
9
      DOUBLE  PRECISION x,xmu
 
10
      INTEGER IH,ipdg
 
11
C
 
12
C     Include
 
13
C
 
14
      include 'pdf.inc'
 
15
C      
 
16
      integer i,j,ihlast(20),ipart,iporg,ireuse,imemlast(20),iset,imem
 
17
     &     ,i_replace,ii,ipartlast(20)
 
18
      double precision xlast(20),xmulast(20),pdflast(20)
 
19
      save ihlast,xlast,xmulast,pdflast,imemlast,ipartlast
 
20
      data ihlast/20*-99/
 
21
      data ipartlast/20*-99/
 
22
      data xlast/20*-99d9/
 
23
      data xmulast/20*-99d9/
 
24
      data pdflast/20*-99d9/
 
25
      data imemlast/20*-99/
 
26
      data i_replace/20/
 
27
 
 
28
c     Make sure we have a reasonable Bjorken x. Note that even though
 
29
c     x=0 is not reasonable, we prefer to simply return pdg2pdf=0
 
30
c     instead of stopping the code, as this might accidentally happen.
 
31
      if (x.eq.0d0) then
 
32
         pdg2pdf=0d0
 
33
         return
 
34
      elseif (x.lt.0d0 .or. x.gt.1d0) then
 
35
         write (*,*) 'PDF not supported for Bjorken x ', x
 
36
         open(unit=26,file='../../../error',status='unknown')
 
37
         write(26,*) 'Error: PDF not supported for Bjorken x ',x
 
38
         stop 1
 
39
      endif
 
40
 
 
41
      ipart=ipdg
 
42
      if(iabs(ipart).eq.21) ipart=0
 
43
      if(iabs(ipart).eq.22) ipart=7
 
44
      iporg=ipart
 
45
 
 
46
c     This will be called for any PDG code, but we only support up to 7
 
47
      if(iabs(ipart).gt.7)then
 
48
         write(*,*) 'PDF not supported for pdg ',ipdg
 
49
         write(*,*) 'For lepton colliders, please set the lpp* '//
 
50
     $    'variables to 0 in the run_card'  
 
51
         open(unit=26,file='../../../error',status='unknown')
 
52
         write(26,*) 'Error: PDF not supported for pdg ',ipdg
 
53
         stop 1
 
54
      endif
 
55
 
 
56
c     Determine the iset used in lhapdf
 
57
      call getnset(iset)
 
58
      if (iset.ne.1) then
 
59
         write (*,*) 'PDF not supported for Bjorken x ', x
 
60
         open(unit=26,file='../../../error',status='unknown')
 
61
         write(26,*) 'Error: PDF not supported for Bjorken x ',x
 
62
         stop 1
 
63
      endif
 
64
 
 
65
c     Determine the member of the set (function of lhapdf)
 
66
      call getnmem(iset,imem)
 
67
 
 
68
      ireuse = 0
 
69
      ii=i_replace
 
70
      do i=1,20
 
71
c     Check if result can be reused since any of last twenty
 
72
c     calls. Start checking with the last call and move back in time
 
73
         if (ih.eq.ihlast(ii)) then
 
74
            if (ipart.eq.ipartlast(ii)) then
 
75
               if (x.eq.xlast(ii)) then
 
76
                  if (xmu.eq.xmulast(ii)) then
 
77
                     if (imem.eq.imemlast(ii)) then
 
78
                        ireuse = ii
 
79
                        exit
 
80
                     endif
 
81
                  endif
 
82
               endif
 
83
            endif
 
84
         endif
 
85
         ii=ii-1
 
86
         if (ii.eq.0) ii=ii+20
 
87
      enddo
 
88
 
 
89
c     Reuse previous result, if possible
 
90
      if (ireuse.gt.0) then
 
91
         if (pdflast(ireuse).ne.-99d9) then
 
92
            pdg2pdf=pdflast(ireuse)
 
93
            return 
 
94
         endif
 
95
      endif
 
96
 
 
97
c Calculated a new value: replace the value computed longest ago
 
98
      i_replace=mod(i_replace,20)+1
 
99
 
 
100
c     Call lhapdf and give the current values to the arrays that should
 
101
c     be saved
 
102
      call evolvepart(ipart,x,xmu,pdg2pdf)
 
103
      pdg2pdf=pdg2pdf/x
 
104
      pdflast(i_replace)=pdg2pdf
 
105
      xlast(i_replace)=x
 
106
      xmulast(i_replace)=xmu
 
107
      ihlast(i_replace)=ih
 
108
      imemlast(i_replace)=imem
 
109
c
 
110
      return
 
111
      end
 
112