1 | //STARTHEADER |
---|
2 | // $Id: ClusterSequenceActiveArea.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/PseudoJet.hh" |
---|
30 | #include "fastjet/ClusterSequence.hh" |
---|
31 | #include "fastjet/ClusterSequenceActiveArea.hh" |
---|
32 | #include "fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh" |
---|
33 | #include<iostream> |
---|
34 | #include<vector> |
---|
35 | #include<sstream> |
---|
36 | #include<algorithm> |
---|
37 | #include<cmath> |
---|
38 | #include<valarray> |
---|
39 | |
---|
40 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
41 | |
---|
42 | |
---|
43 | using namespace std; |
---|
44 | |
---|
45 | |
---|
46 | //int ClusterSequenceActiveArea::_n_seed_warnings = 0; |
---|
47 | //const int _max_seed_warnings = 10; |
---|
48 | |
---|
49 | //---------------------------------------------------------------------- |
---|
50 | /// global routine for running active area |
---|
51 | void ClusterSequenceActiveArea::_initialise_and_run_AA ( |
---|
52 | const JetDefinition & jet_def_in, |
---|
53 | const GhostedAreaSpec & ghost_spec, |
---|
54 | const bool & writeout_combinations) { |
---|
55 | |
---|
56 | bool continue_running; |
---|
57 | _initialise_AA(jet_def_in, ghost_spec, writeout_combinations, continue_running); |
---|
58 | if (continue_running) { |
---|
59 | _run_AA(ghost_spec); |
---|
60 | _postprocess_AA(ghost_spec); |
---|
61 | } |
---|
62 | } |
---|
63 | |
---|
64 | //---------------------------------------------------------------------- |
---|
65 | void ClusterSequenceActiveArea::_resize_and_zero_AA () { |
---|
66 | // initialize our local area information |
---|
67 | _average_area.resize(2*_jets.size()); _average_area = 0.0; |
---|
68 | _average_area2.resize(2*_jets.size()); _average_area2 = 0.0; |
---|
69 | _average_area_4vector.resize(2*_jets.size()); |
---|
70 | _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0); |
---|
71 | _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0; |
---|
72 | } |
---|
73 | |
---|
74 | //---------------------------------a------------------------------------- |
---|
75 | void ClusterSequenceActiveArea::_initialise_AA ( |
---|
76 | const JetDefinition & jet_def_in, |
---|
77 | const GhostedAreaSpec & ghost_spec, |
---|
78 | const bool & writeout_combinations, |
---|
79 | bool & continue_running) |
---|
80 | { |
---|
81 | |
---|
82 | // store this for future use |
---|
83 | _ghost_spec_repeat = ghost_spec.repeat(); |
---|
84 | |
---|
85 | // make sure placeholders are there & zeroed |
---|
86 | _resize_and_zero_AA(); |
---|
87 | |
---|
88 | // for future reference... |
---|
89 | _maxrap_for_area = ghost_spec.ghost_maxrap(); |
---|
90 | _safe_rap_for_area = _maxrap_for_area - jet_def_in.R(); |
---|
91 | |
---|
92 | // Make sure we'll have at least one repetition -- then we can |
---|
93 | // deduce the unghosted clustering sequence from one of the ghosted |
---|
94 | // sequences. If we do not have any repetitions, then get the |
---|
95 | // unghosted sequence from the plain unghosted clustering. |
---|
96 | // |
---|
97 | // NB: all decanting and filling of initial history will then |
---|
98 | // be carried out by base-class routine |
---|
99 | if (ghost_spec.repeat() <= 0) { |
---|
100 | _initialise_and_run(jet_def_in, writeout_combinations); |
---|
101 | continue_running = false; |
---|
102 | return; |
---|
103 | } |
---|
104 | |
---|
105 | // transfer all relevant info into internal variables |
---|
106 | _decant_options(jet_def_in, writeout_combinations); |
---|
107 | |
---|
108 | // set up the history entries for the initial particles (those |
---|
109 | // currently in _jets) |
---|
110 | _fill_initial_history(); |
---|
111 | |
---|
112 | // by default it does not... |
---|
113 | _has_dangerous_particles = false; |
---|
114 | |
---|
115 | continue_running = true; |
---|
116 | } |
---|
117 | |
---|
118 | |
---|
119 | //---------------------------------------------------------------------- |
---|
120 | void ClusterSequenceActiveArea::_run_AA (const GhostedAreaSpec & ghost_spec) { |
---|
121 | // record the input jets as they are currently |
---|
122 | vector<PseudoJet> input_jets(_jets); |
---|
123 | |
---|
124 | // code for testing the unique tree |
---|
125 | vector<int> unique_tree; |
---|
126 | |
---|
127 | // run the clustering multiple times so as to get areas of all the jets |
---|
128 | for (int irepeat = 0; irepeat < ghost_spec.repeat(); irepeat++) { |
---|
129 | |
---|
130 | ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets, |
---|
131 | jet_def(), ghost_spec); |
---|
132 | |
---|
133 | _has_dangerous_particles |= clust_seq.has_dangerous_particles(); |
---|
134 | if (irepeat == 0) { |
---|
135 | // take the non-ghost part of the history and put into our own |
---|
136 | // history. |
---|
137 | _transfer_ghost_free_history(clust_seq); |
---|
138 | // get the "unique" order that will be used for transferring all areas. |
---|
139 | unique_tree = unique_history_order(); |
---|
140 | } |
---|
141 | |
---|
142 | // transfer areas from clust_seq into our object |
---|
143 | _transfer_areas(unique_tree, clust_seq); |
---|
144 | } |
---|
145 | } |
---|
146 | |
---|
147 | |
---|
148 | //---------------------------------------------------------------------- |
---|
149 | /// run the postprocessing for the active area (and derived classes) |
---|
150 | void ClusterSequenceActiveArea::_postprocess_AA (const GhostedAreaSpec & ghost_spec) { |
---|
151 | _average_area /= ghost_spec.repeat(); |
---|
152 | _average_area2 /= ghost_spec.repeat(); |
---|
153 | if (ghost_spec.repeat() > 1) { |
---|
154 | // the VC compiler complains if one puts everything on a single line. |
---|
155 | // An alternative solution would be to use -1.0 (+single line) |
---|
156 | const double tmp = ghost_spec.repeat()-1; |
---|
157 | _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/tmp); |
---|
158 | } else { |
---|
159 | _average_area2 = 0.0; |
---|
160 | } |
---|
161 | |
---|
162 | _non_jet_area /= ghost_spec.repeat(); |
---|
163 | _non_jet_area2 /= ghost_spec.repeat(); |
---|
164 | _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/ |
---|
165 | ghost_spec.repeat()); |
---|
166 | _non_jet_number /= ghost_spec.repeat(); |
---|
167 | |
---|
168 | // following bizarre way of writing things is related to |
---|
169 | // poverty of operations on PseudoJet objects (as well as some confusion |
---|
170 | // in one or two places) |
---|
171 | for (unsigned i = 0; i < _average_area_4vector.size(); i++) { |
---|
172 | _average_area_4vector[i] = (1.0/ghost_spec.repeat()) * _average_area_4vector[i]; |
---|
173 | } |
---|
174 | //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl; |
---|
175 | } |
---|
176 | |
---|
177 | |
---|
178 | // //---------------------------------------------------------------------- |
---|
179 | // void ClusterSequenceActiveArea::_initialise_and_run_AA ( |
---|
180 | // const JetDefinition & jet_def, |
---|
181 | // const GhostedAreaSpec & ghost_spec, |
---|
182 | // const bool & writeout_combinations) |
---|
183 | // { |
---|
184 | // |
---|
185 | // // store this for future use |
---|
186 | // _ghost_spec_repeat = ghost_spec.repeat(); |
---|
187 | // |
---|
188 | // // initialize our local area information |
---|
189 | // _average_area.resize(2*_jets.size()); _average_area = 0.0; |
---|
190 | // _average_area2.resize(2*_jets.size()); _average_area2 = 0.0; |
---|
191 | // _average_area_4vector.resize(2*_jets.size()); |
---|
192 | // _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0); |
---|
193 | // _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0; |
---|
194 | // |
---|
195 | // // for future reference... |
---|
196 | // _maxrap_for_area = ghost_spec.ghost_maxrap(); |
---|
197 | // _safe_rap_for_area = _maxrap_for_area - jet_def.R(); |
---|
198 | // |
---|
199 | // // Make sure we'll have at least one repetition -- then we can |
---|
200 | // // deduce the unghosted clustering sequence from one of the ghosted |
---|
201 | // // sequences. If we do not have any repetitions, then get the |
---|
202 | // // unghosted sequence from the plain unghosted clustering. |
---|
203 | // // |
---|
204 | // // NB: all decanting and filling of initial history will then |
---|
205 | // // be carried out by base-class routine |
---|
206 | // if (ghost_spec.repeat() <= 0) { |
---|
207 | // _initialise_and_run(jet_def, writeout_combinations); |
---|
208 | // return; |
---|
209 | // } |
---|
210 | // |
---|
211 | // // transfer all relevant info into internal variables |
---|
212 | // _decant_options(jet_def, writeout_combinations); |
---|
213 | // |
---|
214 | // // set up the history entries for the initial particles (those |
---|
215 | // // currently in _jets) |
---|
216 | // _fill_initial_history(); |
---|
217 | // |
---|
218 | // // record the input jets as they are currently |
---|
219 | // vector<PseudoJet> input_jets(_jets); |
---|
220 | // |
---|
221 | // // code for testing the unique tree |
---|
222 | // vector<int> unique_tree; |
---|
223 | // |
---|
224 | // |
---|
225 | // |
---|
226 | // |
---|
227 | // // run the clustering multiple times so as to get areas of all the jets |
---|
228 | // for (int irepeat = 0; irepeat < ghost_spec.repeat(); irepeat++) { |
---|
229 | // |
---|
230 | // ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets, |
---|
231 | // jet_def, ghost_spec); |
---|
232 | // |
---|
233 | // if (irepeat == 0) { |
---|
234 | // // take the non-ghost part of the history and put into our own |
---|
235 | // // history. |
---|
236 | // _transfer_ghost_free_history(clust_seq); |
---|
237 | // // get the "unique" order that will be used for transferring all areas. |
---|
238 | // unique_tree = unique_history_order(); |
---|
239 | // } |
---|
240 | // |
---|
241 | // // transfer areas from clust_seq into our object |
---|
242 | // _transfer_areas(unique_tree, clust_seq); |
---|
243 | // } |
---|
244 | // |
---|
245 | // _average_area /= ghost_spec.repeat(); |
---|
246 | // _average_area2 /= ghost_spec.repeat(); |
---|
247 | // if (ghost_spec.repeat() > 1) { |
---|
248 | // _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/ |
---|
249 | // (ghost_spec.repeat()-1)); |
---|
250 | // } else { |
---|
251 | // _average_area2 = 0.0; |
---|
252 | // } |
---|
253 | // |
---|
254 | // _non_jet_area /= ghost_spec.repeat(); |
---|
255 | // _non_jet_area2 /= ghost_spec.repeat(); |
---|
256 | // _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/ |
---|
257 | // ghost_spec.repeat()); |
---|
258 | // _non_jet_number /= ghost_spec.repeat(); |
---|
259 | // |
---|
260 | // // following bizarre way of writing things is related to |
---|
261 | // // poverty of operations on PseudoJet objects (as well as some confusion |
---|
262 | // // in one or two places) |
---|
263 | // for (unsigned i = 0; i < _average_area_4vector.size(); i++) { |
---|
264 | // _average_area_4vector[i] = (1.0/ghost_spec.repeat()) * _average_area_4vector[i]; |
---|
265 | // } |
---|
266 | // //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl; |
---|
267 | // |
---|
268 | // |
---|
269 | // } |
---|
270 | // |
---|
271 | |
---|
272 | |
---|
273 | //---------------------------------------------------------------------- |
---|
274 | double ClusterSequenceActiveArea::pt_per_unit_area( |
---|
275 | mean_pt_strategies strat, double range) const { |
---|
276 | |
---|
277 | vector<PseudoJet> incl_jets = inclusive_jets(); |
---|
278 | vector<double> pt_over_areas; |
---|
279 | |
---|
280 | for (unsigned i = 0; i < incl_jets.size(); i++) { |
---|
281 | if (abs(incl_jets[i].rap()) < _safe_rap_for_area) { |
---|
282 | double this_area; |
---|
283 | if ( strat == median_4vector ) { |
---|
284 | this_area = area_4vector(incl_jets[i]).perp(); |
---|
285 | } else { |
---|
286 | this_area = area(incl_jets[i]); |
---|
287 | } |
---|
288 | pt_over_areas.push_back(incl_jets[i].perp()/this_area); |
---|
289 | } |
---|
290 | } |
---|
291 | |
---|
292 | // there is nothing inside our region, so answer will always be zero |
---|
293 | if (pt_over_areas.size() == 0) {return 0.0;} |
---|
294 | |
---|
295 | // get median (pt/area) [this is the "old" median definition. It considers |
---|
296 | // only the "real" jets in calculating the median, i.e. excluding the |
---|
297 | // only-ghost ones] |
---|
298 | sort(pt_over_areas.begin(), pt_over_areas.end()); |
---|
299 | double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2]; |
---|
300 | |
---|
301 | // new median definition that takes into account non-jet area (i.e. |
---|
302 | // jets composed only of ghosts), and for fractional median position |
---|
303 | // interpolates between the corresponding entries in the pt_over_areas array |
---|
304 | double nj_median_pos = (pt_over_areas.size()-1 - _non_jet_number)/2.0; |
---|
305 | double nj_median_ratio; |
---|
306 | if (nj_median_pos >= 0 && pt_over_areas.size() > 1) { |
---|
307 | int int_nj_median = int(nj_median_pos); |
---|
308 | nj_median_ratio = |
---|
309 | pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos) |
---|
310 | + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median); |
---|
311 | } else { |
---|
312 | nj_median_ratio = 0.0; |
---|
313 | } |
---|
314 | |
---|
315 | |
---|
316 | // get various forms of mean (pt/area) |
---|
317 | double pt_sum = 0.0, pt_sum_with_cut = 0.0; |
---|
318 | double area_sum = _non_jet_area, area_sum_with_cut = _non_jet_area; |
---|
319 | double ratio_sum = 0.0; |
---|
320 | double ratio_n = _non_jet_number; |
---|
321 | for (unsigned i = 0; i < incl_jets.size(); i++) { |
---|
322 | if (abs(incl_jets[i].rap()) < _safe_rap_for_area) { |
---|
323 | double this_area; |
---|
324 | if ( strat == median_4vector ) { |
---|
325 | this_area = area_4vector(incl_jets[i]).perp(); |
---|
326 | } else { |
---|
327 | this_area = area(incl_jets[i]); |
---|
328 | } |
---|
329 | pt_sum += incl_jets[i].perp(); |
---|
330 | area_sum += this_area; |
---|
331 | double ratio = incl_jets[i].perp()/this_area; |
---|
332 | if (ratio < range*nj_median_ratio) { |
---|
333 | pt_sum_with_cut += incl_jets[i].perp(); |
---|
334 | area_sum_with_cut += this_area; |
---|
335 | ratio_sum += ratio; ratio_n++; |
---|
336 | } |
---|
337 | } |
---|
338 | } |
---|
339 | |
---|
340 | if (strat == play) { |
---|
341 | double trunc_sum = 0, trunc_sumsqr = 0; |
---|
342 | vector<double> means(pt_over_areas.size()), sd(pt_over_areas.size()); |
---|
343 | for (unsigned i = 0; i < pt_over_areas.size() ; i++ ) { |
---|
344 | double ratio = pt_over_areas[i]; |
---|
345 | trunc_sum += ratio; |
---|
346 | trunc_sumsqr += ratio*ratio; |
---|
347 | means[i] = trunc_sum / (i+1); |
---|
348 | sd[i] = sqrt(abs(means[i]*means[i] - trunc_sumsqr/(i+1))); |
---|
349 | cerr << "i, means, sd: " <<i<<", "<< means[i] <<", "<<sd[i]<<", "<< |
---|
350 | sd[i]/sqrt(i+1.0)<<endl; |
---|
351 | } |
---|
352 | cout << "-----------------------------------"<<endl; |
---|
353 | for (unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) { |
---|
354 | cout << "Median "<< i <<" = " << pt_over_areas[i]<<endl; |
---|
355 | } |
---|
356 | cout << "Number of non-jets: "<<_non_jet_number<<endl; |
---|
357 | cout << "Area of non-jets: "<<_non_jet_area<<endl; |
---|
358 | cout << "Default median position: " << (pt_over_areas.size()-1)/2.0<<endl; |
---|
359 | cout << "NJ median position: " << nj_median_pos <<endl; |
---|
360 | cout << "NJ median value: " << nj_median_ratio <<endl; |
---|
361 | return 0.0; |
---|
362 | } |
---|
363 | |
---|
364 | switch(strat) { |
---|
365 | case median: |
---|
366 | case median_4vector: |
---|
367 | return nj_median_ratio; |
---|
368 | case non_ghost_median: |
---|
369 | return non_ghost_median_ratio; |
---|
370 | case pttot_over_areatot: |
---|
371 | return pt_sum / area_sum; |
---|
372 | case pttot_over_areatot_cut: |
---|
373 | return pt_sum_with_cut / area_sum_with_cut; |
---|
374 | case mean_ratio_cut: |
---|
375 | return ratio_sum/ratio_n; |
---|
376 | default: |
---|
377 | return nj_median_ratio; |
---|
378 | } |
---|
379 | |
---|
380 | } |
---|
381 | |
---|
382 | |
---|
383 | // The following functionality is now provided by the base class |
---|
384 | // //---------------------------------------------------------------------- |
---|
385 | // // fit a parabola to pt/area as a function of rapidity, using the |
---|
386 | // // formulae of CCN28-36 (which actually fits f = a+b*x^2) |
---|
387 | // void ClusterSequenceActiveArea::parabolic_pt_per_unit_area( |
---|
388 | // double & a, double & b, double raprange, double exclude_above, |
---|
389 | // bool use_area_4vector) const { |
---|
390 | // |
---|
391 | // double this_raprange; |
---|
392 | // if (raprange <= 0) {this_raprange = _safe_rap_for_area;} |
---|
393 | // else {this_raprange = raprange;} |
---|
394 | // |
---|
395 | // int n=0; |
---|
396 | // int n_excluded = 0; |
---|
397 | // double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0; |
---|
398 | // |
---|
399 | // vector<PseudoJet> incl_jets = inclusive_jets(); |
---|
400 | // |
---|
401 | // for (unsigned i = 0; i < incl_jets.size(); i++) { |
---|
402 | // if (abs(incl_jets[i].rap()) < this_raprange) { |
---|
403 | // double this_area; |
---|
404 | // if ( use_area_4vector ) { |
---|
405 | // this_area = area_4vector(incl_jets[i]).perp(); |
---|
406 | // } else { |
---|
407 | // this_area = area(incl_jets[i]); |
---|
408 | // } |
---|
409 | // double f = incl_jets[i].perp()/this_area; |
---|
410 | // if (exclude_above <= 0.0 || f < exclude_above) { |
---|
411 | // double x = incl_jets[i].rap(); double x2 = x*x; |
---|
412 | // mean_f += f; |
---|
413 | // mean_x2 += x2; |
---|
414 | // mean_x4 += x2*x2; |
---|
415 | // mean_fx2 += f*x2; |
---|
416 | // n++; |
---|
417 | // } else { |
---|
418 | // n_excluded++; |
---|
419 | // } |
---|
420 | // } |
---|
421 | // } |
---|
422 | // |
---|
423 | // if (n <= 1) { |
---|
424 | // // meaningful results require at least two jets inside the |
---|
425 | // // area -- mind you if there are empty jets we should be in |
---|
426 | // // any case doing something special... |
---|
427 | // a = 0.0; |
---|
428 | // b = 0.0; |
---|
429 | // } else { |
---|
430 | // mean_f /= n; |
---|
431 | // mean_x2 /= n; |
---|
432 | // mean_x4 /= n; |
---|
433 | // mean_fx2 /= n; |
---|
434 | // |
---|
435 | // b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4); |
---|
436 | // a = mean_f - b*mean_x2; |
---|
437 | // } |
---|
438 | // //cerr << "n_excluded = "<< n_excluded << endl; |
---|
439 | // } |
---|
440 | |
---|
441 | |
---|
442 | //---------------------------------------------------------------------- |
---|
443 | double ClusterSequenceActiveArea::empty_area(const Selector & selector) const { |
---|
444 | // make sure that the selector applies jet by jet |
---|
445 | if (! selector.applies_jet_by_jet()){ |
---|
446 | throw Error("ClusterSequenceActiveArea: empty area can only be computed from selectors applying jet by jet"); |
---|
447 | } |
---|
448 | |
---|
449 | double empty = 0.0; |
---|
450 | // first deal with ghost jets |
---|
451 | for (unsigned i = 0; i < _ghost_jets.size(); i++) { |
---|
452 | if (selector.pass(_ghost_jets[i])) { |
---|
453 | empty += _ghost_jets[i].area; |
---|
454 | } |
---|
455 | } |
---|
456 | // then deal with unclustered ghosts |
---|
457 | for (unsigned i = 0; i < _unclustered_ghosts.size(); i++) { |
---|
458 | if (selector.pass(_unclustered_ghosts[i])) { |
---|
459 | empty += _unclustered_ghosts[i].area; |
---|
460 | } |
---|
461 | } |
---|
462 | empty /= _ghost_spec_repeat; |
---|
463 | return empty; |
---|
464 | } |
---|
465 | |
---|
466 | //---------------------------------------------------------------------- |
---|
467 | double ClusterSequenceActiveArea::n_empty_jets(const Selector & selector) const { |
---|
468 | _check_selector_good_for_median(selector); |
---|
469 | |
---|
470 | double inrange = 0; |
---|
471 | for (unsigned i = 0; i < _ghost_jets.size(); i++) { |
---|
472 | if (selector.pass(_ghost_jets[i])) inrange++; |
---|
473 | } |
---|
474 | inrange /= _ghost_spec_repeat; |
---|
475 | return inrange; |
---|
476 | } |
---|
477 | |
---|
478 | //---------------------------------------------------------------------- |
---|
479 | /// transfer the history (and jet-momenta) from clust_seq to our |
---|
480 | /// own internal structure while removing ghosts |
---|
481 | void ClusterSequenceActiveArea::_transfer_ghost_free_history( |
---|
482 | const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq) { |
---|
483 | |
---|
484 | const vector<history_element> & gs_history = ghosted_seq.history(); |
---|
485 | vector<int> gs2self_hist_map(gs_history.size()); |
---|
486 | |
---|
487 | // first transfer info about strategy used (which isn't necessarily |
---|
488 | // always the one that got asked for...) |
---|
489 | _strategy = ghosted_seq.strategy_used(); |
---|
490 | |
---|
491 | // work our way through to first non-trivial combination |
---|
492 | unsigned igs = 0; |
---|
493 | unsigned iself = 0; |
---|
494 | while (igs < gs_history.size() && gs_history[igs].parent1 == InexistentParent) { |
---|
495 | // record correspondence |
---|
496 | if (!ghosted_seq.is_pure_ghost(igs)) { |
---|
497 | gs2self_hist_map[igs] = iself++; |
---|
498 | } else { |
---|
499 | gs2self_hist_map[igs] = Invalid; |
---|
500 | } |
---|
501 | igs++; |
---|
502 | }; |
---|
503 | |
---|
504 | // make sure the count of non-ghost initial jets is equal to |
---|
505 | // what we already have in terms of initial jets |
---|
506 | assert(iself == _history.size()); |
---|
507 | |
---|
508 | // if there was no clustering in this event (e.g. SISCone passive |
---|
509 | // area with zero input particles, or with a pt cut on stable cones |
---|
510 | // that kills all jets), then don't bother with the rest (which |
---|
511 | // would crash!) |
---|
512 | if (igs == gs_history.size()) return; |
---|
513 | |
---|
514 | // now actually transfer things |
---|
515 | do { |
---|
516 | // if we are a pure ghost, then go on to next round |
---|
517 | if (ghosted_seq.is_pure_ghost(igs)) { |
---|
518 | gs2self_hist_map[igs] = Invalid; |
---|
519 | continue; |
---|
520 | } |
---|
521 | |
---|
522 | const history_element & gs_hist_el = gs_history[igs]; |
---|
523 | |
---|
524 | bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1); |
---|
525 | bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2); |
---|
526 | |
---|
527 | // if exactly one parent is a ghost then maintain info about the |
---|
528 | // non-ghost correspondence for this jet, and then go on to next |
---|
529 | // recombination in the ghosted sequence |
---|
530 | if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) { |
---|
531 | gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2]; |
---|
532 | continue; |
---|
533 | } |
---|
534 | if (!parent1_is_ghost && parent2_is_ghost) { |
---|
535 | gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1]; |
---|
536 | continue; |
---|
537 | } |
---|
538 | |
---|
539 | // no parents are ghosts... |
---|
540 | if (gs_hist_el.parent2 >= 0) { |
---|
541 | // recombination of two non-ghosts |
---|
542 | gs2self_hist_map[igs] = _history.size(); |
---|
543 | // record the recombination in our own sequence |
---|
544 | int newjet_k; // dummy var -- not used |
---|
545 | int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index; |
---|
546 | int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index; |
---|
547 | //cerr << "recombining "<< jet_i << " and "<< jet_j << endl; |
---|
548 | _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k); |
---|
549 | } else { |
---|
550 | // we have a non-ghost that has become a beam-jet |
---|
551 | assert(gs_history[igs].parent2 == BeamJet); |
---|
552 | // record position |
---|
553 | gs2self_hist_map[igs] = _history.size(); |
---|
554 | // record the recombination in our own sequence |
---|
555 | _do_iB_recombination_step( |
---|
556 | _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index, |
---|
557 | gs_hist_el.dij); |
---|
558 | } |
---|
559 | } while (++igs < gs_history.size()); |
---|
560 | |
---|
561 | } |
---|
562 | |
---|
563 | //---------------------------------------------------------------------- |
---|
564 | void ClusterSequenceActiveArea::_transfer_areas( |
---|
565 | const vector<int> & unique_hist_order, |
---|
566 | const ClusterSequenceActiveAreaExplicitGhosts & ghosted_seq ) { |
---|
567 | |
---|
568 | const vector<history_element> & gs_history = ghosted_seq.history(); |
---|
569 | const vector<PseudoJet> & gs_jets = ghosted_seq.jets(); |
---|
570 | vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order(); |
---|
571 | |
---|
572 | const double tolerance = 1e-11; // to decide when two jets are the same |
---|
573 | |
---|
574 | int j = -1; |
---|
575 | int hist_index = -1; |
---|
576 | |
---|
577 | valarray<double> our_areas(_history.size()); |
---|
578 | our_areas = 0.0; |
---|
579 | |
---|
580 | valarray<PseudoJet> our_area_4vectors(_history.size()); |
---|
581 | our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0); |
---|
582 | |
---|
583 | for (unsigned i = 0; i < gs_history.size(); i++) { |
---|
584 | // only consider composite particles |
---|
585 | unsigned gs_hist_index = gs_unique_hist_order[i]; |
---|
586 | if (gs_hist_index < ghosted_seq.n_particles()) continue; |
---|
587 | const history_element & gs_hist = gs_history[gs_unique_hist_order[i]]; |
---|
588 | int parent1 = gs_hist.parent1; |
---|
589 | int parent2 = gs_hist.parent2; |
---|
590 | |
---|
591 | if (parent2 == BeamJet) { |
---|
592 | // need to look at parent to get the actual jet |
---|
593 | const PseudoJet & jet = |
---|
594 | gs_jets[gs_history[parent1].jetp_index]; |
---|
595 | double area_local = ghosted_seq.area(jet); |
---|
596 | PseudoJet ext_area = ghosted_seq.area_4vector(jet); |
---|
597 | |
---|
598 | if (ghosted_seq.is_pure_ghost(parent1)) { |
---|
599 | // record the existence of the pure ghost jet for future use |
---|
600 | _ghost_jets.push_back(GhostJet(jet,area_local)); |
---|
601 | if (abs(jet.rap()) < _safe_rap_for_area) { |
---|
602 | _non_jet_area += area_local; |
---|
603 | _non_jet_area2 += area_local*area_local; |
---|
604 | _non_jet_number += 1; |
---|
605 | } |
---|
606 | } else { |
---|
607 | |
---|
608 | // get next "combined-particle" index in our own history |
---|
609 | // making sure we don't go beyond its bounds (if we do |
---|
610 | // then we're in big trouble anyway...) |
---|
611 | while (++j < int(_history.size())) { |
---|
612 | hist_index = unique_hist_order[j]; |
---|
613 | if (hist_index >= _initial_n) break;} |
---|
614 | |
---|
615 | // sanity checking -- do not overrun |
---|
616 | if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching"); |
---|
617 | |
---|
618 | // sanity check -- make sure we are taking about the same |
---|
619 | // jet in reference and new sequences |
---|
620 | int refjet_index = _history[_history[hist_index].parent1].jetp_index; |
---|
621 | assert(refjet_index >= 0 && refjet_index < int(_jets.size())); |
---|
622 | const PseudoJet & refjet = _jets[refjet_index]; |
---|
623 | |
---|
624 | //cerr << "Inclusive" << endl; |
---|
625 | //cerr << gs_history[parent1].jetp_index << " " << gs_jets.size() << endl; |
---|
626 | //cerr << _history[_history[hist_index].parent1].jetp_index << " " << _jets.size() << endl; |
---|
627 | |
---|
628 | // If pt disagrees check E; if they both disagree there's a |
---|
629 | // problem here... NB: a massive particle with zero pt may |
---|
630 | // have its pt changed when a ghost is added -- this is why we |
---|
631 | // also require the energy to be wrong before complaining |
---|
632 | _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance, |
---|
633 | ghosted_seq); |
---|
634 | |
---|
635 | // set the area at this clustering stage |
---|
636 | our_areas[hist_index] = area_local; |
---|
637 | our_area_4vectors[hist_index] = ext_area; |
---|
638 | |
---|
639 | // update the parent as well -- that way its area is the area |
---|
640 | // immediately before clustering (i.e. resolve an ambiguity in |
---|
641 | // the Cambridge case and ensure in the kt case that the original |
---|
642 | // particles get a correct area) |
---|
643 | our_areas[_history[hist_index].parent1] = area_local; |
---|
644 | our_area_4vectors[_history[hist_index].parent1] = ext_area; |
---|
645 | |
---|
646 | } |
---|
647 | } |
---|
648 | else if (!ghosted_seq.is_pure_ghost(parent1) && |
---|
649 | !ghosted_seq.is_pure_ghost(parent2)) { |
---|
650 | |
---|
651 | // get next "combined-particle" index in our own history |
---|
652 | while (++j < int(_history.size())) { |
---|
653 | hist_index = unique_hist_order[j]; |
---|
654 | if (hist_index >= _initial_n) break;} |
---|
655 | |
---|
656 | // sanity checking -- do not overrun |
---|
657 | if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching"); |
---|
658 | |
---|
659 | // make sure that our reference history entry is also for |
---|
660 | // an exclusive (dij) clustering (otherwise the comparison jet |
---|
661 | // will not exist) |
---|
662 | if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)"); |
---|
663 | |
---|
664 | //cerr << "Exclusive: hist_index,hist_size: " << hist_index << " " << _history.size()<< endl; |
---|
665 | //cerr << gs_hist.jetp_index << " " << gs_jets.size() << endl; |
---|
666 | //cerr << _history[hist_index].jetp_index << " " << _jets.size() << endl; |
---|
667 | |
---|
668 | const PseudoJet & jet = gs_jets[gs_hist.jetp_index]; |
---|
669 | const PseudoJet & refjet = _jets[_history[hist_index].jetp_index]; |
---|
670 | |
---|
671 | // run sanity check |
---|
672 | _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance, |
---|
673 | ghosted_seq); |
---|
674 | |
---|
675 | // update area and our local index (maybe redundant since later |
---|
676 | // the descendants will reupdate it?) |
---|
677 | double area_local = ghosted_seq.area(jet); |
---|
678 | our_areas[hist_index] += area_local; |
---|
679 | |
---|
680 | PseudoJet ext_area = ghosted_seq.area_4vector(jet); |
---|
681 | |
---|
682 | // GPS TMP debugging (jetclu) ----------------------- |
---|
683 | //ext_area = PseudoJet(1e-100,1e-100,1e-100,4e-100); |
---|
684 | //our_area_4vectors[hist_index] = ext_area; |
---|
685 | //cout << "aa " |
---|
686 | // << our_area_4vectors[hist_index].px() << " " |
---|
687 | // << our_area_4vectors[hist_index].py() << " " |
---|
688 | // << our_area_4vectors[hist_index].pz() << " " |
---|
689 | // << our_area_4vectors[hist_index].E() << endl; |
---|
690 | //cout << "bb " |
---|
691 | // << ext_area.px() << " " |
---|
692 | // << ext_area.py() << " " |
---|
693 | // << ext_area.pz() << " " |
---|
694 | // << ext_area.E() << endl; |
---|
695 | //--------------------------------------------------- |
---|
696 | |
---|
697 | _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area); |
---|
698 | |
---|
699 | // now update areas of parents (so that they becomes areas |
---|
700 | // immediately before clustering occurred). This is of use |
---|
701 | // because it allows us to set the areas of the original hard |
---|
702 | // particles in the kt algorithm; for the Cambridge case it |
---|
703 | // means a jet's area will be the area just before it clusters |
---|
704 | // with another hard jet. |
---|
705 | const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index]; |
---|
706 | int our_parent1 = _history[hist_index].parent1; |
---|
707 | our_areas[our_parent1] = ghosted_seq.area(jet1); |
---|
708 | our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1); |
---|
709 | |
---|
710 | const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index]; |
---|
711 | int our_parent2 = _history[hist_index].parent2; |
---|
712 | our_areas[our_parent2] = ghosted_seq.area(jet2); |
---|
713 | our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2); |
---|
714 | } |
---|
715 | |
---|
716 | } |
---|
717 | |
---|
718 | // now add unclustered ghosts to the relevant list so that we can |
---|
719 | // calculate empty area later. |
---|
720 | vector<PseudoJet> unclust = ghosted_seq.unclustered_particles(); |
---|
721 | for (unsigned iu = 0; iu < unclust.size(); iu++) { |
---|
722 | if (ghosted_seq.is_pure_ghost(unclust[iu])) { |
---|
723 | double area_local = ghosted_seq.area(unclust[iu]); |
---|
724 | _unclustered_ghosts.push_back(GhostJet(unclust[iu],area_local)); |
---|
725 | } |
---|
726 | } |
---|
727 | |
---|
728 | /* |
---|
729 | * WARNING: |
---|
730 | * _average_area has explicitly been sized initially to 2*jets().size() |
---|
731 | * which can be bigger than our_areas (of size _history.size() |
---|
732 | * if there are some unclustered particles. |
---|
733 | * So we must take care about boundaries |
---|
734 | */ |
---|
735 | |
---|
736 | for (unsigned int area_index = 0; area_index<our_areas.size(); area_index++){ |
---|
737 | _average_area[area_index] += our_areas[area_index]; |
---|
738 | _average_area2[area_index] += our_areas[area_index]*our_areas[area_index]; |
---|
739 | } |
---|
740 | |
---|
741 | //_average_area_4vector += our_area_4vectors; |
---|
742 | // Use the proper recombination scheme when averaging the area_4vectors |
---|
743 | // over multiple ghost runs (i.e. the repeat stage); |
---|
744 | for (unsigned i = 0; i < our_area_4vectors.size(); i++) { |
---|
745 | _jet_def.recombiner()->plus_equal(_average_area_4vector[i], |
---|
746 | our_area_4vectors[i]); |
---|
747 | } |
---|
748 | } |
---|
749 | |
---|
750 | |
---|
751 | /// check if two jets have the same momentum to within the |
---|
752 | /// tolerance (and if pt's are not the same we're forgiving and |
---|
753 | /// look to see if the energy is the same) |
---|
754 | void ClusterSequenceActiveArea::_throw_unless_jets_have_same_perp_or_E( |
---|
755 | const PseudoJet & jet, |
---|
756 | const PseudoJet & refjet, |
---|
757 | double tolerance, |
---|
758 | const ClusterSequenceActiveAreaExplicitGhosts & jets_ghosted_seq |
---|
759 | ) const { |
---|
760 | |
---|
761 | if (abs(jet.perp2()-refjet.perp2()) > |
---|
762 | tolerance*max(jet.perp2(),refjet.perp2()) |
---|
763 | && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) { |
---|
764 | ostringstream ostr; |
---|
765 | ostr << "Could not match clustering sequence for an inclusive/exclusive jet when reconstructing areas. See FAQ for possible explanations." << endl; |
---|
766 | ostr << " Ref-Jet: " |
---|
767 | << refjet.px() << " " |
---|
768 | << refjet.py() << " " |
---|
769 | << refjet.pz() << " " |
---|
770 | << refjet.E() << endl; |
---|
771 | ostr << " New-Jet: " |
---|
772 | << jet.px() << " " |
---|
773 | << jet.py() << " " |
---|
774 | << jet.pz() << " " |
---|
775 | << jet.E() << endl; |
---|
776 | if (jets_ghosted_seq.has_dangerous_particles()) { |
---|
777 | ostr << " NB: some particles have pt too low wrt ghosts -- this may be the cause" << endl;} |
---|
778 | //ostr << jet.perp() << " " << refjet.perp() << " " |
---|
779 | // << jet.perp() - refjet.perp() << endl; |
---|
780 | throw Error(ostr.str()); |
---|
781 | } |
---|
782 | } |
---|
783 | |
---|
784 | FASTJET_END_NAMESPACE |
---|
785 | |
---|