← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

Viewing changes to vendor/CutTools/src/qcdloop/ffcb2.f

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/CutTools/src/qcdloop/ffcb2.f
 
1
*###[ ffcb2:
 
2
        subroutine ffcb2(cb2p,cb2d,cb1,cb0,ca0i,xp,xma,xmb,piDpj,ier)
 
3
***#[*comment:***********************************************************
 
4
*                                                                       *
 
5
*       Calculate    1   /    d^n Q Q(mu) Q(nu)                         *
 
6
*                 ------ | ------------------------                     *
 
7
*                 i pi^2 / (Q^2-ma^2)((Q+p)^2-mb^2)                     *
 
8
*                                                p mu                   *
 
9
*                   = B2p*p(mu)*p(nu) + B2d*delta    /p^2               *
 
10
*                                                p nu                   *
 
11
*                                                                       *
 
12
*       Input:  cb1        complex      vector twopoint function        *
 
13
*               cb0        complex      scalar twopoint function        *
 
14
*               ca0i(2)    complex      scalar onepoint function with   *
 
15
*                                               m1,m2                   *
 
16
*               xp         complex      p.p in B&D metric               *
 
17
*               xma,2      complex      m_1^2,m_2^2                     *
 
18
*               piDpj(3,3) complex      dotproducts between s1,s2,p     *
 
19
*               ier        integer      digits lost so far              *
 
20
*       Output: cb2p       complex      coefficient of p(mu)*p(nu)      *
 
21
*               cb2d       complex      coefficient of delta()/p^2      *
 
22
*               ier        integer      digits lost                     *
 
23
*                                                                       *
 
24
***#]*comment:***********************************************************
 
25
*  #[ declarations:
 
26
        implicit none
 
27
*
 
28
*       arguments
 
29
*
 
30
        integer ier
 
31
        DOUBLE COMPLEX xp,xma,xmb,piDpj(3,3)
 
32
        DOUBLE COMPLEX cb2p,cb2d,cb1,cb0,ca0i(2)
 
33
*
 
34
*       local variables
 
35
*
 
36
        integer ier0,i,j
 
37
        DOUBLE COMPLEX dmap,dmbp,dmamb,cc
 
38
        DOUBLE PRECISION absc
 
39
        DOUBLE PRECISION rm1,rm2,rp,rpiDpj(3,3),sprec
 
40
*
 
41
*       common blocks
 
42
*
 
43
        include 'ff.h'
 
44
*
 
45
*       statement function
 
46
*
 
47
        absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
 
48
*
 
49
*  #] declarations:
 
50
*  #[ real case:
 
51
        if ( DIMAG(xma).eq.0 .and. DIMAG(xmb).eq.0 ) then
 
52
            rm1 = DBLE(xma)
 
53
            rm2 = DBLE(xmb)
 
54
            rp  = DBLE(xp)
 
55
            do 20 j=1,3
 
56
                do 10 i=1,3
 
57
                    rpiDpj(i,j) = DBLE(piDpj(i,j))
 
58
   10           continue
 
59
   20       continue
 
60
            sprec = precx
 
61
            precx = precc
 
62
            call ffxb2(cb2p,cb2d,cb1,cb0,ca0i,rp,rm1,rm2,rpiDpj,
 
63
     +          ier)
 
64
            precx = sprec
 
65
            return
 
66
        endif
 
67
*  #] real case:
 
68
*  #[ get differences:
 
69
        ier0 = 0
 
70
        dmamb = xma - xmb
 
71
        dmap = xma - xp
 
72
        dmbp = xmb - xp
 
73
        if ( lwarn ) then
 
74
            if ( absc(dmamb) .lt. xloss*absc(xma) .and. xma .ne. xmb )
 
75
     +          call ffwarn(97,ier0,absc(dmamb),absc(xma))
 
76
            if ( absc(dmap) .lt. xloss*absc(xp) .and. xp .ne. xma )
 
77
     +          call ffwarn(98,ier0,absc(dmap),absc(xp))
 
