go-hep.org/x/hep@v0.38.1/fastjet/recombiner.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  	"fmt"
     9  	"math"
    10  
    11  	"go-hep.org/x/hep/fmom"
    12  )
    13  
    14  type Recombiner interface {
    15  	Description() string
    16  	Recombine(j1, j2 *Jet) (Jet, error)
    17  	Preprocess(jet *Jet) error
    18  	Scheme() RecombinationScheme
    19  }
    20  
    21  type DefaultRecombiner struct {
    22  	scheme RecombinationScheme
    23  }
    24  
    25  func NewRecombiner(scheme RecombinationScheme) DefaultRecombiner {
    26  	return DefaultRecombiner{
    27  		scheme: scheme,
    28  	}
    29  }
    30  
    31  func (rec DefaultRecombiner) Description() string {
    32  	str := rec.scheme.String()
    33  	return str + " scheme recombination"
    34  }
    35  
    36  func (rec DefaultRecombiner) Recombine(j1, j2 *Jet) (Jet, error) {
    37  	w1 := 0.0
    38  	w2 := 0.0
    39  
    40  	switch rec.Scheme() {
    41  	case EScheme:
    42  		return NewJet(
    43  			j1.Px()+j2.Px(),
    44  			j1.Py()+j2.Py(),
    45  			j1.Pz()+j2.Pz(),
    46  			j1.E()+j2.E(),
    47  		), nil
    48  
    49  	case PtScheme, EtScheme, BIPtScheme:
    50  		w1 = j1.Pt()
    51  		w2 = j2.Pt()
    52  
    53  	case Pt2Scheme, Et2Scheme, BIPt2Scheme:
    54  		w1 = j1.Pt2()
    55  		w2 = j2.Pt2()
    56  
    57  	default:
    58  		return Jet{}, fmt.Errorf("fastjet.Recombine: invalid recombination scheme (%v)", rec.Scheme())
    59  	}
    60  
    61  	pt := j1.Pt() + j2.Pt()
    62  	if pt != 0.0 {
    63  		y := (w1*j1.Rapidity() + w2*j2.Rapidity()) / (w1 + w2)
    64  		phi1 := j1.Phi()
    65  		phi2 := j2.Phi()
    66  		if phi1-phi2 > math.Pi {
    67  			phi2 += 2 * math.Pi
    68  		}
    69  		if phi1-phi2 < -math.Pi {
    70  			phi2 -= 2 * math.Pi
    71  		}
    72  		phi := (w1*phi1 + w2*phi2) / (w1 + w2)
    73  		return NewJet(
    74  			pt*math.Cos(phi),
    75  			pt*math.Sin(phi),
    76  			pt*math.Sinh(y),
    77  			pt*math.Cosh(y),
    78  		), nil
    79  	}
    80  
    81  	return NewJet(0, 0, 0, 0), nil
    82  }
    83  
    84  func (rec DefaultRecombiner) Preprocess(jet *Jet) error {
    85  
    86  	switch rec.Scheme() {
    87  	case EScheme, BIPtScheme, BIPt2Scheme:
    88  		return nil
    89  
    90  	case PtScheme, Pt2Scheme:
    91  		// these schemes (as in the ktjet impl.) need massless
    92  		// initial 4-vectors, with essentially E=|p|
    93  		pz := jet.Pz()
    94  		e := math.Sqrt(jet.Pt2() + pz*pz)
    95  		jet.PxPyPzE = fmom.NewPxPyPzE(jet.Px(), jet.Py(), jet.Pz(), e)
    96  		return nil
    97  
    98  	case EtScheme, Et2Scheme:
    99  		// these schemes (as in the ktjet impl.) need massless
   100  		// initial 4-vectors, with essentially E=|p|
   101  		pz := jet.Pz()
   102  		rescale := jet.E() / math.Sqrt(jet.Pt2()+pz*pz)
   103  		jet.PxPyPzE = fmom.NewPxPyPzE(
   104  			rescale*jet.Px(),
   105  			rescale*jet.Py(),
   106  			rescale*jet.Pz(),
   107  			jet.E(),
   108  		)
   109  		return nil
   110  
   111  	default:
   112  		return fmt.Errorf("fastjet.Preprocess: invalid recombination scheme (%v)", rec.Scheme())
   113  	}
   114  }
   115  
   116  func (rec DefaultRecombiner) Scheme() RecombinationScheme {
   117  	return rec.scheme
   118  }