← Back to branch summary

~maddevelopers/mg5amcnlo/2.9.4

« back to all changes in this revision

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

pass to v2.0.0

Show diffs side-by-side

added added

removed removed

Lines of Context:
vendor/CutTools/src/qcdloop/aaxdx.f
 
1
*###[ aaxdx :
 
2
        subroutine aaxdx(cbxi,ccxi,cdxi,d0,xmm,xpi,level,ier)
 
3
***#[*comment:***********************************************************
 
4
*                                                                       *
 
5
*       Calculation of four point tensor integrals.  Just a wrapper     *
 
6
*       for ffxdx nowadays, see there for the real description.         *
 
7
*                                                                       *
 
8
*       Input:                                                          *
 
9
*               xpi       the same as in Geert Jan's routines           *
 
10
*               level     rank of tensor(integral)                      *
 
11
*       Output:                                                         *
 
12
*               cbxi(12)  cb0(1),cb1(1),[cb2(2)]         x 6            *
 
13
*               ccxi(28)  cc0(1),cc1(2),cc2(4),[cc3(6)]  x 4            *
 
14
*               cdxi(24)  cd0(1),cd1(3),cd2(7),cd3(13)                  *
 
15
*                                                                       *
 
16
***#]*comment:***********************************************************
 
17
*  #[ declarations :
 
18
        implicit none
 
19
*
 
20
*       arguments
 
21
*
 
22
        integer ier,level
 
23
        DOUBLE PRECISION xpi(13),d0,xmm
 
24
        DOUBLE COMPLEX cbxi(12),ccxi(28),cdxi(24)
 
25
*
 
26
*       local variables
 
27
*
 
28
        DOUBLE COMPLEX caxi(4)
 
29
        DOUBLE PRECISION maxi(4),mbxi(12),mcxi(28),mdxi(24)
 
30
*
 
31
*  #] declarations :
 
32
*  #[ call ffxdx:
 
33
*
 
34
        call ffxdx(caxi,maxi,cbxi,mbxi,ccxi,mcxi,cdxi,mdxi,d0,xmm,xpi,
 
35
     +          level,ier)
 
36
*
 
37
*  #] call ffxdx:
 
38
*###] aaxdx :
 
39
        end
 
40
*###[ ffxdx:
 
41
        subroutine ffxdx(caxi,maxi,cbxi,mbxi,ccxi,mcxi,cdxi,mdxi,d0,xmm,
 
42
     +          xpi,level,ier)
 
43
***#[*comment:***********************************************************
 
44
*                                                                       *
 
45
*       Calculation of four point form factors with more accurate       *
 
46
*       error estimates.  Calls ffxd1, the rest is still here.          *
 
47
*       Definitions:                                                    *
 
48
*       D0 = 1/(i pi^2)*\int d^4 Q 1/((Q^2-m_1^2)((Q+p1)^2-m2^2)        *
 
49
*               ((Q+p1+p2)^2-m3^2))((Q-p4)^2-m4^2))                     *
 
50
*       D1 = 1/(i pi^2)*\int d^n Q Q(mu)/(...)                          *
 
51
*          = D11*p1 + D12*p2 + D13*p3                                   *
 
52
*       D2 = D21*p1*p1 + D22*p2*p2 + D23*p3*p3 + D24*(p1*p2+p2*p1) +    *
 
53
*            D25*(p1*p3+p3*p1) + D26*(p2*p3+p3*p2) + D27*g              *
 
54
*       D3 = D31*p1*p1*p1 + D32*p2*p2*p2 + D33*p3*p3*p3 + D34*(p1*p1*p2+*
 
55
*            p1*p2*p1+p2*p1*p1) + D35*(p1*p1*p3+p1*p3*p1+p3*p1*p1) + D36*
 
56
*            *(p1*p2*p2+p2*p1*p2+p1*p2*p2) + D37*(p1*p3*p3+p3*p1*p3+p1* *
 
57
*            p3*p3) + D38*(p2*p2*p3+p2*p3*p2+p3*p2*p2) + D39*(p2*p3*p3+ *
 
58
*            p3*p2*p3+p2*p3*p3) + D310*(p1*p2*p3+p2*p3*p1+p3*p1*p2+p1*p3*
 
59
*            *p2+p3*p2*p1+p2*p1*p3) + D311*(p1*g+g*p1+'g*p1*g') + D312* *
 
60
*            (p2*g+g*p2+'g*p2*g') + D313*(p3*g+g*p3+'g*p3*g')           *
 
61
*       D4 has not yet been implemented                                 *
 
62
*                                                                       *
 
63
*       Input:  xpi(13)      real       m_i^2 (1:4), p_{i-4}^2 (4:8),s,t*
 
64
*                                       optionally u,v,w (see ffxd0a)   *
 
65
*               d0,xmm       real       renormalisation constants       *
 
66
*               level        integer    rank of tensor (integral)       *
 
67
*       Output: caxi(4)      complex    A0(m_i^2) only when level>3     *
 
68
*               maxi(12)     real       max term in sum to caxi()       *
 
69
*               cbxi(12)     complex    6x(B0,B11,B21,B22)(p_i^2)       *
 
70
*               mbxi(12)     real       max term in sum to cbxi()       *
 
71
*               ccxi(28)     complex    4x(C0,C1(2),C2(4))              *
 
72
*               mcxi(28)     real       max term in sum to ccxi()       *
 
73
*               cdxi(24)     complex    D0,D1(3),D2(7),D3(13)           *
 
74
*               mdxi(24)     real       max term in sum to cdxi()       *
 
75
*       Note that if level<3 some of these are not defined.             *
 
76
*                                                                       *
 
77
***#]*comment:***********************************************************
 
78
*  #[ declarations :
 
79
        implicit none
 
80
*
 
81
*       arguments
 
82
*
 
83
        integer ier,level
 
84
        DOUBLE PRECISION xpi(13),d0,xmm
 
85
        DOUBLE COMPLEX caxi(4),cbxi(12),ccxi(28),cdxi(24)
 
86
        DOUBLE PRECISION maxi(4),mbxi(12),mcxi(28),mdxi(24)
 
87
*
 
88
*       local variables
 
89
*
 
90
        integer i,j,cl,ier0,ier1,iinx(6,4)
 
91
        DOUBLE PRECISION xpi3(6),fdel2i(4),absc,big
 
92
        DOUBLE COMPLEX caxj(12),cbxj(48),ccxj(52),cc,cd0
 
93
        DOUBLE PRECISION maxj(12),mbxj(48),mcxj(52)
 
94
        save iinx
 
95
*
 
96
*       common blocks
 
97
*
 
98
        include 'ff.h'
 
99
        include 'aa.h'
 
100
*
 
101
*       statement function
 
102
*
 
103
        absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
 
104
*
 
105
*       data
 
106
*
 
107
        data iinx /2,3,4,6,7,10,
 
108
     +             1,3,4,9,7,8,
 
109
     +             1,2,4,5,10,8,
 
110
     +             1,2,3,5,6,9/
 
111
*
 
112
*  #] declarations :
 
113
*  #[ initialisations:
 
114
*
 
115
*       initialize to something ridiculous so that one immediately
 
116
*       notices when it is accidentally used...
 
117
*
 
118
        big = 1/(1d20*xclogm)
 
119
        do 8 i=1,4
 
120
            caxi(i) = big
 
121
    8   continue
 
122
        do 9 i=1,12
 
123
            cbxi(i) = big
 
124
    9   continue
 
125
        do 10 i=1,28
 
126
            ccxi(i) = big
 
127
   10   continue
 
128
        do 11 i=1,24
 
129
            cdxi(i) = big
 
130
   11   continue
 
131
*
 
132
*  #] initialisations:
 
133
*  #[ get D0:
 
134
*       D0-function (ff)
 
135
*          futhermore dotpr and determinants are delivered by ff
 
136
        ldot = .TRUE.
 
137
        ier1 = ier
 
138
        call ffxd0(cdxi(1),xpi,ier1)
 