78
            if ( absc(dmbp) .lt. xloss*absc(xp) .and. xp .ne. xmb )
 
79
     +          call ffwarn(99,ier0,absc(dmbp),absc(xp))
 
80
        endif
 
81
*  #] get differences:
 
82
*  #[ call ffcb2a:
 
83
        call ffcb2a(cb2p,cb2d,cb1,cb0,ca0i,xp,xma,xmb,dmap,dmbp,dmamb,
 
84
     +          piDpj,ier)
 
85
*  #] call ffcb2a:
 
86
*###] ffcb2:
 
87
        end
 
88
*###[ ffcb2a:
 
89
        subroutine ffcb2a(cb2p,cb2d,cb1,cb0,ca0i,xp,xma,xmb,
 
90
     +          dmap,dmbp,dmamb,piDpj,ier)
 
91
***#[*comment:***********************************************************
 
92
*                                                                       *
 
93
*       see ffcb2, plus differences.                                    *
 
94
*                                                                       *
 
95
***#]*comment:***********************************************************
 
96
*  #[ declarations:
 
97
        implicit none
 
98
*
 
99
*       arguments
 
100
*
 
101
        integer ier
 
102
        DOUBLE COMPLEX xp,xma,xmb,dmap,dmbp,dmamb,piDpj(3,3)
 
103
        DOUBLE COMPLEX cb2p,cb2d,cb1,cb0,ca0i(2)
 
104
*
 
105
*       local variables
 
106
*
 
107
        integer i,j,ier0,init,ithres
 
108
        logical llogmm,lreal
 
109
        DOUBLE PRECISION absc,xmax,xmxp,rloss
 
110
        DOUBLE PRECISION rm1,rm2,rp,rm1p,rm2p,rm1m2,rpiDpj(3,3),sprec
 
111
        DOUBLE COMPLEX delsp,xlam,xlo3,xlogmm,zfflo1,zfflo3
 
112
        DOUBLE COMPLEX cc,cs(6),cb21,cb22,csom
 
113
        DOUBLE COMPLEX cqi(3),cqiqj(3,3),qiDqj(3,3)
 
114
        save init
 
115
*
 
116
*       common blocks
 
117
*
 
118
        include 'ff.h'
 
119
*
 
120
*       statement function
 
121
*
 
122
        absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
 
123
*
 
124
*  #] declarations:
 
125
*  #[ real case:
 
126
        if ( DIMAG(xma).eq.0 .and. DIMAG(xmb).eq.0 ) then
 
127
            lreal = .TRUE.
 
128
            if ( lwrite ) print *,'ffcb2a: real masses'
 
129
        elseif ( nschem.le.4 ) then
 
130
            lreal = .TRUE.
 
131
            if ( lwrite .or. init.eq.0 ) then
 
132
                init = 1
 
133
                print *,'ffcb2a: nschem <= 4, ignoring complex masses:',
 
134
     +                  nschem
 
135
            endif
 
136
        elseif ( nschem.le.6 ) then
 
137
            if ( lwrite .or. init.eq.0 ) then
 
138
                init = 1
 
139
                print *,'ffcb2a: nschem = 5,6 complex masses near ',
 
140
     +                  'threshold: ',nschem
 
141
            endif
 
142
            cqi(1) = xma
 
143
            cqi(2) = xmb
 
144
            cqi(3) = xp
 
145
            cqiqj(1,2) = dmamb
 
146
            cqiqj(2,1) = -cqiqj(1,2)
 
147
            cqiqj(1,3) = dmap
 
148
            cqiqj(3,1) = -cqiqj(1,3)
 
149
            cqiqj(2,3) = dmbp
 
150
            cqiqj(3,2) = -cqiqj(2,3)
 
151
            cqiqj(1,1) = 0
 
152
            cqiqj(2,2) = 0
 
153
            cqiqj(3,3) = 0
 
154
            call ffthre(ithres,cqi,cqiqj,3,1,2,3)
 
155
            if ( ithres.eq.0 .or. ithres.eq.1 .and. nschem.eq.5 ) then
 
156
                lreal = .TRUE.
 
