33 #include "fastjet/ClusterSequenceAreaBase.hh"
36 FASTJET_BEGIN_NAMESPACE
42 LimitedWarning ClusterSequenceAreaBase::_warnings;
43 LimitedWarning ClusterSequenceAreaBase::_warnings_zero_area;
44 LimitedWarning ClusterSequenceAreaBase::_warnings_empty_area;
55 double ClusterSequenceAreaBase::empty_area(
const Selector & selector)
const {
57 if (has_explicit_ghosts()) {
return 0.0;}
58 else {
return empty_area_from_jets(inclusive_jets(0.0), selector);}
67 double ClusterSequenceAreaBase::empty_area_from_jets(
68 const std::vector<PseudoJet> & all_jets,
70 _check_selector_good_for_median(selector);
72 double empty = selector.
area();
73 for (
unsigned i = 0; i < all_jets.size(); i++) {
74 if (selector.
pass(all_jets[i])) empty -= area(all_jets[i]);
79 double ClusterSequenceAreaBase::median_pt_per_unit_area(
const Selector & selector)
const {
80 return median_pt_per_unit_something(selector,
false);
83 double ClusterSequenceAreaBase::median_pt_per_unit_area_4vector(
const Selector & selector)
const {
84 return median_pt_per_unit_something(selector,
true);
92 double ClusterSequenceAreaBase::median_pt_per_unit_something(
93 const Selector & selector,
bool use_area_4vector)
const {
95 double median, sigma, mean_area;
96 get_median_rho_and_sigma(selector, use_area_4vector, median, sigma, mean_area);
106 void ClusterSequenceAreaBase::parabolic_pt_per_unit_area(
107 double & a,
double & b,
const Selector & selector,
108 double exclude_above,
bool use_area_4vector)
const {
112 _check_selector_good_for_median(selector);
116 double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0;
118 vector<PseudoJet> incl_jets = inclusive_jets();
120 for (
unsigned i = 0; i < incl_jets.size(); i++) {
121 if (selector.
pass(incl_jets[i])) {
123 if ( use_area_4vector ) {
124 this_area = area_4vector(incl_jets[i]).perp();
126 this_area = area(incl_jets[i]);
128 double f = incl_jets[i].perp()/this_area;
129 if (exclude_above <= 0.0 || f < exclude_above) {
130 double x = incl_jets[i].rap();
double x2 = x*x;
154 b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
155 a = mean_f - b*mean_x2;
162 void ClusterSequenceAreaBase::get_median_rho_and_sigma(
163 const Selector & selector,
bool use_area_4vector,
164 double & median,
double & sigma,
double & mean_area)
const {
166 vector<PseudoJet> incl_jets = inclusive_jets();
167 get_median_rho_and_sigma(incl_jets, selector, use_area_4vector,
168 median, sigma, mean_area,
true);
172 void ClusterSequenceAreaBase::get_median_rho_and_sigma(
173 const vector<PseudoJet> & all_jets,
174 const Selector & selector,
bool use_area_4vector,
175 double & median,
double & sigma,
double & mean_area,
176 bool all_are_incl)
const {
178 _check_jet_alg_good_for_median();
183 _check_selector_good_for_median(selector);
185 vector<double> pt_over_areas;
186 double total_area = 0.0;
187 double total_njets = 0;
189 for (
unsigned i = 0; i < all_jets.size(); i++) {
190 if (selector.
pass(all_jets[i])) {
192 if (use_area_4vector) {
193 this_area = area_4vector(all_jets[i]).perp();
195 this_area = area(all_jets[i]);
199 pt_over_areas.push_back(all_jets[i].perp()/this_area);
201 _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).");
204 total_area += this_area;
210 if (pt_over_areas.size() == 0) {
220 sort(pt_over_areas.begin(), pt_over_areas.end());
224 double posn[2] = {0.5, (1.0-0.6827)/2.0};
227 double n_empty, empty_a;
228 if (has_explicit_ghosts()) {
236 }
else if (all_are_incl) {
238 empty_a = empty_area(selector);
239 n_empty = n_empty_jets(selector);
244 empty_a = empty_area_from_jets(all_jets, selector);
245 mean_area = total_area / total_njets;
246 n_empty = empty_a / mean_area;
250 total_njets += n_empty;
251 total_area += empty_a;
256 int pt_over_areas_size = pt_over_areas.size();
257 if (n_empty < -pt_over_areas_size/4.0)
258 _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.");
260 for (
int i = 0; i < 2; i++) {
261 double nj_median_pos =
262 (pt_over_areas_size-1.0 + n_empty)*posn[i] - n_empty;
263 double nj_median_ratio;
264 if (nj_median_pos >= 0 && pt_over_areas_size > 1) {
265 int int_nj_median = int(nj_median_pos);
268 if (int_nj_median+1 > pt_over_areas_size-1){
269 int_nj_median = pt_over_areas_size-2;
270 nj_median_pos = pt_over_areas_size-1;
274 pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
275 + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
277 nj_median_ratio = 0.0;
279 res[i] = nj_median_ratio;
282 double error = res[0] - res[1];
283 mean_area = total_area / total_njets;
284 sigma = error * sqrt(mean_area);
292 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
const double rho,
295 vector<PseudoJet> sub_jets;
296 vector<PseudoJet> jets_local =
sorted_by_pt(inclusive_jets(ptmin));
297 for (
unsigned i=0; i<jets_local.size(); i++) {
298 PseudoJet sub_jet = subtracted_jet(jets_local[i],rho);
299 sub_jets.push_back(sub_jet);
308 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
312 double rho = median_pt_per_unit_area_4vector(selector);
313 return subtracted_jets(rho,ptmin);
319 const double rho)
const {
323 if (rho*area4vect.
perp() < jet.
perp() ) {
324 sub_jet = jet - rho*area4vect;
325 }
else { sub_jet =
PseudoJet(0.0,0.0,0.0,0.0); }
343 double rho = median_pt_per_unit_area_4vector(selector);
344 PseudoJet sub_jet = subtracted_jet(jet, rho);
350 double ClusterSequenceAreaBase::subtracted_pt(
const PseudoJet & jet,
352 bool use_area_4vector)
const {
353 if ( use_area_4vector ) {
354 PseudoJet sub_jet = subtracted_jet(jet,rho);
355 return sub_jet.
perp();
357 return jet.
perp() - rho*area(jet);
365 double ClusterSequenceAreaBase::subtracted_pt(
const PseudoJet & jet,
367 bool use_area_4vector)
const {
368 if ( use_area_4vector ) {
369 PseudoJet sub_jet = subtracted_jet(jet,selector);
370 return sub_jet.
perp();
372 double rho = median_pt_per_unit_area(selector);
373 return subtracted_pt(jet,rho,
false);
379 void ClusterSequenceAreaBase::_check_selector_good_for_median(
const Selector &selector)
const{
382 throw Error(
"ClusterSequenceAreaBase: empty area can only be computed from selectors with a finite area");
387 throw Error(
"ClusterSequenceAreaBase: empty area can only be computed from selectors that apply jet by jet");
393 void ClusterSequenceAreaBase::_check_jet_alg_good_for_median()
const {
394 if (jet_def().jet_algorithm() != kt_algorithm
395 && jet_def().jet_algorithm() != cambridge_algorithm
396 && jet_def().jet_algorithm() != cambridge_for_passive_algorithm) {
397 _warnings.warn(
"ClusterSequenceAreaBase: jet_def being used may not be suitable for estimating diffuse backgrounds (good options are kt, cam)");
403 FASTJET_END_NAMESPACE
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
bool applies_jet_by_jet() const
returns true if this can be applied jet by jet
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...
const SharedPtr< PseudoJetStructureBase > & structure_shared_ptr() const
return a reference to the shared pointer to the PseudoJetStructureBase associated with this PseudoJet...
int cluster_hist_index() const
return the cluster_hist_index, intended to be used by clustering routines.
void set_cluster_hist_index(const int index)
set the cluster_hist_index, intended to be used by clustering routines.
base class corresponding to errors that can be thrown by FastJet
bool has_finite_area() const
returns true if it has a meaningful and finite area (i.e.
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
double area() const
returns the rapidity-phi area associated with the Selector (throws InvalidArea if the area does not m...
int user_index() const
return the user_index,
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure)
set the associated structure
double perp() const
returns the scalar transverse momentum
bool pass(const PseudoJet &jet) const
return true if the jet passes the selection
Class to contain pseudojets, including minimal information of use to jet-clustering routines...