139
        if ( ier1.gt.10 ) then
 
140
            if ( ltest ) then
 
141
                print *,'ffxdx: id = ',id,', nevent = ',nevent
 
142
                print *,'ffxdx: lost ',ier1,' digits in D0 with isgnal '
 
143
     +                  ,isgnal,', trying other roots, isgnal ',-isgnal
 
144
                print *,'       if OK (no further messages) adding this'
 
145
     +                  ,' to your code will improve speed'
 
146
            endif
 
147
            isgnal = -isgnal
 
148
            ier0 = ier
 
149
            call ffxd0(cd0,xpi,ier0)
 
150
            isgnal = -isgnal
 
151
            if ( ier0.lt.ier1 ) then
 
152
                ier1 = ier0
 
153
                cdxi(1) = cd0
 
154
            endif
 
155
        endif
 
156
        if ( awrite ) then
 
157
            print *,'    '
 
158
            print *,'ffxdx : level 0: id,nevent ',id,nevent
 
159
            print *,'D0 =',cdxi(1),ier1
 
160
        endif
 
161
        if ( ier1 .gt. 10 ) then
 
162
            print *,'ffxdx: id = ',id,', nevent = ',nevent
 
163
            print *,'ffxdx: error: D0 not stable, lost ',ier1,' digits'
 
164
            print *,'       please try another permutation or contact',
 
165
     +          ' author (t19@nikhef.nl)'
 
166
            print *,'xpi = ',xpi
 
167
        endif
 
168
*       note that we may have lost another factor xloss**3 or so
 
169
        mdxi(1) = absc(cdxi(1))*DBLE(10)**mod(ier1,50)
 
170
*
 
171
        if (level .eq. 0) goto 990
 
172
*
 
173
*  #] get D0:
 
174
*  #[ need C-functions till c-level=(level-1):
 
175
        if ( level.gt.3 ) then
 
176
            print *,'ffxdx: error: higher than third rank ',
 
177
     +          'not yet implemented'
 
178
            stop
 
179
        endif
 
180
        cl = level-1
 
181
*       go trough the 4 different cancellation patterns
 
182
        if ( awrite ) then
 
183
            print '(a,i1)','####[ C-function output: to level ',cl
 
184
        endif
 
185
        do 100 i=1,4
 
186
            do 60 j=1,6
 
187
                xpi3(j) = xpi(iinx(j,i))
 
188
   60       continue
 
189
            ier0 = ier
 
190
            call ffxcx( caxj(3*i-2),maxj(3*i-2),cbxj(12*i-11),
 
191
     +          mbxj(12*i-11),ccxj(13*i-12),mcxj(13*i-12),d0,xmm,xpi3,
 
192
     +          cl,ier0)
 
193
            ier1 = max(ier1,ier0)
 
194
            fdel2i(i)=fdel2
 
195
  100   continue
 
196
        if ( awrite ) then
 
197
            print '(a)','####] C-function output:'
 
198
        endif
 
199
*  #] need C-functions till c-level=(level-1):
 
200
*  #[ break to let ffzdz tie in:
 
201
*
 
202
*       convert ??xj to ??xi
 
203
*
 
204
        call ffdji(ccxi,mcxi,cbxi,mbxi,caxi,maxi,
 
205
     +          ccxj,mcxj,cbxj,mbxj,caxj,maxj,level)
 
206
*
 
207
*       and call the real routine for the rest
 
208
*
 
209
        call ffxdxp(caxj,maxj,cbxj,mbxj,ccxj,mcxj,cdxi,mdxi,xpi,fdel2i,
 
210
     +          level,ier1)
 
211
*  #] break to let ffzdz tie in:
 
212
  990   ier = ier1
 
213
        end
 
214
        subroutine ffxdxp(caxj,maxj,cbxj,mbxj,ccxj,mcxj,cdxi,mdxi,xpi,
 
215
     +          fdel2i,level,ier)
 
216
*  #[ declarations :
 
217
        implicit none
 
218
*
 
219
*       arguments
 
220
*
 
221
        integer ier,level
 
222
        DOUBLE PRECISION xpi(13)
 
223
        DOUBLE COMPLEX caxj(12),cbxj(48),ccxj(52),cdxi(24)
 
224
        DOUBLE PRECISION maxj(12),mbxj(48),mcxj(52),mdxi(24),fdel2i(4)
 
225
*
 
226
*       local variables
 
227
*
 
228
        integer i,j,ier0,ier1,ier2,iinx(6,4),bij(12)
 
229
        DOUBLE PRECISION xi4(6),f1,f2,f3,absc,xmax,Rm(20:55),
 
230
     +          d11max,d22max,d23max,d24max,d0,xmm
 
231
        DOUBLE PRECISION dl2pij(6,6),del3p
 
232
        DOUBLE PRECISION mc0i(4),mxy(3),mc21i(4),mc22i(4),mc23i(4),
 
233
     +          mc24i(4),mc11i(4),mc12i(4),md1i(3)
 
234
        DOUBLE COMPLEX R20,R21,R22,R30,R31,R32,R33,R34,R35,R36,R37,R38,
 
235
     +          R41,R42,R43,R44,R45,R46,R47,R48,R49,R50,R51,R52,R53,R54,
 
236
     +          R55,cd1i(3),cc0i(4),cc11i(4),cc12i(4),
 
237
     +          cc21i(4),cc22i(4),cc23i(4),cc24i(4),cc,cxy(3)
 
238
        DOUBLE COMPLEX cd4pppp(3,3,3,3),cd4ppdel(3,3),cd4deldel,
 
239
     +          cd3ppp(3,3,3),cd3pdel(3),cd2pp(3,3),cd2del,
 
240
     +          cb0ij(4,4),ca0i(4),cd2(7)
 
241
        save iinx,bij
 
242
*
 
243
*       common blocks
 
244
*
 
245
        include 'ff.h'
 
246
        include 'aa.h'
 
247
*
 
248
*       data
 
249
*
 
250
        data iinx /2,3,4,6,7,10,
 
251
     +             1,3,4,9,7,8,
 
252
     +             1,2,4,5,10,8,
 
253
     +             1,2,3,5,6,9/
 
254
        data bij /1,2,5,6,9,10,17,18,21,22,33,34/
 
255
*
 
256
*       statement function
 
257
*
 
258
        absc(cc) = abs(DBLE(cc)) + abs(DIMAG(cc))
 
259
*
 
260
*  #] declarations :
 
261
*  #[ kinematical quantities for 4pv-red :
 
262
*       if ( abs(fdel3)  .lt. 1.d-6 ) then
 
263
*           print *,'kinematical det = 0, PV-scheme breaks down'
 
264
*           print *,fdel3
 
265
*           goto 990
 
266
*       endif
 
267
        if ( atest ) then
 
268
            del3p =
 
269
     +  - xpi(5)*xpi(5)*xpi(7) + xpi(5)*xpi(6)*xpi(7) + xpi(5)*xpi(6)*
 
270
     + xpi(8) - xpi(5)*xpi(6)*xpi(9) - xpi(5)*xpi(7)*xpi(7) + xpi(5)*
 
271
     + xpi(7)*xpi(8) + xpi(5)*xpi(7)*xpi(9) + xpi(5)*xpi(7)*xpi(10) -
 
272
     + xpi(5)*xpi(8)*xpi(10) + xpi(5)*xpi(9)*xpi(10) - xpi(6)*xpi(6)*
 
273
     + xpi(8) + xpi(6)*xpi(7)*xpi(8) - xpi(6)*xpi(7)*xpi(10) - xpi(6)*
 
274
     + xpi(8)*xpi(8) + xpi(6)*xpi(8)*xpi(9) + xpi(6)*xpi(8)*xpi(10) +
 
275
     + xpi(6)*xpi(9)*xpi(10) - xpi(7)*xpi(8)*xpi(9) + xpi(7)*xpi(9)*
 
276
     + xpi(10) + xpi(8)*xpi(9)*xpi(10) - xpi(9)*xpi(9)*xpi(10) -
 