157
                if ( lwrite ) print *,'ffcb2a: no threshold'
 
158
            else
 
159
                if ( lwrite ) print *,'ffcb2a: found threshold'
 
160
                lreal = .FALSE.
 
161
            endif
 
162
        else
 
163
            lreal = .FALSE.
 
164
        endif
 
165
        if ( lreal ) then
 
166
            rm1 = DBLE(xma)
 
167
            rm2 = DBLE(xmb)
 
168
            rp  = DBLE(xp)
 
169
            rm1p  = DBLE(dmap)
 
170
            rm2p  = DBLE(dmbp)
 
171
            rm1m2 = DBLE(dmamb)
 
172
            do 20 j=1,3
 
173
                do 10 i=1,3
 
174
                    rpiDpj(i,j) = DBLE(piDpj(i,j))
 
175
   10           continue
 
176
   20       continue
 
177
            sprec = precx
 
178
            precx = precc
 
179
            call ffxb2a(cb2p,cb2d,cb1,cb0,ca0i,rp,rm1,rm2,rm1p,rm2p,
 
180
     +          rm1m2,rpiDpj,ier)
 
181
            precx = sprec
 
182
            return
 
183
        endif
 
184
*  #] real case:
 
185
*  #[ test input:
 
186
        if ( ltest ) then
 
187
            ier0 = ier
 
188
            call ffcot2(qiDqj,xp,xma,xmb,dmap,dmbp,dmamb,ier0)
 
189
            rloss = xloss*DBLE(10)**(-mod(ier0,50))
 
190
            do 40 j=1,3
 
191
                do 30 i=1,3
 
192
                    if ( rloss*absc(piDpj(i,j)-qiDqj(i,j)).gt.precc*
 
193
     +                  absc(qiDqj(i,j)) ) print *,'ffcb2a: ',
 
194
     +                  'error: piDpj(',i,j,') wrong: ',piDpj(i,j),
 
195
     +                  qiDqj(i,j),piDpj(i,j)-qiDqj(i,j),ier0
 
196
   30           continue
 
197
   40       continue
 
198
        endif
 
199
*  #] test input:
 
200
*  #[ p^2 != 0:
 
201
        if ( xp .ne. 0 ) then
 
202
*       #[ normal case:
 
203
            call ffclmb(xlam,-xp,-xmb,-xma,dmbp,dmap,dmamb,ier)
 
204
            delsp = -xlam/4
 
205
*
 
206
*           the first one is simple...
 
207
*
 
208
            cs(1) = 2*cb1*piDpj(1,3)
 
209
            cs(2) = ca0i(2)
 
210
            cb2p = cs(1) + cs(2)
 
211
            if ( absc(cb2p) .lt. xloss*absc(cs(2)) ) then
 
212
                if ( lwarn ) call ffwarn(214,ier,absc(cb2p),absc(cs(2)))
 
213
            endif
 
214
*
 
215
*           the next one ain't.
 
216
*
 
217
            cs(1) = ca0i(2)
 
218
            cs(2) = 2*xma*cb0
 
219
            cs(3) = -2*piDpj(1,3)*cb1
 
220
            cs(4) = xma+xmb
 
221
            cs(5) = -xp/3
 
222
            cb2d = cs(1) + cs(2) + cs(3) + cs(4) + cs(5)
 
223
            xmax = max(absc(cs(2)),absc(cs(3)),absc(cs(4)),absc(cs(5)))
 
224
            if ( absc(cb2d) .ge. xloss*xmax ) goto 110
 
225
            if ( lwrite ) then
 
226
                print '(a,2e30.16,e12.4)','cb2d  = ',cb2d,xmax
 
227
                print *,'with cs '
 
228
                print '(i3,2e30.16)',(i,cs(i),i=1,5)
 
229
            endif
 
230
            if ( lwarn ) then
 
231
                call ffwarn(214,ier,absc(cb2d),xmax)
 
232
                if ( lwrite ) then
 
233
                    print *,'xp,xma,xmb = ',xp,xma,xmb
 
234
                endif
 
