Fork me on GitHub

source: git/external/fastjet/tools/Subtractor.cc@ 35cdc46

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 35cdc46 was 35cdc46, checked in by Pavel Demin <demin@…>, 10 years ago

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

  • Property mode set to 100644
File size: 6.9 KB
Line 
1//FJSTARTHEADER
2// $Id: Subtractor.cc 3594 2014-08-12 15:25:11Z 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>
35using namespace std;
36
37FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
38
39const double Subtractor::_invalid_rho = -numeric_limits<double>::infinity();
40
41
42//----------------------------------------------------------------------
43// ctor
44Subtractor::Subtractor(double rho) : _bge(0), _rho(rho) {
45 assert(_rho>0.0);
46 set_defaults();
47}
48
49//----------------------------------------------------------------------
50void Subtractor::set_defaults(){
51 _use_rho_m = false; // likely to change in future releases!!
52 _safe_mass = false; // likely to change in future releases!!
53
54 _sel_known_vertex = Selector();
55 _sel_leading_vertex = Selector();
56}
57
58//----------------------------------------------------------------------
59// perform the subtraction of a given jet
60PseudoJet Subtractor::result(const PseudoJet & jet) const {
61 if (!jet.has_area()){
62 throw Error("Trying to subtract a jet without area support");
63 }
64
65 PseudoJet known_lv, known_pu;
66 PseudoJet unknown = jet;
67 if (_sel_known_vertex.worker()){
68 // separate the jet constituents in 3 groups:
69 // unknown vertex
70 // known vertex, leading vertex
71 // known vertex, non-leading vertex (PU)
72 vector<PseudoJet> constits_unknown, constits_known;
73 _sel_known_vertex.sift(jet.constituents(),
74 constits_known,
75 constits_unknown);
76 vector<PseudoJet> constits_known_lv, constits_known_pu;
77 _sel_leading_vertex.sift(constits_known,
78 constits_known_lv,
79 constits_known_pu);
80
81 // For the parts related to the known vertices (LV or PU), we just
82 // sum the 4-momenta. For the unknown part, we assign it the full
83 // jet area.
84 known_lv = (constits_known_lv.size()!=0)
85 ? SelectorIdentity().sum(constits_known_lv) : 0.0*jet;
86 known_pu = (constits_known_pu.size()!=0)
87 ? SelectorIdentity().sum(constits_known_pu) : 0.0*jet;
88 if (constits_unknown.size()==0){
89 // no need for any form of subtraction!
90 PseudoJet subtracted_jet = jet;
91 subtracted_jet.reset_momentum(known_lv);
92 return subtracted_jet;
93 }
94 unknown = jet; // that keeps all info including area
95 unknown.reset_momentum(SelectorIdentity().sum(constits_unknown));
96 } else {
97 known_lv = jet; // ensures correct rap-phi!
98 known_lv *= 0.0;
99 known_pu = known_lv;
100 }
101
102 // prepare for the subtraction and compute the 4-vector to be
103 // subtracted
104 PseudoJet subtracted_jet = jet;
105 PseudoJet to_subtract = known_pu + _amount_to_subtract(unknown);
106
107 // sanity check for the transverse momentum
108 if (to_subtract.pt2() < jet.pt2() ) {
109 // this subtraction should retain the jet's structural
110 // information
111 subtracted_jet -= to_subtract;
112 } else {
113 // this sets the jet's momentum while maintaining all of the jet's
114 // structural information
115 subtracted_jet.reset_momentum(known_lv);
116 return subtracted_jet;
117 }
118
119 // make sure that in the end the pt is at least the one known to
120 // come from the leading vertex
121 if (subtracted_jet.pt2() < known_lv.pt2()){
122 subtracted_jet.reset_momentum(known_lv);
123 return subtracted_jet;
124 }
125
126 // sanity check for the mass (if needed)
127 if ((_safe_mass) && (subtracted_jet.m2() < known_lv.m2())){
128 // in this case, we keep pt and phi as obtained from the
129 // subtraction above and take rap and m from the part that comes
130 // from the leading vertex (or the original jet if nothing comes
131 // from the leading vertex)
132 subtracted_jet.reset_momentum(PtYPhiM(subtracted_jet.pt(),
133 known_lv.rap(),
134 subtracted_jet.phi(),
135 known_lv.m()));
136 }
137
138 return subtracted_jet;
139}
140
141//----------------------------------------------------------------------
142std::string Subtractor::description() const{
143 if (_bge != 0) {
144 string desc = "Subtractor that uses the following background estimator to determine rho: "+_bge->description();
145 if (use_rho_m()) desc += "; including the rho_m correction";
146 if (safe_mass()) desc += "; including mass safety tests";
147 if (_sel_known_vertex.worker()){
148 desc += "; using known vertex selection: "+_sel_known_vertex.description()+" and leading vertex selection: "+_sel_leading_vertex.description();
149 }
150 return desc;
151 } else if (_rho != _invalid_rho) {
152 ostringstream ostr;
153 ostr << "Subtractor that uses a fixed value of rho = " << _rho;
154 return ostr.str();
155 } else {
156 return "Uninitialised subtractor";
157 }
158}
159
160//----------------------------------------------------------------------
161// compute the 4-vector that should be subtracted from the given
162// jet
163PseudoJet Subtractor::_amount_to_subtract(const PseudoJet &jet) const{
164 // the "transverse momentum" part
165 double rho;
166 if (_bge != 0) {
167 rho = _bge->rho(jet);
168 } else if (_rho != _invalid_rho) {
169 rho = _rho;
170 } else {
171 throw Error("default Subtractor does not have any information about the background, which is needed to perform the subtraction");
172 }
173
174 PseudoJet area = jet.area_4vector();
175 PseudoJet to_subtract = rho*area;
176
177 double const rho_m_warning_threshold = 1e-5;
178
179 // add an optional contribution from the unknown particles masses
180 if (_use_rho_m){
181 assert(_bge != 0); // test done in "set_use_rho_m()"
182 to_subtract += _bge->rho_m(jet) * PseudoJet(0.0, 0.0, area.pz(), area.E());
183 } else if (_bge &&
184 _bge->has_rho_m() &&
185 _bge->rho_m(jet) > rho_m_warning_threshold * rho) {
186 _unused_rho_m_warning.warn("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");
187 }
188
189 return to_subtract;
190}
191
192
193FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.