Fork me on GitHub

source: git/external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.cc@ af2c385

ImprovedOutputFile Timing dual_readout llp
Last change on this file since af2c385 was 1f1f858, checked in by Michele Selvaggi <michele.selvaggi@…>, 10 years ago

added RecursiveTools package from fastjet contribs

  • Property mode set to 100644
File size: 10.1 KB
RevLine 
[1f1f858]1// $Id: RecursiveSymmetryCutBase.cc 700 2014-07-07 12:50:05Z gsoyez $
2//
3// Copyright (c) 2014-, Gavin P. Salam, Gregory Soyez, Jesse Thaler
4//
5//----------------------------------------------------------------------
6// This file is part of FastJet contrib.
7//
8// It is free software; you can redistribute it and/or modify it under
9// the terms of the GNU General Public License as published by the
10// Free Software Foundation; either version 2 of the License, or (at
11// your option) any later version.
12//
13// It is distributed in the hope that it will be useful, but WITHOUT
14// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
16// License for more details.
17//
18// You should have received a copy of the GNU General Public License
19// along with this code. If not, see <http://www.gnu.org/licenses/>.
20//----------------------------------------------------------------------
21
22#include "RecursiveSymmetryCutBase.hh"
23#include "fastjet/JetDefinition.hh"
24#include "fastjet/ClusterSequenceAreaBase.hh"
25#include <algorithm>
26#include <cstdlib>
27
28using namespace std;
29
30FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
31
32namespace contrib{
33
34LimitedWarning RecursiveSymmetryCutBase::_negative_mass_warning;
35LimitedWarning RecursiveSymmetryCutBase::_mu2_gt1_warning;
36//LimitedWarning RecursiveSymmetryCutBase::_nonca_warning;
37LimitedWarning RecursiveSymmetryCutBase::_explicit_ghost_warning;
38
39bool RecursiveSymmetryCutBase::_verbose = false;
40
41//----------------------------------------------------------------------
42PseudoJet RecursiveSymmetryCutBase::result(const PseudoJet & jet) const {
43 // construct the input jet (by default, recluster with C/A)
44
45 if (! jet.has_constituents()){
46 throw Error("RecursiveSymmetryCutBase can only be applied to jets with constituents");
47 }
48
49 PseudoJet j =
50 _do_reclustering
51 ? _recluster ? (*_recluster)(jet)
52 : Recluster(cambridge_algorithm, JetDefinition::max_allowable_R)(jet)
53 : jet;
54
55 // issue a warning if the jet is not obtained through a C/A
56 // clustering
57 // if ((! j.has_associated_cluster_sequence()) ||
58 // (j.validated_cs()->jet_def().jet_algorithm() != cambridge_algorithm))
59 // _nonca_warning.warn("RecursiveSymmetryCutBase is designed to be applied on jets from a Cambridge/Aachen clustering; use it with other algorithms at your own risk.");
60
61 if (! j.has_valid_cluster_sequence()){
62 throw Error("RecursiveSymmetryCutBase can only be applied to jets associated to a (valid) cluster sequence");
63 }
64
65 if (_subtractor) {
66 const ClusterSequenceAreaBase * csab =
67 dynamic_cast<const ClusterSequenceAreaBase *>(j.associated_cs());
68 if (csab == 0 || (!csab->has_explicit_ghosts()))
69 _explicit_ghost_warning.warn("RecursiveSymmetryCutBase: there is no clustering sequence, or it lacks explicit ghosts: subtraction is not guaranteed to function properly");
70 }
71
72 // establish the first subjet and optionally subtract it
73 PseudoJet subjet = j;
74 if (_subtractor && (!_input_jet_is_subtracted)) {
75 subjet = (*_subtractor)(subjet);
76 }
77
78 bool use_mu_cut = (_mu_cut != numeric_limits<double>::infinity());
79
80 // variables for tracking what will happen
81 PseudoJet piece1, piece2;
82
83 // vectors for storing optional verbose structure
84 // these hold the deltaR, symmetry, and mu values of dropped branches
85 std::vector<double> dropped_delta_R;
86 std::vector<double> dropped_symmetry;
87 std::vector<double> dropped_mu;
88
89 // now recurse into the jet's structure
90 while (subjet.has_parents(piece1, piece2)) {
91
92 // first sanity check:
93 // - zero or negative pts are not allowed for the input subjet
94 // - zero or negative masses are not allowed for configurations
95 // in which the mass will effectively appear in a denominator
96 // (The masses will be checked later)
97 if (subjet.pt2() <= 0) return PseudoJet();
98
99 if (_subtractor) {
100 piece1 = (*_subtractor)(piece1);
101 piece2 = (*_subtractor)(piece2);
102 }
103
104 // determine the symmetry parameter
105 double sym;
106
107 if (_symmetry_measure == y) {
108 // the original d_{ij}/m^2 choice from MDT
109 // first make sure the mass is sensible
110 if (subjet.m2() <= 0) {
111 _negative_mass_warning.warn("RecursiveSymmetryCutBase: cannot calculate y, because (sub)jet mass is negative; bailing out");
112 return _result_no_substructure(subjet); //TBC: do we return the hardest parent? A NULL PseudoJet?
113 }
114 sym = piece1.kt_distance(piece2) / subjet.m2();
115
116 } else if (_symmetry_measure == vector_z) {
117 // min(pt1, pt2)/(pt), where the denominator is a vector sum
118 // of the two subjets
119 sym = min(piece1.pt(), piece2.pt()) / subjet.pt();
120
121 } else if (_symmetry_measure == scalar_z) {
122 // min(pt1, pt2)/(pt1+pt2), where the denominator is a scalar sum
123 // of the two subjets
124 double pt1 = piece1.pt();
125 double pt2 = piece2.pt();
126 // make sure denominator is non-zero
127 sym = pt1 + pt2;
128 if (sym == 0) return PseudoJet();
129 sym = min(pt1, pt2) / sym;
130
131 } else {
132 throw Error ("Unrecognized choice of symmetry_measure");
133 }
134
135 // determine the symmetry cut
136 // (This function is specified in the derived classes)
137 double this_symmetry_cut = symmetry_cut_fn(piece1, piece2);
138
139 // and make a first tagging decision based on symmetry cut
140 bool tagged = (sym > this_symmetry_cut);
141
142 // if tagged based on symmetry cut, then check the mu cut (if relevant)
143 // and update the tagging decision. Calculate mu^2 regardless, for cases
144 // of users not cutting on mu2, but still interested in its value.
145 double mu2 = max(piece1.m2(), piece2.m2())/subjet.m2();
146 if (tagged && use_mu_cut) {
147 // first a sanity check -- mu2 won't be sensible if the subjet mass
148 // is negative, so we can't then trust the mu cut - bail out
149 if (subjet.m2() <= 0) {
150 _negative_mass_warning.warn("RecursiveSymmetryCutBase: cannot trust mu, because (sub)jet mass is negative; bailing out");
151 return PseudoJet();
152 }
153 if (mu2 > 1) _mu2_gt1_warning.warn("RecursiveSymmetryCutBase encountered mu^2 value > 1");
154 if (mu2 > pow(_mu_cut,2)) tagged = false;
155 }
156
157
158 // if we've tagged the splitting, return the jet with its substructure
159 if (tagged) {
160 // record relevant information
161 StructureType * structure = new StructureType(subjet);
162 structure->_symmetry = sym;
163 structure->_mu = (mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2);
164 structure->_delta_R = piece1.delta_R(piece2);
165 if (_verbose_structure) {
166 structure->_has_verbose = true;
167 structure->_dropped_symmetry = dropped_symmetry;
168 structure->_dropped_mu = dropped_mu;
169 structure->_dropped_delta_R = dropped_delta_R;
170 }
171 subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
172 return subjet;
173 }
174
175 // if desired, store information about dropped branches before recursing
176 if (_verbose_structure) {
177 dropped_delta_R.push_back(piece1.delta_R(piece2));
178 dropped_symmetry.push_back(sym);
179 dropped_mu.push_back((mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2));
180 }
181
182 // otherwise continue unclustering, allowing for the different
183 // ways of choosing which parent to look into
184 int choice;
185 if (_recursion_choice == larger_mt) {
186 choice = piece1.mt2() > piece2.mt2() ? 1 : 2;
187
188 } else if (_recursion_choice == larger_pt) {
189 choice = piece1.pt2() > piece2.pt2() ? 1 : 2;
190
191 } else if (_recursion_choice == larger_m) {
192 choice = piece1.m2() > piece2.m2() ? 1 : 2;
193
194 } else {
195 throw Error ("Unrecognized value for recursion_choice");
196 }
197 if (_verbose) cout << "choice is " << choice << endl;;
198 subjet = (choice == 1) ? piece1 : piece2;
199 } // (subjet.has_parents(...))
200
201 if (_verbose) cout << "reached end; returning null jet " << endl;
202
203 // decide on tagging versus grooming mode here
204 PseudoJet result = _result_no_substructure(subjet);
205
206 if (result != 0) {
207 // if in grooming mode, add dummy structure information
208 StructureType * structure = new StructureType(result);
209 structure->_symmetry = 0.0;
210 structure->_mu = 0.0;
211 structure->_delta_R = 0.0;
212 if (_verbose_structure) { // still want to store verbose information about dropped branches
213 structure->_has_verbose = true;
214 structure->_dropped_symmetry = dropped_symmetry;
215 structure->_dropped_mu = dropped_mu;
216 structure->_dropped_delta_R = dropped_delta_R;
217 }
218 result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
219 }
220
221 return result;
222}
223
224
225//----------------------------------------------------------------------
226string RecursiveSymmetryCutBase::description() const {
227 ostringstream ostr;
228 ostr << "Recursive " << (_grooming_mode ? "Groomer" : "Tagger") << " with a symmetry cut ";
229
230 switch(_symmetry_measure) {
231 case y:
232 ostr << "y"; break;
233 case scalar_z:
234 ostr << "scalar_z"; break;
235 case vector_z:
236 ostr << "vector_z"; break;
237 default:
238 cerr << "failed to interpret symmetry_measure" << endl; exit(-1);
239 }
240 ostr << " > " << symmetry_cut_description();
241
242 if (_mu_cut != numeric_limits<double>::infinity()) {
243 ostr << ", mass-drop cut mu=max(m1,m2)/m < " << _mu_cut;
244 } else {
245 ostr << ", no mass-drop requirement";
246 }
247
248 ostr << ", recursion into the subjet with larger ";
249 switch(_recursion_choice) {
250 case larger_pt:
251 ostr << "pt"; break;
252 case larger_mt:
253 ostr << "mt(=sqrt(m^2+pt^2))"; break;
254 case larger_m:
255 ostr << "mass"; break;
256 default:
257 cerr << "failed to interpret recursion_choice" << endl; exit(-1);
258 }
259
260 if (_subtractor) {
261 ostr << " and subtractor: " << _subtractor->description();
262 if (_input_jet_is_subtracted) {ostr << " (input jet is assumed already subtracted)";}
263 }
264 return ostr.str();
265}
266
267// decide what to return when no substructure has been found
268PseudoJet RecursiveSymmetryCutBase::_result_no_substructure(const PseudoJet &last_parent) const{
269 if (_grooming_mode){
270 // in grooming mode, return the last parent
271 return last_parent;
272 } else {
273 // in tagging mode, return an empty PseudoJet
274 return PseudoJet();
275 }
276}
277
278
279} // namespace contrib
280
281FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.