Index: contrib/contribs/RecursiveTools/trunk/ChangeLog =================================================================== --- contrib/contribs/RecursiveTools/trunk/ChangeLog (revision 1280) +++ contrib/contribs/RecursiveTools/trunk/ChangeLog (revision 1281) @@ -1,692 +1,699 @@ +2021-08-16 Mon + + * RecursiveSymmetryCutBase.cc: + added potential fix for issue signalled by Pierre-Antoine Delsart + (rare divide by zero in calculation of mu2; now mu2 gets set to -1 + when the parent mass is <= 0) + 2018-11-02 Jesse Thaler * AUTHORS: updated journal for RecursiveSoftDrop 2018-10-30 Gregory Soyez * RecursiveSoftDrop.cc: fixed a few typos in comments * RecursiveSoftDrop.hh: used the native FJ Recluster tool when available (did create conflicts in some cases) 2018-06-18 Jesse Thaler * README Fixed incorrect order of zcut and beta in the README for SoftDrop example. 2018-05-29 Jesse Thaler * VERSION * NEWS Changed to 2.0.0-beta2, noted in news 2018-04-21 Jesse Thaler * AUTHORS: updated arxiv number for RecursiveSoftDrop 2018-04-21 Gavin Salam * README: updated arxiv number for RecursiveSoftDrop & Bottom-up Soft Drop. 2018-04-04 Gregory Soyez * RecursiveSoftDrop.cc (contrib): fixed syntax of calls to structure_of<...> (thanks to Attila Krasznahorkay) 2018-01-24 Gregory Soyez * RecursiveSoftDrop.cc: for the (dafault) dynamical R0 implementation, the dynamical R0 is evolved independently in each branch. 2017-10-10 Jesse Thaler * AUTHORS Added mention of BUSD * README Some tweaks to the wording 2017-10-11 Gregory Soyez * IteratedSoftDrop.hh: IteratedSoftDropInfo::size() and multiplicity() now return an unsigned int instead of a double 2017-10-10 Jesse Thaler * AUTHORS Updated journal reference for ISD * example_isd.{cc,ref}: Included soft drop multiplicity in example * README Added warning at top that documentation has not been updated * VERSION Changed to 2.0.0-beta1 2017-10-10 Gregory Soyez * example_isd.ref: updated to reflect the bugfix below * IteratedSoftDrop.cc: fixed issue in description (was taking sqrt of -ve number when there were no angular cut) * NEWS: drafted for RecursiveTools-2.0.0-beta1 * TODO: updated list in view of a beta release * example_isd.cc: pointed to the right header for IteratedSoftDrop * RecursiveSoftDrop.{hh,cc}: * IteratedSoftDrop.{hh,cc}: moved IteratedSoftDrop to its own file (replacing the old implementation) * example_advanced_usage.ref: updated reference file following the fix below. 2017-09-28 Gregory Soyez * RecursiveSymmetryCutBase.cc: when no substructure is found, keep the _symmetry, _delta_R and _mu structure variables at -1. This for example allows one to trigger on substructure by checking if delta_R>=0. Note: we could set it to 0 (as it was before) and trigger using _delta_R>0 but there might be some genuine substructure exactly collinear. 2017-09-19 Gregory Soyez * example_isd.ref: updated to the latest version of the example 2017-09-18 Gregory Soyez * Makefile: updating make check to use all the examples (failing on ISD as expected, see below) * example_bottomup_softdrop.cc: fixed typo * example_bottomup_softdrop.ref: * example_recursive_softdrop.ref: added reference output for this example * RecursiveSoftDrop.{hh,cc}: * RecursiveSymmetryCutBase.{cc,hh}: moved the "all_prongs" method from the base structure t oa standalone function in RecursiveSoftDrop.hh In practice, this is irrelevant for mMDT and SD (since pieces() gets the job done, and the substructure class does not have (as is) reliable info to get the full structure) * RecursiveSymmetryCutBase.cc: revamped a series of requests for substructure info to better handle possible recursion into deeper jet substructure. * RecursiveSoftDrop.{hh,cc}: updated "description" to reuse the info from the base class * example_isd.cc: updated to use the newer implementation of ISD. Checked that it gives the same results as the former implementation. Note: make check still needs fixing because the example now computes a different set of angularities * RecursiveSoftDrop.hh: added a few helpers to IteratedSoftDropInfo ([] operator and size, meaning it can be used as a vector >) * RecursiveSymmetryCutBase.cc: . fixed bugs in the calculation of the geometric distances for ee coordinates . fixed bug in the computation of the (zg,thetag) pairs [it was returning the groomed ones instead of the ones passing the condition] * example_recursive_softdrop.cc: set the R0 parameter to the original jet radius 2017-09-13 Gregory Soyez * example_recursive_softdrop.cc: tied up a few comments and the code output * RecursiveSymmetryCutBase.{hh,cc}: removed the unneeded _is_composite * RecursiveSoftDrop.cc: fixed issue with "verbose" dropped info on branches with no further substructure 2017-09-11 Gregory Soyez * RecursiveSoftDrop.{hh,cc}: have IteratedSoftDDrop returning by default an object of type IteratedSoftDropInfo; added several helpers 2017-09-08 Gregory Soyez * RecursiveSoftDrop.{hh,cc}: updated IteratedSoftDrop to give it the flexibility of RecursiveSoftDrop * RecursiveSymmetryCutBase.hh: fixed typo in comment * example_mmdt_ee.cc: *** ADDED *** added an example to illustrat usage in ee collisions * example_isd.cc: * BottomUpSoftDrop.cc: * IteratedSoftDrop.cc: * RecursiveSoftDrop.cc: Fixed compilation issues with FJ 3.0 (mostly the usage of features introduced only in FJ3.1) * RecursiveSymmetryCutBase.{hh,cc}: used the internal Recluster class for FJ<3.1.0 and the FJ antive one for FJ>=3.1.0 * BottomUpSoftDrop.{hh,cc}: moved the implementation of global_grooming to the source file and fixed a few trivial compilation errors 2017-09-08 Frédéric Dreyer * BottomUpSoftDrop.hh: added the global_grooming method to process full event 2017-09-07 Gregory Soyez * RecursiveSoftDrop.cc: cleaned (mostly by removing older commented-out code) * RecursiveSoftDrop.{hh,cc}: * RecursiveSymmetryCutBase.{hh,cc}: * SoftDrop.cc: added support for ee coordinates. For that, the symmetry measure has to be set to either theta_E (which uses the 3-vector angle theta) or to cos_theta_E which uses sqrt(2*[1-cos(theta)]) Accordingly, the recursion_choice can be set to larger_E to recurse into the branch with the largest energy. The larger_m mode, recorsing into the larger-mass branch is also possible but not advised (for the same reason as the pp version). * RecursiveSymmetryCutBase.{hh,cc}: switched to the Recluster class provided with FastJet. ASlso included the recluster description to RecursiveSymmetryCutBase when it is user-defined. 2017-09-06 Gregory Soyez * BottomUpSoftDrop.{hh,cc}: . renamed SoftDropStructure -> BottomUpSoftDropStructure SoftDropRecombiner -> BottomUpSoftDropRecombiner SoftDropPlugin -> BottomUpSoftDropPlugin . moved 'description' to source file (instead of header) . kept the "area" information when available (jets will just appear as having a 0 area) . added most class description (main class still missing) * RecursiveSoftDrop.cc: . put internal classes to handle CS history in an internal namespace . replaced the "switch" in the mail loop by a series of if (allows us a few simplificatins/cleaning) . more uniform treatment of issues in the presence of an angular cut (as part of the above reorganisation) * example_advanced_usage.ref: updated reference output file following the bugfix (missing "grooming mode" initialisation in one of the SoftDrop ctors) on 2017-08-01 * RecursiveSymmetryCutBase.cc: removed redundent code 2017-08-10 Gregory Soyez * RecursiveSoftDrop.cc: fixed trivial typo in variable name >>>>>>> .r1071 2017-08-04 Gregory Soyez * RecursiveSoftDrop.cc: do not impose an angular cut in IterativeSD if it is -ve. 2017-08-01 Gregory Soyez * example_recursive_softdrop.cc: added a series of optional flags * RecursiveSoftDrop.cc: fixed a few issues with the fixed depth version * RecursiveSymmetryCutBase.hh: a jet is now considered as havig substructure if deltaR>0 (coherent with released version) * SoftDrop.hh: bugfix: set the "grooming mode" by default in all ctors EDIT: caused issue with make check, fixed on 2017-09-069 (see above) * RecursiveSoftDrop.{hh,cc}: added support for the "same depth" variant * RecursiveSymmetryCutBase.cc: also associate a RecursiveSymmetryCutBase::StructureType structure to the result jet in case it is just a single particle (in grooming mode) 2017-07-31 Gregory Soyez * RecursiveSymmetryCutBase.{hh,cc}: added the option to pass an extra parameter to the symmetry cut function * RecursiveSymmetryCutBase.{hh,cc}: * ModifiedMassDropTagger.hh * SoftDrop.hh: minor adaptions due to the above change + added a few methods to query the class information (symmetry cut, beta, ...) * RecursiveSoftDrop.{hh,cc}: Added support for - a dynamical R0 - recursing only in the hardest branch - imposing a min deltaR cut Added a tentative IterativeSoftDrop class 2017-07-28 Gregory Soyez * RecursiveSoftDrop.cc: adapted to the latest changes in RecursiveSymmetryCutBase * RecursiveSymmetryCutBase.{hh,cc}: reorganised the output of the recursion step (recurse_one_step) using an enum to make iot more readable (and to fix issues where the dropped prong is actually 0, e.g. after subtraction) 2017-07-26 Gregory Soyez * example_recursive_softdrop.cc: *** ADDED *** added a draft example for the use of RecursiveSoftDrop * RecursiveSoftDrop.{hh,cc}: *** ADDED *** added a first pass at an implementation of RecursiveSoftDrop. This is based on Frederic's initial version but keeps the branching structure of the jet. Some of the features, like dynamical R0, direct access to the substructure or the same depth variant, are still unsupported. * SoftDrop.hh: declared _beta, _symmetry_cut and _R0sqr as protected (was private) so they ca n be used by RecursiveSoftDrop * RecursiveSymmetryCutBase.{hh,cc}: extracted from result() the part that performs one step of the resursion (implemented as recurse_one_step()). This is repeatedly called by result(). It has specific conventions to indicate whether or not some substructure has been found or if one ran into some issue. 2017-04-25 Kevin Zhou * IteratedSoftDrop.hh . Added Doxygen documentation * RecursiveSymmetryCutBase.hh . Added references to ISD 2017-04-25 Jesse Thaler * AUTHORS, COPYING: . Added ISD arXiv number * example_isd.{cc,ref} . Added ISD arXiv number . Changing z_cut to be optimal value (with Lambda = 1 GeV) . Tweaked formatting * IteratedSoftDrop.{hh,cc} . Added ISD arXiv number . Assert added if recluster does not return one jet. * README . Added ISD arXiv number and tweaked wording 2017-04-20 Kevin Zhou * example_isd.{cc,ref} ** ADDED ** * IteratedSoftDrop.{cc,hh} ** ADDED ** * Makefile . Added IteratedSoftDrop (ISD) as appropriate * AUTHORS . Added myself to author list . Added placeholder for ISD paper * COPYING . Added placeholder for ISD paper * README . Added description of ISD * TODO . Added tasks to integrate ISD with other classes, e.g. SD, Recluster, and the future RecursiveSoftDrop (RSD) * NEWS . Added dummy for release of 1.1.0 * VERSION: . Switched version number to 1.1.0-dev * example_advanced_usage.cc: . Added $Id$ tag, didn't change anything else 2014-07-30 Gregory Soyez * Recluster.hh: fixed the name of the #define for the header 2014-07-09 Gregory Soyez + Jesse * NEWS: release of RecursiveTools v1.0.0 * VERSION: switched version number to 1.0.0 2014-07-08 Gavin Salam * README (RecursionChoice): added ref to .hh for constness specs of virtual functions (to reduce risk of failed overloading due to wrong constness). 2014-07-08 Gregory Soyez + Jesse * README: partially rewrote the description of set_subtractor 2014-07-07 Gregory Soyez + Jesse * example_advanced_usage.cc: * example_softdrop.cc: a few small fixed in the header of the files * VERSION: switched over to 1.0.0-alpha2-devel * README: Reordered a few things and added a few details. * Makefile (check): Be quiter during "make check" * Recluster.cc (contrib): Documented the "single" ctor argument * RecursiveSymmetryCutBase.cc (contrib): If the user sets himself the reclustering, disable the "non-CA" warning (we assume that he knows what he is doing). Mentioned in the comments that non-CA reclustering has to be used at the user's risk. Also throw when th input jet has no constituents or when there is no cluster sequence after reclustering. * Recluster.cc (contrib): replaced a remaining mention to "filtering" by reclustering 2014-07-04 Jesse Thaler * VERSION . Ready for alpha release 2014-06-17 Jesse Thaler * example_advanced_usage.{cc,ref} ** ADDED ** * Makefile . New example file to test a bunch of soft drop options . Put in makefile as well . Fixed nasty memory bug with pointers to Recluster * RecursiveSymmetryCutBase.cc * example_softdrop.ref . description() now says Groomer vs. Tagger * RecursiveSymmetryCutBase.{cc,hh} . Added optional verbose logging information about kinematics of dropped branches * example_softdrop.cc * example_advanced_usage.cc . Fixed 2014-06-16 Gregory Soyez * Makefile: also install the RecursiveSymmetryuCutBase.hh header 2014-06-13 Jesse Thaler * AUTHORS . Added myself to author list . Put complete bibliographic details on papers * COPYING . Added boilerplate MC/GPLv2 statement * example.cc: ** REMOVED ** renamed to... * example_mmdt.cc ** ADDED ** * Makefile . Made name change for consistency . Made corresponding changes in Makefile * example_mmdt_sub.cc: * example_mmdt.cc: * example_recluster.cc: . light editing of comments * example_softdrop.cc: . light editing of comments . added assert for sdjet != 0, since SoftDrop is a groomer * ModifiedMassDropTagger.hh * Recluster.{cc,hh} * RecursiveSymmetryCutBase.{cc,hh} * SoftDrop.hh . Updated some comments * README . Updated to include basic usage description and some technical details * TODO: . Added some discussion points. 2014-06-13 Gregory Soyez * example_softdrop.{cc,ref}: added an example for SoftDrop * SoftDrop.{hh,cc}: * ModifiedMassDropTagger.{hh,cc}: * RecursiveSymmetryCutBase.{hh,cc}: *** ADDED *** . added a base class for both the mMDT and SoftDrop . made mMDT and SoftDrop inherit from RecursiveSymmetryCutBase . moved the reclustering to the base class. By default, both mMDT and SoftDrop now recluster the jet with C/A . added set_grooming_mode and set_tagging_mode methods to the base class * Merging the development branch 1.0-beta1-softdrop-addition back into the trunk (will correspond to revision 682) * VERSION: switched back to 1.0.0-devel * SoftDrop.{hh,cc}: added support for re-clustering through set_reclustering(bool, Recluster*). By default, reclustering is done with Cambridge/Aachen. * example_recluster.{cc,ref}: *** ADDED *** added an example of reclustering * Recluster.{hh,cc}: added a 'single' ctor argument [true by default]. When true, the hardest jet after reclustering is returned, otherwise, the result is a composite jet with all the subjets as its pieces. 2014-05-15 Gregory Soyez * VERSION: set version number to 1.0-alpha-PUWS14.1 in preparation for a fastjet-contrib release for the pileup-workshop at CERN on May 2014. 2014-04-25 Gregory Soyez * ModifiedMassDropTagger.hh: * ModifiedMassDropTagger.cc: Added comments at various places * AUTHORS: * README: Updated info about what is now included in this contrib * SoftDrop.hh: *** ADDED *** * SoftDrop.cc: *** ADDED *** * Recluster.hh: *** ADDED *** * Recluster.cc; *** ADDED *** Added tools for reclustering and softDrop 2014-04-25 Gregory Soyez branch started at revision 611 to start including SoftDrop in the Recursivetols contrib 2014-04-24 Gregory Soyez * ModifiedMassDropTagger.hh: added a mention of the fact that when result is called in the presence of a subtractor, then the output is a subtracted PseudoJet. * ModifiedMassDropTagger.hh: declared _symmetry_cut as protected (rather than provate) so it can be accessed if symmetry_cut_description() is overloaded. * example.cc: trivial typo fixed in a comment 2014-02-04 Gavin Salam * VERSION: upped it to 1.0-beta1 * example_mmdt_sub.cc (main): added an #if to make sure FJ3.1 features are only used if FJ3.1 is available. (Currently FJ3.1 is only available to FJ developers). 2014-01-26 Gavin Salam * VERSION: renamed to 1.0-beta0 * ModifiedMassDropTagger.hh: * README: added info on author names * example_mmdt_sub.ref: *** ADDED *** added reference output for the pileup test. 2014-01-25 Gavin Salam * example_mmdt_sub.cc: * Makefile: added an extra example illustrating functionality with pileup subtraction. 2014-01-24 Gavin Salam * ModifiedMassDropTagger.cc: Reorganised code so that (sub)jet.m2()>0 check is only used when absolutely necessary: so if using a scalar_z symmetry measure, whenever scalar_z < zcut, then there is no point checking the mu condition. This means that there's no issue if the (sub)jet mass is negative, and one simply recurses down into the jet. (Whereas before it would bail out, reducing the tagging efficiency). Also removed the verbose code. 2014-01-23 Gavin Salam * ModifiedMassDropTagger.cc|hh: * example.cc replaced "asymmetry" with "symmetry" in a number of places; implemented the structural information and added it to the example; added a new simplified constructor; improved doxygen documentation; started renaming -> RecursiveTools * README some tidying * VERSION 1.0-b0 2014-01-22 Gavin Salam * ModifiedMassDropTagger.cc (contrib): -ve mass now bails out also when using the "y" asymmetry measure. Also, default my is now infinite. 2014-01-20 Gavin Salam + Gregory * ModifiedMassDropTagger.cc|hh: introduced a virtual asymmetry_cut_fn (essentially a dummy function returning a constant), to allow for derived classes to do fancier things. added warning about non-C/A clustering. explicitly labelled some (inherited) virtual functions as virtual. 2014-01-20 Gavin Salam * example.ref: * example.cc (main): got a first working example and make check setup. * ModifiedMassDropTagger.cc|hh: improved doxygen comments; added option whereby input jet is assumed already subtracted 2014-01-19 Gavin Salam * ModifiedMassDropTagger.cc|hh: * many other files Initial creation, with basic code for MMDT Index: contrib/contribs/RecursiveTools/trunk/RecursiveSymmetryCutBase.cc =================================================================== --- contrib/contribs/RecursiveTools/trunk/RecursiveSymmetryCutBase.cc (revision 1280) +++ contrib/contribs/RecursiveTools/trunk/RecursiveSymmetryCutBase.cc (revision 1281) @@ -1,642 +1,647 @@ // $Id$ // // Copyright (c) 2014-, Gavin P. Salam, Gregory Soyez, Jesse Thaler // //---------------------------------------------------------------------- // This file is part of FastJet contrib. // // It is free software; you can redistribute it and/or modify it under // the terms of the GNU General Public License as published by the // Free Software Foundation; either version 2 of the License, or (at // your option) any later version. // // It is distributed in the hope that it will be useful, but WITHOUT // ANY WARRANTY; without even the implied warranty of MERCHANTABILITY // or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public // License for more details. // // You should have received a copy of the GNU General Public License // along with this code. If not, see . //---------------------------------------------------------------------- #include "RecursiveSymmetryCutBase.hh" #include "fastjet/JetDefinition.hh" #include "fastjet/ClusterSequenceAreaBase.hh" #include #include #include using namespace std; FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh namespace contrib{ LimitedWarning RecursiveSymmetryCutBase::_negative_mass_warning; LimitedWarning RecursiveSymmetryCutBase::_mu2_gt1_warning; //LimitedWarning RecursiveSymmetryCutBase::_nonca_warning; LimitedWarning RecursiveSymmetryCutBase::_explicit_ghost_warning; bool RecursiveSymmetryCutBase::_verbose = false; //---------------------------------------------------------------------- PseudoJet RecursiveSymmetryCutBase::result(const PseudoJet & jet) const { // construct the input jet (by default, recluster with C/A) if (! jet.has_constituents()){ throw Error("RecursiveSymmetryCutBase can only be applied to jets with constituents"); } PseudoJet j = _recluster_if_needed(jet); // sanity check: the jet must have a valid CS if (! j.has_valid_cluster_sequence()){ throw Error("RecursiveSymmetryCutBase can only be applied to jets associated to a (valid) cluster sequence"); } // check that area information is there in case we have a subtractor // GS: do we really need this since subtraction may not require areas? if (_subtractor) { const ClusterSequenceAreaBase * csab = dynamic_cast(j.associated_cs()); if (csab == 0 || (!csab->has_explicit_ghosts())) _explicit_ghost_warning.warn("RecursiveSymmetryCutBase: there is no clustering sequence, or it lacks explicit ghosts: subtraction is not guaranteed to function properly"); } // establish the first subjet and optionally subtract it PseudoJet subjet = j; if (_subtractor && (!_input_jet_is_subtracted)) { subjet = (*_subtractor)(subjet); } // variables for tracking what will happen PseudoJet piece1, piece2; // vectors for storing optional verbose structure // these hold the deltaR, symmetry, and mu values of dropped branches std::vector dropped_delta_R; std::vector dropped_symmetry; std::vector dropped_mu; double sym, mu2; // now recurse into the jet's structure RecursionStatus status; while ((status=recurse_one_step(subjet, piece1, piece2, sym, mu2)) != recursion_success) { // start with sanity checks: if ((status == recursion_issue) || (status == recursion_no_parents)) { // we should return piece1 by our convention for recurse_one_step PseudoJet result; if (status == recursion_issue){ result = piece1; if (_verbose) cout << "reached end; returning null jet " << endl; } else { result = _result_no_substructure(piece1); if (_verbose) cout << "no parents found; returning last PJ or empty jet" << endl; } if (result != 0) { // if in grooming mode, add dummy structure information StructureType * structure = new StructureType(result); // structure->_symmetry = 0.0; // structure->_mu = 0.0; // structure->_delta_R = 0.0; if (_verbose_structure) { // still want to store verbose information about dropped branches structure->_has_verbose = true; structure->_dropped_symmetry = dropped_symmetry; structure->_dropped_mu = dropped_mu; structure->_dropped_delta_R = dropped_delta_R; } result.set_structure_shared_ptr(SharedPtr(structure)); } return result; } assert(status == recursion_dropped); // if desired, store information about dropped branches before recursing if (_verbose_structure) { dropped_delta_R.push_back(piece1.delta_R(piece2)); dropped_symmetry.push_back(sym); dropped_mu.push_back((mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2)); } subjet = piece1; } // we've tagged the splitting, return the jet with its substructure StructureType * structure = new StructureType(subjet); structure->_symmetry = sym; structure->_mu = (mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2); structure->_delta_R = sqrt(squared_geometric_distance(piece1, piece2)); if (_verbose_structure) { structure->_has_verbose = true; structure->_dropped_symmetry = dropped_symmetry; structure->_dropped_mu = dropped_mu; structure->_dropped_delta_R = dropped_delta_R; } subjet.set_structure_shared_ptr(SharedPtr(structure)); return subjet; } //---------------------------------------------------------------------- // the method below is the one actually performing one step of the // recursion. // // It returns a status code (defined above) // // In case of success, all the information is filled // In case of "no parents", piee1 is the same subjet // In case of trouble, piece2 will be a 0 PJ and piece1 is the PJ we // should return (either 0 itself if the issue was critical, or // non-wero in case of a minor issue just causing the recursion to // stop) RecursiveSymmetryCutBase::RecursionStatus RecursiveSymmetryCutBase::recurse_one_step(const PseudoJet & subjet, PseudoJet &piece1, PseudoJet &piece2, double &sym, double &mu2, void *extra_parameters) const { if (!subjet.has_parents(piece1, piece2)){ piece1 = subjet; piece2 = PseudoJet(); return recursion_no_parents; } // first sanity check: // - zero or negative pts are not allowed for the input subjet // - zero or negative masses are not allowed for configurations // in which the mass will effectively appear in a denominator // (The masses will be checked later) if (subjet.pt2() <= 0){ // this is a critical problem, return an empty PJ piece1 = piece2 = PseudoJet(); return recursion_issue; } if (_subtractor) { piece1 = (*_subtractor)(piece1); piece2 = (*_subtractor)(piece2); } // determine the symmetry parameter if (_symmetry_measure == y) { // the original d_{ij}/m^2 choice from MDT // first make sure the mass is sensible if (subjet.m2() <= 0) { _negative_mass_warning.warn("RecursiveSymmetryCutBase: cannot calculate y, because (sub)jet mass is negative; bailing out"); // since rounding errors can give -ve masses, be a it more // tolerant and consider that no substructure has been found piece1 = _result_no_substructure(subjet); piece2 = PseudoJet(); return recursion_issue; } sym = piece1.kt_distance(piece2) / subjet.m2(); } else if (_symmetry_measure == vector_z) { // min(pt1, pt2)/(pt), where the denominator is a vector sum // of the two subjets sym = min(piece1.pt(), piece2.pt()) / subjet.pt(); } else if (_symmetry_measure == scalar_z) { // min(pt1, pt2)/(pt1+pt2), where the denominator is a scalar sum // of the two subjets double pt1 = piece1.pt(); double pt2 = piece2.pt(); // make sure denominator is non-zero sym = pt1 + pt2; if (sym == 0){ // this is a critical problem, return an empty PJ piece1 = piece2 = PseudoJet(); return recursion_issue; } sym = min(pt1, pt2) / sym; } else if ((_symmetry_measure == theta_E) || (_symmetry_measure == cos_theta_E)){ // min(E1, E2)/(E1+E2) double E1 = piece1.E(); double E2 = piece2.E(); // make sure denominator is non-zero sym = E1 + E2; if (sym == 0){ // this is a critical problem, return an empty PJ piece1 = piece2 = PseudoJet(); return recursion_issue; } sym = min(E1, E2) / sym; } else { throw Error ("Unrecognized choice of symmetry_measure"); } // determine the symmetry cut // (This function is specified in the derived classes) double this_symmetry_cut = symmetry_cut_fn(piece1, piece2, extra_parameters); // and make a first tagging decision based on symmetry cut bool tagged = (sym > this_symmetry_cut); // if tagged based on symmetry cut, then check the mu cut (if relevant) // and update the tagging decision. Calculate mu^2 regardless, for cases // of users not cutting on mu2, but still interested in its value. bool use_mu_cut = (_mu_cut != numeric_limits::infinity()); - mu2 = max(piece1.m2(), piece2.m2())/subjet.m2(); + if (subjet.m2() > 0) { + mu2 = max(piece1.m2(), piece2.m2())/subjet.m2(); + } else { + // use this to signal problems + mu2 = -1.0; + } if (tagged && use_mu_cut) { // first a sanity check -- mu2 won't be sensible if the subjet mass // is negative, so we can't then trust the mu cut - bail out if (subjet.m2() <= 0) { _negative_mass_warning.warn("RecursiveSymmetryCutBase: cannot trust mu, because (sub)jet mass is negative; bailing out"); piece1 = piece2 = PseudoJet(); return recursion_issue; } if (mu2 > 1) _mu2_gt1_warning.warn("RecursiveSymmetryCutBase encountered mu^2 value > 1"); if (mu2 > pow(_mu_cut,2)) tagged = false; } // we'll continue unclustering, allowing for the different // ways of choosing which parent to look into if (_recursion_choice == larger_pt) { if (piece1.pt2() < piece2.pt2()) std::swap(piece1, piece2); } else if (_recursion_choice == larger_mt) { if (piece1.mt2() < piece2.mt2()) std::swap(piece1, piece2); } else if (_recursion_choice == larger_m) { if (piece1.m2() < piece2.m2()) std::swap(piece1, piece2); } else if (_recursion_choice == larger_E) { if (piece1.E() < piece2.E()) std::swap(piece1, piece2); } else { throw Error ("Unrecognized value for recursion_choice"); } return tagged ? recursion_success : recursion_dropped; } //---------------------------------------------------------------------- string RecursiveSymmetryCutBase::description() const { ostringstream ostr; ostr << "Recursive " << (_grooming_mode ? "Groomer" : "Tagger") << " with a symmetry cut "; switch(_symmetry_measure) { case y: ostr << "y"; break; case scalar_z: ostr << "scalar_z"; break; case vector_z: ostr << "vector_z"; break; case theta_E: ostr << "theta_E"; break; case cos_theta_E: ostr << "cos_theta_E"; break; default: cerr << "failed to interpret symmetry_measure" << endl; exit(-1); } ostr << " > " << symmetry_cut_description(); if (_mu_cut != numeric_limits::infinity()) { ostr << ", mass-drop cut mu=max(m1,m2)/m < " << _mu_cut; } else { ostr << ", no mass-drop requirement"; } ostr << ", recursion into the subjet with larger "; switch(_recursion_choice) { case larger_pt: ostr << "pt"; break; case larger_mt: ostr << "mt(=sqrt(m^2+pt^2))"; break; case larger_m: ostr << "mass"; break; case larger_E: ostr << "energy"; break; default: cerr << "failed to interpret recursion_choice" << endl; exit(-1); } if (_subtractor) { ostr << ", subtractor: " << _subtractor->description(); if (_input_jet_is_subtracted) {ostr << " (input jet is assumed already subtracted)";} } if (_recluster) { ostr << " and reclustering using " << _recluster->description(); } return ostr.str(); } //---------------------------------------------------------------------- // helper for handling the reclustering PseudoJet RecursiveSymmetryCutBase::_recluster_if_needed(const PseudoJet &jet) const{ if (! _do_reclustering) return jet; if (_recluster) return (*_recluster)(jet); if (is_ee()){ #if FASTJET_VERSION_NUMBER >= 30100 return Recluster(JetDefinition(ee_genkt_algorithm, JetDefinition::max_allowable_R, 0.0), true)(jet); #else return Recluster(JetDefinition(ee_genkt_algorithm, JetDefinition::max_allowable_R, 0.0))(jet); #endif } return Recluster(cambridge_algorithm, JetDefinition::max_allowable_R)(jet); } //---------------------------------------------------------------------- // decide what to return when no substructure has been found double RecursiveSymmetryCutBase::squared_geometric_distance(const PseudoJet &j1, const PseudoJet &j2) const{ if (_symmetry_measure == theta_E){ double dot_3d = j1.px()*j2.px() + j1.py()*j2.py() + j1.pz()*j2.pz(); double cos_theta = max(-1.0,min(1.0, dot_3d/sqrt(j1.modp2()*j2.modp2()))); double theta = acos(cos_theta); return theta*theta; } else if (_symmetry_measure == cos_theta_E){ double dot_3d = j1.px()*j2.px() + j1.py()*j2.py() + j1.pz()*j2.pz(); return max(0.0, 2*(1-dot_3d/sqrt(j1.modp2()*j2.modp2()))); } return j1.squared_distance(j2); } //---------------------------------------------------------------------- PseudoJet RecursiveSymmetryCutBase::_result_no_substructure(const PseudoJet &last_parent) const{ if (_grooming_mode){ // in grooming mode, return the last parent return last_parent; } else { // in tagging mode, return an empty PseudoJet return PseudoJet(); } } //======================================================================== // implementation of the details of the structure // the number of dropped subjets int RecursiveSymmetryCutBase::StructureType::dropped_count(bool global) const { check_verbose("dropped_count()"); // if this jet has no substructure, just return an empty vector if (!has_substructure()) return _dropped_delta_R.size(); // deal with the non-global case if (!global) return _dropped_delta_R.size(); // for the global case, we've unfolded the recursion (likely more // efficient as it requires less copying) unsigned int count = 0; vector to_parse; to_parse.push_back(this); unsigned int i_parse = 0; while (i_parse_dropped_delta_R.size(); // check if we need to recurse deeper in the substructure // // we can have 2 situations here for the underlying structure (the // one we've wrapped around): // - it's of the clustering type // - it's a composite jet // only in the 2nd case do we have to recurse deeper const CompositeJetStructure *css = dynamic_cast(current->_structure.get()); if (css == 0){ ++i_parse; continue; } vector prongs = css->pieces(PseudoJet()); // argument irrelevant assert(prongs.size() == 2); for (unsigned int i_prong=0; i_prong<2; ++i_prong){ if (prongs[i_prong].has_structure_of()){ RecursiveSymmetryCutBase::StructureType* prong_structure = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr(); if (prong_structure->has_substructure()) to_parse.push_back(prong_structure); } } ++i_parse; } return count; } // the delta_R of all the dropped subjets vector RecursiveSymmetryCutBase::StructureType::dropped_delta_R(bool global) const { check_verbose("dropped_delta_R()"); // if this jet has no substructure, just return an empty vector if (!has_substructure()) return vector(); // deal with the non-global case if (!global) return _dropped_delta_R; // for the global case, we've unfolded the recursion (likely more // efficient as it requires less copying) vector all_dropped; vector to_parse; to_parse.push_back(this); unsigned int i_parse = 0; while (i_parse_dropped_delta_R.begin(), current->_dropped_delta_R.end()); // check if we need to recurse deeper in the substructure // // we can have 2 situations here for the underlying structure (the // one we've wrapped around): // - it's of the clustering type // - it's a composite jet // only in the 2nd case do we have to recurse deeper const CompositeJetStructure *css = dynamic_cast(current->_structure.get()); if (css == 0){ ++i_parse; continue; } vector prongs = css->pieces(PseudoJet()); // argument irrelevant assert(prongs.size() == 2); for (unsigned int i_prong=0; i_prong<2; ++i_prong){ if (prongs[i_prong].has_structure_of()){ RecursiveSymmetryCutBase::StructureType* prong_structure = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr(); if (prong_structure->has_substructure()) to_parse.push_back(prong_structure); } } ++i_parse; } return all_dropped; } // the symmetry of all the dropped subjets vector RecursiveSymmetryCutBase::StructureType::dropped_symmetry(bool global) const { check_verbose("dropped_symmetry()"); // if this jet has no substructure, just return an empty vector if (!has_substructure()) return vector(); // deal with the non-global case if (!global) return _dropped_symmetry; // for the global case, we've unfolded the recursion (likely more // efficient as it requires less copying) vector all_dropped; vector to_parse; to_parse.push_back(this); unsigned int i_parse = 0; while (i_parse_dropped_symmetry.begin(), current->_dropped_symmetry.end()); // check if we need to recurse deeper in the substructure // // we can have 2 situations here for the underlying structure (the // one we've wrapped around): // - it's of the clustering type // - it's a composite jet // only in the 2nd case do we have to recurse deeper const CompositeJetStructure *css = dynamic_cast(current->_structure.get()); if (css == 0){ ++i_parse; continue; } vector prongs = css->pieces(PseudoJet()); // argument irrelevant assert(prongs.size() == 2); for (unsigned int i_prong=0; i_prong<2; ++i_prong){ if (prongs[i_prong].has_structure_of()){ RecursiveSymmetryCutBase::StructureType* prong_structure = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr(); if (prong_structure->has_substructure()) to_parse.push_back(prong_structure); } } ++i_parse; } return all_dropped; } // the mu of all the dropped subjets vector RecursiveSymmetryCutBase::StructureType::dropped_mu(bool global) const { check_verbose("dropped_mu()"); // if this jet has no substructure, just return an empty vector if (!has_substructure()) return vector(); // deal with the non-global case if (!global) return _dropped_mu; // for the global case, we've unfolded the recursion (likely more // efficient as it requires less copying) vector all_dropped; vector to_parse; to_parse.push_back(this); unsigned int i_parse = 0; while (i_parse_dropped_mu.begin(), current->_dropped_mu.end()); // check if we need to recurse deeper in the substructure // // we can have 2 situations here for the underlying structure (the // one we've wrapped around): // - it's of the clustering type // - it's a composite jet // only in the 2nd case do we have to recurse deeper const CompositeJetStructure *css = dynamic_cast(current->_structure.get()); if (css == 0){ ++i_parse; continue; } vector prongs = css->pieces(PseudoJet()); // argument irrelevant assert(prongs.size() == 2); for (unsigned int i_prong=0; i_prong<2; ++i_prong){ if (prongs[i_prong].has_structure_of()){ RecursiveSymmetryCutBase::StructureType* prong_structure = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr(); if (prong_structure->has_substructure()) to_parse.push_back(prong_structure); } } ++i_parse; } return all_dropped; } // the maximum of the symmetry over the dropped subjets double RecursiveSymmetryCutBase::StructureType::max_dropped_symmetry(bool global) const { check_verbose("max_dropped_symmetry()"); // if there is no substructure, just exit if (!has_substructure()){ return 0.0; } // local value of the max_dropped_symmetry double local_max = (_dropped_symmetry.size() == 0) ? 0.0 : *max_element(_dropped_symmetry.begin(),_dropped_symmetry.end()); // recurse down the structure if instructed to do so if (global){ // we can have 2 situations here for the underlying structure (the // one we've wrapped around): // - it's of the clustering type // - it's a composite jet // only in the 2nd case do we have to recurse deeper const CompositeJetStructure *css = dynamic_cast(_structure.get()); if (css == 0) return local_max; vector prongs = css->pieces(PseudoJet()); // argument irrelevant assert(prongs.size() == 2); for (unsigned int i_prong=0; i_prong<2; ++i_prong){ // check if the prong has further substructure if (prongs[i_prong].has_structure_of()){ RecursiveSymmetryCutBase::StructureType* prong_structure = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr(); local_max = max(local_max, prong_structure->max_dropped_symmetry(true)); } } } return local_max; } //------------------------------------------------------------------------ // helper class to sort by decreasing thetag class SortRecursiveSoftDropStructureZgThetagPair{ public: bool operator()(const pair &p1, const pair &p2) const{ return p1.second > p2.second; } }; //------------------------------------------------------------------------ // the (zg,thetag) pairs of all the splitting that were found and passed the SD condition vector > RecursiveSymmetryCutBase::StructureType::sorted_zg_and_thetag() const { //check_verbose("sorted_zg_and_thetag()"); // if this jet has no substructure, just return an empty vector if (!has_substructure()) return vector >(); // otherwise fill a vector with all the prongs (no specific ordering) vector > all; vector to_parse; to_parse.push_back(this); unsigned int i_parse = 0; while (i_parse(current->_symmetry, current->_delta_R)); vector prongs = current->pieces(PseudoJet()); assert(prongs.size() == 2); for (unsigned int i_prong=0; i_prong<2; ++i_prong){ if (prongs[i_prong].has_structure_of()){ RecursiveSymmetryCutBase::StructureType* prong_structure = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr(); if (prong_structure->has_substructure()) to_parse.push_back(prong_structure); } } ++i_parse; } sort(all.begin(), all.end(), SortRecursiveSoftDropStructureZgThetagPair()); return all; } } // namespace contrib FASTJET_END_NAMESPACE