diff --git a/cmd/lhef2hepmc/main.go b/cmd/lhef2hepmc/main.go
index d50a17554..52f0f5eed 100644
--- a/cmd/lhef2hepmc/main.go
+++ b/cmd/lhef2hepmc/main.go
@@ -132,21 +132,21 @@ func main() {
 			Barcode: -1,
 		}
 		p1 := hepmc.Particle{
-			Momentum: fmom.PxPyPzE{
+			Momentum: fmom.NewPxPyPzE(
 				0, 0,
 				dec.Run.EBMUP[0],
 				dec.Run.EBMUP[0],
-			},
+			),
 			PdgID:   dec.Run.IDBMUP[0],
 			Status:  4,
 			Barcode: 1,
 		}
 		p2 := hepmc.Particle{
-			Momentum: fmom.PxPyPzE{
+			Momentum: fmom.NewPxPyPzE(
 				0, 0,
 				dec.Run.EBMUP[1],
 				dec.Run.EBMUP[1],
-			},
+			),
 			PdgID:   dec.Run.IDBMUP[1],
 			Status:  4,
 			Barcode: 2,
@@ -170,12 +170,12 @@ func main() {
 			}
 			nmax += 1
 			vtx.AddParticleOut(&hepmc.Particle{
-				Momentum: fmom.PxPyPzE{
+				Momentum: fmom.NewPxPyPzE(
 					lhevt.PUP[i][0],
 					lhevt.PUP[i][1],
 					lhevt.PUP[i][2],
 					lhevt.PUP[i][3],
-				},
+				),
 				GeneratedMass: lhevt.PUP[i][4],
 				PdgID:         lhevt.IDUP[i],
 				Status:        1,
diff --git a/fads/fastjet_finder.go b/fads/fastjet_finder.go
index 6b51bdce7..b9a6d11ed 100644
--- a/fads/fastjet_finder.go
+++ b/fads/fastjet_finder.go
@@ -149,8 +149,10 @@ func (tsk *FastJetFinder) Process(ctx fwk.Context) error {
 	dphiMax := 0.0
 	output = make([]Candidate, 0, len(outjets))
 	for i := range outjets {
-		jet := &outjets[i]
-		area := fmom.PxPyPzE{0, 0, 0, 0}
+		var (
+			jet  = &outjets[i]
+			area fmom.PxPyPzE
+		)
 		if tsk.areaDef != nil {
 			// FIXME
 			panic("not implemented")
@@ -187,7 +189,7 @@ func (tsk *FastJetFinder) Process(ctx fwk.Context) error {
 			cand.Add(cst)
 		}
 
-		cand.Pos[3] = time / wtime
+		cand.Pos.P4.T = time / wtime
 		cand.Area = area
 		cand.DEta = detaMax
 		cand.DPhi = dphiMax
diff --git a/fads/utils.go b/fads/utils.go
index 65f301628..6d0054799 100644
--- a/fads/utils.go
+++ b/fads/utils.go
@@ -46,5 +46,5 @@ func newPtEtaPhiE(pt, eta, phi, ene float64) fmom.PxPyPzE {
 	py := pt * math.Sin(phi)
 	pz := pt * math.Sinh(eta)
 
-	return fmom.PxPyPzE{px, py, pz, ene}
+	return fmom.NewPxPyPzE(px, py, pz, ene)
 }
diff --git a/fmom/eetaphim.go b/fmom/eetaphim.go
index 794ef394f..df61c0977 100644
--- a/fmom/eetaphim.go
+++ b/fmom/eetaphim.go
@@ -8,10 +8,12 @@ import (
 	"math"
 )
 
-type EEtaPhiM [4]float64
+type EEtaPhiM struct {
+	P4 Vec4
+}
 
-func NewEEtaPhiM(et, eta, phi, m float64) EEtaPhiM {
-	return EEtaPhiM([4]float64{et, eta, phi, m})
+func NewEEtaPhiM(e, eta, phi, m float64) EEtaPhiM {
+	return EEtaPhiM{P4: Vec4{X: e, Y: eta, Z: phi, T: m}}
 }
 
 func (p4 *EEtaPhiM) Clone() P4 {
@@ -20,19 +22,19 @@ func (p4 *EEtaPhiM) Clone() P4 {
 }
 
 func (p4 *EEtaPhiM) E() float64 {
-	return p4[0]
+	return p4.P4.X
 }
 
 func (p4 *EEtaPhiM) Eta() float64 {
-	return p4[1]
+	return p4.P4.Y
 }
 
 func (p4 *EEtaPhiM) Phi() float64 {
-	return p4[2]
+	return p4.P4.Z
 }
 
 func (p4 *EEtaPhiM) M() float64 {
-	return p4[3]
+	return p4.P4.T
 }
 
 func (p4 *EEtaPhiM) M2() float64 {
@@ -143,8 +145,8 @@ func (p4 *EEtaPhiM) Pz() float64 {
 }
 
 func (p4 *EEtaPhiM) Set(p P4) {
-	p4[0] = p.E()
-	p4[1] = p.Eta()
-	p4[2] = p.Phi()
-	p4[3] = p.M()
+	p4.P4.X = p.E()
+	p4.P4.Y = p.Eta()
+	p4.P4.Z = p.Phi()
+	p4.P4.T = p.M()
 }
diff --git a/fmom/etetaphim.go b/fmom/etetaphim.go
index f74b44568..18845e124 100644
--- a/fmom/etetaphim.go
+++ b/fmom/etetaphim.go
@@ -8,10 +8,12 @@ import (
 	"math"
 )
 
-type EtEtaPhiM [4]float64
+type EtEtaPhiM struct {
+	P4 Vec4
+}
 
 func NewEtEtaPhiM(et, eta, phi, m float64) EtEtaPhiM {
-	return EtEtaPhiM([4]float64{et, eta, phi, m})
+	return EtEtaPhiM{P4: Vec4{X: et, Y: eta, Z: phi, T: m}}
 }
 
 func (p4 *EtEtaPhiM) Clone() P4 {
@@ -20,19 +22,19 @@ func (p4 *EtEtaPhiM) Clone() P4 {
 }
 
 func (p4 *EtEtaPhiM) Et() float64 {
-	return p4[0]
+	return p4.P4.X
 }
 
 func (p4 *EtEtaPhiM) Eta() float64 {
-	return p4[1]
+	return p4.P4.Y
 }
 
 func (p4 *EtEtaPhiM) Phi() float64 {
-	return p4[2]
+	return p4.P4.Z
 }
 
 func (p4 *EtEtaPhiM) M() float64 {
-	return p4[3]
+	return p4.P4.T
 }
 
 func (p4 *EtEtaPhiM) M2() float64 {
@@ -143,8 +145,8 @@ func (p4 *EtEtaPhiM) Pz() float64 {
 }
 
 func (p4 *EtEtaPhiM) Set(p P4) {
-	p4[0] = p.Et()
-	p4[1] = p.Eta()
-	p4[2] = p.Phi()
-	p4[3] = p.M()
+	p4.P4.X = p.Et()
+	p4.P4.Y = p.Eta()
+	p4.P4.Z = p.Phi()
+	p4.P4.T = p.M()
 }
diff --git a/fmom/iptcotthphim.go b/fmom/iptcotthphim.go
index 6aeffaa2e..dc61181a0 100644
--- a/fmom/iptcotthphim.go
+++ b/fmom/iptcotthphim.go
@@ -8,10 +8,12 @@ import (
 	"math"
 )
 
-type IPtCotThPhiM [4]float64
+type IPtCotThPhiM struct {
+	P4 Vec4
+}
 
 func NewIPtCotThPhiM(pt, eta, phi, m float64) IPtCotThPhiM {
-	return IPtCotThPhiM([4]float64{pt, eta, phi, m})
+	return IPtCotThPhiM{P4: Vec4{X: pt, Y: eta, Z: phi, T: m}}
 }
 
 func (p4 *IPtCotThPhiM) Clone() P4 {
@@ -20,19 +22,19 @@ func (p4 *IPtCotThPhiM) Clone() P4 {
 }
 
 func (p4 *IPtCotThPhiM) IPt() float64 {
-	return p4[0]
+	return p4.P4.X
 }
 
 func (p4 *IPtCotThPhiM) CotTh() float64 {
-	return p4[1]
+	return p4.P4.Y
 }
 
 func (p4 *IPtCotThPhiM) Phi() float64 {
-	return p4[2]
+	return p4.P4.Z
 }
 
 func (p4 *IPtCotThPhiM) M() float64 {
-	return p4[3]
+	return p4.P4.T
 }
 
 func (p4 *IPtCotThPhiM) Pt() float64 {
@@ -137,8 +139,8 @@ func (p4 *IPtCotThPhiM) SinPhi() float64 {
 }
 
 func (p4 *IPtCotThPhiM) Set(p P4) {
-	p4[0] = p.IPt()
-	p4[1] = p.CotTh()
-	p4[2] = p.Phi()
-	p4[3] = p.M()
+	p4.P4.X = p.IPt()
+	p4.P4.Y = p.CotTh()
+	p4.P4.Z = p.Phi()
+	p4.P4.T = p.M()
 }
diff --git a/fmom/ops.go b/fmom/ops.go
index 6f8900d19..2d3faa51b 100644
--- a/fmom/ops.go
+++ b/fmom/ops.go
@@ -110,10 +110,10 @@ func IAdd(dst, src P4) P4 {
 	default:
 		panic(fmt.Errorf("fmom: invalid P4 concrete value: %#v", dst))
 	}
-	p4[0] += src.Px()
-	p4[1] += src.Py()
-	p4[2] += src.Pz()
-	p4[3] += src.E()
+	p4.P4.X += src.Px()
+	p4.P4.Y += src.Py()
+	p4.P4.Z += src.Pz()
+	p4.P4.T += src.E()
 	sum.Set(p4)
 	return sum
 }
diff --git a/fmom/p4.go b/fmom/p4.go
index fc3c4b842..8326899bf 100644
--- a/fmom/p4.go
+++ b/fmom/p4.go
@@ -30,3 +30,8 @@ type P4 interface {
 	Set(p4 P4)
 	Clone() P4
 }
+
+// Vec4 holds the four components of a Lorentz vector.
+type Vec4 struct {
+	X, Y, Z, T float64
+}
diff --git a/fmom/ptetaphim.go b/fmom/ptetaphim.go
index 465484614..dfc49f8ce 100644
--- a/fmom/ptetaphim.go
+++ b/fmom/ptetaphim.go
@@ -8,10 +8,12 @@ import (
 	"math"
 )
 
-type PtEtaPhiM [4]float64
+type PtEtaPhiM struct {
+	P4 Vec4
+}
 
 func NewPtEtaPhiM(pt, eta, phi, m float64) PtEtaPhiM {
-	return PtEtaPhiM([4]float64{pt, eta, phi, m})
+	return PtEtaPhiM{P4: Vec4{X: pt, Y: eta, Z: phi, T: m}}
 }
 
 func (p4 *PtEtaPhiM) Clone() P4 {
@@ -20,19 +22,19 @@ func (p4 *PtEtaPhiM) Clone() P4 {
 }
 
 func (p4 *PtEtaPhiM) Pt() float64 {
-	return p4[0]
+	return p4.P4.X
 }
 
 func (p4 *PtEtaPhiM) Eta() float64 {
-	return p4[1]
+	return p4.P4.Y
 }
 
 func (p4 *PtEtaPhiM) Phi() float64 {
-	return p4[2]
+	return p4.P4.Z
 }
 
 func (p4 *PtEtaPhiM) M() float64 {
-	return p4[3]
+	return p4.P4.T
 }
 
 func (p4 *PtEtaPhiM) E() float64 {
@@ -150,8 +152,8 @@ func (p4 *PtEtaPhiM) Pz() float64 {
 }
 
 func (p4 *PtEtaPhiM) Set(p P4) {
-	p4[0] = p.Pt()
-	p4[1] = p.Eta()
-	p4[2] = p.Phi()
-	p4[3] = p.M()
+	p4.P4.X = p.Pt()
+	p4.P4.Y = p.Eta()
+	p4.P4.Z = p.Phi()
+	p4.P4.T = p.M()
 }
diff --git a/fmom/pxpypze.go b/fmom/pxpypze.go
index f1c01ef2f..91255787e 100644
--- a/fmom/pxpypze.go
+++ b/fmom/pxpypze.go
@@ -8,10 +8,12 @@ import (
 	"math"
 )
 
-type PxPyPzE [4]float64
+type PxPyPzE struct {
+	P4 Vec4
+}
 
 func NewPxPyPzE(px, py, pz, e float64) PxPyPzE {
-	return PxPyPzE([4]float64{px, py, pz, e})
+	return PxPyPzE{P4: Vec4{X: px, Y: py, Z: pz, T: e}}
 }
 
 func (p4 *PxPyPzE) Clone() P4 {
@@ -19,37 +21,15 @@ func (p4 *PxPyPzE) Clone() P4 {
 	return &pp
 }
 
-func (p4 *PxPyPzE) Px() float64 {
-	return p4[0]
-}
-
-func (p4 *PxPyPzE) Py() float64 {
-	return p4[1]
-}
-
-func (p4 *PxPyPzE) Pz() float64 {
-	return p4[2]
-}
-
-func (p4 *PxPyPzE) E() float64 {
-	return p4[3]
-}
+func (p4 *PxPyPzE) Px() float64 { return p4.P4.X }
+func (p4 *PxPyPzE) Py() float64 { return p4.P4.Y }
+func (p4 *PxPyPzE) Pz() float64 { return p4.P4.Z }
+func (p4 *PxPyPzE) E() float64  { return p4.P4.T }
 
-func (p4 *PxPyPzE) X() float64 {
-	return p4[0]
-}
-
-func (p4 *PxPyPzE) Y() float64 {
-	return p4[1]
-}
-
-func (p4 *PxPyPzE) Z() float64 {
-	return p4[2]
-}
-
-func (p4 *PxPyPzE) T() float64 {
-	return p4[3]
-}
+func (p4 *PxPyPzE) X() float64 { return p4.P4.X }
+func (p4 *PxPyPzE) Y() float64 { return p4.P4.Y }
+func (p4 *PxPyPzE) Z() float64 { return p4.P4.Z }
+func (p4 *PxPyPzE) T() float64 { return p4.P4.T }
 
 func (p4 *PxPyPzE) M2() float64 {
 	px := p4.Px()
@@ -206,8 +186,8 @@ func (p4 *PxPyPzE) Rapidity() float64 {
 }
 
 func (p4 *PxPyPzE) Set(p P4) {
-	p4[0] = p.Px()
-	p4[1] = p.Py()
-	p4[2] = p.Pz()
-	p4[3] = p.E()
+	p4.P4.X = p.Px()
+	p4.P4.Y = p.Py()
+	p4.P4.Z = p.Pz()
+	p4.P4.T = p.E()
 }
diff --git a/hepmc/decoder.go b/hepmc/decoder.go
index a846163b6..101d3041f 100644
--- a/hepmc/decoder.go
+++ b/hepmc/decoder.go
@@ -503,22 +503,22 @@ func (dec *Decoder) decodeVertex(evt *Event, vtx *Vertex, pidxToEndVtx map[int]i
 		return err
 	}
 
-	vtx.Position[0], err = tokens.float64()
+	vtx.Position.P4.X, err = tokens.float64()
 	if err != nil {
 		return err
 	}
 
-	vtx.Position[1], err = tokens.float64()
+	vtx.Position.P4.Y, err = tokens.float64()
 	if err != nil {
 		return err
 	}
 
-	vtx.Position[2], err = tokens.float64()
+	vtx.Position.P4.Z, err = tokens.float64()
 	if err != nil {
 		return err
 	}
 
-	vtx.Position[3], err = tokens.float64()
+	vtx.Position.P4.T, err = tokens.float64()
 	if err != nil {
 		return err
 	}
@@ -602,22 +602,22 @@ func (dec *Decoder) decodeParticle(evt *Event, p *Particle, pidxToEndVtx map[int
 		return err
 	}
 
-	p.Momentum[0], err = tokens.float64()
+	p.Momentum.P4.X, err = tokens.float64()
 	if err != nil {
 		return err
 	}
 
-	p.Momentum[1], err = tokens.float64()
+	p.Momentum.P4.Y, err = tokens.float64()
 	if err != nil {
 		return err
 	}
 
-	p.Momentum[2], err = tokens.float64()
+	p.Momentum.P4.Z, err = tokens.float64()
 	if err != nil {
 		return err
 	}
 
-	p.Momentum[3], err = tokens.float64()
+	p.Momentum.P4.T, err = tokens.float64()
 	if err != nil {
 		return err
 	}
diff --git a/hepmc/hepmc.go b/hepmc/hepmc.go
index a2b211b0e..cc7fcabb5 100644
--- a/hepmc/hepmc.go
+++ b/hepmc/hepmc.go
@@ -506,8 +506,10 @@ func (vtx *Vertex) removeParticleOut(p *Particle) error {
 
 // Print prints the vertex to w in a human-readable format
 func (vtx *Vertex) Print(w io.Writer) error {
-	var err error
-	zero := fmom.PxPyPzE{0, 0, 0, 0}
+	var (
+		err  error
+		zero fmom.PxPyPzE
+	)
 	if vtx.Barcode != 0 {
 		if vtx.Position != zero {
 			_, err = fmt.Fprintf(
diff --git a/hepmc/hepmc_test.go b/hepmc/hepmc_test.go
index 1f79c1a39..6691724b8 100644
--- a/hepmc/hepmc_test.go
+++ b/hepmc/hepmc_test.go
@@ -282,7 +282,7 @@ func ExampleEvent_buildFromScratch() {
 	}
 
 	err = v1.AddParticleIn(&hepmc.Particle{
-		Momentum: fmom.PxPyPzE{0, 0, 7000, 7000},
+		Momentum: fmom.NewPxPyPzE(0, 0, 7000, 7000),
 		PdgID:    2212,
 		Status:   3,
 	})
@@ -297,7 +297,7 @@ func ExampleEvent_buildFromScratch() {
 	}
 
 	err = v2.AddParticleIn(&hepmc.Particle{
-		Momentum: fmom.PxPyPzE{0, 0, -7000, 7000},
+		Momentum: fmom.NewPxPyPzE(0, 0, -7000, 7000),
 		PdgID:    2212,
 		Status:   3,
 		//Barcode:  2,
@@ -308,7 +308,7 @@ func ExampleEvent_buildFromScratch() {
 
 	// create the outgoing particles of v1 and v2
 	p3 := &hepmc.Particle{
-		Momentum: fmom.PxPyPzE{.750, -1.569, 32.191, 32.238},
+		Momentum: fmom.NewPxPyPzE(.750, -1.569, 32.191, 32.238),
 		PdgID:    1,
 		Status:   3,
 		// Barcode: 3,
@@ -319,7 +319,7 @@ func ExampleEvent_buildFromScratch() {
 	}
 
 	p4 := &hepmc.Particle{
-		Momentum: fmom.PxPyPzE{-3.047, -19., -54.629, 57.920},
+		Momentum: fmom.NewPxPyPzE(-3.047, -19., -54.629, 57.920),
 		PdgID:    -2,
 		Status:   3,
 		// Barcode: 4,
@@ -347,7 +347,7 @@ func ExampleEvent_buildFromScratch() {
 	}
 
 	err = v3.AddParticleOut(&hepmc.Particle{
-		Momentum: fmom.PxPyPzE{-3.813, 0.113, -1.833, 4.233},
+		Momentum: fmom.NewPxPyPzE(-3.813, 0.113, -1.833, 4.233),
 		PdgID:    22,
 		Status:   1,
 	})
@@ -356,7 +356,7 @@ func ExampleEvent_buildFromScratch() {
 	}
 
 	p5 := &hepmc.Particle{
-		Momentum: fmom.PxPyPzE{1.517, -20.68, -20.605, 85.925},
+		Momentum: fmom.NewPxPyPzE(1.517, -20.68, -20.605, 85.925),
 		PdgID:    -24,
 		Status:   3,
 	}
@@ -367,7 +367,7 @@ func ExampleEvent_buildFromScratch() {
 
 	// create v4
 	v4 := &hepmc.Vertex{
-		Position: fmom.PxPyPzE{0.12, -0.3, 0.05, 0.004},
+		Position: fmom.NewPxPyPzE(0.12, -0.3, 0.05, 0.004),
 	}
 	err = evt.AddVertex(v4)
 	if err != nil {
@@ -380,7 +380,7 @@ func ExampleEvent_buildFromScratch() {
 	}
 
 	err = v4.AddParticleOut(&hepmc.Particle{
-		Momentum: fmom.PxPyPzE{-2.445, 28.816, 6.082, 29.552},
+		Momentum: fmom.NewPxPyPzE(-2.445, 28.816, 6.082, 29.552),
 		PdgID:    1,
 		Status:   1,
 	})
@@ -389,7 +389,7 @@ func ExampleEvent_buildFromScratch() {
 	}
 
 	err = v4.AddParticleOut(&hepmc.Particle{
-		Momentum: fmom.PxPyPzE{3.962, -49.498, -26.687, 56.373},
+		Momentum: fmom.NewPxPyPzE(3.962, -49.498, -26.687, 56.373),
 		PdgID:    -2,
 		Status:   1,
 	})