Fork me on GitHub

source: git/external/fastjet/JetDefinition.cc@ 0e6b14e

ImprovedOutputFile Timing
Last change on this file since 0e6b14e was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 7 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 18.5 KB
Line 
1//FJSTARTHEADER
2// $Id: JetDefinition.cc 4354 2018-04-22 07:12:37Z salam $
3//
4// Copyright (c) 2005-2018, 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).
95bool JetDefinition::is_spherical() const {
96 if (jet_algorithm() == plugin_algorithm) {
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//----------------------------------------------------------------------
106string JetDefinition::description() const {
107 ostringstream name;
108
109 name << description_no_recombiner();
110
111 if ((jet_algorithm() == plugin_algorithm) || (jet_algorithm() == undefined_jet_algorithm)){
112 return name.str();
113 }
114
115 if (n_parameters_for_algorithm(jet_algorithm()) == 0)
116 name << " with ";
117 else
118 name << " and ";
119 name << recombiner()->description();
120
121 return name.str();
122}
123
124//----------------------------------------------------------------------
125string JetDefinition::description_no_recombiner() const {
126
127 ostringstream name;
128 if (jet_algorithm() == plugin_algorithm) {
129 return plugin()->description();
130 } else if (jet_algorithm() == undefined_jet_algorithm) {
131 return "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
132 }
133
134 name << algorithm_description(jet_algorithm());
135 switch (n_parameters_for_algorithm(jet_algorithm())){
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
142 if (jet_algorithm() == cambridge_for_passive_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//----------------------------------------------------------------------
154string JetDefinition::algorithm_description(const JetAlgorithm jet_alg){
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//----------------------------------------------------------------------
172unsigned int JetDefinition::n_parameters_for_algorithm(const JetAlgorithm jet_alg){
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//----------------------------------------------------------------------
182void JetDefinition::set_recombination_scheme(
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
192void JetDefinition::set_recombiner(const JetDefinition &other_jet_def){
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.reset(other_jet_def._shared_recombiner);
210
211 // NB: it is tempting to go via set_recombiner and then to sort
212 // out the shared part, but this would be dangerous in the
213 // specific (rare?) case where other_jet_def is the same as this
214 // it deletes_recombiner_when_unused. In that case the shared
215 // pointer reset would delete the recombiner.
216}
217
218
219// returns true if the current jet definitions shares the same
220// recombiner as teh one passed as an argument
221bool JetDefinition::has_same_recombiner(const JetDefinition &other_jd) const{
222 // first make sure that they have the same recombination scheme
223 const RecombinationScheme & scheme = recombination_scheme();
224 if (other_jd.recombination_scheme() != scheme) return false;
225
226 // if the scheme is "external", also check that they have the same
227 // recombiner
228 return (scheme != external_scheme)
229 || (recombiner() == other_jd.recombiner());
230}
231
232/// causes the JetDefinition to handle the deletion of the
233/// recombiner when it is no longer used
234void JetDefinition::delete_recombiner_when_unused(){
235 if (_recombiner == 0){
236 throw Error("tried to call JetDefinition::delete_recombiner_when_unused() for a JetDefinition without a user-defined recombination scheme");
237 } else if (_shared_recombiner.get()) {
238 throw Error("Error in JetDefinition::delete_recombiner_when_unused: the recombiner is already scheduled for deletion when unused (or was already set as shared)");
239 }
240
241 _shared_recombiner.reset(_recombiner);
242}
243
244/// allows to let the JetDefinition handle the deletion of the
245/// plugin when it is no longer used
246void JetDefinition::delete_plugin_when_unused(){
247 if (_plugin == 0){
248 throw Error("tried to call JetDefinition::delete_plugin_when_unused() for a JetDefinition without a plugin");
249 }
250
251 _plugin_shared.reset(_plugin);
252}
253
254
255
256string JetDefinition::DefaultRecombiner::description() const {
257 switch(_recomb_scheme) {
258 case E_scheme:
259 return "E scheme recombination";
260 case pt_scheme:
261 return "pt scheme recombination";
262 case pt2_scheme:
263 return "pt2 scheme recombination";
264 case Et_scheme:
265 return "Et scheme recombination";
266 case Et2_scheme:
267 return "Et2 scheme recombination";
268 case BIpt_scheme:
269 return "boost-invariant pt scheme recombination";
270 case BIpt2_scheme:
271 return "boost-invariant pt2 scheme recombination";
272 case WTA_pt_scheme:
273 return "pt-ordered Winner-Takes-All recombination";
274 // Energy-ordering can lead to dangerous situations with particles at
275 // rest. We instead implement the WTA_modp_scheme
276 //
277 // case WTA_E_scheme:
278 // return "energy-ordered Winner-Takes-All recombination";
279 case WTA_modp_scheme:
280 return "|3-momentum|-ordered Winner-Takes-All recombination";
281 default:
282 ostringstream err;
283 err << "DefaultRecombiner: unrecognized recombination scheme "
284 << _recomb_scheme;
285 throw Error(err.str());
286 }
287}
288
289
290void JetDefinition::DefaultRecombiner::recombine(
291 const PseudoJet & pa, const PseudoJet & pb,
292 PseudoJet & pab) const {
293
294 double weighta, weightb;
295
296 switch(_recomb_scheme) {
297 case E_scheme:
298 // a call to reset turns out to be somewhat more efficient
299 // than a sum and assignment
300 //pab = pa + pb;
301 pab.reset(pa.px()+pb.px(),
302 pa.py()+pb.py(),
303 pa.pz()+pb.pz(),
304 pa.E ()+pb.E ());
305 return;
306 // all remaining schemes are massless recombinations and locally
307 // we just set weights, while the hard work is done below...
308 case pt_scheme:
309 case Et_scheme:
310 case BIpt_scheme:
311 weighta = pa.perp();
312 weightb = pb.perp();
313 break;
314 case pt2_scheme:
315 case Et2_scheme:
316 case BIpt2_scheme:
317 weighta = pa.perp2();
318 weightb = pb.perp2();
319 break;
320 case WTA_pt_scheme:{
321 const PseudoJet & phard = (pa.pt2() >= pb.pt2()) ? pa : pb;
322 /// keep y,phi and m from the hardest, sum pt
323 pab.reset_PtYPhiM(pa.pt()+pb.pt(),
324 phard.rap(), phard.phi(), phard.m());
325 return;}
326 // Energy-ordering can lead to dangerous situations with particles at
327 // rest. We instead implement the WTA_modp_scheme
328 //
329 // case WTA_E_scheme:{
330 // const PseudoJet & phard = (pa.E() >= pb.E()) ? pa : pb;
331 // /// keep 3-momentum direction and mass from the hardest, sum energies
332 // ///
333 // /// If the particle with the largest energy is at rest, the sum
334 // /// remains at rest, implying that the mass of the sum is larger
335 // /// than the mass of pa.
336 // double Eab = pa.E() + pb.E();
337 // double scale = (phard.modp2()==0.0)
338 // ? 0.0
339 // : sqrt((Eab*Eab - phard.m2())/phard.modp2());
340 // pab.reset(phard.px()*scale, phard.py()*scale, phard.pz()*scale, Eab);
341 // return;}
342 case WTA_modp_scheme:{
343 // Note: we need to compute both a and b modp. And we need pthard
344 // and its modp. If we want to avoid repeating the test and do
345 // only 2 modp calculations, we'd have to duplicate the code (or
346 // use a pair<const PJ&>). An alternative is to write modp_soft as
347 // modp_ab-modp_hard but this could suffer from larger rounding
348 // errors
349 bool a_hardest = (pa.modp2() >= pb.modp2());
350 const PseudoJet & phard = a_hardest ? pa : pb;
351 const PseudoJet & psoft = a_hardest ? pb : pa;
352 /// keep 3-momentum direction and mass from the hardest, sum modp
353 ///
354 /// If the hardest particle is at rest, the sum remains at rest
355 /// (the energy of the sum is therefore the mass of pa)
356 double modp_hard = phard.modp();
357 double modp_ab = modp_hard + psoft.modp();
358 if (phard.modp2()==0.0){
359 pab.reset(0.0, 0.0, 0.0, phard.m());
360 } else {
361 double scale = modp_ab/modp_hard;
362 pab.reset(phard.px()*scale, phard.py()*scale, phard.pz()*scale,
363 sqrt(modp_ab*modp_ab + phard.m2()));
364 }
365 return;}
366 default:
367 ostringstream err;
368 err << "DefaultRecombiner: unrecognized recombination scheme "
369 << _recomb_scheme;
370 throw Error(err.str());
371 }
372
373 double perp_ab = pa.perp() + pb.perp();
374 if (perp_ab != 0.0) { // weights also non-zero...
375 double y_ab = (weighta * pa.rap() + weightb * pb.rap())/(weighta+weightb);
376
377 // take care with periodicity in phi...
378 double phi_a = pa.phi(), phi_b = pb.phi();
379 if (phi_a - phi_b > pi) phi_b += twopi;
380 if (phi_a - phi_b < -pi) phi_b -= twopi;
381 double phi_ab = (weighta * phi_a + weightb * phi_b)/(weighta+weightb);
382
383 // this is much more efficient...
384 pab.reset_PtYPhiM(perp_ab,y_ab,phi_ab);
385 // pab = PseudoJet(perp_ab*cos(phi_ab),
386 // perp_ab*sin(phi_ab),
387 // perp_ab*sinh(y_ab),
388 // perp_ab*cosh(y_ab));
389 } else { // weights are zero
390 //pab = PseudoJet(0.0,0.0,0.0,0.0);
391 pab.reset(0.0, 0.0, 0.0, 0.0);
392 }
393}
394
395
396void JetDefinition::DefaultRecombiner::preprocess(PseudoJet & p) const {
397 switch(_recomb_scheme) {
398 case E_scheme:
399 case BIpt_scheme:
400 case BIpt2_scheme:
401 case WTA_pt_scheme:
402 //case WTA_E_scheme:
403 case WTA_modp_scheme:
404 break;
405 case pt_scheme:
406 case pt2_scheme:
407 {
408 // these schemes (as in the ktjet implementation) need massless
409 // initial 4-vectors with essentially E=|p|.
410 double newE = sqrt(p.perp2()+p.pz()*p.pz());
411 p.reset_momentum(p.px(), p.py(), p.pz(), newE);
412 // FJ2.x version
413 // int user_index = p.user_index();
414 // p = PseudoJet(p.px(), p.py(), p.pz(), newE);
415 // p.set_user_index(user_index);
416 }
417 break;
418 case Et_scheme:
419 case Et2_scheme:
420 {
421 // these schemes (as in the ktjet implementation) need massless
422 // initial 4-vectors with essentially E=|p|.
423 double rescale = p.E()/sqrt(p.perp2()+p.pz()*p.pz());
424 p.reset_momentum(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
425 // FJ2.x version
426 // int user_index = p.user_index();
427 // p = PseudoJet(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
428 // p.set_user_index(user_index);
429 }
430 break;
431 default:
432 ostringstream err;
433 err << "DefaultRecombiner: unrecognized recombination scheme "
434 << _recomb_scheme;
435 throw Error(err.str());
436 }
437}
438
439void JetDefinition::Plugin::set_ghost_separation_scale(double /*scale*/) const {
440 throw Error("set_ghost_separation_scale not supported");
441}
442
443
444
445//-------------------------------------------------------------------------------
446// helper functions to build a jet made of pieces
447//
448// This is the extended version with support for a user-defined
449// recombination-scheme
450// -------------------------------------------------------------------------------
451
452// build a "CompositeJet" from the vector of its pieces
453//
454// the user passes the reciombination scheme used to "sum" the pieces.
455PseudoJet join(const vector<PseudoJet> & pieces, const JetDefinition::Recombiner & recombiner){
456 // compute the total momentum
457 //--------------------------------------------------
458 PseudoJet result; // automatically initialised to 0
459 if (pieces.size()>0){
460 result = pieces[0];
461 for (unsigned int i=1; i<pieces.size(); i++)
462 recombiner.plus_equal(result, pieces[i]);
463 }
464
465 // attach a CompositeJetStructure to the result
466 //--------------------------------------------------
467 CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces, &recombiner);
468
469 result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct));
470
471 return result;
472}
473
474// build a "CompositeJet" from a single PseudoJet
475PseudoJet join(const PseudoJet & j1,
476 const JetDefinition::Recombiner & recombiner){
477 return join(vector<PseudoJet>(1,j1), recombiner);
478}
479
480// build a "CompositeJet" from two PseudoJet
481PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
482 const JetDefinition::Recombiner & recombiner){
483 vector<PseudoJet> pieces;
484 pieces.push_back(j1);
485 pieces.push_back(j2);
486 return join(pieces, recombiner);
487}
488
489// build a "CompositeJet" from 3 PseudoJet
490PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3,
491 const JetDefinition::Recombiner & recombiner){
492 vector<PseudoJet> pieces;
493 pieces.push_back(j1);
494 pieces.push_back(j2);
495 pieces.push_back(j3);
496 return join(pieces, recombiner);
497}
498
499// build a "CompositeJet" from 4 PseudoJet
500PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4,
501 const JetDefinition::Recombiner & recombiner){
502 vector<PseudoJet> pieces;
503 pieces.push_back(j1);
504 pieces.push_back(j2);
505 pieces.push_back(j3);
506 pieces.push_back(j4);
507 return join(pieces, recombiner);
508}
509
510
511
512
513FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.