235
            endif
 
236
  110       continue
 
237
*       #] give up:
 
238
            cb2p = cb2p*DBLE(1/(2*xp))
 
239
            cb2d = cb2d*(1/DBLE(6))
 
240
*  #] p^2 != 0:
 
241
*  #[ p^2=0:
 
242
        elseif ( dmamb .ne. 0 ) then
 
243
            if ( init.eq.0 ) then
 
244
                init = 1
 
245
                print *,' '
 
246
                print *,'ffcb2a: note: in this case p^2=0 B21 is ',
 
247
     +                  'returned rather than B2p which is undefined'
 
248
                print *,' '
 
249
            endif
 
250
            if ( dmamb .ne. 0 ) then
 
251
*               #[ B21:
 
252
                llogmm = .FALSE.
 
253
*
 
254
*               B21 (see thesis, b21.frm)
 
255
*
 
256
                cs(1) = xma**2/3/dmamb**3*ca0i(1)
 
257
                cs(2) = (-xma**2 + xma*xmb - xmb**2/3)/dmamb**3*ca0i(2)
 
258
                cs(3) = (5*xma**3/18 - xma*xmb**2/2 + 2*xmb**3/9)
 
259
     +                  /dmamb**3
 
260
                cb21 = cs(1)+cs(2)+cs(3)
 
261
                xmax = max(absc(cs(2)),absc(cs(3)))
 
262
                if ( absc(cb21).gt.xloss**2*xmax ) goto 160
 
263
                if ( lwrite ) then
 
264
                    print *,'cb21 = ',cb21,xmax
 
265
                    print *,'with cs '
 
266
                    print '(i3,2e30.16)',(i,cs(i),i=1,3)
 
267
                endif
 
268
*
 
269
*               ma ~ mb
 
270
*
 
271
                if ( absc(dmamb).lt.xloss*absc(xma) ) then
 
272
                    xlogmm = zfflo1(dmamb/xma,ier)
 
273
                else
 
274
                    xlogmm = log(xmb/xma)
 
275
                endif
 
276
                llogmm = .TRUE.
 
277
                cs(1) = (xma/dmamb)/6
 
278
                cs(2) = (xma/dmamb)**2/3
 
279
                cs(3) = (xma/dmamb)**3*xlogmm/3
 
280
                cs(4) = -2/DBLE(9) + ca0i(1)/(3*xma)
 
281
                cs(5) = -xlogmm/3
 
282
                csom = cs(1)+cs(2)+cs(3)+cs(4)+cs(5)
 
283
                xmxp = max(absc(cs(2)),absc(cs(3)),absc(cs(4)),
 
284
     +                  absc(cs(5)))
 
285
                if ( lwrite ) then
 
286
                    print *,'cb21+= ',csom,xmxp
 
287
                    print *,'with cs '
 
288
                    print '(i3,2e30.16)',(i,cs(i),i=1,5)
 
289
                endif
 
290
                if ( xmxp.lt.xmax ) then
 
291
                    xmax = xmxp
 
292
                    cb21 = csom
 
293
                    if ( absc(cb21).gt.xloss**2*xmax ) goto 160
 
294
                endif
 
295
*
 
296
*               and last try
 
297
*
 
298
                xlo3 = zfflo3(dmamb/xma,ier)
 
299
                cs(1) = (dmamb/xma)**2/6
 
300
                cs(2) = (dmamb/xma)/3
 
301
                cs(3) = xlo3/(3*(dmamb/xma)**3)
 
302
*same           cs(4) = -2/DBLE(9) + ca0i(1)/(3*xma)
 
303
                cs(5) = -xlo3/3
 
304
                csom = cs(1)+cs(2)+cs(3)+cs(4)+cs(5)
 
305
                xmxp = max(absc(cs(2)),absc(cs(3)),absc(cs(4)),
 
306
     +                  absc(cs(5)))
 
307
                if ( lwrite ) then
 
308
                    print *,'cb21+= ',csom,xmxp
 
309
                    print *,'with cs '
 
