FastJet  3.0.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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