| 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 |
|
|---|
| 28 | using namespace std;
|
|---|
| 29 |
|
|---|
| 30 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
| 31 |
|
|---|
| 32 | namespace contrib{
|
|---|
| 33 |
|
|---|
| 34 | LimitedWarning RecursiveSymmetryCutBase::_negative_mass_warning;
|
|---|
| 35 | LimitedWarning RecursiveSymmetryCutBase::_mu2_gt1_warning;
|
|---|
| 36 | //LimitedWarning RecursiveSymmetryCutBase::_nonca_warning;
|
|---|
| 37 | LimitedWarning RecursiveSymmetryCutBase::_explicit_ghost_warning;
|
|---|
| 38 |
|
|---|
| 39 | bool RecursiveSymmetryCutBase::_verbose = false;
|
|---|
| 40 |
|
|---|
| 41 | //----------------------------------------------------------------------
|
|---|
| 42 | PseudoJet 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 | //----------------------------------------------------------------------
|
|---|
| 226 | string 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
|
|---|
| 268 | PseudoJet 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 |
|
|---|
| 281 | FASTJET_END_NAMESPACE
|
|---|
