CMS 3D CMS Logo

BackgroundEstimator.cc
Go to the documentation of this file.
2 #include <fastjet/ClusterSequenceAreaBase.hh>
3 #include <iostream>
4 
5 using namespace fastjet;
6 using namespace std;
7 
8 //---------------------------------------------------------------------
9 // \class BackgroundEstimator
10 // Class to estimate the density of the background per unit area
11 //
12 // The default behaviour of this class is to compute the global
13 // properties of the background as it is done in ClusterSequenceArea.
14 // On top of that, we provide methods to specify an explicit set of
15 // jets to use or a list of jets to exclude.
16 // We also provide all sorts of additional information regarding
17 // the background estimation like the jets that have been used or
18 // the number of pure-ghost jets.
19 //---------------------------------------------------------------------
20 
21 // default ctor
22 // - csa the ClusterSequenceArea to use
23 // - range the range over which jets will be considered
24 BackgroundEstimator::BackgroundEstimator(const ClusterSequenceAreaBase &csa, const RangeDefinition &range)
25  : _csa(csa), _range(range) {
26  reset();
27 }
28 
29 // default dtor
31 
32 // reset to default values
33 // set the list of included jets to the inclusive jets and clear the excluded ones
35  // set the list of included jets to the inclusive jets
36  _included_jets = _csa.inclusive_jets();
37  _all_from_inclusive = true;
38  //set_included_jets(_csa.inclusive_jets());
39 
40  // clear the list of explicitly excluded jets
41  _excluded_jets.clear();
42  //set_excluded_jets(vector<PseudoJet>());
43 
44  // set the remaining default parameters
45  set_use_area_4vector(); // true by default
46 
47  // reset the computed values
48  _median_rho = _sigma = _mean_area = 0.0;
50  _empty_area = 0.0;
51  _uptodate = false;
52 }
53 
54 // do the actual job
56  //TODO: check that the alg is OK for median computation
57  //_check_jet_alg_good_for_median();
58 
59  // fill the vector of pt/area with the jets
60  // - in included_jets
61  // - not in excluded_jets
62  // - in the range
63  vector<double> pt_over_areas;
64  double total_area = 0.0;
65 
66  _n_jets_used = 0;
67  _n_jets_excluded = 0;
68 
69  for (unsigned i = 0; i < _included_jets.size(); i++) {
70  const PseudoJet &current_jet = _included_jets[i];
71 
72  // check that the jet is not explicitly excluded
73  // we'll compare them using their cluster_history_index
74  bool excluded = false;
75  int ref_idx = current_jet.cluster_hist_index();
76  for (unsigned int j = 0; j < _excluded_jets.size(); j++)
77  excluded |= (_excluded_jets[j].cluster_hist_index() == ref_idx);
78 
79  // check if the jet is in the range
80  if (_range.is_in_range(current_jet)) {
81  if (excluded) {
82  // keep track of the explicitly excluded jets
84  } else {
85  double this_area = (_use_area_4vector) ? _csa.area_4vector(current_jet).perp() : _csa.area(current_jet);
86 
87  pt_over_areas.push_back(current_jet.perp() / this_area);
88  total_area += this_area;
89  _n_jets_used++;
90  }
91  }
92  }
93 
94  // there is nothing inside our region, so answer will always be zero
95  if (pt_over_areas.empty()) {
96  _median_rho = 0.0;
97  _sigma = 0.0;
98  _mean_area = 0.0;
99  return;
100  }
101 
102  // get median (pt/area) [this is the "old" median definition. It considers
103  // only the "real" jets in calculating the median, i.e. excluding the
104  // only-ghost ones; it will be supplemented with more info below]
105  sort(pt_over_areas.begin(), pt_over_areas.end());
106 
107  // determine the number of empty jets
108  _empty_area = 0.0;
109  _n_empty_jets = 0.0;
110  if (_csa.has_explicit_ghosts()) {
111  _empty_area = 0.0;
112  _n_empty_jets = 0;
113  } else if (_all_from_inclusive) {
114  _empty_area = _csa.empty_area(_range);
115  _n_empty_jets = _csa.n_empty_jets(_range);
116  } else {
117  _empty_area = _csa.empty_area_from_jets(_included_jets, _range);
118  _mean_area = total_area / _n_jets_used; // temporary value
120  }
121 
122  double total_njets = _n_jets_used + _n_empty_jets;
123  total_area += _empty_area;
124 
125  // now get the median & error, accounting for empty jets
126  // define the fractions of distribution at median, median-1sigma
127  double posn[2] = {0.5, (1.0 - 0.6827) / 2.0};
128  double res[2];
129 
130  for (int i = 0; i < 2; i++) {
131  double nj_median_pos = (total_njets - 1) * posn[i] - _n_empty_jets;
132  double nj_median_ratio;
133  if (nj_median_pos >= 0 && pt_over_areas.size() > 1) {
134  int int_nj_median = int(nj_median_pos);
135  nj_median_ratio = pt_over_areas[int_nj_median] * (int_nj_median + 1 - nj_median_pos) +
136  pt_over_areas[int_nj_median + 1] * (nj_median_pos - int_nj_median);
137  } else {
138  nj_median_ratio = 0.0;
139  }
140  res[i] = nj_median_ratio;
141  }
142 
143  // store the results
144  double error = res[0] - res[1];
145  _median_rho = res[0];
146  _mean_area = total_area / total_njets;
148 
149  // record that the computation has been performed
150  _uptodate = true;
151 }
FastTimerService_cff.range
range
Definition: FastTimerService_cff.py:34
BackgroundEstimator.h
mps_fire.i
i
Definition: mps_fire.py:428
fastjet
Definition: BackgroundEstimator.h:8
fastjet::BackgroundEstimator::_n_jets_used
unsigned int _n_jets_used
number of jets used to estimate the background
Definition: BackgroundEstimator.h:159
fastjet::BackgroundEstimator::_sigma
double _sigma
background estimated fluctuations
Definition: BackgroundEstimator.h:157
fastjet::BackgroundEstimator::set_use_area_4vector
void set_use_area_4vector(bool use_it=true)
Definition: BackgroundEstimator.h:112
fastjet::BackgroundEstimator::_use_area_4vector
bool _use_area_4vector
Definition: BackgroundEstimator.h:153
fastjet::BackgroundEstimator::_median_rho
double _median_rho
background estimated density per unit area
Definition: BackgroundEstimator.h:156
relativeConstraints.error
error
Definition: relativeConstraints.py:53
fastjet::BackgroundEstimator::_csa
const ClusterSequenceAreaBase & _csa
cluster sequence to get jets and areas from
Definition: BackgroundEstimator.h:148
fastjet::BackgroundEstimator::_excluded_jets
std::vector< PseudoJet > _excluded_jets
jets to be excluded
Definition: BackgroundEstimator.h:151
fastjet::BackgroundEstimator::_all_from_inclusive
bool _all_from_inclusive
when true, we'll assume that the incl jets are the complete set
Definition: BackgroundEstimator.h:152
mathSSE::sqrt
T sqrt(T t)
Definition: SSEVec.h:19
fastjet::BackgroundEstimator::_compute
void _compute()
do the actual job
Definition: BackgroundEstimator.cc:55
fastjet::BackgroundEstimator::_included_jets
std::vector< PseudoJet > _included_jets
jets to be used
Definition: BackgroundEstimator.h:150
fastjet::BackgroundEstimator::reset
void reset()
Definition: BackgroundEstimator.cc:34
jetUpdater_cfi.sort
sort
Definition: jetUpdater_cfi.py:29
fastjet::BackgroundEstimator::_range
const RangeDefinition & _range
range to compute the background in
Definition: BackgroundEstimator.h:149
createfilelist.int
int
Definition: createfilelist.py:10
fastjet::BackgroundEstimator::_n_jets_excluded
unsigned int _n_jets_excluded
number of jets that have explicitly been excluded
Definition: BackgroundEstimator.h:160
fastjet::BackgroundEstimator::_n_empty_jets
double _n_empty_jets
number of empty (pure-ghost) jets
Definition: BackgroundEstimator.h:161
res
Definition: Electron.h:6
fastjet::BackgroundEstimator::_empty_area
double _empty_area
the empty (pure-ghost/unclustered) area!
Definition: BackgroundEstimator.h:162
std
Definition: JetResolutionObject.h:76
fastjet::BackgroundEstimator::_uptodate
bool _uptodate
true when the background computation is up-to-date
Definition: BackgroundEstimator.h:165
fastjet::BackgroundEstimator::BackgroundEstimator
BackgroundEstimator(const ClusterSequenceAreaBase &csa, const RangeDefinition &range)
Definition: BackgroundEstimator.cc:24
dqmiolumiharvest.j
j
Definition: dqmiolumiharvest.py:66
fastjet::BackgroundEstimator::~BackgroundEstimator
~BackgroundEstimator()
default dtor
Definition: BackgroundEstimator.cc:30
fastjet::BackgroundEstimator::_mean_area
double _mean_area
mean area of the jets used to estimate the background
Definition: BackgroundEstimator.h:158