1 | //STARTHEADER |
---|
2 | // $Id: PseudoJet.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
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 | |
---|
30 | #include "fastjet/Error.hh" |
---|
31 | #include "fastjet/PseudoJet.hh" |
---|
32 | #include "fastjet/ClusterSequence.hh" |
---|
33 | #include "fastjet/ClusterSequenceAreaBase.hh" |
---|
34 | #include "fastjet/CompositeJetStructure.hh" |
---|
35 | #include<valarray> |
---|
36 | #include<iostream> |
---|
37 | #include<sstream> |
---|
38 | #include<cmath> |
---|
39 | #include<algorithm> |
---|
40 | #include <cstdarg> |
---|
41 | |
---|
42 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
43 | |
---|
44 | using namespace std; |
---|
45 | |
---|
46 | |
---|
47 | //---------------------------------------------------------------------- |
---|
48 | // another constructor... |
---|
49 | PseudoJet::PseudoJet(const double px_in, const double py_in, const double pz_in, const double E_in) { |
---|
50 | |
---|
51 | _E = E_in ; |
---|
52 | _px = px_in; |
---|
53 | _py = py_in; |
---|
54 | _pz = pz_in; |
---|
55 | |
---|
56 | this->_finish_init(); |
---|
57 | |
---|
58 | // some default values for the history and user indices |
---|
59 | _reset_indices(); |
---|
60 | |
---|
61 | } |
---|
62 | |
---|
63 | |
---|
64 | //---------------------------------------------------------------------- |
---|
65 | /// do standard end of initialisation |
---|
66 | void PseudoJet::_finish_init () { |
---|
67 | _kt2 = this->px()*this->px() + this->py()*this->py(); |
---|
68 | _phi = pseudojet_invalid_phi; |
---|
69 | } |
---|
70 | |
---|
71 | //---------------------------------------------------------------------- |
---|
72 | void PseudoJet::_set_rap_phi() const { |
---|
73 | |
---|
74 | if (_kt2 == 0.0) { |
---|
75 | _phi = 0.0; } |
---|
76 | else { |
---|
77 | _phi = atan2(this->py(),this->px()); |
---|
78 | } |
---|
79 | if (_phi < 0.0) {_phi += twopi;} |
---|
80 | if (_phi >= twopi) {_phi -= twopi;} // can happen if phi=-|eps<1e-15|? |
---|
81 | if (this->E() == abs(this->pz()) && _kt2 == 0) { |
---|
82 | // Point has infinite rapidity -- convert that into a very large |
---|
83 | // number, but in such a way that different 0-pt momenta will have |
---|
84 | // different rapidities (so as to lift the degeneracy between |
---|
85 | // them) [this can be relevant at parton-level] |
---|
86 | double MaxRapHere = MaxRap + abs(this->pz()); |
---|
87 | if (this->pz() >= 0.0) {_rap = MaxRapHere;} else {_rap = -MaxRapHere;} |
---|
88 | } else { |
---|
89 | // get the rapidity in a way that's modestly insensitive to roundoff |
---|
90 | // error when things pz,E are large (actually the best we can do without |
---|
91 | // explicit knowledge of mass) |
---|
92 | double effective_m2 = max(0.0,m2()); // force non tachyonic mass |
---|
93 | double E_plus_pz = _E + abs(_pz); // the safer of p+, p- |
---|
94 | // p+/p- = (p+ p-) / (p-)^2 = (kt^2+m^2)/(p-)^2 |
---|
95 | _rap = 0.5*log((_kt2 + effective_m2)/(E_plus_pz*E_plus_pz)); |
---|
96 | if (_pz > 0) {_rap = - _rap;} |
---|
97 | } |
---|
98 | |
---|
99 | } |
---|
100 | |
---|
101 | |
---|
102 | //---------------------------------------------------------------------- |
---|
103 | // return a valarray four-momentum |
---|
104 | valarray<double> PseudoJet::four_mom() const { |
---|
105 | valarray<double> mom(4); |
---|
106 | mom[0] = _px; |
---|
107 | mom[1] = _py; |
---|
108 | mom[2] = _pz; |
---|
109 | mom[3] = _E ; |
---|
110 | return mom; |
---|
111 | } |
---|
112 | |
---|
113 | //---------------------------------------------------------------------- |
---|
114 | // Return the component corresponding to the specified index. |
---|
115 | // taken from CLHEP |
---|
116 | double PseudoJet::operator () (int i) const { |
---|
117 | switch(i) { |
---|
118 | case X: |
---|
119 | return px(); |
---|
120 | case Y: |
---|
121 | return py(); |
---|
122 | case Z: |
---|
123 | return pz(); |
---|
124 | case T: |
---|
125 | return e(); |
---|
126 | default: |
---|
127 | ostringstream err; |
---|
128 | err << "PseudoJet subscripting: bad index (" << i << ")"; |
---|
129 | throw Error(err.str()); |
---|
130 | } |
---|
131 | return 0.; |
---|
132 | } |
---|
133 | |
---|
134 | //---------------------------------------------------------------------- |
---|
135 | // return the pseudorapidity |
---|
136 | double PseudoJet::pseudorapidity() const { |
---|
137 | if (px() == 0.0 && py() ==0.0) return MaxRap; |
---|
138 | if (pz() == 0.0) return 0.0; |
---|
139 | |
---|
140 | double theta = atan(perp()/pz()); |
---|
141 | if (theta < 0) theta += pi; |
---|
142 | return -log(tan(theta/2)); |
---|
143 | } |
---|
144 | |
---|
145 | //---------------------------------------------------------------------- |
---|
146 | // return "sum" of two pseudojets |
---|
147 | PseudoJet operator+ (const PseudoJet & jet1, const PseudoJet & jet2) { |
---|
148 | //return PseudoJet(jet1.four_mom()+jet2.four_mom()); |
---|
149 | return PseudoJet(jet1.px()+jet2.px(), |
---|
150 | jet1.py()+jet2.py(), |
---|
151 | jet1.pz()+jet2.pz(), |
---|
152 | jet1.E() +jet2.E() ); |
---|
153 | } |
---|
154 | |
---|
155 | //---------------------------------------------------------------------- |
---|
156 | // return difference of two pseudojets |
---|
157 | PseudoJet operator- (const PseudoJet & jet1, const PseudoJet & jet2) { |
---|
158 | //return PseudoJet(jet1.four_mom()-jet2.four_mom()); |
---|
159 | return PseudoJet(jet1.px()-jet2.px(), |
---|
160 | jet1.py()-jet2.py(), |
---|
161 | jet1.pz()-jet2.pz(), |
---|
162 | jet1.E() -jet2.E() ); |
---|
163 | } |
---|
164 | |
---|
165 | //---------------------------------------------------------------------- |
---|
166 | // return the product, coeff * jet |
---|
167 | PseudoJet operator* (double coeff, const PseudoJet & jet) { |
---|
168 | //return PseudoJet(coeff*jet.four_mom()); |
---|
169 | // the following code is hopefully more efficient |
---|
170 | PseudoJet coeff_times_jet(jet); |
---|
171 | coeff_times_jet *= coeff; |
---|
172 | return coeff_times_jet; |
---|
173 | } |
---|
174 | |
---|
175 | //---------------------------------------------------------------------- |
---|
176 | // return the product, coeff * jet |
---|
177 | PseudoJet operator* (const PseudoJet & jet, double coeff) { |
---|
178 | return coeff*jet; |
---|
179 | } |
---|
180 | |
---|
181 | //---------------------------------------------------------------------- |
---|
182 | // return the ratio, jet / coeff |
---|
183 | PseudoJet operator/ (const PseudoJet & jet, double coeff) { |
---|
184 | return (1.0/coeff)*jet; |
---|
185 | } |
---|
186 | |
---|
187 | //---------------------------------------------------------------------- |
---|
188 | /// multiply the jet's momentum by the coefficient |
---|
189 | void PseudoJet::operator*=(double coeff) { |
---|
190 | _px *= coeff; |
---|
191 | _py *= coeff; |
---|
192 | _pz *= coeff; |
---|
193 | _E *= coeff; |
---|
194 | _kt2*= coeff*coeff; |
---|
195 | // phi and rap are unchanged |
---|
196 | } |
---|
197 | |
---|
198 | //---------------------------------------------------------------------- |
---|
199 | /// divide the jet's momentum by the coefficient |
---|
200 | void PseudoJet::operator/=(double coeff) { |
---|
201 | (*this) *= 1.0/coeff; |
---|
202 | } |
---|
203 | |
---|
204 | |
---|
205 | //---------------------------------------------------------------------- |
---|
206 | /// add the other jet's momentum to this jet |
---|
207 | void PseudoJet::operator+=(const PseudoJet & other_jet) { |
---|
208 | _px += other_jet._px; |
---|
209 | _py += other_jet._py; |
---|
210 | _pz += other_jet._pz; |
---|
211 | _E += other_jet._E ; |
---|
212 | _finish_init(); // we need to recalculate phi,rap,kt2 |
---|
213 | } |
---|
214 | |
---|
215 | |
---|
216 | //---------------------------------------------------------------------- |
---|
217 | /// subtract the other jet's momentum from this jet |
---|
218 | void PseudoJet::operator-=(const PseudoJet & other_jet) { |
---|
219 | _px -= other_jet._px; |
---|
220 | _py -= other_jet._py; |
---|
221 | _pz -= other_jet._pz; |
---|
222 | _E -= other_jet._E ; |
---|
223 | _finish_init(); // we need to recalculate phi,rap,kt2 |
---|
224 | } |
---|
225 | |
---|
226 | //---------------------------------------------------------------------- |
---|
227 | bool operator==(const PseudoJet & a, const PseudoJet & b) { |
---|
228 | if (a.px() != b.px()) return false; |
---|
229 | if (a.py() != b.py()) return false; |
---|
230 | if (a.pz() != b.pz()) return false; |
---|
231 | if (a.E () != b.E ()) return false; |
---|
232 | |
---|
233 | if (a.user_index() != b.user_index()) return false; |
---|
234 | if (a.cluster_hist_index() != b.cluster_hist_index()) return false; |
---|
235 | if (a.user_info_ptr() != b.user_info_ptr()) return false; |
---|
236 | if (a.structure_ptr() != b.structure_ptr()) return false; |
---|
237 | |
---|
238 | return true; |
---|
239 | } |
---|
240 | |
---|
241 | //---------------------------------------------------------------------- |
---|
242 | // check if the jet has zero momentum |
---|
243 | bool operator==(const PseudoJet & jet, const double val) { |
---|
244 | if (val != 0) |
---|
245 | throw Error("comparing a PseudoJet with a non-zero constant (double) is not allowed."); |
---|
246 | return (jet.px() == 0 && jet.py() == 0 && |
---|
247 | jet.pz() == 0 && jet.E() == 0); |
---|
248 | } |
---|
249 | |
---|
250 | |
---|
251 | |
---|
252 | //---------------------------------------------------------------------- |
---|
253 | /// transform this jet (given in lab) into a jet in the rest |
---|
254 | /// frame of prest |
---|
255 | // |
---|
256 | // NB: code adapted from that in herwig f77 (checked how it worked |
---|
257 | // long ago) |
---|
258 | PseudoJet & PseudoJet::boost(const PseudoJet & prest) { |
---|
259 | |
---|
260 | if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0) |
---|
261 | return *this; |
---|
262 | |
---|
263 | double m_local = prest.m(); |
---|
264 | assert(m_local != 0); |
---|
265 | |
---|
266 | double pf4 = ( px()*prest.px() + py()*prest.py() |
---|
267 | + pz()*prest.pz() + E()*prest.E() )/m_local; |
---|
268 | double fn = (pf4 + E()) / (prest.E() + m_local); |
---|
269 | _px += fn*prest.px(); |
---|
270 | _py += fn*prest.py(); |
---|
271 | _pz += fn*prest.pz(); |
---|
272 | _E = pf4; |
---|
273 | |
---|
274 | _finish_init(); // we need to recalculate phi,rap,kt2 |
---|
275 | return *this; |
---|
276 | } |
---|
277 | |
---|
278 | |
---|
279 | //---------------------------------------------------------------------- |
---|
280 | /// transform this jet (given in the rest frame of prest) into a jet |
---|
281 | /// in the lab frame; |
---|
282 | // |
---|
283 | // NB: code adapted from that in herwig f77 (checked how it worked |
---|
284 | // long ago) |
---|
285 | PseudoJet & PseudoJet::unboost(const PseudoJet & prest) { |
---|
286 | |
---|
287 | if (prest.px() == 0.0 && prest.py() == 0.0 && prest.pz() == 0.0) |
---|
288 | return *this; |
---|
289 | |
---|
290 | double m_local = prest.m(); |
---|
291 | assert(m_local != 0); |
---|
292 | |
---|
293 | double pf4 = ( -px()*prest.px() - py()*prest.py() |
---|
294 | - pz()*prest.pz() + E()*prest.E() )/m_local; |
---|
295 | double fn = (pf4 + E()) / (prest.E() + m_local); |
---|
296 | _px -= fn*prest.px(); |
---|
297 | _py -= fn*prest.py(); |
---|
298 | _pz -= fn*prest.pz(); |
---|
299 | _E = pf4; |
---|
300 | |
---|
301 | _finish_init(); // we need to recalculate phi,rap,kt2 |
---|
302 | return *this; |
---|
303 | } |
---|
304 | |
---|
305 | |
---|
306 | //---------------------------------------------------------------------- |
---|
307 | /// returns true if the momenta of the two input jets are identical |
---|
308 | bool have_same_momentum(const PseudoJet & jeta, const PseudoJet & jetb) { |
---|
309 | return jeta.px() == jetb.px() |
---|
310 | && jeta.py() == jetb.py() |
---|
311 | && jeta.pz() == jetb.pz() |
---|
312 | && jeta.E() == jetb.E(); |
---|
313 | } |
---|
314 | |
---|
315 | //---------------------------------------------------------------------- |
---|
316 | void PseudoJet::set_cached_rap_phi(double rap_in, double phi_in) { |
---|
317 | _rap = rap_in; _phi = phi_in; |
---|
318 | if (_phi >= twopi) _phi -= twopi; |
---|
319 | if (_phi < 0) _phi += twopi; |
---|
320 | } |
---|
321 | |
---|
322 | //---------------------------------------------------------------------- |
---|
323 | void PseudoJet::reset_momentum_PtYPhiM(double pt_in, double y_in, double phi_in, double m_in) { |
---|
324 | assert(phi_in < 2*twopi && phi_in > -twopi); |
---|
325 | double ptm = (m_in == 0) ? pt_in : sqrt(pt_in*pt_in+m_in*m_in); |
---|
326 | double exprap = exp(y_in); |
---|
327 | double pminus = ptm/exprap; |
---|
328 | double pplus = ptm*exprap; |
---|
329 | double px_local = pt_in*cos(phi_in); |
---|
330 | double py_local = pt_in*sin(phi_in); |
---|
331 | reset_momentum(px_local,py_local,0.5*(pplus-pminus),0.5*(pplus+pminus)); |
---|
332 | set_cached_rap_phi(y_in,phi_in); |
---|
333 | } |
---|
334 | |
---|
335 | //---------------------------------------------------------------------- |
---|
336 | /// return a pseudojet with the given pt, y, phi and mass |
---|
337 | PseudoJet PtYPhiM(double pt, double y, double phi, double m) { |
---|
338 | assert(phi < 2*twopi && phi > -twopi); |
---|
339 | double ptm = (m == 0) ? pt : sqrt(pt*pt+m*m); |
---|
340 | double exprap = exp(y); |
---|
341 | double pminus = ptm/exprap; |
---|
342 | double pplus = ptm*exprap; |
---|
343 | double px = pt*cos(phi); |
---|
344 | double py = pt*sin(phi); |
---|
345 | PseudoJet mom(px,py,0.5*(pplus-pminus),0.5*(pplus+pminus)); |
---|
346 | mom.set_cached_rap_phi(y,phi); |
---|
347 | return mom; |
---|
348 | //return PseudoJet(pt*cos(phi), pt*sin(phi), ptm*sinh(y), ptm*cosh(y)); |
---|
349 | } |
---|
350 | |
---|
351 | |
---|
352 | //---------------------------------------------------------------------- |
---|
353 | // return kt-distance between this jet and another one |
---|
354 | double PseudoJet::kt_distance(const PseudoJet & other) const { |
---|
355 | //double distance = min(this->kt2(), other.kt2()); |
---|
356 | double distance = min(_kt2, other._kt2); |
---|
357 | double dphi = abs(phi() - other.phi()); |
---|
358 | if (dphi > pi) {dphi = twopi - dphi;} |
---|
359 | double drap = rap() - other.rap(); |
---|
360 | distance = distance * (dphi*dphi + drap*drap); |
---|
361 | return distance; |
---|
362 | } |
---|
363 | |
---|
364 | |
---|
365 | //---------------------------------------------------------------------- |
---|
366 | // return squared cylinder (eta-phi) distance between this jet and another one |
---|
367 | double PseudoJet::plain_distance(const PseudoJet & other) const { |
---|
368 | double dphi = abs(phi() - other.phi()); |
---|
369 | if (dphi > pi) {dphi = twopi - dphi;} |
---|
370 | double drap = rap() - other.rap(); |
---|
371 | return (dphi*dphi + drap*drap); |
---|
372 | } |
---|
373 | |
---|
374 | //---------------------------------------------------------------------- |
---|
375 | /// returns other.phi() - this.phi(), i.e. the phi distance to |
---|
376 | /// other, constrained to be in range -pi .. pi |
---|
377 | double PseudoJet::delta_phi_to(const PseudoJet & other) const { |
---|
378 | double dphi = other.phi() - phi(); |
---|
379 | if (dphi > pi) dphi -= twopi; |
---|
380 | if (dphi < -pi) dphi += twopi; |
---|
381 | return dphi; |
---|
382 | } |
---|
383 | |
---|
384 | |
---|
385 | string PseudoJet::description() const{ |
---|
386 | // the "default" case of a PJ which does not belong to any cluster sequence |
---|
387 | if (!_structure()) |
---|
388 | return "standard PseudoJet (with no associated clustering information)"; |
---|
389 | |
---|
390 | // for all the other cases, the description comes from the structure |
---|
391 | return _structure()->description(); |
---|
392 | } |
---|
393 | |
---|
394 | |
---|
395 | |
---|
396 | //---------------------------------------------------------------------- |
---|
397 | // |
---|
398 | // The following methods access the associated jet structure (if any) |
---|
399 | // |
---|
400 | //---------------------------------------------------------------------- |
---|
401 | |
---|
402 | |
---|
403 | //---------------------------------------------------------------------- |
---|
404 | // check whether this PseudoJet has an associated parent |
---|
405 | // ClusterSequence |
---|
406 | bool PseudoJet::has_associated_cluster_sequence() const{ |
---|
407 | return (_structure()) && (_structure->has_associated_cluster_sequence()); |
---|
408 | } |
---|
409 | |
---|
410 | //---------------------------------------------------------------------- |
---|
411 | // get a (const) pointer to the associated ClusterSequence (NULL if |
---|
412 | // inexistent) |
---|
413 | const ClusterSequence* PseudoJet::associated_cluster_sequence() const{ |
---|
414 | if (! has_associated_cluster_sequence()) return NULL; |
---|
415 | |
---|
416 | return _structure->associated_cluster_sequence(); |
---|
417 | } |
---|
418 | |
---|
419 | |
---|
420 | //---------------------------------------------------------------------- |
---|
421 | // check whether this PseudoJet has an associated parent |
---|
422 | // ClusterSequence that is still valid |
---|
423 | bool PseudoJet::has_valid_cluster_sequence() const{ |
---|
424 | return (_structure()) && (_structure->has_valid_cluster_sequence()); |
---|
425 | } |
---|
426 | |
---|
427 | //---------------------------------------------------------------------- |
---|
428 | // If there is a valid cluster sequence associated with this jet, |
---|
429 | // returns a pointer to it; otherwise throws an Error. |
---|
430 | // |
---|
431 | // Open question: should these errors be upgraded to classes of their |
---|
432 | // own so that they can be caught? [Maybe, but later] |
---|
433 | const ClusterSequence * PseudoJet::validated_cs() const { |
---|
434 | return validated_structure_ptr()->validated_cs(); |
---|
435 | } |
---|
436 | |
---|
437 | |
---|
438 | //---------------------------------------------------------------------- |
---|
439 | // set the associated structure |
---|
440 | void PseudoJet::set_structure_shared_ptr(const SharedPtr<PseudoJetStructureBase> &structure_in){ |
---|
441 | _structure = structure_in; |
---|
442 | } |
---|
443 | |
---|
444 | //---------------------------------------------------------------------- |
---|
445 | // return true if there is some strusture associated with this PseudoJet |
---|
446 | bool PseudoJet::has_structure() const{ |
---|
447 | return _structure(); |
---|
448 | } |
---|
449 | |
---|
450 | //---------------------------------------------------------------------- |
---|
451 | // return a pointer to the structure (of type |
---|
452 | // PseudoJetStructureBase*) associated with this PseudoJet. |
---|
453 | // |
---|
454 | // return NULL if there is no associated structure |
---|
455 | const PseudoJetStructureBase* PseudoJet::structure_ptr() const { |
---|
456 | if (!_structure()) return NULL; |
---|
457 | return _structure(); |
---|
458 | } |
---|
459 | |
---|
460 | //---------------------------------------------------------------------- |
---|
461 | // return a non-const pointer to the structure (of type |
---|
462 | // PseudoJetStructureBase*) associated with this PseudoJet. |
---|
463 | // |
---|
464 | // return NULL if there is no associated structure |
---|
465 | // |
---|
466 | // Only use this if you know what you are doing. In any case, |
---|
467 | // prefer the 'structure_ptr()' (the const version) to this method, |
---|
468 | // unless you really need a write access to the PseudoJet's |
---|
469 | // underlying structure. |
---|
470 | PseudoJetStructureBase* PseudoJet::structure_non_const_ptr(){ |
---|
471 | if (!_structure()) return NULL; |
---|
472 | return _structure(); |
---|
473 | } |
---|
474 | |
---|
475 | //---------------------------------------------------------------------- |
---|
476 | // return a pointer to the structure (of type |
---|
477 | // PseudoJetStructureBase*) associated with this PseudoJet. |
---|
478 | // |
---|
479 | // throw an error if there is no associated structure |
---|
480 | const PseudoJetStructureBase* PseudoJet::validated_structure_ptr() const { |
---|
481 | if (!_structure()) |
---|
482 | throw Error("Trying to access the structure of a PseudoJet which has no associated structure"); |
---|
483 | return _structure(); |
---|
484 | } |
---|
485 | |
---|
486 | //---------------------------------------------------------------------- |
---|
487 | // return a reference to the shared pointer to the |
---|
488 | // PseudoJetStructureBase associated with this PseudoJet |
---|
489 | const SharedPtr<PseudoJetStructureBase> & PseudoJet::structure_shared_ptr() const { |
---|
490 | return _structure; |
---|
491 | } |
---|
492 | |
---|
493 | |
---|
494 | //---------------------------------------------------------------------- |
---|
495 | // check if it has been recombined with another PseudoJet in which |
---|
496 | // case, return its partner through the argument. Otherwise, |
---|
497 | // 'partner' is set to 0. |
---|
498 | // |
---|
499 | // false is also returned if this PseudoJet has no associated |
---|
500 | // ClusterSequence |
---|
501 | bool PseudoJet::has_partner(PseudoJet &partner) const{ |
---|
502 | return validated_structure_ptr()->has_partner(*this, partner); |
---|
503 | } |
---|
504 | |
---|
505 | //---------------------------------------------------------------------- |
---|
506 | // check if it has been recombined with another PseudoJet in which |
---|
507 | // case, return its child through the argument. Otherwise, 'child' |
---|
508 | // is set to 0. |
---|
509 | // |
---|
510 | // false is also returned if this PseudoJet has no associated |
---|
511 | // ClusterSequence, with the child set to 0 |
---|
512 | bool PseudoJet::has_child(PseudoJet &child) const{ |
---|
513 | return validated_structure_ptr()->has_child(*this, child); |
---|
514 | } |
---|
515 | |
---|
516 | //---------------------------------------------------------------------- |
---|
517 | // check if it is the product of a recombination, in which case |
---|
518 | // return the 2 parents through the 'parent1' and 'parent2' |
---|
519 | // arguments. Otherwise, set these to 0. |
---|
520 | // |
---|
521 | // false is also returned if this PseudoJet has no parent |
---|
522 | // ClusterSequence |
---|
523 | bool PseudoJet::has_parents(PseudoJet &parent1, PseudoJet &parent2) const{ |
---|
524 | return validated_structure_ptr()->has_parents(*this, parent1, parent2); |
---|
525 | } |
---|
526 | |
---|
527 | //---------------------------------------------------------------------- |
---|
528 | // check if the current PseudoJet contains the one passed as |
---|
529 | // argument |
---|
530 | // |
---|
531 | // false is also returned if this PseudoJet has no associated |
---|
532 | // ClusterSequence. |
---|
533 | bool PseudoJet::contains(const PseudoJet &constituent) const{ |
---|
534 | return validated_structure_ptr()->object_in_jet(constituent, *this); |
---|
535 | } |
---|
536 | |
---|
537 | //---------------------------------------------------------------------- |
---|
538 | // check if the current PseudoJet is contained the one passed as |
---|
539 | // argument |
---|
540 | // |
---|
541 | // false is also returned if this PseudoJet has no associated |
---|
542 | // ClusterSequence |
---|
543 | bool PseudoJet::is_inside(const PseudoJet &jet) const{ |
---|
544 | return validated_structure_ptr()->object_in_jet(*this, jet); |
---|
545 | } |
---|
546 | |
---|
547 | |
---|
548 | //---------------------------------------------------------------------- |
---|
549 | // returns true if the PseudoJet has constituents |
---|
550 | bool PseudoJet::has_constituents() const{ |
---|
551 | return (_structure()) && (_structure->has_constituents()); |
---|
552 | } |
---|
553 | |
---|
554 | //---------------------------------------------------------------------- |
---|
555 | // retrieve the constituents. |
---|
556 | vector<PseudoJet> PseudoJet::constituents() const{ |
---|
557 | return validated_structure_ptr()->constituents(*this); |
---|
558 | } |
---|
559 | |
---|
560 | |
---|
561 | //---------------------------------------------------------------------- |
---|
562 | // returns true if the PseudoJet has support for exclusive subjets |
---|
563 | bool PseudoJet::has_exclusive_subjets() const{ |
---|
564 | return (_structure()) && (_structure->has_exclusive_subjets()); |
---|
565 | } |
---|
566 | |
---|
567 | //---------------------------------------------------------------------- |
---|
568 | // return a vector of all subjets of the current jet (in the sense |
---|
569 | // of the exclusive algorithm) that would be obtained when running |
---|
570 | // the algorithm with the given dcut. |
---|
571 | // |
---|
572 | // Time taken is O(m ln m), where m is the number of subjets that |
---|
573 | // are found. If m gets to be of order of the total number of |
---|
574 | // constituents in the jet, this could be substantially slower than |
---|
575 | // just getting that list of constituents. |
---|
576 | // |
---|
577 | // an Error is thrown if this PseudoJet has no currently valid |
---|
578 | // associated ClusterSequence |
---|
579 | std::vector<PseudoJet> PseudoJet::exclusive_subjets (const double & dcut) const { |
---|
580 | return validated_structure_ptr()->exclusive_subjets(*this, dcut); |
---|
581 | } |
---|
582 | |
---|
583 | //---------------------------------------------------------------------- |
---|
584 | // return the size of exclusive_subjets(...); still n ln n with same |
---|
585 | // coefficient, but marginally more efficient than manually taking |
---|
586 | // exclusive_subjets.size() |
---|
587 | // |
---|
588 | // an Error is thrown if this PseudoJet has no currently valid |
---|
589 | // associated ClusterSequence |
---|
590 | int PseudoJet::n_exclusive_subjets(const double & dcut) const { |
---|
591 | return validated_structure_ptr()->n_exclusive_subjets(*this, dcut); |
---|
592 | } |
---|
593 | |
---|
594 | //---------------------------------------------------------------------- |
---|
595 | // return the list of subjets obtained by unclustering the supplied |
---|
596 | // jet down to n subjets (or all constituents if there are fewer |
---|
597 | // than n). |
---|
598 | // |
---|
599 | // requires n ln n time |
---|
600 | // |
---|
601 | // an Error is thrown if this PseudoJet has no currently valid |
---|
602 | // associated ClusterSequence |
---|
603 | std::vector<PseudoJet> PseudoJet::exclusive_subjets_up_to (int nsub) const { |
---|
604 | return validated_structure_ptr()->exclusive_subjets_up_to(*this, nsub); |
---|
605 | } |
---|
606 | |
---|
607 | //---------------------------------------------------------------------- |
---|
608 | // Same as exclusive_subjets_up_to but throws an error if there are |
---|
609 | // fewer than nsub particles in the jet |
---|
610 | std::vector<PseudoJet> PseudoJet::exclusive_subjets (int nsub) const { |
---|
611 | vector<PseudoJet> subjets = exclusive_subjets_up_to(nsub); |
---|
612 | if (int(subjets.size()) < nsub) { |
---|
613 | ostringstream err; |
---|
614 | err << "Requested " << nsub << " exclusive subjets, but there were only " |
---|
615 | << subjets.size() << " particles in the jet"; |
---|
616 | throw Error(err.str()); |
---|
617 | } |
---|
618 | return subjets; |
---|
619 | } |
---|
620 | |
---|
621 | //---------------------------------------------------------------------- |
---|
622 | // return the dij that was present in the merging nsub+1 -> nsub |
---|
623 | // subjets inside this jet. |
---|
624 | // |
---|
625 | // an Error is thrown if this PseudoJet has no currently valid |
---|
626 | // associated ClusterSequence |
---|
627 | double PseudoJet::exclusive_subdmerge(int nsub) const { |
---|
628 | return validated_structure_ptr()->exclusive_subdmerge(*this, nsub); |
---|
629 | } |
---|
630 | |
---|
631 | //---------------------------------------------------------------------- |
---|
632 | // return the maximum dij that occurred in the whole event at the |
---|
633 | // stage that the nsub+1 -> nsub merge of subjets occurred inside |
---|
634 | // this jet. |
---|
635 | // |
---|
636 | // an Error is thrown if this PseudoJet has no currently valid |
---|
637 | // associated ClusterSequence |
---|
638 | double PseudoJet::exclusive_subdmerge_max(int nsub) const { |
---|
639 | return validated_structure_ptr()->exclusive_subdmerge_max(*this, nsub); |
---|
640 | } |
---|
641 | |
---|
642 | |
---|
643 | // returns true if a jet has pieces |
---|
644 | // |
---|
645 | // By default a single particle or a jet coming from a |
---|
646 | // ClusterSequence have no pieces and this methos will return false. |
---|
647 | bool PseudoJet::has_pieces() const{ |
---|
648 | return ((_structure()) && (_structure->has_pieces(*this))); |
---|
649 | } |
---|
650 | |
---|
651 | // retrieve the pieces that make up the jet. |
---|
652 | // |
---|
653 | // By default a jet does not have pieces. |
---|
654 | // If the underlying interface supports "pieces" retrieve the |
---|
655 | // pieces from there. |
---|
656 | std::vector<PseudoJet> PseudoJet::pieces() const{ |
---|
657 | return validated_structure_ptr()->pieces(*this); |
---|
658 | // if (!has_pieces()) |
---|
659 | // throw Error("Trying to retrieve the pieces of a PseudoJet that has no support for pieces."); |
---|
660 | // |
---|
661 | // return _structure->pieces(*this); |
---|
662 | } |
---|
663 | |
---|
664 | |
---|
665 | //---------------------------------------------------------------------- |
---|
666 | // the following ones require a computation of the area in the |
---|
667 | // associated ClusterSequence (See ClusterSequenceAreaBase for details) |
---|
668 | //---------------------------------------------------------------------- |
---|
669 | |
---|
670 | //---------------------------------------------------------------------- |
---|
671 | // if possible, return a valid ClusterSequenceAreaBase pointer; otherwise |
---|
672 | // throw an error |
---|
673 | const ClusterSequenceAreaBase * PseudoJet::validated_csab() const { |
---|
674 | const ClusterSequenceAreaBase *csab = dynamic_cast<const ClusterSequenceAreaBase*>(validated_cs()); |
---|
675 | if (csab == NULL) throw Error("you requested jet-area related information, but the PseudoJet does not have associated area information."); |
---|
676 | return csab; |
---|
677 | } |
---|
678 | |
---|
679 | |
---|
680 | //---------------------------------------------------------------------- |
---|
681 | // check if it has a defined area |
---|
682 | bool PseudoJet::has_area() const{ |
---|
683 | //if (! has_associated_cluster_sequence()) return false; |
---|
684 | if (! has_structure()) return false; |
---|
685 | return (validated_structure_ptr()->has_area() != 0); |
---|
686 | } |
---|
687 | |
---|
688 | //---------------------------------------------------------------------- |
---|
689 | // return the jet (scalar) area. |
---|
690 | // throw an Error if there is no support for area in the associated CS |
---|
691 | double PseudoJet::area() const{ |
---|
692 | return validated_structure_ptr()->area(*this); |
---|
693 | } |
---|
694 | |
---|
695 | //---------------------------------------------------------------------- |
---|
696 | // return the error (uncertainty) associated with the determination |
---|
697 | // of the area of this jet. |
---|
698 | // throws an Error if there is no support for area in the associated CS |
---|
699 | double PseudoJet::area_error() const{ |
---|
700 | return validated_structure_ptr()->area_error(*this); |
---|
701 | } |
---|
702 | |
---|
703 | //---------------------------------------------------------------------- |
---|
704 | // return the jet 4-vector area |
---|
705 | // throws an Error if there is no support for area in the associated CS |
---|
706 | PseudoJet PseudoJet::area_4vector() const{ |
---|
707 | return validated_structure_ptr()->area_4vector(*this); |
---|
708 | } |
---|
709 | |
---|
710 | //---------------------------------------------------------------------- |
---|
711 | // true if this jet is made exclusively of ghosts |
---|
712 | // throws an Error if there is no support for area in the associated CS |
---|
713 | bool PseudoJet::is_pure_ghost() const{ |
---|
714 | return validated_structure_ptr()->is_pure_ghost(*this); |
---|
715 | } |
---|
716 | |
---|
717 | |
---|
718 | //---------------------------------------------------------------------- |
---|
719 | // |
---|
720 | // end of the methods accessing the information in the associated |
---|
721 | // Cluster Sequence |
---|
722 | // |
---|
723 | //---------------------------------------------------------------------- |
---|
724 | |
---|
725 | //---------------------------------------------------------------------- |
---|
726 | /// provide a meaningful error message for InexistentUserInfo |
---|
727 | PseudoJet::InexistentUserInfo::InexistentUserInfo() : Error("you attempted to perform a dynamic cast of a PseudoJet's extra info, but the extra info pointer was null") |
---|
728 | {} |
---|
729 | |
---|
730 | |
---|
731 | //---------------------------------------------------------------------- |
---|
732 | // sort the indices so that values[indices[0..n-1]] is sorted |
---|
733 | // into increasing order |
---|
734 | void sort_indices(vector<int> & indices, |
---|
735 | const vector<double> & values) { |
---|
736 | IndexedSortHelper index_sort_helper(&values); |
---|
737 | sort(indices.begin(), indices.end(), index_sort_helper); |
---|
738 | } |
---|
739 | |
---|
740 | |
---|
741 | |
---|
742 | //---------------------------------------------------------------------- |
---|
743 | /// given a vector of values with a one-to-one correspondence with the |
---|
744 | /// vector of objects, sort objects into an order such that the |
---|
745 | /// associated values would be in increasing order |
---|
746 | template<class T> vector<T> objects_sorted_by_values( |
---|
747 | const vector<T> & objects, |
---|
748 | const vector<double> & values) { |
---|
749 | |
---|
750 | assert(objects.size() == values.size()); |
---|
751 | |
---|
752 | // get a vector of indices |
---|
753 | vector<int> indices(values.size()); |
---|
754 | for (size_t i = 0; i < indices.size(); i++) {indices[i] = i;} |
---|
755 | |
---|
756 | // sort the indices |
---|
757 | sort_indices(indices, values); |
---|
758 | |
---|
759 | // copy the objects |
---|
760 | vector<T> objects_sorted(objects.size()); |
---|
761 | |
---|
762 | // place the objects in the correct order |
---|
763 | for (size_t i = 0; i < indices.size(); i++) { |
---|
764 | objects_sorted[i] = objects[indices[i]]; |
---|
765 | } |
---|
766 | |
---|
767 | return objects_sorted; |
---|
768 | } |
---|
769 | |
---|
770 | //---------------------------------------------------------------------- |
---|
771 | /// return a vector of jets sorted into decreasing kt2 |
---|
772 | vector<PseudoJet> sorted_by_pt(const vector<PseudoJet> & jets) { |
---|
773 | vector<double> minus_kt2(jets.size()); |
---|
774 | for (size_t i = 0; i < jets.size(); i++) {minus_kt2[i] = -jets[i].kt2();} |
---|
775 | return objects_sorted_by_values(jets, minus_kt2); |
---|
776 | } |
---|
777 | |
---|
778 | //---------------------------------------------------------------------- |
---|
779 | /// return a vector of jets sorted into increasing rapidity |
---|
780 | vector<PseudoJet> sorted_by_rapidity(const vector<PseudoJet> & jets) { |
---|
781 | vector<double> rapidities(jets.size()); |
---|
782 | for (size_t i = 0; i < jets.size(); i++) {rapidities[i] = jets[i].rap();} |
---|
783 | return objects_sorted_by_values(jets, rapidities); |
---|
784 | } |
---|
785 | |
---|
786 | //---------------------------------------------------------------------- |
---|
787 | /// return a vector of jets sorted into decreasing energy |
---|
788 | vector<PseudoJet> sorted_by_E(const vector<PseudoJet> & jets) { |
---|
789 | vector<double> energies(jets.size()); |
---|
790 | for (size_t i = 0; i < jets.size(); i++) {energies[i] = -jets[i].E();} |
---|
791 | return objects_sorted_by_values(jets, energies); |
---|
792 | } |
---|
793 | |
---|
794 | //---------------------------------------------------------------------- |
---|
795 | /// return a vector of jets sorted into increasing pz |
---|
796 | vector<PseudoJet> sorted_by_pz(const vector<PseudoJet> & jets) { |
---|
797 | vector<double> pz(jets.size()); |
---|
798 | for (size_t i = 0; i < jets.size(); i++) {pz[i] = jets[i].pz();} |
---|
799 | return objects_sorted_by_values(jets, pz); |
---|
800 | } |
---|
801 | |
---|
802 | |
---|
803 | |
---|
804 | //------------------------------------------------------------------------------- |
---|
805 | // helper functions to build a jet made of pieces |
---|
806 | //------------------------------------------------------------------------------- |
---|
807 | |
---|
808 | // build a "CompositeJet" from the vector of its pieces |
---|
809 | // |
---|
810 | // In this case, E-scheme recombination is assumed to compute the |
---|
811 | // total momentum |
---|
812 | PseudoJet join(const vector<PseudoJet> & pieces){ |
---|
813 | // compute the total momentum |
---|
814 | //-------------------------------------------------- |
---|
815 | PseudoJet result; // automatically initialised to 0 |
---|
816 | for (unsigned int i=0; i<pieces.size(); i++) |
---|
817 | result += pieces[i]; |
---|
818 | |
---|
819 | // attach a CompositeJetStructure to the result |
---|
820 | //-------------------------------------------------- |
---|
821 | CompositeJetStructure *cj_struct = new CompositeJetStructure(pieces); |
---|
822 | |
---|
823 | result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct)); |
---|
824 | |
---|
825 | return result; |
---|
826 | } |
---|
827 | |
---|
828 | // build a "CompositeJet" from a single PseudoJet |
---|
829 | PseudoJet join(const PseudoJet & j1){ |
---|
830 | return join(vector<PseudoJet>(1,j1)); |
---|
831 | } |
---|
832 | |
---|
833 | // build a "CompositeJet" from two PseudoJet |
---|
834 | PseudoJet join(const PseudoJet & j1, const PseudoJet & j2){ |
---|
835 | vector<PseudoJet> pieces; |
---|
836 | pieces.push_back(j1); |
---|
837 | pieces.push_back(j2); |
---|
838 | return join(pieces); |
---|
839 | } |
---|
840 | |
---|
841 | // build a "CompositeJet" from 3 PseudoJet |
---|
842 | PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3){ |
---|
843 | vector<PseudoJet> pieces; |
---|
844 | pieces.push_back(j1); |
---|
845 | pieces.push_back(j2); |
---|
846 | pieces.push_back(j3); |
---|
847 | return join(pieces); |
---|
848 | } |
---|
849 | |
---|
850 | // build a "CompositeJet" from 4 PseudoJet |
---|
851 | PseudoJet join(const PseudoJet & j1, const PseudoJet & j2, const PseudoJet & j3, const PseudoJet & j4){ |
---|
852 | vector<PseudoJet> pieces; |
---|
853 | pieces.push_back(j1); |
---|
854 | pieces.push_back(j2); |
---|
855 | pieces.push_back(j3); |
---|
856 | pieces.push_back(j4); |
---|
857 | return join(pieces); |
---|
858 | } |
---|
859 | |
---|
860 | |
---|
861 | |
---|
862 | |
---|
863 | FASTJET_END_NAMESPACE |
---|
864 | |
---|