FastJet 3.4.3
Loading...
Searching...
No Matches
JetDefinition.cc
1//FJSTARTHEADER
2// $Id$
3//
4// Copyright (c) 2005-2024, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5//
6//----------------------------------------------------------------------
7// This file is part of FastJet.
8//
9// FastJet is free software; you can redistribute it and/or modify
10// it under the terms of the GNU General Public License as published by
11// the Free Software Foundation; either version 2 of the License, or
12// (at your option) any later version.
13//
14// The algorithms that underlie FastJet have required considerable
15// development. They are described in the original FastJet paper,
16// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17// FastJet as part of work towards a scientific publication, please
18// quote the version you use and include a citation to the manual and
19// optionally also to hep-ph/0512210.
20//
21// FastJet is distributed in the hope that it will be useful,
22// but WITHOUT ANY WARRANTY; without even the implied warranty of
23// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24// GNU General Public License for more details.
25//
26// You should have received a copy of the GNU General Public License
27// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28//----------------------------------------------------------------------
29//FJENDHEADER
30
31#include "fastjet/JetDefinition.hh"
32#include "fastjet/Error.hh"
33#include "fastjet/CompositeJetStructure.hh"
34#include<sstream>
35
36FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
37
38using namespace std;
39
40const double JetDefinition::max_allowable_R = 1000.0;
41
42//----------------------------------------------------------------------
43// [NB: implementation was getting complex, so in 2.4-devel moved it
44// from .hh to .cc]
45JetDefinition::JetDefinition(JetAlgorithm jet_algorithm_in,
46 double R_in,
47 RecombinationScheme recomb_scheme_in,
48 Strategy strategy_in,
49 int nparameters) :
50 _jet_algorithm(jet_algorithm_in), _Rparam(R_in), _strategy(strategy_in) {
51
52 // set R parameter or ensure its sensibleness, as appropriate
53 if (_jet_algorithm == ee_kt_algorithm) {
54 _Rparam = 4.0; // introduce a fictional R that ensures that
55 // our clustering sequence will not produce
56 // "beam" jets except when only a single particle remains.
57 // Any value > 2 would have done here
58 } else {
59 // We maintain some limit on R because particles with pt=0, m=0
60 // can have rapidities O(100000) and one doesn't want the
61 // clustering to start including them as if their rapidities were
62 // physical.
63 if (R_in > max_allowable_R) {
64 ostringstream oss;
65 oss << "Requested R = " << R_in << " for jet definition is larger than max_allowable_R = " << max_allowable_R;
66 throw Error(oss.str());
67 }
68 }
69
70 // cross-check the number of parameters that were declared in setting up the
71 // algorithm (passed internally from the public constructors)
72 unsigned int nparameters_expected = n_parameters_for_algorithm(jet_algorithm_in);
73 if (nparameters != (int) nparameters_expected){
74 ostringstream oss;
75 oss << "The jet algorithm you requested ("
76 << jet_algorithm_in << ") should be constructed with " << nparameters_expected
77 << " parameter(s) but was called with " << nparameters << " parameter(s)\n";
78 throw Error(oss.str());
79 }
80
81 // make sure the strategy requested is sensible
82 assert (_strategy != plugin_strategy);
83
84 _plugin = NULL;
85 set_recombination_scheme(recomb_scheme_in);
86 set_extra_param(0.0); // make sure it's defined
87}
88
89
90//----------------------------------------------------------------------
91// returns true if the jet definition involves an algorithm
92// intended for use on a spherical geometry (e.g. e+e- algorithms,
93// as opposed to most pp algorithms, which use a cylindrical,
94// rapidity-phi geometry).
97 return plugin()->is_spherical();
98 } else {
99 return (jet_algorithm() == ee_kt_algorithm || // as of 2013-02-14, the two
100 jet_algorithm() == ee_genkt_algorithm // native spherical algorithms
101 );
102 }
103}
104
105//----------------------------------------------------------------------
107 ostringstream name;
108
110
112 return name.str();
113 }
114
116 name << " with ";
117 else
118 name << " and ";
119 name << recombiner()->description();
120
121 return name.str();
122}
123
124//----------------------------------------------------------------------
126
127 ostringstream name;
129 return plugin()->description();
130 } else if (jet_algorithm() == undefined_jet_algorithm) {
131 return "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
132 }
133
136 case 0: name << " (NB: no R)"; break;
137 case 1: name << " with R = " << R(); break; // the parameter is always R
138 case 2:
139 // the 1st parameter is always R
140 name << " with R = " << R();
141 // the 2nd depends on the algorithm
143 name << "and a special hack whereby particles with kt < "
144 << extra_param() << "are treated as passive ghosts";
145 } else {
146 name << ", p = " << extra_param();
147 }
148 };
149
150 return name.str();
151}
152
153//----------------------------------------------------------------------
155 ostringstream name;
156 switch (jet_alg){
157 case plugin_algorithm: return "plugin algorithm";
158 case kt_algorithm: return "Longitudinally invariant kt algorithm";
159 case cambridge_algorithm: return "Longitudinally invariant Cambridge/Aachen algorithm";
160 case antikt_algorithm: return "Longitudinally invariant anti-kt algorithm";
161 case genkt_algorithm: return "Longitudinally invariant generalised kt algorithm";
162 case cambridge_for_passive_algorithm: return "Longitudinally invariant Cambridge/Aachen algorithm";
163 case ee_kt_algorithm: return "e+e- kt (Durham) algorithm (NB: no R)";
164 case ee_genkt_algorithm: return "e+e- generalised kt algorithm";
165 case undefined_jet_algorithm: return "undefined jet algorithm";
166 default:
167 throw Error("JetDefinition::algorithm_description(): unrecognized jet_algorithm");
168 };
169}
170
171//----------------------------------------------------------------------
173 switch (jet_alg) {
174 case ee_kt_algorithm: return 0;
175 case genkt_algorithm:
176 case ee_genkt_algorithm: return 2;
177 default: return 1;
178 };
179}
180
181//----------------------------------------------------------------------
183 RecombinationScheme recomb_scheme) {
184 _default_recombiner = JetDefinition::DefaultRecombiner(recomb_scheme);
185
186 // do not forget to delete the existing recombiner if needed
187 if (_shared_recombiner) _shared_recombiner.reset();
188
189 _recombiner = 0;
190}
191
193 // make sure the "invariants" of the other jet def are sensible
194 assert(other_jet_def._recombiner ||
195 other_jet_def.recombination_scheme() != external_scheme);
196
197 // first treat the situation where we're using the default recombiner
198 if (other_jet_def._recombiner == 0){
199 set_recombination_scheme(other_jet_def.recombination_scheme());
200 return;
201 }
202
203 // in other cases, copy the pointer to the recombiner
204 _recombiner = other_jet_def._recombiner;
205 // set the default recombiner appropriately
206 _default_recombiner = DefaultRecombiner(external_scheme);
207 // and set the _shared_recombiner to the same state
208 // as in the other_jet_def, whatever that was
209 _shared_recombiner = other_jet_def._shared_recombiner;
210 //_shared_recombiner.reset(other_jet_def._shared_recombiner);
211
212 // NB: it is tempting to go via set_recombiner and then to sort
213 // out the shared part, but this would be dangerous in the
214 // specific (rare?) case where other_jet_def is the same as this
215 // it deletes_recombiner_when_unused. In that case the shared
216 // pointer reset would delete the recombiner.
217}
218
219
220// returns true if the current jet definitions shares the same
221// recombiner as teh one passed as an argument
223 // first make sure that they have the same recombination scheme
224 const RecombinationScheme & scheme = recombination_scheme();
225 if (other_jd.recombination_scheme() != scheme) return false;
226
227 // if the scheme is "external", also check that they have the same
228 // recombiner
229 return (scheme != external_scheme)
230 || (recombiner() == other_jd.recombiner());
231}
232
233/// causes the JetDefinition to handle the deletion of the
234/// recombiner when it is no longer used
236 if (_recombiner == 0){
237 throw Error("tried to call JetDefinition::delete_recombiner_when_unused() for a JetDefinition without a user-defined recombination scheme");
238 } else if (_shared_recombiner.get()) {
239 throw Error("Error in JetDefinition::delete_recombiner_when_unused: the recombiner is already scheduled for deletion when unused (or was already set as shared)");
240 }
241
242 _shared_recombiner.reset(_recombiner);
243}
244
245/// allows to let the JetDefinition handle the deletion of the
246/// plugin when it is no longer used
248 if (_plugin == 0){
249 throw Error("tried to call JetDefinition::delete_plugin_when_unused() for a JetDefinition without a plugin");
250 }
251
252 _plugin_shared.reset(_plugin);
253}
254
255
256
258 switch(_recomb_scheme) {
259 case E_scheme:
260 return "E scheme recombination";
261 case pt_scheme:
262 return "pt scheme recombination";
263 case pt2_scheme:
264 return "pt2 scheme recombination";
265 case Et_scheme:
266 return "Et scheme recombination";
267 case Et2_scheme:
268 return "Et2 scheme recombination";
269 case BIpt_scheme:
270 return "boost-invariant pt scheme recombination";
271 case BIpt2_scheme:
272 return "boost-invariant pt2 scheme recombination";
273 case WTA_pt_scheme:
274 return "pt-ordered Winner-Takes-All recombination";
275 // Energy-ordering can lead to dangerous situations with particles at
276 // rest. We instead implement the WTA_modp_scheme
277 //
278 // case WTA_E_scheme:
279 // return "energy-ordered Winner-Takes-All recombination";
280 case WTA_modp_scheme:
281 return "|3-momentum|-ordered Winner-Takes-All recombination";
282 default:
283 ostringstream err;
284 err << "DefaultRecombiner: unrecognized recombination scheme "
285 << _recomb_scheme;
286 throw Error(err.str());
287 }
288}
289
290
292 const PseudoJet & pa, const PseudoJet & pb,
293 PseudoJet & pab) const {
294
295 double weighta, weightb;
296
297 switch(_recomb_scheme) {
298 case E_scheme:
299 // a call to reset turns out to be somewhat more efficient
300 // than a sum and assignment
301 //pab = pa + pb;
302 pab.reset(pa.px()+pb.px(),
303 pa.py()+pb.py(),
304 pa.pz()+pb.pz(),
305 pa.E ()+pb.E ());
306 return;
307 // all remaining schemes are massless recombinations and locally
308 // we just set weights, while the hard work is done below...
309 case pt_scheme:
310 case Et_scheme:
311 case BIpt_scheme:
312 weighta = pa.perp();
313 weightb = pb.perp();
314 break;
315 case pt2_scheme:
316 case Et2_scheme:
317 case BIpt2_scheme:
318 weighta = pa.perp2();
319 weightb = pb.perp2();
320 break;
321 case WTA_pt_scheme:{
322 const PseudoJet & phard = (pa.pt2() >= pb.pt2()) ? pa : pb;
323 /// keep y,phi and m from the hardest, sum pt
324 pab.reset_PtYPhiM(pa.pt()+pb.pt(),
325 phard.rap(), phard.phi(), phard.m());
326 return;}
327 // Energy-ordering can lead to dangerous situations with particles at
328 // rest. We instead implement the WTA_modp_scheme
329 //
330 // case WTA_E_scheme:{
331 // const PseudoJet & phard = (pa.E() >= pb.E()) ? pa : pb;
332 // /// keep 3-momentum direction and mass from the hardest, sum energies
333 // ///
334 // /// If the particle with the largest energy is at rest, the sum
335 // /// remains at rest, implying that the mass of the sum is larger
336 // /// than the mass of pa.
337 // double Eab = pa.E() + pb.E();
338 // double scale = (phard.modp2()==0.0)
339 // ? 0.0
340 // : sqrt((Eab*Eab - phard.m2())/phard.modp2());
341 // pab.reset(phard.px()*scale, phard.py()*scale, phard.pz()*scale, Eab);
342 // return;}
343 case WTA_modp_scheme:{
344 // Note: we need to compute both a and b modp. And we need pthard
345 // and its modp. If we want to avoid repeating the test and do
346 // only 2 modp calculations, we'd have to duplicate the code (or
347 // use a pair<const PJ&>). An alternative is to write modp_soft as
348 // modp_ab-modp_hard but this could suffer from larger rounding
349 // errors
350 bool a_hardest = (pa.modp2() >= pb.modp2());
351 const PseudoJet & phard = a_hardest ? pa : pb;
352 const PseudoJet & psoft = a_hardest ? pb : pa;
353 /// keep 3-momentum direction and mass from the hardest, sum modp
354 ///
355 /// If the hardest particle is at rest, the sum remains at rest
356 /// (the energy of the sum is therefore the mass of pa)
357 double modp_hard = phard.modp();
358 double modp_ab = modp_hard + psoft.modp();
359 if (phard.modp2()==0.0){
360 pab.reset(0.0, 0.0, 0.0, phard.m());
361 } else {
362 double scale = modp_ab/modp_hard;
363 pab.reset(phard.px()*scale, phard.py()*scale, phard.pz()*scale,
364 sqrt(modp_ab*modp_ab + phard.m2()));
365 }
366 return;}
367 default:
368 ostringstream err;
369 err << "DefaultRecombiner: unrecognized recombination scheme "
370 << _recomb_scheme;
371 throw Error(err.str());
372 }
373
374 double perp_ab = pa.perp() + pb.perp();
375 if (perp_ab != 0.0) { // weights also non-zero...
376 double y_ab = (weighta * pa.rap() + weightb * pb.rap())/(weighta+weightb);
377
378 // take care with periodicity in phi...
379 double phi_a = pa.phi(), phi_b = pb.phi();
380 if (phi_a - phi_b > pi) phi_b += twopi;
381 if (phi_a - phi_b < -pi) phi_b -= twopi;
382 double phi_ab = (weighta * phi_a + weightb * phi_b)/(weighta+weightb);
383
384 // this is much more efficient...
385 pab.reset_PtYPhiM(perp_ab,y_ab,phi_ab);
386 // pab = PseudoJet(perp_ab*cos(phi_ab),
387 // perp_ab*sin(phi_ab),
388 // perp_ab*sinh(y_ab),
389 // perp_ab*cosh(y_ab));
390 } else { // weights are zero
391 //pab = PseudoJet(0.0,0.0,0.0,0.0);
392 pab.reset(0.0, 0.0, 0.0, 0.0);
393 }
394}
395
396
398 switch(_recomb_scheme) {
399 case E_scheme:
400 case BIpt_scheme:
401 case BIpt2_scheme:
402 case WTA_pt_scheme:
403 //case WTA_E_scheme:
404 case WTA_modp_scheme:
405 break;
406 case pt_scheme:
407 case pt2_scheme:
408 {
409 // these schemes (as in the ktjet implementation) need massless
410 // initial 4-vectors with essentially E=|p|.
411 double newE = sqrt(p.perp2()+p.pz()*p.pz());
412 p.reset_momentum(p.px(), p.py(), p.pz(), newE);
413 // FJ2.x version
414 // int user_index = p.user_index();
415 // p = PseudoJet(p.px(), p.py(), p.pz(), newE);
416 // p.set_user_index(user_index);
417 }
418 break;
419 case Et_scheme:
420 case Et2_scheme:
421 {
422 // these schemes (as in the ktjet implementation) need massless
423 // initial 4-vectors with essentially E=|p|.
424 double rescale = p.E()/sqrt(p.perp2()+p.pz()*p.pz());
425 p.reset_momentum(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
426 // FJ2.x version
427 // int user_index = p.user_index();
428 // p = PseudoJet(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
429 // p.set_user_index(user_index);
430 }
431 break;
432 default:
433 ostringstream err;
434 err << "DefaultRecombiner: unrecognized recombination scheme "
435 << _recomb_scheme;
436 throw Error(err.str());
437 }
438}
439
441 throw Error("set_ghost_separation_scale not supported");
442}
443
444
445
446//-------------------------------------------------------------------------------
447// helper functions to build a jet made of pieces
448//
449// This is the extended version with support for a user-defined
450// recombination-scheme
451// -------------------------------------------------------------------------------
452
453// build a "CompositeJet" from the vector of its pieces
454//
455// the user passes the reciombination scheme used to "sum" the pieces.
456PseudoJet join(const vector<PseudoJet> & pieces, const JetDefinition::Recombiner & recombiner){
457 // compute the total momentum
458 //--------------------------------------------------
459 PseudoJet result; // automatically initialised to 0
460 if (pieces.size()>0){
461 result = pieces[0];
462 for (unsigned int i=1; i<pieces.size(); i++)
463 recombiner.plus_equal(result, pieces[i]);
464 }
465
466 // attach a CompositeJetStructure to the result
467 //--------------------------------------------------
468 CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces, &recombiner);
469
471
472 return result;
473}
474
475// build a "CompositeJet" from a single PseudoJet
476PseudoJet join(const PseudoJet & j1,
478 return join(vector<PseudoJet>(1,j1), recombiner);
479}
480
481// build a "CompositeJet" from two PseudoJet
482PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
484 vector<PseudoJet> pieces;
485 pieces.push_back(j1);
486 pieces.push_back(j2);
487 return join(pieces, recombiner);
488}
489
490// build a "CompositeJet" from 3 PseudoJet
491PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3,
493 vector<PseudoJet> pieces;
494 pieces.push_back(j1);
495 pieces.push_back(j2);
496 pieces.push_back(j3);
497 return join(pieces, recombiner);
498}
499
500// build a "CompositeJet" from 4 PseudoJet
501PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4,
503 vector<PseudoJet> pieces;
504 pieces.push_back(j1);
505 pieces.push_back(j2);
506 pieces.push_back(j3);
507 pieces.push_back(j4);
508 return join(pieces, recombiner);
509}
510
511
512
513
514FASTJET_END_NAMESPACE
The structure for a jet made of pieces.
base class corresponding to errors that can be thrown by FastJet
Definition Error.hh:52
A class that will provide the recombination scheme facilities and/or allow a user to extend these fac...
virtual void recombine(const PseudoJet &pa, const PseudoJet &pb, PseudoJet &pab) const override
recombine pa and pb and put result into pab
virtual std::string description() const override
return a textual description of the recombination scheme implemented here
virtual void preprocess(PseudoJet &p) const override
routine called to preprocess each input jet (to make all input jets compatible with the scheme requir...
virtual void set_ghost_separation_scale(double scale) const
set the ghost separation scale for passive area determinations in future runs (strictly speaking that...
virtual bool is_spherical() const
returns true if the plugin implements an algorithm intended for use on a spherical geometry (e....
virtual std::string description() const =0
return a textual description of the jet-definition implemented in this plugin
An abstract base class that will provide the recombination scheme facilities and/or allow a user to e...
virtual std::string description() const =0
return a textual description of the recombination scheme implemented here
void plus_equal(PseudoJet &pa, const PseudoJet &pb) const
pa += pb in the given recombination scheme.
class that is intended to hold a full definition of the jet clusterer
std::string description() const
return a textual description of the current jet definition
bool is_spherical() const
returns true if the jet definition involves an algorithm intended for use on a spherical geometry (e....
const Plugin * plugin() const
return a pointer to the plugin
JetAlgorithm jet_algorithm() const
return information about the definition...
static std::string algorithm_description(const JetAlgorithm jet_alg)
a short textual description of the algorithm jet_alg
const Recombiner * recombiner() const
returns a pointer to the currently defined recombiner.
static unsigned int n_parameters_for_algorithm(const JetAlgorithm jet_alg)
the number of parameters associated to a given jet algorithm
void delete_recombiner_when_unused()
calling this tells the JetDefinition to handle the deletion of the recombiner when it is no longer us...
void set_extra_param(double xtra_param)
(re)set the general purpose extra parameter
bool has_same_recombiner(const JetDefinition &other_jd) const
returns true if the current jet definitions shares the same recombiner as the one passed as an argume...
std::string description_no_recombiner() const
returns a description not including the recombiner information
static const double max_allowable_R
R values larger than max_allowable_R are not allowed.
void set_recombiner(const Recombiner *recomb)
set the recombiner class to the one provided
void set_recombination_scheme(RecombinationScheme)
set the recombination scheme to the one provided
void delete_plugin_when_unused()
calling this causes the JetDefinition to handle the deletion of the plugin when it is no longer used
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition PseudoJet.hh:68
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition PseudoJet.hh:138
double modp() const
return the 3-vector modulus = sqrt(px^2+py^2+pz^2)
Definition PseudoJet.hh:180
void reset_momentum(double px, double py, double pz, double E)
reset the 4-momentum according to the supplied components but leave all other information (indices,...
double phi() const
returns phi (in the range 0..2pi)
Definition PseudoJet.hh:123
double perp() const
returns the scalar transverse momentum
Definition PseudoJet.hh:158
double perp2() const
returns the squared transverse momentum
Definition PseudoJet.hh:156
double m2() const
returns the squared invariant mass // like CLHEP
Definition PseudoJet.hh:163
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure_in)
set the associated structure
Definition PseudoJet.cc:561
void reset_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in=0.0)
reset the PseudoJet according to the specified pt, rapidity, azimuth and mass (also resetting indices...
Definition PseudoJet.hh:333
double pt() const
returns the scalar transverse momentum
Definition PseudoJet.hh:154
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP)
double modp2() const
return the squared 3-vector modulus = px^2+py^2+pz^2
Definition PseudoJet.hh:178
void reset(double px, double py, double pz, double E)
reset the 4-momentum according to the supplied components and put the user and history indices back t...
double pt2() const
returns the squared transverse momentum
Definition PseudoJet.hh:152
An implementation of shared pointers that is broadly similar to C++11 shared_ptr (https://en....
Definition SharedPtr.hh:341
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
@ plugin_strategy
the plugin has been used...
RecombinationScheme
The various recombination schemes.
@ E_scheme
summing the 4-momenta
@ BIpt2_scheme
pt^2 weighted recombination of y,phi (and summing of pt's) no preprocessing
@ pt2_scheme
pt^2 weighted recombination of y,phi (and summing of pt's) with preprocessing to make things massless...
@ pt_scheme
pt weighted recombination of y,phi (and summing of pt's) with preprocessing to make things massless b...
@ external_scheme
for the user's external scheme
@ Et_scheme
pt weighted recombination of y,phi (and summing of pt's) with preprocessing to make things massless b...
@ WTA_pt_scheme
pt-based Winner-Takes-All (WTA) recombination: the result of the recombination has the rapidity,...
@ WTA_modp_scheme
mod-p-based Winner-Takes-All (WTA) recombination: the result of the recombination gets the 3-vector d...
@ BIpt_scheme
pt weighted recombination of y,phi (and summing of pt's), with no preprocessing
@ Et2_scheme
pt^2 weighted recombination of y,phi (and summing of pt's) with preprocessing to make things massless...
JetAlgorithm
the various families of jet-clustering algorithm
@ ee_genkt_algorithm
the e+e- genkt algorithm (R > 2 and p=1 gives ee_kt)
@ genkt_algorithm
like the k_t but with distance measures dij = min(kti^{2p},ktj^{2p}) Delta R_{ij}^2 / R^2 diB = 1/kti...
@ undefined_jet_algorithm
the value for the jet algorithm in a JetDefinition for which no algorithm has yet been defined
@ plugin_algorithm
any plugin algorithm supplied by the user
@ ee_kt_algorithm
the e+e- kt algorithm
@ cambridge_algorithm
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
@ cambridge_for_passive_algorithm
a version of cambridge with a special distance measure for particles whose pt is < extra_param(); thi...
@ antikt_algorithm
like the k_t but with distance measures dij = min(1/kti^2,1/ktj^2) Delta R_{ij}^2 / R^2 diB = 1/kti^2
@ kt_algorithm
the longitudinally invariant kt algorithm