29 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
30 #include <fastjet/ClusterSequenceArea.hh>
31 #include <fastjet/ClusterSequenceStructure.hh>
34 FASTJET_BEGIN_NAMESPACE
38 double BackgroundJetScalarPtDensity::result(
const PseudoJet & jet)
const {
39 std::vector<PseudoJet> constituents = jet.
constituents();
41 for (
unsigned i = 0; i < constituents.size(); i++) {
42 scalar_pt += pow(constituents[i].perp(), _pt_power);
44 return scalar_pt / jet.
area();
49 double BackgroundRescalingYPolynomial::result(
const PseudoJet & jet)
const {
52 double rescaling = _a0 + _a1*y + _a2*y2 + _a3*y2*y + _a4*y2*y2;
59 LimitedWarning JetMedianBackgroundEstimator::_warnings_preliminary;
71 JetMedianBackgroundEstimator::JetMedianBackgroundEstimator(
const Selector &rho_range,
74 : _rho_range(rho_range), _jet_def(jet_def), _area_def(area_def) {
80 _check_jet_alg_good_for_median();
107 throw Error(
"JetMedianBackgroundEstimator::set_particles can only be called if you set the jet (and area) definition explicitly through the class constructor");
153 _check_jet_alg_good_for_median();
157 throw Error(
"JetMedianBackgroundEstimator: either an area with explicit ghosts (recommended) or a Selector with finite area is needed (to allow for the computation of the empty area)");
174 throw Error(
"JetMedianBackgroundEstimator::JetMedianBackgroundEstimator: At least one jet is needed to compute the background properties");
179 if (! (jets[0].has_associated_cluster_sequence()) && (jets[0].has_area()))
180 throw Error(
"JetMedianBackgroundEstimator::JetMedianBackgroundEstimator: the jets used to estimate the background properties must be associated with a valid ClusterSequenceAreaBase");
182 _csi = jets[0].structure_shared_ptr();
186 for (
unsigned int i=1;i<jets.size(); i++){
187 if (! jets[i].has_associated_cluster_sequence())
188 throw Error(
"JetMedianBackgroundEstimator::set_jets(...): the jets used to estimate the background properties must be associated with a valid ClusterSequenceAreaBase");
190 if (jets[i].structure_shared_ptr().get() != _csi.
get())
191 throw Error(
"JetMedianBackgroundEstimator::set_jets(...): all the jets used to estimate the background properties must share the same ClusterSequence");
195 _check_jet_alg_good_for_median();
199 throw Error(
"JetMedianBackgroundEstimator: either an area with explicit ghosts (recommended) or a Selector with finite area is needed (to allow for the computation of the empty area)");
204 _included_jets = jets;
218 throw Error(
"The background estimation is obtained from a selector that takes a reference jet. rho(PseudoJet) should be used in that case");
219 _recompute_if_needed();
226 throw Error(
"The background estimation is obtained from a selector that takes a reference jet. rho(PseudoJet) should be used in that case");
227 _recompute_if_needed();
238 _recompute_if_needed(jet);
239 double our_rho = _rho;
240 if (_rescaling_class != 0) {
241 our_rho *= (*_rescaling_class)(jet);
253 _recompute_if_needed(jet);
254 double our_sigma = _sigma;
255 if (_rescaling_class != 0) {
256 our_sigma *= (*_rescaling_class)(jet);
275 _n_jets_used = _n_empty_jets = 0;
276 _empty_area = _mean_area = 0.0;
278 _included_jets.clear();
280 _jet_density_class = 0;
281 _rescaling_class = 0;
291 _warnings_preliminary.
warn(
"JetMedianBackgroundEstimator::set_jet_density_class: density classes are still preliminary in FastJet 3.0. Their interface may differ in future releases (without guaranteeing backward compatibility).");
292 _jet_density_class = jet_density_class_in;
303 desc <<
"JetMedianBackgroundEstimator, using " << _jet_def.
description()
305 <<
" and selecting jets with " << _rho_range.
description();
317 void JetMedianBackgroundEstimator::_recompute_if_needed(
const PseudoJet &jet){
322 if (jet == _current_reference)
return;
330 _recompute_if_needed();
334 void JetMedianBackgroundEstimator::_compute()
const {
340 vector<double> vector_for_median;
341 double total_area = 0.0;
345 vector<PseudoJet> selected_jets = _rho_range(_included_jets);
348 for (
unsigned i = 0; i < selected_jets.size(); i++) {
349 const PseudoJet & current_jet = selected_jets[i];
351 double this_area = (_use_area_4vector) ? current_jet.area_4vector().perp() : current_jet.area();
355 if (_jet_density_class == 0) {
356 median_input = current_jet.perp()/this_area;
358 median_input = (*_jet_density_class)(current_jet);
360 if (_rescaling_class != 0) {
361 median_input /= (*_rescaling_class)(current_jet);
363 vector_for_median.push_back(median_input);
364 total_area += this_area;
367 _warnings_zero_area.
warn(
"JetMedianBackgroundEstimator::_compute(...): 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).");
373 if (vector_for_median.size() == 0) {
381 const ClusterSequenceAreaBase * csab = (
dynamic_cast<ClusterSequenceStructure*
>(_csi()))->validated_csab();
382 if (csab->has_explicit_ghosts()) {
386 _empty_area = csab->empty_area(_rho_range);
387 _n_empty_jets = csab->n_empty_jets(_rho_range);
390 double total_njets = _n_jets_used + _n_empty_jets;
391 total_area += _empty_area;
398 _mean_area = total_area / total_njets;
399 _sigma = stand_dev * sqrt(_mean_area);
409 void JetMedianBackgroundEstimator::_check_csa_alive()
const{
410 ClusterSequenceStructure* csa =
dynamic_cast<ClusterSequenceStructure*
>(_csi());
412 throw Error(
"JetMedianBackgroundEstimator: there is no cluster sequence associated with the JetMedianBackgroundEstimator");
414 if (! dynamic_cast<ClusterSequenceStructure*>(_csi())->has_associated_cluster_sequence())
415 throw Error(
"JetMedianBackgroundEstimator: modifications are no longer possible as the underlying ClusterSequence has gone out of scope");
422 void JetMedianBackgroundEstimator::_check_jet_alg_good_for_median()
const{
423 const JetDefinition * jet_def = &_jet_def;
428 const ClusterSequence * cs =
dynamic_cast<ClusterSequenceStructure*
>(_csi())->validated_cs();
429 jet_def = &(cs->jet_def());
435 _warnings.
warn(
"JetMedianBackgroundEstimator: jet_def being used may not be suitable for estimating diffuse backgrounds (good alternatives are kt, cam)");
442 FASTJET_END_NAMESPACE