#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 std::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 | |
| std::valarray< double > | _average_area |
| child classes benefit from having these at their disposal | |
| std::valarray< double > | _average_area2 |
| child classes benefit from having these at their disposal | |
| std::valarray< PseudoJet > | _average_area_4vector |
Private Member Functions | |
| void | _extract_tree (std::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, std::valarray< bool > &, const std::valarray< int > &, std::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, std::valarray< bool > &, const std::valarray< int > &, std::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.
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};
|
|
|
default constructor
Definition at line 60 of file ClusterSequenceActiveArea.hh. 00060 {}
|
|
||||||||||||||||||||||||
|
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 }
|
|
|
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
|
|
||||||||||||||||||||
|
do the part of the extraction associated with pos, working through its children and their parents
Reimplemented from fastjet::ClusterSequence. |
|
||||||||||||||||||||
|
do the part of the extraction associated with the parents of pos.
Reimplemented from fastjet::ClusterSequence. |
|
||||||||||||||||||||
|
Definition at line 77 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::JetDefinition::R(), and fastjet::GhostedAreaSpec::repeat(). Referenced by fastjet::ClusterSequence1GhostPassiveArea::_initialise_and_run_1GPA(), and _initialise_and_run_AA(). 00082 {
00083
00084 // store this for future use
00085 _area_spec_repeat = area_spec.repeat();
00086
00087 // make sure placeholders are there & zeroed
00088 _resize_and_zero_AA();
00089
00090 // for future reference...
00091 _maxrap_for_area = area_spec.ghost_maxrap();
00092 _safe_rap_for_area = _maxrap_for_area - jet_def.R();
00093
00094 // Make sure we'll have at least one repetition -- then we can
00095 // deduce the unghosted clustering sequence from one of the ghosted
00096 // sequences. If we do not have any repetitions, then get the
00097 // unghosted sequence from the plain unghosted clustering.
00098 //
00099 // NB: all decanting and filling of initial history will then
00100 // be carried out by base-class routine
00101 if (area_spec.repeat() <= 0) {
00102 _initialise_and_run(jet_def, writeout_combinations);
00103 continue_running = false;
00104 return;
00105 }
00106
00107 // transfer all relevant info into internal variables
00108 _decant_options(jet_def, writeout_combinations);
00109
00110 // set up the history entries for the initial particles (those
00111 // currently in _jets)
00112 _fill_initial_history();
00113
00114 // by default it does not...
00115 _has_dangerous_particles = false;
00116
00117 continue_running = true;
00118 }
|
|
||||||||||||||||
|
global routine for running active area
Definition at line 53 of file ClusterSequenceActiveArea.cc. References _initialise_AA(), _postprocess_AA(), and _run_AA(). Referenced by fastjet::ClusterSequencePassiveArea::_initialise_and_run_PA(). 00056 {
00057
00058 bool continue_running;
00059 _initialise_AA(jet_def, area_spec, writeout_combinations, continue_running);
00060 if (continue_running) {
00061 _run_AA(area_spec);
00062 _postprocess_AA(area_spec);
00063 }
00064 }
|
|
|
run the postprocessing for the active area (and derived classes)
Definition at line 152 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(). 00152 {
00153 _average_area /= area_spec.repeat();
00154 _average_area2 /= area_spec.repeat();
00155 if (area_spec.repeat() > 1) {
00156 // the VC compiler complains if one puts everything on a single line.
00157 // An alternative solution would be to use -1.0 (+single line)
00158 const double tmp = area_spec.repeat()-1;
00159 _average_area2 = sqrt(abs(_average_area2 - _average_area*_average_area)/tmp);
00160 } else {
00161 _average_area2 = 0.0;
00162 }
00163
00164 _non_jet_area /= area_spec.repeat();
00165 _non_jet_area2 /= area_spec.repeat();
00166 _non_jet_area2 = sqrt(abs(_non_jet_area2 - _non_jet_area*_non_jet_area)/
00167 area_spec.repeat());
00168 _non_jet_number /= area_spec.repeat();
00169
00170 // following bizarre way of writing things is related to
00171 // poverty of operations on PseudoJet objects (as well as some confusion
00172 // in one or two places)
00173 for (unsigned i = 0; i < _average_area_4vector.size(); i++) {
00174 _average_area_4vector[i] = (1.0/area_spec.repeat()) * _average_area_4vector[i];
00175 }
00176 //cerr << "Non-jet area = " << _non_jet_area << " +- " << _non_jet_area2<<endl;
00177 }
|
|
|
Definition at line 67 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(). 00067 {
00068 // initialize our local area information
00069 _average_area.resize(2*_jets.size()); _average_area = 0.0;
00070 _average_area2.resize(2*_jets.size()); _average_area2 = 0.0;
00071 _average_area_4vector.resize(2*_jets.size());
00072 _average_area_4vector = PseudoJet(0.0,0.0,0.0,0.0);
00073 _non_jet_area = 0.0; _non_jet_area2 = 0.0; _non_jet_number=0.0;
00074 }
|
|
|
Definition at line 122 of file ClusterSequenceActiveArea.cc. References _has_dangerous_particles, fastjet::ClusterSequence::_jets, _transfer_areas(), _transfer_ghost_free_history(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::has_dangerous_particles(), fastjet::ClusterSequence::jet_def(), fastjet::GhostedAreaSpec::repeat(), and fastjet::ClusterSequence::unique_history_order(). Referenced by _initialise_and_run_AA(). 00122 {
00123 // record the input jets as they are currently
00124 vector<PseudoJet> input_jets(_jets);
00125
00126 // code for testing the unique tree
00127 vector<int> unique_tree;
00128
00129 // run the clustering multiple times so as to get areas of all the jets
00130 for (int irepeat = 0; irepeat < area_spec.repeat(); irepeat++) {
00131
00132 ClusterSequenceActiveAreaExplicitGhosts clust_seq(input_jets,
00133 jet_def(), area_spec);
00134
00135 _has_dangerous_particles |= clust_seq.has_dangerous_particles();
00136 if (irepeat == 0) {
00137 // take the non-ghost part of the history and put into our own
00138 // history.
00139 _transfer_ghost_free_history(clust_seq);
00140 // get the "unique" order that will be used for transferring all areas.
00141 unique_tree = unique_history_order();
00142 }
00143
00144 // transfer areas from clust_seq into our object
00145 _transfer_areas(unique_tree, clust_seq);
00146 }
00147 }
|
|
||||||||||||||||||||
|
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 742 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(). 00747 {
00748
00749 if (abs(jet.perp2()-refjet.perp2()) >
00750 tolerance*max(jet.perp2(),refjet.perp2())
00751 && abs(jet.E()-refjet.E()) > tolerance*max(jet.E(),refjet.E())) {
00752 ostringstream ostr;
00753 ostr << "Could not match clustering sequence for an inclusive/exclusive jet when reconstructing areas" << endl;
00754 ostr << " Ref-Jet: "
00755 << refjet.px() << " "
00756 << refjet.py() << " "
00757 << refjet.pz() << " "
00758 << refjet.E() << endl;
00759 ostr << " New-Jet: "
00760 << jet.px() << " "
00761 << jet.py() << " "
00762 << jet.pz() << " "
00763 << jet.E() << endl;
00764 if (jets_ghosted_seq.has_dangerous_particles()) {
00765 ostr << " NB: some particles have pt too low wrt ghosts -- this may be the cause" << endl;}
00766 //ostr << jet.perp() << " " << refjet.perp() << " "
00767 // << jet.perp() - refjet.perp() << endl;
00768 throw Error(ostr.str());
00769 }
00770 }
|
|
||||||||||||
|
transfer areas from the ClusterSequenceActiveAreaExplicitGhosts object into our internal area bookkeeping. .. Definition at line 552 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(). 00554 {
00555
00556 const vector<history_element> & gs_history = ghosted_seq.history();
00557 const vector<PseudoJet> & gs_jets = ghosted_seq.jets();
00558 vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order();
00559
00560 const double tolerance = 1e-11; // to decide when two jets are the same
00561
00562 int j = -1;
00563 int hist_index = -1;
00564
00565 valarray<double> our_areas(_history.size());
00566 our_areas = 0.0;
00567
00568 valarray<PseudoJet> our_area_4vectors(_history.size());
00569 our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0);
00570
00571 for (unsigned i = 0; i < gs_history.size(); i++) {
00572 // only consider composite particles
00573 unsigned gs_hist_index = gs_unique_hist_order[i];
00574 if (gs_hist_index < ghosted_seq.n_particles()) continue;
00575 const history_element & gs_hist = gs_history[gs_unique_hist_order[i]];
00576 int parent1 = gs_hist.parent1;
00577 int parent2 = gs_hist.parent2;
00578
00579 if (parent2 == BeamJet) {
00580 // need to look at parent to get the actual jet
00581 const PseudoJet & jet =
00582 gs_jets[gs_history[parent1].jetp_index];
00583 double area = ghosted_seq.area(jet);
00584 PseudoJet ext_area = ghosted_seq.area_4vector(jet);
00585
00586 if (ghosted_seq.is_pure_ghost(parent1)) {
00587 // record the existence of the pure ghost jet for future use
00588 _ghost_jets.push_back(GhostJet(jet,area));
00589 if (abs(jet.rap()) < _safe_rap_for_area) {
00590 _non_jet_area += area;
00591 _non_jet_area2 += area*area;
00592 _non_jet_number += 1;
00593 }
00594 } else {
00595
00596 // get next "combined-particle" index in our own history
00597 // making sure we don't go beyond its bounds (if we do
00598 // then we're in big trouble anyway...)
00599 while (++j < int(_history.size())) {
00600 hist_index = unique_hist_order[j];
00601 if (hist_index >= _initial_n) break;}
00602
00603 // sanity checking -- do not overrun
00604 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching");
00605
00606 // sanity check -- make sure we are taking about the same
00607 // jet in reference and new sequences
00608 int refjet_index = _history[_history[hist_index].parent1].jetp_index;
00609 assert(refjet_index >= 0 && refjet_index < int(_jets.size()));
00610 const PseudoJet & refjet = _jets[refjet_index];
00611
00612 //cerr << "Inclusive" << endl;
00613 //cerr << gs_history[parent1].jetp_index << " " << gs_jets.size() << endl;
00614 //cerr << _history[_history[hist_index].parent1].jetp_index << " " << _jets.size() << endl;
00615
00616 // If pt disagrees check E; if they both disagree there's a
00617 // problem here... NB: a massive particle with zero pt may
00618 // have its pt changed when a ghost is added -- this is why we
00619 // also require the energy to be wrong before complaining
00620 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
00621 ghosted_seq);
00622
00623 // set the area at this clustering stage
00624 our_areas[hist_index] = area;
00625 our_area_4vectors[hist_index] = ext_area;
00626
00627 // update the parent as well -- that way its area is the area
00628 // immediately before clustering (i.e. resolve an ambiguity in
00629 // the Cambridge case and ensure in the kt case that the original
00630 // particles get a correct area)
00631 our_areas[_history[hist_index].parent1] = area;
00632 our_area_4vectors[_history[hist_index].parent1] = ext_area;
00633
00634 }
00635 }
00636 else if (!ghosted_seq.is_pure_ghost(parent1) &&
00637 !ghosted_seq.is_pure_ghost(parent2)) {
00638
00639 // get next "combined-particle" index in our own history
00640 while (++j < int(_history.size())) {
00641 hist_index = unique_hist_order[j];
00642 if (hist_index >= _initial_n) break;}
00643
00644 // sanity checking -- do not overrun
00645 if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching");
00646
00647 // make sure that our reference history entry is also for
00648 // an exclusive (dij) clustering (otherwise the comparison jet
00649 // will not exist)
00650 if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)");
00651
00652 //cerr << "Exclusive: hist_index,hist_size: " << hist_index << " " << _history.size()<< endl;
00653 //cerr << gs_hist.jetp_index << " " << gs_jets.size() << endl;
00654 //cerr << _history[hist_index].jetp_index << " " << _jets.size() << endl;
00655
00656 const PseudoJet & jet = gs_jets[gs_hist.jetp_index];
00657 const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];
00658
00659 // run sanity check
00660 _throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
00661 ghosted_seq);
00662
00663 // update area and our local index (maybe redundant since later
00664 // the descendants will reupdate it?)
00665 double area = ghosted_seq.area(jet);
00666 our_areas[hist_index] += area;
00667
00668 PseudoJet ext_area = ghosted_seq.area_4vector(jet);
00669
00670 // GPS TMP debugging (jetclu) -----------------------
00671 //ext_area = PseudoJet(1e-100,1e-100,1e-100,4e-100);
00672 //our_area_4vectors[hist_index] = ext_area;
00673 //cout << "aa "
00674 // << our_area_4vectors[hist_index].px() << " "
00675 // << our_area_4vectors[hist_index].py() << " "
00676 // << our_area_4vectors[hist_index].pz() << " "
00677 // << our_area_4vectors[hist_index].E() << endl;
00678 //cout << "bb "
00679 // << ext_area.px() << " "
00680 // << ext_area.py() << " "
00681 // << ext_area.pz() << " "
00682 // << ext_area.E() << endl;
00683 //---------------------------------------------------
00684
00685 _jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area);
00686
00687 // now update areas of parents (so that they becomes areas
00688 // immediately before clustering occurred). This is of use
00689 // because it allows us to set the areas of the original hard
00690 // particles in the kt algorithm; for the Cambridge case it
00691 // means a jet's area will be the area just before it clusters
00692 // with another hard jet.
00693 const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index];
00694 int our_parent1 = _history[hist_index].parent1;
00695 our_areas[our_parent1] = ghosted_seq.area(jet1);
00696 our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1);
00697
00698 const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index];
00699 int our_parent2 = _history[hist_index].parent2;
00700 our_areas[our_parent2] = ghosted_seq.area(jet2);
00701 our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2);
00702 }
00703
00704 }
00705
00706 // now add unclustered ghosts to the relevant list so that we can
00707 // calculate empty area later.
00708 vector<PseudoJet> unclust = ghosted_seq.unclustered_particles();
00709 for (unsigned iu = 0; iu < unclust.size(); iu++) {
00710 if (ghosted_seq.is_pure_ghost(unclust[iu])) {
00711 double area = ghosted_seq.area(unclust[iu]);
00712 _unclustered_ghosts.push_back(GhostJet(unclust[iu],area));
00713 }
00714 }
00715
00716 /*
00717 * WARNING:
00718 * _average_area has explicitly been sized initially to 2*jets().size()
00719 * which can be bigger than our_areas (of size _history.size()
00720 * if there are some unclustered particles.
00721 * So we must take care about boundaries
00722 */
00723
00724 for (unsigned int area_index = 0; area_index<our_areas.size(); area_index++){
00725 _average_area[area_index] += our_areas[area_index];
00726 _average_area2[area_index] += our_areas[area_index]*our_areas[area_index];
00727 }
00728
00729 //_average_area_4vector += our_area_4vectors;
00730 // Use the proper recombination scheme when averaging the area_4vectors
00731 // over multiple ghost runs (i.e. the repeat stage);
00732 for (unsigned i = 0; i < our_area_4vectors.size(); i++) {
00733 _jet_def.recombiner()->plus_equal(_average_area_4vector[i],
00734 our_area_4vectors[i]);
00735 }
00736 }
|
|
|
transfer the history (and jet-momenta) from clust_seq to our own internal structure while removing ghosts
Definition at line 476 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::InexistentParent, and fastjet::ClusterSequence::Invalid. Referenced by fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), and _run_AA(). 00477 {
00478
00479 const vector<history_element> & gs_history = ghosted_seq.history();
00480 vector<int> gs2self_hist_map(gs_history.size());
00481
00482 // work our way through to first non-trivial combination
00483 unsigned igs = 0;
00484 unsigned iself = 0;
00485 while (gs_history[igs].parent1 == InexistentParent) {
00486 // record correspondence
00487 if (!ghosted_seq.is_pure_ghost(igs)) {
00488 gs2self_hist_map[igs] = iself++;
00489 } else {
00490 gs2self_hist_map[igs] = Invalid;
00491 }
00492 igs++;
00493 };
00494
00495 // make sure the count of non-ghost initial jets is equal to
00496 // what we already have in terms of initial jets
00497 assert(iself == _history.size());
00498
00499 // now actually transfer things
00500 do {
00501 // if we are a pure ghost, then go on to next round
00502 if (ghosted_seq.is_pure_ghost(igs)) {
00503 gs2self_hist_map[igs] = Invalid;
00504 continue;
00505 }
00506
00507 const history_element & gs_hist_el = gs_history[igs];
00508
00509 bool parent1_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent1);
00510 bool parent2_is_ghost = ghosted_seq.is_pure_ghost(gs_hist_el.parent2);
00511
00512 // if exactly one parent is a ghost then maintain info about the
00513 // non-ghost correspondence for this jet, and then go on to next
00514 // recombination in the ghosted sequence
00515 if (parent1_is_ghost && !parent2_is_ghost && gs_hist_el.parent2 >= 0) {
00516 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent2];
00517 continue;
00518 }
00519 if (!parent1_is_ghost && parent2_is_ghost) {
00520 gs2self_hist_map[igs] = gs2self_hist_map[gs_hist_el.parent1];
00521 continue;
00522 }
00523
00524 // no parents are ghosts...
00525 if (gs_hist_el.parent2 >= 0) {
00526 // recombination of two non-ghosts
00527 gs2self_hist_map[igs] = _history.size();
00528 // record the recombination in our own sequence
00529 int newjet_k; // dummy var -- not used
00530 int jet_i = _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index;
00531 int jet_j = _history[gs2self_hist_map[gs_hist_el.parent2]].jetp_index;
00532 //cerr << "recombining "<< jet_i << " and "<< jet_j << endl;
00533 _do_ij_recombination_step(jet_i, jet_j, gs_hist_el.dij, newjet_k);
00534 } else {
00535 // we have a non-ghost that has become a beam-jet
00536 assert(gs_history[igs].parent2 == BeamJet);
00537 // record position
00538 gs2self_hist_map[igs] = _history.size();
00539 // record the recombination in our own sequence
00540 _do_iB_recombination_step(
00541 _history[gs2self_hist_map[gs_hist_el.parent1]].jetp_index,
00542 gs_hist_el.dij);
00543 }
00544 } while (++igs < gs_history.size());
00545
00546 // finally transfer info about strategy used (which isn't necessarily
00547 // always the one that got asked for...)
00548 _strategy = ghosted_seq.strategy_used();
00549 }
|
|
|
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(), and pt_per_unit_area(). 00069 {
00070 return _average_area[jet.cluster_hist_index()];};
|
|
|
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()];};
|
|
|
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()];};
|
|
|
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 445 of file ClusterSequenceActiveArea.cc. References _area_spec_repeat, _ghost_jets, _unclustered_ghosts, and fastjet::RangeDefinition::is_in_range(). Referenced by fastjet::ClusterSequencePassiveArea::empty_area(). 00445 {
00446 double empty = 0.0;
00447 // first deal with ghost jets
00448 for (unsigned i = 0; i < _ghost_jets.size(); i++) {
00449 if (range.is_in_range(_ghost_jets[i])) {
00450 empty += _ghost_jets[i].area;
00451 }
00452 }
00453 // then deal with unclustered ghosts
00454 for (unsigned i = 0; i < _unclustered_ghosts.size(); i++) {
00455 if (range.is_in_range(_unclustered_ghosts[i])) {
00456 empty += _unclustered_ghosts[i].area;
00457 }
00458 }
00459 empty /= _area_spec_repeat;
00460 return empty;
00461 }
|
|
|
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;}
|
|
|
return the true number of empty jets (replaces ClusterSequenceAreaBase::n_empty_jets(. ..)) Reimplemented from fastjet::ClusterSequenceAreaBase. Reimplemented in fastjet::ClusterSequence1GhostPassiveArea. Definition at line 464 of file ClusterSequenceActiveArea.cc. References _area_spec_repeat, and _ghost_jets. 00464 {
00465 double inrange = 0;
00466 for (unsigned i = 0; i < _ghost_jets.size(); i++) {
00467 if (range.is_in_range(_ghost_jets[i])) inrange++;
00468 }
00469 inrange /= _area_spec_repeat;
00470 return inrange;
00471 }
|
|
||||||||||||
|
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 276 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, play, pttot_over_areatot, and pttot_over_areatot_cut. 00277 {
00278
00279 vector<PseudoJet> incl_jets = inclusive_jets();
00280 vector<double> pt_over_areas;
00281
00282 for (unsigned i = 0; i < incl_jets.size(); i++) {
00283 if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
00284 double this_area;
00285 if ( strat == median_4vector ) {
00286 this_area = area_4vector(incl_jets[i]).perp();
00287 } else {
00288 this_area = area(incl_jets[i]);
00289 }
00290 pt_over_areas.push_back(incl_jets[i].perp()/this_area);
00291 }
00292 }
00293
00294 // there is nothing inside our region, so answer will always be zero
00295 if (pt_over_areas.size() == 0) {return 0.0;}
00296
00297 // get median (pt/area) [this is the "old" median definition. It considers
00298 // only the "real" jets in calculating the median, i.e. excluding the
00299 // only-ghost ones]
00300 sort(pt_over_areas.begin(), pt_over_areas.end());
00301 double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2];
00302
00303 // new median definition that takes into account non-jet area (i.e.
00304 // jets composed only of ghosts), and for fractional median position
00305 // interpolates between the corresponding entries in the pt_over_areas array
00306 double nj_median_pos = (pt_over_areas.size()-1 - _non_jet_number)/2.0;
00307 double nj_median_ratio;
00308 if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
00309 int int_nj_median = int(nj_median_pos);
00310 nj_median_ratio =
00311 pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
00312 + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
00313 } else {
00314 nj_median_ratio = 0.0;
00315 }
00316
00317
00318 // get various forms of mean (pt/area)
00319 double pt_sum = 0.0, pt_sum_with_cut = 0.0;
00320 double area_sum = _non_jet_area, area_sum_with_cut = _non_jet_area;
00321 double ratio_sum = 0.0;
00322 double ratio_n = _non_jet_number;
00323 for (unsigned i = 0; i < incl_jets.size(); i++) {
00324 if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
00325 double this_area;
00326 if ( strat == median_4vector ) {
00327 this_area = area_4vector(incl_jets[i]).perp();
00328 } else {
00329 this_area = area(incl_jets[i]);
00330 }
00331 pt_sum += incl_jets[i].perp();
00332 area_sum += this_area;
00333 double ratio = incl_jets[i].perp()/this_area;
00334 if (ratio < range*nj_median_ratio) {
00335 pt_sum_with_cut += incl_jets[i].perp();
00336 area_sum_with_cut += this_area;
00337 ratio_sum += ratio; ratio_n++;
00338 }
00339 }
00340 }
00341
00342 if (strat == play) {
00343 double trunc_sum = 0, trunc_sumsqr = 0;
00344 vector<double> means(pt_over_areas.size()), sd(pt_over_areas.size());
00345 for (unsigned i = 0; i < pt_over_areas.size() ; i++ ) {
00346 double ratio = pt_over_areas[i];
00347 trunc_sum += ratio;
00348 trunc_sumsqr += ratio*ratio;
00349 means[i] = trunc_sum / (i+1);
00350 sd[i] = sqrt(abs(means[i]*means[i] - trunc_sumsqr/(i+1)));
00351 cerr << "i, means, sd: " <<i<<", "<< means[i] <<", "<<sd[i]<<", "<<
00352 sd[i]/sqrt(i+1.0)<<endl;
00353 }
00354 cout << "-----------------------------------"<<endl;
00355 for (unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) {
00356 cout << "Median "<< i <<" = " << pt_over_areas[i]<<endl;
00357 }
00358 cout << "Number of non-jets: "<<_non_jet_number<<endl;
00359 cout << "Area of non-jets: "<<_non_jet_area<<endl;
00360 cout << "Default median position: " << (pt_over_areas.size()-1)/2.0<<endl;
00361 cout << "NJ median position: " << nj_median_pos <<endl;
00362 cout << "NJ median value: " << nj_median_ratio <<endl;
00363 return 0.0;
00364 }
00365
00366 switch(strat) {
00367 case median:
00368 case median_4vector:
00369 return nj_median_ratio;
00370 case non_ghost_median:
00371 return non_ghost_median_ratio;
00372 case pttot_over_areatot:
00373 return pt_sum / area_sum;
00374 case pttot_over_areatot_cut:
00375 return pt_sum_with_cut / area_sum_with_cut;
00376 case mean_ratio_cut:
00377 return ratio_sum/ratio_n;
00378 default:
00379 return nj_median_ratio;
00380 }
00381
00382 }
|
|
|
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(). |
|
|
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(). |
|
|
child classes benefit from having these at their disposal
Definition at line 143 of file ClusterSequenceActiveArea.hh. Referenced by _postprocess_AA(), _resize_and_zero_AA(), fastjet::ClusterSequence1GhostPassiveArea::_run_1GPA(), and _transfer_areas(). |
|
|
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(). |
|
|
Definition at line 196 of file ClusterSequenceActiveArea.hh. Referenced by _transfer_areas(), empty_area(), and n_empty_jets(). |
|
|
Definition at line 159 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_AA(), and _run_AA(). |
|
|
Definition at line 156 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_AA(). |
|
|
Definition at line 154 of file ClusterSequenceActiveArea.hh. Referenced by _postprocess_AA(), _resize_and_zero_AA(), _transfer_areas(), and pt_per_unit_area(). |
|
|
Definition at line 154 of file ClusterSequenceActiveArea.hh. Referenced by _postprocess_AA(), _resize_and_zero_AA(), and _transfer_areas(). |
|
|
Definition at line 154 of file ClusterSequenceActiveArea.hh. Referenced by _postprocess_AA(), _resize_and_zero_AA(), _transfer_areas(), and pt_per_unit_area(). |
|
|
Definition at line 157 of file ClusterSequenceActiveArea.hh. Referenced by _initialise_AA(), _transfer_areas(), and pt_per_unit_area(). |
|
|
Definition at line 197 of file ClusterSequenceActiveArea.hh. Referenced by _transfer_areas(), and empty_area(). |
1.4.2