277
     + xpi(9)*xpi(10)*xpi(10)
 
278
            del3p = del3p/4
 
279
            xmax = max(abs(xpi(5)),abs(xpi(6)),abs(xpi(7)),abs(xpi(8)),
 
280
     +          abs(xpi(9)),abs(xpi(10)))**3
 
281
            if ( abs(del3p-fdel3).gt.1d-12*xmax ) then
 
282
                print *,'ffxdxp: fdel3 wrong: ',fdel3,del3p,fdel3-del3p,
 
283
     +                  xmax
 
284
            endif
 
285
        endif
 
286
*
 
287
*       inverse kinematical matrix xi4  (3X3)
 
288
*
 
289
        ier2 = ier
 
290
        call aaxi4(xi4,ier2)
 
291
        ier2 = ier2 - ier
 
292
*
 
293
*           f-functions:
 
294
        f1 = 2*fpij4(1,5)
 
295
        f2 = 2*fpij4(1,6)
 
296
        f3 = 2*fpij4(1,7)
 
297
*
 
298
*  #] kinematical quantities for 4pv-red :
 
299
*  #[ level 1 : D11,D12,D13,C0(I)
 
300
*               need 4 diff C0(I)-functions,I=1,2,3
 
301
        cc0i(1)=ccxj(1)
 
302
        cc0i(2)=ccxj(14)
 
303
        cc0i(3)=ccxj(27)
 
304
        cc0i(4)=ccxj(40)
 
305
        mc0i(1)=mcxj(1)
 
306
        mc0i(2)=mcxj(14)
 
307
        mc0i(3)=mcxj(27)
 
308
        mc0i(4)=mcxj(40)
 
309
        ier1 = ier
 
310
        call ffxd1a(cdxi(2),mdxi(2),cdxi(1),mdxi(1),cc0i,mc0i,
 
311
     +          xpi,fpij4,fdel3,fdel2i,ier1)
 
312
        if ( awrite ) then
 
313
            print *,'GEERT JANs-scheme: id,nevent ',id,nevent
 
314
            print *,'D11=',cdxi(2),mdxi(2),ier1
 
315
            print *,'D12=',cdxi(3),mdxi(3),ier1
 
316
            print *,'D13=',cdxi(4),mdxi(4),ier1
 
317
        endif
 
318
*
 
319
        if ( lwarn .and. atest ) then
 
320
*           PV-reduction
 
321
            R20 = ( f1*cdxi(1)+cc0i(2)-cc0i(1) )/2
 
322
            R21 = ( f2*cdxi(1)+cc0i(3)-cc0i(2) )/2
 
323
            R22 = ( f3*cdxi(1)+cc0i(4)-cc0i(3) )/2
 
324
            Rm(20) = ( max(abs(f1)*mdxi(1),mc0i(2),mc0i(1)) )/2
 
325
            Rm(21) = ( max(abs(f2)*mdxi(1),mc0i(3),mc0i(2)) )/2
 
326
            Rm(22) = ( max(abs(f3)*mdxi(1),mc0i(4),mc0i(3)) )/2
 
327
            cd1i(1)=xi4(1)*R20+xi4(4)*R21+xi4(5)*R22
 
328
            cd1i(2)=xi4(4)*R20+xi4(2)*R21+xi4(6)*R22
 
329
            cd1i(3)=xi4(5)*R20+xi4(6)*R21+xi4(3)*R22
 
330
            md1i(1)=abs(xi4(1))*Rm(20)+abs(xi4(4))*Rm(21)+
 
331
     +          abs(xi4(5))*Rm(22)
 
332
            md1i(2)=abs(xi4(4))*Rm(20)+abs(xi4(2))*Rm(21)+
 
333
     +          abs(xi4(6))*Rm(22)
 
334
            md1i(3)=abs(xi4(5))*Rm(20)+abs(xi4(6))*Rm(21)+
 
335
     +          abs(xi4(3))*Rm(22)
 
336
            do 139 i=1,3
 
337
                if ( xloss**2*absc(cdxi(i+1)-cd1i(i)).gt.precc*max(
 
338
     +              mdxi(i+1),md1i(i)) ) print *,'ffxdx: error: FF D1',
 
339
     +              i,' disagrees with PV:',cdxi(i+1),cd1i(i),
 
340
     +              cdxi(i+1)-cd1i(i),max(mdxi(i+1),md1i(i))
 
341
  139       continue
 
342
            if (awrite) then
 
343
                print *,' '
 
344
                print *,'ffxdx : level 1: id,nevent ',id,nevent
 
345
                print *,'D11=',cd1i(1)
 
346
                print *,'D12=',cd1i(2)
 
347
                print *,'D13=',cd1i(3)
 
348
            endif
 
349
        endif
 
350
*
 
351
        if ( level.eq.1 ) then
 
352
            if ( lwarn ) then
 
353
                xmax = 0
 
354
                do i=1,4
 
355
                    if ( absc(cdxi(i)).ne.0 ) then
 
356
                        xmax = max(xmax,mdxi(i)/absc(cdxi(i)))
 
357
                    elseif ( mdxi(i).ne.0 ) then
 
358
                        xmax = max(xmax,1/precc)
 
359
                    endif
 
360
                enddo
 
361
                ier1 = int(log10(xmax))
 
362
                if ( awrite ) print *,'ier = ',ier1
 
363
            else
 
364
                ier1 = 0
 
365
            endif
 
366
            goto 990
 
367
        endif
 
368
*
 
369
*  #] level 1 :
 
370
*  #[ level 2 : D21,D22,D23,D24,D25,D26,D27,C11(I),C12(I)
 
371
*           need 4 diff C1-functions
 
372
        ier = ier1
 
373
        do 14 i=1,4
 
374
            j = 2 + (i-1)*13
 
375
            cc11i(i) = ccxj(j)
 
376
            cc12i(i) = ccxj(j+1)
 
377
            mc11i(i) = mcxj(j)
 
378
            mc12i(i) = mcxj(j+1)
 
379
   14   continue
 
380
*               PV-reduction
 
381
        cdxi(11)=-( f1*cdxi(2)+f2*cdxi(3)+f3*cdxi(4)-cc0i(1) )/2
 
382
     +           +xpi(1)*cdxi(1)
 
383
        if ( lwarn ) then
 
384
***         d11max = max(absc(f1*cdxi(2)),absc(f2*cdxi(3)),
 
385
***     +               absc(f3*cdxi(4)),absc(cc0i(1)),2*absc(xpi(1)*cdxi(1)))/2
 
386
***         if ( absc(cdxi(11)).lt.xloss*d11max ) then
 
387
***             ier0 = ier
 
388
***             call ffwarn(283,ier0,absc(cdxi(11)),d11max)
 
389
***             ier1 = max(ier1,ier0)
 
390
***         endif
 
391
            mdxi(11) = max(abs(f1)*mdxi(2),abs(f2)*mdxi(3),
 
392
     +          abs(f3)*mdxi(4),mc0i(1),2*xpi(1)*mdxi(1))/2
 
393
        endif
 
394
        R30=( f1*cdxi(2) + cc11i(2) + cc0i(1)  )/2 - cdxi(11)
 
395
        R31=( f2*cdxi(2) + cc11i(3) - cc11i(2) )/2
 
396
        R32=( f3*cdxi(2) + cc11i(4) - cc11i(3) )/2
 
397
        R33=( f1*cdxi(3) + cc11i(2) - cc11i(1) )/2
 
398
        R34=( f2*cdxi(3) + cc12i(3) - cc11i(2) )/2 - cdxi(11)
 
399
        R35=( f3*cdxi(3) + cc12i(4) - cc12i(3) )/2
 
400
        R36=( f1*cdxi(4) + cc12i(2) - cc12i(1) )/2
 
401
        R37=( f2*cdxi(4) + cc12i(3) - cc12i(2) )/2
 
402
        R38=( f3*cdxi(4)            - cc12i(3) )/2 - cdxi(11)
 
