JHUGen MELA  JHUGen v7.5.6, MELA v2.4.2
Matrix element calculations as used in JHUGen.
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