Fork me on GitHub

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

ImprovedOutputFile Timing
Last change on this file since d4f8a8c was b7b836a, checked in by Pavel Demin <pavel-demin@…>, 6 years ago

update FastJet library to 3.3.1 and FastJet Contrib library to 1.036

  • Property mode set to 100644
File size: 19.2 KB
Line 
1// $Id: RecursiveSoftDrop.cc 1111 2018-04-04 10:06:11Z 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 // initialize 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 // initialise to the full ca_jet
166 priority_queue<internal_recursive_softdrop::RSDHistoryElement*, vector<internal_recursive_softdrop::RSDHistoryElement*>, internal_recursive_softdrop::OrderRSDHistoryElements> active_branches;
167 active_branches.push(& (history[0]));
168
169 PseudoJet parent, piece1, piece2;
170 double sym, mu2;
171
172 // which R0 to use
173 //double R0sqr = _R0sqr;
174
175 // loop over C/A tree until we reach the appropriate number of subjets
176 while ((continue_grooming(n_tagged)) && (active_branches.size())) {
177 // get the element corresponding to the max dR
178 // and the associated PJ
179 internal_recursive_softdrop::RSDHistoryElement * elm = active_branches.top();
180 PseudoJet parent = cs_jets[cs_history[elm->current_in_ca_tree].jetp_index];
181
182 // do one step of SD
183 RecursionStatus status = recurse_one_step(parent, piece1, piece2, sym, mu2, &elm->R0_squared);
184
185 // check if we passed the SD condition
186 if (status==recursion_success){
187 // check for the optional angular cut
188 if ((_min_dR2 > 0) && (squared_geometric_distance(piece1,piece2) < _min_dR2))
189 break;
190
191 // both subjets are kept in the list for potential further de-clustering
192 elm->child1_in_history = history.size();
193 elm->child2_in_history = history.size()+1;
194 elm->symmetry = sym;
195 elm->mu2 = mu2;
196 active_branches.pop();
197
198 // update the history
199 double next_R0_squared = (_dynamical_R0)
200 ? piece1.squared_distance(piece2) : elm->R0_squared;
201
202 internal_recursive_softdrop::RSDHistoryElement elm1(piece1, this, next_R0_squared);
203 history.push_back(elm1);
204 active_branches.push(&(history.back()));
205 internal_recursive_softdrop::RSDHistoryElement elm2(piece2, this, next_R0_squared);
206 history.push_back(elm2);
207 if (!_hardest_branch_only){
208 active_branches.push(&(history.back()));
209 }
210
211 ++n_tagged;
212 } else if (status==recursion_dropped){
213 // check for the optional angular cut
214 if ((_min_dR2 > 0) && (squared_geometric_distance(piece1,piece2) < _min_dR2))
215 break;
216
217 active_branches.pop();
218 // tagging failed and the softest branch should be dropped
219 // keep track of what has een groomed away
220 max_njet -= piece2.constituents().size();
221 elm->dropped_delta_R .push_back((elm->theta_squared >= 0) ? sqrt(elm->theta_squared) : -sqrt(elm->theta_squared));
222 elm->dropped_symmetry.push_back(sym);
223 elm->dropped_mu .push_back((mu2>=0) ? sqrt(mu2) : -sqrt(mu2));
224
225 // keep the hardest bhanch in the recursion
226 elm->reset(piece1, this);
227 active_branches.push(elm);
228 } else if (status==recursion_no_parents){
229 if (_min_dR2 > 0) break;
230 active_branches.pop();
231 // nothing specific to do: we just keep the curent jet as a "leaf"
232 } else { // recursion_issue
233 active_branches.pop();
234 // we've met an issue
235 // if the piece2 is null as well, it means we've had a critical problem.
236 // In that case, return an empty PseudoJet
237 if (piece2 == 0) return PseudoJet();
238
239 // otherwise, we should consider "piece2" as a final particle
240 // not to be recursed into
241 if (_min_dR2 > 0) break;
242 max_njet -= (piece2.constituents().size()-1);
243 break;
244 }
245
246 // If the missing number of tags is exactly the number of objects
247 // we have left in the recursion, stop
248 if (n_tagged == max_njet) break;
249 }
250
251 // now we have a bunch of history elements that we can use to build the final jet
252 vector<PseudoJet> mapped_to_history(history.size());
253 unsigned int history_index = history.size();
254 do {
255 --history_index;
256 const internal_recursive_softdrop::RSDHistoryElement & elm = history[history_index];
257
258 // two kinds of events: either just a final leave, poteitially with grooming
259 // or a brandhing (also with potential grooming at the end)
260 if (elm.child1_in_history<0){
261 // this is a leaf, i.e. with no further sustructure
262 PseudoJet & subjet = mapped_to_history[history_index]
263 = cs_jets[cs_history[elm.current_in_ca_tree].jetp_index];
264
265 StructureType * structure = new StructureType(subjet);
266 if (has_verbose_structure()){
267 structure->set_verbose(true);
268 structure->set_dropped_delta_R (elm.dropped_delta_R);
269 structure->set_dropped_symmetry(elm.dropped_symmetry);
270 structure->set_dropped_mu (elm.dropped_mu);
271 }
272 subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
273 } else {
274 PseudoJet & subjet = mapped_to_history[history_index]
275 = join(mapped_to_history[elm.child1_in_history], mapped_to_history[elm.child2_in_history]);
276 StructureType * structure = new StructureType(subjet, sqrt(elm.theta_squared), elm.symmetry, sqrt(elm.mu2));
277 if (has_verbose_structure()){
278 structure->set_verbose(true);
279 structure->set_dropped_delta_R (elm.dropped_delta_R);
280 structure->set_dropped_symmetry(elm.dropped_symmetry);
281 structure->set_dropped_mu (elm.dropped_mu);
282 }
283 subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
284 }
285 } while (history_index>0);
286
287 return mapped_to_history[0];
288}
289
290// this routine applies the Soft Drop criterion recursively on the
291// CA tree, recursing into all the branches found during the previous iteration
292// until n layers have been found (or until it converges)
293PseudoJet RecursiveSoftDrop::result_fixed_depth(const PseudoJet &jet) const {
294 // start by reclustering jet with C/A algorithm
295 PseudoJet ca_jet = _recluster_if_needed(jet);
296
297 if (! ca_jet.has_valid_cluster_sequence()){
298 throw Error("RecursiveSoftDrop can only be applied to jets associated to a (valid) cluster sequence");
299 }
300
301 const ClusterSequence *cs = ca_jet.validated_cluster_sequence();
302 const vector<ClusterSequence::history_element> &cs_history = cs->history();
303 const vector<PseudoJet> &cs_jets = cs->jets();
304
305 // initialize counter to 1 subjet (i.e. the full ca_jet)
306 int n_depth = 0;
307 int max_njet = ca_jet.constituents().size();
308
309 // create the list of branches
310 unsigned int max_history_size = 2*max_njet;
311 //if ((_n>0) && (_n<max_njet-1)){ max_history_size = 2*(_n+1); }
312
313 // we need to pre-allocate the space for the vector so that the
314 // pointers are not invalidated
315 vector<internal_recursive_softdrop::RSDHistoryElement> history;
316 history.reserve(max_history_size); // could be one shorter
317 history.push_back(internal_recursive_softdrop::RSDHistoryElement(ca_jet, this, _R0sqr));
318 history.back().theta_squared = _R0sqr;
319
320 // create a priority queue containing the subjets and a comparison definition
321 // initialize counter to 1 subjet (i.e. the full ca_jet)
322 list<internal_recursive_softdrop::RSDHistoryElement*> active_branches;
323 active_branches.push_back(& (history[0]));
324
325 PseudoJet parent, piece1, piece2;
326
327 while ((continue_grooming(n_depth)) && (active_branches.size())) {
328 // loop over all the branches and look for substructure
329 list<internal_recursive_softdrop::RSDHistoryElement*>::iterator hist_it=active_branches.begin();
330 while (hist_it!=active_branches.end()){
331 // get the element corresponding to the max dR
332 // and the associated PJ
333 internal_recursive_softdrop::RSDHistoryElement * elm = (*hist_it);
334 PseudoJet parent = cs_jets[cs_history[elm->current_in_ca_tree].jetp_index];
335
336 // we need to iterate this branch until we find some substructure
337 PseudoJet result_sd;
338 if (_dynamical_R0){
339 SoftDrop sd(_beta, _symmetry_cut, symmetry_measure(), sqrt(elm->theta_squared),
340 mu_cut(), recursion_choice(), subtractor());
341 sd.set_reclustering(false);
342 sd.set_verbose_structure(has_verbose_structure());
343 result_sd = sd(parent);
344 } else {
345 result_sd = SoftDrop::result(parent);
346 }
347
348 // if we had an empty PJ, that means we ran into some problems.
349 // just return an empty PJ ourselves
350 if (result_sd == 0) return PseudoJet();
351
352 // update the history element to reflect our iteration
353 elm->current_in_ca_tree = result_sd.cluster_hist_index();
354
355 if (has_verbose_structure()){
356 elm->dropped_delta_R = result_sd.structure_of<SoftDrop>().dropped_delta_R();
357 elm->dropped_symmetry = result_sd.structure_of<SoftDrop>().dropped_symmetry();
358 elm->dropped_mu = result_sd.structure_of<SoftDrop>().dropped_mu();
359 }
360
361 // if some substructure was found:
362 if (result_sd.structure_of<SoftDrop>().has_substructure()){
363 // update the history element to reflect our iteration
364 elm->child1_in_history = history.size();
365 elm->child2_in_history = history.size()+1;
366 elm->theta_squared = result_sd.structure_of<SoftDrop>().delta_R();
367 elm->theta_squared *= elm->theta_squared;
368 elm->symmetry = result_sd.structure_of<SoftDrop>().symmetry();
369 elm->mu2 = result_sd.structure_of<SoftDrop>().mu();
370 elm->mu2 *= elm->mu2;
371
372 // the next iteration will have to handle 2 new history
373 // elements (the R0squared argument here is unused)
374 result_sd.has_parents(piece1, piece2);
375 internal_recursive_softdrop::RSDHistoryElement elm1(piece1, this, _R0sqr);
376 history.push_back(elm1);
377 // insert it in the active branches if needed
378 if (elm1.theta_squared>0)
379 active_branches.insert(hist_it,&(history.back())); // insert just before
380
381 internal_recursive_softdrop::RSDHistoryElement elm2(piece2, this, _R0sqr);
382 history.push_back(elm2);
383 if ((!_hardest_branch_only) && (elm2.theta_squared>0)){
384 active_branches.insert(hist_it,&(history.back())); // insert just before
385 }
386 }
387 // otherwise we've just reached the end of the recursion the
388 // history information has been updated above
389 //
390 // we just need to make sure that we do not recurse into that
391 // element any longer
392
393 list<internal_recursive_softdrop::RSDHistoryElement*>::iterator current = hist_it;
394 ++hist_it;
395 active_branches.erase(current);
396 } // loop over branches at current depth
397 ++n_depth;
398 } // loop over depth
399
400 // now we have a bunch of history elements that we can use to build the final jet
401 vector<PseudoJet> mapped_to_history(history.size());
402 unsigned int history_index = history.size();
403 do {
404 --history_index;
405 const internal_recursive_softdrop::RSDHistoryElement & elm = history[history_index];
406
407 // two kinds of events: either just a final leave, poteitially with grooming
408 // or a brandhing (also with potential grooming at the end)
409 if (elm.child1_in_history<0){
410 // this is a leaf, i.e. with no further sustructure
411 PseudoJet & subjet = mapped_to_history[history_index]
412 = cs_jets[cs_history[elm.current_in_ca_tree].jetp_index];
413
414 StructureType * structure = new StructureType(subjet);
415 if (has_verbose_structure()){
416 structure->set_verbose(true);
417 }
418 subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
419 } else {
420 PseudoJet & subjet = mapped_to_history[history_index]
421 = join(mapped_to_history[elm.child1_in_history], mapped_to_history[elm.child2_in_history]);
422 StructureType * structure = new StructureType(subjet, sqrt(elm.theta_squared), elm.symmetry, sqrt(elm.mu2));
423 if (has_verbose_structure()){
424 structure->set_verbose(true);
425 structure->set_dropped_delta_R (elm.dropped_delta_R);
426 structure->set_dropped_symmetry(elm.dropped_symmetry);
427 structure->set_dropped_mu (elm.dropped_mu);
428 }
429 subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
430 }
431 } while (history_index>0);
432
433 return mapped_to_history[0];
434}
435
436
437//========================================================================
438// implementation of the helpers
439//========================================================================
440
441// helper to get all the prongs in a jet that has been obtained using
442// RecursiveSoftDrop (instead of recursively parsing the 1->2
443// composite jet structure)
444vector<PseudoJet> recursive_soft_drop_prongs(const PseudoJet & rsd_jet){
445 // make sure that the jet has the appropriate RecursiveSoftDrop structure
446 if (!rsd_jet.has_structure_of<RecursiveSoftDrop>())
447 return vector<PseudoJet>();
448
449 // if this jet has no substructure, just return a 1-prong object
450 if (!rsd_jet.structure_of<RecursiveSoftDrop>().has_substructure())
451 return vector<PseudoJet>(1, rsd_jet);
452
453 // otherwise fill a vector with all the prongs (no specific ordering)
454 vector<PseudoJet> prongs;
455
456 // parse the list of PseudoJet we still need to deal with
457 vector<PseudoJet> to_parse = rsd_jet.pieces(); // valid both for a C/A recombination step or a RSD join
458 unsigned int i_parse = 0;
459 while (i_parse<to_parse.size()){
460 const PseudoJet &current = to_parse[i_parse];
461
462 if ((current.has_structure_of<RecursiveSymmetryCutBase>()) &&
463 (current.structure_of<RecursiveSymmetryCutBase>().has_substructure())){
464 // if this has some deeper substructure, add it to the list of
465 // things to further process
466 vector<PseudoJet> pieces = current.pieces();
467 assert(pieces.size() == 2);
468 to_parse[i_parse] = pieces[0];
469 to_parse.push_back(pieces[1]);
470 } else {
471 // no further substructure, just add this as a branch
472 prongs.push_back(current);
473 ++i_parse;
474 }
475 }
476
477 return prongs;
478}
479
480}
481
482FASTJET_END_NAMESPACE
Note: See TracBrowser for help on using the repository browser.