403
        cdxi(5) = xi4(1)*R30+xi4(4)*R31+xi4(5)*R32
 
404
        cdxi(6) = xi4(4)*R33+xi4(2)*R34+xi4(6)*R35
 
405
        cdxi(7) = xi4(5)*R36+xi4(6)*R37+xi4(3)*R38
 
406
        cdxi(8) = xi4(4)*R30+xi4(2)*R31+xi4(6)*R32
 
407
        cdxi(9) = xi4(5)*R30+xi4(6)*R31+xi4(3)*R32
 
408
        cdxi(10)= xi4(5)*R33+xi4(6)*R34+xi4(3)*R35
 
409
        if ( lwarn ) then
 
410
***         Rm(30)=max(absc(f1*cdxi(2)),absc(cc11i(2)),absc(cc0i(1)),
 
411
***     +               2*d11max)/2
 
412
***         Rm(31)=max(absc(f2*cdxi(2)),absc(cc11i(3)),absc(cc11i(2)))/2
 
413
***         Rm(32)=max(absc(f3*cdxi(2)),absc(cc11i(4)),absc(cc11i(3)))/2
 
414
***         Rm(33)=max(absc(f1*cdxi(3)),absc(cc11i(2)),absc(cc11i(1)))/2
 
415
***         Rm(34)=max(absc(f2*cdxi(3)),absc(cc12i(3)),absc(cc11i(2)),
 
416
***     +               2*d11max)/2
 
417
***         Rm(35)=max(absc(f3*cdxi(3)),absc(cc12i(4)),absc(cc12i(3)))/2
 
418
***         Rm(36)=max(absc(f1*cdxi(4)),absc(cc12i(2)),absc(cc12i(1)))/2
 
419
***         Rm(37)=max(absc(f2*cdxi(4)),absc(cc12i(3)),absc(cc12i(2)))/2
 
420
***         Rm(38)=max(absc(f3*cdxi(4)),absc(cc12i(3)),2*d11max)/2
 
421
***         xmax = max(abs(xi4(1))*Rm(30),abs(xi4(4))*Rm(31),abs(xi4(5))
 
422
***     +               *Rm(32))
 
423
***         if ( absc(cdxi(5)).lt.xloss*xmax ) then
 
424
***             ier0 = ier
 
425
***             call ffwarn(282,ier0,absc(cdxi(5)),xmax)
 
426
***             ier1 = max(ier1,ier0)
 
427
***         endif
 
428
***         xmax = max(abs(xi4(4))*Rm(33),abs(xi4(2))*Rm(34),abs(xi4(6))
 
429
***     +               *Rm(35))
 
430
***         if ( absc(cdxi(6)).lt.xloss*xmax ) then
 
431
***             ier0 = ier
 
432
***             call ffwarn(281,ier0,absc(cdxi(6)),xmax)
 
433
***             ier1 = max(ier1,ier0)
 
434
***         endif
 
435
***         xmax = max(abs(xi4(5))*Rm(36),abs(xi4(6))*Rm(37),abs(xi4(3))
 
436
***     +               *Rm(38))
 
437
***         if ( absc(cdxi(7)).lt.xloss*xmax ) then
 
438
***             ier0 = ier
 
439
***             call ffwarn(280,ier0,absc(cdxi(7)),xmax)
 
440
***             ier1 = max(ier1,ier0)
 
441
***         endif
 
442
***         xmax = max(abs(xi4(4))*Rm(30),abs(xi4(2))*Rm(31),abs(xi4(6))
 
443
***     +               *Rm(32))
 
444
***         if ( absc(cdxi(8)).lt.xloss*xmax ) then
 
445
***             ier0 = ier
 
446
***             call ffwarn(279,ier0,absc(cdxi(8)),xmax)
 
447
***             ier1 = max(ier1,ier0)
 
448
***         endif
 
449
***         xmax = max(abs(xi4(5))*Rm(30),abs(xi4(6))*Rm(31),abs(xi4(3))
 
450
***     +               *Rm(32))
 
451
***         if ( absc(cdxi(9)).lt.xloss*xmax ) then
 
452
***             ier0 = ier
 
453
***             call ffwarn(278,ier0,absc(cdxi(9)),xmax)
 
454
***             ier1 = max(ier1,ier0)
 
455
***         endif
 
456
***         xmax = max(abs(xi4(5))*Rm(33),abs(xi4(6))*Rm(34),abs(xi4(3))
 
457
***     +               *Rm(35))
 
458
***         if ( absc(cdxi(10)).lt.xloss*xmax ) then
 
459
***             ier0 = ier
 
460
***             call ffwarn(277,ier0,absc(cdxi(10)),xmax)
 
461
***             ier1 = max(ier1,ier0)
 
462
***         endif
 
463
*
 
464
*           the maximum values of the whole value (not only in this step)
 
465
*
 
466
            Rm(30) = max(abs(f1)*mdxi(2),mc11i(2),mc0i(1),2*mdxi(11))/2
 
467
            Rm(31) = max(abs(f2)*mdxi(2),mc11i(3),mc11i(2))/2
 
468
            Rm(32) = max(abs(f3)*mdxi(2),mc11i(4),mc11i(3))/2
 
469
            Rm(33) = max(abs(f1)*mdxi(3),mc11i(2),mc11i(1))/2
 
470
            Rm(34) = max(abs(f2)*mdxi(3),mc12i(3),mc11i(2),2*mdxi(11))/2
 
471
            Rm(35) = max(abs(f3)*mdxi(3),mc12i(4),mc12i(3))/2
 
472
            Rm(36) = max(abs(f1)*mdxi(4),mc12i(2),mc12i(1))/2
 
473
            Rm(37) = max(abs(f2)*mdxi(4),mc12i(3),mc12i(2))/2
 
474
            Rm(38) = max(abs(f3)*mdxi(4),mc12i(3),2*mdxi(11))/2
 
475
            mdxi(5) = max(abs(xi4(1))*Rm(30),abs(xi4(4))*Rm(31),
 
476
     +          abs(xi4(5))*Rm(32))
 
477
            mdxi(6) = max(abs(xi4(4))*Rm(33),abs(xi4(2))*Rm(34),
 
478
     +          abs(xi4(6))*Rm(35))
 
479
            mdxi(7) = max(abs(xi4(5))*Rm(36),abs(xi4(6))*Rm(37),
 
480
     +          abs(xi4(3))*Rm(38))
 
481
            mdxi(8) = max(abs(xi4(4))*Rm(30),abs(xi4(2))*Rm(31),
 
482
     +          abs(xi4(6))*Rm(32))
 
483
            mdxi(9) = max(abs(xi4(5))*Rm(30),abs(xi4(6))*Rm(31),
 
484
     +          abs(xi4(3))*Rm(32))
 
485
            mdxi(10)= max(abs(xi4(5))*Rm(33),abs(xi4(6))*Rm(34),
 
486
     +          abs(xi4(3))*Rm(35))
 
487
        endif
 
488
*       redundancy check
 
489
        if ( lwarn .and. atest ) then
 
490
            cxy(1) = xi4(1)*R33+xi4(4)*R34+xi4(5)*R35
 
491
            cxy(2) = xi4(1)*R36+xi4(4)*R37+xi4(5)*R38
 
492
            cxy(3) = xi4(4)*R36+xi4(2)*R37+xi4(6)*R38
 
493
            mxy(1) = abs(xi4(1))*Rm(33)+abs(xi4(4))*Rm(34)+abs(xi4(5))*
 
494
     +          Rm(35)
 
495
            mxy(2) = abs(xi4(1))*Rm(36)+abs(xi4(4))*Rm(37)+abs(xi4(5))*
 
496
     +          Rm(38)
 
497
            mxy(3) = abs(xi4(4))*Rm(36)+abs(xi4(2))*Rm(37)+abs(xi4(6))*
 
498
     +          Rm(38)
 
