← Back to branch summary

~maddevelopers/mg5amcnlo/3.0.2-alpha0

« back to all changes in this revision

Viewing changes to Template/Source/is-sud.f

Added Template and HELAS into bzr

Show diffs side-by-side

added added

removed removed

Lines of Context:
Template/Source/is-sud.f
 
1
      function calcissud(kkfl,xx1,xx2,pt2min,pt2,nfnevl,relerr,ifail)
 
2
 
 
3
c...arguments
 
4
      integer kkfl
 
5
      double precision calcissud,xx1,xx2,pt2min,pt2
 
6
c...global variables
 
7
      double precision bthres,alam5,alam4,s,x1,x2
 
8
      integer kfl
 
9
      common/sudpars/bthres,alam5,alam4,s,x1,x2,kfl
 
10
      data alam5,alam4/2.11384241947499864d-2,3.6864E-002/
 
11
      data bthres,s/22.3109,196d6/
 
12
c...local variables
 
13
      integer n,minpts,maxpts,iwk,nfnevl,ifail
 
14
      parameter(n=2,minpts=5,maxpts=50000,iwk=10298) !iwk=7*(1+maxpts/34)
 
15
      double precision fsud,a(n),b(n),eps,wk(iwk),result,relerr
 
16
      parameter(eps=0.001)
 
17
      external fsud
 
18
 
 
19
      x1=xx1
 
20
      x2=xx2
 
21
      kfl=kkfl
 
22
 
 
23
c...limits for pt2
 
24
      a(1)=log(pt2min)
 
25
      if(abs(kfl).eq.5) a(1)=log(max(pt2min,bthres))
 
26
      b(1)=log(pt2)
 
27
c...limits for z
 
28
      a(2)=x1
 
29
      b(2)=1
 
30
 
 
31
      if(pt2.le.exp(a(1)))then
 
32
        calcissud=1
 
33
        nfnevl=0
 
34
        relerr=0
 
35
        ifail=0
 
36
        return
 
37
      endif
 
38
 
 
39
      if(abs(kfl).eq.5.and.x1.gt.0.6)then
 
40
        calcissud=0d0
 
41
        nfnevl=0
 
42
        relerr=0
 
43
        ifail=0
 
44
        return
 
45
      endif
 
46
 
 
47
c...perform integration
 
48
      call dadmul(fsud,n,a,b,minpts,maxpts,eps,wk,iwk,
 
49
     $   result,relerr,nfnevl,ifail)
 
50
 
 
51
      calcissud=exp(-result)
 
52
      return
 
53
      end
 
54
 
 
55
 
 
56
C...FSUD is the integrand called by dadmul 
 
57
      double precision function fsud(n,x)
 
58
      implicit none
 
59
 
 
60
c...arguments
 
61
      integer n
 
62
      double precision x(*)
 
63
 
 
64
c...global variables
 
65
      double precision bthres,alam5,alam4,s,x1,x2
 
66
      integer kfl
 
67
      common/sudpars/bthres,alam5,alam4,s,x1,x2,kfl
 
68
 
 
69
c...functions
 
70
      double precision pdg2pdf
 
71
      external pdg2pdf
 
72
      
 
73
c...local variables
 
74
      double precision b5,b4
 
75
      data b5,b4/3.8333333333333335,4.1666666666666/
 
76
      double precision pt2,pt,z,zmax,zint,s12
 
77
      DOUBLE PRECISION XPQ1(-7:7),XPQ2(-7:7)
 
78
      DATA XPQ1,XPQ2/30*0d0/
 
79
      double precision qm2(5),qmrel
 
80
      data qm2/3*0,2.2725,22.3109/
 
81
      double precision ptold,xiold,kflold,xpqold
 
82
      integer kfla,i
 
83
 
 
84
      pt2=exp(x(1))
 
85
      if(pt2.gt.s) then
 
86
        print *,'Error: Probably forgot log() for pt2 limits!'
 
87
        stop
 
88
      endif
 
89
      pt=dsqrt(pt2)
 
90
      z=x(2)
 
91
      s12=s*x1*x2
 
92
      zmax=1d0+pt2/(2d0*s12)-sqrt(pt2**2/(4d0*s12**2)+pt2/s12)
 
93
      kfla=iabs(kfl)
 
94
 
 
95
      fsud=0d0
 
96
      if(z.gt.zmax.or.z.le.x1) return
 
97
 
 
98
c   PT integral: 1/pt2*alpha_s(pt)/2pi
 
99
      if(pt2.gt.bthres) then
 
100
        fsud=1d0/(b5*log(pt2/alam5))
 
101
      else! if(pt2.gt.2.2725000000000000)then
 
102
        fsud=1d0/(b4*log(pt2/alam4))
 
103
      endif
 
104
 
 
105
c   Color coherence suppression
 
106
      if(.false.)then
 
107
        fsud=fsud*min(1d0,s12/4d0/pt2)
 
108
      endif
 
109
c   Z integral: need pdf ratio and splitting function
 
110
      zint=0
 
111
c      CALL PFTOPDG(1,x1/z,DSQRT(pt2),XPQ2)
 
112
      if(kfla.lt.6.and.kfl.ne.0)then
 
113
        XPQ1(kfl)=max(pdg2pdf(1,kfl,x1,pt),1d-20)
 
114
        XPQ2(kfl)=pdg2pdf(1,kfl,x1/z,pt)
 
115
        XPQ2(0)=pdg2pdf(1,21,x1/z,pt)
 
116
        qmrel=0
 
117
        if(qm2(kfla).gt.0)
 
118
     $     qmrel=2*z*(1-z)*qm2(kfla)/pt2
 
