1
* $Id: ffxe0.f,v 1.4 1996/01/10 15:36:51 gj Exp $
3
subroutine ffxe0(ce0,cd0i,xpi,ier)
4
***#[*comment:***********************************************************
9
* e0= -----\dq |(q^2-m_1^2)((q+p_1)^2-m_2^2)...((q-p_5)^2-m_5^2| *
12
* following the five four-point-function method in .... *
13
* As an extra the five fourpoint function Di are also returned *
14
* if ( ldot ) the dotproducts are left behind in fpij5(15,15) in *
15
* /ffdot/ and the external determinants fdel4 and fdl3i(5) in *
18
* Input: xpi = m_i^2 (real) i=1,5 *
19
* xpi = p_i.p_i (real) i=6,10 (note: B&D metric) *
20
* xpi = (p_i+p_{i+1})^2 (r) i=11,15 *
21
* xpi = (p_i+p_{i+2})^2 (r) i=16,20 OR 0 *
23
* Output: ce0 (complex) *
24
* cd0i(5) (complex) D0 with s_i missing *
25
* ier (integr) 0=ok 1=inaccurate 2=error *
27
***#]*comment:***********************************************************
33
DOUBLE PRECISION xpi(20)
34
DOUBLE COMPLEX ce0,cd0i(5)
39
integer i,j,NMIN,NMAX,ier0,i6,i7,i8,i9
40
parameter(NMIN=15,NMAX=20)
41
DOUBLE PRECISION dpipj(NMIN,NMAX),xmax
42
logical lp5(NMAX-NMIN)
50
* simulate the differences in the masses etc..
52
if ( lwrite ) print *,'ffxe0: input xpi: ',xpi
57
if ( xpi(i+15) .eq. 0 ) then
60
if ( i7 .ge. 11 ) i7 = 6
62
if ( i8 .ge. 11 ) i8 = 6
64
if ( i9 .ge. 11 ) i9 = 6
65
xpi(i+15) = xpi(i6)+xpi(i7)+xpi(i8)-xpi(i6+5)-xpi(i7+5)+
67
xmax = max(abs(xpi(i6)),abs(xpi(i7)),abs(xpi(i8)),abs(
68
+ xpi(i6+5)),abs(xpi(i7+5)),abs(xpi(i9+5)))
69
if ( abs(xpi(i+15)) .lt. xloss*xmax )
70
+ call ffwarn(168,ier,xpi(i+15),xmax)
77
* next the differences
82
if ( i .le. NMIN ) dpipj(i,i) = 0
83
do 10 j=1,min(i-1,NMIN)
84
dpipj(j,i) = xpi(j) - xpi(i)
85
if ( i .le. NMIN ) then
86
dpipj(i,j) = -dpipj(j,i)
88
* we do not need the differences of the u-like variables accurately
89
if ( i.gt.10 .and. j.gt.10 ) goto 10
90
if ( abs(dpipj(j,i)) .lt. xloss*abs(xpi(i))
91
+ .and. xpi(i) .ne. xpi(j) ) then
92
call ffwarn(158,ier0,dpipj(j,i),xpi(i))
93
if ( lwrite ) print *,'between xpi(',i,
101
dpipj(j,i) = xpi(j) - xpi(i)
105
* #] get differences:
107
call ffxe0a(ce0,cd0i,xpi,dpipj,ier)
119
subroutine ffxe0a(ce0,cd0i,xpi,dpipj,ier)
120
***#[*comment:***********************************************************
125
* e0= -----\dq |(q^2-m_1^2)((q+p_1)^2-m_2^2)...((q-p_5)^2-m_5^2| *
128
* following the five four-point-function method in .... *
129
* As an extra the five fourpoint function Di are also returned *
130
* if ( ldot ) the dotproducts are left behind in fpij5(15,15) in *
131
* /ffdot/ and the external determinants fdel4 and fdl3i(5) in *
134
* Input: xpi = m_i^2 (real) i=1,5 *
135
* xpi = p_i.p_i (real) i=6,10 (note: B&D metric) *
136
* xpi = (p_i+p_{i+1})^2 (r) i=11,15 *
137
* xpi = (p_i+p_{i+2})^2 (r) i=16,20 *
138
* dpipj(15,20) (real) = pi(i) - pi(j) *
140
* Output: ce0 (complex) *
141
* cd0i(5) (complex) D0 with s_i missing *
142
* ier (integer) <50:lost # digits 100=error *
144
***#]*comment:***********************************************************
151
DOUBLE COMPLEX ce0,cd0i(5)
152
DOUBLE PRECISION xpi(20),dpipj(15,20)
156
integer i,j,ii(10),ii4(6),ieri(5),ier0,imin,itype,ndiv,idone,
159
DOUBLE COMPLEX c,cfac,cs,csum
160
DOUBLE PRECISION dl5s,dl4p,xpi4(13),dpipj4(10,13),piDpj4(10,10),
161
+ absc,xmax,piDpj(15,15),xqi4(13),dqiqj4(10,13),
162
+ qiDqj4(10,10),del2s,xmx5(5),dl4ri(5)
171
absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
175
data ii4 /5,6,7,8,9,10/
178
* #[ initialisations:
185
* #] initialisations:
189
call ffdot5(piDpj,xpi,dpipj,ier)
193
fpij5(j,i) = piDpj(j,i)
201
call ffdl4p(dl4p,xpi,piDpj,15,ii,ier0)
202
* if ( dl4p .lt. 0 ) then
208
call ffdel5(dl5s,xpi,piDpj,15,ier)
210
print *,'ffxe0: dl5s = ',dl5s
218
if ( lwrite ) print *,'ffxe0a: fourpoint function nr ',i
220
* get the coefficient determinant
223
call ffdl4r(dl4ri(i),xpi,piDpj,15,i,ieri(i))
225
* get four-point momenta
227
call ffpi54(xpi4,dpipj4,piDpj4,xpi,dpipj,piDpj,i,ieri(i))
229
* first try IR divergent function to avoid error messages from ffrot4
232
call ffxdir(cs,cfac,idone,xpi4,dpipj4,6,ndiv,ier1)
233
if ( idone .gt. 0 ) then
235
xmax = abs(cs)*10d0**(-mod((ier1-ieri(i)),50))
238
* rotate to calculable posistion
240
call ffrot4(irota4,del2s,xqi4,dqiqj4,qiDqj4,xpi4,dpipj4,
241
+ piDpj4,5,itype,ieri(i))
242
if ( itype .lt. 0 ) then
243
print *,'ffxe0: error: Cannot handle this ',
244
+ ' 4point masscombination yet:'
245
print *,(xpi(j),j=1,20)
248
if ( itype .eq. 1 ) then
251
print *,'ffxe0a: Cannot handle del2s = 0 yet'
256
if ( itype .eq. 2 ) then
257
print *,'ffxe0a: no doubly IR divergent yet'
261
* get fourpoint function
266
call ffxd0e(cs,cfac,xmax, .FALSE.,ndiv,xqi4,dqiqj4,
267
+ qiDqj4,del2s,ldel2s,ieri(i))
268
if ( ieri(i).gt.10 ) then
270
print *,'ffxe0: id = ',id,', nevent = ',nevent
271
print *,'ffxe0: lost ',ieri(i),
272
+ ' digits in D0 with isgnal ',isgnal,
273
+ ', trying other roots, isgnal ',-isgnal
277
call ffxd0e(cs,cfac,xmax, .TRUE.,ndiv,xqi4,dqiqj4,
278
+ qiDqj4,del2s,ldel2s,ieri(i))
287
xmx5(i) = xmax*absc(cfac)
289
call ffdl3p(fdl3i(i),piDpj4,10,ii4,ii4,ieri(i))
290
* let's hope tha tthese have been set by ffxd0e...
294
call ffdel4(fdel4s,xpi4,piDpj4,10,ier0)
295
if ( xloss*10d0**(-ier0-1)*abs(fdl4si(i)-fdel4s)
296
+ .gt. precx*abs(fdel4s) ) then
297
print *,'ffxe0a: error: Del4s was not correct',
298
+ fdl4si(i),fdel4s,fdl4si(i)-fdel4s,ier0
301
if ( lwrite ) print *,'ffxe0: fdel4s = ',fdel4s
313
csum = csum + imin*DBLE(dl4ri(i))*cd0i(i)
314
if ( ieri(i) .gt. 50 ) then
315
ieri(i) = mod(ieri(i),50)
317
xmax = max(xmax,dl4ri(i)*xmx5(i)*DBLE(10)**mod(ieri(i),50))
320
* Check for cancellations in the final adding up
322
if ( lwarn .and. 2*absc(csum) .lt. xloss*xmax )
323
+ call ffwarn(161,ier,absc(csum),xmax)
325
* Check for a sum close to the minimum of the range (underflow
328
if ( lwarn .and. absc(csum).lt.xalogm/precc .and. csum.ne.0 )
329
+ call ffwarn(162,ier,absc(csum),xalogm/precc)
331
* If the imaginary part is very small it most likely is zero
332
* (can be removed, just esthetically more pleasing)
334
if ( abs(DIMAG(csum)) .lt. precc*abs(DBLE(csum)) )
335
+ csum = DCMPLX(DBLE(csum))
339
ce0 = csum*(1/DBLE(2*dl5s))
343
print '(a,5e16.8,i6)','cs,del4r,D0 = ',
344
+ DBLE(dl4ri(i))*cd0i(i)*(1/DBLE(2*dl5s)),
345
+ dl4ri(i)/DBLE(2*dl5s),cd0i(i),ieri(i)
347
print '(a,2e24.16,i6)','ffxe0a: ce0 = ',ce0,ier
353
subroutine ffxe00(ce0,cd0i,dl4ri,xpi,piDpj,ier)
354
***#[*comment:***********************************************************
359
* e0= -----\dq |(q^2-m_1^2)((q+p_1)^2-m_2^2)...((q-p_5)^2-m_5^2| *
362
* following the five four-point-function method in .... *
363
* The four five fourpoint function Di are input in this version. *
365
* Input: cd0i(5) (complex) D0 with s_i missing *
366
* dl4ri(5) (real) coeff of D0 *
367
* xpi = m_i^2 (real) i=1,5 *
368
* xpi = p_i.p_i (real) i=6,10 (note: B&D metric) *
369
* xpi = (p_i+p_{i+1})^2 (r) i=11,15 *
370
* xpi = (p_i+p_{i+2})^2 (r) i=16,20 *
371
* piDpj(15,15) (real) pi.pj *
373
* Output: ce0 (complex) *
374
* ier (integer) <50:lost # digits 100=error *
376
***#]*comment:***********************************************************
383
DOUBLE COMPLEX ce0,cd0i(5)
384
DOUBLE PRECISION dl4ri(5),xpi(20),piDpj(15,15)
388
integer i,ii(10),imin,ier0
389
DOUBLE COMPLEX c,csum
390
DOUBLE PRECISION dl5s,dl4p,absc,xmax
398
absc(c) = abs(DBLE(c)) + abs(DIMAG(c))
400
* #[ initialisations:
405
print *,'ffxe00: input:'
406
print *,' cd0i = ',cd0i
407
print *,' dl4ri = ',dl4ri
408
print *,' xpi = ',xpi
411
* #] initialisations:
420
call ffdl4p(dl4p,xpi,piDpj,15,ii,ier0)
424
call ffdel5(dl5s,xpi,piDpj,15,ier)
426
print *,'ffxe00: dl5s = ',dl5s
437
csum = csum + imin*DBLE(dl4ri(i))*cd0i(i)
438
xmax = max(xmax,abs(dl4ri(i))*absc(cd0i(i)))
441
* Check for cancellations in the final adding up
443
if ( lwarn .and. 2*absc(csum) .lt. xloss*xmax )
444
+ call ffwarn(161,ier,absc(csum),xmax)
446
* Check for a sum close to the minimum of the range (underflow
449
if ( lwarn .and. absc(csum).lt.xalogm/precc .and. csum.ne.0 )
450
+ call ffwarn(162,ier,absc(csum),xalogm/precc)
452
* If the imaginary part is very small it most likely is zero
453
* (can be removed, just esthetically more pleasing)
455
if ( abs(DIMAG(csum)) .lt. precc*abs(DBLE(csum)) )
456
+ csum = DCMPLX(DBLE(csum))
460
ce0 = csum*(1/DBLE(2*dl5s))
466
subroutine ffdot5(piDpj,xpi,dpipj,ier)
467
***#[*comment:***********************************************************
469
* calculate the dotproducts pi.pj with *
471
* xpi(i) = s_i i=1,5 *
472
* xpi(i) = p_i i=6,10 *
473
* xpi(i) = p_i+p_{i+1} i=11,15 *
475
***#]*comment:***********************************************************
482
DOUBLE PRECISION xpi(20),dpipj(15,20),piDpj(15,15)
486
integer is1,is2,is3,is4,ip6,ip7,ip8,ip11,ip12,ip14,i,j,
487
+ igehad(15,15),itel,i1,i2,i3,i4,i5,i6,ierin,ier0
488
* werkt niet bij Absoft
489
* parameter (locwrt=.TRUE.)
491
DOUBLE PRECISION xheck,xmax
499
data locwrt /.FALSE./
502
if ( ltest ) call ffxhck(xpi,dpipj,15,ier)
515
if ( is2 .eq. 6 ) is2 = 1
517
if ( is3 .eq. 6 ) is3 = 1
522
* we have now defined a 3point function
532
* #[ all in one vertex:
536
piDpj(is1,is1) = xpi(is1)
537
piDpj(ip6,ip6) = xpi(ip6)
538
piDpj(ip11,ip11) = xpi(ip11)
540
igehad(is1,is1) = igehad(is1,is1) + 1
541
igehad(ip6,ip6) = igehad(ip6,ip6) + 1
542
igehad(ip11,ip11) = igehad(ip11,ip11) + 1
547
if ( xpi(is2) .le. xpi(is1) ) then
548
piDpj(is1,is2) = (dpipj(is1,ip6) + xpi(is2))/2
550
piDpj(is1,is2) = (dpipj(is2,ip6) + xpi(is1))/2
552
piDpj(is2,is1) = piDpj(is1,is2)
554
igehad(is1,is2) = igehad(is1,is2) + 1
555
igehad(is2,is1) = igehad(is2,is1) + 1
557
if ( lwarn .and. abs(piDpj(is1,is2)) .lt.
558
+ xloss*min(xpi(is1),xpi(is2)) ) then
560
call ffwarn(105,ier0,piDpj(is1,is2),min(xpi(is1),
567
if ( xpi(is1) .le. xpi(is3) ) then
568
piDpj(is3,is1) = (dpipj(is3,ip11) + xpi(is1))/2
570
piDpj(is3,is1) = (dpipj(is1,ip11) + xpi(is3))/2
572
piDpj(is1,is3) = piDpj(is3,is1)
574
igehad(is1,is3) = igehad(is1,is3) + 1
575
igehad(is3,is1) = igehad(is3,is1) + 1
577
if ( lwarn .and. abs(piDpj(is1,is3)) .lt.
578
+ xloss*min(xpi(is1),xpi(is3)) ) then
580
call ffwarn(106,ier0,
581
+ piDpj(is1,is3),min(xpi(is1),xpi(is3)))
587
if ( abs(xpi(ip6)) .le. xpi(is1) ) then
588
piDpj(ip6,is1) = (dpipj(is2,is1) - xpi(ip6))/2
590
piDpj(ip6,is1) = (dpipj(is2,ip6) - xpi(is1))/2
592
piDpj(is1,ip6) = piDpj(ip6,is1)
594
igehad(is1,ip6) = igehad(is1,ip6) + 1
595
igehad(ip6,is1) = igehad(ip6,is1) + 1
597
if ( lwarn .and. abs(piDpj(ip6,is1)) .lt.
598
+ xloss*min(abs(xpi(ip6)),xpi(is1))) then
600
call ffwarn(107,ier0,
601
+ piDpj(ip6,is1),min( abs(xpi(ip6)),xpi(is1)))
607
if ( abs(xpi(ip6)) .le. xpi(is2) ) then
608
piDpj(ip6,is2) = (dpipj(is2,is1) + xpi(ip6))/2
610
piDpj(ip6,is2) = (dpipj(ip6,is1) + xpi(is2))/2
613
igehad(is2,ip6) = igehad(is2,ip6) + 1
614
igehad(ip6,is2) = igehad(ip6,is2) + 1
616
piDpj(is2,ip6) = piDpj(ip6,is2)
617
if ( lwarn .and. abs(piDpj(ip6,is2)) .lt.
618
+ xloss*min(abs(xpi(ip6)),xpi(is2))) then
620
call ffwarn(108,ier0,
621
+ piDpj(ip6,is2),min(abs(xpi(ip6)),xpi (is2)))
627
if ( abs(xpi(ip11)) .le. xpi(is1) ) then
628
piDpj(ip11,is1) = -(dpipj(is1,is3) + xpi(ip11))/2
630
piDpj(ip11,is1) = -(dpipj(ip11,is3) + xpi(is1))/2
632
piDpj(is1,ip11) = piDpj(ip11,is1)
634
igehad(is1,ip11) = igehad(is1,ip11) + 1
635
igehad(ip11,is1) = igehad(ip11,is1) + 1
637
if ( lwarn .and. abs(piDpj(ip11,is1)) .lt.
638
+ xloss*min(abs(xpi(ip11)),xpi(is1))) then
640
call ffwarn(109,ier0,
641
+ piDpj(ip11,is1),min(abs(xpi(ip11)),xpi(is1)))
647
if ( abs(xpi(ip11)) .le. xpi(is3) ) then
648
piDpj(ip11,is3) = -(dpipj(is1,is3) - xpi(ip11))/2
650
piDpj(ip11,is3) = -(dpipj(is1,ip11) - xpi(is3))/2
652
piDpj(is3,ip11) = piDpj(ip11,is3)
654
igehad(is3,ip11) = igehad(is3,ip11) + 1
655
igehad(ip11,is3) = igehad(ip11,is3) + 1
657
if ( lwarn .and. abs(piDpj(ip11,is3)) .lt.
658
+ xloss*min(abs(xpi(ip11)),xpi(is3))) then
660
call ffwarn(109,ier0,
661
+ piDpj(ip11,is3),min(abs(xpi(ip11)),xpi(is3)))
664
* #] all in one vertex:
665
* #[ all in one 3point:
669
if ( min(abs(dpipj(is2,is1)),abs(dpipj(ip11,ip7))) .le.
670
+ min(abs(dpipj(ip11,is1)),abs(dpipj(is2,ip7))) ) then
671
piDpj(ip6,is3) = (dpipj(ip11,ip7) + dpipj(is2,is1))/2
673
piDpj(ip6,is3) = (dpipj(ip11,is1) + dpipj(is2,ip7))/2
675
piDpj(is3,ip6) = piDpj(ip6,is3)
677
igehad(is3,ip6) = igehad(is3,ip6) + 1
678
igehad(ip6,is3) = igehad(ip6,is3) + 1
680
if ( lwarn .and. abs(piDpj(ip6,is3)) .lt.
681
+ xloss*min(abs(dpipj(ip11,ip7)),abs(dpipj(ip11,is1))) )
684
call ffwarn(110,ier0,piDpj(ip6,is3),
685
+ min(abs(dpipj(ip11,ip7)),abs(dpipj(ip11,is1))))
691
if ( min(abs(dpipj(is3,is2)),abs(dpipj(ip6,ip11))) .le.
692
+ min(abs(dpipj(ip6,is2)),abs(dpipj(is3,ip11))) ) then
693
piDpj(ip7,is1) = (dpipj(ip6,ip11) + dpipj(is3,is2))/2
695
piDpj(ip7,is1) = (dpipj(ip6,is2) + dpipj(is3,ip11))/2
697
piDpj(is1,ip7) = piDpj(ip7,is1)
699
igehad(is1,ip7) = igehad(is1,ip7) + 1
700
igehad(ip7,is1) = igehad(ip7,is1) + 1
702
if ( lwarn .and. abs(piDpj(ip7,is1)) .lt.
703
+ xloss*min(abs(dpipj(ip6,ip11)),abs(dpipj(ip6,is2))) )
706
call ffwarn(111,ier0,piDpj(ip7,is1),
707
+ min(abs(dpipj(ip6,ip11)),abs(dpipj(ip6,is2))))
713
if ( min(abs(dpipj(is1,is3)),abs(dpipj(ip7,ip6))) .le.
714
+ min(abs(dpipj(ip7,is3)),abs(dpipj(is1,ip6))) ) then
715
piDpj(ip11,is2) = -(dpipj(ip7,ip6) + dpipj(is1,is3))/2
717
piDpj(ip11,is2) = -(dpipj(ip7,is3) + dpipj(is1,ip6))/2
719
piDpj(is2,ip11) = piDpj(ip11,is2)
721
igehad(is2,ip11) = igehad(is2,ip11) + 1
722
igehad(ip11,is2) = igehad(ip11,is2) + 1
724
if ( lwarn .and. abs(piDpj(ip11,is2)) .lt.
725
+ xloss*min(abs(dpipj(ip7,ip6)),abs(dpipj(ip7,is3))) )
728
call ffwarn(112,ier0,piDpj(ip11,is2),
729
+ min(abs(dpipj(ip7,ip6)),abs(dpipj(ip7,is3))))
732
* #] all in one 3point:
733
* #[ all external 3point:
737
if ( abs(xpi(ip7)) .le. abs(xpi(ip6)) ) then
738
piDpj(ip6,ip7) = (dpipj(ip11,ip6) - xpi(ip7))/2
740
piDpj(ip6,ip7) = (dpipj(ip11,ip7) - xpi(ip6))/2
742
piDpj(ip7,ip6) = piDpj(ip6,ip7)
744
igehad(ip7,ip6) = igehad(ip7,ip6) + 1
745
igehad(ip6,ip7) = igehad(ip6,ip7) + 1
747
if ( lwarn .and. abs(piDpj(ip6,ip7)) .lt.
748
+ xloss*min(abs(xpi(ip6)),abs(xpi(ip7))) ) then
750
call ffwarn(113,ier0,piDpj(ip6,ip7),
751
+ min(abs(xpi(ip6)),abs(xpi(ip7))))
757
if ( abs(xpi(ip11)) .le. abs(xpi(ip7)) ) then
758
piDpj(ip7,ip11) = -(dpipj(ip6,ip7) - xpi(ip11))/2
760
piDpj(ip7,ip11) = -(dpipj(ip6,ip11) - xpi(ip7))/2
762
piDpj(ip11,ip7) = piDpj(ip7,ip11)
764
igehad(ip11,ip7) = igehad(ip11,ip7) + 1
765
igehad(ip7,ip11) = igehad(ip7,ip11) + 1
767
if ( lwarn .and. abs(piDpj(ip7,ip11)) .lt.
768
+ xloss*min(abs(xpi(ip7)),abs(xpi(ip11))) ) then
770
call ffwarn(114,ier0,piDpj(ip7,ip11),
771
+ min(abs(xpi(ip7)),abs(xpi(ip11))))
777
if ( abs(xpi(ip6)) .le. abs(xpi(ip11)) ) then
778
piDpj(ip11,ip6) = -(dpipj(ip7,ip11) - xpi(ip6))/2
780
piDpj(ip11,ip6) = -(dpipj(ip7,ip6) - xpi(ip11))/2
782
piDpj(ip6,ip11) = piDpj(ip11,ip6)
784
igehad(ip6,ip11) = igehad(ip6,ip11) + 1
785
igehad(ip11,ip6) = igehad(ip11,ip6) + 1
787
if ( lwarn .and. abs(piDpj(ip11,ip6)) .lt.
788
+ xloss*min(abs(xpi(ip11)),abs(xpi(ip6))) ) then
790
call ffwarn(115,ier0,piDpj(ip11,ip6),
791
+ min(abs(xpi(ip11)),abs(xpi(ip6))))
794
* #] all external 3point:
795
* #[ the other 3point:
797
if ( is4 .eq. 6 ) is4 = 1
801
* we now work with the threepoint configuration
813
if ( itel .eq. 1 ) then
820
elseif ( itel .eq. 2 ) then
836
* in one go: the opposite sides
838
if ( min(abs(dpipj(i3,i2)),abs(dpipj(i4,i6))) .le.
839
+ min(abs(dpipj(i4,i2)),abs(dpipj(i3,i6))) ) then
840
piDpj(i5,i1) = (dpipj(i3,i2) + dpipj(i4,i6))/2
842
piDpj(i5,i1) = (dpipj(i4,i2) + dpipj(i3,i6))/2
844
piDpj(i1,i5) = piDpj(i5,i1)
846
igehad(i1,i5) = igehad(i1,i5) + 1
847
igehad(i5,i1) = igehad(i5,i1) + 1
849
if ( lwarn .and. abs(piDpj(i5,i1)) .lt. xloss*
850
+ min(abs(dpipj(i4,i6)),abs(dpipj(i4,i2))) ) then
852
call ffwarn(111,ier0,piDpj(i5,i1),
853
+ min(abs(dpipj(i4,i6)),abs(dpipj(i4,i2))))
857
* and the remaining external ones
859
if ( abs(xpi(i5)) .le. abs(xpi(i4)) ) then
860
piDpj(i4,i5) = (dpipj(i6,i4) - xpi(i5))/2
862
piDpj(i4,i5) = (dpipj(i6,i5) - xpi(i4))/2
864
piDpj(i5,i4) = piDpj(i4,i5)
866
igehad(i5,i4) = igehad(i5,i4) + 1
867
igehad(i4,i5) = igehad(i4,i5) + 1
869
if ( lwarn .and. abs(piDpj(i4,i5)) .lt.
870
+ xloss*min(abs(xpi(i4)),abs(xpi(i5))) ) then
872
call ffwarn(113,ier0,piDpj(i4,i5),
873
+ min(abs(xpi(i4)),abs(xpi(i5))))
877
* #] the other 3point:
881
* we now have the fourpoint configuration
894
if ( itel .eq. 1 ) then
905
if ( min(abs(dpipj(i3,ip11)),abs(dpipj(i4,ip12))) .le.
906
+ min(abs(dpipj(i4,ip11)),abs(dpipj(i3,ip12))) ) then
907
piDpj(i1,i2) = (dpipj(i3,ip11) + dpipj(i4,ip12))/2
909
piDpj(i1,i2) = (dpipj(i4,ip11) + dpipj(i3,ip12))/2
911
piDpj(i2,i1) = piDpj(i1,i2)
913
igehad(i1,i2) = igehad(i1,i2) + 1
914
igehad(i2,i1) = igehad(i2,i1) + 1
916
if ( lwarn .and. abs(piDpj(i2,i1)) .lt. xloss*
917
+ min(abs(dpipj(i4,ip12)),abs(dpipj(i4,ip11))))
920
call ffwarn(111,ier0,piDpj(i2,i1),
921
+ min(abs(dpipj(i4,ip12)),abs(dpipj(i4,ip11))))
926
* we are only left with p11.p12 etc.
928
if ( min(abs(dpipj(ip14,ip8)),abs(dpipj(ip7,ip6))) .le.
929
+ min(abs(dpipj(ip7,ip8)),abs(dpipj(ip14,ip6))) ) then
930
piDpj(ip11,ip12) = (dpipj(ip7,ip6) + dpipj(ip14,ip8))/2
932
piDpj(ip11,ip12) = (dpipj(ip7,ip8) + dpipj(ip14,ip6))/2
934
piDpj(ip12,ip11) = piDpj(ip11,ip12)
936
igehad(ip12,ip11) = igehad(ip12,ip11) + 1
937
igehad(ip11,ip12) = igehad(ip11,ip12) + 1
939
if ( lwarn .and. abs(piDpj(ip11,ip12)) .lt. xloss*
940
+ min(abs(dpipj(ip7,ip6)),abs(dpipj(ip7,ip8))) ) then
942
call ffwarn(112,ier0,piDpj(ip11,ip12),
943
+ min(abs(dpipj(ip7,ip6)),abs(dpipj(ip7,ip8))))
950
print *,'We hebben gehad:'
951
print '(15i2)',igehad
956
* sum over all (incoming) momenta => 0
961
xheck = xheck + piDpj(j,i)
962
xmax = max(abs(piDpj(j,i)),xmax)
964
if ( xloss*abs(xheck) .gt. precx*xmax ) print *,
965
+ 'ffdot5: error: dotproducts with p(',i,
966
+ ') wrong: (som(.p(i))<>0) ',
967
+ (piDpj(i,j),j=6,10),xheck
969
* sum over all (incoming) momentum pairs => 0
974
xheck = xheck + piDpj(j,i)
975
xmax = max(abs(piDpj(j,i)),xmax)
977
if ( xloss*abs(xheck) .gt. precx*xmax ) print *,
978
+ 'ffdot5: error: dotproducts with p(',i,
979
+ ') wrong: (som(.(p(i)+p(i+1)))<>0) ',
980
+ (piDpj(i,j),j=11,15),xheck
985
if ( piDpj(i,j) .ne. piDpj(j,i) ) print *,
986
+ 'ffdot5: error: piDpj(',i,j,') <> piDpj',j,i,')'
991
if ( piDpj(i,i) .ne. xpi(i) ) print *,'ffdot5: error: ',
992
+ 'piDpj(',i,i,') <> xpi(',i,')'
997
* see if indeed pi+p(i+1) = p(i+5)
1001
if ( i1 .eq. 11 ) i1 = 6
1004
* check that si+p(i+5) = s(i+1)
1009
xheck = piDpj(j,i)+piDpj(i1,i)-piDpj(i2,i)
1010
xmax = max(abs(piDpj(j,i)),abs(piDpj(i2,i)),
1012
if ( xloss*abs(xheck) .gt. precx*xmax ) print *,
1013
+ 'ffdot5: error: piDpj(',j,i,')+piDpj(',
1014
+ i2,i,')-piDpj(',i1,i,') <> 0',xmax,xheck
1023
subroutine ffpi54(xpi4,dpipj4,piDpj4,xpi,dpipj,piDpj,inum,ier)
1024
***#[*comment:***********************************************************
1026
* Gets the dotproducts pertaining to the fourpoint function with *
1027
* s_i missing out of the five point function dotproduct array. *
1029
* Input: xpi real(20) si.si,pi.pi *
1030
* dpipj real(15,20) xpi(i) - xpi(j) *
1031
* piDpj real(15,15) pi(i).pi(j) *
1032
* inum integer 1--5 *
1034
* Output: xpi4 real(13) *
1035
* dpipj4 real(10,13) *
1036
* piDpj4 real(10,10) *
1039
***#]*comment:***********************************************************
1046
DOUBLE PRECISION xpi(20),dpipj(15,20),piDpj(15,15),xpi4(13),
1047
+ dpipj4(10,13),piDpj4(10,10),qDq(10,10)
1051
integer i,j,iplace(11,5),isigns(11,5),ier0
1053
DOUBLE PRECISION xmax
1062
+ 2,3,4,5, 07,08,09,15, 12,13, 17,
1063
+ 1,3,4,5, 11,08,09,10, 14,13, 18,
1064
+ 1,2,4,5, 06,12,09,10, 14,15, 19,
1065
+ 1,2,3,5, 06,07,13,10, 11,15, 20,
1066
+ 1,2,3,4, 06,07,08,14, 11,12, 16/
1069
+ +1,+1,+1,+1, +1,+1,+1,+1, -1,+1, +1,
1070
+ +1,+1,+1,+1, +1,+1,+1,+1, +1,+1, +1,
1071
+ +1,+1,+1,+1, +1,+1,+1,+1, +1,-1, +1,
1072
+ +1,+1,+1,+1, +1,+1,+1,+1, -1,-1, +1,
1073
+ +1,+1,+1,+1, +1,+1,+1,+1, -1,+1, +1/
1080
xpi4(i) = xpi(iplace(i,inum))
1082
dpipj4(j,i) = dpipj(iplace(j,inum),iplace(i,inum))
1086
* these cannot be simply copied I think
1088
xpi4(12) = -xpi4(5)+xpi4(6)-xpi4(7)+xpi4(8)+xpi4(9)+xpi4(10)
1090
xmax = max(abs(xpi4(5)),abs(xpi4(6)),abs(xpi4(7)),
1091
+ abs(xpi4(8)),abs(xpi4(9)),abs(xpi4(10)))
1092
if ( abs(xpi4(12)) .lt. xloss*xmax )
1093
+ call ffwarn(154,ier,xpi4(12),xmax)
1095
xpi4(13) = xpi4(5)-xpi4(6)+xpi4(7)-xpi4(8)+xpi4(9)+xpi4(10)
1097
xmax = max(abs(xpi4(5)),abs(xpi4(6)),abs(xpi4(7)),
1098
+ abs(xpi4(8)),abs(xpi4(9)),abs(xpi4(10)))
1099
if ( abs(xpi4(13)) .lt. xloss*xmax )
1100
+ call ffwarn(155,ier,xpi4(13),xmax)
1103
* and the differences
1107
dpipj4(j,i) = xpi4(j) - xpi4(i)
1111
* copy the dotproducts (watch the signs of p9,p10!)
1115
piDpj4(j,i) = isigns(j,inum)*isigns(i,inum)*
1116
+ piDpj(iplace(j,inum),iplace(i,inum))
1122
print *,'ffpi54: xpi4 = ',xpi4
1126
call ffxhck(xpi4,dpipj4,10,ier0)
1127
call ffxuvw(xpi4,dpipj4,ier0)
1128
if ( ier0 .ne. 0 ) print *,'ffpi54: error detected'
1133
call ffdot4(qDq,xpi4,dpipj4,10,ier0)
1136
if ( xloss*10d0**(-mod(ier0,50))*abs(qDq(j,i)-
1137
+ piDpj4(j,i)) .gt. precx*abs(qDq(j,i)) ) print *,
1138
+ 'ffpi54: error: piDpj4(',j,i,') not correct: ',
1139
+ piDpj4(j,i),qDq(j,i),piDpj4(j,i)-qDq(j,i),ier0
1147
subroutine ffxe0r(ce0,cd0i,xpi,ier)
1148
***#[*comment:***********************************************************
1150
* Tries all 12 permutations of the 5pointfunction *
1152
***#]*comment:***********************************************************
1157
DOUBLE PRECISION xpi(20),xqi(20)
1158
DOUBLE COMPLEX ce0,cd0i(5),ce0p,cd0ip(5),cd0ipp(5)
1159
integer inew(20,nrot),irota,ier1,i,j,k,icon,ialsav,init
1165
+ /1,2,3,4,5, 6,7,8,9,10,11,12,13,14,15, 16,17,18,19,20,
1166
+ 2,1,3,4,5, 6,11,8,9,15,7,14,13,12,10, 16,18,17,19,-20,
1167
+ 1,3,2,4,5, 11,7,12,9,10,6,8,15,14,13, -16,17,19,18,20,
1168
+ 1,2,4,3,5, 6,12,8,13,10,14,7,9,11,15, 16,-17,18,20,19,
1169
+ 1,2,3,5,4, 6,7,13,9,14,11,15,8,10,12, 20,17,-18,19,16,
1170
+ 5,2,3,4,1, 15,7,8,14,10,13,12,11,9,6, 17,16,18,-19,20,
1171
+ 2,1,4,3,5, 6,14,8,13,15,12,11,9,7,10, 16,-18,17,20,-19,
1172
+ 1,3,2,5,4, 11,7,15,9,14,6,13,12,10,8, -20,17,-19,18,16,
1173
+ 5,2,4,3,1, 15,12,8,11,10,9,7,14,13,6, 17,-16,18,-20,19,
1174
+ 2,1,3,5,4, 6,11,13,9,12,7,10,8,15,14, 20,18,-17,19,-16,
1175
+ 5,3,2,4,1, 13,7,12,14,10,15,8,6,9,11, -17,16,19,-18,20,
1176
+ 1,3,5,2,4, 11,13,15,12,14,10,7,9,6,8,-20,-17,-19,-16,-18/
1179
* #[ open console for some activity on screen:
1180
if ( init .eq. 0 ) then
1183
open(icon,file='CON:',status='old',err=11)
1191
* #] open console for some activity on screen:
1199
if ( inew(i,irota) .lt. 0 ) then
1200
xqi(-inew(i,irota)) = 0
1202
xqi(inew(i,irota)) = xpi(i)
1205
print '(a,i2,a,i2)','---#[ rotation ',irota,': isgnal ',
1207
if (lcon) write(icon,'(a,i2,a,i2)')'rotation ',irota,',
1213
call ffxe0(ce0p,cd0ip,xqi,ier1)
1215
print '(a,i1,a,i2)','---#] rotation ',irota,': isgnal ',
1217
print '(a,2g28.16,i3)','e0 = ',ce0p,ier1
1219
cd0ipp(k) = cd0ip(inew(k,irota))
1220
print '(a,2g28.16,i3)','d0 = ',cd0ipp(k),k
1222
if (lcon) write(icon,'(a,2g28.16,i3)')'e0 = ',ce0p,ier1
1223
if ( ier1 .lt. ier ) then
added
removed
vendor/CutTools/src/qcdloop/ffxe0.f
