JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
mod_HiggsJJ.F90
Go to the documentation of this file.
1 module modhiggsjj
2  use modparameters
3  use modmisc
4  implicit none
5  private
6 
7  public :: evalamp_wbfh_unsymm_sa,evalamp_wbfh_unsymm_sa_select,evalamp_wbfh_unsymm_sa_select_exact
8  public :: evalamp_sbfh_unsymm_sa,evalamp_sbfh_unsymm_sa_select,evalamp_sbfh_unsymm_sa_select_exact
9  !public :: wrapHVV
10 
11  !-- general definitions, to be merged with Markus final structure
12  real(dp), public, parameter :: tag1 = 1.0_dp
13  real(dp), public, parameter :: tag2 = 1.0_dp
14  real(dp), public, parameter :: tagbot = 1.0_dp
15 
16 
17  CONTAINS
18 
19 
20 
21  !-- SM: |g2| = alphas/(six*pi), so multiply ME**2 by (|g2|*gs)**2
22  !-- g3 not supported yet
23  subroutine evalamp_sbfh_unsymm_sa(p,res)
24  real(dp), intent(in) :: p(4,5)
25  real(dp), intent(out) :: res(-5:5,-5:5)
26  real(dp) :: sprod(4,4)
27  complex(dp) :: za(4,4), zb(4,4)
28  real(dp) :: restmp, restmpid
29  integer :: i, j
30 
31  res = zero
32 
33  call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
34 
35  !-- gg -> gg
36  call me2_ggggh(1,2,3,4,za,zb,sprod,restmp)
37  restmp = restmp * avegg * symmfac
38  res(0,0) = res(0,0) + restmp
39 
40  !-- gg -> qqb
41  call me2_qbqggh(4,3,1,2,za,zb,sprod,restmp)
42  restmp = restmp * avegg
43  res(0,0) = res(0,0) + restmp * nf
44 
45  !-- gq -> gq
46  call me2_qbqggh(2,4,1,3,za,zb,sprod,restmp)
47  restmp = restmp * aveqg
48  do i = 1,5
49  res(0,i) = res(0,i) + restmp
50  res(0,-i) = res(0,-i) + restmp
51  enddo
52  call me2_qbqggh(1,4,2,3,za,zb,sprod,restmp)
53  restmp = restmp * aveqg
54  do i = 1,5
55  res(i,0) = res(i,0) + restmp
56  res(-i,0) = res(-i,0) + restmp
57  enddo
58 
59  !-- qqb -> gg
60  call me2_qbqggh(1,2,3,4,za,zb,sprod,restmp)
61  restmp = restmp * aveqq * symmfac
62  do i = 1,5
63  res(i,-i) = res(i,-i) + restmp
64  res(-i,i) = res(-i,i) + restmp
65  enddo
66 
67  !-- qqb -> qqb
68  call me2_qbqqbq(1,2,4,3,za,zb,sprod,restmp,restmpid)
69  restmp = restmpid * aveqq
70  do i = 1,5
71  res(i,-i) =res(i,-i) + restmp
72  enddo
73  call me2_qbqqbq(2,1,4,3,za,zb,sprod,restmp,restmpid)
74  restmp = restmpid * aveqq
75  do i = 1,5
76  res(-i,i) =res(-i,i) + restmp
77  enddo
78 
79  !-- qqb -> rrb
80  call me2_qbqqbq(1,2,4,3,za,zb,sprod,restmp,restmpid)
81  restmp = restmp * aveqq
82  do i = 1,5
83  res(i,-i) =res(i,-i) + restmp * (nf-1.0_dp)
84  res(-i,i) = res(-i,i) + restmp * (nf-1.0_dp)
85  enddo
86 
87  !-- qrb -> qrb
88  call me2_qbqqbq(1,3,4,2,za,zb,sprod,restmp,restmpid)
89  restmp = restmp * aveqq
90  do i = 1,5
91  do j = 1,5
92  if (i.ne.j) res(i,-j) = res(i,-j) + restmp
93  enddo
94  enddo
95  call me2_qbqqbq(2,3,4,1,za,zb,sprod,restmp,restmpid)
96  restmp = restmp * aveqq
97  do i = 1,5
98  do j = 1,5
99  if (i.ne.j) res(-j,i) = res(-j,i) + restmp
100  enddo
101  enddo
102 
103  !-- qq -> qq
104  call me2_qbqqbq(1,3,2,4,za,zb,sprod,restmp,restmpid)
105  restmp = restmpid * aveqq * symmfac
106  do i = 1,5
107  res(i,i) = res(i,i) + restmp
108  res(-i,-i) = res(-i,-i) + restmp
109  enddo
110 
111  !-- qr -> qr
112  call me2_qbqqbq(1,3,2,4,za,zb,sprod,restmp,restmpid)
113  restmp = restmp * aveqq
114  do i = 1,5
115  do j = 1,5
116  if (i.ne.j) then
117  res(i,j) = res(i,j) + restmp
118  res(-i,-j) = res(-i,-j) + restmp
119  endif
120  enddo
121  enddo
122 
123  res(:,:) = res(:,:) * (2d0/3d0*alphas**2)**2
124  return
125 
126  end subroutine evalamp_sbfh_unsymm_sa
127 
128 
129 
130 
131 
132  !-- SM: |g2| = alphas/(six*pi), so multiply ME**2 by (|g2|*gs)**2
133  !-- g3 not supported yet
134  subroutine evalamp_sbfh_unsymm_sa_select(p,iSel,jSel,flav_tag,iflip,res)
135  real(dp), intent(in) :: p(4,5)
136  integer, intent(in) :: isel,jsel,flav_tag
137  real(dp), intent(out) :: res(-5:5,-5:5)
138  real(dp) :: sprod(4,4)
139  complex(dp) :: za(4,4), zb(4,4)
140  real(dp) :: restmp, restmpid
141  integer :: i, j,iflip
142 
143  res = zero
144 
145  call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
146  iflip = 1
147 
148  !-- gg -> gg
149  if( isel.eq.pdfglu_ .and. jsel.eq.pdfglu_ .and. flav_tag.eq.1 ) then
150  call me2_ggggh(1,2,3,4,za,zb,sprod,restmp)
151  restmp = restmp * avegg * symmfac
152  restmp = restmp * (2d0/3d0*alphas**2)**2
153  res(isel,jsel) = restmp
154  return
155  endif
156 
157 
158  !-- gg -> qqb
159  if( isel.eq.pdfglu_ .and. jsel.eq.pdfglu_ .and. flav_tag.eq.2 ) then
160  call me2_qbqggh(4,3,1,2,za,zb,sprod,restmp)
161  restmp = restmp * avegg
162  restmp = restmp * (2d0/3d0*alphas**2)**2
163  res(isel,jsel) = restmp * nf
164  return
165  endif
166 
167 
168 
169  !-- gq -> gq
170  if( isel.eq.pdfglu_ .and. jsel.ne.0 ) then
171  call me2_qbqggh(2,4,1,3,za,zb,sprod,restmp)
172  restmp = restmp * aveqg
173  restmp = restmp * (2d0/3d0*alphas**2)**2
174  res(isel,jsel) = restmp
175  return
176  endif
177 
178  if( jsel.eq.pdfglu_ .and. isel.ne.0 ) then
179  iflip = 2
180  call me2_qbqggh(1,4,2,3,za,zb,sprod,restmp)
181  restmp = restmp * aveqg
182  restmp = restmp * (2d0/3d0*alphas**2)**2
183  res(isel,jsel) = restmp
184  return
185  endif
186 
187 
188 
189  !-- qqb -> gg
190  if( isel.ne.0 .and. jsel.eq.-isel .and. flav_tag.eq.1 ) then
191  call me2_qbqggh(1,2,3,4,za,zb,sprod,restmp)
192  restmp = restmp * aveqq * symmfac
193  restmp = restmp * (2d0/3d0*alphas**2)**2
194  res(isel,jsel) = restmp
195  return
196  endif
197 
198 
199  !-- qqb -> qqb
200  if( isel.gt.0 .and. jsel.eq.-isel .and. flav_tag.eq.2 ) then
201  call me2_qbqqbq(1,2,4,3,za,zb,sprod,restmp,restmpid)
202  restmp = restmpid * aveqq
203  restmp = restmp * (2d0/3d0*alphas**2)**2
204  res(isel,jsel) = restmp
205  return
206  endif
207 
208  if( isel.lt.0 .and. jsel.eq.-isel .and. flav_tag.eq.2 ) then
209  iflip = 2
210  call me2_qbqqbq(2,1,4,3,za,zb,sprod,restmp,restmpid)
211  restmp = restmpid * aveqq
212  restmp = restmp * (2d0/3d0*alphas**2)**2
213  res(isel,jsel) = restmp
214  return
215  endif
216 
217 
218  !-- qqb -> rrb
219  if( isel.ne.0 .and. jsel.eq.-isel .and. flav_tag.eq.3 ) then
220  call me2_qbqqbq(1,2,4,3,za,zb,sprod,restmp,restmpid)
221  restmp = restmp * aveqq
222  restmp = restmp * (2d0/3d0*alphas**2)**2
223  res(isel,jsel) = restmp * (nf-1.0_dp)
224  return
225  endif
226 
227 
228 
229  !-- qrb -> qrb
230  if( isel.gt.0 .and. jsel.lt.0 .and. isel.ne.jsel ) then
231  call me2_qbqqbq(1,3,4,2,za,zb,sprod,restmp,restmpid)
232  restmp = restmp * aveqq
233  restmp = restmp * (2d0/3d0*alphas**2)**2
234  res(isel,jsel) = restmp
235  return
236  endif
237 
238  if( isel.lt.0 .and. jsel.gt.0 .and. isel.ne.jsel ) then
239  iflip = 2
240  call me2_qbqqbq(2,3,4,1,za,zb,sprod,restmp,restmpid)
241  restmp = restmp * aveqq
242  restmp = restmp * (2d0/3d0*alphas**2)**2
243  res(isel,jsel) = restmp
244  return
245  endif
246 
247 
248 
249 
250  !-- qq -> qq
251  if( isel.ne.0 .and. isel.eq.jsel ) then
252  call me2_qbqqbq(1,3,2,4,za,zb,sprod,restmp,restmpid)
253  restmp = restmpid * aveqq * symmfac
254  restmp = restmp * (2d0/3d0*alphas**2)**2
255  res(isel,jsel) = restmp
256  return
257  endif
258 
259 
260 
261  !-- qr -> qr
262  if( isel.gt.0 .and. jsel.gt.0 .and. isel.ne.jsel ) then
263  call me2_qbqqbq(1,3,2,4,za,zb,sprod,restmp,restmpid)
264  restmp = restmp * aveqq
265  restmp = restmp * (2d0/3d0*alphas**2)**2
266  res(isel,jsel) = restmp
267  return
268  endif
269 
270  if( isel.lt.0 .and. jsel.lt.0 .and. isel.ne.jsel ) then
271  call me2_qbqqbq(1,3,2,4,za,zb,sprod,restmp,restmpid)
272  restmp = restmp * aveqq
273  restmp = restmp * (2d0/3d0*alphas**2)**2
274  res(isel,jsel) = restmp
275  return
276  endif
277 
278 
279 
280  return
281 
282  end subroutine
283 
284 
285  subroutine evalamp_sbfh_unsymm_sa_select_exact(p,iSel,jSel,rSel,sSel,res)
286  real(dp), intent(in) :: p(4,5)
287  integer, intent(in) :: isel,jsel,rsel,ssel!,flavor_tag ! flavor_tag for TEST
288  real(dp), intent(out) :: res
289  real(dp) :: sprod(4,4)
290  complex(dp) :: za(4,4), zb(4,4)
291  real(dp) :: restmp, restmpid
292  integer :: i, j
293  integer :: j1,j2,k1,k2
294  logical :: isgggg,isqqgg,isqqqq,isqqqq_idq
295 
296  restmp=0.0_dp
297  restmpid=0.0_dp
298  res=0.0_dp
299  isgggg=.false.
300  isqqgg=.false.
301  isqqqq=.false.
302  isqqqq_idq=.false.
303  j1 = 0
304  j2 = 0
305  k1 = 0
306  k2 = 0
307  if( (abs(isel).eq.pdftop_ .or. abs(jsel).eq.pdftop_) .or. (abs(rsel).eq.pdftop_ .or. abs(ssel).eq.pdftop_) ) return
308 
309  !print *, "Begin EvalAmp_SBFH_UnSymm_SA_Select_exact"
310  !print *, "iSel: ",iSel,", jSel: ",jSel," rSel: ",rSel," sSel: ",sSel
311 
312  if(isel.eq.pdfglu_) then ! 1==g
313  if(jsel.eq.pdfglu_) then ! gg
314  if(rsel.eq.pdfglu_ .and. ssel .eq. pdfglu_) then ! gg->gg
315  isgggg=.true.
316  j1=1
317  j2=2
318  k1=3
319  k2=4
320 
321  elseif(ssel.eq.(-rsel)) then ! gg->qqb/qbq
322  isqqgg=.true.
323  if(rsel.gt.0) then ! gg->qqb
324  j1=4
325  j2=3
326  else ! gg->qbq
327  j1=3
328  j2=4
329  endif
330  k1=1
331  k2=2
332  endif
333 
334  elseif(abs(jsel).eq.pdfup_ .or. abs(jsel).eq.pdfchm_ .or. abs(jsel).eq.pdfdn_ .or. abs(jsel).eq.pdfstr_ .or. abs(jsel).eq.pdfbot_) then ! gq/gqb
335  isqqgg=.true.
336  k1=1
337  if(jsel.gt.0) then ! gq
338  j1=2
339  if(rsel.eq.isel .and. ssel.eq.jsel) then ! gq->gq
340  k2=3
341  j2=4
342  elseif(ssel.eq.isel .and. rsel.eq.jsel) then ! gq->qg
343  j2=3
344  k2=4
345  endif
346  else ! gqb
347  j2=2
348  if(rsel.eq.isel .and. ssel.eq.jsel) then ! gqb->gqb
349  k2=3
350  j1=4
351  elseif(ssel.eq.isel .and. rsel.eq.jsel) then ! gqb->qbg
352  j1=3
353  k2=4
354  endif
355  endif
356  endif
357  elseif(abs(isel).eq.pdfup_ .or. abs(isel).eq.pdfchm_ .or. abs(isel).eq.pdfdn_ .or. abs(isel).eq.pdfstr_ .or. abs(isel).eq.pdfbot_) then ! q/qb?->?
358  if(jsel.eq.pdfglu_) then ! qg->?
359  isqqgg=.true.
360  k1=2
361  if(isel.gt.0) then
362  j1=1
363  if(rsel.eq.isel .and. ssel.eq.jsel) then ! qg->qg
364  k2=4
365  j2=3
366  elseif(ssel.eq.isel .and. rsel.eq.jsel) then ! qg->gq
367  j2=4
368  k2=3
369  endif
370  else ! qbg
371  j2=1
372  if(rsel.eq.isel .and. ssel.eq.jsel) then ! qbg->qbg
373  k2=4
374  j1=3
375  elseif(ssel.eq.isel .and. rsel.eq.jsel) then ! qbg->gqb
376  j1=4
377  k2=3
378  endif
379  endif
380 
381  elseif(abs(jsel).eq.pdfup_ .or. abs(jsel).eq.pdfchm_ .or. abs(jsel).eq.pdfdn_ .or. abs(jsel).eq.pdfstr_ .or. abs(jsel).eq.pdfbot_) then ! qq/qqb/qbq/qbqb->?
382  if(isel.gt.0) then ! qq'/qqb'->?
383  j1=1
384  else ! qbq'/qbqb'->?
385  j2=1
386  endif
387 
388  if(jsel.eq.(-isel)) then ! qqb/qbq->?
389  if(jsel.gt.0) then ! qbq->?
390  j1=2
391  else ! qqb->?
392  j2=2
393  endif
394 
395  if(rsel.eq.pdfglu_ .and. ssel .eq. pdfglu_) then ! qqb/qbq->gg
396  isqqgg=.true.
397  k1=3
398  k2=4
399  elseif(ssel.eq.(-rsel)) then ! qqb->q'qb'/qb'q'
400  isqqqq=.true.
401  if(rsel.gt.0) then ! ->q'qb'
402  if(rsel .eq. isel .or. rsel.eq.jsel) isqqqq_idq=.true.
403  k1=4
404  k2=3
405  else ! ->qb'q'
406  if(rsel .eq. isel .or. rsel.eq.jsel) isqqqq_idq=.true.
407  k1=3
408  k2=4
409  endif
410  endif
411 
412  elseif(jsel.eq.isel) then ! qq/qbqb->?
413  isqqqq=.true.
414  if(jsel.gt.0) then ! qq
415  k1=2
416  else ! qbqb
417  k2=2
418  endif
419  if(rsel.gt.0 .and. rsel.eq.ssel .and. rsel.eq.isel) then ! qq->qq
420  isqqqq_idq=.true.
421  j2=3
422  k2=4
423  elseif(rsel.lt.0 .and. rsel.eq.ssel .and. rsel.eq.isel) then ! qbqb->qbqb
424  isqqqq_idq=.true.
425  j1=3
426  k1=4
427  endif
428 
429  elseif(jsel.eq.sign(jsel,isel)) then ! qq'/qbqb'->?
430  isqqqq=.true.
431  if(jsel.gt.0) then ! qq'
432  k1=2 ! j1=1
433  else ! qbqb'
434  k2=2 ! j2=1
435  endif
436  if(rsel.gt.0 .and. rsel.ne.ssel .and. rsel.eq.sign(rsel,ssel) .and. rsel.eq.isel) then ! ->qq'
437  j2=3
438  k2=4
439  elseif(rsel.lt.0 .and. rsel.ne.ssel .and. rsel.eq.sign(rsel,ssel) .and. rsel.eq.isel) then ! ->qbqb'
440  j1=3
441  k1=4
442  elseif(rsel.gt.0 .and. rsel.ne.ssel .and. rsel.eq.sign(rsel,ssel) .and. rsel.eq.jsel) then ! ->q'q
443  j2=4
444  k2=3
445  elseif(rsel.lt.0 .and. rsel.ne.ssel .and. rsel.eq.sign(rsel,ssel) .and. rsel.eq.jsel) then ! ->qb'qb
446  j1=4
447  k1=3
448  endif
449 
450  elseif(jsel.eq.(-sign(jsel,isel))) then ! qqb'/qbq'->?
451  isqqqq=.true.
452  if(jsel.gt.0) then ! qbq'
453  k1=2 ! j2=1
454  else ! qqb'
455  k2=2 ! j1=1
456  endif
457  if(rsel.gt.0 .and. rsel.ne.ssel .and. rsel.ne.sign(rsel,ssel) .and. rsel.eq.isel) then ! qqb'->qqb'
458  j2=3
459  k1=4
460  elseif(rsel.lt.0 .and. rsel.ne.ssel .and. rsel.ne.sign(rsel,ssel) .and. rsel.eq.isel) then ! qbq'->qbq'
461  j1=3
462  k2=4
463  elseif(rsel.lt.0 .and. rsel.ne.ssel .and. rsel.ne.sign(rsel,ssel) .and. rsel.eq.jsel) then ! qqb'->qb'q
464  j2=4
465  k1=3
466  elseif(rsel.gt.0 .and. rsel.ne.ssel .and. rsel.ne.sign(rsel,ssel) .and. rsel.eq.jsel) then ! qbq'->q'qb
467  j1=4
468  k2=3
469  endif
470  endif
471 
472  endif
473  endif
474 
475  if(j1.eq.0 .or. j2.eq.0 .or. k1.eq.0 .or. k2.eq.0) then
476  print *,"Could not recognize the incoming flavors!"
477  print *,"iSel: ",isel,", jSel: ",jsel," rSel: ",rsel," sSel: ",ssel
478  print *,"j1: ",j1,", j2: ",j2," k1: ",k1," k2: ",k2
479  if(isgggg) then
480  print *,"is gggg"
481  endif
482  if(isqqgg) then
483  print *,"is qbqgg"
484  endif
485  if(isqqqq) then
486  print *,"is qbqQBQ"
487  endif
488  if(isqqqq_idq) then
489  print *,"is qbqQBQ id"
490  endif
491  endif
492 
493  call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
494 
495  if(isgggg) then
496  call me2_ggggh(j1,j2,k1,k2,za,zb,sprod,restmp)
497  restmp = restmp * avegg * symmfac
498  elseif(isqqqq) then
499  call me2_qbqqbq(j1,j2,k1,k2,za,zb,sprod,restmp,restmpid)
500  if(isqqqq_idq) then
501  restmp = restmpid * aveqq
502  if(isel.eq.jsel) restmp = restmp * symmfac
503  else
504  restmp = restmp * aveqq
505  if(isel .eq. (-jsel) .and. rsel .eq. (-ssel) .and. abs(isel) .ne. abs(rsel)) restmp = restmp * (nf-1.0_dp)
506  endif
507  elseif(isqqgg) then
508  call me2_qbqggh(j1,j2,k1,k2,za,zb,sprod,restmp)
509  if( isel.eq.pdfglu_ .and. jsel.eq.pdfglu_) then
510  restmp = restmp * avegg * nf
511  elseif( isel.ne.pdfglu_ .and. jsel.ne.pdfglu_) then
512  restmp = restmp * aveqq * symmfac
513  else
514  restmp = restmp * aveqg
515  endif
516  endif
517 
518  restmp = restmp * (2d0/3d0*alphas**2)**2
519  res = restmp
520  return
521  end subroutine
522 
523 
524 
525  subroutine evalamp_wbfh_unsymm_sa(p,res)
526  real(dp), intent(in) :: p(4,5)
527  real(dp), intent(out) :: res(-5:5,-5:5)
528  complex(dp) :: amp_z(-1:1,-1:1), amp_z_b(-1:1,-1:1)
529  complex(dp) :: amp_w(-1:1,-1:1)
530  real(dp) :: sprod(4,4)
531  complex(dp) :: za(4,4), zb(4,4)
532  real(dp) :: lu, ru
533  real(dp) :: ld, rd
534  real(dp) :: couplz
535  real(dp) :: couplw
536  real(dp) :: restmp=0d0
537  integer :: i, j, j1, j2, iflip, pdfindex(2)
538 
539  lu = al_qup**2
540  ru = ar_qup**2
541  ld = al_qdn**2
542  rd = ar_qdn**2
543  couplz = gwsq * xw/twosc**2
544  couplw = gwsq/two
545 
546  res = zero
547 
548  call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
549 
550  !-- qq->qq, up
551  amp_z = a0_vv_4f(4,1,3,2,za,zb,sprod,m_z,ga_z)
552  amp_z_b = -a0_vv_4f(3,1,4,2,za,zb,sprod,m_z,ga_z)
553 
554 ! ! adding contract terms
555 ! iprop12 = sprod(4,1) - mv**2 + ci * mv * ga_v
556 ! iprop34 = sprod(3,2) - mv**2 + ci * mv * ga_v
557 !
558 ! amp_z(-1,-1) = amp_z(-1,-1) + amp_z(-1,-1)*iprop12/(aL_QUp*couplz)*ehz_L_U &
559 ! + amp_z(-1,-1)*iprop34/(aL_QUp*couplz)*ehz_L_U &
560 ! + amp_z(-1,-1)*iprop12/(aL_QUp*couplz)*ehz_L_U &
561 ! *iprop34/(aL_QUp*couplz)*ehz_L_U
562 !
563 ! amp_z(+1,-1) = amp_z(+1,-1) + amp_z(+1,-1)*iprop12/(aR_QUp*couplz)*ehz_R_U &
564 ! + amp_z(+1,-1)*iprop34/(aL_QUp*couplz)*ehz_L_U &
565 ! + amp_z(+1,-1)*iprop12/(aR_QUp*couplz)*ehz_R_U &
566 ! *iprop34/(aL_QUp*couplz)*ehz_L_U
567 !
568 ! amp_z(-1,+1) = amp_z(-1,+1) + amp_z(-1,+1)*iprop12/(aL_QUp*couplz)*ehz_L_U &
569 ! + amp_z(-1,+1)*iprop34/(aR_QUp*couplz)*ehz_R_U &
570 ! + amp_z(-1,+1)*iprop12/(aL_QUp*couplz)*ehz_L_U &
571 ! *iprop34/(aR_QUp*couplz)*ehz_R_U
572 !
573 ! amp_z(+1,+1) = amp_z(+1,+1) + amp_z(+1,+1)*iprop12/(aR_QUp*couplz)*ehz_R_U &
574 ! + amp_z(+1,+1)*iprop34/(aR_QUp*couplz)*ehz_R_U &
575 ! + amp_z(+1,+1)*iprop12/(aR_QUp*couplz)*ehz_R_U &
576 ! *iprop34/(aR_QUp*couplz)*ehz_R_U
577 
578 
579 ! Lu = aL_QUp + iprop12/couplz*ehz_L_U
580 ! Ru = aR_QUp + iprop12/couplz*ehz_R_U
581 ! amp_z(-1,-1) = amp_z(-1,-1)*Lu*Lu
582 ! amp_z(-1,+1) = amp_z(-1,+1)*Lu*Ru
583 ! amp_z(+1,-1) = amp_z(+1,-1)*Ru*Lu
584 ! amp_z(+1,+1) = amp_z(+1,+1)*Ru*Ru
585 !
586 ! amp_z_b(-1,-1) = amp_z_b(-1,-1)*Lu*Lu
587 ! amp_z_b(-1,+1) = amp_z_b(-1,+1)*Lu*Ru
588 ! amp_z_b(+1,-1) = amp_z_b(+1,-1)*Ru*Lu
589 ! amp_z_b(+1,+1) = amp_z_b(+1,+1)*Ru*Ru
590 
591  restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) * lu**2 + &
592  (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) * lu * ru + &
593  (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) * lu * ru + &
594  (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2) * ru**2) * xn**2
595 
596  restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) * lu**2 + &
597  two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp) * ru**2) * xn
598 
599  restmp = restmp * symmfac * aveqq * couplz**2
600 
601  res(pdfup_,pdfup_) = restmp
602  res(pdfchm_,pdfchm_) = restmp
603 
604 
605 
606 
607 
608  !-- qq->qq, down
609  restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) * ld**2 + &
610  (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) * ld * rd + &
611  (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) * ld * rd + &
612  (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2) * rd**2) * xn**2
613 
614  restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) * ld**2 + &
615  two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp) * rd**2) * xn
616 
617  restmp = restmp * symmfac * aveqq * couplz**2
618 
619  res(pdfdn_,pdfdn_) = restmp
620  res(pdfstr_,pdfstr_) = restmp
621  res(pdfbot_,pdfbot_) = restmp * tagbot
622 
623  !-- qbqb->qbqb, aup
624  amp_z = a0_vv_4f(1,4,2,3,za,zb,sprod,m_z,ga_z)
625  amp_z_b = -a0_vv_4f(1,3,2,4,za,zb,sprod,m_z,ga_z)
626 
627  restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) * lu**2 + &
628  (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) * lu * ru + &
629  (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) * lu * ru + &
630  (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2) * ru**2) * xn**2
631 
632  restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) * lu**2 + &
633  two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp) * ru**2) * xn
634 
635  restmp = restmp * symmfac * aveqq * couplz**2
636 
637  res(pdfaup_,pdfaup_) = restmp
638  res(pdfachm_,pdfachm_) = restmp
639 
640  !-- qbqb->qbqb, adn
641  restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) * ld**2 + &
642  (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) * ld * rd + &
643  (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) * ld * rd + &
644  (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2) * rd**2) * xn**2
645 
646  restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) * ld**2 + &
647  two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp) * rd**2) * xn
648 
649  restmp = restmp * symmfac * aveqq * couplz**2
650 
651  res(pdfadn_,pdfadn_) = restmp
652  res(pdfastr_,pdfastr_) = restmp
653  res(pdfabot_,pdfabot_) = restmp * tagbot
654 
655  !-- W/Z interference
656  j1 = 1
657  j2 = 2
658  do iflip = 1, 2
659  !-- ud -> ud
660  amp_z = a0_vv_4f(4,j2,3,j1,za,zb,sprod,m_z,ga_z)
661  amp_w = -a0_vv_4f(4,j1,3,j2,za,zb,sprod,m_w,ga_w,usewwcoupl=.true.)
662 
663  restmp = ((abs(amp_z(-1,-1))**2) * ld * lu + &
664  (abs(amp_z(-1,+1))**2) * ld * ru + &
665  (abs(amp_z(+1,-1))**2) * rd * lu + &
666  (abs(amp_z(+1,+1))**2) * rd * ru) * couplz**2 * xn**2
667 
668  restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * xn**2
669 
670  restmp = restmp + two * real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) * &
671  al_qup * al_qdn * couplz * couplw * xn
672 
673  restmp = restmp * aveqq
674 
675  pdfindex = flip(iflip,pdfup_,pdfdn_)
676  res(pdfindex(1),pdfindex(2)) = restmp
677  pdfindex = flip(iflip,pdfchm_,pdfstr_)
678  res(pdfindex(1),pdfindex(2)) = restmp
679 
680  !-- ubdb -> ubdb
681  amp_z = a0_vv_4f(j2,4,j1,3,za,zb,sprod,m_z,ga_z)
682  amp_w = -a0_vv_4f(j1,4,j2,3,za,zb,sprod,m_w,ga_w,usewwcoupl=.true.)
683 
684  restmp = ((abs(amp_z(-1,-1))**2) * ld * lu + &
685  (abs(amp_z(-1,+1))**2) * ld * ru + &
686  (abs(amp_z(+1,-1))**2) * rd * lu + &
687  (abs(amp_z(+1,+1))**2) * rd * ru) * couplz**2 * xn**2
688 
689  restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * xn**2
690 
691  restmp = restmp + two * real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) * &
692  al_qup * al_qdn * couplz * couplw * xn
693 
694  restmp = restmp * aveqq
695 
696  pdfindex = flip(iflip,pdfaup_,pdfadn_)
697  res(pdfindex(1),pdfindex(2)) = restmp
698  pdfindex = flip(iflip,pdfachm_,pdfastr_)
699  res(pdfindex(1),pdfindex(2)) = restmp
700 
701  j1 = 2
702  j2 = 1
703  enddo
704 
705  !-- non-interfering diagrams below
706  j1 = 1
707  j2 = 2
708  do iflip = 1,2
709 
710  !-- qqb processes
711 
712  !--uub -> uub // ddb
713  amp_z = a0_vv_4f(3,j1,j2,4,za,zb,sprod,m_z,ga_z)
714  amp_w = a0_vv_4f(3,j1,j2,4,za,zb,sprod,m_w,ga_w,usewwcoupl=.true.)
715 
716  restmp = ((abs(amp_z(-1,-1))**2) * lu * lu + &
717  (abs(amp_z(-1,+1))**2) * lu * ru + &
718  (abs(amp_z(+1,-1))**2) * ru * lu + &
719  (abs(amp_z(+1,+1))**2) * ru * ru) * couplz**2 * spinavg * tag1
720  restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * spinavg * tag2
721 
722  pdfindex = flip(iflip,pdfup_,pdfaup_)
723  res(pdfindex(1),pdfindex(2)) = restmp
724  pdfindex = flip(iflip,pdfchm_,pdfachm_)
725  res(pdfindex(1),pdfindex(2)) = restmp
726  pdfindex = flip(iflip,pdfup_,pdfachm_)
727  res(pdfindex(1),pdfindex(2)) = restmp
728  pdfindex = flip(iflip,pdfchm_,pdfaup_)
729  res(pdfindex(1),pdfindex(2)) = restmp
730 
731  !--udb -> udb
732  restmp = ((abs(amp_z(-1,-1))**2) * lu * ld + &
733  (abs(amp_z(-1,+1))**2) * lu * rd + &
734  (abs(amp_z(+1,-1))**2) * ru * ld + &
735  (abs(amp_z(+1,+1))**2) * ru * rd) * couplz**2 * spinavg
736 
737  pdfindex = flip(iflip,pdfup_,pdfadn_)
738  res(pdfindex(1),pdfindex(2)) = restmp
739  pdfindex = flip(iflip,pdfchm_,pdfastr_)
740  res(pdfindex(1),pdfindex(2)) = restmp
741  pdfindex = flip(iflip,pdfup_,pdfastr_)
742  res(pdfindex(1),pdfindex(2)) = restmp
743  pdfindex = flip(iflip,pdfchm_,pdfadn_)
744  res(pdfindex(1),pdfindex(2)) = restmp
745  pdfindex = flip(iflip,pdfup_,pdfabot_)
746  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
747  pdfindex = flip(iflip,pdfchm_,pdfabot_)
748  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
749 
750  !--dub -> dub
751  restmp = ((abs(amp_z(-1,-1))**2) * ld * lu + &
752  (abs(amp_z(-1,+1))**2) * ld * ru + &
753  (abs(amp_z(+1,-1))**2) * rd * lu + &
754  (abs(amp_z(+1,+1))**2) * rd * ru) * couplz**2 * spinavg
755 
756  pdfindex = flip(iflip,pdfdn_,pdfaup_)
757  res(pdfindex(1),pdfindex(2)) = restmp
758  pdfindex = flip(iflip,pdfstr_,pdfachm_)
759  res(pdfindex(1),pdfindex(2)) = restmp
760  pdfindex = flip(iflip,pdfdn_,pdfachm_)
761  res(pdfindex(1),pdfindex(2)) = restmp
762  pdfindex = flip(iflip,pdfstr_,pdfaup_)
763  res(pdfindex(1),pdfindex(2)) = restmp
764  pdfindex = flip(iflip,pdfbot_,pdfaup_)
765  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
766  pdfindex = flip(iflip,pdfbot_,pdfachm_)
767  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
768 
769  !--ddb -> uub/ddb
770  restmp = ((abs(amp_z(-1,-1))**2) * ld * ld + &
771  (abs(amp_z(-1,+1))**2) * ld * rd + &
772  (abs(amp_z(+1,-1))**2) * rd * ld + &
773  (abs(amp_z(+1,+1))**2) * rd * rd) * couplz**2 * spinavg * tag1
774 
775  pdfindex = flip(iflip,pdfbot_,pdfabot_)
776  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
777  pdfindex = flip(iflip,pdfdn_,pdfabot_)
778  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
779  pdfindex = flip(iflip,pdfstr_,pdfabot_)
780  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
781  pdfindex = flip(iflip,pdfbot_,pdfadn_)
782  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
783  pdfindex = flip(iflip,pdfbot_,pdfastr_)
784  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
785 
786  restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * spinavg * tag2
787 
788  pdfindex = flip(iflip,pdfdn_,pdfadn_)
789  res(pdfindex(1),pdfindex(2)) = restmp
790  pdfindex = flip(iflip,pdfstr_,pdfastr_)
791  res(pdfindex(1),pdfindex(2)) = restmp
792  pdfindex = flip(iflip,pdfdn_,pdfastr_)
793  res(pdfindex(1),pdfindex(2)) = restmp
794  pdfindex = flip(iflip,pdfstr_,pdfadn_)
795  res(pdfindex(1),pdfindex(2)) = restmp
796 
797  !-- non-symmetric qq processes
798  amp_z = a0_vv_4f(3,j1,4,j2,za,zb,sprod,m_z,ga_z)
799  amp_w = a0_vv_4f(3,j2,4,j1,za,zb,sprod,m_w,ga_w,usewwcoupl=.true.)
800 
801  !--uc -> uc
802  restmp = ((abs(amp_z(-1,-1))**2) * lu * lu + &
803  (abs(amp_z(-1,+1))**2) * lu * ru + &
804  (abs(amp_z(+1,-1))**2) * ru * lu + &
805  (abs(amp_z(+1,+1))**2) * ru * ru) * couplz**2 * spinavg
806 
807  pdfindex = flip(iflip,pdfup_,pdfchm_)
808  res(pdfindex(1),pdfindex(2)) = restmp
809 
810  !--us -> us/cd
811  restmp = ((abs(amp_z(-1,-1))**2) * lu * ld + &
812  (abs(amp_z(-1,+1))**2) * lu * rd + &
813  (abs(amp_z(+1,-1))**2) * ru * ld + &
814  (abs(amp_z(+1,+1))**2) * ru * rd) * couplz**2 * spinavg * tag1
815 
816  pdfindex = flip(iflip,pdfup_,pdfbot_)
817  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
818  pdfindex = flip(iflip,pdfchm_,pdfbot_)
819  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
820 
821  restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * spinavg * tag2
822 
823  pdfindex = flip(iflip,pdfup_,pdfstr_)
824  res(pdfindex(1),pdfindex(2)) = restmp
825  pdfindex = flip(iflip,pdfchm_,pdfdn_)
826  res(pdfindex(1),pdfindex(2)) = restmp
827 
828  !--ds -> ds
829  restmp = ((abs(amp_z(-1,-1))**2) * ld * ld + &
830  (abs(amp_z(-1,+1))**2) * ld * rd + &
831  (abs(amp_z(+1,-1))**2) * rd * ld + &
832  (abs(amp_z(+1,+1))**2) * rd * rd) * couplz**2 * spinavg
833 
834  pdfindex = flip(iflip,pdfdn_,pdfstr_)
835  res(pdfindex(1),pdfindex(2)) = restmp
836  pdfindex = flip(iflip,pdfdn_,pdfbot_)
837  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
838  pdfindex = flip(iflip,pdfstr_,pdfbot_)
839  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
840 
841  !-- qbqb processes
842  amp_z = a0_vv_4f(j1,3,j2,4,za,zb,sprod,m_z,ga_z)
843  amp_w = a0_vv_4f(j1,4,j2,3,za,zb,sprod,m_w,ga_w,usewwcoupl=.true.)
844 
845  !--ubcb -> ubcb
846  restmp = ((abs(amp_z(-1,-1))**2) * lu * lu + &
847  (abs(amp_z(-1,+1))**2) * lu * ru + &
848  (abs(amp_z(+1,-1))**2) * ru * lu + &
849  (abs(amp_z(+1,+1))**2) * ru * ru) * couplz**2 * spinavg
850 
851  pdfindex = flip(iflip,pdfaup_,pdfachm_)
852  res(pdfindex(1),pdfindex(2)) = restmp
853 
854  !--ubsb -> ubsb//cbdb
855  restmp = ((abs(amp_z(-1,-1))**2) * lu * ld + &
856  (abs(amp_z(-1,+1))**2) * lu * rd + &
857  (abs(amp_z(+1,-1))**2) * ru * ld + &
858  (abs(amp_z(+1,+1))**2) * ru * rd) * couplz**2 * spinavg * tag1
859 
860  pdfindex = flip(iflip,pdfaup_,pdfabot_)
861  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
862  pdfindex = flip(iflip,pdfachm_,pdfabot_)
863  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
864 
865  restmp = restmp + abs(amp_w(-1,-1))**2 * couplw**2 * spinavg * tag2
866 
867  pdfindex = flip(iflip,pdfaup_,pdfastr_)
868  res(pdfindex(1),pdfindex(2)) = restmp
869  pdfindex = flip(iflip,pdfachm_,pdfadn_)
870  res(pdfindex(1),pdfindex(2)) = restmp
871 
872  !--dbsb -> dbsb
873  restmp = ((abs(amp_z(-1,-1))**2) * ld * ld + &
874  (abs(amp_z(-1,+1))**2) * ld * rd + &
875  (abs(amp_z(+1,-1))**2) * rd * ld + &
876  (abs(amp_z(+1,+1))**2) * rd * rd) * couplz**2 * spinavg
877 
878  pdfindex = flip(iflip,pdfadn_,pdfastr_)
879  res(pdfindex(1),pdfindex(2)) = restmp
880  pdfindex = flip(iflip,pdfadn_,pdfabot_)
881  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
882  pdfindex = flip(iflip,pdfastr_,pdfabot_)
883  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
884 
885  j1 = 2
886  j2 = 1
887 
888  enddo
889 
890  return
891 
892  end subroutine evalamp_wbfh_unsymm_sa
893 
894 
895 
896 
897  SUBROUTINE evalamp_wbfh_unsymm_sa_select(p,iSel,jSel,zz_fusion,iflip,res)
898  implicit none
899  real(dp), intent(in) :: p(4,5)
900  real(dp), intent(out) :: res(-5:5,-5:5)
901  logical, intent(in) :: zz_fusion
902  integer, intent(in) :: isel,jsel
903  complex(dp) :: amp_z(-1:1,-1:1), amp_z_b(-1:1,-1:1)
904  complex(dp) :: amp_w(-1:1,-1:1)
905  real(dp) :: sprod(4,4)
906  complex(dp) :: za(4,4), zb(4,4)
907  real(dp) :: restmp
908  integer :: i, j, j1, j2, iflip, pdfindex(2)
909 
910  call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
911  iflip = 1
912 
913  !-- qq->qq, up
914  if( (isel.eq.pdfup_ .and. jsel.eq.pdfup_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfchm_) ) then
915 
916  amp_z = a0_zz_4f(4,1,3,2,za,zb,sprod,isel,jsel)
917  amp_z_b = -a0_zz_4f(3,1,4,2,za,zb,sprod,isel,jsel)
918 
919  restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) + &
920  (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) + &
921  (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) + &
922  (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2)) * xn**2
923 
924  restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) + &
925  two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp)) * xn
926 
927  restmp = restmp * symmfac * aveqq
928 
929  if( .not. zz_fusion ) restmp = 0.0d0
930 
931  res(pdfup_,pdfup_) = restmp
932  res(pdfchm_,pdfchm_) = restmp
933 
934  return
935  endif
936 
937 
938  !-- qq->qq, down
939  if( (isel.eq.pdfdn_ .and. jsel.eq.pdfdn_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfstr_) .or. (isel.eq.pdfbot_ .and. jsel.eq.pdfbot_) ) then
940  amp_z = a0_zz_4f(4,1,3,2,za,zb,sprod,isel,jsel)
941  amp_z_b = -a0_zz_4f(3,1,4,2,za,zb,sprod,isel,jsel)
942 
943  restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) + &
944  (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) + &
945  (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) + &
946  (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2)) * xn**2
947 
948  restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) + &
949  two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp)) * xn
950 
951  restmp = restmp * symmfac * aveqq
952 
953  if( .not. zz_fusion ) restmp = 0.0d0
954 
955  res(pdfdn_,pdfdn_) = restmp
956  res(pdfstr_,pdfstr_) = restmp
957  res(pdfbot_,pdfbot_) = restmp * tagbot
958  return
959  endif
960 
961 
962  !-- qbqb->qbqb, aup
963  if( (isel.eq.pdfaup_ .and. jsel.eq.pdfaup_) .or. (isel.eq.pdfachm_ .and. jsel.eq.pdfachm_) ) then
964  amp_z = a0_zz_4f(1,4,2,3,za,zb,sprod,isel,jsel)
965  amp_z_b = -a0_zz_4f(1,3,2,4,za,zb,sprod,isel,jsel)
966 
967  restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) + &
968  (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) + &
969  (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) + &
970  (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2)) * xn**2
971 
972  restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) + &
973  two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp)) * xn
974 
975  restmp = restmp * symmfac * aveqq
976 
977  if( .not. zz_fusion ) restmp = 0.0d0
978 
979  res(pdfaup_,pdfaup_) = restmp
980  res(pdfachm_,pdfachm_) = restmp
981  return
982  endif
983 
984 
985  !-- qbqb->qbqb, adn
986  if( (isel.eq.pdfadn_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfastr_ .and. jsel.eq.pdfastr_) .or. (isel.eq.pdfabot_ .and. jsel.eq.pdfabot_) ) then
987  amp_z = a0_zz_4f(1,4,2,3,za,zb,sprod,isel,jsel)
988  amp_z_b = -a0_zz_4f(1,3,2,4,za,zb,sprod,isel,jsel)
989  restmp = ((abs(amp_z(-1,-1))**2+abs(amp_z_b(-1,-1))**2) + &
990  (abs(amp_z(-1,+1))**2+abs(amp_z_b(-1,+1))**2) + &
991  (abs(amp_z(+1,-1))**2+abs(amp_z_b(+1,-1))**2) + &
992  (abs(amp_z(+1,+1))**2+abs(amp_z_b(+1,+1))**2)) * xn**2
993 
994  restmp = restmp + (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) + &
995  two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp)) * xn
996 
997  restmp = restmp * symmfac * aveqq
998 
999  if( .not. zz_fusion ) restmp = 0.0d0
1000 
1001  res(pdfadn_,pdfadn_) = restmp
1002  res(pdfastr_,pdfastr_) = restmp
1003  res(pdfabot_,pdfabot_) = restmp * tagbot
1004  return
1005  endif
1006 
1007 
1008 
1009  ! NEED TO ADJSUT J1,J2 in W AMPS HERE
1010  !-- W/Z interference
1011  j1 = 1
1012  j2 = 2
1013  iflip = 1
1014 
1015  !-- ud -> ud
1016  if( (isel.eq.pdfup_ .and. jsel.eq.pdfdn_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfstr_) ) then
1017  amp_z = a0_zz_4f(4,j2,3,j1,za,zb,sprod,jsel,isel)
1018  amp_w = -a0_ww_4f(4,j1,3,j2,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.true.)
1019  if( zz_fusion ) then
1020  restmp = ((abs(amp_z(-1,-1))**2) + &
1021  (abs(amp_z(-1,+1))**2) + &
1022  (abs(amp_z(+1,-1))**2) + &
1023  (abs(amp_z(+1,+1))**2)) * xn**2
1024  else
1025  restmp= ( &
1026  abs(amp_w(-1,-1))**2 + &
1027  abs(amp_w(-1,+1))**2 + &
1028  abs(amp_w(+1,-1))**2 + &
1029  abs(amp_w(+1,+1))**2 &
1030  ) * xn**2
1031  restmp = restmp + two * ( &
1032  real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) + &
1033  real(amp_z(-1,+1)*conjg(amp_w(-1,+1)),kind=dp) + &
1034  real(amp_z(+1,-1)*conjg(amp_w(+1,-1)),kind=dp) + &
1035  real(amp_z(+1,+1)*conjg(amp_w(+1,+1)),kind=dp) &
1036  ) * xn
1037 
1038  endif
1039 
1040  restmp = restmp * aveqq
1041 
1042  pdfindex = flip(iflip,pdfup_,pdfdn_)
1043  res(pdfindex(1),pdfindex(2)) = restmp
1044 
1045  pdfindex = flip(iflip,pdfchm_,pdfstr_)
1046  res(pdfindex(1),pdfindex(2)) = restmp
1047 
1048  return
1049  endif
1050 
1051 
1052  !-- ubdb -> ubdb
1053  if( (isel.eq.pdfaup_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfachm_ .and. jsel.eq.pdfastr_) ) then
1054  amp_z = a0_zz_4f(j2,4,j1,3,za,zb,sprod,jsel,isel)
1055  amp_w = -a0_ww_4f(j1,4,j2,3,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.false.)
1056 
1057  if( zz_fusion ) then
1058  restmp = ((abs(amp_z(-1,-1))**2) + &
1059  (abs(amp_z(-1,+1))**2) + &
1060  (abs(amp_z(+1,-1))**2) + &
1061  (abs(amp_z(+1,+1))**2)) * xn**2
1062  else
1063  restmp= ( &
1064  abs(amp_w(-1,-1))**2 + &
1065  abs(amp_w(-1,+1))**2 + &
1066  abs(amp_w(+1,-1))**2 + &
1067  abs(amp_w(+1,+1))**2 &
1068  ) * xn**2
1069  restmp = restmp + two * ( &
1070  real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) + &
1071  real(amp_z(-1,+1)*conjg(amp_w(-1,+1)),kind=dp) + &
1072  real(amp_z(+1,-1)*conjg(amp_w(+1,-1)),kind=dp) + &
1073  real(amp_z(+1,+1)*conjg(amp_w(+1,+1)),kind=dp) &
1074  ) * xn
1075 
1076  endif
1077 
1078  restmp = restmp * aveqq
1079 
1080  pdfindex = flip(iflip,pdfaup_,pdfadn_)
1081  res(pdfindex(1),pdfindex(2)) = restmp
1082  pdfindex = flip(iflip,pdfachm_,pdfastr_)
1083  res(pdfindex(1),pdfindex(2)) = restmp
1084 
1085  return
1086  endif
1087 
1088 
1089  ! NEED TO ADJSUT J1,J2 in W AMPS HERE
1090  !-- W/Z interference
1091  j1 = 2
1092  j2 = 1
1093  iflip = 2
1094 
1095  !-- ud -> ud
1096  if( (isel.eq.pdfdn_ .and. jsel.eq.pdfup_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfchm_) ) then
1097  amp_z = a0_zz_4f(4,j2,3,j1,za,zb,sprod,isel,jsel)
1098  amp_w = -a0_ww_4f(4,j1,3,j2,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.true.)
1099  if( zz_fusion ) then
1100  restmp = ((abs(amp_z(-1,-1))**2) + &
1101  (abs(amp_z(-1,+1))**2) + &
1102  (abs(amp_z(+1,-1))**2) + &
1103  (abs(amp_z(+1,+1))**2)) * xn**2
1104  else
1105  restmp= ( &
1106  abs(amp_w(-1,-1))**2 + &
1107  abs(amp_w(-1,+1))**2 + &
1108  abs(amp_w(+1,-1))**2 + &
1109  abs(amp_w(+1,+1))**2 &
1110  ) * xn**2
1111  restmp = restmp + two * ( &
1112  real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) + &
1113  real(amp_z(-1,+1)*conjg(amp_w(-1,+1)),kind=dp) + &
1114  real(amp_z(+1,-1)*conjg(amp_w(+1,-1)),kind=dp) + &
1115  real(amp_z(+1,+1)*conjg(amp_w(+1,+1)),kind=dp) &
1116  ) * xn
1117 
1118  endif
1119 
1120  restmp = restmp * aveqq
1121 
1122  pdfindex = flip(iflip,pdfup_,pdfdn_)
1123  res(pdfindex(1),pdfindex(2)) = restmp
1124 
1125  pdfindex = flip(iflip,pdfchm_,pdfstr_)
1126  res(pdfindex(1),pdfindex(2)) = restmp
1127 
1128  return
1129  endif
1130 
1131 
1132  !-- ubdb -> ubdb
1133  if( (isel.eq.pdfadn_ .and. jsel.eq.pdfaup_) .or. (isel.eq.pdfastr_ .and. jsel.eq.pdfachm_) ) then
1134  amp_z = a0_zz_4f(j2,4,j1,3,za,zb,sprod,isel,jsel)
1135  amp_w = -a0_ww_4f(j1,4,j2,3,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.false.)
1136 
1137  if( zz_fusion ) then
1138  restmp = ((abs(amp_z(-1,-1))**2) + &
1139  (abs(amp_z(-1,+1))**2) + &
1140  (abs(amp_z(+1,-1))**2) + &
1141  (abs(amp_z(+1,+1))**2)) * xn**2
1142  else
1143  restmp= ( &
1144  abs(amp_w(-1,-1))**2 + &
1145  abs(amp_w(-1,+1))**2 + &
1146  abs(amp_w(+1,-1))**2 + &
1147  abs(amp_w(+1,+1))**2 &
1148  ) * xn**2
1149  restmp = restmp + two * ( &
1150  real(amp_z(-1,-1)*conjg(amp_w(-1,-1)),kind=dp) + &
1151  real(amp_z(-1,+1)*conjg(amp_w(-1,+1)),kind=dp) + &
1152  real(amp_z(+1,-1)*conjg(amp_w(+1,-1)),kind=dp) + &
1153  real(amp_z(+1,+1)*conjg(amp_w(+1,+1)),kind=dp) &
1154  ) * xn
1155  endif
1156 
1157  restmp = restmp * aveqq
1158 
1159  pdfindex = flip(iflip,pdfaup_,pdfadn_)
1160  res(pdfindex(1),pdfindex(2)) = restmp
1161 
1162  pdfindex = flip(iflip,pdfachm_,pdfastr_)
1163  res(pdfindex(1),pdfindex(2)) = restmp
1164 
1165  return
1166  endif
1167 
1168 
1169  !-- non-interfering diagrams below
1170  j1 = 1
1171  j2 = 2
1172  iflip = 1
1173 
1174  !-- qqb processes
1175 
1176  !--uub -> uub // ddb
1177  if( (isel.eq.pdfup_ .and. jsel.eq.pdfaup_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfachm_) .or. (isel.eq.pdfup_ .and. jsel.eq.pdfachm_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfaup_) ) then
1178 
1179  if( zz_fusion ) then
1180  amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,isel,jsel)
1181  restmp = ((abs(amp_z(-1,-1))**2) + &
1182  (abs(amp_z(-1,+1))**2) + &
1183  (abs(amp_z(+1,-1))**2) + &
1184  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1185  else
1186  amp_w = a0_ww_4f(4,j1,j2,3,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.true.)
1187  restmp= ( &
1188  abs(amp_w(-1,-1))**2 + &
1189  abs(amp_w(-1,+1))**2 + &
1190  abs(amp_w(+1,-1))**2 + &
1191  abs(amp_w(+1,+1))**2 &
1192  ) * spinavg * tag2
1193  endif
1194 
1195  pdfindex = flip(iflip,pdfup_,pdfaup_)
1196  res(pdfindex(1),pdfindex(2)) = restmp
1197 
1198  pdfindex = flip(iflip,pdfchm_,pdfachm_)
1199  res(pdfindex(1),pdfindex(2)) = restmp
1200 
1201  pdfindex = flip(iflip,pdfup_,pdfachm_)
1202  res(pdfindex(1),pdfindex(2)) = restmp
1203 
1204  pdfindex = flip(iflip,pdfchm_,pdfaup_)
1205  res(pdfindex(1),pdfindex(2)) = restmp
1206 
1207  return
1208  endif
1209 
1210 
1211  !--udb -> udb
1212  if( (isel.eq.pdfup_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfastr_) .or. (isel.eq.pdfup_ .and. jsel.eq.pdfastr_) .or. &
1213  (isel.eq.pdfchm_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfup_ .and. jsel.eq.pdfabot_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfabot_) ) then
1214  amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,isel,jsel)
1215  restmp = ((abs(amp_z(-1,-1))**2) + &
1216  (abs(amp_z(-1,+1))**2) + &
1217  (abs(amp_z(+1,-1))**2) + &
1218  (abs(amp_z(+1,+1))**2)) * spinavg
1219 
1220  if( .not. zz_fusion ) restmp = 0.0d0
1221 
1222  pdfindex = flip(iflip,pdfup_,pdfadn_)
1223  res(pdfindex(1),pdfindex(2)) = restmp
1224 
1225  pdfindex = flip(iflip,pdfchm_,pdfastr_)
1226  res(pdfindex(1),pdfindex(2)) = restmp
1227 
1228  pdfindex = flip(iflip,pdfup_,pdfastr_)
1229  res(pdfindex(1),pdfindex(2)) = restmp
1230 
1231  pdfindex = flip(iflip,pdfchm_,pdfadn_)
1232  res(pdfindex(1),pdfindex(2)) = restmp
1233 
1234  pdfindex = flip(iflip,pdfup_,pdfabot_)
1235  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1236 
1237  pdfindex = flip(iflip,pdfchm_,pdfabot_)
1238  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1239  return
1240  endif
1241 
1242 
1243  !--dub -> dub
1244  if( (isel.eq.pdfdn_ .and. jsel.eq.pdfaup_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfachm_) .or. (isel.eq.pdfdn_ .and. jsel.eq.pdfachm_) .or. &
1245  (isel.eq.pdfstr_ .and. jsel.eq.pdfaup_) .or. (isel.eq.pdfbot_ .and. jsel.eq.pdfaup_) .or. (isel.eq.pdfbot_ .and. jsel.eq.pdfachm_) ) then
1246  amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,isel,jsel)
1247  restmp = ((abs(amp_z(-1,-1))**2) + &
1248  (abs(amp_z(-1,+1))**2) + &
1249  (abs(amp_z(+1,-1))**2) + &
1250  (abs(amp_z(+1,+1))**2)) * spinavg
1251 
1252  if( .not. zz_fusion ) restmp = 0.0d0
1253 
1254  pdfindex = flip(iflip,pdfdn_,pdfaup_)
1255  res(pdfindex(1),pdfindex(2)) = restmp
1256 
1257  pdfindex = flip(iflip,pdfstr_,pdfachm_)
1258  res(pdfindex(1),pdfindex(2)) = restmp
1259 
1260  pdfindex = flip(iflip,pdfdn_,pdfachm_)
1261  res(pdfindex(1),pdfindex(2)) = restmp
1262 
1263  pdfindex = flip(iflip,pdfstr_,pdfaup_)
1264  res(pdfindex(1),pdfindex(2)) = restmp
1265 
1266  pdfindex = flip(iflip,pdfbot_,pdfaup_)
1267  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1268 
1269  pdfindex = flip(iflip,pdfbot_,pdfachm_)
1270  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1271  return
1272  endif
1273 
1274 
1275  !--ddb -> uub/ddb
1276  if( (isel.eq.pdfbot_ .and. jsel.eq.pdfabot_) .or. (isel.eq.pdfdn_ .and. jsel.eq.pdfabot_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfabot_) .or.&
1277  (isel.eq.pdfbot_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfbot_ .and. jsel.eq.pdfastr_) ) then
1278  amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,isel,jsel)
1279 
1280  restmp = ((abs(amp_z(-1,-1))**2) + &
1281  (abs(amp_z(-1,+1))**2) + &
1282  (abs(amp_z(+1,-1))**2) + &
1283  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1284 
1285  if( .not. zz_fusion ) restmp = 0.0d0
1286 
1287  pdfindex = flip(iflip,pdfbot_,pdfabot_)
1288  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1289 
1290  pdfindex = flip(iflip,pdfdn_,pdfabot_)
1291  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1292 
1293  pdfindex = flip(iflip,pdfstr_,pdfabot_)
1294  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1295 
1296  pdfindex = flip(iflip,pdfbot_,pdfadn_)
1297  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1298 
1299  pdfindex = flip(iflip,pdfbot_,pdfastr_)
1300  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1301  return
1302  endif
1303 
1304 
1305  if( (isel.eq.pdfdn_ .and. jsel.eq.pdfadn_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfastr_) .or. (isel.eq.pdfdn_ .and. jsel.eq.pdfastr_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfadn_) ) then
1306  if( zz_fusion ) then
1307  amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,isel,jsel)
1308  restmp = ((abs(amp_z(-1,-1))**2) + &
1309  (abs(amp_z(-1,+1))**2) + &
1310  (abs(amp_z(+1,-1))**2) + &
1311  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1312  else
1313  amp_w = a0_ww_4f(4,j1,j2,3,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.false.)
1314  restmp= ( &
1315  abs(amp_w(-1,-1))**2 + &
1316  abs(amp_w(-1,+1))**2 + &
1317  abs(amp_w(+1,-1))**2 + &
1318  abs(amp_w(+1,+1))**2 &
1319  ) * spinavg * tag2
1320  endif
1321 
1322  pdfindex = flip(iflip,pdfdn_,pdfadn_)
1323  res(pdfindex(1),pdfindex(2)) = restmp
1324 
1325  pdfindex = flip(iflip,pdfstr_,pdfastr_)
1326  res(pdfindex(1),pdfindex(2)) = restmp
1327 
1328  pdfindex = flip(iflip,pdfdn_,pdfastr_)
1329  res(pdfindex(1),pdfindex(2)) = restmp
1330 
1331  pdfindex = flip(iflip,pdfstr_,pdfadn_)
1332  res(pdfindex(1),pdfindex(2)) = restmp
1333 
1334  return
1335  endif
1336 
1337 
1338  if( (isel.eq.pdfup_ .and. jsel.eq.pdfchm_) ) then
1339 
1340  !-- non-symmetric qq processes
1341  amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,isel,jsel)
1342  !--uc -> uc
1343  restmp = ((abs(amp_z(-1,-1))**2) + &
1344  (abs(amp_z(-1,+1))**2) + &
1345  (abs(amp_z(+1,-1))**2) + &
1346  (abs(amp_z(+1,+1))**2)) * spinavg
1347 
1348  if( .not. zz_fusion ) restmp = 0.0d0
1349 
1350  pdfindex = flip(iflip,pdfup_,pdfchm_)
1351  res(pdfindex(1),pdfindex(2)) = restmp
1352 
1353  return
1354  endif
1355 
1356 
1357  !--us -> us/cd
1358  if( (isel.eq.pdfup_ .and. jsel.eq.pdfbot_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfbot_) ) then
1359 
1360  amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,isel,jsel)
1361  restmp = ((abs(amp_z(-1,-1))**2) + &
1362  (abs(amp_z(-1,+1))**2) + &
1363  (abs(amp_z(+1,-1))**2) + &
1364  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1365 
1366  if( .not. zz_fusion ) restmp = 0.0d0
1367 
1368  pdfindex = flip(iflip,pdfup_,pdfbot_)
1369  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1370 
1371  pdfindex = flip(iflip,pdfchm_,pdfbot_)
1372  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1373  return
1374  endif
1375 
1376 
1377  if( (isel.eq.pdfup_ .and. jsel.eq.pdfstr_) .or. (isel.eq.pdfchm_ .and. jsel.eq.pdfdn_) ) then
1378 
1379  if( zz_fusion ) then
1380  amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,isel,jsel)
1381  restmp = ((abs(amp_z(-1,-1))**2) + &
1382  (abs(amp_z(-1,+1))**2) + &
1383  (abs(amp_z(+1,-1))**2) + &
1384  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1385  else
1386  amp_w = a0_ww_4f(3,j2,4,j1,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.false.)
1387  restmp= ( &
1388  abs(amp_w(-1,-1))**2 + &
1389  abs(amp_w(-1,+1))**2 + &
1390  abs(amp_w(+1,-1))**2 + &
1391  abs(amp_w(+1,+1))**2 &
1392  ) * spinavg * tag2
1393  endif
1394 
1395  pdfindex = flip(iflip,pdfup_,pdfstr_)
1396  res(pdfindex(1),pdfindex(2)) = restmp
1397 
1398  pdfindex = flip(iflip,pdfchm_,pdfdn_)
1399  res(pdfindex(1),pdfindex(2)) = restmp
1400 
1401  return
1402  endif
1403 
1404 
1405  !--ds -> ds
1406  if( (isel.eq.pdfdn_ .and. jsel.eq.pdfstr_) .or. (isel.eq.pdfdn_ .and. jsel.eq.pdfbot_) .or. (isel.eq.pdfstr_ .and. jsel.eq.pdfbot_) ) then
1407 
1408  amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,isel,jsel)
1409  restmp = ((abs(amp_z(-1,-1))**2) + &
1410  (abs(amp_z(-1,+1))**2) + &
1411  (abs(amp_z(+1,-1))**2) + &
1412  (abs(amp_z(+1,+1))**2)) * spinavg
1413 
1414  if( .not. zz_fusion ) restmp = 0.0d0
1415 
1416  pdfindex = flip(iflip,pdfdn_,pdfstr_)
1417  res(pdfindex(1),pdfindex(2)) = restmp
1418 
1419  pdfindex = flip(iflip,pdfdn_,pdfbot_)
1420  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1421 
1422  pdfindex = flip(iflip,pdfstr_,pdfbot_)
1423  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1424  return
1425  endif
1426 
1427 
1428  !-- qbqb processes
1429  if( (isel.eq.pdfaup_ .and. jsel.eq.pdfachm_) ) then
1430  amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,isel,jsel)
1431  !--ubcb -> ubcb
1432  restmp = ((abs(amp_z(-1,-1))**2) + &
1433  (abs(amp_z(-1,+1))**2) + &
1434  (abs(amp_z(+1,-1))**2) + &
1435  (abs(amp_z(+1,+1))**2)) * spinavg
1436 
1437  if( .not. zz_fusion ) restmp = 0.0d0
1438 
1439 
1440  pdfindex = flip(iflip,pdfaup_,pdfachm_)
1441  res(pdfindex(1),pdfindex(2)) = restmp
1442  return
1443  endif
1444 
1445 
1446  !--ubsb -> ubsb//cbdb
1447  if( (isel.eq.pdfaup_ .and. jsel.eq.pdfabot_) .or. (isel.eq.pdfachm_ .and. jsel.eq.pdfabot_) ) then
1448  amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,isel,jsel)
1449  restmp = ((abs(amp_z(-1,-1))**2) + &
1450  (abs(amp_z(-1,+1))**2) + &
1451  (abs(amp_z(+1,-1))**2) + &
1452  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1453 
1454  if( .not. zz_fusion ) restmp = 0.0d0
1455 
1456  pdfindex = flip(iflip,pdfaup_,pdfabot_)
1457  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1458 
1459  pdfindex = flip(iflip,pdfachm_,pdfabot_)
1460  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1461  return
1462  endif
1463 
1464 
1465  if( (isel.eq.pdfaup_ .and. jsel.eq.pdfastr_) .or. (isel.eq.pdfachm_ .and. jsel.eq.pdfadn_) ) then
1466 
1467  if( zz_fusion ) then
1468  amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,isel,jsel)
1469  restmp = ((abs(amp_z(-1,-1))**2) + &
1470  (abs(amp_z(-1,+1))**2) + &
1471  (abs(amp_z(+1,-1))**2) + &
1472  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1473  else
1474  amp_w = a0_ww_4f(j1,4,j2,3,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.false.)
1475  restmp= ( &
1476  abs(amp_w(-1,-1))**2 + &
1477  abs(amp_w(-1,+1))**2 + &
1478  abs(amp_w(+1,-1))**2 + &
1479  abs(amp_w(+1,+1))**2 &
1480  ) * spinavg * tag2
1481  endif
1482 
1483  pdfindex = flip(iflip,pdfaup_,pdfastr_)
1484  res(pdfindex(1),pdfindex(2)) = restmp
1485 
1486  pdfindex = flip(iflip,pdfachm_,pdfadn_)
1487  res(pdfindex(1),pdfindex(2)) = restmp
1488 
1489 
1490  return
1491  endif
1492 
1493 
1494  !--dbsb -> dbsb
1495  if( (isel.eq.pdfadn_ .and. jsel.eq.pdfastr_) .or. (isel.eq.pdfadn_ .and. jsel.eq.pdfabot_) .or. (isel.eq.pdfastr_ .and. jsel.eq.pdfabot_) ) then
1496 
1497  amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,isel,jsel)
1498  restmp = ((abs(amp_z(-1,-1))**2) + &
1499  (abs(amp_z(-1,+1))**2) + &
1500  (abs(amp_z(+1,-1))**2) + &
1501  (abs(amp_z(+1,+1))**2)) * spinavg
1502 
1503  if( .not. zz_fusion ) restmp = 0.0d0
1504 
1505 
1506  pdfindex = flip(iflip,pdfadn_,pdfastr_)
1507  res(pdfindex(1),pdfindex(2)) = restmp
1508 
1509  pdfindex = flip(iflip,pdfadn_,pdfabot_)
1510  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1511 
1512  pdfindex = flip(iflip,pdfastr_,pdfabot_)
1513  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1514  return
1515  endif
1516 
1517  !-------------
1518 
1519  j1 = 2
1520  j2 = 1
1521  iflip = 2
1522  !-- qqb processes
1523 
1524  !--uub -> uub // ddb
1525  if( (jsel.eq.pdfup_ .and. isel.eq.pdfaup_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfachm_) .or. (jsel.eq.pdfup_ .and. isel.eq.pdfachm_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfaup_) ) then
1526 
1527  if( zz_fusion ) then
1528  amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,jsel,isel)
1529  restmp = ((abs(amp_z(-1,-1))**2) + &
1530  (abs(amp_z(-1,+1))**2) + &
1531  (abs(amp_z(+1,-1))**2) + &
1532  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1533  else
1534  amp_w = a0_ww_4f(4,j1,j2,3,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.true.)
1535  restmp= ( &
1536  abs(amp_w(-1,-1))**2 + &
1537  abs(amp_w(-1,+1))**2 + &
1538  abs(amp_w(+1,-1))**2 + &
1539  abs(amp_w(+1,+1))**2 &
1540  ) * spinavg * tag2
1541  endif
1542 
1543  pdfindex = flip(iflip,pdfup_,pdfaup_)
1544  res(pdfindex(1),pdfindex(2)) = restmp
1545 
1546  pdfindex = flip(iflip,pdfchm_,pdfachm_)
1547  res(pdfindex(1),pdfindex(2)) = restmp
1548 
1549  pdfindex = flip(iflip,pdfup_,pdfachm_)
1550  res(pdfindex(1),pdfindex(2)) = restmp
1551 
1552  pdfindex = flip(iflip,pdfchm_,pdfaup_)
1553  res(pdfindex(1),pdfindex(2)) = restmp
1554 
1555 
1556  return
1557  endif
1558 
1559 
1560  !--udb -> udb
1561  if( (jsel.eq.pdfup_ .and. isel.eq.pdfadn_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfastr_) .or. (jsel.eq.pdfup_ .and. isel.eq.pdfastr_) .or. &
1562  (jsel.eq.pdfchm_ .and. isel.eq.pdfadn_) .or. (jsel.eq.pdfup_ .and. isel.eq.pdfabot_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfabot_) ) then
1563  amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,jsel,isel)
1564  restmp = ((abs(amp_z(-1,-1))**2) + &
1565  (abs(amp_z(-1,+1))**2) + &
1566  (abs(amp_z(+1,-1))**2) + &
1567  (abs(amp_z(+1,+1))**2)) * spinavg
1568 
1569  if( .not. zz_fusion ) restmp = 0.0d0
1570 
1571 
1572  pdfindex = flip(iflip,pdfup_,pdfadn_)
1573  res(pdfindex(1),pdfindex(2)) = restmp
1574 
1575  pdfindex = flip(iflip,pdfchm_,pdfastr_)
1576  res(pdfindex(1),pdfindex(2)) = restmp
1577 
1578  pdfindex = flip(iflip,pdfup_,pdfastr_)
1579  res(pdfindex(1),pdfindex(2)) = restmp
1580 
1581  pdfindex = flip(iflip,pdfchm_,pdfadn_)
1582  res(pdfindex(1),pdfindex(2)) = restmp
1583 
1584  pdfindex = flip(iflip,pdfup_,pdfabot_)
1585  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1586 
1587  pdfindex = flip(iflip,pdfchm_,pdfabot_)
1588  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1589  return
1590  endif
1591 
1592 
1593  !--dub -> dub
1594  if( (jsel.eq.pdfdn_ .and. isel.eq.pdfaup_) .or. (jsel.eq.pdfstr_ .and. isel.eq.pdfachm_) .or. (jsel.eq.pdfdn_ .and. isel.eq.pdfachm_) .or. &
1595  (jsel.eq.pdfstr_ .and. isel.eq.pdfaup_) .or. (jsel.eq.pdfbot_ .and. isel.eq.pdfaup_) .or. (jsel.eq.pdfbot_ .and. isel.eq.pdfachm_) ) then
1596  amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,jsel,isel)
1597  restmp = ((abs(amp_z(-1,-1))**2) + &
1598  (abs(amp_z(-1,+1))**2) + &
1599  (abs(amp_z(+1,-1))**2) + &
1600  (abs(amp_z(+1,+1))**2)) * spinavg
1601 
1602  if( .not. zz_fusion ) restmp = 0.0d0
1603 
1604 
1605  pdfindex = flip(iflip,pdfdn_,pdfaup_)
1606  res(pdfindex(1),pdfindex(2)) = restmp
1607 
1608  pdfindex = flip(iflip,pdfstr_,pdfachm_)
1609  res(pdfindex(1),pdfindex(2)) = restmp
1610 
1611  pdfindex = flip(iflip,pdfdn_,pdfachm_)
1612  res(pdfindex(1),pdfindex(2)) = restmp
1613 
1614  pdfindex = flip(iflip,pdfstr_,pdfaup_)
1615  res(pdfindex(1),pdfindex(2)) = restmp
1616 
1617  pdfindex = flip(iflip,pdfbot_,pdfaup_)
1618  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1619 
1620  pdfindex = flip(iflip,pdfbot_,pdfachm_)
1621  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1622  return
1623  endif
1624 
1625 
1626  !--ddb -> uub/ddb
1627  if( (jsel.eq.pdfbot_ .and. isel.eq.pdfabot_) .or. (jsel.eq.pdfdn_ .and. isel.eq.pdfabot_) .or. (jsel.eq.pdfstr_ .and. isel.eq.pdfabot_) .or.&
1628  (jsel.eq.pdfbot_ .and. isel.eq.pdfadn_) .or. (jsel.eq.pdfbot_ .and. isel.eq.pdfastr_) ) then
1629  amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,jsel,isel)
1630 
1631  restmp = ((abs(amp_z(-1,-1))**2) + &
1632  (abs(amp_z(-1,+1))**2) + &
1633  (abs(amp_z(+1,-1))**2) + &
1634  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1635 
1636  if( .not. zz_fusion ) restmp = 0.0d0
1637 
1638 
1639  pdfindex = flip(iflip,pdfbot_,pdfabot_)
1640  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1641 
1642  pdfindex = flip(iflip,pdfdn_,pdfabot_)
1643  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1644 
1645  pdfindex = flip(iflip,pdfstr_,pdfabot_)
1646  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1647 
1648  pdfindex = flip(iflip,pdfbot_,pdfadn_)
1649  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1650 
1651  pdfindex = flip(iflip,pdfbot_,pdfastr_)
1652  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1653  return
1654  endif
1655 
1656 
1657 
1658 
1659  if( (jsel.eq.pdfdn_ .and. isel.eq.pdfadn_) .or. (jsel.eq.pdfstr_ .and. isel.eq.pdfastr_) .or. (jsel.eq.pdfdn_ .and. isel.eq.pdfastr_) .or. (jsel.eq.pdfstr_ .and. isel.eq.pdfadn_) ) then
1660  if( zz_fusion ) then
1661  amp_z = a0_zz_4f(3,j1,j2,4,za,zb,sprod,jsel,isel)
1662  restmp = ((abs(amp_z(-1,-1))**2) + &
1663  (abs(amp_z(-1,+1))**2) + &
1664  (abs(amp_z(+1,-1))**2) + &
1665  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1666  else
1667  amp_w = a0_ww_4f(4,j1,j2,3,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.false.)
1668  restmp= ( &
1669  abs(amp_w(-1,-1))**2 + &
1670  abs(amp_w(-1,+1))**2 + &
1671  abs(amp_w(+1,-1))**2 + &
1672  abs(amp_w(+1,+1))**2 &
1673  ) * spinavg * tag2
1674  endif
1675 
1676  pdfindex = flip(iflip,pdfdn_,pdfadn_)
1677  res(pdfindex(1),pdfindex(2)) = restmp
1678 
1679  pdfindex = flip(iflip,pdfstr_,pdfastr_)
1680  res(pdfindex(1),pdfindex(2)) = restmp
1681 
1682  pdfindex = flip(iflip,pdfdn_,pdfastr_)
1683  res(pdfindex(1),pdfindex(2)) = restmp
1684 
1685  pdfindex = flip(iflip,pdfstr_,pdfadn_)
1686  res(pdfindex(1),pdfindex(2)) = restmp
1687 
1688  return
1689  endif
1690 
1691 
1692  if( (jsel.eq.pdfup_ .and. isel.eq.pdfchm_) ) then
1693 
1694  !-- non-symmetric qq processes
1695  amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,jsel,isel)
1696  !--uc -> uc
1697  restmp = ((abs(amp_z(-1,-1))**2) + &
1698  (abs(amp_z(-1,+1))**2) + &
1699  (abs(amp_z(+1,-1))**2) + &
1700  (abs(amp_z(+1,+1))**2)) * spinavg
1701 
1702  if( .not. zz_fusion ) restmp = 0.0d0
1703 
1704  pdfindex = flip(iflip,pdfup_,pdfchm_)
1705  res(pdfindex(1),pdfindex(2)) = restmp
1706  return
1707  endif
1708 
1709 
1710  !--us -> us/cd
1711  if( (jsel.eq.pdfup_ .and. isel.eq.pdfbot_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfbot_) ) then
1712 
1713  amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,jsel,isel)
1714  restmp = ((abs(amp_z(-1,-1))**2) + &
1715  (abs(amp_z(-1,+1))**2) + &
1716  (abs(amp_z(+1,-1))**2) + &
1717  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1718 
1719  if( .not. zz_fusion ) restmp = 0.0d0
1720 
1721  pdfindex = flip(iflip,pdfup_,pdfbot_)
1722  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1723 
1724  pdfindex = flip(iflip,pdfchm_,pdfbot_)
1725  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1726  return
1727  endif
1728 
1729 
1730  if( (jsel.eq.pdfup_ .and. isel.eq.pdfstr_) .or. (jsel.eq.pdfchm_ .and. isel.eq.pdfdn_) ) then
1731 
1732  if( zz_fusion ) then
1733  amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,jsel,isel)
1734  restmp = ((abs(amp_z(-1,-1))**2) + &
1735  (abs(amp_z(-1,+1))**2) + &
1736  (abs(amp_z(+1,-1))**2) + &
1737  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1738  else
1739  amp_w = a0_ww_4f(3,j2,4,j1,za,zb,sprod,isel,jsel,usewwcoupl=.true.,wpm_flip=.false.)
1740  restmp= ( &
1741  abs(amp_w(-1,-1))**2 + &
1742  abs(amp_w(-1,+1))**2 + &
1743  abs(amp_w(+1,-1))**2 + &
1744  abs(amp_w(+1,+1))**2 &
1745  ) * spinavg * tag2
1746  endif
1747 
1748  pdfindex = flip(iflip,pdfup_,pdfstr_)
1749  res(pdfindex(1),pdfindex(2)) = restmp
1750 
1751  pdfindex = flip(iflip,pdfchm_,pdfdn_)
1752  res(pdfindex(1),pdfindex(2)) = restmp
1753 
1754  return
1755  endif
1756 
1757 
1758  !--ds -> ds
1759  if( (jsel.eq.pdfdn_ .and. isel.eq.pdfstr_) .or. (jsel.eq.pdfdn_ .and. isel.eq.pdfbot_) .or. (jsel.eq.pdfstr_ .and. isel.eq.pdfbot_) ) then
1760 
1761  amp_z = a0_zz_4f(3,j1,4,j2,za,zb,sprod,jsel,isel)
1762  restmp = ((abs(amp_z(-1,-1))**2) + &
1763  (abs(amp_z(-1,+1))**2) + &
1764  (abs(amp_z(+1,-1))**2) + &
1765  (abs(amp_z(+1,+1))**2)) * spinavg
1766 
1767  if( .not. zz_fusion ) restmp = 0.0d0
1768 
1769  pdfindex = flip(iflip,pdfdn_,pdfstr_)
1770  res(pdfindex(1),pdfindex(2)) = restmp
1771 
1772  pdfindex = flip(iflip,pdfdn_,pdfbot_)
1773  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1774 
1775  pdfindex = flip(iflip,pdfstr_,pdfbot_)
1776  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1777  return
1778  endif
1779 
1780 
1781  !-- qbqb processes
1782  if( (jsel.eq.pdfaup_ .and. isel.eq.pdfachm_) ) then
1783  amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,jsel,isel)
1784  !--ubcb -> ubcb
1785  restmp = ((abs(amp_z(-1,-1))**2)+ &
1786  (abs(amp_z(-1,+1))**2) + &
1787  (abs(amp_z(+1,-1))**2) + &
1788  (abs(amp_z(+1,+1))**2)) * spinavg
1789 
1790  if( .not. zz_fusion ) restmp = 0.0d0
1791 
1792 
1793  pdfindex = flip(iflip,pdfaup_,pdfachm_)
1794  res(pdfindex(1),pdfindex(2)) = restmp
1795  return
1796  endif
1797 
1798 
1799  !--ubsb -> ubsb//cbdb
1800  if( (jsel.eq.pdfaup_ .and. isel.eq.pdfabot_) .or. (jsel.eq.pdfachm_ .and. isel.eq.pdfabot_) ) then
1801  amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,jsel,isel)
1802  restmp = ((abs(amp_z(-1,-1))**2) + &
1803  (abs(amp_z(-1,+1))**2) + &
1804  (abs(amp_z(+1,-1))**2) + &
1805  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1806 
1807  if( .not. zz_fusion ) restmp = 0.0d0
1808 
1809 
1810  pdfindex = flip(iflip,pdfaup_,pdfabot_)
1811  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1812 
1813  pdfindex = flip(iflip,pdfachm_,pdfabot_)
1814  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1815  return
1816  endif
1817 
1818 
1819  if( (jsel.eq.pdfaup_ .and. isel.eq.pdfastr_) .or. (jsel.eq.pdfachm_ .and. isel.eq.pdfadn_) ) then
1820  if( zz_fusion ) then
1821  amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,jsel,isel)
1822  restmp = ((abs(amp_z(-1,-1))**2) + &
1823  (abs(amp_z(-1,+1))**2) + &
1824  (abs(amp_z(+1,-1))**2) + &
1825  (abs(amp_z(+1,+1))**2)) * spinavg * tag1
1826  else
1827  amp_w = a0_ww_4f(j1,4,j2,3,za,zb,sprod,jsel,isel,usewwcoupl=.true.,wpm_flip=.false.)
1828  restmp= ( &
1829  abs(amp_w(-1,-1))**2 + &
1830  abs(amp_w(-1,+1))**2 + &
1831  abs(amp_w(+1,-1))**2 + &
1832  abs(amp_w(+1,+1))**2 &
1833  ) * spinavg * tag2
1834  endif
1835 
1836  pdfindex = flip(iflip,pdfaup_,pdfastr_)
1837  res(pdfindex(1),pdfindex(2)) = restmp
1838 
1839  pdfindex = flip(iflip,pdfachm_,pdfadn_)
1840  res(pdfindex(1),pdfindex(2)) = restmp
1841 
1842  return
1843  endif
1844 
1845 
1846  !--dbsb -> dbsb
1847  if( (jsel.eq.pdfadn_ .and. isel.eq.pdfastr_) .or. (jsel.eq.pdfadn_ .and. isel.eq.pdfabot_) .or. (jsel.eq.pdfastr_ .and. isel.eq.pdfabot_) ) then
1848 
1849  amp_z = a0_zz_4f(j1,3,j2,4,za,zb,sprod,jsel,isel)
1850  restmp = ((abs(amp_z(-1,-1))**2) + &
1851  (abs(amp_z(-1,+1))**2) + &
1852  (abs(amp_z(+1,-1))**2) + &
1853  (abs(amp_z(+1,+1))**2)) * spinavg
1854 
1855  if( .not. zz_fusion ) restmp = 0.0d0
1856 
1857  pdfindex = flip(iflip,pdfadn_,pdfastr_)
1858  res(pdfindex(1),pdfindex(2)) = restmp
1859 
1860  pdfindex = flip(iflip,pdfadn_,pdfabot_)
1861  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1862 
1863  pdfindex = flip(iflip,pdfastr_,pdfabot_)
1864  res(pdfindex(1),pdfindex(2)) = restmp * tagbot
1865  return
1866  endif
1867 
1868  RETURN
1869  END SUBROUTINE evalamp_wbfh_unsymm_sa_select
1870 
1871  SUBROUTINE evalamp_wbfh_unsymm_sa_select_exact(p,iSel,jSel,rSel,sSel,res)
1872  implicit none
1873  real(dp), intent(in) :: p(4,5)
1874  real(dp), intent(out) :: res
1875  integer, intent(in) :: isel,jsel,rsel,ssel
1876  complex(dp) :: amp_z(-1:1,-1:1), amp_z_b(-1:1,-1:1)
1877  complex(dp) :: amp_w(-1:1,-1:1)
1878  real(dp) :: sprod(4,4)
1879  complex(dp) :: za(4,4), zb(4,4)
1880  real(dp) :: restmp
1881  real(dp) :: ckm_wfactor
1882  real(dp) :: kfactor_z,kfactor_w
1883  integer :: i, j, jz1, jz2, jw1, jw2, pdfindex(2)
1884  integer :: line_i,line_j
1885  integer :: kz1, kz2, kw1, kw2
1886  logical :: zz_fusion,ww_fusion
1887 
1888  !print *, "Begin EvalAmp_WBFH_UnSymm_SA_Select_exact"
1889  !print *, "iSel: ",iSel,", jSel: ",jSel," rSel: ",rSel," sSel: ",sSel
1890 
1891  res=0.0_dp
1892 
1893  zz_fusion=.false.
1894  ww_fusion=.false.
1895 
1896  ckm_wfactor=1.0_dp
1897  kfactor_z=1.0_dp
1898  kfactor_w=1.0_dp
1899  restmp=0.0_dp
1900 
1901  line_i=0
1902  line_j=0
1903  jz1 = 1
1904  jz2 = 2
1905  jw1 = jz1
1906  jw2 = jz2
1907 
1908  if( (abs(isel).eq.pdftop_ .or. abs(jsel).eq.pdftop_) .or. (abs(rsel).eq.pdftop_ .or. abs(ssel).eq.pdftop_) ) return
1909  if( (abs(isel).eq.pdfglu_ .or. abs(jsel).eq.pdfglu_) .or. (abs(rsel).eq.pdfglu_ .or. abs(ssel).eq.pdfglu_) ) return
1910 
1911  if( &
1912  (isel.eq.rsel .and. jsel.eq.ssel) &
1913  .or. &
1914  (isel.eq.ssel .and. jsel.eq.rsel) &
1915  ) then
1916  zz_fusion=.true.
1917  if( isel.eq.ssel .and. jsel.eq.rsel ) then
1918  kz1 = 4
1919  kz2 = 3
1920  kfactor_z = 1.0_dp ! dsqrt(ScaleFactor(iSel,sSel)*ScaleFactor(jSel,rSel))
1921  !print *, "EvalAmp_WBFH_UnSymm_SA_Select_exact: isZZ and is-jr"
1922  else
1923  kz1 = 3
1924  kz2 = 4
1925  kfactor_z = 1.0_dp ! dsqrt(ScaleFactor(iSel,rSel)*ScaleFactor(jSel,sSel))
1926  !print *, "EvalAmp_WBFH_UnSymm_SA_Select_exact: isZZ and ir-js"
1927  endif
1928  endif
1929 
1930  if( &
1931  ( (sign(isel,rsel).eq.isel .and. sign(jsel,ssel).eq.jsel) .and. (abs(isel-rsel).eq.1 .or. abs(isel-rsel).eq.3 .or. abs(isel-rsel).eq.5) .and. (abs(jsel-ssel).eq.1 .or. abs(jsel-ssel).eq.3 .or. abs(jsel-ssel).eq.5) ) &
1932  .or. &
1933  ( (sign(isel,ssel).eq.isel .and. sign(jsel,rsel).eq.jsel) .and. (abs(isel-ssel).eq.1 .or. abs(isel-ssel).eq.3 .or. abs(isel-ssel).eq.5) .and. (abs(jsel-rsel).eq.1 .or. abs(jsel-rsel).eq.3 .or. abs(jsel-rsel).eq.5) ) &
1934  ) then
1935  ww_fusion=.true.
1936  ! W_is W_jr fusion
1937  if( (sign(isel,ssel).eq.isel .and. sign(jsel,rsel).eq.jsel) .and. (abs(isel-ssel).eq.1 .or. abs(isel-ssel).eq.3 .or. abs(isel-ssel).eq.5) .and. (abs(jsel-rsel).eq.1 .or. abs(jsel-rsel).eq.3 .or. abs(jsel-rsel).eq.5) ) then
1938  kw1 = 4
1939  kw2 = 3
1940  kfactor_w = 1.0_dp ! dsqrt(ScaleFactor(iSel,sSel)*ScaleFactor(jSel,rSel))
1941 
1942  !if(ZZ_fusion) then ! Special treatment for WW+ZZ interference, not included through phasespace
1943  ckm_wfactor = ckmbare(isel,ssel)*ckmbare(jsel,rsel)
1944  !print *, "iSel, sSel: ",iSel," ",sSel, "; jSel, rSel: ",jSel," ",rSel, ", ckm: ",ckm_wfactor
1945  !endif
1946  !print *, "EvalAmp_WBFH_UnSymm_SA_Select_exact: isWW and is-jr"
1947  else ! W_ir W_js fusion
1948  kw1 = 3
1949  kw2 = 4
1950  kfactor_w = 1.0_dp ! dsqrt(ScaleFactor(iSel,rSel)*ScaleFactor(jSel,sSel))
1951 
1952  ckm_wfactor = ckmbare(isel,rsel)*ckmbare(jsel,ssel)
1953  !print *, "iSel, rSel: ",iSel," ",rSel, "; jSel, sSel: ",jSel," ",sSel, ", ckm: ",ckm_wfactor
1954  !print *, "EvalAmp_WBFH_UnSymm_SA_Select_exact: isWW and ir-js"
1955  endif
1956  endif
1957 
1958 
1959  if( .not.(zz_fusion .or. ww_fusion) ) return
1960  if(isel.lt.0) then
1961  if(zz_fusion) call swapi(jz1,kz1)
1962  if(ww_fusion) call swapi(jw1,kw1)
1963  endif
1964  if(jsel.lt.0) then
1965  if(zz_fusion) call swapi(jz2,kz2)
1966  if(ww_fusion) call swapi(jw2,kw2)
1967  endif
1968 
1969  line_i = isel
1970  line_j = jsel
1971 
1972  if(ww_fusion) then
1973  if (isel.eq.pdfup_ .or. isel.eq.pdfchm_ .or. isel.eq.pdfadn_ .or. isel.eq.pdfastr_ .or. isel.eq.pdfabot_) then ! W+ should be passed as the second set of partons
1974  call swapi(jw1,jw2)
1975  call swapi(kw1,kw2)
1976  call swapi(line_i,line_j)
1977  if(zz_fusion) then ! If also ZZ fusion, swap it as well
1978  call swapi(jz1,jz2)
1979  call swapi(kz1,kz2)
1980  endif
1981  endif
1982  endif
1983 
1984  call spinoru2(4,(/-p(:,1),-p(:,2),p(:,3),p(:,4)/),za,zb,sprod)
1985 
1986  !if(ZZ_fusion) print *, "jz1: ",jz1,", jz2: ",jz2," kz1: ",kz1," kz2: ",kz2
1987  !if(WW_fusion) print *, "jw1: ",jw1,", jw2: ",jw2," kw1: ",kw1," kw2: ",kw2
1988 
1989  !ckm_wfactor = 1.0_dp ! TEST!!!
1990  if(zz_fusion) then
1991 
1992  amp_z = a0_zz_4f(kz1,jz1,kz2,jz2,za,zb,sprod,line_i,line_j)
1993 
1994  restmp = restmp + &
1995  ((abs(amp_z(-1,-1))**2) + &
1996  (abs(amp_z(-1,+1))**2) + &
1997  (abs(amp_z(+1,-1))**2) + &
1998  (abs(amp_z(+1,+1))**2)) * xn**2
1999 
2000  if(isel.eq.jsel) then
2001  amp_z_b = -a0_zz_4f(kz2,jz1,kz1,jz2,za,zb,sprod,line_i,line_j)
2002  restmp = restmp + &
2003  ((abs(amp_z_b(-1,-1))**2) + &
2004  (abs(amp_z_b(-1,+1))**2) + &
2005  (abs(amp_z_b(+1,-1))**2) + &
2006  (abs(amp_z_b(+1,+1))**2)) * xn**2
2007  restmp = restmp + &
2008  (two * real(amp_z(-1,-1)*conjg(amp_z_b(-1,-1)),kind=dp) + &
2009  two * real(amp_z(+1,+1)*conjg(amp_z_b(+1,+1)),kind=dp)) * xn
2010 
2011  restmp = restmp * symmfac
2012  endif
2013  restmp = restmp * kfactor_z**2
2014  endif
2015 
2016  if(ww_fusion) then
2017  amp_w = -a0_ww_4f(kw1,jw1,kw2,jw2,za,zb,sprod,line_i,line_j,usewwcoupl=.true.)
2018  restmp = restmp + &
2019  ( &
2020  abs(amp_w(-1,-1))**2 + &
2021  abs(amp_w(-1,+1))**2 + &
2022  abs(amp_w(+1,-1))**2 + &
2023  abs(amp_w(+1,+1))**2 &
2024  ) * xn**2 * kfactor_w**2 * abs(ckm_wfactor)**2
2025  if(zz_fusion) then
2026  restmp = restmp + &
2027  two * ( &
2028  real(amp_w(-1,-1)*ckm_wfactor*conjg(amp_z(-1,-1)),kind=dp) + &
2029  real(amp_w(-1,+1)*ckm_wfactor*conjg(amp_z(-1,+1)),kind=dp) + &
2030  real(amp_w(+1,-1)*ckm_wfactor*conjg(amp_z(+1,-1)),kind=dp) + &
2031  real(amp_w(+1,+1)*ckm_wfactor*conjg(amp_z(+1,+1)),kind=dp) &
2032  ) * xn * kfactor_z * kfactor_w
2033  endif
2034  endif
2035 
2036  restmp = restmp * aveqq
2037  if(abs(isel).eq.pdfbot_ .or. abs(jsel).eq.pdfbot_) restmp = restmp * tagbot
2038  res = restmp
2039 
2040  RETURN
2041  END SUBROUTINE evalamp_wbfh_unsymm_sa_select_exact
2042 
2043 
2044 
2045 
2046 !-- TEST FUNCTION THAT WRAPS THE HVV AMPLITUDE
2047 !-- FOR IT TO WORK, COMMENT OUT THE LINES BELOW,
2048 !-- AND THE QUOTED LINE IN MODHIGGS::CALCHELAMP2:
2049 !-- "if ((q_q).lt.-0.1d0 .or. (q3_q3).lt.-0.1d0 .or. (q4_q4).lt.-0.1d0) return ! if negative invariant masses return zero"
2050 
2051 ! SUBROUTINE wrapHVV(p,iSel,jSel,rSel,sSel,res)
2052 ! use ModHiggs
2053 ! implicit none
2054 ! real(dp), intent(in) :: p(4,5)
2055 ! real(dp), intent(out) :: res
2056 ! integer, intent(in) :: iSel,jSel,rSel,sSel
2057 ! integer :: i, j, jz1, jz2, jw1, jw2, pdfindex(2)
2058 ! integer :: kz1, kz2, kw1, kw2
2059 ! logical :: ZZ_fusion,WW_fusion
2060 ! complex(dp) :: A_VV(1:8)
2061 ! integer, parameter :: ZZMode=00,ZgsMode=01,gsZMode=02,gsgsMode=03
2062 ! integer, parameter :: WWMode=10
2063 ! integer, parameter :: ggMode=20
2064 ! integer, parameter :: ZgMode=30,gsgMode=31
2065 ! integer :: MY_IDUP(3:6),i3,i4
2066 ! real(dp) :: pUsed(4,6)
2067 
2068 ! if( (abs(iSel).eq.pdfTop_ .or. abs(jSel).eq.pdfTop_) .or. (abs(rSel).eq.pdfTop_ .or. abs(sSel).eq.pdfTop_) ) return
2069 
2070 ! res=0.0_dp
2071 
2072 ! ZZ_fusion=.false.
2073 ! WW_fusion=.false.
2074 
2075 ! jz1 = 1
2076 ! jz2 = 2
2077 ! jw1 = jz1
2078 ! jw2 = jz2
2079 
2080 ! if( &
2081 ! (iSel.eq.rSel .and. jSel.eq.sSel) &
2082 ! .or. &
2083 ! (iSel.eq.sSel .and. jSel.eq.rSel) &
2084 ! ) then
2085 ! ZZ_fusion=.true.
2086 ! if( iSel.eq.sSel .and. jSel.eq.rSel ) then
2087 ! kz1 = 4
2088 ! kz2 = 3
2089 ! else
2090 ! kz1 = 3
2091 ! kz2 = 4
2092 ! endif
2093 ! endif
2094 
2095 ! if( &
2096 ! ( (sign(iSel,rSel).eq.iSel .and. sign(jSel,sSel).eq.jSel) .and. (abs(iSel-rSel).eq.1 .or. abs(iSel-rSel).eq.3 .or. abs(iSel-rSel).eq.5) .and. (abs(jSel-sSel).eq.1 .or. abs(jSel-sSel).eq.3 .or. abs(jSel-sSel).eq.5) ) &
2097 ! .or. &
2098 ! ( (sign(iSel,sSel).eq.iSel .and. sign(jSel,rSel).eq.jSel) .and. (abs(iSel-sSel).eq.1 .or. abs(iSel-sSel).eq.3 .or. abs(iSel-sSel).eq.5) .and. (abs(jSel-rSel).eq.1 .or. abs(jSel-rSel).eq.3 .or. abs(jSel-rSel).eq.5) ) &
2099 ! ) then
2100 ! WW_fusion=.true.
2101 ! ! W_is W_jr fusion
2102 ! if( (sign(iSel,sSel).eq.iSel .and. sign(jSel,rSel).eq.jSel) .and. (abs(iSel-sSel).eq.1 .or. abs(iSel-sSel).eq.3 .or. abs(iSel-sSel).eq.5) .and. (abs(jSel-rSel).eq.1 .or. abs(jSel-rSel).eq.3 .or. abs(jSel-rSel).eq.5) ) then
2103 ! kw1 = 4
2104 ! kw2 = 3
2105 ! else ! W_ir W_js fusion
2106 ! kw1 = 3
2107 ! kw2 = 4
2108 ! endif
2109 ! endif
2110 
2111 
2112 ! if( .not.(ZZ_fusion .or. WW_fusion) ) return
2113 ! if(iSel.lt.0) then
2114 ! if(ZZ_fusion) call swapi(jz1,kz1)
2115 ! if(WW_fusion) call swapi(jw1,kw1)
2116 ! endif
2117 ! if(jSel.lt.0) then
2118 ! if(ZZ_fusion) call swapi(jz2,kz2)
2119 ! if(WW_fusion) call swapi(jw2,kw2)
2120 ! endif
2121 
2122 ! if(WW_fusion) then
2123 ! if (iSel.eq.pdfUp_ .or. iSel.eq.pdfChm_ .or. iSel.eq.pdfADn_ .or. iSel.eq.pdfAStr_ .or. iSel.eq.pdfABot_) then ! W+ should be passed as the second set of partons
2124 ! call swapi(jw1,jw2)
2125 ! call swapi(kw1,kw2)
2126 ! if(ZZ_fusion) then ! If also ZZ fusion, swap it as well
2127 ! call swapi(jz1,jz2)
2128 ! call swapi(kz1,kz2)
2129 ! endif
2130 ! endif
2131 ! endif
2132 
2133 ! if(jz1.eq.1) then
2134 ! MY_IDUP(4) = -convertFromPartIndex(iSel)
2135 ! pUsed(:,4) = -p(:,1)
2136 ! else if(jz1.eq.2) then
2137 ! MY_IDUP(4) = -convertFromPartIndex(jSel)
2138 ! pUsed(:,4) = -p(:,2)
2139 ! else if(jz1.eq.3) then
2140 ! MY_IDUP(4) = convertFromPartIndex(rSel)
2141 ! pUsed(:,4) = p(:,3)
2142 ! else if(jz1.eq.4) then
2143 ! MY_IDUP(4) = convertFromPartIndex(sSel)
2144 ! pUsed(:,4) = p(:,4)
2145 ! endif
2146 ! if(kz1.eq.1) then
2147 ! MY_IDUP(3) = -convertFromPartIndex(iSel)
2148 ! pUsed(:,3) = -p(:,1)
2149 ! else if(kz1.eq.2) then
2150 ! MY_IDUP(3) = -convertFromPartIndex(jSel)
2151 ! pUsed(:,3) = -p(:,2)
2152 ! else if(kz1.eq.3) then
2153 ! MY_IDUP(3) = convertFromPartIndex(rSel)
2154 ! pUsed(:,3) = p(:,3)
2155 ! else if(kz1.eq.4) then
2156 ! MY_IDUP(3) = convertFromPartIndex(sSel)
2157 ! pUsed(:,3) = p(:,4)
2158 ! endif
2159 ! if(jz2.eq.1) then
2160 ! MY_IDUP(6) = -convertFromPartIndex(iSel)
2161 ! pUsed(:,6) = -p(:,1)
2162 ! else if(jz2.eq.2) then
2163 ! MY_IDUP(6) = -convertFromPartIndex(jSel)
2164 ! pUsed(:,6) = -p(:,2)
2165 ! else if(jz2.eq.3) then
2166 ! MY_IDUP(6) = convertFromPartIndex(rSel)
2167 ! pUsed(:,6) = p(:,3)
2168 ! else if(jz2.eq.4) then
2169 ! MY_IDUP(6) = convertFromPartIndex(sSel)
2170 ! pUsed(:,6) = p(:,4)
2171 ! endif
2172 ! if(kz2.eq.1) then
2173 ! MY_IDUP(5) = -convertFromPartIndex(iSel)
2174 ! pUsed(:,5) = -p(:,1)
2175 ! else if(kz2.eq.2) then
2176 ! MY_IDUP(5) = -convertFromPartIndex(jSel)
2177 ! pUsed(:,5) = -p(:,2)
2178 ! else if(kz2.eq.3) then
2179 ! MY_IDUP(5) = convertFromPartIndex(rSel)
2180 ! pUsed(:,5) = p(:,3)
2181 ! else if(kz2.eq.4) then
2182 ! MY_IDUP(5) = convertFromPartIndex(sSel)
2183 ! pUsed(:,5) = p(:,4)
2184 ! endif
2185 
2186 ! pUSed(:,1) = pUSed(:,3) + pUSed(:,4) + pUSed(:,5) + pUSed(:,6)
2187 ! pUSed(:,2) = zero
2188 
2189 ! if(ZZ_fusion) then
2190 ! print *, "iSel: ",iSel,", jSel: ",jSel," rSel: ",rSel," sSel: ",sSel
2191 ! print *, "jz1: ",jz1,", jz2: ",jz2," kz1: ",kz1," kz2: ",kz2
2192 ! print *,"MY_IDUP:",MY_IDUP
2193 ! print *,"p1:",p(:,1)
2194 ! print *,"p2:",p(:,2)
2195 ! print *,"p3:",p(:,3)
2196 ! print *,"p4:",p(:,4)
2197 ! print *,"p5:",p(:,5)
2198 ! print *,"pUsed1:",pUsed(:,1)
2199 ! print *,"pUsed2:",pUsed(:,2)
2200 ! print *,"pUsed3:",pUsed(:,3)
2201 ! print *,"pUsed4:",pUsed(:,4)
2202 ! print *,"pUsed5:",pUsed(:,5)
2203 ! print *,"pUsed6:",pUsed(:,6)
2204 
2205 ! A_VV(:) = 0d0
2206 ! do i3=1,2; do i4=1,2! sum over helicities
2207 ! call calcHelAmp2((/3,4,5,6/),ZZMode,MY_IDUP,pUsed,i3,i4,A_VV(1))
2208 ! if( includeGammaStar ) then
2209 ! call calcHelAmp2((/3,4,5,6/),ZgsMode,MY_IDUP,pUsed,i3,i4,A_VV(3))
2210 ! call calcHelAmp2((/3,4,5,6/),gsZMode,MY_IDUP,pUsed,i3,i4,A_VV(5))
2211 ! call calcHelAmp2((/3,4,5,6/),gsgsMode,MY_IDUP,pUsed,i3,i4,A_VV(7))
2212 ! endif
2213 
2214 ! if( (MY_IDUP(3).eq.MY_IDUP(5)) .and. (MY_IDUP(4).eq.MY_IDUP(6)) .and. iSel.eq.jSel ) then ! iSel.eq.jSel to avoid ZH diagram
2215 ! call calcHelAmp2((/5,4,3,6/),ZZMode,MY_IDUP,pUsed,i3,i4,A_VV(2))
2216 ! if( includeGammaStar ) then
2217 ! call calcHelAmp2((/5,4,3,6/),ZgsMode,MY_IDUP,pUsed,i3,i4,A_VV(4))
2218 ! call calcHelAmp2((/5,4,3,6/),gsZMode,MY_IDUP,pUsed,i3,i4,A_VV(6))
2219 ! call calcHelAmp2((/5,4,3,6/),gsgsMode,MY_IDUP,pUsed,i3,i4,A_VV(8))
2220 ! endif
2221 ! A_VV(2) = -A_VV(2)! minus from Fermi statistics
2222 ! A_VV(4) = -A_VV(4)
2223 ! A_VV(6) = -A_VV(6)
2224 ! A_VV(8) = -A_VV(8)
2225 ! endif
2226 
2227 ! res = res + (A_VV(1)+A_VV(3)+A_VV(5)+A_VV(7))*dconjg(A_VV(1)+A_VV(3)+A_VV(5)+A_VV(7))! interfere the 3456 pieces
2228 ! res = res + (A_VV(2)+A_VV(4)+A_VV(6)+A_VV(8))*dconjg(A_VV(2)+A_VV(4)+A_VV(6)+A_VV(8))! interfere the 5436 pieces
2229 ! if( (MY_IDUP(3).eq.MY_IDUP(5)) .and. (MY_IDUP(4).eq.MY_IDUP(6)) .and. (i3.eq.i4) .and. iSel.eq.jSel ) then! interfere the 3456 with 5436 pieces
2230 ! res = res + 2d0/3d0*dreal( A_VV(1)*dconjg( A_VV(2)+A_VV(4)+A_VV(6)+A_VV(8) ) )
2231 ! res = res + 2d0/3d0*dreal( A_VV(3)*dconjg( A_VV(2)+A_VV(4)+A_VV(6)+A_VV(8) ) )
2232 ! res = res + 2d0/3d0*dreal( A_VV(5)*dconjg( A_VV(2)+A_VV(4)+A_VV(6)+A_VV(8) ) )
2233 ! res = res + 2d0/3d0*dreal( A_VV(7)*dconjg( A_VV(2)+A_VV(4)+A_VV(6)+A_VV(8) ) )
2234 ! endif
2235 
2236 ! print *,"i3 i4 res: ",i3,i4,res
2237 ! enddo; enddo
2238 ! if( (MY_IDUP(3).eq.MY_IDUP(5)) .and. (MY_IDUP(4).eq.MY_IDUP(6)) .and. iSel.eq.jSel ) res = res/2d0
2239 
2240 ! endif
2241 ! res = res * aveqq
2242 
2243 ! RETURN
2244 ! END SUBROUTINE
2245 
2246 
2247 
2248 
2249 
2250 
2251  function flip(i,a1,a2)
2252  integer :: flip(2)
2253  integer :: i, a1, a2
2254 
2255  if (i .eq. 1) flip = (/a1,a2/)
2256  if (i .eq. 2) flip = (/a2,a1/)
2257 
2258  return
2259 
2260  end function flip
2261 
2262 
2263  !-------------------------------------------------------------------------
2264  !-- amplitudes below
2265  !-- Older WBF code uses A0_VV_4f, newer subroutines use A0_ZZ_4f or A0_WW_4f -- U.S.: Fully checked with decay amplitudes through a wrapper
2266  function a0_vv_4f(j1,j2,j3,j4,za,zb,sprod,mv,ga_v,useWWcoupl,Wpm_flip)
2267  use modmisc
2268  implicit none
2269  complex(dp) :: a0_vv_4f(-1:1,-1:1)
2270  integer :: j1,j2,j3,j4
2271  complex(dp) :: za(4,4), zb(4,4)
2272  real(dp) :: mv, ga_v
2273  logical,optional :: usewwcoupl,wpm_flip
2274  real(dp) :: sprod(4,4),q2wplus,q2wminus
2275  real(dp) :: mhsq, q1q2, kcoupl
2276  complex(dp) :: a1, a2, a3, struc1, struc2, struc3
2277  complex(dp) :: zab2
2278  complex(dp) :: iprop12, iprop34
2279  complex(dp) :: vvcoupl_prime(4)
2280  integer :: vv_it
2281  integer :: i,j,k,l
2282 
2283  zab2(j1,j2,j3,j4) = za(j1,j2)*zb(j2,j4) + za(j1,j3)*zb(j3,j4)
2284  !The previous line works, and assigns the whole zab2 correctly.
2285  !I have no idea how it works.
2286  !If you don't believe me, please uncomment the following lines:
2287  !do i=1,4
2288  ! do j=1,4
2289  ! do k=1,4
2290  ! do l=1,4
2291  ! print *,i,j,k,l,zab2(i,j,k,l),za(i,j)*zb(j,l) + za(i,k)*zb(k,l)
2292  ! enddo
2293  ! enddo
2294  ! enddo
2295  !enddo
2296  !print *,j1,j2,j3,j4
2297  !pause
2298 
2299  a0_vv_4f = czero
2300 
2301  q1q2 = (sprod(j1,j3)+sprod(j1,j4)+sprod(j2,j3)+sprod(j2,j4))/two
2302  mhsq = two * q1q2 + sprod(j1,j2) + sprod(j3,j4)
2303 
2304  kcoupl = q1q2/lambda**2
2305 
2306  q2wplus = sprod(j1,j2)
2307  q2wminus = sprod(j3,j4)
2308  if( present(wpm_flip) ) then
2309  if( wpm_flip ) call swapr(q2wplus,q2wminus)
2310  endif
2311  if( .not.present(usewwcoupl) ) then
2312  do vv_it=1,4
2313  vvcoupl_prime(vv_it) = hvvspinzerodynamiccoupling(vv_it,q2wplus,q2wminus,mhsq)
2314  enddo
2315  else
2316  do vv_it=1,4
2317  vvcoupl_prime(vv_it) = hvvspinzerodynamiccoupling(vv_it,q2wplus,q2wminus,mhsq,trywwcoupl=usewwcoupl)
2318  enddo
2319  endif
2320 
2321 
2322  a1 = vvcoupl_prime(1) * mv**2/mhsq + vvcoupl_prime(2) * two * q1q2/mhsq + vvcoupl_prime(3) * kcoupl * q1q2/mhsq
2323  a2 = -two * vvcoupl_prime(2) - kcoupl * vvcoupl_prime(3)
2324  a3 = -two * vvcoupl_prime(4)
2325 
2326  struc1 = two * (a1 * mhsq - ci * a3 * q1q2)
2327  struc2 = a2 + ci * a3
2328  struc3 = two * ci * a3
2329 
2330 
2331  a0_vv_4f(-1,-1) = za(j1,j3)*zb(j4,j2) * struc1 + &
2332  zab2(j1,j3,j4,j2)*zab2(j3,j1,j2,j4) * struc2 + &
2333  za(j1,j2)*za(j3,j4)*zb(j4,j2)**2 * struc3
2334 
2335  a0_vv_4f(-1,+1) = za(j1,j4)*zb(j3,j2) * struc1 + &
2336  zab2(j1,j3,j4,j2)*zab2(j4,j1,j2,j3) * struc2 + &
2337  za(j1,j2)*za(j4,j3)*zb(j3,j2)**2 * struc3
2338 
2339  a0_vv_4f(+1,-1) = za(j2,j3)*zb(j4,j1) * struc1 + &
2340  zab2(j2,j3,j4,j1)*zab2(j3,j1,j2,j4) * struc2 + &
2341  za(j2,j1)*za(j3,j4)*zb(j4,j1)**2 * struc3
2342 
2343  a0_vv_4f(+1,+1) = za(j2,j4)*zb(j3,j1) * struc1 + &
2344  zab2(j2,j3,j4,j1)*zab2(j4,j1,j2,j3) * struc2 + &
2345  za(j2,j1)*za(j4,j3)*zb(j3,j1)**2 * struc3
2346 
2347  iprop12 = sprod(j1,j2) - mv**2 + ci * mv * ga_v
2348  iprop34 = sprod(j3,j4) - mv**2 + ci * mv * ga_v
2349 
2350  a0_vv_4f = a0_vv_4f/vev /iprop12/iprop34
2351 
2352  return
2353 
2354  end function a0_vv_4f
2355 
2356  !-- line = 1/2 --> down/up couplings included
2357  function a0_zz_4f(j1,j2,j3,j4,za,zb,sprod,iSel,jSel)
2358  use modmisc
2359  implicit none
2360  real(dp), dimension(5) :: lz,rz,lct,rct
2361  real(dp), parameter, dimension(5) :: la = (/qdl,qul,qdl,qul,qdl/)
2362  real(dp), parameter, dimension(5) :: ra = (/qdr,qur,qdr,qur,qdr/)
2363  complex(dp) :: a0_zz_4f(-1:1,-1:1)
2364  integer :: j1,j2,j3,j4,isel,jsel,line1,line2
2365  complex(dp) :: za(4,4),zb(4,4)
2366  real(dp) :: sprod(4,4)
2367  real(dp) :: mhsq,q1q2,kcoupl,q12sq,q34sq
2368  complex(dp) :: a1_zz,a2_zz,a3_zz,struc_zz(3),fac_zz(-1:1,-1:1)
2369  complex(dp) :: a1_az,a2_az,a3_az,struc_az(3),fac_az(-1:1,-1:1)
2370  complex(dp) :: a1_za,a2_za,a3_za,struc_za(3),fac_za(-1:1,-1:1)
2371  complex(dp) :: a1_aa,a2_aa,a3_aa,struc_aa(3),fac_aa(-1:1,-1:1)
2372  complex(dp) :: a1_zzp,a2_zzp,a3_zzp,struc_zzp(3),fac_zzp(-1:1,-1:1)
2373  complex(dp) :: a1_zpz,a2_zpz,a3_zpz,struc_zpz(3),fac_zpz(-1:1,-1:1)
2374  complex(dp) :: a1_azp,a2_azp,a3_azp,struc_azp(3),fac_azp(-1:1,-1:1)
2375  complex(dp) :: a1_zpa,a2_zpa,a3_zpa,struc_zpa(3),fac_zpa(-1:1,-1:1)
2376  complex(dp) :: a1_zpzp,a2_zpzp,a3_zpzp,struc_zpzp(3),fac_zpzp(-1:1,-1:1)
2377  complex(dp) :: helcoup(1:3,-1:1,-1:1)
2378  complex(dp) :: zab2
2379  complex(dp) :: iprop12,iprop34, zpprop12,zpprop34
2380  complex(dp) :: vvcoupl_prime_zz(4),vvcoupl_prime_az(4),vvcoupl_prime_za(4),vvcoupl_prime_aa(2:4)
2381  complex(dp) :: vvcoupl_prime_zzp(4),vvcoupl_prime_zpz(4),vvcoupl_prime_zpzp(4)
2382  complex(dp) :: vvcoupl_prime_azp(4),vvcoupl_prime_zpa(4)
2383  integer :: vv_it
2384  integer :: i,j,k,l
2385 
2386  zab2(j1,j2,j3,j4) = za(j1,j2)*zb(j2,j4) + za(j1,j3)*zb(j3,j4)
2387 
2388  a0_zz_4f = czero
2389  helcoup(:,:,:)=czero
2390  struc_zz(:)=czero
2391  struc_az(:)=czero
2392  struc_za(:)=czero
2393  struc_aa(:)=czero
2394  struc_zzp(:)=czero
2395  struc_zpz(:)=czero
2396  struc_azp(:)=czero
2397  struc_zpa(:)=czero
2398  struc_zpzp(:)=czero
2399  fac_zz(:,:)=czero
2400  fac_az(:,:)=czero
2401  fac_za(:,:)=czero
2402  fac_aa(:,:)=czero
2403  fac_zzp(:,:)=czero
2404  fac_zpz(:,:)=czero
2405  fac_azp(:,:)=czero
2406  fac_zpa(:,:)=czero
2407  fac_zpzp(:,:)=czero
2408 
2409  lz = (/al_qdn,al_qup,al_qdn,al_qup,al_qdn/)
2410  rz = (/ar_qdn,ar_qup,ar_qdn,ar_qup,ar_qdn/)
2411 
2412  if( abs(isel).eq.pdftop_ .or. abs(jsel).eq.pdftop_ ) return
2413  if( abs(isel).eq.pdfglu_ .or. abs(jsel).eq.pdfglu_ ) return
2414 
2415  line1=abs(isel)
2416  line2=abs(jsel)
2417 
2418  q1q2 = (sprod(j1,j3)+sprod(j1,j4)+sprod(j2,j3)+sprod(j2,j4))/two
2419  mhsq = two * q1q2 + sprod(j1,j2) + sprod(j3,j4)
2420 
2421  kcoupl = q1q2/lambda**2
2422 
2423  q12sq = sprod(j1,j2)
2424  q34sq = sprod(j3,j4)
2425 
2426  iprop12 = q12sq - m_z**2 + ci * m_z * ga_z
2427  iprop34 = q34sq - m_z**2 + ci * m_z * ga_z
2428 
2429  !-- set up couplings
2430  do vv_it=1,4
2431  vvcoupl_prime_zz(vv_it) = hvvspinzerodynamiccoupling(vv_it,q12sq,q34sq,mhsq)
2432  enddo
2433 
2434  a1_zz = vvcoupl_prime_zz(1) * m_z**2/mhsq + vvcoupl_prime_zz(2) * two * q1q2/mhsq + vvcoupl_prime_zz(3) * kcoupl * q1q2/mhsq
2435  a2_zz = -two * vvcoupl_prime_zz(2) - kcoupl * vvcoupl_prime_zz(3)
2436  a3_zz = -two * vvcoupl_prime_zz(4)
2437 
2438  struc_zz(1) = two * (a1_zz * mhsq - ci * a3_zz * q1q2)
2439  struc_zz(2) = (a2_zz + ci * a3_zz)
2440  struc_zz(3) = two * ci * a3_zz
2441  struc_zz(:) = struc_zz(:) * couplzffsq
2442 
2443  fac_zz(-1,-1) = lz(line1) * lz(line2)/iprop12/iprop34
2444  fac_zz(-1,+1) = lz(line1) * rz(line2)/iprop12/iprop34
2445  fac_zz(+1,-1) = rz(line1) * lz(line2)/iprop12/iprop34
2446  fac_zz(+1,+1) = rz(line1) * rz(line2)/iprop12/iprop34
2447 
2448  if( includegammastar ) then
2449 
2450  do vv_it=1,4
2451  vvcoupl_prime_az(vv_it) = hvvspinzerodynamiccoupling(vv_it+4,0d0,q12sq,mhsq)
2452  vvcoupl_prime_za(vv_it) = hvvspinzerodynamiccoupling(vv_it+4,0d0,q34sq,mhsq)
2453  enddo
2454  do vv_it=1,3
2455  vvcoupl_prime_aa(vv_it+1) = hvvspinzerodynamiccoupling(vv_it+8,q12sq,q34sq,mhsq)
2456  enddo
2457 
2458  a1_aa = vvcoupl_prime_aa(2) * two * q1q2/mhsq + vvcoupl_prime_aa(3) * kcoupl * q1q2/mhsq
2459  a2_aa = -two * vvcoupl_prime_aa(2) - kcoupl * vvcoupl_prime_aa(3)
2460  a3_aa = -two * vvcoupl_prime_aa(4)
2461 
2462  a1_za = vvcoupl_prime_za(1) * m_z**2/mhsq + vvcoupl_prime_za(2) * two * q1q2/mhsq + vvcoupl_prime_za(3) * kcoupl * q1q2/mhsq
2463  a2_za = -two * vvcoupl_prime_za(2) - kcoupl * vvcoupl_prime_za(3)
2464  a3_za = -two * vvcoupl_prime_za(4)
2465 
2466  a1_az = vvcoupl_prime_az(1) * m_z**2/mhsq + vvcoupl_prime_az(2) * two * q1q2/mhsq + vvcoupl_prime_az(3) * kcoupl * q1q2/mhsq
2467  a2_az = -two * vvcoupl_prime_az(2) - kcoupl * vvcoupl_prime_az(3)
2468  a3_az = -two * vvcoupl_prime_az(4)
2469 
2470  struc_aa(1) = two * (a1_aa * mhsq - ci * a3_aa * q1q2)
2471  struc_aa(2) = (a2_aa + ci * a3_aa)
2472  struc_aa(3) = two * ci * a3_aa
2473  struc_aa(:) = struc_aa(:) * couplaffsq
2474 
2475  struc_az(1) = two * (a1_az * mhsq - ci * a3_az * q1q2)
2476  struc_az(2) = (a2_az + ci * a3_az)
2477  struc_az(3) = two * ci * a3_az
2478  struc_az(:) = struc_az(:) * couplazff
2479 
2480  struc_za(1) = two * (a1_za * mhsq - ci * a3_za * q1q2)
2481  struc_za(2) = (a2_za + ci * a3_za)
2482  struc_za(3) = two * ci * a3_za
2483  struc_za(:) = struc_za(:) * couplazff
2484 
2485  fac_aa(-1,-1) = la(line1) * la(line2)/q12sq/q34sq
2486  fac_aa(-1,+1) = la(line1) * ra(line2)/q12sq/q34sq
2487  fac_aa(+1,-1) = ra(line1) * la(line2)/q12sq/q34sq
2488  fac_aa(+1,+1) = ra(line1) * ra(line2)/q12sq/q34sq
2489 
2490  fac_az(-1,-1) = la(line1) * lz(line2)/q12sq/iprop34
2491  fac_az(-1,+1) = la(line1) * rz(line2)/q12sq/iprop34
2492  fac_az(+1,-1) = ra(line1) * lz(line2)/q12sq/iprop34
2493  fac_az(+1,+1) = ra(line1) * rz(line2)/q12sq/iprop34
2494 
2495  fac_za(-1,-1) = lz(line1) * la(line2)/iprop12/q34sq
2496  fac_za(-1,+1) = lz(line1) * ra(line2)/iprop12/q34sq
2497  fac_za(+1,-1) = rz(line1) * la(line2)/iprop12/q34sq
2498  fac_za(+1,+1) = rz(line1) * ra(line2)/iprop12/q34sq
2499 
2500  endif
2501 
2502  if ( includevprime ) then
2503 
2504  lct = (/ezp_dn_left,ezp_up_left,ezp_str_left,ezp_chm_left,ezp_bot_left/)
2505  rct = (/ezp_dn_right,ezp_up_right,ezp_str_right,ezp_chm_right,ezp_bot_right/)
2506 
2507  do vv_it=1,4
2508  vvcoupl_prime_zzp(vv_it) = hvvspinzerodynamiccoupling(11+vv_it,q12sq,q34sq,mhsq)
2509  vvcoupl_prime_zpz(vv_it) = hvvspinzerodynamiccoupling(11+vv_it,q34sq,q12sq,mhsq)
2510  vvcoupl_prime_zpzp(vv_it) = hvvspinzerodynamiccoupling(15+vv_it,q12sq,q34sq,mhsq)
2511  enddo
2512 
2513  if (m_zprime.gt.0d0) then
2514  zpprop12 = q12sq - m_zprime**2 + ci * m_zprime * ga_zprime
2515  zpprop34 = q34sq - m_zprime**2 + ci * m_zprime * ga_zprime
2516  elseif (m_zprime.eq.0d0) then
2517  zpprop12 = q12sq
2518  zpprop34 = q34sq
2519  else
2520  zpprop12 = m_z**2
2521  zpprop34 = m_z**2
2522  endif
2523 
2524 
2525  a1_zzp = vvcoupl_prime_zzp(1) * m_z**2/mhsq + vvcoupl_prime_zzp(2) * two * q1q2/mhsq + vvcoupl_prime_zzp(3) * kcoupl * q1q2/mhsq
2526  a2_zzp = -two * vvcoupl_prime_zzp(2) - kcoupl * vvcoupl_prime_zzp(3)
2527  a3_zzp = -two * vvcoupl_prime_zzp(4)
2528  a1_zpz = vvcoupl_prime_zpz(1) * m_z**2/mhsq + vvcoupl_prime_zpz(2) * two * q1q2/mhsq + vvcoupl_prime_zpz(3) * kcoupl * q1q2/mhsq
2529  a2_zpz = -two * vvcoupl_prime_zpz(2) - kcoupl * vvcoupl_prime_zpz(3)
2530  a3_zpz = -two * vvcoupl_prime_zpz(4)
2531  a1_zpzp = vvcoupl_prime_zpzp(1) * m_z**2/mhsq + vvcoupl_prime_zpzp(2) * two * q1q2/mhsq + vvcoupl_prime_zpzp(3) * kcoupl * q1q2/mhsq
2532  a2_zpzp = -two * vvcoupl_prime_zpzp(2) - kcoupl * vvcoupl_prime_zpzp(3)
2533  a3_zpzp = -two * vvcoupl_prime_zpzp(4)
2534 
2535  struc_zzp(1) = two * (a1_zzp * mhsq - ci * a3_zzp * q1q2)
2536  struc_zzp(2) = (a2_zzp + ci * a3_zzp)
2537  struc_zzp(3) = two * ci * a3_zzp
2538  struc_zpz(1) = two * (a1_zpz * mhsq - ci * a3_zpz * q1q2)
2539  struc_zpz(2) = (a2_zpz + ci * a3_zpz)
2540  struc_zpz(3) = two * ci * a3_zpz
2541  struc_zpzp(1) = two * (a1_zpzp * mhsq - ci * a3_zpzp * q1q2)
2542  struc_zpzp(2) = (a2_zpzp + ci * a3_zpzp)
2543  struc_zpzp(3) = two * ci * a3_zpzp
2544 
2545  struc_zzp(:) = struc_zzp(:) * couplzffsq
2546  struc_zpz(:) = struc_zpz(:) * couplzffsq
2547  struc_zpzp(:) = struc_zpzp(:) * couplzffsq
2548 
2549  fac_zpzp(-1,-1) = lct(line1) * lct(line2)/zpprop12/zpprop34
2550  fac_zpzp(-1,+1) = lct(line1) * rct(line2)/zpprop12/zpprop34
2551  fac_zpzp(+1,-1) = rct(line1) * lct(line2)/zpprop12/zpprop34
2552  fac_zpzp(+1,+1) = rct(line1) * rct(line2)/zpprop12/zpprop34
2553 
2554  fac_zpz(-1,-1) = lct(line1) * lz(line2)/zpprop12/iprop34
2555  fac_zpz(-1,+1) = lct(line1) * rz(line2)/zpprop12/iprop34
2556  fac_zpz(+1,-1) = rct(line1) * lz(line2)/zpprop12/iprop34
2557  fac_zpz(+1,+1) = rct(line1) * rz(line2)/zpprop12/iprop34
2558 
2559  fac_zzp(-1,-1) = lz(line1) * lct(line2)/iprop12/zpprop34
2560  fac_zzp(-1,+1) = lz(line1) * rct(line2)/iprop12/zpprop34
2561  fac_zzp(+1,-1) = rz(line1) * lct(line2)/iprop12/zpprop34
2562  fac_zzp(+1,+1) = rz(line1) * rct(line2)/iprop12/zpprop34
2563 
2564 
2565  if( includegammastar ) then
2566  do vv_it=1,4
2567  vvcoupl_prime_azp(vv_it) = hvvspinzerodynamiccoupling(19+vv_it,q12sq,q34sq,mhsq)
2568  vvcoupl_prime_zpa(vv_it) = hvvspinzerodynamiccoupling(19+vv_it,q34sq,q12sq,mhsq)
2569  enddo
2570 
2571  a1_azp = vvcoupl_prime_azp(1) * m_z**2/mhsq + vvcoupl_prime_azp(2) * two * q1q2/mhsq + vvcoupl_prime_azp(3) * kcoupl * q1q2/mhsq
2572  a2_azp = -two * vvcoupl_prime_azp(2) - kcoupl * vvcoupl_prime_azp(3)
2573  a3_azp = -two * vvcoupl_prime_azp(4)
2574 
2575  a1_zpa = vvcoupl_prime_zpa(1) * m_z**2/mhsq + vvcoupl_prime_zpa(2) * two * q1q2/mhsq + vvcoupl_prime_zpa(3) * kcoupl * q1q2/mhsq
2576  a2_zpa = -two * vvcoupl_prime_zpa(2) - kcoupl * vvcoupl_prime_zpa(3)
2577  a3_zpa = -two * vvcoupl_prime_zpa(4)
2578 
2579  struc_azp(1) = two * (a1_azp * mhsq - ci * a3_azp * q1q2)
2580  struc_azp(2) = (a2_azp + ci * a3_azp)
2581  struc_azp(3) = two * ci * a3_azp
2582  struc_azp(:) = struc_azp(:) * couplazff
2583 
2584  struc_zpa(1) = two * (a1_zpa * mhsq - ci * a3_zpa * q1q2)
2585  struc_zpa(2) = (a2_zpa + ci * a3_zpa)
2586  struc_zpa(3) = two * ci * a3_zpa
2587  struc_zpa(:) = struc_zpa(:) * couplazff
2588 
2589  fac_zpa(-1,-1) = lct(line1) * la(line2)/zpprop12/q34sq
2590  fac_zpa(-1,+1) = lct(line1) * ra(line2)/zpprop12/q34sq
2591  fac_zpa(+1,-1) = rct(line1) * la(line2)/zpprop12/q34sq
2592  fac_zpa(+1,+1) = rct(line1) * ra(line2)/zpprop12/q34sq
2593 
2594  fac_azp(-1,-1) = la(line1) * lct(line2)/q12sq/zpprop34
2595  fac_azp(-1,+1) = la(line1) * rct(line2)/q12sq/zpprop34
2596  fac_azp(+1,-1) = ra(line1) * lct(line2)/q12sq/zpprop34
2597  fac_azp(+1,+1) = ra(line1) * rct(line2)/q12sq/zpprop34
2598  endif
2599 
2600  endif
2601 
2602  do i=-1,1,2
2603  do j=-1,1,2
2604  helcoup(1:3,i,j) = struc_zz(1:3) * fac_zz(i,j) + &
2605  struc_az(1:3) * fac_az(i,j) + &
2606  struc_za(1:3) * fac_za(i,j) + &
2607  struc_zpz(1:3) * fac_zpz(i,j) + &
2608  struc_zzp(1:3) * fac_zzp(i,j) + &
2609  struc_zpa(1:3) * fac_zpa(i,j) + &
2610  struc_azp(1:3) * fac_azp(i,j) + &
2611  struc_aa(1:3) * fac_aa(i,j) + &
2612  struc_zpzp(1:3) * fac_zpzp(i,j)
2613  enddo
2614  enddo
2615 
2616 
2617  a0_zz_4f(-1,-1) = za(j1,j3)*zb(j4,j2) * helcoup(1,-1,-1) + &
2618  zab2(j1,j3,j4,j2)*zab2(j3,j1,j2,j4) * helcoup(2,-1,-1) + &
2619  za(j1,j2)*za(j3,j4)*zb(j4,j2)**2 * helcoup(3,-1,-1)
2620 
2621  a0_zz_4f(-1,+1) = za(j1,j4)*zb(j3,j2) * helcoup(1,-1,+1) + &
2622  zab2(j1,j3,j4,j2)*zab2(j4,j1,j2,j3) * helcoup(2,-1,+1) + &
2623  za(j1,j2)*za(j4,j3)*zb(j3,j2)**2 * helcoup(3,-1,+1)
2624 
2625  a0_zz_4f(+1,-1) = za(j2,j3)*zb(j4,j1) * helcoup(1,+1,-1) + &
2626  zab2(j2,j3,j4,j1)*zab2(j3,j1,j2,j4) * helcoup(2,+1,-1) + &
2627  za(j2,j1)*za(j3,j4)*zb(j4,j1)**2 * helcoup(3,+1,-1)
2628 
2629  a0_zz_4f(+1,+1) = za(j2,j4)*zb(j3,j1) * helcoup(1,+1,+1) + &
2630  zab2(j2,j3,j4,j1)*zab2(j4,j1,j2,j3) * helcoup(2,+1,+1) + &
2631  za(j2,j1)*za(j4,j3)*zb(j3,j1)**2 * helcoup(3,+1,+1)
2632 
2633  a0_zz_4f = a0_zz_4f/vev
2634 
2635  return
2636 
2637  end function a0_zz_4f
2638 
2639  function a0_ww_4f(j1,j2,j3,j4,za,zb,sprod,iSel,jSel,useWWcoupl,Wpm_flip)
2640  use modmisc
2641  implicit none
2642  complex(dp) :: a0_ww_4f(-1:1,-1:1)
2643  integer :: j1,j2,j3,j4,isel,jsel,line1,line2
2644  complex(dp) :: za(4,4), zb(4,4)
2645  logical,optional :: usewwcoupl,wpm_flip
2646  real(dp) :: sprod(4,4),q2wplus,q2wminus
2647  real(dp) :: mhsq, q1q2, kcoupl
2648  complex(dp) :: vvcoupl_prime_ww(4)
2649  complex(dp) :: a1_ww, a2_ww, a3_ww, struc_ww(1:3)
2650  complex(dp) :: vvcoupl_prime_wwp(4)
2651  complex(dp) :: a1_wwp, a2_wwp, a3_wwp, struc_wwp(1:3)
2652  complex(dp) :: vvcoupl_prime_wpw(4)
2653  complex(dp) :: a1_wpw, a2_wpw, a3_wpw, struc_wpw(1:3)
2654  complex(dp) :: vvcoupl_prime_wpwp(4)
2655  complex(dp) :: a1_wpwp, a2_wpwp, a3_wpwp, struc_wpwp(1:3)
2656  complex(dp) :: zab2
2657  complex(dp) :: iprop12, iprop34, wpprop12, wpprop34
2658  complex(dp) :: helcoup(1:3,-1:1,-1:1)
2659  real(dp), dimension(5) :: lct,rct
2660  integer :: vv_it
2661  integer :: i,j,k,l
2662 
2663  zab2(j1,j2,j3,j4) = za(j1,j2)*zb(j2,j4) + za(j1,j3)*zb(j3,j4)
2664 
2665  a0_ww_4f = czero
2666  helcoup(:,:,:)=czero
2667 
2668  if( abs(isel).eq.pdftop_ .or. abs(jsel).eq.pdftop_ ) return
2669  if( abs(isel).eq.pdfglu_ .or. abs(jsel).eq.pdfglu_ ) return
2670 
2671  line1=abs(isel)
2672  line2=abs(jsel)
2673 
2674  q1q2 = (sprod(j1,j3)+sprod(j1,j4)+sprod(j2,j3)+sprod(j2,j4))/two
2675  mhsq = two * q1q2 + sprod(j1,j2) + sprod(j3,j4)
2676 
2677  kcoupl = q1q2/lambda**2
2678 
2679  q2wplus = sprod(j1,j2)
2680  q2wminus = sprod(j3,j4)
2681  if( present(wpm_flip) ) then
2682  if( wpm_flip ) then
2683  call swap(q2wplus,q2wminus)
2684  call swap(line1,line2)
2685  endif
2686  endif
2687  if( .not.present(usewwcoupl) ) then
2688  do vv_it=1,4
2689  vvcoupl_prime_ww(vv_it) = hvvspinzerodynamiccoupling(vv_it,q2wplus,q2wminus,mhsq)
2690  enddo
2691  else
2692  do vv_it=1,4
2693  vvcoupl_prime_ww(vv_it) = hvvspinzerodynamiccoupling(vv_it,q2wplus,q2wminus,mhsq,trywwcoupl=usewwcoupl)
2694  enddo
2695  endif
2696 
2697  ! SM contribution
2698  iprop12 = q2wplus - m_w**2 + ci * m_w * ga_w
2699  iprop34 = q2wminus - m_w**2 + ci * m_w * ga_w
2700 
2701  a1_ww = vvcoupl_prime_ww(1) * m_w**2/mhsq + vvcoupl_prime_ww(2) * two * q1q2/mhsq + vvcoupl_prime_ww(3) * kcoupl * q1q2/mhsq
2702  a2_ww = -two * vvcoupl_prime_ww(2) - kcoupl * vvcoupl_prime_ww(3)
2703  a3_ww = -two * vvcoupl_prime_ww(4)
2704 
2705  struc_ww(1) = two * (a1_ww * mhsq - ci * a3_ww * q1q2)
2706  struc_ww(2) = a2_ww + ci * a3_ww
2707  struc_ww(3) = two * ci * a3_ww
2708  struc_ww(:) = struc_ww(:) * couplwffsq
2709 
2710  helcoup(1:3,-1,-1) = struc_ww(1:3)/iprop12/iprop34
2711 
2712  if ( includevprime ) then
2713 
2714  lct = (/ewp_up_left,ewp_up_left,ewp_chm_left,ewp_chm_left,ewp_top_left/)
2715  rct = (/ewp_up_right,ewp_up_right,ewp_chm_right,ewp_chm_right,ewp_top_right/)
2716  ! Adjust couplings convention
2717  lct(:) = lct(:)/bl
2718  rct(:) = rct(:)/bl
2719 
2720  do vv_it=1,4
2721  vvcoupl_prime_wwp(vv_it) = hvvspinzerodynamiccoupling(11+vv_it,q2wplus,q2wminus,mhsq,trywwcoupl=usewwcoupl)
2722  vvcoupl_prime_wpw(vv_it) = hvvspinzerodynamiccoupling(11+vv_it,q2wminus,q2wplus,mhsq,trywwcoupl=usewwcoupl)
2723  vvcoupl_prime_wpwp(vv_it) = hvvspinzerodynamiccoupling(15+vv_it,q2wplus,q2wminus,mhsq,trywwcoupl=usewwcoupl)
2724  enddo
2725 
2726 
2727  if (m_wprime.gt.0d0) then
2728  wpprop12 = q2wplus - m_wprime**2 + ci * m_wprime * ga_wprime
2729  wpprop34 = q2wminus - m_wprime**2 + ci * m_wprime * ga_wprime
2730  elseif (m_wprime.eq.0d0) then
2731  wpprop12 = q2wplus
2732  wpprop34 = q2wminus
2733  else
2734  wpprop12 = m_w**2
2735  wpprop34 = m_w**2
2736  endif
2737 
2738 
2739  a1_wwp = vvcoupl_prime_wwp(1) * m_w**2/mhsq + vvcoupl_prime_wwp(2) * two * q1q2/mhsq + vvcoupl_prime_wwp(3) * kcoupl * q1q2/mhsq
2740  a2_wwp = -two * vvcoupl_prime_wwp(2) - kcoupl * vvcoupl_prime_wwp(3)
2741  a3_wwp = -two * vvcoupl_prime_wwp(4)
2742  a1_wpw = vvcoupl_prime_wpw(1) * m_w**2/mhsq + vvcoupl_prime_wpw(2) * two * q1q2/mhsq + vvcoupl_prime_wpw(3) * kcoupl * q1q2/mhsq
2743  a2_wpw = -two * vvcoupl_prime_wpw(2) - kcoupl * vvcoupl_prime_wpw(3)
2744  a3_wpw = -two * vvcoupl_prime_wpw(4)
2745  a1_wpwp = vvcoupl_prime_wpwp(1) * m_w**2/mhsq + vvcoupl_prime_wpwp(2) * two * q1q2/mhsq + vvcoupl_prime_wpwp(3) * kcoupl * q1q2/mhsq
2746  a2_wpwp = -two * vvcoupl_prime_wpwp(2) - kcoupl * vvcoupl_prime_wpwp(3)
2747  a3_wpwp = -two * vvcoupl_prime_wpwp(4)
2748 
2749  struc_wwp(1) = two * (a1_wwp * mhsq - ci * a3_wwp * q1q2)
2750  struc_wwp(2) = (a2_wwp + ci * a3_wwp)
2751  struc_wwp(3) = two * ci * a3_wwp
2752  struc_wpw(1) = two * (a1_wpw * mhsq - ci * a3_wpw * q1q2)
2753  struc_wpw(2) = (a2_wpw + ci * a3_wpw)
2754  struc_wpw(3) = two * ci * a3_wpw
2755  struc_wpwp(1) = two * (a1_wpwp * mhsq - ci * a3_wpwp * q1q2)
2756  struc_wpwp(2) = (a2_wpwp + ci * a3_wpwp)
2757  struc_wpwp(3) = two * ci * a3_wpwp
2758 
2759  struc_wwp(:) = struc_wwp(:) * couplwffsq
2760  struc_wpw(:) = struc_wpw(:) * couplwffsq
2761  struc_wpwp(:) = struc_wpwp(:) * couplwffsq
2762 
2763  helcoup(1:3,-1,-1) = helcoup(1:3,-1,-1) + &
2764  struc_wpwp(1:3) * lct(line1) * lct(line2)/wpprop12/wpprop34 + &
2765  struc_wpw(1:3) * lct(line1) /wpprop12/iprop34 + &
2766  struc_wwp(1:3) * lct(line2)/iprop12/wpprop34
2767 
2768  helcoup(1:3,-1,+1) = helcoup(1:3,-1,+1) + &
2769  struc_wpwp(1:3) * lct(line1) * rct(line2)/wpprop12/wpprop34 + &
2770  struc_wpw(1:3) * lct(line1) * czero /wpprop12/iprop34 + &
2771  struc_wwp(1:3) * rct(line2)/iprop12/wpprop34
2772 
2773  helcoup(1:3,+1,-1) = helcoup(1:3,+1,-1) + &
2774  struc_wpwp(1:3) * rct(line1) * lct(line2)/wpprop12/wpprop34 + &
2775  struc_wpw(1:3) * rct(line1) /wpprop12/iprop34 + &
2776  struc_wwp(1:3) * czero * lct(line2)/iprop12/wpprop34
2777 
2778  helcoup(1:3,+1,+1) = helcoup(1:3,+1,+1) + &
2779  struc_wpwp(1:3) * rct(line1) * rct(line2)/wpprop12/wpprop34 + &
2780  struc_wpw(1:3) * rct(line1) * czero /wpprop12/iprop34 + &
2781  struc_wwp(1:3) * czero * rct(line2)/iprop12/wpprop34
2782  endif
2783 
2784 
2785  a0_ww_4f(-1,-1) = za(j1,j3)*zb(j4,j2) * helcoup(1,-1,-1) + &
2786  zab2(j1,j3,j4,j2)*zab2(j3,j1,j2,j4) * helcoup(2,-1,-1) + &
2787  za(j1,j2)*za(j3,j4)*zb(j4,j2)**2 * helcoup(3,-1,-1)
2788 
2789  if ( includevprime ) then
2790  a0_ww_4f(-1,+1) = za(j1,j4)*zb(j3,j2) * helcoup(1,-1,+1) + &
2791  zab2(j1,j3,j4,j2)*zab2(j4,j1,j2,j3) * helcoup(2,-1,+1) + &
2792  za(j1,j2)*za(j4,j3)*zb(j3,j2)**2 * helcoup(3,-1,+1)
2793 
2794  a0_ww_4f(+1,-1) = za(j2,j3)*zb(j4,j1) * helcoup(1,+1,-1) + &
2795  zab2(j2,j3,j4,j1)*zab2(j3,j1,j2,j4) * helcoup(2,+1,-1) + &
2796  za(j2,j1)*za(j3,j4)*zb(j4,j1)**2 * helcoup(3,+1,-1)
2797 
2798  a0_ww_4f(+1,+1) = za(j2,j4)*zb(j3,j1) * helcoup(1,+1,+1) + &
2799  zab2(j2,j3,j4,j1)*zab2(j4,j1,j2,j3) * helcoup(2,+1,+1) + &
2800  za(j2,j1)*za(j4,j3)*zb(j3,j1)**2 * helcoup(3,+1,+1)
2801  endif
2802 
2803 
2804  a0_ww_4f = a0_ww_4f/vev
2805 
2806 
2807  return
2808 
2809  end function a0_ww_4f
2810 
2811 
2812  !-- QCD amplitudes squared below
2813  subroutine me2_ggggh(j1,j2,j3,j4,za,zb,sprod,res)
2814  integer, intent(in) :: j1,j2,j3,j4
2815  complex(dp), intent(in) :: za(4,4), zb(4,4)
2816  real(dp), intent(in) :: sprod(4,4)
2817  real(dp), intent(out) :: res
2818  complex(dp) :: a1234(-1:1,-1:1,-1:1,-1:1), a1324(-1:1,-1:1,-1:1,-1:1)
2819  complex(dp) :: aphi1234(1:2,-1:1,-1:1,-1:1,-1:1), aphi1324(1:2,-1:1,-1:1,-1:1,-1:1)
2820  real(dp), parameter :: col = 8.0_dp * ca**3 * cf
2821  integer :: i1,i2,i3,i4
2822  complex(dp) :: scalar, pseudo
2823 
2824  res = zero
2825 
2826  scalar = ghg2
2827  pseudo = -ghg4
2828 
2829  aphi1234 = a0phigggg_xxxx(j1,j2,j3,j4,za,zb,sprod)
2830  aphi1324 = a0phigggg_xxxx(j1,j3,j2,j4,za,zb,sprod)
2831 
2832  a1234(:,:,:,:) = scalar * (aphi1234(1,:,:,:,:) + aphi1234(2,:,:,:,:)) + &
2833  pseudo * (-ci) * (aphi1234(1,:,:,:,:) - aphi1234(2,:,:,:,:))
2834 
2835  a1324(:,:,:,:) = scalar * (aphi1324(1,:,:,:,:) + aphi1324(2,:,:,:,:)) + &
2836  pseudo * (-ci) * (aphi1324(1,:,:,:,:) - aphi1324(2,:,:,:,:))
2837 
2838  do i1 = -1,1,2
2839  do i2 = -1,1,2
2840  do i3 = -1,1,2
2841  do i4 = -1,1,2
2842 
2843  res = res + real(a1234(i1,i2,i3,i4)*conjg(a1234(i1,i2,i3,i4)),kind=dp)
2844  res = res + real(a1324(i1,i2,i3,i4)*conjg(a1324(i1,i2,i3,i4)),kind=dp)
2845  res = res + real(a1234(i1,i2,i3,i4)*conjg(a1324(i1,i3,i2,i4)),kind=dp)
2846 
2847  enddo
2848  enddo
2849  enddo
2850  enddo
2851 
2852  res = res * col / vev**2
2853 
2854  return
2855 
2856  end subroutine me2_ggggh
2857 
2858  subroutine me2_qbqggh(j1,j2,j3,j4,za,zb,sprod,res)
2859  integer, intent(in) :: j1,j2,j3,j4
2860  complex(dp), intent(in) :: za(4,4), zb(4,4)
2861  real(dp), intent(in) :: sprod(4,4)
2862  real(dp), intent(out) :: res
2863  complex(dp) :: a1234(-1:1,-1:1,-1:1), a1243(-1:1,-1:1,-1:1)
2864  complex(dp) :: aphi1234(1:2,-1:1,-1:1,-1:1), aphi1243(1:2,-1:1,-1:1,-1:1)
2865  real(dp), parameter :: col1 = 4.0_dp * xn * cf**2
2866  real(dp), parameter :: col2 = 2.0_dp * xn * cf * (2.0_dp * cf - ca)
2867  integer :: i12,i3,i4
2868  complex(dp) :: scalar, pseudo
2869 
2870  res = zero
2871 
2872  scalar = ghg2
2873  pseudo = -ghg4
2874 
2875  aphi1234 = a0phiqbqgg_xxx(j1,j2,j3,j4,za,zb,sprod)
2876  aphi1243 = a0phiqbqgg_xxx(j1,j2,j4,j3,za,zb,sprod)
2877 
2878  a1234(:,:,:) = scalar * (aphi1234(1,:,:,:) + aphi1234(2,:,:,:)) + &
2879  pseudo * (-ci) * (aphi1234(1,:,:,:) - aphi1234(2,:,:,:))
2880 
2881  a1243(:,:,:) = scalar * (aphi1243(1,:,:,:) + aphi1243(2,:,:,:)) + &
2882  pseudo * (-ci) * (aphi1243(1,:,:,:) - aphi1243(2,:,:,:))
2883 
2884  do i12 = -1,1,2
2885  do i3 = -1,1,2
2886  do i4 = -1,1,2
2887 
2888  res = res + real(a1234(i12,i3,i4)*conjg(a1234(i12,i3,i4)),kind=dp) * col1
2889  res = res + real(a1243(i12,i4,i3)*conjg(a1243(i12,i4,i3)),kind=dp) * col1
2890  res = res + two * real(a1234(i12,i3,i4)*conjg(a1243(i12,i4,i3)),kind=dp) * col2
2891 
2892  enddo
2893  enddo
2894  enddo
2895 
2896  res = res / vev**2
2897 
2898  return
2899 
2900  end subroutine me2_qbqggh
2901 
2902  subroutine me2_qbqqbq(j1,j2,j3,j4,za,zb,sprod,res_diff,res_id)
2903  integer, intent(in) :: j1,j2,j3,j4
2904  complex(dp), intent(in) :: za(4,4), zb(4,4)
2905  real(dp), intent(in) :: sprod(4,4)
2906  real(dp), intent(out) :: res_diff, res_id
2907  complex(dp) :: amp_a(-1:1,-1:1), amp_b(-1:1,-1:1)
2908  complex(dp) :: aphi_a(1:2,-1:1,-1:1), aphi_b(1:2,-1:1,-1:1)
2909  real(dp), parameter :: col = xn**2 - one
2910  integer :: h12,h34
2911  complex(dp) :: scalar, pseudo
2912 
2913  res_diff = zero
2914  res_id = zero
2915 
2916  scalar = ghg2
2917  pseudo = -ghg4
2918 
2919  aphi_a = a0phiqbqqbq_xx(j1,j2,j3,j4,za,zb,sprod)
2920  aphi_b = -a0phiqbqqbq_xx(j1,j4,j3,j2,za,zb,sprod)
2921 
2922  amp_a(:,:) = scalar * (aphi_a(1,:,:) + aphi_a(2,:,:)) + &
2923  pseudo * (-ci) * (aphi_a(1,:,:) - aphi_a(2,:,:))
2924 
2925  amp_b(:,:) = scalar * (aphi_b(1,:,:) + aphi_b(2,:,:)) + &
2926  pseudo * (-ci) * (aphi_b(1,:,:) - aphi_b(2,:,:))
2927 
2928  do h12 = -1,1,2
2929  do h34 = -1,1,2
2930  res_diff = res_diff + real(amp_a(h12,h34)*conjg(amp_a(h12,h34)),kind=dp)
2931  res_id = res_id + real(amp_a(h12,h34)*conjg(amp_a(h12,h34)),kind=dp)
2932  res_id = res_id + real(amp_b(h12,h34)*conjg(amp_b(h12,h34)),kind=dp)
2933  enddo
2934  res_id = res_id - two/xn * real(amp_a(h12,h12)*conjg(amp_b(h12,h12)),kind=dp)
2935  enddo
2936 
2937  res_id = res_id * col / vev**2
2938  res_diff = res_diff * col / vev**2
2939 
2940  return
2941 
2942  end subroutine me2_qbqqbq
2943 
2944 
2945  function a0phigggg_xxxx(j1,j2,j3,j4,za,zb,sprod)
2946  complex(dp) :: a0phigggg_xxxx(1:2,-1:1,-1:1,-1:1,-1:1)
2947  integer :: j1,j2,j3,j4
2948  complex(dp) :: za(4,4),zb(4,4)
2949  real(dp) :: sprod(4,4)
2950 
2951  a0phigggg_xxxx = czero
2952 
2953  a0phigggg_xxxx(1,+1,+1,-1,-1) = a0phiggggmmpp(j3,j4,j1,j2,za,zb,sprod)
2954  a0phigggg_xxxx(1,+1,-1,+1,-1) = a0phiggggmpmp(j2,j3,j4,j1,za,zb,sprod)
2955  a0phigggg_xxxx(1,+1,-1,-1,+1) = a0phiggggmmpp(j2,j3,j4,j1,za,zb,sprod)
2956  a0phigggg_xxxx(1,+1,-1,-1,-1) = a0phiggggpmmm(j1,j2,j3,j4,za,zb,sprod)
2957  a0phigggg_xxxx(1,-1,+1,+1,-1) = a0phiggggmmpp(j4,j1,j2,j3,za,zb,sprod)
2958  a0phigggg_xxxx(1,-1,+1,-1,+1) = a0phiggggmpmp(j1,j2,j3,j4,za,zb,sprod)
2959  a0phigggg_xxxx(1,-1,+1,-1,-1) = a0phiggggpmmm(j2,j3,j4,j1,za,zb,sprod)
2960  a0phigggg_xxxx(1,-1,-1,+1,+1) = a0phiggggmmpp(j1,j2,j3,j4,za,zb,sprod)
2961  a0phigggg_xxxx(1,-1,-1,+1,-1) = a0phiggggpmmm(j3,j4,j1,j2,za,zb,sprod)
2962  a0phigggg_xxxx(1,-1,-1,-1,+1) = a0phiggggpmmm(j4,j1,j2,j3,za,zb,sprod)
2963  a0phigggg_xxxx(1,-1,-1,-1,-1) = a0phiggggmmmm(j1,j2,j3,j4,za,zb,sprod)
2964 
2965  a0phigggg_xxxx(2,-1,-1,+1,+1) = conjg(a0phigggg_xxxx(1,+1,+1,-1,-1))
2966  a0phigggg_xxxx(2,-1,+1,-1,+1) = conjg(a0phigggg_xxxx(1,+1,-1,+1,-1))
2967  a0phigggg_xxxx(2,-1,+1,+1,-1) = conjg(a0phigggg_xxxx(1,+1,-1,-1,+1))
2968  a0phigggg_xxxx(2,-1,+1,+1,+1) = conjg(a0phigggg_xxxx(1,+1,-1,-1,-1))
2969  a0phigggg_xxxx(2,+1,-1,-1,+1) = conjg(a0phigggg_xxxx(1,-1,+1,+1,-1))
2970  a0phigggg_xxxx(2,+1,-1,+1,-1) = conjg(a0phigggg_xxxx(1,-1,+1,-1,+1))
2971  a0phigggg_xxxx(2,+1,-1,+1,+1) = conjg(a0phigggg_xxxx(1,-1,+1,-1,-1))
2972  a0phigggg_xxxx(2,+1,+1,-1,-1) = conjg(a0phigggg_xxxx(1,-1,-1,+1,+1))
2973  a0phigggg_xxxx(2,+1,+1,-1,+1) = conjg(a0phigggg_xxxx(1,-1,-1,+1,-1))
2974  a0phigggg_xxxx(2,+1,+1,+1,-1) = conjg(a0phigggg_xxxx(1,-1,-1,-1,+1))
2975  a0phigggg_xxxx(2,+1,+1,+1,+1) = conjg(a0phigggg_xxxx(1,-1,-1,-1,-1))
2976 
2977  return
2978 
2979  end function a0phigggg_xxxx
2980 
2981  function a0phiqbqgg_xxx(j1,j2,j3,j4,za,zb,sprod)
2982  complex(dp) :: a0phiqbqgg_xxx(1:2,-1:1,-1:1,-1:1)
2983  integer :: j1,j2,j3,j4
2984  complex(dp) :: za(4,4), zb(4,4)
2985  real(dp) :: sprod(4,4)
2986  real(dp) :: s3
2987  complex(dp) :: zab2
2988 
2989  s3(j1,j2,j3) = sprod(j1,j2) + sprod(j1,j3) + sprod(j2,j3)
2990  zab2(j1,j2,j3,j4) = za(j1,j2)*zb(j2,j4) + za(j1,j3)*zb(j3,j4)
2991 
2992  a0phiqbqgg_xxx = czero
2993 
2994  a0phiqbqgg_xxx(1,-1,+1,-1) = -za(j1,j4)**2*za(j2,j4)/(za(j1,j2)*za(j2,j3)*za(j3,j4))
2995  a0phiqbqgg_xxx(1,-1,-1,+1) = za(j1,j3)**3/(za(j1,j2)*za(j3,j4)*za(j4,j1))
2996  a0phiqbqgg_xxx(1,-1,-1,-1) = -((za(j1,j3)*zab2(j4,j1,j3,j2)**2)/ &
2997  (s3(j1,j2,j3)*sprod(j1,j2)*zb(j2,j3))) - &
2998  ((1/sprod(j1,j2) + 1/sprod(j4,j1))*za(j4,j1)*zab2(j3,j1,j4,j2)**2)/ &
2999  (s3(j4,j1,j2)*zb(j2,j4)) + &
3000  zab2(j1,j3,j4,j2)**2/(za(j1,j2)*zb(j2,j3)*zb(j2,j4)*zb(j3,j4))
3001 
3002  a0phiqbqgg_xxx(2,-1,+1,-1) = -zb(j1,j3)*zb(j2,j3)**2/(zb(j1,j2)*zb(j3,j4)*zb(j4,j1))
3003  a0phiqbqgg_xxx(2,-1,-1,+1) = zb(j2,j4)**3/(zb(j1,j2)*zb(j2,j3)*zb(j3,j4))
3004  a0phiqbqgg_xxx(2,-1,+1,+1) = -(zab2(j1,j3,j4,j2)**2/(za(j1,j3)*za(j1,j4)*za(j3,j4)*zb(j1,j2))) - &
3005  ((1/sprod(j1,j2) + 1/sprod(j2,j3))*zab2(j1,j2,j3,j4)**2*zb(j2,j3))/ &
3006  (s3(j1,j2,j3)*za(j1,j3)) + &
3007  (zab2(j1,j2,j4,j3)**2*zb(j2,j4))/ &
3008  (s3(j4,j1,j2)*sprod(j1,j2)*za(j1,j4))
3009 
3010  a0phiqbqgg_xxx(1,+1,+1,-1) = conjg(a0phiqbqgg_xxx(2,-1,-1,+1))
3011  a0phiqbqgg_xxx(1,+1,-1,+1) = conjg(a0phiqbqgg_xxx(2,-1,+1,-1))
3012  a0phiqbqgg_xxx(1,+1,-1,-1) = conjg(a0phiqbqgg_xxx(2,-1,+1,+1))
3013 
3014  a0phiqbqgg_xxx(2,+1,+1,-1) = conjg(a0phiqbqgg_xxx(1,-1,-1,+1))
3015  a0phiqbqgg_xxx(2,+1,-1,+1) = conjg(a0phiqbqgg_xxx(1,-1,+1,-1))
3016  a0phiqbqgg_xxx(2,+1,+1,+1) = conjg(a0phiqbqgg_xxx(1,-1,-1,-1))
3017 
3018  return
3019 
3020  end function a0phiqbqgg_xxx
3021 
3022  function a0phiqbqqbq_xx(j1,j2,j3,j4,za,zb,sprod)
3023  complex(dp) :: a0phiqbqqbq_xx(1:2,-1:1,-1:1)
3024  integer :: j1,j2,j3,j4
3025  complex(dp) :: za(4,4), zb(4,4)
3026  real(dp) :: sprod(4,4)
3027 
3028  a0phiqbqqbq_xx = czero
3029 
3030  a0phiqbqqbq_xx(1,-1,+1) = za(j1,j4)**2/(za(j1,j2)*za(j3,j4))
3031  a0phiqbqqbq_xx(1,-1,-1) = za(j1,j3)**2/(za(j1,j2)*za(j4,j3))
3032 
3033  a0phiqbqqbq_xx(2,-1,+1) = zb(j2,j3)**2/(zb(j1,j2)*zb(j3,j4))
3034  a0phiqbqqbq_xx(2,-1,-1) = zb(j2,j4)**2/(zb(j1,j2)*zb(j4,j3))
3035 
3036  a0phiqbqqbq_xx(1,+1,-1) = conjg(a0phiqbqqbq_xx(2,-1,+1))
3037  a0phiqbqqbq_xx(1,+1,+1) = conjg(a0phiqbqqbq_xx(2,-1,-1))
3038 
3039  a0phiqbqqbq_xx(2,+1,-1) = conjg(a0phiqbqqbq_xx(1,-1,+1))
3040  a0phiqbqqbq_xx(2,+1,+1) = conjg(a0phiqbqqbq_xx(1,-1,-1))
3041 
3042  return
3043 
3044  end function a0phiqbqqbq_xx
3045 
3046 
3047  function a0phiggggpmmm(j1,j2,j3,j4,za,zb,sprod)
3048  integer :: j1,j2,j3,j4
3049  complex(dp) :: za(4,4), zb(4,4)
3050  real(dp) :: sprod(4,4)
3051  complex(dp) :: a0phiggggpmmm
3052 
3053  real(dp) :: s3
3054  complex(dp) :: zab2
3055 
3056  s3(j1,j2,j3)=sprod(j1,j2)+sprod(j2,j3)+sprod(j3,j1)
3057  zab2(j1,j2,j3,j4)=za(j1,j2)*zb(j2,j4)+za(j1,j3)*zb(j3,j4)
3058 
3059  a0phiggggpmmm = &
3060  +(zab2(j3,j2,j4,j1)*za(j2,j4))**2/(s3(j1,j2,j4)*sprod(j1,j2)*sprod(j1,j4)) &
3061  +(zab2(j4,j2,j3,j1)*za(j2,j3))**2/(s3(j1,j2,j3)*sprod(j1,j2)*sprod(j2,j3)) &
3062  +(zab2(j2,j3,j4,j1)*za(j3,j4))**2/(s3(j1,j3,j4)*sprod(j1,j4)*sprod(j3,j4)) &
3063  -za(j2,j4)/(za(j1,j2)*zb(j2,j3)*zb(j3,j4)*za(j4,j1)) &
3064  *(-sprod(j2,j3)*zab2(j2,j3,j4,j1)/zb(j4,j1) &
3065  -sprod(j3,j4)*zab2(j4,j2,j3,j1)/zb(j1,j2) &
3066  -s3(j2,j3,j4)*za(j2,j4))
3067 
3068  return
3069 
3070  end function a0phiggggpmmm
3071 
3072 
3073  function a0phiggggmpmp(j1,j2,j3,j4,za,zb,sprod)
3074  integer :: j1,j2,j3,j4
3075  complex(dp) :: za(4,4), zb(4,4)
3076  real(dp) :: sprod(4,4)
3077  complex(dp) :: a0phiggggmpmp
3078 
3079  a0phiggggmpmp = za(j1,j3)**4 &
3080  /(za(j1,j2)*za(j2,j3)*za(j3,j4)*za(j4,j1))
3081 
3082  return
3083 
3084  end function a0phiggggmpmp
3085 
3086 
3087  function a0phiggggmmpp(j1,j2,j3,j4,za,zb,sprod)
3088  integer :: j1,j2,j3,j4
3089  complex(dp) :: za(4,4), zb(4,4)
3090  real(dp) :: sprod(4,4)
3091  complex(dp) :: a0phiggggmmpp
3092 
3093  a0phiggggmmpp = za(j1,j2)**4 &
3094  /(za(j1,j2)*za(j2,j3)*za(j3,j4)*za(j4,j1))
3095 
3096  return
3097 
3098  end function a0phiggggmmpp
3099 
3100 
3101  function a0phiggggmmmm(j1,j2,j3,j4,za,zb,sprod)
3102  integer :: j1,j2,j3,j4
3103  complex(dp) :: za(4,4), zb(4,4)
3104  real(dp) :: sprod(4,4),qsq
3105  complex(dp) :: a0phiggggmmmm
3106 
3107  qsq = sprod(j1,j2)+sprod(j1,j3)+sprod(j1,j4)+sprod(j2,j3)+sprod(j2,j4)+sprod(j3,j4)
3108  a0phiggggmmmm=qsq**2/(zb(j1,j2)*zb(j2,j3)*zb(j3,j4)*zb(j4,j1))
3109 
3110  return
3111 
3112  end function a0phiggggmmmm
3113 
3114  !-------------------------------------------------------------------------
3115  !-- generic functions below
3116  !- MCFM spinors
3117  subroutine spinoru2(n,p,za,zb,s)
3118  implicit none
3119  integer, intent(in) :: n
3120  real(dp), intent(in) :: p(4,n)
3121  complex(dp), intent(out) :: za(n,n), zb(n,n)
3122  real(dp), intent(out) :: s(n,n)
3123  integer :: i,j
3124  complex(dp) :: c23(n), f(n)
3125  real(dp) :: rt(n)
3126 
3127  !---if one of the vectors happens to be zero this routine fails.
3128  do j=1,n
3129  za(j,j)=czero
3130  zb(j,j)=za(j,j)
3131 
3132  !-----positive energy case
3133  if (p(1,j) .gt. zero) then
3134  rt(j)=sqrt(abs(p(2,j)+p(1,j)))
3135  c23(j)=cmplx(p(4,j),-p(3,j),kind=dp)
3136  f(j)=(one,zero)
3137  else
3138  !-----negative energy case
3139  rt(j)=sqrt(abs(-p(1,j)-p(2,j)))
3140  c23(j)=cmplx(-p(4,j),p(3,j),kind=dp)
3141  f(j)=ci
3142  endif
3143  enddo
3144 
3145  do i=2,n
3146 
3147  do j=1,i-1
3148  s(i,j)=two*scr(p(:,i),p(:,j))
3149  za(i,j)=f(i)*f(j) * ( c23(i)*cmplx(rt(j)/(rt(i)+1d-16),kind=dp)-c23(j)*cmplx(rt(i)/(rt(j)+1d-16),kind=dp) )
3150 
3151  if (abs(s(i,j)).lt.1d-5) then
3152  zb(i,j)=-(f(i)*f(j))**2*conjg(za(i,j))
3153  else
3154  zb(i,j)=-cmplx(s(i,j),kind=dp)/(za(i,j)+1d-16)
3155  endif
3156 
3157  za(j,i)=-za(i,j)
3158  zb(j,i)=-zb(i,j)
3159  s(j,i)=s(i,j)
3160 
3161  enddo
3162 
3163  enddo
3164 
3165  return
3166 
3167  end subroutine spinoru2
3168 
3169 
3170 
3171 end module modhiggsjj
modhiggsjj::a0_zz_4f
complex(dp) function, dimension(-1:1,-1:1) a0_zz_4f(j1, j2, j3, j4, za, zb, sprod, iSel, jSel)
Definition: mod_HiggsJJ.F90:2358
modparameters::includegammastar
logical, public includegammastar
Definition: mod_Parameters.F90:213
modparameters::pdfstr_
integer, parameter, public pdfstr_
Definition: mod_Parameters.F90:1127
modhiggsjj::a0phiqbqgg_xxx
complex(dp) function, dimension(1:2,-1:1,-1:1,-1:1) a0phiqbqgg_xxx(j1, j2, j3, j4, za, zb, sprod)
Definition: mod_HiggsJJ.F90:2982
modparameters::vev
real(8), public vev
Definition: mod_Parameters.F90:249
modparameters::twosc
real(8), public twosc
Definition: mod_Parameters.F90:255
modparameters::couplazff
real(dp), public couplazff
Definition: mod_Parameters.F90:1080
modparameters::pdfastr_
integer, parameter, public pdfastr_
Definition: mod_Parameters.F90:1133
modparameters::ewp_chm_right
complex(8), public ewp_chm_right
Definition: mod_Parameters.F90:708
modhiggsjj::evalamp_wbfh_unsymm_sa_select
subroutine, public evalamp_wbfh_unsymm_sa_select(p, iSel, jSel, zz_fusion, iflip, res)
Definition: mod_HiggsJJ.F90:898
modhiggsjj::a0phiqbqqbq_xx
complex(dp) function, dimension(1:2,-1:1,-1:1) a0phiqbqqbq_xx(j1, j2, j3, j4, za, zb, sprod)
Definition: mod_HiggsJJ.F90:3023
modhiggsjj::tagbot
real(dp), parameter, public tagbot
Definition: mod_HiggsJJ.F90:14
modparameters::hvvspinzerodynamiccoupling
complex(8) function hvvspinzerodynamiccoupling(index, sWplus, sWminus, sWW, tryWWcoupl)
Definition: mod_Parameters.F90:1161
modparameters::ghg2
complex(8), public ghg2
Definition: mod_Parameters.F90:371
modparameters::pdfglu_
integer, parameter, public pdfglu_
Definition: mod_Parameters.F90:1124
modparameters::pdftop_
integer, parameter, public pdftop_
Definition: mod_Parameters.F90:1130
modparameters::ezp_up_left
complex(8), public ezp_up_left
Definition: mod_Parameters.F90:606
modparameters::aveqq
real(dp), parameter, public aveqq
Definition: mod_Parameters.F90:1021
modparameters::ga_wprime
real(8), public ga_wprime
Definition: mod_Parameters.F90:713
modparameters::a2
complex(8), public a2
Definition: mod_Parameters.F90:924
modhiggsjj::a0_vv_4f
complex(dp) function, dimension(-1:1,-1:1) a0_vv_4f(j1, j2, j3, j4, za, zb, sprod, mv, ga_v, useWWcoupl, Wpm_flip)
Definition: mod_HiggsJJ.F90:2267
modparameters::qul
real(8), parameter, public qul
Definition: mod_Parameters.F90:1036
modparameters::ga_w
real(8), public ga_w
Definition: mod_Parameters.F90:229
modmisc::swapr
subroutine swapr(i, j)
Definition: mod_Misc.F90:403
modparameters::pdfdn_
integer, parameter, public pdfdn_
Definition: mod_Parameters.F90:1125
modhiggsjj::evalamp_sbfh_unsymm_sa_select_exact
subroutine, public evalamp_sbfh_unsymm_sa_select_exact(p, iSel, jSel, rSel, sSel, res)
Definition: mod_HiggsJJ.F90:286
modparameters::cf
real(dp), parameter, public cf
Definition: mod_Parameters.F90:1017
modparameters::one
real(8), parameter, public one
Definition: mod_Parameters.F90:83
modparameters::m_zprime
real(8), public m_zprime
Definition: mod_Parameters.F90:619
modparameters::ezp_bot_right
complex(8), public ezp_bot_right
Definition: mod_Parameters.F90:615
modparameters::ewp_top_right
complex(8), public ewp_top_right
Definition: mod_Parameters.F90:710
modparameters::two
real(8), parameter, public two
Definition: mod_Parameters.F90:84
modparameters::couplaffsq
real(dp), public couplaffsq
Definition: mod_Parameters.F90:1081
modparameters::includevprime
logical, public includevprime
Definition: mod_Parameters.F90:214
modparameters::ci
complex(8), parameter, public ci
Definition: mod_Parameters.F90:88
modparameters::avegg
real(dp), parameter, public avegg
Definition: mod_Parameters.F90:1019
modparameters::lambda
real(8), parameter, public lambda
Definition: mod_Parameters.F90:354
modhiggsjj::tag1
real(dp), parameter, public tag1
Definition: mod_HiggsJJ.F90:12
modparameters::zero
real(8), parameter, public zero
Definition: mod_Parameters.F90:85
modparameters::gwsq
real(8), public gwsq
Definition: mod_Parameters.F90:250
modparameters::a3
complex(8), public a3
Definition: mod_Parameters.F90:925
modparameters::ezp_chm_right
complex(8), public ezp_chm_right
Definition: mod_Parameters.F90:609
modparameters::nf
real(dp), parameter, public nf
Definition: mod_Parameters.F90:1014
modparameters::m_wprime
real(8), public m_wprime
Definition: mod_Parameters.F90:712
modparameters::ar_qdn
real(8), public ar_qdn
Definition: mod_Parameters.F90:1064
modparameters::ckmbare
real(8) function ckmbare(id1in, id2in)
Definition: mod_Parameters.F90:1563
modhiggsjj::evalamp_sbfh_unsymm_sa
subroutine, public evalamp_sbfh_unsymm_sa(p, res)
Definition: mod_HiggsJJ.F90:24
modparameters::bl
real(8), public bl
Definition: mod_Parameters.F90:1066
modparameters::ezp_chm_left
complex(8), public ezp_chm_left
Definition: mod_Parameters.F90:608
modparameters::ca
real(dp), parameter, public ca
Definition: mod_Parameters.F90:1016
modhiggsjj::tag2
real(dp), parameter, public tag2
Definition: mod_HiggsJJ.F90:13
modparameters::xn
real(dp), parameter, public xn
Definition: mod_Parameters.F90:1015
modhiggsjj::a0phigggg_xxxx
complex(dp) function, dimension(1:2,-1:1,-1:1,-1:1,-1:1) a0phigggg_xxxx(j1, j2, j3, j4, za, zb, sprod)
Definition: mod_HiggsJJ.F90:2946
modparameters::couplwffsq
real(dp), public couplwffsq
Definition: mod_Parameters.F90:1078
modparameters::xw
real(8), public xw
Definition: mod_Parameters.F90:253
modparameters::al_qdn
real(8), public al_qdn
Definition: mod_Parameters.F90:1065
modparameters::ga_zprime
real(8), public ga_zprime
Definition: mod_Parameters.F90:620
modhiggsjj::spinoru2
subroutine spinoru2(n, p, za, zb, s)
Definition: mod_HiggsJJ.F90:3118
modparameters::ewp_up_right
complex(8), public ewp_up_right
Definition: mod_Parameters.F90:706
modhiggsjj::flip
integer function, dimension(2) flip(i, a1, a2)
Definition: mod_HiggsJJ.F90:2252
modparameters::ezp_str_right
complex(8), public ezp_str_right
Definition: mod_Parameters.F90:613
modparameters::czero
complex(8), parameter, public czero
Definition: mod_Parameters.F90:86
modparameters::pdfup_
integer, parameter, public pdfup_
Definition: mod_Parameters.F90:1126
modparameters::couplzffsq
real(dp), public couplzffsq
Definition: mod_Parameters.F90:1079
modparameters::ezp_up_right
complex(8), public ezp_up_right
Definition: mod_Parameters.F90:607
modparameters::aveqg
real(dp), parameter, public aveqg
Definition: mod_Parameters.F90:1020
modhiggsjj::me2_qbqggh
subroutine me2_qbqggh(j1, j2, j3, j4, za, zb, sprod, res)
Definition: mod_HiggsJJ.F90:2859
modparameters::pdfaup_
integer, parameter, public pdfaup_
Definition: mod_Parameters.F90:1132
modparameters::qdl
real(8), parameter, public qdl
Definition: mod_Parameters.F90:1038
modparameters::ga_z
real(8), public ga_z
Definition: mod_Parameters.F90:227
modparameters::m_w
real(8), public m_w
Definition: mod_Parameters.F90:228
modparameters::spinavg
real(8), parameter, public spinavg
Definition: mod_Parameters.F90:91
modparameters::ewp_chm_left
complex(8), public ewp_chm_left
Definition: mod_Parameters.F90:707
modparameters::al_qup
real(8), public al_qup
Definition: mod_Parameters.F90:1063
modparameters::ga_v
real(8), public ga_v
Definition: mod_Parameters.F90:78
modparameters::ewp_up_left
complex(8), public ewp_up_left
Definition: mod_Parameters.F90:705
modparameters::pdfchm_
integer, parameter, public pdfchm_
Definition: mod_Parameters.F90:1128
modhiggsjj::me2_ggggh
subroutine me2_ggggh(j1, j2, j3, j4, za, zb, sprod, res)
Definition: mod_HiggsJJ.F90:2814
modparameters::ezp_dn_left
complex(8), public ezp_dn_left
Definition: mod_Parameters.F90:610
modparameters
Definition: mod_Parameters.F90:1
modhiggsjj::evalamp_wbfh_unsymm_sa
subroutine, public evalamp_wbfh_unsymm_sa(p, res)
Definition: mod_HiggsJJ.F90:526
modhiggsjj::a0_ww_4f
complex(dp) function, dimension(-1:1,-1:1) a0_ww_4f(j1, j2, j3, j4, za, zb, sprod, iSel, jSel, useWWcoupl, Wpm_flip)
Definition: mod_HiggsJJ.F90:2640
modmisc
Definition: mod_Misc.F90:1
modparameters::pdfbot_
integer, parameter, public pdfbot_
Definition: mod_Parameters.F90:1129
modhiggsjj::evalamp_wbfh_unsymm_sa_select_exact
subroutine, public evalamp_wbfh_unsymm_sa_select_exact(p, iSel, jSel, rSel, sSel, res)
Definition: mod_HiggsJJ.F90:1872
modparameters::ezp_dn_right
complex(8), public ezp_dn_right
Definition: mod_Parameters.F90:611
modparameters::a1
complex(8), public a1
Definition: mod_Parameters.F90:923
modparameters::symmfac
real(8), parameter, public symmfac
Definition: mod_Parameters.F90:91
modparameters::qur
real(8), parameter, public qur
Definition: mod_Parameters.F90:1037
modparameters::pdfachm_
integer, parameter, public pdfachm_
Definition: mod_Parameters.F90:1134
modparameters::m_z
real(8), public m_z
Definition: mod_Parameters.F90:226
modparameters::dp
integer, parameter, public dp
Definition: mod_Parameters.F90:11
modmisc::scr
double precision function scr(p1, p2)
Definition: mod_Misc.F90:135
modparameters::alphas
real(dp), public alphas
Definition: mod_Parameters.F90:269
modhiggsjj::a0phiggggmpmp
complex(dp) function a0phiggggmpmp(j1, j2, j3, j4, za, zb, sprod)
Definition: mod_HiggsJJ.F90:3074
modparameters::ewp_top_left
complex(8), public ewp_top_left
Definition: mod_Parameters.F90:709
modparameters::ghg4
complex(8), public ghg4
Definition: mod_Parameters.F90:373
modhiggsjj::evalamp_sbfh_unsymm_sa_select
subroutine, public evalamp_sbfh_unsymm_sa_select(p, iSel, jSel, flav_tag, iflip, res)
Definition: mod_HiggsJJ.F90:135
modhiggsjj::a0phiggggmmmm
complex(dp) function a0phiggggmmmm(j1, j2, j3, j4, za, zb, sprod)
Definition: mod_HiggsJJ.F90:3102
modhiggsjj::a0phiggggmmpp
complex(dp) function a0phiggggmmpp(j1, j2, j3, j4, za, zb, sprod)
Definition: mod_HiggsJJ.F90:3088
modparameters::ezp_str_left
complex(8), public ezp_str_left
Definition: mod_Parameters.F90:612
modhiggsjj::a0phiggggpmmm
complex(dp) function a0phiggggpmmm(j1, j2, j3, j4, za, zb, sprod)
Definition: mod_HiggsJJ.F90:3048
modhiggsjj::me2_qbqqbq
subroutine me2_qbqqbq(j1, j2, j3, j4, za, zb, sprod, res_diff, res_id)
Definition: mod_HiggsJJ.F90:2903
modmisc::swapi
subroutine swapi(i, j)
Definition: mod_Misc.F90:393
modparameters::qdr
real(8), parameter, public qdr
Definition: mod_Parameters.F90:1039
modparameters::pdfadn_
integer, parameter, public pdfadn_
Definition: mod_Parameters.F90:1131
modparameters::pdfabot_
integer, parameter, public pdfabot_
Definition: mod_Parameters.F90:1135
modparameters::ezp_bot_left
complex(8), public ezp_bot_left
Definition: mod_Parameters.F90:614
modmisc::swap
Definition: mod_Misc.F90:5
modhiggsjj
Definition: mod_HiggsJJ.F90:1
modparameters::ar_qup
real(8), public ar_qup
Definition: mod_Parameters.F90:1062