119
        zint=zint+XPQ2(kfl)/(z*XPQ1(kfl))*4d0/3d0*
 
120
     $     ((1+z**2)/(1-z)-qmrel)
 
121
        zint=zint+XPQ2(0)/(z*XPQ1(kfl))*0.5d0*(z**2+(1-z)**2+qmrel)
 
122
      elseif(kfl.eq.21.or.kfl.eq.0)then
 
123
        XPQ1(0)=pdg2pdf(1,0,x1,pt)
 
124
c        CALL PFTOPDG(1,x1/z,pt,XPQ2)
 
125
        do i=-2,3
 
126
          XPQ2(i)=pdg2pdf(1,i,x1/z,pt)
 
127
          if(i.ne.0)
 
128
     $       zint=zint+XPQ2(i)
 
129
        enddo
 
130
        zint=zint+XPQ2(3) ! for anti-s dist
 
131
        zint=zint/(z*XPQ1(0))*4d0/3d0*(1+(1-z)**2)/z
 
132
        zint=zint+XPQ2(0)/(z*XPQ1(0))*6*(1-z*(1-z))**2/(z*(1-z))
 
133
      endif
 
134
      fsud=fsud*zint
 
135
 
 
136
      return
 
137
      end
 
138
 
 
139
 
 
140
      double precision function pdg2pdf(ih,ipdg,x,xmu)
 
141
c***************************************************************************
 
142
c     Based on pdf.f, wrapper for calling the pdf of MCFM
 
143
c***************************************************************************
 
144
      implicit none
 
145
c
 
146
c     Arguments
 
147
c
 
148
      DOUBLE  PRECISION x,xmu
 
149
      INTEGER IH,ipdg,ipart
 
150
C
 
151
C     Include
 
152
C
 
153
      include 'PDF/pdf.inc'
 
154
C      
 
155
      double precision Ctq3df,Ctq4Fn,Ctq5Pdf,Ctq6Pdf,Ctq5L
 
156
      integer mode,Irt
 
157
 
 
158
      ipart=ipdg
 
159
      if(ipart.eq.21) ipart=0
 
160
 
 
161
      if (pdlabel(1:5) .eq. 'cteq3') then
 
162
C     
 
163
         if (pdlabel .eq. 'cteq3_m') then
 
164
            mode=1
 
165
         elseif (pdlabel .eq. 'cteq3_l') then
 
166
            mode=2
 
167
         elseif (pdlabel .eq. 'cteq3_d') then
 
168
            mode=3
 
169
         endif
 
170
 
 
171
         
 
172
         if(iabs(ipart).ge.1.and.iabs(ipart).le.2)
 
173
     $      ipart=sign(3-iabs(ipart),ipart)
 
174
 
 
175
         pdg2pdf=Ctq3df(mode,ipart,x,xmu,Irt)/x
 
176
 
 
177
         if(ipdg.ge.1.and.ipdg.le.2)
 
178
     $      pdg2pdf=pdg2pdf+Ctq3df(mode,-ipart,x,xmu,Irt)/x
 
179
 
 
180
C     
 
181
      elseif (pdlabel(1:5) .eq. 'cteq4') then
 
182
C     
 
183
         if (pdlabel .eq. 'cteq4_m') then
 
184
            mode=1
 
185
         elseif (pdlabel .eq. 'cteq4_d') then
 
186
            mode=2
 
187
         elseif (pdlabel .eq. 'cteq4_l') then
 
188
            mode=3
 
189
         elseif (pdlabel .eq. 'cteq4a1') then
 
190
            mode=4
 
191
         elseif (pdlabel .eq. 'cteq4a2') then
 
192
            mode=5
 
193
         elseif (pdlabel .eq. 'cteq4a3') then
 
194
            mode=6
 
195
         elseif (pdlabel .eq. 'cteq4a4') then
 
196
            mode=7
 
197
         elseif (pdlabel .eq. 'cteq4a5') then
 
198
            mode=8
 
199
         elseif (pdlabel .eq. 'cteq4hj') then
 
200
            mode=9
 
201
         elseif (pdlabel .eq. 'cteq4lq') then
 
202
            mode=10
 
203
         endif
 
204
         
 
205
         if(iabs(ipart).ge.1.and.iabs(ipart).le.2)
 
206
     $      ipart=sign(3-iabs(ipart),ipart)
 
207
 
 
208
         pdg2pdf=Ctq4Fn(mode,ipart,x,xmu)
 
209
C
 
210
      elseif (pdlabel .eq. 'cteq5l1') then
 
211
C
 
212
         if(iabs(ipart).ge.1.and.iabs(ipart).le.2)
 
213
     $      ipart=sign(3-iabs(ipart),ipart)
 
214
 
 
215
         pdg2pdf=Ctq5L(ipart,x,xmu)
 
216
C         
 
217
      elseif ((pdlabel(1:5) .eq. 'cteq5') .or. 
 
218
     .        (pdlabel(1:4) .eq. 'ctq5')) then
 
219
C         
 
220
         if(iabs(ipart).ge.1.and.iabs(ipart).le.2)
 
221
     $      ipart=sign(3-iabs(ipart),ipart)
 
222
 
 
223
         pdg2pdf=Ctq5Pdf(ipart,x,xmu)
 
224
C                  
 
225
      elseif (pdlabel(1:5) .eq. 'cteq6') then
 
226
C         
 
227
         if(iabs(ipart).ge.1.and.iabs(ipart).le.2)
 
228
     $      ipart=sign(3-iabs(ipart),ipart)
 
229
 
 
230
         pdg2pdf=Ctq6Pdf(ipart,x,xmu)
 
231
 
 
232
      endif      
 
233
 
 
234
      end
 
235