FastJet  3.0.6
JetDefinition.cc
1 //STARTHEADER
2 // $Id: JetDefinition.cc 2712 2011-11-16 22:38:44Z salam $
3 //
4 // Copyright (c) 2005-2011, 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 and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 #include "fastjet/JetDefinition.hh"
30 #include "fastjet/Error.hh"
31 #include "fastjet/CompositeJetStructure.hh"
32 #include<sstream>
33 
34 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
35 
36 using namespace std;
37 
38 const double JetDefinition::max_allowable_R = 1000.0;
39 
40 //----------------------------------------------------------------------
41 // [NB: implementation was getting complex, so in 2.4-devel moved it
42 // from .hh to .cc]
43 JetDefinition::JetDefinition(JetAlgorithm jet_algorithm_in,
44  double R_in,
45  Strategy strategy_in,
46  RecombinationScheme recomb_scheme_in,
47  int nparameters) :
48  _jet_algorithm(jet_algorithm_in), _Rparam(R_in), _strategy(strategy_in) {
49 
50  // set R parameter or ensure its sensibleness, as appropriate
51  if (_jet_algorithm == ee_kt_algorithm) {
52  _Rparam = 4.0; // introduce a fictional R that ensures that
53  // our clustering sequence will not produce
54  // "beam" jets except when only a single particle remains.
55  // Any value > 2 would have done here
56  } else {
57  // We maintain some limit on R because particles with pt=0, m=0
58  // can have rapidities O(100000) and one doesn't want the
59  // clustering to start including them as if their rapidities were
60  // physical.
61  if (R_in > max_allowable_R) {
62  ostringstream oss;
63  oss << "Requested R = " << R_in << " for jet definition is larger than max_allowable_R = " << max_allowable_R;
64  throw Error(oss.str());
65  }
66  }
67 
68  // cross-check the number of parameters that were declared in setting up the
69  // algorithm (passed internally from the public constructors)
70  switch (jet_algorithm_in) {
71  case ee_kt_algorithm:
72  if (nparameters != 0) {
73  ostringstream oss;
74  oss << "ee_kt_algorithm should be constructed with 0 parameters but was called with "
75  << nparameters << " parameter(s)\n";
76  throw Error(oss.str());
77  }
78  break;
79  case genkt_algorithm:
80  case ee_genkt_algorithm:
81  if (nparameters != 2) {
82  ostringstream oss;
83  oss << "(ee_)genkt_algorithm should be constructed with 2 parameters but was called with "
84  << nparameters << " parameter(s)\n";
85  throw Error(oss.str());
86  }
87  break;
88  default:
89  if (nparameters != 1) {
90  ostringstream oss;
91  oss << "The jet algorithm you requested ("
92  << jet_algorithm_in << ") should be constructed with 1 parameter but was called with "
93  << nparameters << " parameter(s)\n";
94  throw Error(oss.str());
95  }
96  }
97 
98  // make sure the strategy requested is sensible
99  assert (_strategy != plugin_strategy);
100 
101  _plugin = NULL;
102  set_recombination_scheme(recomb_scheme_in);
103  set_extra_param(0.0); // make sure it's defined
104 }
105 
106 
107 //----------------------------------------------------------------------
109  ostringstream name;
110  if (jet_algorithm() == plugin_algorithm) {
111  return plugin()->description();
112  } else if (jet_algorithm() == kt_algorithm) {
113  name << "Longitudinally invariant kt algorithm with R = " << R();
114  name << " and " << recombiner()->description();
115  } else if (jet_algorithm() == cambridge_algorithm) {
116  name << "Longitudinally invariant Cambridge/Aachen algorithm with R = "
117  << R() ;
118  name << " and " << recombiner()->description();
119  } else if (jet_algorithm() == antikt_algorithm) {
120  name << "Longitudinally invariant anti-kt algorithm with R = "
121  << R() ;
122  name << " and " << recombiner()->description();
123  } else if (jet_algorithm() == genkt_algorithm) {
124  name << "Longitudinally invariant generalised kt algorithm with R = "
125  << R() << ", p = " << extra_param();
126  name << " and " << recombiner()->description();
128  name << "Longitudinally invariant Cambridge/Aachen algorithm with R = "
129  << R() << "and a special hack whereby particles with kt < "
130  << extra_param() << "are treated as passive ghosts";
131  } else if (jet_algorithm() == ee_kt_algorithm) {
132  name << "e+e- kt (Durham) algorithm (NB: no R)";
133  name << " with " << recombiner()->description();
134  } else if (jet_algorithm() == ee_genkt_algorithm) {
135  name << "e+e- generalised kt algorithm with R = "
136  << R() << ", p = " << extra_param();
137  name << " and " << recombiner()->description();
138  } else if (jet_algorithm() == undefined_jet_algorithm) {
139  name << "uninitialised JetDefinition (jet_algorithm=undefined_jet_algorithm)" ;
140  } else {
141  throw Error("JetDefinition::description(): unrecognized jet_algorithm");
142  }
143  return name.str();
144 }
145 
146 
148  RecombinationScheme recomb_scheme) {
149  _default_recombiner = JetDefinition::DefaultRecombiner(recomb_scheme);
150 
151  // do not forget to delete the existing recombiner if needed
152  if (_recombiner_shared()) _recombiner_shared.reset();
153 
154  _recombiner = 0;
155 }
156 
157 
158 // returns true if the current jet definitions shares the same
159 // recombiner as teh one passed as an argument
161  // first make sure that they have the same recombination scheme
162  const RecombinationScheme & scheme = recombination_scheme();
163  if (other_jd.recombination_scheme() != scheme) return false;
164 
165  // if the scheme is "external", also check that they ahve the same
166  // recombiner
167  return (scheme != external_scheme)
168  || (recombiner() == other_jd.recombiner());
169 }
170 
171 /// allows to let the JetDefinition handle the deletion of the
172 /// recombiner when it is no longer used
174  if (_recombiner == 0){
175  throw Error("tried to call JetDefinition::delete_recombiner_when_unused() for a JetDefinition without a user-defined recombination scheme");
176  }
177 
178  _recombiner_shared.reset(_recombiner);
179 }
180 
181 /// allows to let the JetDefinition handle the deletion of the
182 /// plugin when it is no longer used
184  if (_plugin == 0){
185  throw Error("tried to call JetDefinition::delete_plugin_when_unused() for a JetDefinition without a plugin");
186  }
187 
188  _plugin_shared.reset(_plugin);
189 }
190 
191 
192 
194  switch(_recomb_scheme) {
195  case E_scheme:
196  return "E scheme recombination";
197  case pt_scheme:
198  return "pt scheme recombination";
199  case pt2_scheme:
200  return "pt2 scheme recombination";
201  case Et_scheme:
202  return "Et scheme recombination";
203  case Et2_scheme:
204  return "Et2 scheme recombination";
205  case BIpt_scheme:
206  return "boost-invariant pt scheme recombination";
207  case BIpt2_scheme:
208  return "boost-invariant pt2 scheme recombination";
209  default:
210  ostringstream err;
211  err << "DefaultRecombiner: unrecognized recombination scheme "
212  << _recomb_scheme;
213  throw Error(err.str());
214  }
215 }
216 
217 
219  const PseudoJet & pa, const PseudoJet & pb,
220  PseudoJet & pab) const {
221 
222  double weighta, weightb;
223 
224  switch(_recomb_scheme) {
225  case E_scheme:
226  // a call to reset turns out to be somewhat more efficient
227  // than a sum and assignment
228  //pab = pa + pb;
229  pab.reset(pa.px()+pb.px(),
230  pa.py()+pb.py(),
231  pa.pz()+pb.pz(),
232  pa.E ()+pb.E ());
233  return;
234  // all remaining schemes are massless recombinations and locally
235  // we just set weights, while the hard work is done below...
236  case pt_scheme:
237  case Et_scheme:
238  case BIpt_scheme:
239  weighta = pa.perp();
240  weightb = pb.perp();
241  break;
242  case pt2_scheme:
243  case Et2_scheme:
244  case BIpt2_scheme:
245  weighta = pa.perp2();
246  weightb = pb.perp2();
247  break;
248  default:
249  ostringstream err;
250  err << "DefaultRecombiner: unrecognized recombination scheme "
251  << _recomb_scheme;
252  throw Error(err.str());
253  }
254 
255  double perp_ab = pa.perp() + pb.perp();
256  if (perp_ab != 0.0) { // weights also non-zero...
257  double y_ab = (weighta * pa.rap() + weightb * pb.rap())/(weighta+weightb);
258 
259  // take care with periodicity in phi...
260  double phi_a = pa.phi(), phi_b = pb.phi();
261  if (phi_a - phi_b > pi) phi_b += twopi;
262  if (phi_a - phi_b < -pi) phi_b -= twopi;
263  double phi_ab = (weighta * phi_a + weightb * phi_b)/(weighta+weightb);
264 
265  // this is much more efficient...
266  pab.reset_PtYPhiM(perp_ab,y_ab,phi_ab);
267  // pab = PseudoJet(perp_ab*cos(phi_ab),
268  // perp_ab*sin(phi_ab),
269  // perp_ab*sinh(y_ab),
270  // perp_ab*cosh(y_ab));
271  } else { // weights are zero
272  //pab = PseudoJet(0.0,0.0,0.0,0.0);
273  pab.reset(0.0, 0.0, 0.0, 0.0);
274  }
275 }
276 
277 
279  switch(_recomb_scheme) {
280  case E_scheme:
281  case BIpt_scheme:
282  case BIpt2_scheme:
283  break;
284  case pt_scheme:
285  case pt2_scheme:
286  {
287  // these schemes (as in the ktjet implementation) need massless
288  // initial 4-vectors with essentially E=|p|.
289  double newE = sqrt(p.perp2()+p.pz()*p.pz());
290  p.reset_momentum(p.px(), p.py(), p.pz(), newE);
291  // FJ2.x version
292  // int user_index = p.user_index();
293  // p = PseudoJet(p.px(), p.py(), p.pz(), newE);
294  // p.set_user_index(user_index);
295  }
296  break;
297  case Et_scheme:
298  case Et2_scheme:
299  {
300  // these schemes (as in the ktjet implementation) need massless
301  // initial 4-vectors with essentially E=|p|.
302  double rescale = p.E()/sqrt(p.perp2()+p.pz()*p.pz());
303  p.reset_momentum(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
304  // FJ2.x version
305  // int user_index = p.user_index();
306  // p = PseudoJet(rescale*p.px(), rescale*p.py(), rescale*p.pz(), p.E());
307  // p.set_user_index(user_index);
308  }
309  break;
310  default:
311  ostringstream err;
312  err << "DefaultRecombiner: unrecognized recombination scheme "
313  << _recomb_scheme;
314  throw Error(err.str());
315  }
316 }
317 
319  throw Error("set_ghost_separation_scale not supported");
320 }
321 
322 
323 
324 //-------------------------------------------------------------------------------
325 // helper functions to build a jet made of pieces
326 //
327 // This is the extended version with support for a user-defined
328 // recombination-scheme
329 // -------------------------------------------------------------------------------
330 
331 // build a "CompositeJet" from the vector of its pieces
332 //
333 // the user passes the reciombination scheme used to "sum" the pieces.
334 PseudoJet join(const vector<PseudoJet> & pieces, const JetDefinition::Recombiner & recombiner){
335  // compute the total momentum
336  //--------------------------------------------------
337  PseudoJet result; // automatically initialised to 0
338  if (pieces.size()>0){
339  result = pieces[0];
340  for (unsigned int i=1; i<pieces.size(); i++)
341  recombiner.plus_equal(result, pieces[i]);
342  }
343 
344  // attach a CompositeJetStructure to the result
345  //--------------------------------------------------
346  CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces, &recombiner);
347 
349 
350  return result;
351 }
352 
353 // build a "CompositeJet" from a single PseudoJet
354 PseudoJet join(const PseudoJet & j1,
355  const JetDefinition::Recombiner & recombiner){
356  return join(vector<PseudoJet>(1,j1), recombiner);
357 }
358 
359 // build a "CompositeJet" from two PseudoJet
360 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2,
361  const JetDefinition::Recombiner & recombiner){
362  vector<PseudoJet> pieces;
363  pieces.push_back(j1);
364  pieces.push_back(j2);
365  return join(pieces, recombiner);
366 }
367 
368 // build a "CompositeJet" from 3 PseudoJet
369 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3,
370  const JetDefinition::Recombiner & recombiner){
371  vector<PseudoJet> pieces;
372  pieces.push_back(j1);
373  pieces.push_back(j2);
374  pieces.push_back(j3);
375  return join(pieces, recombiner);
376 }
377 
378 // build a "CompositeJet" from 4 PseudoJet
379 PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4,
380  const JetDefinition::Recombiner & recombiner){
381  vector<PseudoJet> pieces;
382  pieces.push_back(j1);
383  pieces.push_back(j2);
384  pieces.push_back(j3);
385  pieces.push_back(j4);
386  return join(pieces, recombiner);
387 }
388 
389 
390 
391 
392 FASTJET_END_NAMESPACE
pt weighted recombination of y,phi (and summing of pt's), with no preprocessing
void set_recombination_scheme(RecombinationScheme)
set the recombination scheme to the one provided
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:121
a version of cambridge with a special distance measure for particles whose pt is < extra_param() ...
like the k_t but with distance measures dij = min(kti^{2p},ktj^{2p}) Delta R_{ij}^2 / R^2 diB = 1/kti...
void delete_recombiner_when_unused()
calling this tells the JetDefinition to handle the deletion of the recombiner when it is no longer us...
summing the 4-momenta
the value for the jet algorithm in a JetDefinition for which no algorithm has yet been defined ...
pt^2 weighted recombination of y,phi (and summing of pt's) no preprocessing
JetAlgorithm jet_algorithm() const
return information about the definition...
void delete_plugin_when_unused()
allows to let the JetDefinition handle the deletion of the plugin when it is no longer used ...
pt^2 weighted recombination of y,phi (and summing of pt's) with preprocessing to make things massless...
the e+e- genkt algorithm (R > 2 and p=1 gives ee_kt)
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:274
the longitudinally invariant kt algorithm
The structure for a jet made of pieces.
const Plugin * plugin() const
return a pointer to the plugin
void set_extra_param(double xtra_param)
(re)set the general purpose extra parameter
any plugin algorithm supplied by the user
the plugin has been used...
A class that will provide the recombination scheme facilities and/or allow a user to extend these fac...
virtual std::string description() const
return a textual description of the recombination scheme implemented here
std::string description() const
return a textual description of the current jet definition
pt^2 weighted recombination of y,phi (and summing of pt's) with preprocessing to make things massless...
virtual void preprocess(PseudoJet &p) const
routine called to preprocess each input jet (to make all input jets compatible with the scheme requir...
virtual std::string description() const =0
return a textual description of the jet-definition implemented in this plugin
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:41
RecombinationScheme
the various recombination schemes
virtual void set_ghost_separation_scale(double scale) const
set the ghost separation scale for passive area determinations in future runs (strictly speaking that...
void plus_equal(PseudoJet &pa, const PseudoJet &pb) const
pa += pb in the given recombination scheme.
static const double max_allowable_R
R values larger than max_allowable_R are not allowed.
an implementation of C++0x shared pointers (or boost's)
Definition: SharedPtr.hh:114
pt weighted recombination of y,phi (and summing of pt's) with preprocessing to make things massless b...
virtual std::string description() const =0
return a textual description of the recombination scheme implemented here
pt weighted recombination of y,phi (and summing of pt's) with preprocessing to make things massless b...
const Recombiner * recombiner() const
return a pointer to the currently defined recombiner.
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...
Definition: PseudoJet.hh:915
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:106
virtual void recombine(const PseudoJet &pa, const PseudoJet &pb, PseudoJet &pab) const
recombine pa and pb and put result into pab
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...
the e+e- kt algorithm
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure)
set the associated structure
Definition: PseudoJet.cc:450
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:141
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
An abstract base class that will provide the recombination scheme facilities and/or allow a user to e...
JetAlgorithm
the various families of jet-clustering algorithm
Strategy
the various options for the algorithmic strategy to adopt in clustering events with kt and cambridge ...
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:65
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...
Definition: PseudoJet.hh:924
for the user's external scheme
double perp2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:139
class that is intended to hold a full definition of the jet clusterer
bool has_same_recombiner(const JetDefinition &other_jd) const
returns true if the current jet definitions shares the same recombiner as teh one passed as an argume...