1 | //STARTHEADER |
---|
2 | // $Id: JetDefinition.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 | #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 | //---------------------------------------------------------------------- |
---|
108 | string JetDefinition::description() const { |
---|
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(); |
---|
127 | } else if (jet_algorithm() == cambridge_for_passive_algorithm) { |
---|
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 | |
---|
147 | void JetDefinition::set_recombination_scheme( |
---|
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 |
---|
160 | bool JetDefinition::has_same_recombiner(const JetDefinition &other_jd) const{ |
---|
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 |
---|
173 | void JetDefinition::delete_recombiner_when_unused(){ |
---|
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 |
---|
183 | void JetDefinition::delete_plugin_when_unused(){ |
---|
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 | |
---|
193 | string JetDefinition::DefaultRecombiner::description() const { |
---|
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 | |
---|
218 | void JetDefinition::DefaultRecombiner::recombine( |
---|
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 | |
---|
278 | void JetDefinition::DefaultRecombiner::preprocess(PseudoJet & p) const { |
---|
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 | |
---|
318 | void JetDefinition::Plugin::set_ghost_separation_scale(double /*scale*/) const { |
---|
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 | |
---|
348 | result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(cj_struct)); |
---|
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 |
---|