JHUGen MELA  v2.4.1
Matrix element calculations as used in JHUGen. MELA is an important tool that was used for the Higgs boson discovery and for precise measurements of its structure and interactions. Please see the website https://spin.pha.jhu.edu/ and papers cited there for more details, and kindly cite those papers when using this code.
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
Generated by   doxygen 1.8.17