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) {
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{
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{
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
const JetDefinition & jet_def() const
return a reference to the jet definition
a version of cambridge with a special distance measure for particles whose pt is < extra_param() ...
JetAlgorithm jet_algorithm() const
return information about the definition...
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
General class for user to obtain ClusterSequence with additional area information.
the value for the jet algorithm in a JetDefinition for which no algorithm has yet been defined ...
void warn(const std::string &warning)
outputs a warning to standard error (or the user's default warning stream if set) ...
void _median_and_stddev(const std::vector< double > &quantity_vector, double n_empty_jets, double &median, double &stand_dev_if_gaussian, bool do_fj2_calculation=false) const
given a quantity in a vector (e.g.
virtual const ClusterSequenceAreaBase * validated_csab() const
if the jet has valid area information then return a pointer to the associated ClusterSequenceAreaBase...
bool has_finite_area() const
returns true if it has a meaningful and finite area (i.e.
the longitudinally invariant kt algorithm
virtual double area() const
return the jet (scalar) area.
class that holds a generic area definition
virtual double n_empty_jets(const Selector &selector) const
return something similar to the number of pure ghost jets in the given selector's range in an active ...
bool takes_reference() const
returns true if this can be applied jet by jet
Contains any information related to the clustering that should be directly accessible to PseudoJet...
class to provide facilities for giving warnings up to some maximum number of times and to provide glo...
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
double perp() const
returns the scalar transverse momentum
base class that sets interface for extensions of ClusterSequence that provide information about the a...
base class corresponding to errors that can be thrown by FastJet
const SharedPtr< PseudoJetStructureBase > & structure_shared_ptr() const
retrieve a shared pointer to the wrapper to this ClusterSequence
double rap() const
returns the rapidity or some large value when the rapidity is infinite
const Selector & set_reference(const PseudoJet &reference)
set the reference jet for this Selector
std::string description() const
return a description of the current area definition
virtual bool has_explicit_ghosts() const
returns true if ghosts are explicitly included within jets for this ClusterSequence; ...
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
virtual double empty_area(const Selector &selector) const
return the total area, corresponding to the given Selector, that is free of jets, in general based on...
std::string description() const
returns a textual description of the selector
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
std::vector< PseudoJet > inclusive_jets(const double &ptmin=0.0) const
return a vector of all jets (in the sense of the inclusive algorithm) with pt >= ptmin.
std::string description() const
return a textual description of the current jet definition
void delete_self_when_unused()
by calling this routine you tell the ClusterSequence to delete itself when all the Pseudojets associa...
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
class that is intended to hold a full definition of the jet clusterer