JetProducers/src/BackgroundEstimator.cc

0001 #include "RecoJets/JetProducers/interface/BackgroundEstimator.h"
0002 #include <fastjet/ClusterSequenceAreaBase.hh>
0003 #include <iostream>
0004
0005 using namespace fastjet;
0006 using namespace std;
0007
0008 //---------------------------------------------------------------------
0009 // \class BackgroundEstimator
0010 // Class to estimate the density of the background per unit area
0011 //
0012 // The default behaviour of this class is to compute the global
0013 // properties of the background as it is done in ClusterSequenceArea.
0014 // On top of that, we provide methods to specify an explicit set of
0015 // jets to use or a list of jets to exclude.
0016 // We also provide all sorts of additional information regarding
0017 // the background estimation like the jets that have been used or
0018 // the number of pure-ghost jets.
0019 //---------------------------------------------------------------------
0020
0021 // default ctor
0022 //  - csa      the ClusterSequenceArea to use
0023 //  - range    the range over which jets will be considered
0024 BackgroundEstimator::BackgroundEstimator(const ClusterSequenceAreaBase &csa, const RangeDefinition &range)
0025     : _csa(csa), _range(range) {
0026   reset();
0027 }
0028
0029 // default dtor
0030 BackgroundEstimator::~BackgroundEstimator() {}
0031
0032 // reset to default values
0033 // set the list of included jets to the inclusive jets and clear the excluded ones
0034 void BackgroundEstimator::reset() {
0035   // set the list of included jets to the inclusive jets
0036   _included_jets = _csa.inclusive_jets();
0037   _all_from_inclusive = true;
0038   //set_included_jets(_csa.inclusive_jets());
0039
0040   // clear the list of explicitly excluded jets
0041   _excluded_jets.clear();
0042   //set_excluded_jets(vector<PseudoJet>());
0043
0044   // set the remaining default parameters
0045   set_use_area_4vector();  // true by default
0046
0047   // reset the computed values
0048   _median_rho = _sigma = _mean_area = 0.0;
0049   _n_jets_used = _n_jets_excluded = _n_empty_jets = 0;
0050   _empty_area = 0.0;
0051   _uptodate = false;
0052 }
0053
0054 // do the actual job
0055 void BackgroundEstimator::_compute() {
0056   //TODO: check that the alg is OK for median computation
0057   //_check_jet_alg_good_for_median();
0058
0059   // fill the vector of pt/area with the jets
0060   //  - in included_jets
0061   //  - not in excluded_jets
0062   //  - in the range
0063   vector<double> pt_over_areas;
0064   double total_area = 0.0;
0065
0066   _n_jets_used = 0;
0067   _n_jets_excluded = 0;
0068
0069   for (unsigned i = 0; i < _included_jets.size(); i++) {
0070     const PseudoJet &current_jet = _included_jets[i];
0071
0072     // check that the jet is not explicitly excluded
0073     // we'll compare them using their cluster_history_index
0074     bool excluded = false;
0075     int ref_idx = current_jet.cluster_hist_index();
0076     for (unsigned int j = 0; j < _excluded_jets.size(); j++)
0077       excluded |= (_excluded_jets[j].cluster_hist_index() == ref_idx);
0078
0079     // check if the jet is in the range
0080     if (_range.is_in_range(current_jet)) {
0081       if (excluded) {
0082         // keep track of the explicitly excluded jets
0083         _n_jets_excluded++;
0084       } else {
0085         double this_area = (_use_area_4vector) ? _csa.area_4vector(current_jet).perp() : _csa.area(current_jet);
0086
0087         pt_over_areas.push_back(current_jet.perp() / this_area);
0088         total_area += this_area;
0089         _n_jets_used++;
0090       }
0091     }
0092   }
0093
0094   // there is nothing inside our region, so answer will always be zero
0095   if (pt_over_areas.empty()) {
0096     _median_rho = 0.0;
0097     _sigma = 0.0;
0098     _mean_area = 0.0;
0099     return;
0100   }
0101
0102   // get median (pt/area) [this is the "old" median definition. It considers
0103   // only the "real" jets in calculating the median, i.e. excluding the
0104   // only-ghost ones; it will be supplemented with more info below]
0105   sort(pt_over_areas.begin(), pt_over_areas.end());
0106
0107   // determine the number of empty jets
0108   _empty_area = 0.0;
0109   _n_empty_jets = 0.0;
0110   if (_csa.has_explicit_ghosts()) {
0111     _empty_area = 0.0;
0112     _n_empty_jets = 0;
0113   } else if (_all_from_inclusive) {
0114     _empty_area = _csa.empty_area(_range);
0115     _n_empty_jets = _csa.n_empty_jets(_range);
0116   } else {
0117     _empty_area = _csa.empty_area_from_jets(_included_jets, _range);
0118     _mean_area = total_area / _n_jets_used;  // temporary value
0119     _n_empty_jets = _empty_area / _mean_area;
0120   }
0121
0122   double total_njets = _n_jets_used + _n_empty_jets;
0123   total_area += _empty_area;
0124
0125   // now get the median & error, accounting for empty jets
0126   // define the fractions of distribution at median, median-1sigma
0127   double posn[2] = {0.5, (1.0 - 0.6827) / 2.0};
0128   double res[2];
0129
0130   for (int i = 0; i < 2; i++) {
0131     double nj_median_pos = (total_njets - 1) * posn[i] - _n_empty_jets;
0132     double nj_median_ratio;
0133     if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
0134       int int_nj_median = int(nj_median_pos);
0135       nj_median_ratio = pt_over_areas[int_nj_median] * (int_nj_median + 1 - nj_median_pos) +
0136                         pt_over_areas[int_nj_median + 1] * (nj_median_pos - int_nj_median);
0137     } else {
0138       nj_median_ratio = 0.0;
0139     }
0140     res[i] = nj_median_ratio;
0141   }
0142
0143   // store the results
0144   double error = res[0] - res[1];
0145   _median_rho = res[0];
0146   _mean_area = total_area / total_njets;
0147   _sigma = error * sqrt(_mean_area);
0148
0149   // record that the computation has been performed
0150   _uptodate = true;
0151 }