Fork me on GitHub

source: svn/trunk/Utilities/Fastjet/src/ClusterSequenceActiveArea.cc@ 598

Last change on this file since 598 was 11, checked in by severine ovyn, 16 years ago

Fastjet added; CDFCones directory has been changed

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