Loading...
Searching...
No Matches
Subtractor.cc
//FJSTARTHEADER
// $Id$
//
// Copyright (c) 2005-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// The algorithms that underlie FastJet have required considerable
// development. They are described in the original FastJet paper,
// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
// FastJet as part of work towards a scientific publication, please
// quote the version you use and include a citation to the manual and
// optionally also to hep-ph/0512210.
//
// FastJet is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//FJENDHEADER
#include "fastjet/tools/Subtractor.hh"
#include <cassert>
#include <sstream>
#include <limits>
using namespace std;
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
const double Subtractor::_invalid_rho = -numeric_limits<double>::infinity();
LimitedWarning Subtractor::_unused_rho_m_warning;
//----------------------------------------------------------------------
// ctor
Subtractor::Subtractor(double rho) : _bge(0), _rho(rho) {
if (_rho<0.0) throw Error("Subtractor(rho) was passed a negative rho value; rho should be >= 0");
set_defaults();
}
//----------------------------------------------------------------------
// ctor
Subtractor::Subtractor(double rho, double rho_m) : _bge(0), _rho(rho) {
if (_rho<0.0) throw Error("Subtractor(rho, rho_m) was passed a negative rho value; rho should be >= 0");
if (rho_m<0.0) throw Error("Subtractor(rho, rho_m) was passed a negative rho_m value; rho_m should be >= 0");
set_defaults();
_rho_m = rho_m;
set_use_rho_m(true);
}
//----------------------------------------------------------------------
void Subtractor::set_defaults(){
_rho_m = _invalid_rho;
_use_rho_m = false; // likely to change in future releases!!
_safe_mass = false; // likely to change in future releases!!
_sel_known_vertex = Selector();
_sel_leading_vertex = Selector();
}
//----------------------------------------------------------------------
// perform the subtraction of a given jet
PseudoJet Subtractor::result(const PseudoJet & jet) const {
if (!jet.has_area()){
throw Error("Subtractor::result(...): Trying to subtract a jet without area support");
}
PseudoJet known_lv, known_pu;
PseudoJet unknown = jet;
if (_sel_known_vertex.worker()){
// separate the jet constituents in 3 groups:
// unknown vertex
// known vertex, leading vertex
// known vertex, non-leading vertex (PU)
vector<PseudoJet> constits_unknown, constits_known;
_sel_known_vertex.sift(jet.constituents(),
constits_known,
constits_unknown);
vector<PseudoJet> constits_known_lv, constits_known_pu;
_sel_leading_vertex.sift(constits_known,
constits_known_lv,
constits_known_pu);
// For the parts related to the known vertices (LV or PU), we just
// sum the 4-momenta. For the unknown part, we assign it the full
// jet area.
known_lv = (constits_known_lv.size()!=0)
? SelectorIdentity().sum(constits_known_lv) : 0.0*jet;
known_pu = (constits_known_pu.size()!=0)
? SelectorIdentity().sum(constits_known_pu) : 0.0*jet;
if (constits_unknown.size()==0){
// no need for any form of subtraction!
PseudoJet subtracted_jet = jet;
subtracted_jet.reset_momentum(known_lv);
return subtracted_jet;
}
unknown = jet; // that keeps all info including area
unknown.reset_momentum(SelectorIdentity().sum(constits_unknown));
} else {
known_lv = jet; // ensures correct rap-phi!
known_lv *= 0.0;
known_pu = known_lv;
}
// prepare for the subtraction and compute the 4-vector to be
// subtracted
PseudoJet subtracted_jet = jet;
PseudoJet to_subtract = known_pu + _amount_to_subtract(unknown);
// sanity check for the transverse momentum
if (to_subtract.pt2() < jet.pt2() ) {
// this subtraction should retain the jet's structural
// information
subtracted_jet -= to_subtract;
} else {
// this sets the jet's momentum while maintaining all of the jet's
// structural information
subtracted_jet.reset_momentum(known_lv);
return subtracted_jet;
}
// make sure that in the end the pt is at least the one known to
// come from the leading vertex
if (subtracted_jet.pt2() < known_lv.pt2()){
subtracted_jet.reset_momentum(known_lv);
return subtracted_jet;
}
// sanity check for the mass (if needed)
if ((_safe_mass) && (subtracted_jet.m2() < known_lv.m2())){
// in this case, we keep pt and phi as obtained from the
// subtraction above and take rap and m from the part that comes
// from the leading vertex (or the original jet if nothing comes
// from the leading vertex)
subtracted_jet.reset_momentum(PtYPhiM(subtracted_jet.pt(),
known_lv.rap(),
subtracted_jet.phi(),
known_lv.m()));
}
return subtracted_jet;
}
//----------------------------------------------------------------------
std::string Subtractor::description() const{
if (_bge != 0) {
string desc = "Subtractor that uses the following background estimator to determine rho: "+_bge->description();
if (use_rho_m()) desc += "; including the rho_m correction";
if (safe_mass()) desc += "; including mass safety tests";
if (_sel_known_vertex.worker()){
desc += "; using known vertex selection: "+_sel_known_vertex.description()+" and leading vertex selection: "+_sel_leading_vertex.description();
}
return desc;
} else if (_rho != _invalid_rho) {
ostringstream ostr;
ostr << "Subtractor that uses a fixed value of rho = " << _rho;
if (use_rho_m()) ostr << " and rho_m = " << _rho_m;
return ostr.str();
} else {
return "Uninitialised subtractor";
}
}
//----------------------------------------------------------------------
// compute the 4-vector that should be subtracted from the given
// jet
PseudoJet Subtractor::_amount_to_subtract(const PseudoJet &jet) const{
// the "transverse momentum" part
BackgroundEstimate bg_estimate;
double rho;
if (_bge != 0) {
bg_estimate = _bge->estimate(jet);
rho = bg_estimate.rho();
} else if (_rho != _invalid_rho) {
rho = _rho;
} else {
throw Error("Subtractor::_amount_to_subtract(...): default Subtractor does not have any information about the background, needed to perform the subtraction");
}
PseudoJet area = jet.area_4vector();
PseudoJet to_subtract = rho*area;
double const rho_m_warning_threshold = 1e-5;
// add an optional contribution from the unknown particles masses
if (_use_rho_m) {
double rho_m;
if (_bge != 0) {
//if (!_bge->has_rho_m()) throw Error("Subtractor::_amount_to_subtract(...): requested subtraction with rho_m from a background estimator, but the estimator does not have rho_m support");
if (!bg_estimate.has_rho_m()) throw Error("Subtractor::_amount_to_subtract(...): requested subtraction with rho_m from a background estimator, but the estimator does not have rho_m support");
rho_m = bg_estimate.rho_m();
} else if (_rho_m != _invalid_rho) {
rho_m = _rho_m;
} else {
throw Error("Subtractor::_amount_to_subtract(...): default Subtractor does not have any information about the background rho_m, needed to perform the rho_m subtraction");
}
to_subtract += rho_m * PseudoJet(0.0, 0.0, area.pz(), area.E());
} else if (_bge &&
bg_estimate.has_rho_m() &&
bg_estimate.rho_m() > rho_m_warning_threshold * rho) {
_unused_rho_m_warning.warn("Subtractor::_amount_to_subtract(...): Background estimator indicates non-zero rho_m, but use_rho_m()==false in subtractor; consider calling set_use_rho_m(true) to include the rho_m information");
}
return to_subtract;
}
FASTJET_END_NAMESPACE
fastjet::BackgroundEstimate
/// a class that holds the result of the calculation
Definition BackgroundEstimatorBase.hh:53
fastjet::BackgroundEstimate::has_rho_m
bool has_rho_m() const
true if this background estimate has a determination of rho_m.
Definition BackgroundEstimatorBase.hh:86
fastjet::BackgroundEstimate::rho
double rho() const
background density per unit area
Definition BackgroundEstimatorBase.hh:66
fastjet::BackgroundEstimate::rho_m
double rho_m() const
purely longitudinal (particle-mass-induced) component of the background density per unit area
Definition BackgroundEstimatorBase.hh:78
fastjet::Error
base class corresponding to errors that can be thrown by FastJet
Definition Error.hh:52
fastjet::PseudoJet
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition PseudoJet.hh:68
fastjet::PseudoJet::constituents
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition PseudoJet.cc:692
fastjet::PseudoJet::rap
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition PseudoJet.hh:138
fastjet::PseudoJet::reset_momentum
void reset_momentum(double px, double py, double pz, double E)
reset the 4-momentum according to the supplied components but leave all other information (indices,...
Definition PseudoJet.hh:1076
fastjet::PseudoJet::phi
double phi() const
returns phi (in the range 0..2pi)
Definition PseudoJet.hh:123
fastjet::PseudoJet::has_area
virtual bool has_area() const
check if it has a defined area
Definition PseudoJet.cc:820
fastjet::PseudoJet::m2
double m2() const
returns the squared invariant mass // like CLHEP
Definition PseudoJet.hh:163
fastjet::PseudoJet::pt
double pt() const
returns the scalar transverse momentum
Definition PseudoJet.hh:154
fastjet::PseudoJet::m
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP)
Definition PseudoJet.hh:1060
fastjet::PseudoJet::area_4vector
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
Definition PseudoJet.cc:844
fastjet::PseudoJet::pt2
double pt2() const
returns the squared transverse momentum
Definition PseudoJet.hh:152
fastjet::Selector
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition Selector.hh:149
fastjet::Subtractor::use_rho_m
bool use_rho_m() const
returns whether or not the rho_m component is used
Definition Subtractor.hh:115
fastjet::Subtractor::_amount_to_subtract
PseudoJet _amount_to_subtract(const PseudoJet &jet) const
compute the 4-vector that should be subtracted from the given jet
Definition Subtractor.cc:176
fastjet::Subtractor::Subtractor
Subtractor(const BackgroundEstimatorBase *bge)
define a subtractor based on a BackgroundEstimator
Definition Subtractor.hh:65
fastjet::Subtractor::Subtractor
Subtractor()
default constructor
Definition Subtractor.hh:77
fastjet::Subtractor::_invalid_rho
static const double _invalid_rho
a value of rho that is used as a default to label that the stored rho is not valid for subtraction.
Definition Subtractor.hh:208
fastjet::Subtractor::_rho
double _rho
the fixed value of rho and/or rho_m to use if the user has selected that option
Definition Subtractor.hh:189
fastjet::Subtractor::description
virtual std::string description() const
class description
Definition Subtractor.cc:154
fastjet::Subtractor::_safe_mass
bool _safe_mass
ensures that the subtracted mass is +ve
Definition Subtractor.hh:193
fastjet::Subtractor::safe_mass
bool safe_mass() const
returns whether or not safety tests on the mass are included
Definition Subtractor.hh:130
fastjet::Subtractor::set_use_rho_m
void set_use_rho_m(bool use_rho_m_in=true)
when 'use_rho_m' is true, include in the subtraction the correction from rho_m, the purely longitudin...
Definition Subtractor.hh:107
fastjet::Subtractor::_bge
const BackgroundEstimatorBase * _bge
the tool used to estimate the background if has to be mutable in case its underlying selector takes a...
Definition Subtractor.hh:187
fastjet::Subtractor::result
virtual PseudoJet result(const PseudoJet &jet) const
returns a jet that's subtracted
Definition Subtractor.cc:72
fastjet::Subtractor::_use_rho_m
bool _use_rho_m
include the rho_m correction
Definition Subtractor.hh:192
fastjet::Subtractor::_sel_leading_vertex
Selector _sel_leading_vertex
amongst the particles with a known vertex origin, select those coming from the leading vertex
Definition Subtractor.hh:197
fastjet::Subtractor::_sel_known_vertex
Selector _sel_known_vertex
selects the particles with a known vertex origin
Definition Subtractor.hh:195
fastjet::Subtractor::set_defaults
void set_defaults()
reset all parameters to default values
Definition Subtractor.cc:61
fastjet::PtYPhiM
PseudoJet PtYPhiM(double pt, double y, double phi, double m)
return a pseudojet with the given pt, y, phi and mass
Definition PseudoJet.cc:469