35#include "fastjet/ClusterSequenceAreaBase.hh"
38FASTJET_BEGIN_NAMESPACE
44LimitedWarning ClusterSequenceAreaBase::_warnings;
45LimitedWarning ClusterSequenceAreaBase::_warnings_zero_area;
46LimitedWarning ClusterSequenceAreaBase::_warnings_empty_area;
57double ClusterSequenceAreaBase::empty_area(
const Selector & selector)
const {
59 if (has_explicit_ghosts()) {
return 0.0;}
60 else {
return empty_area_from_jets(inclusive_jets(0.0), selector);}
69double ClusterSequenceAreaBase::empty_area_from_jets(
70 const std::vector<PseudoJet> & all_jets,
72 _check_selector_good_for_median(selector);
74 double empty = selector.
area();
75 for (
unsigned i = 0; i < all_jets.size(); i++) {
76 if (selector.
pass(all_jets[i])) empty -= area(all_jets[i]);
86double ClusterSequenceAreaBase::median_pt_per_unit_area(
const Selector & selector)
const {
87 return _median_pt_per_unit_area(selector);
91double ClusterSequenceAreaBase::_median_pt_per_unit_area(
const Selector & selector)
const {
92 return _median_pt_per_unit_something(selector,
false);
102double ClusterSequenceAreaBase::median_pt_per_unit_area_4vector(
const Selector & selector)
const {
103 return _median_pt_per_unit_area_4vector(selector);
107double ClusterSequenceAreaBase::_median_pt_per_unit_area_4vector(
const Selector & selector)
const {
108 return _median_pt_per_unit_something(selector,
true);
118double ClusterSequenceAreaBase::median_pt_per_unit_something(
119 const Selector & selector,
bool use_area_4vector)
const {
120 return _median_pt_per_unit_something(selector, use_area_4vector);
126double ClusterSequenceAreaBase::_median_pt_per_unit_something(
127 const Selector & selector,
bool use_area_4vector)
const {
128 double median, sigma, mean_area;
129 _get_median_rho_and_sigma(selector, use_area_4vector, median, sigma, mean_area);
138void ClusterSequenceAreaBase::parabolic_pt_per_unit_area(
139 double & a,
double & b,
const Selector & selector,
140 double exclude_above,
bool use_area_4vector)
const {
141 return _parabolic_pt_per_unit_area(a, b, selector, exclude_above, use_area_4vector);
144void ClusterSequenceAreaBase::_parabolic_pt_per_unit_area(
145 double & a,
double & b,
const Selector & selector,
146 double exclude_above,
bool use_area_4vector)
const {
150 _check_selector_good_for_median(selector);
153 double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0;
155 vector<PseudoJet> incl_jets = inclusive_jets();
157 for (
unsigned i = 0; i < incl_jets.size(); i++) {
158 if (selector.
pass(incl_jets[i])) {
160 if ( use_area_4vector ) {
161 this_area = area_4vector(incl_jets[i]).perp();
163 this_area = area(incl_jets[i]);
165 double f = incl_jets[i].perp()/this_area;
166 if (exclude_above <= 0.0 || f < exclude_above) {
167 double x = incl_jets[i].rap();
double x2 = x*x;
189 b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
190 a = mean_f - b*mean_x2;
196void ClusterSequenceAreaBase::get_median_rho_and_sigma(
197 const Selector & selector,
bool use_area_4vector,
198 double & median,
double & sigma,
double & mean_area)
const {
199 _get_median_rho_and_sigma(selector, use_area_4vector, median, sigma, mean_area);
202void ClusterSequenceAreaBase::_get_median_rho_and_sigma(
203 const Selector & selector,
bool use_area_4vector,
204 double & median,
double & sigma,
double & mean_area)
const {
206 vector<PseudoJet> incl_jets = inclusive_jets();
207 _get_median_rho_and_sigma(incl_jets, selector, use_area_4vector,
208 median, sigma, mean_area,
true);
211void ClusterSequenceAreaBase::get_median_rho_and_sigma(
212 const vector<PseudoJet> & all_jets,
213 const Selector & selector,
bool use_area_4vector,
214 double & median,
double & sigma,
double & mean_area,
215 bool all_are_incl)
const {
216 _get_median_rho_and_sigma(all_jets, selector, use_area_4vector,
217 median, sigma, mean_area, all_are_incl);
220void ClusterSequenceAreaBase::_get_median_rho_and_sigma(
221 const vector<PseudoJet> & all_jets,
222 const Selector & selector,
bool use_area_4vector,
223 double & median,
double & sigma,
double & mean_area,
224 bool all_are_incl)
const {
226 _check_jet_alg_good_for_median();
231 _check_selector_good_for_median(selector);
233 vector<double> pt_over_areas;
234 double total_area = 0.0;
235 double total_njets = 0;
237 for (
unsigned i = 0; i < all_jets.size(); i++) {
238 if (selector.pass(all_jets[i])) {
240 if (use_area_4vector) {
241 this_area = area_4vector(all_jets[i]).perp();
243 this_area = area(all_jets[i]);
247 pt_over_areas.push_back(all_jets[i].perp()/this_area);
249 _warnings_zero_area.warn(
"ClusterSequenceAreaBase::get_median_rho_and_sigma(...): discarded jet with zero area. Zero-area jets may be due to (i) too large a ghost area (ii) a jet being outside the ghost range (iii) the computation not being done using an appropriate algorithm (kt;C/A).");
252 total_area += this_area;
258 if (pt_over_areas.size() == 0) {
268 sort(pt_over_areas.begin(), pt_over_areas.end());
272 double posn[2] = {0.5, (1.0-0.6827)/2.0};
275 double n_empty, empty_a;
276 if (has_explicit_ghosts()) {
284 }
else if (all_are_incl) {
286 empty_a = empty_area(selector);
287 n_empty = n_empty_jets(selector);
292 empty_a = empty_area_from_jets(all_jets, selector);
293 mean_area = total_area / total_njets;
294 n_empty = empty_a / mean_area;
298 total_njets += n_empty;
299 total_area += empty_a;
304 int pt_over_areas_size = pt_over_areas.size();
305 if (n_empty < -pt_over_areas_size/4.0) {
306 _warnings_empty_area.warn(
"ClusterSequenceAreaBase::get_median_rho_and_sigma(...): the estimated empty area is suspiciously large and negative and may lead to an over-estimation of rho. This may be due to (i) a rare statistical fluctuation or (ii) too small a range used to estimate the background properties.");
309 for (
int i = 0; i < 2; i++) {
310 double nj_median_pos =
311 (pt_over_areas_size-1.0 + n_empty)*posn[i] - n_empty;
312 double nj_median_ratio;
313 if (nj_median_pos >= 0 && pt_over_areas_size > 1) {
314 int int_nj_median = int(nj_median_pos);
317 if (int_nj_median+1 > pt_over_areas_size-1){
318 int_nj_median = pt_over_areas_size-2;
319 nj_median_pos = pt_over_areas_size-1;
323 pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
324 + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
326 nj_median_ratio = 0.0;
328 res[i] = nj_median_ratio;
331 double error = res[0] - res[1];
332 mean_area = total_area / total_njets;
333 sigma = error * sqrt( max(0.0, mean_area) );
341vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
const double rho,
344 return _subtracted_jets(rho,ptmin);
347vector<PseudoJet> ClusterSequenceAreaBase::_subtracted_jets(
const double rho,
350 vector<PseudoJet> sub_jets;
351 vector<PseudoJet> jets_local =
sorted_by_pt(inclusive_jets(ptmin));
352 for (
unsigned i=0; i<jets_local.size(); i++) {
353 PseudoJet sub_jet = _subtracted_jet(jets_local[i],rho);
354 sub_jets.push_back(sub_jet);
363vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
367 double rho = _median_pt_per_unit_area_4vector(selector);
368 return _subtracted_jets(rho,ptmin);
374 const double rho)
const {
375 return _subtracted_jet(jet, rho);
379 const double rho)
const {
383 if (rho*area4vect.
perp() < jet.
perp() ) {
384 sub_jet = jet - rho*area4vect;
385 }
else { sub_jet = PseudoJet(0.0,0.0,0.0,0.0); }
403 return _subtracted_jet(jet, selector);
408 double rho = _median_pt_per_unit_area_4vector(selector);
409 PseudoJet sub_jet = _subtracted_jet(jet, rho);
415double ClusterSequenceAreaBase::subtracted_pt(
const PseudoJet & jet,
417 bool use_area_4vector)
const {
418 return _subtracted_pt(jet, rho, use_area_4vector);
421double ClusterSequenceAreaBase::_subtracted_pt(
const PseudoJet & jet,
423 bool use_area_4vector)
const {
424 if ( use_area_4vector ) {
425 PseudoJet sub_jet = _subtracted_jet(jet,rho);
426 return sub_jet.
perp();
428 return jet.
perp() - rho*area(jet);
436double ClusterSequenceAreaBase::subtracted_pt(
const PseudoJet & jet,
438 bool use_area_4vector)
const {
439 if ( use_area_4vector ) {
440 PseudoJet sub_jet = _subtracted_jet(jet,selector);
441 return sub_jet.
perp();
443 double rho = _median_pt_per_unit_area(selector);
444 return _subtracted_pt(jet,rho,
false);
450void ClusterSequenceAreaBase::_check_selector_good_for_median(
const Selector &selector)
const{
453 throw Error(
"ClusterSequenceAreaBase: empty area can only be computed from selectors with a finite area");
458 throw Error(
"ClusterSequenceAreaBase: empty area can only be computed from selectors that apply jet by jet");
464void ClusterSequenceAreaBase::_check_jet_alg_good_for_median()
const {
465 if (jet_def().jet_algorithm() != kt_algorithm
466 && jet_def().jet_algorithm() != cambridge_algorithm
467 && jet_def().jet_algorithm() != cambridge_for_passive_algorithm) {
468 _warnings.warn(
"ClusterSequenceAreaBase: jet_def being used may not be suitable for estimating diffuse backgrounds (good options are kt, cam)");
base class corresponding to errors that can be thrown by FastJet
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
void set_user_index(const int index)
set the user_index, intended to allow the user to add simple identifying information to a particle/je...
double perp() const
returns the scalar transverse momentum
void set_cluster_hist_index(const int index)
set the cluster_hist_index, intended to be used by clustering routines.
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure_in)
set the associated structure
int user_index() const
return the user_index,
const SharedPtr< PseudoJetStructureBase > & structure_shared_ptr() const
return a reference to the shared pointer to the PseudoJetStructureBase associated with this PseudoJet
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
bool has_finite_area() const
returns true if it has a meaningful and finite area (i.e.
bool pass(const PseudoJet &jet) const
return true if the jet passes the selection
bool applies_jet_by_jet() const
returns true if this can be applied jet by jet
double area() const
returns the rapidity-phi area associated with the Selector (throws InvalidArea if the area does not m...
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2