FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
Subtractor.cc
1 //FJSTARTHEADER
2 // $Id: Subtractor.cc 3670 2014-09-08 14:17:59Z soyez $
3 //
4 // Copyright (c) 2005-2014, 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 
42 //----------------------------------------------------------------------
43 // ctor
44 Subtractor::Subtractor(double rho) : _bge(0), _rho(rho) {
45  if (_rho<0.0) throw Error("Subtractor(rho) was passed a negative rho value; rho should be >= 0");
46  set_defaults();
47 }
48 
49 //----------------------------------------------------------------------
50 // ctor
51 Subtractor::Subtractor(double rho, double rho_m) : _bge(0), _rho(rho) {
52  if (_rho<0.0) throw Error("Subtractor(rho, rho_m) was passed a negative rho value; rho should be >= 0");
53  if (rho_m<0.0) throw Error("Subtractor(rho, rho_m) was passed a negative rho_m value; rho_m should be >= 0");
54  set_defaults();
55  _rho_m = rho_m;
56  set_use_rho_m(true);
57 }
58 
59 //----------------------------------------------------------------------
61  _rho_m = _invalid_rho;
62  _use_rho_m = false; // likely to change in future releases!!
63  _safe_mass = false; // likely to change in future releases!!
64 
67 }
68 
69 //----------------------------------------------------------------------
70 // perform the subtraction of a given jet
72  if (!jet.has_area()){
73  throw Error("Subtractor::result(...): Trying to subtract a jet without area support");
74  }
75 
76  PseudoJet known_lv, known_pu;
77  PseudoJet unknown = jet;
79  // separate the jet constituents in 3 groups:
80  // unknown vertex
81  // known vertex, leading vertex
82  // known vertex, non-leading vertex (PU)
83  vector<PseudoJet> constits_unknown, constits_known;
85  constits_known,
86  constits_unknown);
87  vector<PseudoJet> constits_known_lv, constits_known_pu;
88  _sel_leading_vertex.sift(constits_known,
89  constits_known_lv,
90  constits_known_pu);
91 
92  // For the parts related to the known vertices (LV or PU), we just
93  // sum the 4-momenta. For the unknown part, we assign it the full
94  // jet area.
95  known_lv = (constits_known_lv.size()!=0)
96  ? SelectorIdentity().sum(constits_known_lv) : 0.0*jet;
97  known_pu = (constits_known_pu.size()!=0)
98  ? SelectorIdentity().sum(constits_known_pu) : 0.0*jet;
99  if (constits_unknown.size()==0){
100  // no need for any form of subtraction!
101  PseudoJet subtracted_jet = jet;
102  subtracted_jet.reset_momentum(known_lv);
103  return subtracted_jet;
104  }
105  unknown = jet; // that keeps all info including area
106  unknown.reset_momentum(SelectorIdentity().sum(constits_unknown));
107  } else {
108  known_lv = jet; // ensures correct rap-phi!
109  known_lv *= 0.0;
110  known_pu = known_lv;
111  }
112 
113  // prepare for the subtraction and compute the 4-vector to be
114  // subtracted
115  PseudoJet subtracted_jet = jet;
116  PseudoJet to_subtract = known_pu + _amount_to_subtract(unknown);
117 
118  // sanity check for the transverse momentum
119  if (to_subtract.pt2() < jet.pt2() ) {
120  // this subtraction should retain the jet's structural
121  // information
122  subtracted_jet -= to_subtract;
123  } else {
124  // this sets the jet's momentum while maintaining all of the jet's
125  // structural information
126  subtracted_jet.reset_momentum(known_lv);
127  return subtracted_jet;
128  }
129 
130  // make sure that in the end the pt is at least the one known to
131  // come from the leading vertex
132  if (subtracted_jet.pt2() < known_lv.pt2()){
133  subtracted_jet.reset_momentum(known_lv);
134  return subtracted_jet;
135  }
136 
137  // sanity check for the mass (if needed)
138  if ((_safe_mass) && (subtracted_jet.m2() < known_lv.m2())){
139  // in this case, we keep pt and phi as obtained from the
140  // subtraction above and take rap and m from the part that comes
141  // from the leading vertex (or the original jet if nothing comes
142  // from the leading vertex)
143  subtracted_jet.reset_momentum(PtYPhiM(subtracted_jet.pt(),
144  known_lv.rap(),
145  subtracted_jet.phi(),
146  known_lv.m()));
147  }
148 
149  return subtracted_jet;
150 }
151 
152 //----------------------------------------------------------------------
153 std::string Subtractor::description() const{
154  if (_bge != 0) {
155  string desc = "Subtractor that uses the following background estimator to determine rho: "+_bge->description();
156  if (use_rho_m()) desc += "; including the rho_m correction";
157  if (safe_mass()) desc += "; including mass safety tests";
158  if (_sel_known_vertex.worker()){
159  desc += "; using known vertex selection: "+_sel_known_vertex.description()+" and leading vertex selection: "+_sel_leading_vertex.description();
160  }
161  return desc;
162  } else if (_rho != _invalid_rho) {
163  ostringstream ostr;
164  ostr << "Subtractor that uses a fixed value of rho = " << _rho;
165  if (use_rho_m()) ostr << " and rho_m = " << _rho_m;
166  return ostr.str();
167  } else {
168  return "Uninitialised subtractor";
169  }
170 }
171 
172 //----------------------------------------------------------------------
173 // compute the 4-vector that should be subtracted from the given
174 // jet
176  // the "transverse momentum" part
177  double rho;
178  if (_bge != 0) {
179  rho = _bge->rho(jet);
180  } else if (_rho != _invalid_rho) {
181  rho = _rho;
182  } else {
183  throw Error("Subtractor::_amount_to_subtract(...): default Subtractor does not have any information about the background, needed to perform the subtraction");
184  }
185 
186  PseudoJet area = jet.area_4vector();
187  PseudoJet to_subtract = rho*area;
188 
189  double const rho_m_warning_threshold = 1e-5;
190 
191  // add an optional contribution from the unknown particles masses
192  if (_use_rho_m) {
193  double rho_m;
194 
195  if (_bge != 0) {
196  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");
197  rho_m = _bge->rho_m(jet);
198  } else if (_rho_m != _invalid_rho) {
199  rho_m = _rho_m;
200  } else {
201  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");
202  }
203  to_subtract += rho_m * PseudoJet(0.0, 0.0, area.pz(), area.E());
204  } else if (_bge &&
205  _bge->has_rho_m() &&
206  _bge->rho_m(jet) > rho_m_warning_threshold * rho) {
207  _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");
208  }
209 
210  return to_subtract;
211 }
212 
213 
214 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:361
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:123
void set_defaults()
reset all parameters to default values
Definition: Subtractor.cc:60
Selector _sel_known_vertex
selects the particles with a known vertex origin
Definition: Subtractor.hh:190
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
virtual PseudoJet result(const PseudoJet &jet) const
returns a jet that's subtracted
Definition: Subtractor.cc:71
double _rho
the fixed value of rho and/or rho_m to use if the user has selected that option
Definition: Subtractor.hh:184
virtual bool has_rho_m() const
Returns true if this background estimator has support for determination of rho_m. ...
Subtractor()
default constructor
Definition: Subtractor.hh:77
virtual std::string description() const =0
returns a textual description of the background estimator
virtual std::string description() const
class description
Definition: Subtractor.cc:153
double m2() const
returns the squared invariant mass // like CLHEP
Definition: PseudoJet.hh:148
void warn(const char *warning)
outputs a warning to standard error (or the user's default warning stream if set) ...
bool safe_mass() const
returns whether or not safety tests on the mass are included
Definition: Subtractor.hh:125
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 'use_rho_m' is true, include in the subtraction the correction from rho_m, the purely longitudin...
Definition: Subtractor.hh:102
std::string description() const
returns a textual description of the selector
Definition: Selector.hh:237
bool _use_rho_m
include the rho_m correction
Definition: Subtractor.hh:187
const SharedPtr< SelectorWorker > & worker() const
returns a (reference to) the underlying worker's shared pointer
Definition: Selector.hh:277
virtual double rho_m() const
returns rho_m, the purely longitudinal, particle-mass-induced component of the background density per...
Selector _sel_leading_vertex
amongst the particles with a known vertex origin, select those coming from the leading vertex ...
Definition: Subtractor.hh:192
Class that encodes information about cuts and other selection criteria that can be applied to PseudoJ...
Definition: Selector.hh:149
bool use_rho_m() const
returns whether or not the rho_m component is used
Definition: Subtractor.hh:110
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:108
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
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 bool has_area() const
check if it has a defined area
Definition: PseudoJet.cc:708
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) ...
Definition: PseudoJet.hh:920
double pt() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:139
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition: PseudoJet.cc:580
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:935
PseudoJet _amount_to_subtract(const PseudoJet &jet) const
compute the 4-vector that should be subtracted from the given jet
Definition: Subtractor.cc:175
double pt2() const
returns the squared transverse momentum
Definition: PseudoJet.hh:137
bool _safe_mass
ensures that the subtracted mass is +ve
Definition: Subtractor.hh:188
virtual PseudoJet area_4vector() const
return the jet 4-vector area.
Definition: PseudoJet.cc:732
virtual double rho() const =0
get rho, the background density per unit area