499
            if ( xloss*absc(cxy(1)-cdxi(8)) .gt.precc*max(mxy(1),
 
500
     +                  mdxi(8))
 
501
     +      .or. xloss*absc(cxy(2)-cdxi(9)) .gt.precc*max(mxy(2),
 
502
     +                  mdxi(9))
 
503
     +      .or. xloss*absc(cxy(3)-cdxi(10)).gt.precc*max(mxy(3),
 
504
     +                  mdxi(10)) ) then
 
505
                print *,'ffxdx: error: id/nevent ',id,'/',nevent
 
506
                print *,'redundancy check at level 2 failed: '
 
507
                print *,cxy(1),cdxi(8),absc(cxy(1)-cdxi(8)),
 
508
     +                  max(mxy(1),mdxi(8))
 
509
                print *,cxy(2),cdxi(9),absc(cxy(2)-cdxi(9)),
 
510
     +                  max(mxy(2),mdxi(9))
 
511
                print *,cxy(3),cdxi(10),absc(cxy(3)-cdxi(20)),
 
512
     +                  max(mxy(3),mdxi(10))
 
513
            endif
 
514
        endif
 
515
        if ( awrite ) then
 
516
            print *,' '
 
517
            print *,'ffxdx : level 2: id,nevent ',id,nevent
 
518
            print *,'D21=',cdxi(5),mdxi(5),ier1
 
519
            print *,'D22=',cdxi(6),mdxi(6),ier1
 
520
            print *,'D23=',cdxi(7),mdxi(7),ier1
 
521
            print *,'D24=',cdxi(8),mdxi(8),ier1
 
522
            print *,'    ',cxy(1),mxy(1)
 
523
            print *,'D25=',cdxi(9),mdxi(9),ier1
 
524
            print *,'    ',cxy(2),mxy(2)
 
525
            print *,'D26=',cdxi(10),mdxi(10),ier1
 
526
            print *,'    ',cxy(3),mxy(3)
 
527
            print *,'D27=',cdxi(11),mdxi(11),ier1
 
528
        endif
 
529
*       this goes wrong in the case of complex masses - no way out yet...
 
530
        if ( awrite .and. .FALSE. ) then
 
531
            d0 = 0
 
532
            xmm = 0
 
533
            if ( awrite ) print *,'calling ffxdi with ier = ',ier
 
534
*           the order of the B0s can be deduced from the C0 -> B0 chain
 
535
            cb0ij(1,2) = cbxj(33)
 
536
            cb0ij(1,3) = cbxj(21)
 
537
            cb0ij(1,4) = cbxj(17)
 
538
            cb0ij(2,1) = cbxj(33)
 
539
            cb0ij(2,3) = cbxj( 9)
 
540
            cb0ij(2,4) = cbxj( 5)
 
541
            cb0ij(3,1) = cbxj(21)
 
542
            cb0ij(3,2) = cbxj( 9)
 
543
            cb0ij(3,4) = cbxj( 1)
 
544
            cb0ij(4,1) = cbxj(17)
 
545
            cb0ij(4,2) = cbxj( 5)
 
546
            cb0ij(4,3) = cbxj( 1)
 
547
*           the A0s are not used for the moment
 
548
            call ffxdi(cd4pppp,cd4ppdel,cd4deldel, cd3ppp,cd3pdel,
 
549
     +          cd2pp,cd2del, cd1i, dl2pij, cdxi(1),cc0i,cb0ij,ca0i,
 
550
     +          fdel4s,fdel3,fdel2i, xpi,fpij4, d0,xmm, 2, ier)
 
551
*  #[ convert to PV conventions:
 
552
*
 
553
        ier1 = ier
 
554
        cd2(1) = cd2pp(1,1) - DBLE(fdel2i(1))*cd2del
 
555
        if ( lwarn .and. absc(cd2(1)).lt.xloss*absc(cd2pp(1,1)) ) then
 
556
            call ffwarn(229,ier1,absc(cd2(1)),absc(cd2pp(1,1)))
 
557
        endif
 
558
        cd2(2) = cd2pp(1,2) + DBLE(dl2pij(2,4))*cd2del
 
559
        if ( lwarn .and. absc(cd2(2)).lt.xloss*absc(cd2pp(1,2)) ) then
 
560
            ier0 = ier
 
561
            call ffwarn(229,ier0,absc(cd2(2)),absc(cd2pp(1,2)))
 
562
            ier1 = max(ier1,ier0)
 
563
        endif
 
564
        cd2(3) = cd2pp(1,3) - DBLE(dl2pij(1,4))*cd2del
 
565
        if ( lwarn .and. absc(cd2(3)).lt.xloss*absc(cd2pp(1,3)) ) then
 
566
            ier0 = ier
 
567
            call ffwarn(229,ier0,absc(cd2(3)),absc(cd2pp(1,3)))
 
568
            ier1 = max(ier1,ier0)
 
569
        endif
 
570
        cd2(4) = cd2pp(2,2) - DBLE(xpi(5)*xpi(7)-fpij4(5,7)**2)*cd2del
 
571
        if ( lwarn .and. absc(cd2(4)).lt.xloss*absc(cd2pp(2,2)) ) then
 
572
            ier0 = ier
 
573
            call ffwarn(229,ier0,absc(cd2(4)),absc(cd2pp(2,2)))
 
574
            ier1 = max(ier1,ier0)
 
575
        endif
 
576
        cd2(5) = cd2pp(2,3) + DBLE(dl2pij(1,2))*cd2del
 
577
        if ( lwarn .and. absc(cd2(5)).lt.xloss*absc(cd2pp(2,3)) ) then
 
578
            ier0 = ier
 
579
            call ffwarn(229,ier0,absc(cd2(5)),absc(cd2pp(2,3)))
 
580
            ier1 = max(ier1,ier0)
 
581
        endif
 
582
        cd2(6) = cd2pp(3,3) - DBLE(fdel2i(4))*cd2del
 
583
        if ( lwarn .and. absc(cd2(6)).lt.xloss*absc(cd2pp(3,3)) ) then
 
584
            ier0 = ier
 
585
            call ffwarn(229,ier0,absc(cd2(6)),absc(cd2pp(3,3)))
 
586
            ier1 = max(ier1,ier0)
 
587
        endif
 
588
        cd2(7) = DBLE(fdel3)*cd2del
 
589
*
 
590
*  #] convert to PV conventions:
 
591
            if ( awrite ) then
 
592
                print *,'ffxdi gives'
 
593
                print *,'D11 = ',cd1i(1),ier1
 
594
                print *,'D12 = ',cd1i(2),ier1
 
595
                print *,'D13 = ',cd1i(3),ier1
 
596
                print *,'D21 = ',cd2(1),ier1
 
597
                print *,'D22 = ',cd2(4),ier1
 
598
                print *,'D23 = ',cd2(6),ier1
 
599
                print *,'D24 = ',cd2(2),ier1
 
600
                print *,'D25 = ',cd2(3),ier1
 
601
                print *,'D26 = ',cd2(5),ier1
 
602
                print *,'D27 = ',cd2(7),ier1
 
603
            endif
 
604
        endif
 
605
*
 
606
        if ( level.eq.2 ) then
 
607
            if ( lwarn ) then
 
608
                xmax = 0
 
609
                do i=1,11
 
610
                    if ( absc(cdxi(i)).ne.0 ) then
 
611
                        xmax = max(xmax,mdxi(i)/absc(cdxi(i)))
 
612
                    elseif ( mdxi(i).ne.0 ) then
 
613
                        xmax = max(xmax,1/precc)
 
614
                    endif
 
615
                enddo
 
616
                ier1 = int(log10(xmax))
 
617
                if ( awrite ) print *,'ier = ',ier1
 
618
            else
 
619
                ier1 = 0
 
620
            endif
 
621
            goto 990
 
622
        endif
 
623
*
 
624
*  #] level 2 :
 
625
*  #[ level 3 : D31,D32,D33,D34,D35,D36,D37,D38,D39,D310,D311,D312,D313
 
626
*               C21(I),C22(I),C23(I),C11(I),C12(I)
 
627
*               need 4 diff C2-functions
 
