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_TopDecay.F90
Go to the documentation of this file.
1 MODULE modtopdecay
2 implicit none
3 
4 public :: topdecay, wdecay
5 public :: ubarspi_dirac, vspi_dirac
6 private
7 
8  CONTAINS
9 
10 
11 
12 
13 
14 SUBROUTINE topdecay(Flavor,Mom,Spinor,TopHel)
15 use modmisc
16 use modparameters
17 implicit none
18 real(8) :: mom(1:4,1:3)
19 integer :: flavor
20 integer,optional :: tophel
21 complex(8) :: spinor(1:4)
22 real(8) :: zeros(1:6),topmom(1:4),nwafactor_top
23 complex(8) :: spi(1:4),barspi(1:4),botspi(1:4),wcurr(1:4)
24 real(8) :: nwafactor_w
25 complex(8) :: wprop
26 real(8),parameter :: nc=xn,nflav=2
27 
28 nwafactor_w = 1d0/dsqrt(2d0*ga_w*m_w)
29 wprop = (0d0,-1d0)*nwafactor_w
30 
31 
32 
33 ! zeros(1:4) = dble(TopMom(1:4)) - Mom(1:4,1)-Mom(1:4,2)-Mom(1:4,3)
34 ! if( any(abs(zeros(1:4)/dble(TopMom(1))).gt.1d-8) ) then
35 ! print *, "ERROR: energy-momentum violation in SUBROUTINE TopDecay(): ",zeros(1:4)
36 ! print *, "momentum dump:"
37 ! print *, Mom(1:4,1:12)
38 ! endif
39 !
40 ! zeros(1) = dble(TopMom(1:4).dot.TopMom(1:4)) - m_Top**2
41 ! zeros(2) = Mom(1:4,1).dot.Mom(1:4,1)
42 ! zeros(3) = Mom(1:4,2).dot.Mom(1:4,2)
43 ! zeros(4) = Mom(1:4,3).dot.Mom(1:4,3)
44 ! if( any(abs(zeros(1:6)/dble(TopMom(1))**2).gt.1d-8) ) then
45 ! print *, "ERROR: onshell-ness violation in SUBROUTINE TopDecay(): ",zeros(1:6)
46 ! print *, "momentum dump:"
47 ! print *, Mom(1:4,1:12)
48 ! endif
49 
50 
51 
52  topmom(1:4) = mom(1:4,1)+mom(1:4,2)+mom(1:4,3)
53  if( topdecays.eq.0 ) then
54 ! if( Flavor.eq.Top_ ) call ubarSpi_Dirac(dcmplx(TopMom(1:4)),M_Top,TopHel,Spinor(1:4))
55 ! if( Flavor.eq.ATop_ ) call vSpi_Dirac(dcmplx(TopMom(1:4)),M_Top,TopHel,Spinor(1:4))
56  spinor(1:4) = 0d0
57  RETURN
58  endif
59 ! elseif( TOPDECAYS.eq.1 ) then
60  nwafactor_top = 1d0/dsqrt(2d0*ga_top*m_top)
61 ! elseif( TOPDECAYS.eq.2 ) then
62 ! NWAFactor_Top = 1d0/dsqrt(2d0*Ga_Top*m_Top)
63 ! NWAFactor_Top = NWAFactor_Top * dsqrt(dsqrt(Nc*NFlav)**2)
64 ! elseif( TOPDECAYS.eq.3 .or. TOPDECAYS.eq.4 ) then
65 ! NWAFactor_Top = 1d0/dsqrt(2d0*Ga_Top*m_Top)
66 ! NWAFactor_Top = NWAFactor_Top * dsqrt(dsqrt(Nc*NFlav))
67 ! endif
68 
69 
70  if( flavor.eq.top_ ) then ! Top quark decay
71 ! assemble lepton current
72  call ubarspi_weyl(dcmplx(mom(1:4,1)),-1,botspi(1:4)) ! bot
73  call vspi_weyl(dcmplx(mom(1:4,2)),+1,spi(1:4)) ! l+ or dn_bar
74  call ubarspi_weyl(dcmplx(mom(1:4,3)),-1,barspi(1:4)) ! nu or up
75  wcurr(1:4) = vbqq_weyl(barspi(1:4),spi(1:4)) * wprop * gwsq ! vbqq introduces -i/Sqrt(2)
76 
77 ! connect to quark current
78  barspi(1:4) = botspi(1:4)
79  spinor(1:4) = vgq_weyl( wcurr(1:4),barspi(1:4) ) ! vgq introduces -i/Sqrt(2)
80  spinor(1:4) =( spb2_weyl(spinor(1:4),dcmplx(topmom(1:4))) + m_top*spinor(1:4) ) * nwafactor_top
81  spinor(1:4) = weyltodirac(spinor(1:4))
82  elseif( flavor.eq.atop_ ) then ! Anti-Top quark decay
83 ! assemble lepton current
84  call vspi_weyl(dcmplx(mom(1:4,1)),+1,botspi(1:4)) ! Abot
85  call ubarspi_weyl(dcmplx(mom(1:4,2)),-1,barspi(1:4)) ! l- or dn
86  call vspi_weyl(dcmplx(mom(1:4,3)),+1,spi(1:4)) ! nubar or up_bar
87  wcurr(1:4) = vbqq_weyl(barspi(1:4),spi(1:4)) * wprop * gwsq ! vbqq introduces -i/Sqrt(2)
88 
89 ! connect to quark current:
90  spi(1:4) = botspi(1:4)
91  spinor(1:4) = vbqg_weyl( spi(1:4),wcurr(1:4) )! vbqg introduces -i/Sqrt(2)
92  spinor(1:4) = ( spi2_weyl(dcmplx(topmom(1:4)),spinor(1:4)) - m_top*spinor(1:4) ) * nwafactor_top
93  spinor(1:4) = weyltodirac(spinor(1:4))
94  endif
95 
96 
97 RETURN
98 END SUBROUTINE
99 
100 
101 
102 
103 
104 
105 SUBROUTINE wdecay(Charge,Mom,WCurr)
106 use modmisc
107 use modparameters
108 implicit none
109 real(8) :: mom(1:4,1:2)
110 integer :: charge
111 real(8) :: wmom(1:4)
112 complex(8) :: spi(1:4),barspi(1:4),botspi(1:4),wcurr(1:4)
113 real(8) :: nwafactor_w
114 complex(8) :: wprop
115 real(8),parameter :: nc=xn,nflav=2
116 
117 nwafactor_w = 1d0/dsqrt(2d0*ga_w*m_w)
118 wprop = (0d0,-1d0)*dsqrt(gwsq)*nwafactor_w
119 
120 
121 
122  wmom(1:4) = mom(1:4,1)+mom(1:4,2)
123  if( topdecays.eq.0 ) then
124  wcurr(1:4) = 0d0
125  RETURN
126  endif
127 
128 
129  if( charge.eq.wp_ ) then ! W+ decay
130 ! assemble lepton current
131  call vspi_weyl(dcmplx(mom(1:4,1)),+1,spi(1:4)) ! l+ or dn_bar
132  call ubarspi_weyl(dcmplx(mom(1:4,2)),-1,barspi(1:4)) ! nu or up
133  wcurr(1:4) = vbqq_weyl(barspi(1:4),spi(1:4)) * wprop ! vbqq introduces -i/Sqrt(2)
134 
135  elseif( charge.eq.wm_ ) then ! W- decay
136 ! assemble lepton current
137  call ubarspi_weyl(dcmplx(mom(1:4,1)),-1,barspi(1:4)) ! l- or dn
138  call vspi_weyl(dcmplx(mom(1:4,2)),+1,spi(1:4)) ! nubar or up_bar
139  wcurr(1:4) = vbqq_weyl(barspi(1:4),spi(1:4)) * wprop ! vbqq introduces -i/Sqrt(2)
140 
141  endif
142 
143 
144 RETURN
145 END SUBROUTINE
146 
147 
148 
149 
150  SUBROUTINE ubarspi_weyl(p,i,ubarSpi) ! i=+1 is ES to Chir_Weyl(.false.), i=-1 is ES to Chir_Weyl(.true.)
151  use modmisc
152  implicit none
153  integer, intent(in):: i
154  complex(8), intent(in) :: p(4)
155  complex(8) :: ubarSpi(4)
156  complex(8) :: ephi
157  real(8) :: p0,px,py,pz
158  complex(8) :: fc, fc2
159 
160  p0=dreal(p(1))
161  px=dreal(p(2))
162  py=dreal(p(3))
163  pz=dreal(p(4))
164 
165  fc2 = p0 + pz
166  fc=cdsqrt(fc2)
167 
168  if (cdabs(fc2).gt.1d-15) then
169 
170  if (i.eq.1) then
171  ubarspi(1)=(0d0,0d0)
172  ubarspi(2)=(0d0,0d0)
173  ubarspi(3)=fc
174  ubarspi(4)=(px-(0d0,1d0)*py)/fc
175  elseif (i.eq.-1) then
176  ubarspi(1)=(px+(0d0,1d0)*py)/fc
177  ubarspi(2)=-fc
178  ubarspi(3)=(0d0,0d0)
179  ubarspi(4)=(0d0,0d0)
180  else
181  call error("wrong helicity setting in ubarSpi_Weyl")
182  endif
183 
184  else
185 
186  if (i.eq.1) then
187  ubarspi(1) = (0d0,0d0)
188  ubarspi(2) = (0d0,0d0)
189  ubarspi(3) = (0d0,0d0)
190  ubarspi(4) = dsqrt(2d0*p0)
191  elseif (i.eq.-1) then
192  ubarspi(1) = dsqrt(2d0*p0)
193  ubarspi(2) = (0d0,0d0)
194  ubarspi(3) = (0d0,0d0)
195  ubarspi(4) = (0d0,0d0)
196  else
197  call error("wrong helicity setting in ubarSpi_Weyl")
198  endif
199 
200  endif
201  return
202  END SUBROUTINE
203 
204 
205 
206 
207 
208  SUBROUTINE vspi_weyl(p,i,vSpi) ! i=+1 is ES to Chir_Weyl(.false.), i=-1 is ES to Chir_Weyl(.true.)
209  use modmisc
210  implicit none
211  integer, intent(in):: i
212  complex(8), intent(in) :: p(4)
213  complex(8) :: vSpi(4)
214  complex(8) :: ephi
215  real(8) :: p0,px,py,pz
216  real(8) :: nx,ny,nz,theta,phi
217  real(8) :: ct,ct2,st,st2,cphi,sphi
218  complex(8) :: fc2, fc
219 
220  p0=dreal(p(1))
221  px=dreal(p(2))
222  py=dreal(p(3))
223  pz=dreal(p(4))
224 
225  fc2 = p0 + pz
226  fc=cdsqrt(fc2)
227 
228  if (cdabs(fc2).gt.1d-15) then
229 
230  if (i.eq.1) then
231  vspi(1)=(0d0,0d0)
232  vspi(2)=(0d0,0d0)
233  vspi(3)=(px-(0d0,1d0)*py)/fc
234  vspi(4)=-fc
235  elseif (i.eq.-1) then
236  vspi(1)=fc
237  vspi(2)=(px+(0d0,1d0)*py)/fc
238  vspi(3)=(0d0,0d0)
239  vspi(4)=(0d0,0d0)
240  else
241  call error("wrong helicity setting in vSpi_Weyl")
242  endif
243 
244  else
245 
246  if (i.eq.1) then
247  vspi(1)=(0d0,0d0)
248  vspi(2)=(0d0,0d0)
249  vspi(3)=dsqrt(2d0*p0)
250  vspi(4)=(0d0,0d0)
251  elseif (i.eq.-1) then
252  vspi(1)=(0d0,0d0)
253  vspi(2)=dsqrt(2d0*p0)
254  vspi(3)=(0d0,0d0)
255  vspi(4)=(0d0,0d0)
256  else
257  call error("wrong helicity setting in vSpi_Weyl")
258  endif
259 
260  endif
261 
262  RETURN
263  END SUBROUTINE
264 
265 
266 
267 
268 
269 
270 
271 
272  FUNCTION vbqg_weyl(sp,e1)
273  use modmisc
274  implicit none
275  complex(8), intent(in) :: e1(:)
276  complex(8), intent(in) :: sp(:)
277  complex(8) :: vbqg_weyl(size(sp))
278  real(8), parameter :: sqrt2 = 1.4142135623730950488016887242096980786d0
279 
280  vbqg_weyl = (0d0,-1d0)/sqrt2*spi2_weyl(e1,sp)
281 
282  END FUNCTION
283 
284 
285 
286  function vgq_weyl(e1,sp)
287  use modmisc
288  implicit none
289  complex(8), intent(in) :: e1(:)
290  complex(8), intent(in) :: sp(:)
291  complex(8) :: vgq_weyl(size(sp))
292  real(8), parameter :: sqrt2 = 1.4142135623730950488016887242096980786d0
293 
294  vgq_weyl = (0d0,-1d0)/sqrt2*spb2_weyl(sp,e1)
295 
296  end function
297 
298 
299 
300  function weyltodirac(sp) ! unitary transformation U to convert a Weyl spinor into the Dirac representation
301  implicit none ! sp can be spinor or bar-spinor, i.e. U^dagger.sp = barsp.U
302  double complex :: sp(1:4)
303  double complex :: weyltodirac(1:4)
304  double precision,parameter :: sqrtfac=1d0/dsqrt(2d0)
305 
306  weyltodirac(1) = sqrtfac*(sp(1)+sp(3))
307  weyltodirac(2) = sqrtfac*(sp(2)+sp(4))
308  weyltodirac(3) = sqrtfac*(sp(1)-sp(3))
309  weyltodirac(4) = sqrtfac*(sp(2)-sp(4))
310  return
311  end function
312 
313 
314 
315 
316  function psp1_(sp1,sp2) result(res)
317  implicit none
318  complex(8), intent(in) :: sp1(:)
319  complex(8), intent(in) :: sp2(:)
320  complex(8) :: res
321 
322  res = sum(sp1(1:)*sp2(1:))
323 
324  end function
325 
326 
327 
328 
329 
330  subroutine ubarspi_dirac(p,m,i,f)
331  use modmisc
332  implicit none
333  integer i
334  real(8) m
335  complex(8) p(4)
336  complex(8) f(4),fc
337  real(8) p0,px,py,pz,fc2
338 
339  p0=dreal(p(1))
340  px=dreal(p(2))
341  py=dreal(p(3))
342  pz=dreal(p(4))
343 
344  fc2=p0+m
345  fc=cdsqrt( dcmplx(fc2))
346 ! fc=dsqrt(fc2)
347 
348  if (i.eq.1) then
349  f(1)=fc
350  f(2)=dcmplx(0d0,0d0)
351  f(3)=-1d0*pz*fc/fc2
352  f(4)=-(px-(0d0,1d0)*py)*fc/fc2
353  elseif (i.eq.-1) then
354  f(1)=dcmplx(0d0,0d0)
355  f(2)=fc
356  f(3)=-(px+(0d0,1d0)*py)*fc/fc2
357  f(4)=pz*fc/fc2
358  else
359  print *, "wrong helicity setting in ubarSpi"
360  stop
361  endif
362 
363  return
364  end subroutine
365 
366 
367 
368 
369 
370 
371 
372 
373  subroutine vspi_dirac(p,m,i,f)
374  use modmisc
375  implicit none
376  integer i
377  real(8) m
378  complex(8) p(4)
379  complex(8) f(4),fc
380  real(8) p0,px,py,pz,fc2
381 
382  p0=dreal(p(1))
383  px=dreal(p(2))
384  py=dreal(p(3))
385  pz=dreal(p(4))
386 
387  fc2 = p0+m
388  fc=cdsqrt(dcmplx(fc2))
389 ! fc=dsqrt(fc2)
390 
391  if (i.eq.1) then
392  f(1)=pz*fc/fc2
393  f(2)=(px+(0d0,1d0)*py)*fc/fc2
394  f(3)=fc
395  f(4)=dcmplx(0d0,0d0)
396  elseif (i.eq.-1) then
397  f(1)=(px-(0d0,1d0)*py)*fc/fc2
398  f(2)=-pz*fc/fc2
399  f(3)=dcmplx(0d0,0d0)
400  f(4)=fc
401  else
402  print *, "wrong helicity setting in vSpi"
403  endif
404 
405  return
406  end SUBROUTINE
407 
408 
409 
410 
411 
412 
413 END MODULE
modtopdecay::topdecay
subroutine, public topdecay(Flavor, Mom, Spinor, TopHel)
Definition: mod_TopDecay.F90:15
modtopdecay::vbqg_weyl
complex(8) function, dimension(size(sp)) vbqg_weyl(sp, e1)
Definition: mod_TopDecay.F90:273
modtopdecay::wdecay
subroutine, public wdecay(Charge, Mom, WCurr)
Definition: mod_TopDecay.F90:106
modparameters::ga_w
real(8), public ga_w
Definition: mod_Parameters.F90:229
modtopdecay
Definition: mod_TopDecay.F90:1
modparameters::topdecays
integer, public topdecays
Definition: mod_Parameters.F90:17
modtopdecay::weyltodirac
double complex function, dimension(1:4) weyltodirac(sp)
Definition: mod_TopDecay.F90:301
modmisc::error
subroutine error(Message, ErrNum)
Definition: mod_Misc.F90:366
modparameters::wm_
integer, target, public wm_
Definition: mod_Parameters.F90:1115
modparameters::gwsq
real(8), public gwsq
Definition: mod_Parameters.F90:250
modtopdecay::vspi_dirac
subroutine, public vspi_dirac(p, m, i, f)
Definition: mod_TopDecay.F90:374
modparameters::m_top
real(8), public m_top
Definition: mod_Parameters.F90:224
modtopdecay::vspi_weyl
subroutine vspi_weyl(p, i, vSpi)
Definition: mod_TopDecay.F90:209
modmisc::spi2_weyl
complex(8) function, dimension(4) spi2_weyl(v, sp)
Definition: mod_Misc.F90:711
modtopdecay::vgq_weyl
complex(8) function, dimension(size(sp)) vgq_weyl(e1, sp)
Definition: mod_TopDecay.F90:287
modparameters::xn
real(dp), parameter, public xn
Definition: mod_Parameters.F90:1015
modmisc::vbqq_weyl
complex(8) function, dimension(dv) vbqq_weyl(sp1, sp2)
Definition: mod_Misc.F90:626
modparameters::m_w
real(8), public m_w
Definition: mod_Parameters.F90:228
modparameters::atop_
integer, target, public atop_
Definition: mod_Parameters.F90:1110
modtopdecay::ubarspi_weyl
subroutine ubarspi_weyl(p, i, ubarSpi)
Definition: mod_TopDecay.F90:151
modparameters
Definition: mod_Parameters.F90:1
modmisc
Definition: mod_Misc.F90:1
modtopdecay::psp1_
complex(8) function psp1_(sp1, sp2)
Definition: mod_TopDecay.F90:317
modparameters::wp_
integer, target, public wp_
Definition: mod_Parameters.F90:1096
modparameters::ga_top
real(8), public ga_top
Definition: mod_Parameters.F90:225
modtopdecay::ubarspi_dirac
subroutine, public ubarspi_dirac(p, m, i, f)
Definition: mod_TopDecay.F90:331
modparameters::top_
integer, target, public top_
Definition: mod_Parameters.F90:1088
modmisc::spb2_weyl
complex(8) function, dimension(4) spb2_weyl(sp, v)
Definition: mod_Misc.F90:657
Generated by   doxygen 1.8.17