|
FastJet 3.0alpha3
|
|
FastJet 3.0alpha3
|
Like ClusterSequence with computation of the active jet area. More...
#include <ClusterSequenceActiveArea.hh>
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 &ghost_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 Selector &selector) const |
| rewrite the empty area from the parent class, so as to use all info at our disposal return the total area, corresponding to a given Selector, that consists of ghost jets or unclustered ghosts | |
| virtual double | n_empty_jets (const Selector &selector) 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 &ghost_spec, const bool &writeout_combinations, bool &continue_running) |
| void | _run_AA (const GhostedAreaSpec &ghost_spec) |
| void | _postprocess_AA (const GhostedAreaSpec &ghost_spec) |
| run the postprocessing for the active area (and derived classes) | |
| void | _initialise_and_run_AA (const JetDefinition &jet_def, const GhostedAreaSpec &ghost_spec, const bool &writeout_combinations=false) |
| does the initialisation and running specific to the active areas class | |
| 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 |
| std::valarray< PseudoJet > | _average_area_4vector |
Like ClusterSequence with computation of the active jet area.
Class that behaves essentially like ClusterSequence except that it also provides access to the area of a jet (which will be a random quantity... Figure out what to do about seeds later...)
This class should not be used directly. Rather use ClusterSequenceArea with the appropriate AreaDefinition
Definition at line 63 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 87 of file ClusterSequenceActiveArea.hh.
{median=0, non_ghost_median, pttot_over_areatot,
pttot_over_areatot_cut, mean_ratio_cut, play,
median_4vector};
| virtual double fastjet::ClusterSequenceActiveArea::area | ( | const PseudoJet & | ) | const [inline, virtual] |
return the area associated with the given jet; this base class returns 0.
Reimplemented from fastjet::ClusterSequenceAreaBase.
Definition at line 76 of file ClusterSequenceActiveArea.hh.
References fastjet::PseudoJet::cluster_hist_index().
{
return _average_area[jet.cluster_hist_index()];};
| virtual double fastjet::ClusterSequenceActiveArea::area_error | ( | const PseudoJet & | ) | 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 78 of file ClusterSequenceActiveArea.hh.
References fastjet::PseudoJet::cluster_hist_index().
{
return _average_area2[jet.cluster_hist_index()];};
| virtual PseudoJet fastjet::ClusterSequenceActiveArea::area_4vector | ( | const PseudoJet & | ) | 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 81 of file ClusterSequenceActiveArea.hh.
References fastjet::PseudoJet::cluster_hist_index().
{
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 276 of file ClusterSequenceActiveArea.cc.
{
vector<PseudoJet> incl_jets = inclusive_jets();
vector<double> pt_over_areas;
for (unsigned i = 0; i < incl_jets.size(); i++) {
if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
double this_area;
if ( strat == median_4vector ) {
this_area = area_4vector(incl_jets[i]).perp();
} else {
this_area = area(incl_jets[i]);
}
pt_over_areas.push_back(incl_jets[i].perp()/this_area);
}
}
// there is nothing inside our region, so answer will always be zero
if (pt_over_areas.size() == 0) {return 0.0;}
// get median (pt/area) [this is the "old" median definition. It considers
// only the "real" jets in calculating the median, i.e. excluding the
// only-ghost ones]
sort(pt_over_areas.begin(), pt_over_areas.end());
double non_ghost_median_ratio = pt_over_areas[pt_over_areas.size()/2];
// new median definition that takes into account non-jet area (i.e.
// jets composed only of ghosts), and for fractional median position
// interpolates between the corresponding entries in the pt_over_areas array
double nj_median_pos = (pt_over_areas.size()-1 - _non_jet_number)/2.0;
double nj_median_ratio;
if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
int int_nj_median = int(nj_median_pos);
nj_median_ratio =
pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
+ pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
} else {
nj_median_ratio = 0.0;
}
// get various forms of mean (pt/area)
double pt_sum = 0.0, pt_sum_with_cut = 0.0;
double area_sum = _non_jet_area, area_sum_with_cut = _non_jet_area;
double ratio_sum = 0.0;
double ratio_n = _non_jet_number;
for (unsigned i = 0; i < incl_jets.size(); i++) {
if (abs(incl_jets[i].rap()) < _safe_rap_for_area) {
double this_area;
if ( strat == median_4vector ) {
this_area = area_4vector(incl_jets[i]).perp();
} else {
this_area = area(incl_jets[i]);
}
pt_sum += incl_jets[i].perp();
area_sum += this_area;
double ratio = incl_jets[i].perp()/this_area;
if (ratio < range*nj_median_ratio) {
pt_sum_with_cut += incl_jets[i].perp();
area_sum_with_cut += this_area;
ratio_sum += ratio; ratio_n++;
}
}
}
if (strat == play) {
double trunc_sum = 0, trunc_sumsqr = 0;
vector<double> means(pt_over_areas.size()), sd(pt_over_areas.size());
for (unsigned i = 0; i < pt_over_areas.size() ; i++ ) {
double ratio = pt_over_areas[i];
trunc_sum += ratio;
trunc_sumsqr += ratio*ratio;
means[i] = trunc_sum / (i+1);
sd[i] = sqrt(abs(means[i]*means[i] - trunc_sumsqr/(i+1)));
cerr << "i, means, sd: " <<i<<", "<< means[i] <<", "<<sd[i]<<", "<<
sd[i]/sqrt(i+1.0)<<endl;
}
cout << "-----------------------------------"<<endl;
for (unsigned i = 0; i <= pt_over_areas.size()/2 ; i++ ) {
cout << "Median "<< i <<" = " << pt_over_areas[i]<<endl;
}
cout << "Number of non-jets: "<<_non_jet_number<<endl;
cout << "Area of non-jets: "<<_non_jet_area<<endl;
cout << "Default median position: " << (pt_over_areas.size()-1)/2.0<<endl;
cout << "NJ median position: " << nj_median_pos <<endl;
cout << "NJ median value: " << nj_median_ratio <<endl;
return 0.0;
}
switch(strat) {
case median:
case median_4vector:
return nj_median_ratio;
case non_ghost_median:
return non_ghost_median_ratio;
case pttot_over_areatot:
return pt_sum / area_sum;
case pttot_over_areatot_cut:
return pt_sum_with_cut / area_sum_with_cut;
case mean_ratio_cut:
return ratio_sum/ratio_n;
default:
return nj_median_ratio;
}
}
| double fastjet::ClusterSequenceActiveArea::empty_area | ( | const Selector & | selector | ) | const [virtual] |
rewrite the empty area from the parent class, so as to use all info at our disposal return the total area, corresponding to a given Selector, that consists of ghost jets or unclustered ghosts
The selector passed as an argument needs to apply jet by jet.
Reimplemented from fastjet::ClusterSequenceAreaBase.
Reimplemented in fastjet::ClusterSequencePassiveArea.
Definition at line 445 of file ClusterSequenceActiveArea.cc.
References fastjet::Selector::applies_jet_by_jet(), and fastjet::Selector::pass().
{
// make sure that the selector applies jet by jet
if (! selector.applies_jet_by_jet()){
throw Error("ClusterSequenceActiveArea: empty area can only be computed from selectors applying jet by jet");
}
double empty = 0.0;
// first deal with ghost jets
for (unsigned i = 0; i < _ghost_jets.size(); i++) {
if (selector.pass(_ghost_jets[i])) {
empty += _ghost_jets[i].area;
}
}
// then deal with unclustered ghosts
for (unsigned i = 0; i < _unclustered_ghosts.size(); i++) {
if (selector.pass(_unclustered_ghosts[i])) {
empty += _unclustered_ghosts[i].area;
}
}
empty /= _ghost_spec_repeat;
return empty;
}
| void fastjet::ClusterSequenceActiveArea::_initialise_and_run_AA | ( | const JetDefinition & | jet_def, |
| const GhostedAreaSpec & | ghost_spec, | ||
| const bool & | writeout_combinations = false |
||
| ) | [protected] |
does the initialisation and running specific to the active areas class
global routine for running active area
Definition at line 53 of file ClusterSequenceActiveArea.cc.
{
bool continue_running;
_initialise_AA(jet_def, ghost_spec, writeout_combinations, continue_running);
if (continue_running) {
_run_AA(ghost_spec);
_postprocess_AA(ghost_spec);
}
}
| void fastjet::ClusterSequenceActiveArea::_transfer_areas | ( | const std::vector< int > & | unique_hist_order, |
| const ClusterSequenceActiveAreaExplicitGhosts & | |||
| ) | [protected] |
transfer areas from the ClusterSequenceActiveAreaExplicitGhosts object into our internal area bookkeeping...
Definition at line 566 of file ClusterSequenceActiveArea.cc.
References fastjet::ClusterSequenceActiveAreaExplicitGhosts::area(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::area_4vector(), fastjet::ClusterSequence::history(), fastjet::ClusterSequenceActiveAreaExplicitGhosts::is_pure_ghost(), fastjet::ClusterSequence::history_element::jetp_index, fastjet::ClusterSequence::jets(), fastjet::ClusterSequence::n_particles(), fastjet::ClusterSequence::history_element::parent2, fastjet::PseudoJet::rap(), fastjet::ClusterSequence::unclustered_particles(), and fastjet::ClusterSequence::unique_history_order().
{
const vector<history_element> & gs_history = ghosted_seq.history();
const vector<PseudoJet> & gs_jets = ghosted_seq.jets();
vector<int> gs_unique_hist_order = ghosted_seq.unique_history_order();
const double tolerance = 1e-11; // to decide when two jets are the same
int j = -1;
int hist_index = -1;
valarray<double> our_areas(_history.size());
our_areas = 0.0;
valarray<PseudoJet> our_area_4vectors(_history.size());
our_area_4vectors = PseudoJet(0.0,0.0,0.0,0.0);
for (unsigned i = 0; i < gs_history.size(); i++) {
// only consider composite particles
unsigned gs_hist_index = gs_unique_hist_order[i];
if (gs_hist_index < ghosted_seq.n_particles()) continue;
const history_element & gs_hist = gs_history[gs_unique_hist_order[i]];
int parent1 = gs_hist.parent1;
int parent2 = gs_hist.parent2;
if (parent2 == BeamJet) {
// need to look at parent to get the actual jet
const PseudoJet & jet =
gs_jets[gs_history[parent1].jetp_index];
double area = ghosted_seq.area(jet);
PseudoJet ext_area = ghosted_seq.area_4vector(jet);
if (ghosted_seq.is_pure_ghost(parent1)) {
// record the existence of the pure ghost jet for future use
_ghost_jets.push_back(GhostJet(jet,area));
if (abs(jet.rap()) < _safe_rap_for_area) {
_non_jet_area += area;
_non_jet_area2 += area*area;
_non_jet_number += 1;
}
} else {
// get next "combined-particle" index in our own history
// making sure we don't go beyond its bounds (if we do
// then we're in big trouble anyway...)
while (++j < int(_history.size())) {
hist_index = unique_hist_order[j];
if (hist_index >= _initial_n) break;}
// sanity checking -- do not overrun
if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in diB matching");
// sanity check -- make sure we are taking about the same
// jet in reference and new sequences
int refjet_index = _history[_history[hist_index].parent1].jetp_index;
assert(refjet_index >= 0 && refjet_index < int(_jets.size()));
const PseudoJet & refjet = _jets[refjet_index];
//cerr << "Inclusive" << endl;
//cerr << gs_history[parent1].jetp_index << " " << gs_jets.size() << endl;
//cerr << _history[_history[hist_index].parent1].jetp_index << " " << _jets.size() << endl;
// If pt disagrees check E; if they both disagree there's a
// problem here... NB: a massive particle with zero pt may
// have its pt changed when a ghost is added -- this is why we
// also require the energy to be wrong before complaining
_throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
ghosted_seq);
// set the area at this clustering stage
our_areas[hist_index] = area;
our_area_4vectors[hist_index] = ext_area;
// update the parent as well -- that way its area is the area
// immediately before clustering (i.e. resolve an ambiguity in
// the Cambridge case and ensure in the kt case that the original
// particles get a correct area)
our_areas[_history[hist_index].parent1] = area;
our_area_4vectors[_history[hist_index].parent1] = ext_area;
}
}
else if (!ghosted_seq.is_pure_ghost(parent1) &&
!ghosted_seq.is_pure_ghost(parent2)) {
// get next "combined-particle" index in our own history
while (++j < int(_history.size())) {
hist_index = unique_hist_order[j];
if (hist_index >= _initial_n) break;}
// sanity checking -- do not overrun
if (j >= int(_history.size())) throw Error("ClusterSequenceActiveArea: overran reference array in dij matching");
// make sure that our reference history entry is also for
// an exclusive (dij) clustering (otherwise the comparison jet
// will not exist)
if (_history[hist_index].parent2 == BeamJet) throw Error("ClusterSequenceActiveArea: could not match clustering sequences (encountered dij matched with diB)");
//cerr << "Exclusive: hist_index,hist_size: " << hist_index << " " << _history.size()<< endl;
//cerr << gs_hist.jetp_index << " " << gs_jets.size() << endl;
//cerr << _history[hist_index].jetp_index << " " << _jets.size() << endl;
const PseudoJet & jet = gs_jets[gs_hist.jetp_index];
const PseudoJet & refjet = _jets[_history[hist_index].jetp_index];
// run sanity check
_throw_unless_jets_have_same_perp_or_E(jet, refjet, tolerance,
ghosted_seq);
// update area and our local index (maybe redundant since later
// the descendants will reupdate it?)
double area = ghosted_seq.area(jet);
our_areas[hist_index] += area;
PseudoJet ext_area = ghosted_seq.area_4vector(jet);
// GPS TMP debugging (jetclu) -----------------------
//ext_area = PseudoJet(1e-100,1e-100,1e-100,4e-100);
//our_area_4vectors[hist_index] = ext_area;
//cout << "aa "
// << our_area_4vectors[hist_index].px() << " "
// << our_area_4vectors[hist_index].py() << " "
// << our_area_4vectors[hist_index].pz() << " "
// << our_area_4vectors[hist_index].E() << endl;
//cout << "bb "
// << ext_area.px() << " "
// << ext_area.py() << " "
// << ext_area.pz() << " "
// << ext_area.E() << endl;
//---------------------------------------------------
_jet_def.recombiner()->plus_equal(our_area_4vectors[hist_index], ext_area);
// now update areas of parents (so that they becomes areas
// immediately before clustering occurred). This is of use
// because it allows us to set the areas of the original hard
// particles in the kt algorithm; for the Cambridge case it
// means a jet's area will be the area just before it clusters
// with another hard jet.
const PseudoJet & jet1 = gs_jets[gs_history[parent1].jetp_index];
int our_parent1 = _history[hist_index].parent1;
our_areas[our_parent1] = ghosted_seq.area(jet1);
our_area_4vectors[our_parent1] = ghosted_seq.area_4vector(jet1);
const PseudoJet & jet2 = gs_jets[gs_history[parent2].jetp_index];
int our_parent2 = _history[hist_index].parent2;
our_areas[our_parent2] = ghosted_seq.area(jet2);
our_area_4vectors[our_parent2] = ghosted_seq.area_4vector(jet2);
}
}
// now add unclustered ghosts to the relevant list so that we can
// calculate empty area later.
vector<PseudoJet> unclust = ghosted_seq.unclustered_particles();
for (unsigned iu = 0; iu < unclust.size(); iu++) {
if (ghosted_seq.is_pure_ghost(unclust[iu])) {
double area = ghosted_seq.area(unclust[iu]);
_unclustered_ghosts.push_back(GhostJet(unclust[iu],area));
}
}
/*
* WARNING:
* _average_area has explicitly been sized initially to 2*jets().size()
* which can be bigger than our_areas (of size _history.size()
* if there are some unclustered particles.
* So we must take care about boundaries
*/
for (unsigned int area_index = 0; area_index<our_areas.size(); area_index++){
_average_area[area_index] += our_areas[area_index];
_average_area2[area_index] += our_areas[area_index]*our_areas[area_index];
}
//_average_area_4vector += our_area_4vectors;
// Use the proper recombination scheme when averaging the area_4vectors
// over multiple ghost runs (i.e. the repeat stage);
for (unsigned i = 0; i < our_area_4vectors.size(); i++) {
_jet_def.recombiner()->plus_equal(_average_area_4vector[i],
our_area_4vectors[i]);
}
}
1.7.4