628
        do 15 i=1,4
 
629
            j = 4 +(i-1)*13
 
630
            cc21i(i)=ccxj(j)
 
631
            cc22i(i)=ccxj(j+1)
 
632
            cc23i(i)=ccxj(j+2)
 
633
            cc24i(i)=ccxj(j+3)
 
634
            mc21i(i)=mcxj(j)
 
635
            mc22i(i)=mcxj(j+1)
 
636
            mc23i(i)=mcxj(j+2)
 
637
            mc24i(i)=mcxj(j+3)
 
638
   15   continue
 
639
*               PV-reduction
 
640
        R53=( f1*cdxi(11) + cc24i(2) - cc24i(1) )/2
 
641
        R54=( f2*cdxi(11) + cc24i(3) - cc24i(2) )/2
 
642
        R55=( f3*cdxi(11) + cc24i(4) - cc24i(3) )/2
 
643
        cdxi(22) = xi4(1)*R53+xi4(4)*R54+xi4(5)*R55
 
644
        cdxi(23) = xi4(4)*R53+xi4(2)*R54+xi4(6)*R55
 
645
        cdxi(24) = xi4(5)*R53+xi4(6)*R54+xi4(3)*R55
 
646
        if ( lwarn ) then
 
647
***         Rm(53)=max(abs(f1)*d11max,absc(cc24i(2)),absc(cc24i(1)))/2
 
648
***         Rm(54)=max(abs(f2)*d11max,absc(cc24i(3)),absc(cc24i(2)))/2
 
649
***         Rm(55)=max(abs(f3)*d11max,absc(cc24i(4)),absc(cc24i(3)))/2
 
650
***         d22max = max(abs(xi4(1))*Rm(53),abs(xi4(4))*Rm(54),
 
651
***     +               abs(xi4(5))*Rm(55))
 
652
***         if ( absc(cdxi(22)).lt.xloss*d22max ) then
 
653
***             ier0 = ier
 
654
***             call ffwarn(276,ier0,absc(cdxi(22)),d22max)
 
655
***             ier1 = max(ier1,ier0)
 
656
***         endif
 
657
***         d23max = max(abs(xi4(4))*Rm(53),abs(xi4(2))*Rm(54),
 
658
***     +               abs(xi4(6))*Rm(55))
 
659
***         if ( absc(cdxi(23)).lt.xloss*d23max ) then
 
660
***             ier0 = ier
 
661
***             call ffwarn(275,ier0,absc(cdxi(23)),d23max)
 
662
***             ier1 = max(ier1,ier0)
 
663
***         endif
 
664
***         d24max = max(abs(xi4(5))*Rm(53),abs(xi4(6))*Rm(54),
 
665
***     +               abs(xi4(3))*Rm(55))
 
666
***         if ( absc(cdxi(24)).lt.xloss*xmax ) then
 
667
***             ier0 = ier
 
668
***             call ffwarn(274,ier0,absc(cdxi(24)),d24max)
 
669
***             ier1 = max(ier1,ier0)
 
670
***         endif
 
671
*
 
672
*           and again for the whole thing
 
673
*
 
674
            Rm(53)=max(abs(f1)*mdxi(11),mc24i(2),mc24i(1))/2
 
675
            Rm(54)=max(abs(f2)*mdxi(11),mc24i(3),mc24i(2))/2
 
676
            Rm(55)=max(abs(f3)*mdxi(11),mc24i(4),mc24i(3))/2
 
677
            mdxi(22) = max(abs(xi4(1))*Rm(53),abs(xi4(4))*Rm(54),
 
678
     +          abs(xi4(5))*Rm(55))
 
679
            mdxi(23) = max(abs(xi4(4))*Rm(53),abs(xi4(2))*Rm(54),
 
680
     +          abs(xi4(6))*Rm(55))
 
681
            mdxi(24) = max(abs(xi4(5))*Rm(53),abs(xi4(6))*Rm(54),
 
682
     +          abs(xi4(3))*Rm(55))
 
683
        endif
 
684
*
 
685
        R41=( f1*cdxi(5) + cc21i(2) - cc0i(1)  )/2-2*cdxi(22)
 
686
        R42=( f2*cdxi(5) + cc21i(3) - cc21i(2) )/2
 
687
        R43=( f3*cdxi(5) + cc21i(4) - cc21i(3) )/2
 
688
        R44=( f1*cdxi(6) + cc21i(2) - cc21i(1) )/2
 
689
        R45=( f2*cdxi(6) + cc22i(3) - cc21i(2) )/2-2*cdxi(23)
 
690
        R46=( f3*cdxi(6) + cc22i(4) - cc22i(3) )/2
 
691
        R47=( f1*cdxi(7) + cc22i(2) - cc22i(1) )/2
 
692
        R48=( f2*cdxi(7) + cc22i(3) - cc22i(2) )/2
 
693
        R49=( f3*cdxi(7)            - cc22i(3) )/2-2*cdxi(24)
 
694
        R50=( f1*cdxi(8) + cc21i(2) + cc11i(1) )/2-cdxi(23)
 
695
        R51=( f2*cdxi(8) + cc23i(3) - cc21i(2) )/2-cdxi(22)
 
696
        R52=( f3*cdxi(8) + cc23i(4) - cc23i(3) )/2
 
697
        cdxi(12) = xi4(1)*R41+xi4(4)*R42+xi4(5)*R43
 
698
        cdxi(13) = xi4(4)*R44+xi4(2)*R45+xi4(6)*R46
 
699
        cdxi(14) = xi4(5)*R47+xi4(6)*R48+xi4(3)*R49
 
700
        cdxi(15) = xi4(4)*R41+xi4(2)*R42+xi4(6)*R43
 
701
        cdxi(16) = xi4(5)*R41+xi4(6)*R42+xi4(3)*R43
 
702
        cdxi(17) = xi4(1)*R44+xi4(4)*R45+xi4(5)*R46
 
703
        cdxi(18) = xi4(1)*R47+xi4(4)*R48+xi4(5)*R49
 
704
        cdxi(19) = xi4(5)*R44+xi4(6)*R45+xi4(3)*R46
 
705
        cdxi(20) = xi4(4)*R47+xi4(2)*R48+xi4(6)*R49
 
706
        cdxi(21) = xi4(5)*R50+xi4(6)*R51+xi4(3)*R52
 
707
        if ( lwarn ) then
 
708
***         Rm(41)=max(absc(f1*cdxi(5)),absc(cc21i(2)),absc(cc0i(1)),
 
709
***     +               4*d22max)/2
 
710
***         Rm(42)=max(absc(f2*cdxi(5)),absc(cc21i(3)),absc(cc21i(2)))/2
 
711
***         Rm(43)=max(absc(f3*cdxi(5)),absc(cc21i(4)),absc(cc21i(3)))/2
 
712
***         Rm(44)=max(absc(f1*cdxi(6)),absc(cc21i(2)),absc(cc21i(1)))/2
 
713
***         Rm(45)=max(absc(f2*cdxi(6)),absc(cc22i(3)),absc(cc21i(2)),
 
714
***     +               4*d23max)/2
 
715
***         Rm(46)=max(absc(f3*cdxi(6)),absc(cc22i(4)),absc(cc22i(3)))/2
 
716
***         Rm(47)=max(absc(f1*cdxi(7)),absc(cc22i(2)),absc(cc22i(1)))/2
 
717
***         Rm(48)=max(absc(f2*cdxi(7)),absc(cc22i(3)),absc(cc22i(2)))/2
 
718
***         Rm(49)=max(absc(f3*cdxi(7)),absc(cc22i(3)),4*d24max)/2
 
719
***         Rm(50)=max(absc(f1*cdxi(8)),absc(cc21i(2)),absc(cc11i(1)),
 
720
***     +               2*d23max)/2
 
721
***         Rm(51)=max(absc(f2*cdxi(8)),absc(cc23i(3)),absc(cc21i(2)),
 
722
***     +               2*d22max)/2
 
