1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403
|
//STARTHEADER
// $Id: ClusterSequenceAreaBase.cc 2687 2011-11-14 11:17:51Z soyez $
//
// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// The algorithms that underlie FastJet have required considerable
// development and are described in hep-ph/0512210. If you use
// FastJet as part of work towards a scientific publication, please
// include a citation to the FastJet paper.
//
// FastJet is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//ENDHEADER
#include "fastjet/ClusterSequenceAreaBase.hh"
#include <algorithm>
FASTJET_BEGIN_NAMESPACE
using namespace std;
/// allow for warnings
LimitedWarning ClusterSequenceAreaBase::_warnings;
LimitedWarning ClusterSequenceAreaBase::_warnings_zero_area;
LimitedWarning ClusterSequenceAreaBase::_warnings_empty_area;
//----------------------------------------------------------------------
/// return the total area, within the selector's range, that is free
/// of jets.
///
/// Calculate this as (range area) - \sum_{i in range} A_i
///
/// for ClusterSequences with explicit ghosts, assume that there will
/// never be any empty area, i.e. it is always filled in by pure
/// ghosts jets. This holds for seq.rec. algorithms
double ClusterSequenceAreaBase::empty_area(const Selector & selector) const {
if (has_explicit_ghosts()) {return 0.0;}
else { return empty_area_from_jets(inclusive_jets(0.0), selector);}
}
//----------------------------------------------------------------------
/// return the total area, within range, that is free of jets.
///
/// Calculate this as (range area) - \sum_{i in range} A_i
///
double ClusterSequenceAreaBase::empty_area_from_jets(
const std::vector<PseudoJet> & all_jets,
const Selector & selector) const {
_check_selector_good_for_median(selector);
double empty = selector.area();
for (unsigned i = 0; i < all_jets.size(); i++) {
if (selector.pass(all_jets[i])) empty -= area(all_jets[i]);
}
return empty;
}
double ClusterSequenceAreaBase::median_pt_per_unit_area(const Selector & selector) const {
return median_pt_per_unit_something(selector,false);
}
double ClusterSequenceAreaBase::median_pt_per_unit_area_4vector(const Selector & selector) const {
return median_pt_per_unit_something(selector,true);
}
//----------------------------------------------------------------------
/// the median of (pt/area) for jets contained within range, counting
/// the empty area as if it were made up of a collection of empty
/// jets each of area (0.55 * pi R^2).
double ClusterSequenceAreaBase::median_pt_per_unit_something(
const Selector & selector, bool use_area_4vector) const {
double median, sigma, mean_area;
get_median_rho_and_sigma(selector, use_area_4vector, median, sigma, mean_area);
return median;
}
//----------------------------------------------------------------------
/// fits a form pt_per_unit_area(y) = a + b*y^2 for jets in range.
/// exclude_above allows one to exclude large values of pt/area from fit.
/// use_area_4vector = true uses the 4vector areas.
void ClusterSequenceAreaBase::parabolic_pt_per_unit_area(
double & a, double & b, const Selector & selector,
double exclude_above, bool use_area_4vector) const {
// sanity check on the selector: we require a finite area and that
// it applies jet by jet (see BackgroundEstimator for more advanced
// usage)
_check_selector_good_for_median(selector);
int n=0;
int n_excluded = 0;
double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0;
vector<PseudoJet> incl_jets = inclusive_jets();
for (unsigned i = 0; i < incl_jets.size(); i++) {
if (selector.pass(incl_jets[i])) {
double this_area;
if ( use_area_4vector ) {
this_area = area_4vector(incl_jets[i]).perp();
} else {
this_area = area(incl_jets[i]);
}
double f = incl_jets[i].perp()/this_area;
if (exclude_above <= 0.0 || f < exclude_above) {
double x = incl_jets[i].rap(); double x2 = x*x;
mean_f += f;
mean_x2 += x2;
mean_x4 += x2*x2;
mean_fx2 += f*x2;
n++;
} else {
n_excluded++;
}
}
}
if (n <= 1) {
// meaningful results require at least two jets inside the
// area -- mind you if there are empty jets we should be in
// any case doing something special...
a = 0.0;
b = 0.0;
} else {
mean_f /= n;
mean_x2 /= n;
mean_x4 /= n;
mean_fx2 /= n;
b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
a = mean_f - b*mean_x2;
}
//cerr << "n_excluded = "<< n_excluded << endl;
}
void ClusterSequenceAreaBase::get_median_rho_and_sigma(
const Selector & selector, bool use_area_4vector,
double & median, double & sigma, double & mean_area) const {
vector<PseudoJet> incl_jets = inclusive_jets();
get_median_rho_and_sigma(incl_jets, selector, use_area_4vector,
median, sigma, mean_area, true);
}
void ClusterSequenceAreaBase::get_median_rho_and_sigma(
const vector<PseudoJet> & all_jets,
const Selector & selector, bool use_area_4vector,
double & median, double & sigma, double & mean_area,
bool all_are_incl) const {
_check_jet_alg_good_for_median();
// sanity check on the selector: we require a finite area and that
// it applies jet by jet (see BackgroundEstimator for more advanced
// usage)
_check_selector_good_for_median(selector);
vector<double> pt_over_areas;
double total_area = 0.0;
double total_njets = 0;
for (unsigned i = 0; i < all_jets.size(); i++) {
if (selector.pass(all_jets[i])) {
double this_area;
if (use_area_4vector) {
this_area = area_4vector(all_jets[i]).perp();
} else {
this_area = area(all_jets[i]);
}
if (this_area>0) {
pt_over_areas.push_back(all_jets[i].perp()/this_area);
} else {
_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).");
}
total_area += this_area;
total_njets += 1.0;
}
}
// there is nothing inside our region, so answer will always be zero
if (pt_over_areas.size() == 0) {
median = 0.0;
sigma = 0.0;
mean_area = 0.0;
return;
}
// 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; it will be supplemented with more info below]
sort(pt_over_areas.begin(), pt_over_areas.end());
// now get the median & error, accounting for empty jets
// define the fractions of distribution at median, median-1sigma
double posn[2] = {0.5, (1.0-0.6827)/2.0};
double res[2];
double n_empty, empty_a;
if (has_explicit_ghosts()) {
// NB: the following lines of code are potentially incorrect in cases
// where there are unclustered particles (empty_area would do a better job,
// at least for active areas). This is not an issue with kt or C/A, or other
// algorithms that cluster all particles (and the median estimation should in
// any case only be done with kt or C/A!)
empty_a = 0.0;
n_empty = 0;
} else if (all_are_incl) {
// the default case
empty_a = empty_area(selector);
n_empty = n_empty_jets(selector);
} else {
// this one is intended to be used when e.g. one runs C/A, then looks at its
// exclusive jets in order to get an effective smaller R value, and passes those
// to this routine.
empty_a = empty_area_from_jets(all_jets, selector);
mean_area = total_area / total_njets; // temporary value
n_empty = empty_a / mean_area;
}
//cout << "*** tot_area = " << total_area << ", empty_a = " << empty_a << endl;
//cout << "*** n_empty = " << n_empty << ", ntotal = " << total_njets << endl;
total_njets += n_empty;
total_area += empty_a;
// we need an int (rather than an unsigned int) with the size of the
// pt_over_areas array, because we'll often be doing subtraction of
// -1, negating it, etc. All of these operations go crazy with unsigned ints.
int pt_over_areas_size = pt_over_areas.size();
if (n_empty < -pt_over_areas_size/4.0)
_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.");
for (int i = 0; i < 2; i++) {
double nj_median_pos =
(pt_over_areas_size-1.0 + n_empty)*posn[i] - n_empty;
double nj_median_ratio;
if (nj_median_pos >= 0 && pt_over_areas_size > 1) {
int int_nj_median = int(nj_median_pos);
// avoid potential overflow issues
if (int_nj_median+1 > pt_over_areas_size-1){
int_nj_median = pt_over_areas_size-2;
nj_median_pos = pt_over_areas_size-1;
}
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;
}
res[i] = nj_median_ratio;
}
median = res[0];
double error = res[0] - res[1];
mean_area = total_area / total_njets;
sigma = error * sqrt(mean_area);
}
/// return a vector of all subtracted jets, using area_4vector, given rho.
/// Only inclusive_jets above ptmin are subtracted and returned.
/// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()),
/// i.e. not necessarily ordered in pt once subtracted
vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(const double rho,
const double ptmin)
const {
vector<PseudoJet> sub_jets;
vector<PseudoJet> jets_local = sorted_by_pt(inclusive_jets(ptmin));
for (unsigned i=0; i<jets_local.size(); i++) {
PseudoJet sub_jet = subtracted_jet(jets_local[i],rho);
sub_jets.push_back(sub_jet);
}
return sub_jets;
}
/// return a vector of subtracted jets, using area_4vector.
/// Only inclusive_jets above ptmin are subtracted and returned.
/// the ordering is the same as that of sorted_by_pt(cs.inclusive_jets()),
/// i.e. not necessarily ordered in pt once subtracted
vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
const Selector & selector,
const double ptmin)
const {
double rho = median_pt_per_unit_area_4vector(selector);
return subtracted_jets(rho,ptmin);
}
/// return a subtracted jet, using area_4vector, given rho
PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
const double rho) const {
PseudoJet area4vect = area_4vector(jet);
PseudoJet sub_jet;
// sanity check
if (rho*area4vect.perp() < jet.perp() ) {
sub_jet = jet - rho*area4vect;
} else { sub_jet = PseudoJet(0.0,0.0,0.0,0.0); }
// make sure the subtracted jet has the same index (cluster, user, csw)
// (i.e. "looks like") the original jet
sub_jet.set_cluster_hist_index(jet.cluster_hist_index());
sub_jet.set_user_index(jet.user_index());
// do not use CS::_set_structure_shared_ptr here, which should
// only be called to maintain the tally during construction
sub_jet.set_structure_shared_ptr(jet.structure_shared_ptr());
return sub_jet;
}
/// return a subtracted jet, using area_4vector; note that this is
/// potentially inefficient if repeatedly used for many different
/// jets, because rho will be recalculated each time around.
PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
const Selector & selector) const {
double rho = median_pt_per_unit_area_4vector(selector);
PseudoJet sub_jet = subtracted_jet(jet, rho);
return sub_jet;
}
/// return the subtracted pt, given rho
double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
const double rho,
bool use_area_4vector) const {
if ( use_area_4vector ) {
PseudoJet sub_jet = subtracted_jet(jet,rho);
return sub_jet.perp();
} else {
return jet.perp() - rho*area(jet);
}
}
/// return the subtracted pt; note that this is
/// potentially inefficient if repeatedly used for many different
/// jets, because rho will be recalculated each time around.
double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
const Selector & selector,
bool use_area_4vector) const {
if ( use_area_4vector ) {
PseudoJet sub_jet = subtracted_jet(jet,selector);
return sub_jet.perp();
} else {
double rho = median_pt_per_unit_area(selector);
return subtracted_pt(jet,rho,false);
}
}
// check the selector is suited for the computations i.e. applies jet
// by jet and has a finite area
void ClusterSequenceAreaBase::_check_selector_good_for_median(const Selector &selector) const{
// make sure the selector has a finite area
if ((! has_explicit_ghosts()) && (! selector.has_finite_area())){
throw Error("ClusterSequenceAreaBase: empty area can only be computed from selectors with a finite area");
}
// make sure the selector applies jet by jet
if (! selector.applies_jet_by_jet()){
throw Error("ClusterSequenceAreaBase: empty area can only be computed from selectors that apply jet by jet");
}
}
/// check the jet algorithm is suitable (and if not issue a warning)
void ClusterSequenceAreaBase::_check_jet_alg_good_for_median() const {
if (jet_def().jet_algorithm() != kt_algorithm
&& jet_def().jet_algorithm() != cambridge_algorithm
&& jet_def().jet_algorithm() != cambridge_for_passive_algorithm) {
_warnings.warn("ClusterSequenceAreaBase: jet_def being used may not be suitable for estimating diffuse backgrounds (good options are kt, cam)");
}
}
FASTJET_END_NAMESPACE
|