go-hep.org/x/hep@v0.38.1/fastjet/jet.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 fastjet
     6  
     7  import (
     8  	"math"
     9  
    10  	"go-hep.org/x/hep/fmom"
    11  )
    12  
    13  const (
    14  	// Used to protect against parton-level events where pt can be zero
    15  	// for some partons, giving rapidity=infinity. KtJet fails in those cases.
    16  	MaxRap = 1e5
    17  )
    18  
    19  // UserInfo holds extra user information in a Jet
    20  type UserInfo any
    21  
    22  // Jet holds minimal information of use for jet-clustering routines
    23  type Jet struct {
    24  	fmom.PxPyPzE
    25  
    26  	UserInfo  UserInfo // holds extra user information for this Jet
    27  	hidx      int      // cluster sequence history index
    28  	structure JetStructure
    29  
    30  	// -- cache
    31  
    32  	rap float64
    33  	pt2 float64
    34  	phi float64
    35  }
    36  
    37  func NewJet(px, py, pz, e float64) Jet {
    38  	jet := Jet{
    39  		PxPyPzE: fmom.NewPxPyPzE(px, py, pz, e),
    40  		hidx:    -1,
    41  	}
    42  	jet.setupCache()
    43  	return jet
    44  }
    45  
    46  func (jet *Jet) setupCache() {
    47  	pt := jet.Pt()
    48  	jet.pt2 = pt * pt
    49  
    50  	var rap float64
    51  	if jet.E() == math.Abs(jet.Pz()) && jet.Pt2() == 0 {
    52  		rap = MaxRap + math.Abs(jet.Pz())
    53  		if jet.Pz() < 0 {
    54  			rap = -rap
    55  		}
    56  	} else {
    57  		m := jet.M()
    58  		m2 := math.Max(0, m*m) // effective mass - force non-tachyonic mass
    59  		e := jet.E() + math.Abs(jet.Pz())
    60  		rap = 0.5 * math.Log((jet.Pt2()+m2)/(e*e))
    61  		if jet.Pz() > 0 {
    62  			rap = -rap
    63  		}
    64  	}
    65  	jet.rap = rap
    66  	jet.phi = jet.PxPyPzE.Phi()
    67  }
    68  
    69  func (jet *Jet) Pt2() float64 {
    70  	return jet.pt2
    71  }
    72  
    73  func (jet *Jet) Phi() float64 {
    74  	return jet.phi
    75  }
    76  
    77  func (jet *Jet) Rapidity() float64 {
    78  	return jet.rap
    79  }
    80  
    81  // Constituents returns the list of constituents for this jet.
    82  func (jet *Jet) Constituents() []Jet {
    83  	subjets, err := jet.structure.Constituents(jet)
    84  	if err != nil {
    85  		panic(err)
    86  	}
    87  	return subjets
    88  }
    89  
    90  // Distance returns the squared cylinder (rapidity-phi) distance between 2 jets
    91  func Distance(j1, j2 *Jet) float64 {
    92  	dphi := deltaPhi(j1, j2)
    93  	drap := deltaRap(j1, j2)
    94  	return dphi*dphi + drap*drap
    95  }
    96  
    97  func deltaPhi(j1, j2 *Jet) float64 {
    98  	dphi := math.Abs(j1.Phi() - j2.Phi())
    99  	if dphi > math.Pi {
   100  		dphi = 2*math.Pi - dphi
   101  	}
   102  	return dphi
   103  }
   104  
   105  func deltaRap(j1, j2 *Jet) float64 {
   106  	return j1.Rapidity() - j2.Rapidity()
   107  }