FastJet  3.3.3
Subtractor.cc
1 //FJSTARTHEADER
2 // $Id: Subtractor.cc 4420 2019-11-29 09:28:20Z soyez $
3 //
4 // Copyright (c) 2005-2019, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet.
8 //
9 // FastJet is free software; you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation; either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // The algorithms that underlie FastJet have required considerable
15 // development. They are described in the original FastJet paper,
16 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17 // FastJet as part of work towards a scientific publication, please
18 // quote the version you use and include a citation to the manual and
19 // optionally also to hep-ph/0512210.
20 //
21 // FastJet is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28 //----------------------------------------------------------------------
29 //FJENDHEADER
30 
31 #include "fastjet/tools/Subtractor.hh"
32 #include <cassert>
33 #include <sstream>
34 #include <limits>
35 using namespace std;
36 
37 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
38 
39 const double Subtractor::_invalid_rho = -numeric_limits<double>::infinity();
40 
41 LimitedWarning Subtractor::_unused_rho_m_warning;
42 
43 //----------------------------------------------------------------------
44 // ctor
45 Subtractor::Subtractor(double rho) : _bge(0), _rho(rho) {
46  if (_rho<0.0) throw Error("Subtractor(rho) was passed a negative rho value; rho should be >= 0");
47  set_defaults();
48 }
49 
50 //----------------------------------------------------------------------
51 // ctor
52 Subtractor::Subtractor(double rho, double rho_m) : _bge(0), _rho(rho) {
53  if (_rho<0.0) throw Error("Subtractor(rho, rho_m) was passed a negative rho value; rho should be >= 0");
54  if (rho_m<0.0) throw Error("Subtractor(rho, rho_m) was passed a negative rho_m value; rho_m should be >= 0");
55  set_defaults();
56  _rho_m = rho_m;
57  set_use_rho_m(true);
58 }
59 
60 //----------------------------------------------------------------------
62  _rho_m = _invalid_rho;
63  _use_rho_m = false; // likely to change in future releases!!
64  _safe_mass = false; // likely to change in future releases!!
65 
68 }
69 
70 //----------------------------------------------------------------------
71 // perform the subtraction of a given jet
73  if (!jet.has_area()){
74  throw Error("Subtractor::result(...): Trying to subtract a jet without area support");
75  }
76 
77  PseudoJet known_lv, known_pu;
78  PseudoJet unknown = jet;
80  // separate the jet constituents in 3 groups:
81  // unknown vertex
82  // known vertex, leading vertex
83  // known vertex, non-leading vertex (PU)
84  vector<PseudoJet> constits_unknown, constits_known;
86  constits_known,
87  constits_unknown);
88  vector<PseudoJet> constits_known_lv, constits_known_pu;
89  _sel_leading_vertex.sift(constits_known,
90  constits_known_lv,
91  constits_known_pu);
92 
93  // For the parts related to the known vertices (LV or PU), we just
94  // sum the 4-momenta. For the unknown part, we assign it the full
95  // jet area.
96  known_lv = (constits_known_lv.size()!=0)
97  ? SelectorIdentity().sum(constits_known_lv) : 0.0*jet;
98  known_pu = (constits_known_pu.size()!=0)
99  ? SelectorIdentity().sum(constits_known_pu) : 0.0*jet;
100  if (constits_unknown.size()==0){
101  // no need for any form of subtraction!
102  PseudoJet subtracted_jet = jet;
103  subtracted_jet.reset_momentum(known_lv);
104  return subtracted_jet;
105  }
106  unknown = jet; // that keeps all info including area
107  unknown.reset_momentum(SelectorIdentity().sum(constits_unknown));
108  } else {
109  known_lv = jet; // ensures correct rap-phi!
110  known_lv *= 0.0;
111  known_pu = known_lv;
112  }
113 
114  // prepare for the subtraction and compute the 4-vector to be
115  // subtracted
116  PseudoJet subtracted_jet = jet;
117  PseudoJet to_subtract = known_pu + _amount_to_subtract(unknown);
118 
119  // sanity check for the transverse momentum
120  if (to_subtract.pt2() < jet.pt2() ) {
121  // this subtraction should retain the jet's structural
122  // information
123  subtracted_jet -= to_subtract;
124  } else {
125  // this sets the jet's momentum while maintaining all of the jet's
126  // structural information
127  subtracted_jet.reset_momentum(known_lv);
128  return subtracted_jet;
129  }
130 
131  // make sure that in the end the pt is at least the one known to
132  // come from the leading vertex
133  if (subtracted_jet.pt2() < known_lv.pt2()){
134  subtracted_jet.reset_momentum(known_lv);
135  return subtracted_jet;
136  }
137 
138  // sanity check for the mass (if needed)
139  if ((_safe_mass) && (subtracted_jet.m2() < known_lv.m2())){
140  // in this case, we keep pt and phi as obtained from the
141  // subtraction above and take rap and m from the part that comes
142  // from the leading vertex (or the original jet if nothing comes
143  // from the leading vertex)
144  subtracted_jet.reset_momentum(PtYPhiM(subtracted_jet.pt(),
145  known_lv.rap(),
146  subtracted_jet.phi(),
147  known_lv.m()));
148  }
149 
150  return subtracted_jet;
151 }
152 
153 //----------------------------------------------------------------------
154 std::string Subtractor::description() const{
155  if (_bge != 0) {
156  string desc = "Subtractor that uses the following background estimator to determine rho: "+_bge->description();
157  if (use_rho_m()) desc += "; including the rho_m correction";
158  if (safe_mass()) desc += "; including mass safety tests";
159  if (_sel_known_vertex.worker()){
160  desc += "; using known vertex selection: "+_sel_known_vertex.description()+" and leading vertex selection: "+_sel_leading_vertex.description();
161  }
162  return desc;
163  } else if (_rho != _invalid_rho) {
164  ostringstream ostr;
165  ostr << "Subtractor that uses a fixed value of rho = " << _rho;
166  if (use_rho_m()) ostr << " and rho_m = " << _rho_m;
167  return ostr.str();
168  } else {
169  return "Uninitialised subtractor";
170  }
171 }
172 
173 //----------------------------------------------------------------------
174 // compute the 4-vector that should be subtracted from the given
175 // jet
177  // the "transverse momentum" part
178  double rho;
179  if (_bge != 0) {
180  rho = _bge->rho(jet);
181  } else if (_rho != _invalid_rho) {
182  rho = _rho;
183  } else {
184  throw Error("Subtractor::_amount_to_subtract(...): default Subtractor does not have any information about the background, needed to perform the subtraction");
185  }
186 
187  PseudoJet area = jet.area_4vector();
188  PseudoJet to_subtract = rho*area;
189 
190  double const rho_m_warning_threshold = 1e-5;
191 
192  // add an optional contribution from the unknown particles masses
193  if (_use_rho_m) {
194  double rho_m;
195 
196  if (_bge != 0) {
197  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");
198  rho_m = _bge->rho_m(jet);
199  } else if (_rho_m != _invalid_rho) {
200  rho_m = _rho_m;
201  } else {
202  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");
203  }
204  to_subtract += rho_m * PseudoJet(0.0, 0.0, area.pz(), area.E());
205  } else if (_bge &&
206  _bge->has_rho_m() &&
207  _bge->rho_m(jet) > rho_m_warning_threshold * rho) {
208  _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");
209  }
210 
211  return to_subtract;
212 }
213 
214 
215 FASTJET_END_NAMESPACE
PseudoJet PtYPhiM(double pt, double y, double phi, double m)
return a pseudojet with the given pt, y, phi and mass
Definition: PseudoJet.cc:365
double m2() const
returns the squared invariant mass // like CLHEP
Definition: PseudoJet.hh:150
void set_defaults()
reset all parameters to default values
Definition: Subtractor.cc:61
Selector _sel_known_vertex
selects the particles with a known vertex origin
Definition: Subtractor.hh:190
virtual bool has_rho_m() const
Returns true if this background estimator has support for determination of rho_m. ...
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
Definition: PseudoJet.cc:736
virtual bool has_area() const
check if it has a defined area
Definition: PseudoJet.cc:712
const SharedPtr< SelectorWorker > & worker() const
returns a (reference to) the underlying worker&#39;s shared pointer
Definition: Selector.hh:277
double _rho
the fixed value of rho and/or rho_m to use if the user has selected that option
Definition: Subtractor.hh:184
Subtractor()
default constructor
Definition: Subtractor.hh:77
virtual std::string description() const =0
returns a textual description of the background estimator
virtual double rho_m() const
returns rho_m, the purely longitudinal, particle-mass-induced component of the background density per...
double pt() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:141
PseudoJet _amount_to_subtract(const PseudoJet &jet) const
compute the 4-vector that should be subtracted from the given jet
Definition: Subtractor.cc:176
void warn(const char *warning)
outputs a warning to standard error (or the user&#39;s default warning stream if set) ...
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition: PseudoJet.cc:584
bool use_rho_m() const
returns whether or not the rho_m component is used
Definition: Subtractor.hh:110
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
void set_use_rho_m(bool use_rho_m_in=true)
when &#39;use_rho_m&#39; is true, include in the subtraction the correction from rho_m, the purely longitudin...
Definition: Subtractor.hh:102
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) ...
Definition: PseudoJet.hh:970
bool _use_rho_m
include the rho_m correction
Definition: Subtractor.hh:187
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:110
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:125
bool safe_mass() const
returns whether or not safety tests on the mass are included
Definition: Subtractor.hh:125
Selector _sel_leading_vertex
amongst the particles with a known vertex origin, select those coming from the leading vertex ...
Definition: Subtractor.hh:192
void sift(const std::vector< PseudoJet > &jets, std::vector< PseudoJet > &jets_that_pass, std::vector< PseudoJet > &jets_that_fail) const
sift the input jets into two vectors – those that pass the selector and those that do not ...
Definition: Selector.cc:153
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:149
virtual std::string description() const
class description
Definition: Subtractor.cc:154
std::string description() const
returns a textual description of the selector
Definition: Selector.hh:237
BackgroundEstimatorBase * _bge
the tool used to estimate the background if has to be mutable in case its underlying selector takes a...
Definition: Subtractor.hh:182
double pt2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:139
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:203
PseudoJet sum(const std::vector< PseudoJet > &jets) const
Return the 4-vector sum of the objects that pass the selection.
Definition: Selector.cc:101
virtual PseudoJet result(const PseudoJet &jet) const
returns a jet that&#39;s subtracted
Definition: Subtractor.cc:72
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
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:985
bool _safe_mass
ensures that the subtracted mass is +ve
Definition: Subtractor.hh:188
virtual double rho() const =0
get rho, the background density per unit area