#include <ClusterSequenceActiveArea.hh>
Inheritance diagram for fastjet::ClusterSequenceActiveArea:
Public Types | |
enum | mean_pt_strategies { median = 0, non_ghost_median, pttot_over_areatot, pttot_over_areatot_cut, mean_ratio_cut, play, median_4vector } |
enum providing a variety of tentative strategies for estimating the background (e.g. More... | |
Public Member Functions | |
ClusterSequenceActiveArea () | |
default constructor | |
template<class L> | |
ClusterSequenceActiveArea (const std::vector< L > &pseudojets, const JetDefinition &jet_def, const GhostedAreaSpec &area_spec, const bool &writeout_combinations=false) | |
constructor based on JetDefinition and GhostedAreaSpec | |
virtual double | area (const PseudoJet &jet) const |
return the area associated with the given jet; this base class returns 0. | |
virtual double | area_error (const PseudoJet &jet) const |
return the error (uncertainty) associated with the determination of the area of this jet; this base class returns 0. | |
virtual PseudoJet | area_4vector (const PseudoJet &jet) const |
return a PseudoJet whose 4-vector is defined by the following integral | |
double | pt_per_unit_area (mean_pt_strategies strat=median, double range=2.0) const |
return the transverse momentum per unit area according to one of the above strategies; for some strategies (those with "cut" in their name) the parameter "range" allows one to exclude a subset of the jets for the background estimation, those that have pt/area > median(pt/area)*range. | |
virtual double | empty_area (const RangeDefinition &range) const |
rewrite the empty area from the parent class, so as to use all info at our disposal return the total area, in the given y-phi range, that consists of ghost jets or unclustered ghosts | |
virtual double | n_empty_jets (const RangeDefinition &range) const |
return the true number of empty jets (replaces ClusterSequenceAreaBase::n_empty_jets(. | |
Protected Member Functions | |
void | _resize_and_zero_AA () |
void | _initialise_AA (const JetDefinition &jet_def, const GhostedAreaSpec &area_spec, const bool &writeout_combinations, bool &continue_running) |
void | _run_AA (const GhostedAreaSpec &area_spec) |
void | _postprocess_AA (const GhostedAreaSpec &area_spec) |
run the postprocessing for the active area (and derived classes) | |
void | _initialise_and_run_AA (const JetDefinition &jet_def, const GhostedAreaSpec &area_spec, const bool &writeout_combinations=false) |
global routine for running active area | |
void | _transfer_ghost_free_history (const ClusterSequenceActiveAreaExplicitGhosts &clust_seq) |
transfer the history (and jet-momenta) from clust_seq to our own internal structure while removing ghosts | |
void | _transfer_areas (const vector< int > &unique_hist_order, const ClusterSequenceActiveAreaExplicitGhosts &) |
transfer areas from the ClusterSequenceActiveAreaExplicitGhosts object into our internal area bookkeeping. | |
bool | has_dangerous_particles () const |
returns true if there are any particles whose transverse momentum if so low that there's a risk of the ghosts having modified the clustering sequence | |
Protected Attributes | |
valarray< double > | _average_area |
child classes benefit from having these at their disposal | |
valarray< double > | _average_area2 |
valarray< PseudoJet > | _average_area_4vector |
Private Member Functions | |
void | _extract_tree (vector< int > &) const |
routine for extracting the tree in an order that will be independent of any degeneracies in the recombination sequence that don't affect the composition of the final jets | |
void | _extract_tree_children (int pos, valarray< bool > &, const valarray< int > &, vector< int > &) const |
do the part of the extraction associated with pos, working through its children and their parents | |
void | _extract_tree_parents (int pos, valarray< bool > &, const valarray< int > &, vector< int > &) const |
do the part of the extraction associated with the parents of pos. | |
void | _throw_unless_jets_have_same_perp_or_E (const PseudoJet &jet, const PseudoJet &refjet, double tolerance, const ClusterSequenceActiveAreaExplicitGhosts &jets_ghosted_seq) const |
check if two jets have the same momentum to within the tolerance (and if pt's are not the same we're forgiving and look to see if the energy is the same) | |
Private Attributes | |
double | _non_jet_area |
double | _non_jet_area2 |
double | _non_jet_number |
double | _maxrap_for_area |
double | _safe_rap_for_area |
bool | _has_dangerous_particles |
int | _area_spec_repeat |
since we are playing nasty games with seeds, we should warn the user a few times | |
std::vector< GhostJet > | _ghost_jets |
std::vector< GhostJet > | _unclustered_ghosts |
Classes | |
class | GhostJet |
a class for our internal storage of ghost jets More... |
.. Figure out what to do about seeds later...)
Definition at line 56 of file ClusterSequenceActiveArea.hh.
enum providing a variety of tentative strategies for estimating the background (e.g.
non-jet) activity in a highly populated event; the one that has been most extensively tested is median.
median | |
non_ghost_median | |
pttot_over_areatot | |
pttot_over_areatot_cut | |
mean_ratio_cut | |
play | |
median_4vector |
Definition at line 80 of file ClusterSequenceActiveArea.hh.
00080 {median=0, non_ghost_median, pttot_over_areatot, 00081 pttot_over_areatot_cut, mean_ratio_cut, play, 00082 median_4vector};
fastjet::ClusterSequenceActiveArea::ClusterSequenceActiveArea | ( | ) | [inline] |
fastjet::ClusterSequenceActiveArea::ClusterSequenceActiveArea | ( | const std::vector< L > & | pseudojets, | |
const JetDefinition & | jet_def, | |||
const GhostedAreaSpec & | area_spec, | |||
const bool & | writeout_combinations = false | |||
) | [inline] |
constructor based on JetDefinition and GhostedAreaSpec
Definition at line 204 of file ClusterSequenceActiveArea.hh.
00207 { 00208 00209 // transfer the initial jets (type L) into our own array 00210 _transfer_input_jets(pseudojets); 00211 00212 // run the clustering for active areas 00213 _initialise_and_run_AA(jet_def, area_spec, writeout_combinations); 00214 00215 }
virtual double fastjet::ClusterSequenceActiveArea::area | ( | const PseudoJet & | jet | ) | const [inline, virtual] |
return the area associated with the given jet; this base class returns 0.
Reimplemented from fastjet::ClusterSequenceAreaBase.
Definition at line 69 of file ClusterSequenceActiveArea.hh.
References fastjet::PseudoJet::cluster_hist_index().
Referenced by _transfer_areas(), empty_area(), and pt_per_unit_area().
00069 { 00070 return _average_area[jet.cluster_hist_index()];};
virtual double fastjet::ClusterSequenceActiveArea::area_error | ( | const PseudoJet & | jet | ) | const [inline, virtual] |
return the error (uncertainty) associated with the determination of the area of this jet; this base class returns 0.
Reimplemented from fastjet::ClusterSequenceAreaBase.
Definition at line 71 of file ClusterSequenceActiveArea.hh.
References fastjet::PseudoJet::cluster_hist_index().
00071 { 00072 return _average_area2[jet.cluster_hist_index()];};
virtual PseudoJet fastjet::ClusterSequenceActiveArea::area_4vector | ( | const PseudoJet & | jet | ) | const [inline, virtual] |
return a PseudoJet whose 4-vector is defined by the following integral
drap d PseudoJet("rap,phi,pt=one") * Theta("rap,phi inside jet boundary")
where PseudoJet("rap,phi,pt=one") is a 4-vector with the given rapidity (rap), azimuth (phi) and pt=1, while Theta("rap,phi inside jet boundary") is a function that is 1 when rap,phi define a direction inside the jet boundary and 0 otherwise.
This base class returns a null 4-vector.
Reimplemented from fastjet::ClusterSequenceAreaBase.
Definition at line 74 of file ClusterSequenceActiveArea.hh.
References fastjet::PseudoJet::cluster_hist_index().
Referenced by pt_per_unit_area().
00074 { 00075 return _average_area_4vector[jet.cluster_hist_index()];};
double fastjet::ClusterSequenceActiveArea::pt_per_unit_area | ( | mean_pt_strategies | strat = median , |
|
double | range = 2.0 | |||
) | const |
return the transverse momentum per unit area according to one of the above strategies; for some strategies (those with "cut" in their name) the parameter "range" allows one to exclude a subset of the jets for the background estimation, those that have pt/area > median(pt/area)*range.
NB: This call is OBSOLETE; use media_pt_per_unit_area from the
Definition at line 271 of file ClusterSequenceActiveArea.cc.
References _non_jet_area, _non_jet_number, _safe_rap_for_area, area(), area_4vector(), fastjet::ClusterSequence::inclusive_jets(), mean_ratio_cut, median, median_4vector, non_ghost_median, fastjet::PseudoJet::perp(), play, pttot_over_areatot, and pttot_over_areatot_cut.
00272 { 00273 00274 vector<PseudoJet> incl_jets = inclusive_jets(); 00275 vector<double> pt_over_areas; 00276 00277 for (unsigned i = 0; i < incl_jets.size(); i++) { 00278 if (abs(incl_jets[i].rap()) < _safe_rap_for_area) { 00279 double this_area; 00280 if ( strat == median_4vector ) { 00281 this_area = area_4vector(incl_jets[i]).perp(); 00282 } else { 00283 this_area = area(incl_jets[i]); 00284 } 00285 pt_over_areas.push_back(incl_jets[i].perp()/this_area); 00286 } 00287 } 00288 00289 // there is nothing inside our region, so answer will always be zero 00290 if (pt_over_areas.size() == 0) {return 0.0;} 00291 00292 // get median (pt/area) [this is the "old" median definition. It considers 00293 // only the "real" jets in calculating the median, i.e. excluding the 00294 // only-ghost ones] 00295 sort(pt_over_areas.begin(), pt_over_areas.end()); 00296 double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2]; 00297 00298 // new median definition that takes into account non-jet area (i.e. 00299 // jets composed only of ghosts), and for fractional median position 00300 // interpolates between the corresponding entries in the pt_over_areas array 00301 double nj_median_pos = (pt_over_areas.size()-1 - _non_jet_number)/2.0; 00302 double nj_median_ratio; 00303 if (nj_median_pos >= 0 && pt_over_areas.size() > 1) { 00304 int int_nj_median = int(nj_median_pos); 00305 nj_median_ratio = 00306 pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos) 00307 + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median); 00308 } else { 00309 nj_median_ratio = 0.0; 00310 } 00311 00312 00313 // get various forms of mean (pt/area) 00314 double pt_sum = 0.0, pt_sum_with_cut = 0.0; 00315 double area_sum = _non_jet_area, area_sum_with_cut = _non_jet_area; 00316 double ratio_sum = 0.0; 00317 double ratio_n = _non_jet_number; 00318 for (unsigned i = 0; i < incl_jets.size(); i++) { 00319 if (abs(incl_jets[i].rap()) < _safe_rap_for_area) { 00320 double this_area; 00321 if ( strat == median_4vector ) { 00322 this_area = area_4vector(incl_jets[i]).perp(); 00323 } else { 00324 this_area = area(incl_jets[i]); 00325 } 00326 pt_sum += incl_jets[i].perp(); 00327 area_sum += this_area; 00328 double ratio = incl_jets[i].perp()/this_area; 00329 if (ratio < range*nj_median_ratio) { 00330 pt_sum_with_cut += incl_jets[i].perp(); 00331 area_sum_with_cut += this_area; 00332 ratio_sum += ratio; ratio_n++; 00333 } 00334 } 00335 } 00336 00337 if (strat == play) { 00338 double trunc_sum = 0, trunc_sumsqr = 0; 00339 vector<double> means(pt_over_areas.size()), sd(pt_over_areas.size()); 00340 for (unsigned i = 0; i < pt_over_areas.size() ; i++ ) { 00341 double ratio = pt_over_areas[i]; 00342 trunc_sum += ratio; 00343 trunc_sumsqr += ratio*ratio; 00344 means[i] = trunc_sum / (i+1); 00345 sd[i] = sqrt(abs(means[i]*means[i] - trunc_sumsqr/(i+1))); 00346 cerr << "i, means, sd: " <<i<<", "<< means[i] <<", "<<sd[i]<<", "<< 00347 sd[i]/sqrt(i+1.0)<<endl; 00348 } 00349 cout << "-----------------------------------"<<endl; 00350 for (unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) { 00351 cout << "Median "<< i <<" = " << pt_over_areas[i]<<endl; 00352 } 00353 cout << "Number of non-jets: "<<_non_jet_number<<endl; 00354 cout << "Area of non-jets: "<<_non_jet_area<<endl; 00355 cout << "Default median position: " << (pt_over_areas.size()-1)/2.0<<endl; 00356 cout << "NJ median position: " << nj_median_pos <<endl; 00357 cout << "NJ median value: " << nj_median_ratio <<endl; 00358 return 0.0; 00359 } 00360 00361 switch(strat) { 00362 case median: 00363 case median_4vector: 00364 return nj_median_ratio; 00365 case non_ghost_median: 00366 return non_ghost_median_ratio; 00367 case pttot_over_areatot: 00368 return pt_sum / area_sum; 00369 case pttot_over_areatot_cut: 00370 return pt_sum_with_cut / area_sum_with_cut; 00371 case mean_ratio_cut: 00372 return ratio_sum/ratio_n; 00373 default: 00374 return nj_median_ratio; 00375 } 00376 00377 }
double fastjet::ClusterSequenceActiveArea::empty_area | ( | const RangeDefinition & | range | ) | const [virtual] |
rewrite the empty area from the parent class, so as to use all info at our disposal return the total area, in the given y-phi range, that consists of ghost jets or unclustered ghosts
Reimplemented from fastjet::ClusterSequenceAreaBase.
Reimplemented in fastjet::ClusterSequencePassiveArea.
Definition at line 440 of file ClusterSequenceActiveArea.cc.
References _area_spec_repeat, _ghost_jets, _unclustered_ghosts, area(), and fastjet::RangeDefinition::is_in_range().
Referenced by fastjet::ClusterSequencePassiveArea::empty_area().
00440 { 00441 double empty = 0.0; 00442 // first deal with ghost jets 00443 for (unsigned i = 0; i < _ghost_jets.size(); i++) { 00444 if (range.is_in_range(_ghost_jets[i])) { 00445 empty += _ghost_jets[i].area; 00446 } 00447 } 00448 // then deal with unclustered ghosts 00449 for (unsigned i = 0; i < _unclustered_ghosts.size(); i++) { 00450 if (range.is_in_range(_unclustered_ghosts[i])) { 00451 empty += _unclustered_ghosts[i].area; 00452 } 00453 } 00454 empty /= _area_spec_repeat; 00455 return empty; 00456 }
double fastjet::ClusterSequenceActiveArea::n_empty_jets | ( | const RangeDefinition & | range | ) | const [virtual] |
return the true number of empty jets (replaces ClusterSequenceAreaBase::n_empty_jets(.
..))
Reimplemented from fastjet::ClusterSequenceAreaBase.
Reimplemented in fastjet::ClusterSequence1GhostPassiveArea.
Definition at line 459 of file ClusterSequenceActiveArea.cc.
References _area_spec_repeat, _ghost_jets, and fastjet::RangeDefinition::is_in_range().
00459 { 00460 double inrange = 0; 00461 for (unsigned i = 0; i < _ghost_jets.size(); i++) { 00462 if (range.is_in_range(_ghost_jets[i])) inrange++; 00463 } 00464 inrange /= _area_spec_repeat; 00465 return inrange; 00466 }
void fastjet::ClusterSequenceActiveArea::_resize_and_zero_AA | ( | ) | [protected] |
Definition at line 64 of file ClusterSequenceActiveArea.cc.
References _average_area, _average_area2, _average_area_4vector, fastjet::ClusterSequence::_jets, _non_jet_area, _non_jet_area2, and _non_jet_number.
Referenced by _initialise_AA(), and fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA().
00064 { 00065 // initialize our local area information 00066 _average_area.resize(2*_jets.size()); _average_area = 0.0; 00067 _average_area2.resize(2*_jets.size()); _average_area2 = 0.0; 00068 _average_area_4vector.resize(2*_jets.size()); 00069 _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0); 00070 _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0; 00071 }
void fastjet::ClusterSequenceActiveArea::_initialise_AA | ( | const JetDefinition & | jet_def, | |
const GhostedAreaSpec & | area_spec, | |||
const bool & | writeout_combinations, | |||
bool & | continue_running | |||
) | [protected] |
Definition at line 74 of file ClusterSequenceActiveArea.cc.
References _area_spec_repeat, fastjet::ClusterSequence::_decant_options(), fastjet::ClusterSequence::_fill_initial_history(), _has_dangerous_particles, fastjet::ClusterSequence::_initialise_and_run(), _maxrap_for_area, _resize_and_zero_AA(), _safe_rap_for_area, fastjet::GhostedAreaSpec::ghost_maxrap(), fastjet::ClusterSequence::jet_def(), fastjet::JetDefinition::R(), and fastjet::GhostedAreaSpec::repeat().
Referenced by fastjet::ClusterSequence1GhostPassiveArea::_initialise_and_run_1GPA(), and _initialise_and_run_AA().
00079 { 00080 00081 // store this for future use 00082 _area_spec_repeat = area_spec.repeat(); 00083 00084 // make sure placeholders are there & zeroed 00085 _resize_and_zero_AA(); 00086 00087 // for future reference... 00088 _maxrap_for_area = area_spec.ghost_maxrap(); 00089 _safe_rap_for_area = _maxrap_for_area - jet_def.R(); 00090 00091 // Make sure we'll have at least one repetition -- then we can 00092 // deduce the unghosted clustering sequence from one of the ghosted 00093 // sequences. If we do not have any repetitions, then get the 00094 // unghosted sequence from the plain unghosted clustering. 00095 // 00096 // NB: all decanting and filling of initial history will then 00097 // be carried out by base-class routine 00098 if (area_spec.repeat() <= 0) { 00099 _initialise_and_run(jet_def, writeout_combinations); 00100 continue_running = false; 00101 return; 00102 } 00103 00104 // transfer all relevant info into internal variables 00105 _decant_options(jet_def, writeout_combinations); 00106 00107 // set up the history entries for the initial particles (those 00108 // currently in _jets) 00109 _fill_initial_history(); 00110 00111 // by default it does not... 00112 _has_dangerous_particles = false; 00113 00114 continue_running = true; 00115 }
void fastjet::ClusterSequenceActiveArea::_run_AA | ( | const GhostedAreaSpec & | area_spec | ) | [protected] |
Definition at line 119 of file ClusterSequenceActiveArea.cc.
References _has_dangerous_particles, fastjet::ClusterSequence::_jets, _transfer_areas(), _transfer_ghost_free_history(), fastjet::ClusterSequence::jet_def(), fastjet::GhostedAreaSpec::repeat(), and fastjet::ClusterSequence::unique_history_order().
Referenced by _initialise_and_run_AA().
00119 { 00120 // record the input jets as they are currently 00121 vector<PseudoJet> input_jets(_jets); 00122 00123 // code for testing the unique tree 00124 vector<int> unique_tree; 00125 00126 // run the clustering multiple times so as to get areas of all the jets 00127 for (int irepeat = 0; irepeat < area_spec.repeat(); irepeat++) { 00128 00129 ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets, 00130 jet_def(), area_spec); 00131 00132 _has_dangerous_particles |= clust_seq.has_dangerous_particles(); 00133 if (irepeat == 0) { 00134 // take the non-ghost part of the history and put into our own 00135 // history. 00136 _transfer_ghost_free_history(clust_seq); 00137 // get the "unique" order that will be used for transferring all areas. 00138 unique_tree = unique_history_order(); 00139 } 00140 00141 // transfer areas from clust_seq into our object 00142 _transfer_areas(unique_tree, clust_seq); 00143 } 00144 }
void fastjet::ClusterSequenceActiveArea::_postprocess_AA | ( | const GhostedAreaSpec & | area_spec | ) | [protected] |
run the postprocessing for the active area (and derived classes)
Definition at line 149 of file ClusterSequenceActiveArea.cc.
References _average_area, _average_area2, _average_area_4vector, _non_jet_area, _non_jet_area2, _non_jet_number, and fastjet::GhostedAreaSpec::repeat().
Referenced by fastjet::ClusterSequence1GhostPassiveArea::_initialise_and_run_1GPA(), and _initialise_and_run_AA().
00149 { 00150 _average_area /= area_spec.repeat(); 00151 _average_area2 /= area_spec.repeat(); 00152 if (area_spec.repeat() > 1) { 00153 _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/ 00154 (area_spec.repeat()-1)); 00155 } else { 00156 _average_area2 = 0.0; 00157 } 00158 00159 _non_jet_area /= area_spec.repeat(); 00160 _non_jet_area2 /= area_spec.repeat(); 00161 _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/ 00162 area_spec.repeat()); 00163 _non_jet_number /= area_spec.repeat(); 00164 00165 // following bizarre way of writing things is related to 00166 // poverty of operations on PseudoJet objects (as well as some confusion 00167 // in one or two places) 00168 for (unsigned i = 0; i < _average_area_4vector.size(); i++) { 00169 _average_area_4vector[i] = (1.0/area_spec.repeat()) * _average_area_4vector[i]; 00170 } 00171 //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl; 00172 }
void fastjet::ClusterSequenceActiveArea::_initialise_and_run_AA | ( | const JetDefinition & | jet_def, | |
const GhostedAreaSpec & | area_spec, | |||
const bool & | writeout_combinations = false | |||
) | [protected] |
global routine for running active area
Definition at line 50 of file ClusterSequenceActiveArea.cc.
References _initialise_AA(), _postprocess_AA(), _run_AA(), and fastjet::ClusterSequence::jet_def().
Referenced by fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA().
00053 { 00054 00055 bool continue_running; 00056 _initialise_AA(jet_def, area_spec, writeout_combinations, continue_running); 00057 if (continue_running) { 00058 _run_AA(area_spec); 00059 _postprocess_AA(area_spec); 00060 } 00061 }
void fastjet::ClusterSequenceActiveArea::_transfer_ghost_free_history | ( | const ClusterSequenceActiveAreaExplicitGhosts & | clust_seq | ) | [protected] |
transfer the history (and jet-momenta) from clust_seq to our own internal structure while removing ghosts
Definition at line 471 of file ClusterSequenceActiveArea.cc.
References fastjet::ClusterSequence::_do_iB_recombination_step(), fastjet::ClusterSequence::_do_ij_recombination_step(), fastjet::ClusterSequence::_history, fastjet::ClusterSequence::_strategy, fastjet::ClusterSequence::BeamJet, fastjet::ClusterSequence::history(), fastjet::ClusterSequence::InexistentParent, fastjet::ClusterSequence::Invalid, fastjet::ClusterSequenceActiveAreaExplicitGhosts::is_pure_ghost(), and fastjet::ClusterSequence::strategy_used().
Referenced by fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), and _run_AA().
00472 { 00473 00474 const vector<history_element> & gs_history = ghosted_seq.history(); 00475 vector<int> gs2self_hist_map(gs_history.size()); 00476 00477 // work our way through to first non-trivial combination 00478 unsigned igs = 0; 00479 unsigned iself = 0; 00480 while (gs_history[igs].parent1 == InexistentParent) { 00481 // record correspondence 00482 if (!ghosted_seq.is_pure_ghost(igs)) { 00483 gs2self_hist_map[igs] = iself++; 00484 } else { 00485 gs2self_hist_map[igs] = Invalid; 00486 } 00487 igs++; 00488 }; 00489 00490 // make sure the count of non-ghost initial jets is equal to 00491 // what we already have in terms of initial jets 00492 assert(iself == _history.size()); 00493 00494 // now actually transfer things 00495 do { 00496 // if we are a pure ghost, then go on to next round 00497 if (ghosted_seq.is_pure_ghost(igs)) { 00498 gs2self_hist_map[igs] = Invalid; 00499 continue; 00500 } 00501 00502 const history_element & gs_hist_el = gs_history[igs]; 00503 00504 bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1); 00505 bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2); 00506 00507 // if exactly one parent is a ghost then maintain info about the 00508 // non-ghost correspondence for this jet, and then go on to next 00509 // recombination in the ghosted sequence 00510 if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) { 00511 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2]; 00512 continue; 00513 } 00514 if (!parent1_is_ghost && parent2_is_ghost) { 00515 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1]; 00516 continue; 00517 } 00518 00519 // no parents are ghosts... 00520 if (gs_hist_el.parent2 >= 0) { 00521 // recombination of two non-ghosts 00522 gs2self_hist_map[igs] = _history.size(); 00523 // record the recombination in our own sequence 00524 int newjet_k; // dummy var -- not used 00525 int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index; 00526 int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index; 00527 //cerr << "recombining "<< jet_i << " and "<< jet_j << endl; 00528 _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k); 00529 } else { 00530 // we have a non-ghost that has become a beam-jet 00531 assert(gs_history[igs].parent2 == BeamJet); 00532 // record position 00533 gs2self_hist_map[igs] = _history.size(); 00534 // record the recombination in our own sequence 00535 _do_iB_recombination_step( 00536 _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index, 00537 gs_hist_el.dij); 00538 } 00539 } while (++igs < gs_history.size()); 00540 00541 // finally transfer info about strategy used (which isn't necessarily 00542 // always the one that got asked for...) 00543 _strategy = ghosted_seq.strategy_used(); 00544 }
void fastjet::ClusterSequenceActiveArea::_transfer_areas | ( | const vector< int > & | unique_hist_order, | |
const ClusterSequenceActiveAreaExplicitGhosts & | ||||
) | [protected] |
transfer areas from the ClusterSequenceActiveAreaExplicitGhosts object into our internal area bookkeeping.
..
Definition at line 547 of file ClusterSequenceActiveArea.cc.
References _average_area, _average_area2, _average_area_4vector, _ghost_jets, fastjet::ClusterSequence::_history, fastjet::ClusterSequence::_initial_n, fastjet::ClusterSequence::_jet_def, fastjet::ClusterSequence::_jets, _non_jet_area, _non_jet_area2, _non_jet_number, _safe_rap_for_area, _throw_unless_jets_have_same_perp_or_E(), _unclustered_ghosts, fastjet::ClusterSequenceActiveAreaExplicitGhosts::area(), area(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::area_4vector(), fastjet::ClusterSequence::BeamJet, fastjet::ClusterSequence::history(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::is_pure_ghost(), fastjet::ClusterSequence::jets(), fastjet::ClusterSequence::n_particles(), fastjet::JetDefinition::recombiner(), fastjet::ClusterSequence::unclustered_particles(), and fastjet::ClusterSequence::unique_history_order().
Referenced by fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), and _run_AA().
00549 { 00550 00551 const vector<history_element> & gs_history = ghosted_seq.history(); 00552 const vector<PseudoJet> & gs_jets = ghosted_seq.jets(); 00553 vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order(); 00554 00555 const double tolerance = 1e-11; // to decide when two jets are the same 00556 00557 int j = -1; 00558 int hist_index = -1; 00559 00560 valarray<double> our_areas(_history.size()); 00561 our_areas = 0.0; 00562 00563 valarray<PseudoJet> our_area_4vectors(_history.size()); 00564 our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0); 00565 00566 for (unsigned i = 0; i < gs_history.size(); i++) { 00567 // only consider composite particles 00568 unsigned gs_hist_index = gs_unique_hist_order[i]; 00569 if (gs_hist_index < ghosted_seq.n_particles()) continue; 00570 const history_element & gs_hist = gs_history[gs_unique_hist_order[i]]; 00571 int parent1 = gs_hist.parent1; 00572 int parent2 = gs_hist.parent2; 00573 00574 if (parent2 == BeamJet) { 00575 // need to look at parent to get the actual jet 00576 const PseudoJet & jet = 00577 gs_jets[gs_history[parent1].jetp_index]; 00578 double area = ghosted_seq.area(jet); 00579 PseudoJet ext_area = ghosted_seq.area_4vector(jet); 00580 00581 if (ghosted_seq.is_pure_ghost(parent1)) { 00582 // record the existence of the pure ghost jet for future use 00583 _ghost_jets.push_back(GhostJet(jet,area)); 00584 if (abs(jet.rap()) < _safe_rap_for_area) { 00585 _non_jet_area += area; 00586 _non_jet_area2 += area*area; 00587 _non_jet_number += 1; 00588 } 00589 } else { 00590 00591 // get next "combined-particle" index in our own history 00592 // making sure we don't go beyond its bounds (if we do 00593 // then we're in big trouble anyway...) 00594 while (++j < int(_history.size())) { 00595 hist_index = unique_hist_order[j]; 00596 if (hist_index >= _initial_n) break;} 00597 00598 // sanity checking -- do not overrun 00599 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching"); 00600 00601 // sanity check -- make sure we are taking about the same 00602 // jet in reference and new sequences 00603 int refjet_index = _history[_history[hist_index].parent1].jetp_index; 00604 assert(refjet_index >= 0 && refjet_index < int(_jets.size())); 00605 const PseudoJet & refjet = _jets[refjet_index]; 00606 00607 //cerr << "Inclusive" << endl; 00608 //cerr << gs_history[parent1].jetp_index << " " << gs_jets.size() << endl; 00609 //cerr << _history[_history[hist_index].parent1].jetp_index << " " << _jets.size() << endl; 00610 00611 // If pt disagrees check E; if they both disagree there's a 00612 // problem here... NB: a massive particle with zero pt may 00613 // have its pt changed when a ghost is added -- this is why we 00614 // also require the energy to be wrong before complaining 00615 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance, 00616 ghosted_seq); 00617 00618 // set the area at this clustering stage 00619 our_areas[hist_index] = area; 00620 our_area_4vectors[hist_index] = ext_area; 00621 00622 // update the parent as well -- that way its area is the area 00623 // immediately before clustering (i.e. resolve an ambiguity in 00624 // the Cambridge case and ensure in the kt case that the original 00625 // particles get a correct area) 00626 our_areas[_history[hist_index].parent1] = area; 00627 our_area_4vectors[_history[hist_index].parent1] = ext_area; 00628 00629 } 00630 } 00631 else if (!ghosted_seq.is_pure_ghost(parent1) && 00632 !ghosted_seq.is_pure_ghost(parent2)) { 00633 00634 // get next "combined-particle" index in our own history 00635 while (++j < int(_history.size())) { 00636 hist_index = unique_hist_order[j]; 00637 if (hist_index >= _initial_n) break;} 00638 00639 // sanity checking -- do not overrun 00640 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching"); 00641 00642 // make sure that our reference history entry is also for 00643 // an exclusive (dij) clustering (otherwise the comparison jet 00644 // will not exist) 00645 if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)"); 00646 00647 //cerr << "Exclusive: hist_index,hist_size: " << hist_index << " " << _history.size()<< endl; 00648 //cerr << gs_hist.jetp_index << " " << gs_jets.size() << endl; 00649 //cerr << _history[hist_index].jetp_index << " " << _jets.size() << endl; 00650 00651 const PseudoJet & jet = gs_jets[gs_hist.jetp_index]; 00652 const PseudoJet & refjet = _jets[_history[hist_index].jetp_index]; 00653 00654 // run sanity check 00655 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance, 00656 ghosted_seq); 00657 00658 // update area and our local index (maybe redundant since later 00659 // the descendants will reupdate it?) 00660 double area = ghosted_seq.area(jet); 00661 our_areas[hist_index] += area; 00662 00663 PseudoJet ext_area = ghosted_seq.area_4vector(jet); 00664 00665 // GPS TMP debugging (jetclu) ----------------------- 00666 //ext_area = PseudoJet(1e-100,1e-100,1e-100,4e-100); 00667 //our_area_4vectors[hist_index] = ext_area; 00668 //cout << "aa " 00669 // << our_area_4vectors[hist_index].px() << " " 00670 // << our_area_4vectors[hist_index].py() << " " 00671 // << our_area_4vectors[hist_index].pz() << " " 00672 // << our_area_4vectors[hist_index].E() << endl; 00673 //cout << "bb " 00674 // << ext_area.px() << " " 00675 // << ext_area.py() << " " 00676 // << ext_area.pz() << " " 00677 // << ext_area.E() << endl; 00678 //--------------------------------------------------- 00679 00680 _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area); 00681 00682 // now update areas of parents (so that they becomes areas 00683 // immediately before clustering occurred). This is of use 00684 // because it allows us to set the areas of the original hard 00685 // particles in the kt algorithm; for the Cambridge case it 00686 // means a jet's area will be the area just before it clusters 00687 // with another hard jet. 00688 const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index]; 00689 int our_parent1 = _history[hist_index].parent1; 00690 our_areas[our_parent1] = ghosted_seq.area(jet1); 00691 our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1); 00692 00693 const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index]; 00694 int our_parent2 = _history[hist_index].parent2; 00695 our_areas[our_parent2] = ghosted_seq.area(jet2); 00696 our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2); 00697 } 00698 00699 } 00700 00701 // now add unclustered ghosts to the relevant list so that we can 00702 // calculate empty area later. 00703 vector<PseudoJet> unclust = ghosted_seq.unclustered_particles(); 00704 for (unsigned iu = 0; iu < unclust.size(); iu++) { 00705 if (ghosted_seq.is_pure_ghost(unclust[iu])) { 00706 double area = ghosted_seq.area(unclust[iu]); 00707 _unclustered_ghosts.push_back(GhostJet(unclust[iu],area)); 00708 } 00709 } 00710 00711 _average_area += our_areas; 00712 _average_area2 += our_areas*our_areas; 00713 00714 //_average_area_4vector += our_area_4vectors; 00715 // Use the proper recombination scheme when averaging the area_4vectors 00716 // over multiple ghost runs (i.e. the repeat stage); 00717 for (unsigned i = 0; i < _average_area_4vector.size(); i++) { 00718 _jet_def.recombiner()->plus_equal(_average_area_4vector[i], 00719 our_area_4vectors[i]); 00720 } 00721 }
bool fastjet::ClusterSequenceActiveArea::has_dangerous_particles | ( | ) | const [inline, protected] |
returns true if there are any particles whose transverse momentum if so low that there's a risk of the ghosts having modified the clustering sequence
Definition at line 149 of file ClusterSequenceActiveArea.hh.
00149 {return _has_dangerous_particles;}
void fastjet::ClusterSequenceActiveArea::_extract_tree | ( | vector< int > & | ) | const [private] |
routine for extracting the tree in an order that will be independent of any degeneracies in the recombination sequence that don't affect the composition of the final jets
void fastjet::ClusterSequenceActiveArea::_extract_tree_children | ( | int | pos, | |
valarray< bool > & | , | |||
const valarray< int > & | , | |||
vector< int > & | ||||
) | const [private] |
do the part of the extraction associated with pos, working through its children and their parents
void fastjet::ClusterSequenceActiveArea::_extract_tree_parents | ( | int | pos, | |
valarray< bool > & | , | |||
const valarray< int > & | , | |||
vector< int > & | ||||
) | const [private] |
do the part of the extraction associated with the parents of pos.
void fastjet::ClusterSequenceActiveArea::_throw_unless_jets_have_same_perp_or_E | ( | const PseudoJet & | jet, | |
const PseudoJet & | refjet, | |||
double | tolerance, | |||
const ClusterSequenceActiveAreaExplicitGhosts & | jets_ghosted_seq | |||
) | const [private] |
check if two jets have the same momentum to within the tolerance (and if pt's are not the same we're forgiving and look to see if the energy is the same)
Definition at line 727 of file ClusterSequenceActiveArea.cc.
References fastjet::PseudoJet::E(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::has_dangerous_particles(), fastjet::PseudoJet::perp2(), fastjet::PseudoJet::px(), fastjet::PseudoJet::py(), and fastjet::PseudoJet::pz().
Referenced by _transfer_areas().
00732 { 00733 00734 if (abs(jet.perp2()-refjet.perp2()) > 00735 tolerance*max(jet.perp2(),refjet.perp2()) 00736 && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) { 00737 ostringstream ostr; 00738 ostr << "Could not match clustering sequence for an inclusive/exclusive jet when reconstructing areas" << endl; 00739 ostr << " Ref-Jet: " 00740 << refjet.px() << " " 00741 << refjet.py() << " " 00742 << refjet.pz() << " " 00743 << refjet.E() << endl; 00744 ostr << " New-Jet: " 00745 << jet.px() << " " 00746 << jet.py() << " " 00747 << jet.pz() << " " 00748 << jet.E() << endl; 00749 if (jets_ghosted_seq.has_dangerous_particles()) { 00750 ostr << " NB: some particles have pt too low wrt ghosts -- this may be the cause" << endl;} 00751 //ostr << jet.perp() << " " << refjet.perp() << " " 00752 // << jet.perp() - refjet.perp() << endl; 00753 throw Error(ostr.str()); 00754 } 00755 }
valarray<double> fastjet::ClusterSequenceActiveArea::_average_area [protected] |
child classes benefit from having these at their disposal
Definition at line 143 of file ClusterSequenceActiveArea.hh.
Referenced by fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA(), _postprocess_AA(), _resize_and_zero_AA(), fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), and _transfer_areas().
valarray<double> fastjet::ClusterSequenceActiveArea::_average_area2 [protected] |
Definition at line 143 of file ClusterSequenceActiveArea.hh.
Referenced by _postprocess_AA(), _resize_and_zero_AA(), fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), and _transfer_areas().
valarray<PseudoJet> fastjet::ClusterSequenceActiveArea::_average_area_4vector [protected] |
Definition at line 144 of file ClusterSequenceActiveArea.hh.
Referenced by fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA(), _postprocess_AA(), _resize_and_zero_AA(), and _transfer_areas().
double fastjet::ClusterSequenceActiveArea::_non_jet_area [private] |
Definition at line 154 of file ClusterSequenceActiveArea.hh.
Referenced by _postprocess_AA(), _resize_and_zero_AA(), _transfer_areas(), and pt_per_unit_area().
double fastjet::ClusterSequenceActiveArea::_non_jet_area2 [private] |
Definition at line 154 of file ClusterSequenceActiveArea.hh.
Referenced by _postprocess_AA(), _resize_and_zero_AA(), and _transfer_areas().
double fastjet::ClusterSequenceActiveArea::_non_jet_number [private] |
Definition at line 154 of file ClusterSequenceActiveArea.hh.
Referenced by _postprocess_AA(), _resize_and_zero_AA(), _transfer_areas(), and pt_per_unit_area().
double fastjet::ClusterSequenceActiveArea::_maxrap_for_area [private] |
double fastjet::ClusterSequenceActiveArea::_safe_rap_for_area [private] |
Definition at line 157 of file ClusterSequenceActiveArea.hh.
Referenced by _initialise_AA(), _transfer_areas(), and pt_per_unit_area().
bool fastjet::ClusterSequenceActiveArea::_has_dangerous_particles [private] |
Definition at line 159 of file ClusterSequenceActiveArea.hh.
Referenced by _initialise_AA(), and _run_AA().
int fastjet::ClusterSequenceActiveArea::_area_spec_repeat [private] |
since we are playing nasty games with seeds, we should warn the user a few times
Definition at line 187 of file ClusterSequenceActiveArea.hh.
Referenced by _initialise_AA(), empty_area(), and n_empty_jets().
std::vector<GhostJet> fastjet::ClusterSequenceActiveArea::_ghost_jets [private] |
Definition at line 196 of file ClusterSequenceActiveArea.hh.
Referenced by _transfer_areas(), empty_area(), and n_empty_jets().
std::vector<GhostJet> fastjet::ClusterSequenceActiveArea::_unclustered_ghosts [private] |
Definition at line 197 of file ClusterSequenceActiveArea.hh.
Referenced by _transfer_areas(), and empty_area().