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
00037 FASTJET_BEGIN_NAMESPACE
00038
00039 using namespace std;
00040
00041
00043 LimitedWarning ClusterSequenceAreaBase::_warnings;
00044
00045
00050 double ClusterSequenceAreaBase::empty_area(const RangeDefinition & range) const {
00051 double empty = range.area();
00052 vector<PseudoJet> incl_jets(inclusive_jets(0.0));
00053 for (unsigned i = 0; i < incl_jets.size(); i++) {
00054 if (range.is_in_range(incl_jets[i])) empty -= area(incl_jets[i]);
00055 }
00056 return empty;
00057 }
00058
00059 double ClusterSequenceAreaBase::median_pt_per_unit_area(const RangeDefinition & range) const {
00060 return median_pt_per_unit_something(range,false);
00061 }
00062
00063 double ClusterSequenceAreaBase::median_pt_per_unit_area_4vector(const RangeDefinition & range) const {
00064 return median_pt_per_unit_something(range,true);
00065 }
00066
00067
00068
00072 double ClusterSequenceAreaBase::median_pt_per_unit_something(
00073 const RangeDefinition & range, bool use_area_4vector) const {
00074
00075 double median, sigma, mean_area;
00076 get_median_rho_and_sigma(range, use_area_4vector, median, sigma, mean_area);
00077 return median;
00078
00079 }
00080
00081
00082
00086 void ClusterSequenceAreaBase::parabolic_pt_per_unit_area(
00087 double & a, double & b, const RangeDefinition & range,
00088 double exclude_above, bool use_area_4vector) const {
00089
00090 int n=0;
00091 int n_excluded = 0;
00092 double mean_f=0, mean_x2=0, mean_x4=0, mean_fx2=0;
00093
00094 vector<PseudoJet> incl_jets = inclusive_jets();
00095
00096 for (unsigned i = 0; i < incl_jets.size(); i++) {
00097 if (range.is_in_range(incl_jets[i])) {
00098 double this_area;
00099 if ( use_area_4vector ) {
00100 this_area = area_4vector(incl_jets[i]).perp();
00101 } else {
00102 this_area = area(incl_jets[i]);
00103 }
00104 double f = incl_jets[i].perp()/this_area;
00105 if (exclude_above <= 0.0 || f < exclude_above) {
00106 double x = incl_jets[i].rap(); double x2 = x*x;
00107 mean_f += f;
00108 mean_x2 += x2;
00109 mean_x4 += x2*x2;
00110 mean_fx2 += f*x2;
00111 n++;
00112 } else {
00113 n_excluded++;
00114 }
00115 }
00116 }
00117
00118 if (n <= 1) {
00119
00120
00121
00122 a = 0.0;
00123 b = 0.0;
00124 } else {
00125 mean_f /= n;
00126 mean_x2 /= n;
00127 mean_x4 /= n;
00128 mean_fx2 /= n;
00129
00130 b = (mean_f*mean_x2 - mean_fx2)/(mean_x2*mean_x2 - mean_x4);
00131 a = mean_f - b*mean_x2;
00132 }
00133
00134 }
00135
00136
00137
00138
00139 void ClusterSequenceAreaBase::get_median_rho_and_sigma(
00140 const RangeDefinition & range, bool use_area_4vector,
00141 double & median, double & sigma, double & mean_area) const {
00142
00143 _check_jet_alg_good_for_median();
00144
00145 vector<double> pt_over_areas;
00146 vector<PseudoJet> incl_jets = inclusive_jets();
00147 double total_area = 0.0;
00148 double total_njets = 0;
00149
00150 for (unsigned i = 0; i < incl_jets.size(); i++) {
00151 if (range.is_in_range(incl_jets[i])) {
00152 double this_area;
00153 if (use_area_4vector) {
00154 this_area = area_4vector(incl_jets[i]).perp();
00155 } else {
00156 this_area = area(incl_jets[i]);
00157 }
00158 pt_over_areas.push_back(incl_jets[i].perp()/this_area);
00159 total_area += this_area;
00160 total_njets += 1.0;
00161 }
00162 }
00163
00164
00165 if (pt_over_areas.size() == 0) {
00166 median = 0.0;
00167 sigma = 0.0;
00168 mean_area = 0.0;
00169 return;
00170 }
00171
00172
00173
00174
00175 sort(pt_over_areas.begin(), pt_over_areas.end());
00176
00177
00178
00179 double posn[2] = {0.5, (1.0-0.6827)/2.0};
00180 double res[2];
00181
00182 double n_empty = n_empty_jets(range);
00183 total_njets += n_empty;
00184 total_area += empty_area(range);
00185
00186 for (int i = 0; i < 2; i++) {
00187 double nj_median_pos =
00188 (pt_over_areas.size()-1 + n_empty)*posn[i] - n_empty;
00189 double nj_median_ratio;
00190 if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
00191 int int_nj_median = int(nj_median_pos);
00192 nj_median_ratio =
00193 pt_over_areas[int_nj_median] * (int_nj_median+1-nj_median_pos)
00194 + pt_over_areas[int_nj_median+1] * (nj_median_pos - int_nj_median);
00195 } else {
00196 nj_median_ratio = 0.0;
00197 }
00198 res[i] = nj_median_ratio;
00199 }
00200 median = res[0];
00201 double error = res[0] - res[1];
00202 mean_area = total_area / total_njets;
00203 sigma = error * sqrt(mean_area);
00204 }
00205
00206
00211 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(const double rho,
00212 const double ptmin)
00213 const {
00214 vector<PseudoJet> sub_jets;
00215 vector<PseudoJet> jets = sorted_by_pt(inclusive_jets(ptmin));
00216 for (unsigned i=0; i<jets.size(); i++) {
00217 PseudoJet sub_jet = subtracted_jet(jets[i],rho);
00218 sub_jets.push_back(sub_jet);
00219 }
00220 return sub_jets;
00221 }
00222
00227 vector<PseudoJet> ClusterSequenceAreaBase::subtracted_jets(
00228 const RangeDefinition & range,
00229 const double ptmin)
00230 const {
00231 double rho = median_pt_per_unit_area_4vector(range);
00232 return subtracted_jets(rho,ptmin);
00233 }
00234
00235
00237 PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
00238 const double rho) const {
00239 PseudoJet area4vect = area_4vector(jet);
00240 PseudoJet sub_jet;
00241
00242 if (rho*area4vect.perp() < jet.perp() ) {
00243 sub_jet = jet - rho*area4vect;
00244 } else { sub_jet = PseudoJet(0.0,0.0,0.0,0.0); }
00245 return sub_jet;
00246 }
00247
00248
00252 PseudoJet ClusterSequenceAreaBase::subtracted_jet(const PseudoJet & jet,
00253 const RangeDefinition & range) const {
00254 double rho = median_pt_per_unit_area_4vector(range);
00255 PseudoJet sub_jet = subtracted_jet(jet, rho);
00256 return sub_jet;
00257 }
00258
00259
00261 double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
00262 const double rho,
00263 bool use_area_4vector) const {
00264 if ( use_area_4vector ) {
00265 PseudoJet sub_jet = subtracted_jet(jet,rho);
00266 return sub_jet.perp();
00267 } else {
00268 return jet.perp() - rho*area(jet);
00269 }
00270 }
00271
00272
00276 double ClusterSequenceAreaBase::subtracted_pt(const PseudoJet & jet,
00277 const RangeDefinition & range,
00278 bool use_area_4vector) const {
00279 if ( use_area_4vector ) {
00280 PseudoJet sub_jet = subtracted_jet(jet,range);
00281 return sub_jet.perp();
00282 } else {
00283 double rho = median_pt_per_unit_area(range);
00284 return subtracted_pt(jet,rho,false);
00285 }
00286 }
00287
00288
00290 void ClusterSequenceAreaBase::_check_jet_alg_good_for_median() const {
00291 if (jet_def().jet_algorithm() != kt_algorithm
00292 && jet_def().jet_algorithm() != cambridge_algorithm
00293 && jet_def().jet_algorithm() != cambridge_for_passive_algorithm) {
00294 _warnings.warn("ClusterSequenceAreaBase: jet_def being used may not be suitable for estimating diffuse backgrounds (good options are kt, cam)");
00295 }
00296 }
00297
00298
00299
00300 FASTJET_END_NAMESPACE