00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035 #include "fastjet/ClusterSequenceAreaBase.hh"
00036 #include <algorithm>
00037
00038 FASTJET_BEGIN_NAMESPACE
00039
00040 using namespace std;
00041
00042
00044 LimitedWarning ClusterSequenceAreaBase::_warnings;
00045 LimitedWarning ClusterSequenceAreaBase::_warnings_zero_area;
00046
00047
00055 double ClusterSequenceAreaBase::empty_area(const RangeDefinition & range) const {
00056
00057 if (has_explicit_ghosts()) {return 0.0;}
00058 else { return empty_area_from_jets(inclusive_jets(0.0), range);}
00059
00060 }
00061
00062
00067 double ClusterSequenceAreaBase::empty_area_from_jets(
00068 const std::vector<PseudoJet> & all_jets,
00069 const RangeDefinition & range) const {
00070
00071 double empty = range.area();
00072 for (unsigned i = 0; i < all_jets.size(); i++) {
00073 if (range.is_in_range(all_jets[i])) empty -= area(all_jets[i]);
00074 }
00075 return empty;
00076 }
00077
00078 double ClusterSequenceAreaBase::median_pt_per_unit_area(const RangeDefinition & range) const {
00079 return median_pt_per_unit_something(range,false);
00080 }
00081
00082 double ClusterSequenceAreaBase::median_pt_per_unit_area_4vector(const RangeDefinition & range) const {
00083 return median_pt_per_unit_something(range,true);
00084 }
00085
00086
00087
00091 double ClusterSequenceAreaBase::median_pt_per_unit_something(
00092 const RangeDefinition & range, bool use_area_4vector) const {
00093
00094 double median, sigma, mean_area;
00095 get_median_rho_and_sigma(range, use_area_4vector, median, sigma, mean_area);
00096 return median;
00097
00098 }
00099
00100
00101
00105 void ClusterSequenceAreaBase::parabolic_pt_per_unit_area(
00106 double & a, double & b, const RangeDefinition & range,
00107 double exclude_above, bool use_area_4vector) const {
00108
00109 int n=0;
00110 int n_excluded = 0;
00111 double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0;
00112
00113 vector<PseudoJet> incl_jets = inclusive_jets();
00114
00115 for (unsigned i = 0; i < incl_jets.size(); i++) {
00116 if (range.is_in_range(incl_jets[i])) {
00117 double this_area;
00118 if ( use_area_4vector ) {
00119 this_area = area_4vector(incl_jets[i]).perp();
00120 } else {
00121 this_area = area(incl_jets[i]);
00122 }
00123 double f = incl_jets[i].perp()/this_area;
00124 if (exclude_above <= 0.0 || f < exclude_above) {
00125 double x = incl_jets[i].rap(); double x2 = x*x;
00126 mean_f += f;
00127 mean_x2 += x2;
00128 mean_x4 += x2*x2;
00129 mean_fx2 += f*x2;
00130 n++;
00131 } else {
00132 n_excluded++;
00133 }
00134 }
00135 }
00136
00137 if (n <= 1) {
00138
00139
00140
00141 a = 0.0;
00142 b = 0.0;
00143 } else {
00144 mean_f /= n;
00145 mean_x2 /= n;
00146 mean_x4 /= n;
00147 mean_fx2 /= n;
00148
00149 b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
00150 a = mean_f - b*mean_x2;
00151 }
00152
00153 }
00154
00155
00156
00157 void ClusterSequenceAreaBase::get_median_rho_and_sigma(
00158 const RangeDefinition & range, bool use_area_4vector,
00159 double & median, double & sigma, double & mean_area) const {
00160
00161 vector<PseudoJet> incl_jets = inclusive_jets();
00162 get_median_rho_and_sigma(incl_jets, range, use_area_4vector,
00163 median, sigma, mean_area, true);
00164 }
00165
00166
00167 void ClusterSequenceAreaBase::get_median_rho_and_sigma(
00168 const vector<PseudoJet> & all_jets,
00169 const RangeDefinition & range, bool use_area_4vector,
00170 double & median, double & sigma, double & mean_area,
00171 bool all_are_incl) const {
00172
00173 _check_jet_alg_good_for_median();
00174
00175 vector<double> pt_over_areas;
00176 double total_area = 0.0;
00177 double total_njets = 0;
00178
00179 for (unsigned i = 0; i < all_jets.size(); i++) {
00180 if (range.is_in_range(all_jets[i])) {
00181 double this_area;
00182 if (use_area_4vector) {
00183 this_area = area_4vector(all_jets[i]).perp();
00184 } else {
00185 this_area = area(all_jets[i]);
00186 }
00187
00188 if (this_area>0) {
00189 pt_over_areas.push_back(all_jets[i].perp()/this_area);
00190 } else {
00191 _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).");
00192 }
00193
00194 total_area += this_area;
00195 total_njets += 1.0;
00196 }
00197 }
00198
00199
00200 if (pt_over_areas.size() == 0) {
00201 median = 0.0;
00202 sigma = 0.0;
00203 mean_area = 0.0;
00204 return;
00205 }
00206
00207
00208
00209
00210 sort(pt_over_areas.begin(), pt_over_areas.end());
00211
00212
00213
00214 double posn[2] = {0.5, (1.0-0.6827)/2.0};
00215 double res[2];
00216
00217 double n_empty, empty_a;
00218 if (has_explicit_ghosts()) {
00219
00220
00221
00222
00223
00224 empty_a = 0.0;
00225 n_empty = 0;
00226 } else if (all_are_incl) {
00227
00228 empty_a = empty_area(range);
00229 n_empty = n_empty_jets(range);
00230 } else {
00231
00232
00233
00234 empty_a = empty_area_from_jets(all_jets, range);
00235 mean_area = total_area / total_njets;
00236 n_empty = empty_a / mean_area;
00237 }
00238
00239
00240 total_njets += n_empty;
00241 total_area += empty_a;
00242
00243 for (int i = 0; i < 2; i++) {
00244 double nj_median_pos =
00245 (pt_over_areas.size()-1 + n_empty)*posn[i] - n_empty;
00246 double nj_median_ratio;
00247 if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
00248 int int_nj_median = int(nj_median_pos);
00249 nj_median_ratio =
00250 pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
00251 + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
00252 } else {
00253 nj_median_ratio = 0.0;
00254 }
00255 res[i] = nj_median_ratio;
00256 }
00257 median = res[0];
00258 double error = res[0] - res[1];
00259 mean_area = total_area / total_njets;
00260 sigma = error * sqrt(mean_area);
00261 }
00262
00263
00268 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(const double rho,
00269 const double ptmin)
00270 const {
00271 vector<PseudoJet> sub_jets;
00272 vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
00273 for (unsigned i=0; i<jets.size(); i++) {
00274 PseudoJet sub_jet = subtracted_jet(jets[i],rho);
00275 sub_jets.push_back(sub_jet);
00276 }
00277 return sub_jets;
00278 }
00279
00284 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
00285 const RangeDefinition & range,
00286 const double ptmin)
00287 const {
00288 double rho = median_pt_per_unit_area_4vector(range);
00289 return subtracted_jets(rho,ptmin);
00290 }
00291
00292
00294 PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
00295 const double rho) const {
00296 PseudoJet area4vect = area_4vector(jet);
00297 PseudoJet sub_jet;
00298
00299 if (rho*area4vect.perp() < jet.perp() ) {
00300 sub_jet = jet - rho*area4vect;
00301 } else { sub_jet = PseudoJet(0.0,0.0,0.0,0.0); }
00302
00303
00304
00305 sub_jet.set_cluster_hist_index(jet.cluster_hist_index());
00306 sub_jet.set_user_index(jet.user_index());
00307
00308 return sub_jet;
00309 }
00310
00311
00315 PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
00316 const RangeDefinition & range) const {
00317 double rho = median_pt_per_unit_area_4vector(range);
00318 PseudoJet sub_jet = subtracted_jet(jet, rho);
00319 return sub_jet;
00320 }
00321
00322
00324 double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
00325 const double rho,
00326 bool use_area_4vector) const {
00327 if ( use_area_4vector ) {
00328 PseudoJet sub_jet = subtracted_jet(jet,rho);
00329 return sub_jet.perp();
00330 } else {
00331 return jet.perp() - rho*area(jet);
00332 }
00333 }
00334
00335
00339 double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
00340 const RangeDefinition & range,
00341 bool use_area_4vector) const {
00342 if ( use_area_4vector ) {
00343 PseudoJet sub_jet = subtracted_jet(jet,range);
00344 return sub_jet.perp();
00345 } else {
00346 double rho = median_pt_per_unit_area(range);
00347 return subtracted_pt(jet,rho,false);
00348 }
00349 }
00350
00351
00353 void ClusterSequenceAreaBase::_check_jet_alg_good_for_median() const {
00354 if (jet_def().jet_algorithm() != kt_algorithm
00355 && jet_def().jet_algorithm() != cambridge_algorithm
00356 && jet_def().jet_algorithm() != cambridge_for_passive_algorithm) {
00357 _warnings.warn("ClusterSequenceAreaBase: jet_def being used may not be suitable for estimating diffuse backgrounds (good options are kt, cam)");
00358 }
00359 }
00360
00361
00362
00363 FASTJET_END_NAMESPACE