go-hep.org/x/hep@v0.38.1/fmom/pxpypze.go (about)

     1  // Copyright ©2017 The go-hep Authors. All rights reserved.
     2  // Use of this source code is governed by a BSD-style
     3  // license that can be found in the LICENSE file.
     4  
     5  package fmom
     6  
     7  import (
     8  	"fmt"
     9  	"math"
    10  )
    11  
    12  type PxPyPzE struct {
    13  	P4 Vec4
    14  }
    15  
    16  func NewPxPyPzE(px, py, pz, e float64) PxPyPzE {
    17  	return PxPyPzE{P4: Vec4{X: px, Y: py, Z: pz, T: e}}
    18  }
    19  
    20  func (p4 PxPyPzE) String() string {
    21  	return fmt.Sprintf(
    22  		"fmom.P4{Px:%v, Py:%v, Pz:%v, E:%v}",
    23  		p4.Px(), p4.Py(), p4.Pz(), p4.E(),
    24  	)
    25  }
    26  
    27  func (p4 *PxPyPzE) Clone() P4 {
    28  	pp := *p4
    29  	return &pp
    30  }
    31  
    32  func (p4 *PxPyPzE) Px() float64 { return p4.P4.X }
    33  func (p4 *PxPyPzE) Py() float64 { return p4.P4.Y }
    34  func (p4 *PxPyPzE) Pz() float64 { return p4.P4.Z }
    35  func (p4 *PxPyPzE) E() float64  { return p4.P4.T }
    36  
    37  func (p4 *PxPyPzE) X() float64 { return p4.P4.X }
    38  func (p4 *PxPyPzE) Y() float64 { return p4.P4.Y }
    39  func (p4 *PxPyPzE) Z() float64 { return p4.P4.Z }
    40  func (p4 *PxPyPzE) T() float64 { return p4.P4.T }
    41  
    42  func (p4 *PxPyPzE) M2() float64 {
    43  	px := p4.Px()
    44  	py := p4.Py()
    45  	pz := p4.Pz()
    46  	e := p4.E()
    47  
    48  	m2 := e*e - (px*px + py*py + pz*pz)
    49  	return m2
    50  }
    51  
    52  func (p4 *PxPyPzE) M() float64 {
    53  	m2 := p4.M2()
    54  	if m2 < 0.0 {
    55  		return -math.Sqrt(-m2)
    56  	}
    57  	return +math.Sqrt(+m2)
    58  }
    59  
    60  func (p4 *PxPyPzE) Eta() float64 {
    61  	px := p4.Px()
    62  	py := p4.Py()
    63  	pz := p4.Pz()
    64  	e := p4.E()
    65  
    66  	// FIXME: should we use a more underflow-friendly formula:
    67  	//  sqrt(a**2 + b**2)
    68  	//   => y.sqrt(1+(x/y)**2) where y=max(|a|,|b|) and x=min(|a|,|b|)
    69  	//
    70  	p := math.Sqrt(px*px + py*py + pz*pz)
    71  	switch p {
    72  	case 0.0:
    73  		return 0
    74  	case +pz:
    75  		return math.Inf(+1)
    76  	case -pz:
    77  		return math.Inf(-1)
    78  	}
    79  	// flip if negative e
    80  	sign := 1.0
    81  	if e < 0 {
    82  		sign = -1.0
    83  	}
    84  	return sign * 0.5 * math.Log((p+pz)/(p-pz))
    85  }
    86  
    87  func (p4 *PxPyPzE) Phi() float64 {
    88  	e := p4.E()
    89  	// flip if negative e
    90  	sign := 1.0
    91  	if e < 0 {
    92  		sign = -1.0
    93  	}
    94  	px := sign * p4.Px()
    95  	py := sign * p4.Py()
    96  	if px == 0.0 && py == 0.0 {
    97  		return 0
    98  	}
    99  	return math.Atan2(py, px)
   100  }
   101  
   102  func (p4 *PxPyPzE) P2() float64 {
   103  	px := p4.Px()
   104  	py := p4.Py()
   105  	pz := p4.Pz()
   106  
   107  	return px*px + py*py + pz*pz
   108  }
   109  
   110  func (p4 *PxPyPzE) P() float64 {
   111  	e := p4.E()
   112  	// flip if negative e
   113  	sign := 1.0
   114  	if e < 0 {
   115  		sign = -1.0
   116  	}
   117  	p2 := p4.P2()
   118  	return sign * math.Sqrt(p2)
   119  }
   120  
   121  func (p4 *PxPyPzE) CosPhi() float64 {
   122  	px := p4.Px()
   123  	ipt := p4.IPt()
   124  	return px * ipt
   125  }
   126  
   127  func (p4 *PxPyPzE) SinPhi() float64 {
   128  	py := p4.Py()
   129  	ipt := p4.IPt()
   130  	return py * ipt
   131  }
   132  
   133  func (p4 *PxPyPzE) TanTh() float64 {
   134  	pt := p4.Pt()
   135  	pz := p4.Pz()
   136  	return pt / pz
   137  }
   138  
   139  func (p4 *PxPyPzE) CotTh() float64 {
   140  	pt := p4.Pt()
   141  	pz := p4.Pz()
   142  	return pz / pt
   143  }
   144  
   145  func (p4 *PxPyPzE) CosTh() float64 {
   146  	pz := p4.Pz()
   147  	p := p4.P()
   148  	return pz / p
   149  }
   150  
   151  func (p4 *PxPyPzE) SinTh() float64 {
   152  	pt := p4.Pt()
   153  	p := p4.P()
   154  	return pt / p
   155  }
   156  
   157  func (p4 *PxPyPzE) Pt() float64 {
   158  	e := p4.E()
   159  	px := p4.Px()
   160  	py := p4.Py()
   161  	// flip if negative e
   162  	sign := 1.0
   163  	if e < 0 {
   164  		sign = -1.0
   165  	}
   166  	return sign * math.Sqrt(px*px+py*py)
   167  }
   168  
   169  func (p4 *PxPyPzE) Et() float64 {
   170  	// to be improved
   171  	e := p4.E()
   172  	sinth := p4.SinTh()
   173  	return e * sinth
   174  }
   175  
   176  func (p4 *PxPyPzE) IPt() float64 {
   177  	pt := p4.Pt()
   178  	return 1.0 / pt
   179  }
   180  
   181  func (p4 *PxPyPzE) Rapidity() float64 {
   182  	e := p4.E()
   183  	pz := p4.Pz()
   184  	switch e {
   185  	case 0.0:
   186  		return 0.0
   187  	case +pz:
   188  		return math.Inf(+1)
   189  	case -pz:
   190  		return math.Inf(-1)
   191  	}
   192  	// invariant under flipping of 4-mom with negative energy
   193  	return 0.5 * math.Log((e+pz)/(e-pz))
   194  }
   195  
   196  func (p4 *PxPyPzE) Set(p P4) {
   197  	p4.P4.X = p.Px()
   198  	p4.P4.Y = p.Py()
   199  	p4.P4.Z = p.Pz()
   200  	p4.P4.T = p.E()
   201  }
   202  
   203  func (p4 *PxPyPzE) SetPtEtaPhiM(pt, eta, phi, m float64) {
   204  	sin, cos := math.Sincos(phi)
   205  	pt = math.Abs(pt)
   206  	p4.P4.X = pt * cos
   207  	p4.P4.Y = pt * sin
   208  	p4.P4.Z = pt * math.Sinh(eta)
   209  
   210  	p2 := p4.P4.X*p4.P4.X + p4.P4.Y*p4.P4.Y + p4.P4.Z*p4.P4.Z
   211  	m2 := m * m
   212  	switch {
   213  	case m >= 0:
   214  		p4.P4.T = math.Sqrt(p2 + m2)
   215  	default:
   216  		p4.P4.T = math.Sqrt(math.Max(0, p2-m2))
   217  	}
   218  }
   219  
   220  func (p4 *PxPyPzE) SetPtEtaPhiE(pt, eta, phi, e float64) {
   221  	sin, cos := math.Sincos(phi)
   222  	pt = math.Abs(pt)
   223  	p4.P4.X = pt * cos
   224  	p4.P4.Y = pt * sin
   225  	p4.P4.Z = pt * math.Sinh(eta)
   226  	p4.P4.T = e
   227  }
   228  
   229  var (
   230  	_ P4           = (*PxPyPzE)(nil)
   231  	_ fmt.Stringer = (*PxPyPzE)(nil)
   232  )