723
***         Rm(52)=max(absc(f3*cdxi(8)),absc(cc23i(4)),absc(cc23i(3)))/2
 
724
***         xmax=max(abs(xi4(1))*Rm(41),abs(xi4(4))*Rm(42),
 
725
***     +               abs(xi4(5))*Rm(43))
 
726
***         if ( absc(cdxi(12)).lt.xloss*xmax ) then
 
727
***             ier0 = ier
 
728
***             call ffwarn(273,ier0,absc(cdxi(12)),xmax)
 
729
***             ier1 = max(ier1,ier0)
 
730
***         endif
 
731
***         xmax=max(abs(xi4(4))*Rm(44),abs(xi4(2))*Rm(45),
 
732
***     +               abs(xi4(6))*Rm(46))
 
733
***         if ( absc(cdxi(13)).lt.xloss*xmax ) then
 
734
***             ier0 = ier
 
735
***             call ffwarn(272,ier0,absc(cdxi(13)),xmax)
 
736
***             ier1 = max(ier1,ier0)
 
737
***         endif
 
738
***         xmax=max(abs(xi4(5))*Rm(47),abs(xi4(6))*Rm(48),
 
739
***     +               abs(xi4(3))*Rm(49))
 
740
***         if ( absc(cdxi(14)).lt.xloss*xmax ) then
 
741
***             ier0 = ier
 
742
***             call ffwarn(271,ier0,absc(cdxi(14)),xmax)
 
743
***             if ( awrite ) then
 
744
***             print *,'xi4(5)*R47,xi4(6)*R48,xi4(3)*R49,cdxi(14) = '
 
745
***             print *,xi4(5)*R47,xi4(6)*R48,xi4(3)*R49,cdxi(14)
 
746
***             print *,xi4(5)*Rm(47),xi4(6)*Rm(48),xi4(3)*Rm(49),xmax
 
747
***             endif
 
748
***             ier1 = max(ier1,ier0)
 
749
***         endif
 
750
***         xmax=max(abs(xi4(4))*Rm(41),abs(xi4(2))*Rm(42),
 
751
***     +               abs(xi4(6))*Rm(43))
 
752
***         if ( absc(cdxi(15)).lt.xloss*xmax ) then
 
753
***             ier0 = ier
 
754
***             call ffwarn(270,ier0,absc(cdxi(15)),xmax)
 
755
***             ier1 = max(ier1,ier0)
 
756
***         endif
 
757
***         xmax=max(abs(xi4(5))*Rm(41),abs(xi4(6))*Rm(42),
 
758
***     +               abs(xi4(3))*Rm(43))
 
759
***         if ( absc(cdxi(16)).lt.xloss*xmax ) then
 
760
***             ier0 = ier
 
761
***             call ffwarn(269,ier0,absc(cdxi(16)),xmax)
 
762
***             ier1 = max(ier1,ier0)
 
763
***         endif
 
764
***         xmax=max(abs(xi4(1))*Rm(44),abs(xi4(4))*Rm(45),
 
765
***     +               abs(xi4(5))*Rm(46))
 
766
***         if ( absc(cdxi(17)).lt.xloss*xmax ) then
 
767
***             ier0 = ier
 
768
***             call ffwarn(268,ier0,absc(cdxi(17)),xmax)
 
769
***             ier1 = max(ier1,ier0)
 
770
***         endif
 
771
***         xmax=max(abs(xi4(1))*Rm(47),abs(xi4(4))*Rm(48),
 
772
***     +               abs(xi4(5))*Rm(49))
 
773
***         if ( absc(cdxi(18)).lt.xloss*xmax ) then
 
774
***             ier0 = ier
 
775
***             call ffwarn(267,ier0,absc(cdxi(18)),xmax)
 
776
***             ier1 = max(ier1,ier0)
 
777
***         endif
 
778
***         xmax=max(abs(xi4(5))*Rm(44),abs(xi4(6))*Rm(45),
 
779
***     +               abs(xi4(3))*Rm(46))
 
780
***         if ( absc(cdxi(19)).lt.xloss*xmax ) then
 
781
***             ier0 = ier
 
782
***             call ffwarn(266,ier0,absc(cdxi(19)),xmax)
 
783
***             ier1 = max(ier1,ier0)
 
784
***         endif
 
785
***         xmax=max(abs(xi4(4))*Rm(47),abs(xi4(2))*Rm(48),
 
786
***     +               abs(xi4(6))*Rm(49))
 
787
***         if ( absc(cdxi(20)).lt.xloss*xmax ) then
 
788
***             ier0 = ier
 
789
***             call ffwarn(265,ier0,absc(cdxi(20)),xmax)
 
790
***             if ( awrite ) then
 
791
***             print *,'xi4(4)*R47,xi4(2)*R48,xi4(6)*R49,cdxi(20) = '
 
792
***             print *,xi4(4)*R47,xi4(2)*R48,xi4(6)*R49,cdxi(20)
 
793
***             print *,xi4(4)*Rm(47),xi4(2)*Rm(48),xi4(6)*Rm(49),xmax
 
794
***             endif
 
795
***             ier1 = max(ier1,ier0)
 
796
***         endif
 
797
***         xmax=max(abs(xi4(5))*Rm(50),abs(xi4(6))*Rm(51),
 
798
***     +               abs(xi4(3))*Rm(52))
 
799
***         if ( absc(cdxi(21)).lt.xloss*xmax ) then
 
800
***             ier0 = ier
 
801
***             call ffwarn(264,ier0,absc(cdxi(21)),xmax)
 
802
***             ier1 = max(ier1,ier0)
 
803
***         endif
 
804
*
 
805
*           again the whole thing, not just this step
 
806
*
 
807
            Rm(41) = max(abs(f1)*mdxi(5),mc21i(2),mc0i(1),4*mdxi(22))/2
 
808
            Rm(42) = max(abs(f2)*mdxi(5),mc21i(3),mc21i(2))/2
 
809
            Rm(43) = max(abs(f3)*mdxi(5),mc21i(4),mc21i(3))/2
 
810
            Rm(44) = max(abs(f1)*mdxi(6),mc21i(2),mc21i(1))/2
 
811
            Rm(45) = max(abs(f2)*mdxi(6),mc22i(3),mc21i(2),4*mdxi(23))/2
 
812
            Rm(46) = max(abs(f3)*mdxi(6),mc22i(4),mc22i(3))/2
 
813
            Rm(47) = max(abs(f1)*mdxi(7),mc22i(2),mc22i(1))/2
 
814
            Rm(48) = max(abs(f2)*mdxi(7),mc22i(3),mc22i(2))/2
 
815
            Rm(49) = max(abs(f3)*mdxi(7),mc22i(3),4*mdxi(24))/2
 
816
            Rm(50) = max(abs(f1)*mdxi(8),mc21i(2),mc11i(1),2*mdxi(23))/2
 
817
            Rm(51) = max(abs(f2)*mdxi(8),mc23i(3),mc21i(2),2*mdxi(22))/2
 
818
            Rm(52) = max(abs(f3)*mdxi(8),mc23i(4),mc23i(3))/2
 
819
            mdxi(12) = max(abs(xi4(1))*Rm(41),abs(xi4(4))*Rm(42),
 
820
     +          abs(xi4(5))*Rm(43))
 
821
            mdxi(13) = max(abs(xi4(4))*Rm(44),abs(xi4(2))*Rm(45),
 
822
     +          abs(xi4(6))*Rm(46))
 
823
            mdxi(14) = max(abs(xi4(5))*Rm(47),abs(xi4(6))*Rm(48),
 
824
     +          abs(xi4(3))*Rm(49))
 
825
            mdxi(15) = max(abs(xi4(4))*Rm(41),abs(xi4(2))*Rm(42),
 
826
     +          abs(xi4(6))*Rm(43))
 
827
            mdxi(16) = max(abs(xi4(5))*Rm(41),abs(xi4(6))*Rm(42),
 
828
     +          abs(xi4(3))*Rm(43))
 
