31 #include "fastjet/tools/JetMedianBackgroundEstimator.hh"
32 #include "fastjet/ClusterSequenceArea.hh"
33 #include "fastjet/ClusterSequenceStructure.hh"
39 FASTJET_BEGIN_NAMESPACE
41 double BackgroundJetScalarPtDensity::result(
const PseudoJet & jet)
const {
46 for (
unsigned i = 0; i < constituents.size(); i++) {
47 scalar_pt += pow(constituents[i].perp(), _pt_power);
49 return scalar_pt / jet.
area();
53 std::string BackgroundJetScalarPtDensity::description()
const {
55 oss <<
"BackgroundScalarJetPtDensity";
56 if (_pt_power != 1.0) oss <<
" with pt_power = " << _pt_power;
62 double BackgroundRescalingYPolynomial::result(
const PseudoJet & jet)
const {
65 double rescaling = _a0 + _a1*y + _a2*y2 + _a3*y2*y + _a4*y2*y2;
72 LimitedWarning JetMedianBackgroundEstimator::_warnings_preliminary;
84 JetMedianBackgroundEstimator::JetMedianBackgroundEstimator(
const Selector &rho_range,
87 : _rho_range(rho_range), _jet_def(jet_def), _area_def(area_def){
93 _check_jet_alg_good_for_median();
131 throw Error(
"JetMedianBackgroundEstimator::set_particles can only be called if you set the jet (and area) definition explicitly through the class constructor");
145 if (seed.size() == 0) {
160 _set_cache_unavailable();
188 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)");
198 _check_jet_alg_good_for_median();
204 _set_cache_unavailable();
217 throw Error(
"JetMedianBackgroundEstimator::JetMedianBackgroundEstimator: At least one jet is needed to compute the background properties");
222 if (! (jets[0].has_associated_cluster_sequence()) && (jets[0].has_area()))
223 throw Error(
"JetMedianBackgroundEstimator::JetMedianBackgroundEstimator: the jets used to estimate the background properties must be associated with a valid ClusterSequenceAreaBase");
229 for (
unsigned int i=1;i<jets.size(); i++){
230 if (! jets[i].has_associated_cluster_sequence())
231 throw Error(
"JetMedianBackgroundEstimator::set_jets(...): the jets used to estimate the background properties must be associated with a valid ClusterSequenceAreaBase");
233 if (jets[i].structure_shared_ptr().get() != csi_shared.
get())
234 throw Error(
"JetMedianBackgroundEstimator::set_jets(...): all the jets used to estimate the background properties must share the same ClusterSequence");
239 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)");
249 _check_jet_alg_good_for_median();
252 _included_jets = jets;
255 _set_cache_unavailable();
272 throw Error(
"The background estimation is obtained from a selector that takes a reference jet. estimate(PseudoJet) should be used in that case");
274 if (!_cache_available) _compute_and_cache_no_overwrite();
282 double rescaling_factor = (_rescaling_class != 0)
283 ? (*_rescaling_class)(jet) : 1.0;
289 local_estimate = _compute(jet);
293 if (!_cache_available) _compute_and_cache_no_overwrite();
297 return local_estimate;
305 throw Error(
"The background estimation is obtained from a selector that takes a reference jet. rho(PseudoJet) should be used in that case");
309 if (!_cache_available) _compute_and_cache_no_overwrite();
316 throw Error(
"The background estimation is obtained from a selector that takes a reference jet. sigma(PseudoJet) should be used in that case");
317 if (!_cache_available) _compute_and_cache_no_overwrite();
329 double rescaling_factor = (_rescaling_class != 0)
330 ? (*_rescaling_class)(jet) : 1.0;
341 if (!_cache_available) _compute_and_cache_no_overwrite();
353 double rescaling_factor = (_rescaling_class != 0)
354 ? (*_rescaling_class)(jet) : 1.0;
362 if (!_cache_available) _compute_and_cache_no_overwrite();
371 throw Error(
"JetMediamBackgroundEstimator: rho_m requested but rho_m calculation is disabled (either eplicitly or due to the presence of a jet density class).");
374 throw Error(
"The background estimation is obtained from a selector that takes a reference jet. rho_m(PseudoJet) should be used in that case");
376 if (!_cache_available) _compute_and_cache_no_overwrite();
387 throw Error(
"JetMediamBackgroundEstimator: sigma_m requested but rho_m/sigma_m calculation is disabled (either explicitly or due to the presence of a jet density class).");
390 throw Error(
"The background estimation is obtained from a selector that takes a reference jet. sigma_m(PseudoJet) should be used in that case");
391 if (!_cache_available) _compute_and_cache_no_overwrite();
400 double rescaling_factor = (_rescaling_class != 0)
401 ? (*_rescaling_class)(jet) : 1.0;
409 if (!_cache_available) _compute_and_cache_no_overwrite();
419 double rescaling_factor = (_rescaling_class != 0)
420 ? (*_rescaling_class)(jet) : 1.0;
428 if (!_cache_available) _compute_and_cache_no_overwrite();
442 if (!_cache_available){
444 throw Error(
"Calls to JetMedianBackgroundEstimator::mean_area() in cases where the background estimation uses a selector that takes a reference jet need to call a method that fills the cached estimate (rho(jet), sigma(jet), ...).");
450 if (!_cache_available) _compute_and_cache_no_overwrite();
462 if (!_cache_available){
464 throw Error(
"Calls to JetMedianBackgroundEstimator::n_jets_used() in cases where the background estimation uses a selector that takes a reference jet need to call a method that fills the cached estimate (rho(jet), sigma(jet), ...).");
470 if (!_cache_available) _compute_and_cache_no_overwrite();
477 vector<PseudoJet> tmp_jets;
483 if (!_cache_available){
485 throw Error(
"Calls to JetMedianBackgroundEstimator::jets_used() in cases where the background estimation uses a selector that takes a reference jet need to call a method that fills the cached estimate (rho(jet), sigma(jet), ...).");
489 Selector local_rho_range = _rho_range;
491 tmp_jets = local_rho_range(_included_jets);
493 if (!_cache_available) _compute_and_cache_no_overwrite();
494 tmp_jets = _rho_range(_included_jets);
497 std::vector<PseudoJet> used_jets;
498 for (
unsigned int i=0; i<tmp_jets.size(); i++){
499 if (tmp_jets[i].area()>0) used_jets.push_back(tmp_jets[i]);
523 if (!_cache_available){
525 throw Error(
"Calls to JetMedianBackgroundEstimator::empty_area() in cases where the background estimation uses a selector that takes a reference jet need to call a method that fills the cached estimate (rho(jet), sigma(jet), ...).");
531 if (!_cache_available) _compute_and_cache_no_overwrite();
552 if (!_cache_available){
554 throw Error(
"Calls to JetMedianBackgroundEstimator::n_empty_jets() in cases where the background estimation uses a selector that takes a reference jet need to call a method that fills the cached estimate (rho(jet), sigma(jet), ...).");
560 if (!_cache_available) _compute_and_cache_no_overwrite();
576 _enable_rho_m =
true;
579 _included_jets.clear();
581 _jet_density_class = 0;
582 _rescaling_class = 0;
584 _set_cache_unavailable();
592 _warnings_preliminary.
warn(
"JetMedianBackgroundEstimator::set_jet_density_class: density classes are still preliminary in FastJet 3.1. Their interface may differ in future releases (without guaranteeing backward compatibility). Note that since FastJet 3.1, rho_m and sigma_m are accessible direclty in JetMedianBackgroundEstimator and GridMedianBackgroundEstimator(with no need for a density class).");
593 _jet_density_class = jet_density_class_in;
594 _set_cache_unavailable();
604 desc <<
"JetMedianBackgroundEstimator, using " << _jet_def.
description()
606 <<
" and selecting jets with " << _rho_range.
description();
624 local_estimate.set_has_sigma(
has_sigma());
625 local_estimate.set_has_rho_m(
has_rho_m());
629 Extras * extras =
new Extras;
631 extras->set_reference_jet(jet);
635 vector<double> vector_for_median_pt;
636 vector<double> vector_for_median_dt;
637 double total_area = 0.0;
640 vector<PseudoJet> selected_jets;
642 Selector local_rho_range = _rho_range;
643 selected_jets = local_rho_range.
set_reference(jet)(_included_jets);
645 selected_jets = _rho_range(_included_jets);
649 double median_input_pt, median_input_dt=0.0;
650 BackgroundJetPtMDensity m_density;
652 unsigned int njets_used = 0;
653 for (
unsigned i = 0; i < selected_jets.size(); i++) {
654 const PseudoJet & current_jet = selected_jets[i];
656 double this_area = (_use_area_4vector) ? current_jet.area_4vector().perp() : current_jet.area();
660 if (_jet_density_class == 0) {
661 median_input_pt = current_jet.perp()/this_area;
663 median_input_pt = (*_jet_density_class)(current_jet);
670 median_input_dt = m_density(current_jet);
673 if (_rescaling_class != 0) {
674 double resc = (*_rescaling_class)(current_jet);
675 median_input_pt /= resc;
676 median_input_dt /= resc;
680 vector_for_median_pt.push_back(median_input_pt);
682 vector_for_median_dt.push_back(median_input_dt);
684 total_area += this_area;
687 _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).");
692 if (vector_for_median_pt.size() == 0) {
694 return local_estimate;
699 const ClusterSequenceAreaBase * csab = (
dynamic_cast<ClusterSequenceStructure*
>(_csi.get()))->validated_csab();
700 if (! (csab->has_explicit_ghosts())) {
702 Selector local_rho_range = _rho_range;
704 extras->set_empty_area (csab->empty_area(local_rho_range));
705 extras->set_n_empty_jets(csab->n_empty_jets(local_rho_range));
707 extras->set_empty_area (csab->empty_area(_rho_range));
708 extras->set_n_empty_jets(csab->n_empty_jets(_rho_range));
712 extras->set_n_jets_used(njets_used);
713 double total_njets = extras->n_jets_used() + extras->n_empty_jets();
714 total_area += extras->empty_area();
716 double rho_tmp, stand_dev;
718 rho_tmp, stand_dev, _provide_fj2_sigma);
719 local_estimate.set_rho(rho_tmp);
722 local_estimate.set_mean_area(total_area / total_njets);
723 local_estimate.set_sigma(stand_dev * sqrt(local_estimate.
mean_area()));
730 local_estimate.set_rho_m(rho_tmp);
731 local_estimate.set_sigma_m(stand_dev * sqrt(local_estimate.
mean_area()));
734 return local_estimate;
738 void JetMedianBackgroundEstimator::_cache_no_overwrite(
const BackgroundEstimate &estimate)
const {
753 if (!_cache_available){
755 _cache_available =
true;
762 void JetMedianBackgroundEstimator::_compute_and_cache_no_overwrite()
const {
768 _cache_no_overwrite(_compute(PseudoJet()));
771 void JetMedianBackgroundEstimator::_cache(
const BackgroundEstimate &estimate)
const {
785 _cache_available =
true;
791 BackgroundEstimate JetMedianBackgroundEstimator::_compute_and_cache_if_needed(
const PseudoJet &jet)
const {
795 BackgroundEstimate local_estimate;
801 return local_estimate;
805 local_estimate = _compute(jet);
806 _cache(local_estimate);
807 return local_estimate;
812 void JetMedianBackgroundEstimator::_check_csa_alive()
const{
813 ClusterSequenceStructure* csa =
dynamic_cast<ClusterSequenceStructure*
>(_csi.get());
815 throw Error(
"JetMedianBackgroundEstimator: there is no cluster sequence associated with the JetMedianBackgroundEstimator");
817 if (! csa->has_associated_cluster_sequence())
818 throw Error(
"JetMedianBackgroundEstimator: modifications are no longer possible as the underlying ClusterSequence has gone out of scope");
825 void JetMedianBackgroundEstimator::_check_jet_alg_good_for_median()
const{
826 const JetDefinition * jet_def = &_jet_def;
831 const ClusterSequence * cs =
dynamic_cast<ClusterSequenceStructure*
>(_csi.get())->validated_cs();
832 jet_def = &(cs->jet_def());
838 _warnings.
warn(
"JetMedianBackgroundEstimator: jet_def being used may not be suitable for estimating diffuse backgrounds (good alternatives are kt, cam)");
843 FASTJET_END_NAMESPACE