Fork me on GitHub

source: git/external/fastjet/contribs/RecursiveTools/RecursiveSoftDrop.cc@ 4c9835d

Last change on this file since 4c9835d was cb80e6f, checked in by Pavel Demin <pavel.demin@…>, 4 years ago

update FastJet library to 3.3.4 and FastJet Contrib library to 1.045

  • Property mode set to 100644
File size: 19.0 KB
Line 
1// $Id: RecursiveSoftDrop.cc 1192 2018-10-30 16:08:36Z gsoyez $
2//
3// Copyright (c) 2017-, Gavin P. Salam, Gregory Soyez, Jesse Thaler,
4// Kevin Zhou, Frederic Dreyer
5//
6//----------------------------------------------------------------------
7// This file is part of FastJet contrib.
8//
9// It is free software; you can redistribute it and/or modify it under
10// the terms of the GNU General Public License as published by the
11// Free Software Foundation; either version 2 of the License, or (at
12// your option) any later version.
13//
14// It is distributed in the hope that it will be useful, but WITHOUT
15// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
17// License for more details.
18//
19// You should have received a copy of the GNU General Public License
20// along with this code. If not, see <http://www.gnu.org/licenses/>.
21//----------------------------------------------------------------------
22
23#include "RecursiveSoftDrop.hh"
24#include "fastjet/ClusterSequence.hh"
25
26using namespace std;
27
28FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
29
30namespace contrib{
31
32namespace internal_recursive_softdrop{
33
34//========================================================================
35/// \class RSDHistoryElement
36/// a helper class to help keeping track od the RSD tree
37///
38/// The element is created at the top of a branch and updated each
39/// time one grooms something away.
40class RSDHistoryElement{
41public:
42 RSDHistoryElement(const PseudoJet &jet, const RecursiveSoftDrop* rsd_ptr, double R0sqr) :
43 R0_squared(R0sqr),
44 child1_in_history(-1), child2_in_history(-1), symmetry(-1.0), mu2(-1.0){
45 reset(jet, rsd_ptr);
46 }
47
48 void reset(const PseudoJet &jet, const RecursiveSoftDrop* rsd_ptr){
49 current_in_ca_tree = jet.cluster_hist_index();
50 PseudoJet piece1, piece2;
51 theta_squared = (jet.has_parents(piece1, piece2))
52 ? rsd_ptr->squared_geometric_distance(piece1,piece2) : 0.0;
53 }
54
55 int current_in_ca_tree; ///< (history) index of the current particle in the C/A tree
56 double theta_squared; ///< squared angle at which this decays
57 double R0_squared; ///< squared angle at the top of the branch
58 ///< (used for RSD with dynamic_R0)
59 int child1_in_history; ///< hardest of the 2 decay products (-1 if untagged)
60 int child2_in_history; ///< softest of the 2 decay products (-1 if untagged)
61
62 // info about what has been dropped and the local substructure
63 vector<double> dropped_delta_R;
64 vector<double> dropped_symmetry;
65 vector<double> dropped_mu;
66 double symmetry, mu2;
67};
68
69
70/// \class OrderRSDHistoryElements
71/// angular ordering of (pointers to) the history elements
72///
73/// our priority queue will use pointers to these elements that are
74/// ordered in angle (of the objects they point to)
75class OrderRSDHistoryElements{
76public:
77 bool operator()(const RSDHistoryElement *e1, const RSDHistoryElement *e2) const {
78 return e1->theta_squared < e2->theta_squared;
79 }
80};
81
82} // internal_recursive_softdrop
83
84//========================================================================
85
86// initialise all the flags and parameters to their default value
87void RecursiveSoftDrop::set_defaults(){
88 set_fixed_depth_mode(false);
89 set_dynamical_R0(false);
90 set_hardest_branch_only(false);
91 set_min_deltaR_squared(-1.0);
92}
93
94// description of the tool
95string RecursiveSoftDrop::description() const{
96 ostringstream res;
97 res << "recursive application of ["
98 << RecursiveSymmetryCutBase::description()
99 << "]";
100
101 if (_fixed_depth){
102 res << ", recursively applied down to a maximal depth of N=";
103 if (_n==-1) res << "infinity"; else res << _n;
104 } else {
105 res << ", applied N=";
106 if (_n==-1) res << "infinity"; else res << _n;
107 res << " times";
108 }
109
110 if (_dynamical_R0)
111 res << ", with R0 dynamically scaled";
112 else
113 res << ", with R0 kept fixed";
114
115 if (_hardest_branch_only)
116 res << ", following only the hardest branch";
117
118 if (_min_dR2>0)
119 res << ", with minimal angle (squared) = " << _min_dR2;
120
121 return res.str();
122}
123
124
125// action on a single jet with RecursiveSoftDrop.
126//
127// uses "result_fixed_tags" by default (i.e. recurse from R0 to
128// smaller angles until n SD conditions have been met), or
129// "result_fixed_depth" where each of the previous SD branches are
130// recirsed into down to a depth of n.
131PseudoJet RecursiveSoftDrop::result(const PseudoJet &jet) const{
132 return _fixed_depth ? result_fixed_depth(jet) : result_fixed_tags(jet);
133}
134
135// this routine applies the Soft Drop criterion recursively on the
136// CA tree until we find n subjets (or until it converges), and
137// adds them together into a groomed PseudoJet
138PseudoJet RecursiveSoftDrop::result_fixed_tags(const PseudoJet &jet) const {
139 // start by reclustering jet with C/A algorithm
140 PseudoJet ca_jet = _recluster_if_needed(jet);
141
142 if (! ca_jet.has_valid_cluster_sequence()){
143 throw Error("RecursiveSoftDrop can only be applied to jets associated to a (valid) cluster sequence");
144 }
145
146 const ClusterSequence *cs = ca_jet.validated_cluster_sequence();
147 const vector<ClusterSequence::history_element> &cs_history = cs->history();
148 const vector<PseudoJet> &cs_jets = cs->jets();
149
150 // initialise counter to 1 subjet (i.e. the full ca_jet)
151 int n_tagged = 0;
152 int max_njet = ca_jet.constituents().size();
153
154 // create the list of branches
155 unsigned int max_history_size = 2*max_njet;
156 if ((_n>0) && (_n<max_njet-1)){ max_history_size = 2*(_n+1); }
157
158 // we need to pre-allocate the space for the vector so that the
159 // pointers are not invalidated
160 vector<internal_recursive_softdrop::RSDHistoryElement> history;
161 history.reserve(max_history_size); // could be one shorter
162 history.push_back(internal_recursive_softdrop::RSDHistoryElement(ca_jet, this, _R0sqr));
163
164 // create a priority queue containing the subjets and a comparison definition
165 priority_queue<internal_recursive_softdrop::RSDHistoryElement*, vector<internal_recursive_softdrop::RSDHistoryElement*>, internal_recursive_softdrop::OrderRSDHistoryElements> active_branches;
166 active_branches.push(& (history[0]));
167
168 PseudoJet parent, piece1, piece2;
169 double sym, mu2;
170
171 // loop over C/A tree until we reach the appropriate number of subjets
172 while ((continue_grooming(n_tagged)) && (active_branches.size())) {
173 // get the element corresponding to the max dR and the associated PJ
174 internal_recursive_softdrop::RSDHistoryElement * elm = active_branches.top();
175 PseudoJet parent = cs_jets[cs_history[elm->current_in_ca_tree].jetp_index];
176
177 // do one step of SD
178 RecursionStatus status = recurse_one_step(parent, piece1, piece2, sym, mu2, &elm->R0_squared);
179
180 // check if we passed the SD condition
181 if (status==recursion_success){
182 // check for the optional angular cut
183 if ((_min_dR2 > 0) && (squared_geometric_distance(piece1,piece2) < _min_dR2))
184 break;
185
186 // both subjets are kept in the list for potential further de-clustering
187 elm->child1_in_history = history.size();
188 elm->child2_in_history = history.size()+1;
189 elm->symmetry = sym;
190 elm->mu2 = mu2;
191 active_branches.pop();
192
193 // update the history
194 double next_R0_squared = (_dynamical_R0)
195 ? piece1.squared_distance(piece2) : elm->R0_squared;
196
197 internal_recursive_softdrop::RSDHistoryElement elm1(piece1, this, next_R0_squared);
198 history.push_back(elm1);
199 active_branches.push(&(history.back()));
200 internal_recursive_softdrop::RSDHistoryElement elm2(piece2, this, next_R0_squared);
201 history.push_back(elm2);
202 if (!_hardest_branch_only){
203 active_branches.push(&(history.back()));
204 }
205
206 ++n_tagged;
207 } else if (status==recursion_dropped){
208 // check for the optional angular cut
209 if ((_min_dR2 > 0) && (squared_geometric_distance(piece1,piece2) < _min_dR2))
210 break;
211
212 active_branches.pop();
213 // tagging failed and the softest branch should be dropped
214 // keep track of what has been groomed away
215 max_njet -= piece2.constituents().size();
216 elm->dropped_delta_R .push_back((elm->theta_squared >= 0) ? sqrt(elm->theta_squared) : -sqrt(elm->theta_squared));
217 elm->dropped_symmetry.push_back(sym);
218 elm->dropped_mu .push_back((mu2>=0) ? sqrt(mu2) : -sqrt(mu2));
219
220 // keep the hardest branch in the recursion
221 elm->reset(piece1, this);
222 active_branches.push(elm);
223 } else if (status==recursion_no_parents){
224 if (_min_dR2 > 0) break;
225 active_branches.pop();
226 // nothing specific to do: we just keep the curent jet as a "leaf"
227 } else { // recursion_issue
228 active_branches.pop();
229 // we've met an issue
230 // if the piece2 is null as well, it means we've had a critical problem.
231 // In that case, return an empty PseudoJet
232 if (piece2 == 0) return PseudoJet();
233
234 // otherwise, we should consider "piece2" as a final particle
235 // not to be recursed into
236 if (_min_dR2 > 0) break;
237 max_njet -= (piece2.constituents().size()-1);
238 break;
239 }
240
241 // If the missing number of tags is exactly the number of objects
242 // we have left in the recursion, stop
243 if (n_tagged == max_njet) break;
244 }
245
246 // now we have a bunch of history elements that we can use to build the final jet
247 vector<PseudoJet> mapped_to_history(history.size());
248 unsigned int history_index = history.size();
249 do {
250 --history_index;
251 const internal_recursive_softdrop::RSDHistoryElement & elm = history[history_index];
252
253 // two kinds of events: either just a final leave, potentially with grooming
254 // or a brandhing (also with potential grooming at the end)
255 if (elm.child1_in_history<0){
256 // this is a leaf, i.e. with no further substructure
257 PseudoJet & subjet = mapped_to_history[history_index]
258 = cs_jets[cs_history[elm.current_in_ca_tree].jetp_index];
259
260 StructureType * structure = new StructureType(subjet);
261 if (has_verbose_structure()){
262 structure->set_verbose(true);
263 structure->set_dropped_delta_R (elm.dropped_delta_R);
264 structure->set_dropped_symmetry(elm.dropped_symmetry);
265 structure->set_dropped_mu (elm.dropped_mu);
266 }
267 subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
268 } else {
269 PseudoJet & subjet = mapped_to_history[history_index]
270 = join(mapped_to_history[elm.child1_in_history], mapped_to_history[elm.child2_in_history]);
271 StructureType * structure = new StructureType(subjet, sqrt(elm.theta_squared), elm.symmetry, sqrt(elm.mu2));
272 if (has_verbose_structure()){
273 structure->set_verbose(true);
274 structure->set_dropped_delta_R (elm.dropped_delta_R);
275 structure->set_dropped_symmetry(elm.dropped_symmetry);
276 structure->set_dropped_mu (elm.dropped_mu);
277 }
278 subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
279 }
280 } while (history_index>0);
281
282 return mapped_to_history[0];
283}
284
285// this routine applies the Soft Drop criterion recursively on the
286// CA tree, recursing into all the branches found during the previous iteration
287// until n layers have been found (or until it converges)
288PseudoJet RecursiveSoftDrop::result_fixed_depth(const PseudoJet &jet) const {
289 // start by reclustering jet with C/A algorithm
290 PseudoJet ca_jet = _recluster_if_needed(jet);
291
292 if (! ca_jet.has_valid_cluster_sequence()){
293 throw Error("RecursiveSoftDrop can only be applied to jets associated to a (valid) cluster sequence");
294 }
295
296 const ClusterSequence *cs = ca_jet.validated_cluster_sequence();
297 const vector<ClusterSequence::history_element> &cs_history = cs->history();
298 const vector<PseudoJet> &cs_jets = cs->jets();
299
300 // initialise counter to 1 subjet (i.e. the full ca_jet)
301 int n_depth = 0;
302 int max_njet = ca_jet.constituents().size();
303
304 // create the list of branches
305 unsigned int max_history_size = 2*max_njet;
306 //if ((_n>0) && (_n<max_njet-1)){ max_history_size = 2*(_n+1); }
307
308 // we need to pre-allocate the space for the vector so that the
309 // pointers are not invalidated
310 vector<internal_recursive_softdrop::RSDHistoryElement> history;
311 history.reserve(max_history_size); // could be one shorter
312 history.push_back(internal_recursive_softdrop::RSDHistoryElement(ca_jet, this, _R0sqr));
313 history.back().theta_squared = _R0sqr;
314
315 // create a priority queue containing the subjets and a comparison definition
316 list<internal_recursive_softdrop::RSDHistoryElement*> active_branches;
317 active_branches.push_back(& (history[0]));
318
319 PseudoJet parent, piece1, piece2;
320
321 while ((continue_grooming(n_depth)) && (active_branches.size())) {
322 // loop over all the branches and look for substructure
323 list<internal_recursive_softdrop::RSDHistoryElement*>::iterator hist_it=active_branches.begin();
324 while (hist_it!=active_branches.end()){
325 // get the element corresponding to the max dR and the associated PJ
326 internal_recursive_softdrop::RSDHistoryElement * elm = (*hist_it);
327 PseudoJet parent = cs_jets[cs_history[elm->current_in_ca_tree].jetp_index];
328
329 // we need to iterate this branch until we find some substructure
330 PseudoJet result_sd;
331 if (_dynamical_R0){
332 SoftDrop sd(_beta, _symmetry_cut, symmetry_measure(), sqrt(elm->theta_squared),
333 mu_cut(), recursion_choice(), subtractor());
334 sd.set_reclustering(false);
335 sd.set_verbose_structure(has_verbose_structure());
336 result_sd = sd(parent);
337 } else {
338 result_sd = SoftDrop::result(parent);
339 }
340
341 // if we had an empty PJ, that means we ran into some problems.
342 // just return an empty PJ ourselves
343 if (result_sd == 0) return PseudoJet();
344
345 // update the history element to reflect our iteration
346 elm->current_in_ca_tree = result_sd.cluster_hist_index();
347
348 if (has_verbose_structure()){
349 elm->dropped_delta_R = result_sd.structure_of<SoftDrop>().dropped_delta_R();
350 elm->dropped_symmetry = result_sd.structure_of<SoftDrop>().dropped_symmetry();
351 elm->dropped_mu = result_sd.structure_of<SoftDrop>().dropped_mu();
352 }
353
354 // if some substructure was found:
355 if (result_sd.structure_of<SoftDrop>().has_substructure()){
356 // update the history element to reflect our iteration
357 elm->child1_in_history = history.size();
358 elm->child2_in_history = history.size()+1;
359 elm->theta_squared = result_sd.structure_of<SoftDrop>().delta_R();
360 elm->theta_squared *= elm->theta_squared;
361 elm->symmetry = result_sd.structure_of<SoftDrop>().symmetry();
362 elm->mu2 = result_sd.structure_of<SoftDrop>().mu();
363 elm->mu2 *= elm->mu2;
364
365 // the next iteration will have to handle 2 new history
366 // elements (the R0squared argument here is unused)
367 result_sd.has_parents(piece1, piece2);
368 internal_recursive_softdrop::RSDHistoryElement elm1(piece1, this, _R0sqr);
369 history.push_back(elm1);
370 // insert it in the active branches if needed
371 if (elm1.theta_squared>0)
372 active_branches.insert(hist_it,&(history.back())); // insert just before
373
374 internal_recursive_softdrop::RSDHistoryElement elm2(piece2, this, _R0sqr);
375 history.push_back(elm2);
376 if ((!_hardest_branch_only) && (elm2.theta_squared>0)){
377 active_branches.insert(hist_it,&(history.back())); // insert just before
378 }
379 }
380 // otherwise we've just reached the end of the recursion the
381 // history information has been updated above
382 //
383 // we just need to make sure that we do not recurse into that
384 // element any longer
385
386 list<internal_recursive_softdrop::RSDHistoryElement*>::iterator current = hist_it;
387 ++hist_it;
388 active_branches.erase(current);
389 } // loop over branches at current depth
390 ++n_depth;
391 } // loop over depth
392
393 // now we have a bunch of history elements that we can use to build the final jet
394 vector<PseudoJet> mapped_to_history(history.size());
395 unsigned int history_index = history.size();
396 do {
397 --history_index;
398 const internal_recursive_softdrop::RSDHistoryElement & elm = history[history_index];
399
400 // two kinds of events: either just a final leave, poteitially with grooming
401 // or a brandhing (also with potential grooming at the end)
402 if (elm.child1_in_history<0){
403 // this is a leaf, i.e. with no further sustructure
404 PseudoJet & subjet = mapped_to_history[history_index]
405 = cs_jets[cs_history[elm.current_in_ca_tree].jetp_index];
406
407 StructureType * structure = new StructureType(subjet);
408 if (has_verbose_structure()){
409 structure->set_verbose(true);
410 }
411 subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
412 } else {
413 PseudoJet & subjet = mapped_to_history[history_index]
414 = join(mapped_to_history[elm.child1_in_history], mapped_to_history[elm.child2_in_history]);
415 StructureType * structure = new StructureType(subjet, sqrt(elm.theta_squared), elm.symmetry, sqrt(elm.mu2));
416 if (has_verbose_structure()){
417 structure->set_verbose(true);
418 structure->set_dropped_delta_R (elm.dropped_delta_R);
419 structure->set_dropped_symmetry(elm.dropped_symmetry);
420 structure->set_dropped_mu (elm.dropped_mu);
421 }
422 subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
423 }
424 } while (history_index>0);
425
426 return mapped_to_history[0];
427}
428
429
430//========================================================================
431// implementation of the helpers
432//========================================================================
433
434// helper to get all the prongs in a jet that has been obtained using
435// RecursiveSoftDrop (instead of recursively parsing the 1->2
436// composite jet structure)
437vector<PseudoJet> recursive_soft_drop_prongs(const PseudoJet & rsd_jet){
438 // make sure that the jet has the appropriate RecursiveSoftDrop structure
439 if (!rsd_jet.has_structure_of<RecursiveSoftDrop>())
440 return vector<PseudoJet>();
441
442 // if this jet has no substructure, just return a 1-prong object
443 if (!rsd_jet.structure_of<RecursiveSoftDrop>().has_substructure())
444 return vector<PseudoJet>(1, rsd_jet);
445
446 // otherwise fill a vector with all the prongs (no specific ordering)
447 vector<PseudoJet> prongs;
448
449 // parse the list of PseudoJet we still need to deal with
450 vector<PseudoJet> to_parse = rsd_jet.pieces(); // valid both for a C/A recombination step or a RSD join
451 unsigned int i_parse = 0;
452 while (i_parse<to_parse.size()){
453 const PseudoJet &current = to_parse[i_parse];
454
455 if ((current.has_structure_of<RecursiveSymmetryCutBase>()) &&
456 (current.structure_of<RecursiveSymmetryCutBase>().has_substructure())){
457 // if this has some deeper substructure, add it to the list of
458 // things to further process
459 vector<PseudoJet> pieces = current.pieces();
460 assert(pieces.size() == 2);
461 to_parse[i_parse] = pieces[0];
462 to_parse.push_back(pieces[1]);
463 } else {
464 // no further substructure, just add this as a branch
465 prongs.push_back(current);
466 ++i_parse;
467 }
468 }
469
470 return prongs;
471}
472
473}
474
475FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.