829
            mdxi(17) = max(abs(xi4(1))*Rm(44),abs(xi4(4))*Rm(45),
 
830
     +          abs(xi4(5))*Rm(46))
 
831
            mdxi(18) = max(abs(xi4(1))*Rm(47),abs(xi4(4))*Rm(48),
 
832
     +          abs(xi4(5))*Rm(49))
 
833
            mdxi(19) = max(abs(xi4(5))*Rm(44),abs(xi4(6))*Rm(45),
 
834
     +          abs(xi4(3))*Rm(46))
 
835
            mdxi(20) = max(abs(xi4(4))*Rm(47),abs(xi4(2))*Rm(48),
 
836
     +          abs(xi4(6))*Rm(49))
 
837
            mdxi(21) = max(abs(xi4(5))*Rm(50),abs(xi4(6))*Rm(51),
 
838
     +          abs(xi4(3))*Rm(52))
 
839
        endif
 
840
*       redundancy check
 
841
        if ( lwarn .and. atest ) then
 
842
            cxy(1) = xi4(1)*R50+xi4(4)*R51+xi4(5)*R52
 
843
            cxy(2) = xi4(4)*R50+xi4(2)*R51+xi4(6)*R52
 
844
            mxy(1) = abs(xi4(1))*Rm(50)+abs(xi4(4))*Rm(51)+ abs(xi4(5))*
 
845
     +          Rm(52)
 
846
            mxy(2) = abs(xi4(4))*Rm(50)+abs(xi4(2))*Rm(51)+ abs(xi4(6))*
 
847
     +          Rm(52)
 
848
            if ( xloss*absc(cxy(1)-cdxi(15)).gt.precc*max(mxy(1),
 
849
     +                  mdxi(15))
 
850
     +      .or. xloss*absc(cxy(2)-cdxi(17)).gt.precc*max(mxy(2),
 
851
     +                  mdxi(17))  ) then
 
852
                print *,'ffxdx: error: id/nevent ',id,'/',nevent
 
853
                print *,'redundancy check at level 3 failed: '
 
854
                print *,cxy(1),cdxi(15),absc(cxy(1)-cdxi(15)),
 
855
     +                  max(mxy(1),mdxi(15))
 
856
                print *,cxy(2),cdxi(17),absc(cxy(2)-cdxi(17)),
 
857
     +                  max(mxy(2),mdxi(17))
 
858
            endif
 
859
        endif
 
860
        if ( awrite ) then
 
861
            print *,' '
 
862
            print *,'ffxdx : level 3: id,nevent ',id,nevent
 
863
            print *,'D31 =',cdxi(12),mdxi(12),ier1
 
864
            print *,'D32 =',cdxi(13),mdxi(13),ier1
 
865
            print *,'D33 =',cdxi(14),mdxi(14),ier1
 
866
            print *,'D34 =',cdxi(15),mdxi(15),ier1
 
867
            print *,'     ',cxy(1)  ,mxy(1)
 
868
            print *,'D35 =',cdxi(16),mdxi(16),ier1
 
869
            print *,'D36 =',cdxi(17),mdxi(17),ier1
 
870
            print *,'     ',cxy(2)  ,mxy(2)
 
871
            print *,'D37 =',cdxi(18),mdxi(18),ier1
 
872
            print *,'D38 =',cdxi(19),mdxi(19),ier1
 
873
            print *,'D39 =',cdxi(20),mdxi(20),ier1
 
874
            print *,'D310=',cdxi(21),mdxi(21),ier1
 
875
            print *,'D311=',cdxi(22),mdxi(22),ier1
 
876
            print *,'D312=',cdxi(23),mdxi(23),ier1
 
877
            print *,'D313=',cdxi(24),mdxi(24),ier1
 
878
        endif
 
879
*
 
880
        if ( level.eq.3 ) then
 
881
            if ( lwarn ) then
 
882
                xmax = 0
 
883
                do i=1,24
 
884
                    if ( absc(cdxi(i)).ne.0 ) then
 
885
                        xmax = max(xmax,mdxi(i)/absc(cdxi(i)))
 
886
                    elseif ( mdxi(i).ne.0 ) then
 
887
                        xmax = max(xmax,1/precc)
 
888
                    endif
 
889
                enddo
 
890
                ier1 = int(log10(xmax))
 
891
                if ( awrite ) print *,'ier = ',ier1
 
892
            else
 
893
                ier1 = 0
 
894
            endif
 
895
            goto 990
 
896
        endif
 
897
*
 
898
*  #] level 3 :
 
899
*  #[ end:
 
900
        print *,'ffxdx: level ',level,' not supported.'
 
901
        stop
 
902
  990   continue
 
903
        ier = ier1 + ier2
 
904
*  #] end:
 
905
*###] ffxdx:
 
906
        end
 
907
*###[ ffdji:
 
908
        subroutine ffdji(ccxi,mcxi,cbxi,mbxi,caxi,maxi,
 
909
     +          ccxj,mcxj,cbxj,mbxj,caxj,maxj,level)
 
910
***#[*comment:***********************************************************
 
911
*                                                                       *
 
912
*       Renumber the [mc][abc]xj arrays into the [mc][abc]xi arrays.    *
 
913
*       Note: the A's are not yet used and not yet renumbered!          *
 
914
*                                                                       *
 
915
***#]*comment:***********************************************************
 
916
*  #[ declarations:
 
917
        implicit none
 
918
*
 
919
*       arguments
 
920
*
 
921
        integer level
 
922
        DOUBLE PRECISION mcxi(28),mbxi(12),maxi(4),
 
923
     +          mcxj(52),mbxj(48),maxj(12)
 
924
        DOUBLE COMPLEX ccxi(28),cbxi(12),caxi(4),
 
925
     +          ccxj(52),cbxj(48),caxj(12)
 
926
*
 
927
*       local variables
 
928
*
 
929
        integer i,j,k,bij(12),beq(6,2)
 
930
        save bij,beq
 
931
*
 
932
*       common
 
933
*
 
934
        include 'aa.h'
 
935
        include 'ff.h'
 
936
*
 
937
*       data
 
938
*
 
939
        data bij /1,2,5,6,9,10,17,18,21,22,33,34/
 
940
        data beq / 0, 4, 8,16,20,32,
 
941
     +            12,24,36,28,40,44/
 
942
*
 
943
*  #] declarations:
 
944
*  #[ renumber:
 
945
*       output preparation
 
946
*          1)C-output: reduce the array ccxj(4*13) to ccxi(4*7)
 
947
*                      c's are calculated only to (level-1)
 
948
        do 130 j=1,4
 
949
            do 131 i=1,7
 
950
                ccxi(i+(j-1)*7) = ccxj(i+(j-1)*13)
 
951
                mcxi(i+(j-1)*7) = mcxj(i+(j-1)*13)
 
952
 131        continue
 
953
 130    continue
 
954
*          2)B-output: reduce the array cbxj(12*4) to cbxi(6*2)
 
955
*                      b's are calculated only to (level-2)
 
956
        do i=1,12
 
957
            cbxi(i) = cbxj(bij(i))
 
958
            mbxi(i) = mbxj(bij(i))
 
959
        enddo
 
960
*       check the symmetry in B0(i,j)
 
961
        if ( atest ) then
 
962
            do 13 i=1,4
 
963
                do 12 j=1,6
 
964
                    if ( xloss*abs(cbxj(i+beq(j,1))-cbxj(i+beq(j,2)))
 
965
     +                          .gt. precc*abs(cbxj(i+beq(j,1))) ) then
 
966
                        print *,'ffxdji: cbxj(',i+beq(j,1),') != cbxj(',
 
967
     +                          i+beq(j,2),'): ',cbxj(i+beq(j,1)),
 
968
     +                          cbxj(i+beq(j,2)),cbxj(i+beq(j,1))-
 
969
     +                          cbxj(i+beq(j,2))
 
970
                    endif
 
971
   12           continue
 
972
   13       continue
 
973
        endif
 
974
*  #] renumber:
 
975
*###] ffdji:
 
976
        end