310
                    print '(i3,2e30.16)',(i,cs(i),i=1,5)
 
311
                endif
 
312
                if ( xmxp.lt.xmax ) then
 
313
                    xmax = xmxp
 
314
                    cb21 = csom
 
315
                    if ( absc(cb21).gt.xloss**2*xmax ) goto 160
 
316
                endif
 
317
*
 
318
*               give up
 
319
*
 
320
                if ( lwarn ) then
 
321
                    call ffwarn(227,ier,absc(cb21),xmax)
 
322
                    if ( lwrite ) then
 
323
                        print *,'xp,xma,xmb = ',xp,xma,xmb
 
324
                    endif
 
325
                endif
 
326
  160           continue
 
327
*               #] B21:
 
328
*               #[ B22:
 
329
*
 
330
*               B22
 
331
*
 
332
                cs(1) = +xma/(4*dmamb)*ca0i(1)
 
333
                cs(2) = -xmb/(4*dmamb)*ca0i(2)
 
334
                cs(3) = (xma+xmb)/8
 
335
                cb22 = cs(1) + cs(2) + cs(3)
 
336
                xmax = max(absc(cs(2)),absc(cs(3)))
 
337
                if ( absc(cb22).gt.xloss*xmax ) goto 210
 
338
                if ( lwrite ) then
 
339
                    print *,'cb22 = ',cb22,xmax
 
340
                    print *,'with cs '
 
341
                    print '(i3,2e30.16)',(i,cs(i),i=1,3)
 
342
                endif
 
343
*
 
344
*               second try, close together
 
345
*
 
346
                if ( .not.llogmm ) then
 
347
                    if ( abs(dmamb).lt.xloss*absc(xma) ) then
 
348
                        xlogmm = zfflo1(dmamb/xma,ier)
 
349
                    else
 
350
                        xlogmm = log(xmb/xma)
 
351
                    endif
 
352
                endif
 
353
                cs(1) = dmamb*( -1/DBLE(8) - ca0i(1)/(4*xma) )
 
354
                cs(2) = dmamb*xlogmm/4
 
355
                cs(3) = xma*(xma/dmamb)/4*xlogmm
 
356
                cs(4) = xma*( 1/DBLE(4) + ca0i(1)/(2*xma) )
 
357
                cs(5) = -xma*xlogmm/2
 
358
                csom = cs(1) + cs(2) + cs(3) + cs(4) + cs(5)
 
359
                xmxp = max(absc(cs(2)),absc(cs(3)),absc(cs(4)),
 
360
     +                  absc(cs(5))) 
 
361
                if ( lwrite ) then
 
362
                    print *,'cb22+= ',csom,xmxp
 
363
                    print *,'with cs '
 
364
                    print '(i3,2e30.16)',(i,cs(i),i=1,2)
 
365
                endif
 
366
                if ( xmxp.lt.xmax ) then
 
367
                    xmax = xmxp
 
368
                    cb22 = csom
 
369
                endif
 
370
                if ( absc(cb22).gt.xloss*xmax ) goto 210
 
371
*
 
372
*               give up
 
373
*
 
374
                if ( lwarn ) then
 
375
                    call ffwarn(214,ier,absc(cb22),xmax)
 
376
                    if ( lwrite ) then
 
377
                        print *,'xp,xma,xmb = ',xp,xma,xmb
 
378
                    endif
 
379
                endif
 
380
  210           continue
 
381
*               #] B22:
 
382
            else
 
383
*
 
384
*               ma=mb: simple
 
385
*
 
386
                cb21 = cb0/3
 
387
                cb22 = xma/2*(cb0 + 1)
 
388
            endif
 
389
            cb2d = cb22
 
390
            cb2p = cb21
 
391
        endif
 
392
*  #] p^2=0:
 
393
*  #[ debug output:
 
394
        if ( lwrite ) then
 
395
            print *,'ffcb2: cb2p = ',cb2p,ier
 
396
            print *,'       cb2d = ',cb2d,ier
 
397
        endif
 
398
*  #] debug output:
 
399
*###] ffcb2a:
 
400
        end