~maddevelopers/mg5amcnlo/laser_beam_option

« back to all changes in this revision

Viewing changes to Template/Source/PDF/pdg2pdf_lhapdf.f

  • Committer: olivier Mattelaer
  • Date: 2013-03-05 18:52:27 UTC
  • mfrom: (239.1.31 mg5_1_5_8)
  • Revision ID: olivier.mattelaer@uclouvain.be-20130305185227-h1358ux804s2wc26
Pass to 1.5.8

Show diffs side-by-side

added added

removed removed

Lines of Context:
7
7
c     Arguments
8
8
c
9
9
      DOUBLE  PRECISION x,xmu
10
 
      INTEGER IH,ipdg,ipart
 
10
      INTEGER IH,ipdg
11
11
C
12
12
C     Include
13
13
C
14
14
      include 'pdf.inc'
15
15
C      
16
 
      double precision Ctq3df,Ctq4Fn,Ctq5Pdf,Ctq6Pdf,Ctq5L
17
 
      integer mode,Irt,ihlast
18
 
      double precision xlast,xmulast,pdflast(-7:7)
19
 
      save ihlast,xlast,xmulast,pdflast
 
16
      integer i,j,ihlast(2),ipart,iporg,ireuse,imemlast(2),iset,imem
 
17
      double precision xlast(2),xmulast(2),pdflast(-7:7,2)
 
18
      save ihlast,xlast,xmulast,pdflast,imemlast
 
19
      data ihlast/2*-99/
 
20
      data xlast/2*-99d9/
 
21
      data xmulast/2*-99d9/
 
22
      data pdflast/30*-99d9/
 
23
      data imemlast/2*-99/
 
24
 
 
25
c     Make sure we have a reasonable Bjorken x. Note that even though
 
26
c     x=0 is not reasonable, we prefer to simply return pdg2pdf=0
 
27
c     instead of stopping the code, as this might accidentally happen.
 
28
      if (x.eq.0d0) then
 
29
         pdg2pdf=0d0
 
30
         return
 
31
      elseif (x.lt.0d0 .or. x.gt.1d0) then
 
32
         write (*,*) 'PDF not supported for Bjorken x ', x
 
33
         open(unit=26,file='../../../error',status='unknown')
 
34
         write(26,*) 'Error: PDF not supported for Bjorken x ',x
 
35
         stop 1
 
36
      endif
20
37
 
21
38
      ipart=ipdg
22
39
      if(ipart.eq.21) ipart=0
23
40
      if(iabs(ipart).eq.22) ipart=7
 
41
      iporg=ipart
24
42
 
25
43
c     This will be called for any PDG code, but we only support up to 7
26
44
      if(iabs(ipart).gt.7)then
27
 
         pdg2pdf=0d0
 
45
         write(*,*) 'PDF not supported for pdg ',ipdg
 
46
         open(unit=26,file='../../../error',status='unknown')
 
47
         write(26,*) 'Error: PDF not supported for pdg ',ipdg
 
48
         stop 1
 
49
      endif
 
50
 
 
51
c     Determine the iset used in lhapdf
 
52
      call getnset(iset)
 
53
      if (iset.ne.1) then
 
54
         write (*,*) 'PDF not supported for Bjorken x ', x
 
55
         open(unit=26,file='../../../error',status='unknown')
 
56
         write(26,*) 'Error: PDF not supported for Bjorken x ',x
 
57
         stop 1
 
58
      endif
 
59
 
 
60
c     Determine the member of the set (function of lhapdf)
 
61
      call getnmem(iset,imem)
 
62
 
 
63
      ireuse = 0
 
64
      do i=1,2
 
65
c     Check if result can be reused since any of last two calls
 
66
         if (x.eq.xlast(i) .and. xmu.eq.xmulast(i) .and.
 
67
     $        imem.eq.imemlast(i) .and. ih.eq.ihlast(i)) then
 
68
            ireuse = i
 
69
         endif
 
70
      enddo
 
71
 
 
72
c     Reuse previous result, if possible
 
73
      if (ireuse.gt.0.and.pdflast(iporg,ireuse).ne.-99d9) then
 
74
         pdg2pdf=pdflast(iporg,ireuse)
28
75
         return 
29
76
      endif
30
77
 
31
 
      if(ih.eq.ihlast.and.x.eq.xlast.and.xmu.eq.xmulast)then
32
 
         pdg2pdf=pdflast(ipart);
33
 
      else
34
 
         call pftopdglha(ih,x,xmu,pdflast)
35
 
         ihlast=ih
36
 
         xlast=x
37
 
         xmulast=xmu
38
 
         pdg2pdf=pdflast(ipart);
 
78
c     Bjorken x and/or facrorization scale and/or PDF set are not
 
79
c     identical to the saved values: this means a new event and we
 
80
c     should reset everything to compute new PDF values. Also, determine
 
81
c     if we should fill ireuse=1 or ireuse=2.
 
82
      if (ireuse.eq.0.and.xlast(1).ne.-99d9.and.xlast(2).ne.-99d9)then
 
83
         do i=1,2
 
84
            xlast(i)=-99d9
 
85
            xmulast(i)=-99d9
 
86
            do j=-7,7
 
87
               pdflast(j,i)=-99d9
 
88
            enddo
 
89
            imemlast(i)=-99
 
90
            ihlast(i)=-99
 
91
         enddo
 
92
c     everything has been reset. Now set ireuse=1 to fill the first
 
93
c     arrays of saved values below
 
94
         ireuse=1
 
95
      else if(ireuse.eq.0.and.xlast(1).ne.-99d9)then
 
96
c     This is first call after everything has been reset, so the first
 
97
c     arrays are already filled with the saved values (hence
 
98
c     xlast(1).ne.-99d9). Fill the second arrays of saved values (done
 
99
c     below) by setting ireuse=2
 
100
         ireuse=2
 
101
      else if(ireuse.eq.0)then
 
102
c     Special: only used for the very first call to this function:
 
103
c     xlast(i) are initialized as data statements to be equal to -99d9
 
104
         ireuse=1
39
105
      endif
40
106
 
 
107
c     Call lhapdf and give the current values to the arrays that should
 
108
c     be saved
 
109
      call pftopdglha(ih,x,xmu,pdflast(-7,ireuse))
 
110
      xlast(ireuse)=x
 
111
      xmulast(ireuse)=xmu
 
112
      ihlast(ireuse)=ih
 
113
      imemlast(ireuse)=imem
 
114
c
 
115
      pdg2pdf=pdflast(ipart,ireuse);
 
116
      return
41
117
      end
42
118