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