Context Navigation

-              r95e6b7a
+              rb7b836a
 // $Id: RecursiveSymmetryCutBase.cc 700 2014-07-07 12:50:05Z gsoyez $
+// $Id: RecursiveSymmetryCutBase.cc 1080 2017-09-28 07:51:37Z gsoyez $
 //
 // Copyright (c) 2014-, Gavin P. Salam, Gregory Soyez, Jesse Thaler
 …
 #include "fastjet/JetDefinition.hh"
 #include "fastjet/ClusterSequenceAreaBase.hh"
+#include <cassert>
 #include <algorithm>
 #include <cstdlib>
 …
 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 =
+    _do_reclustering
+      ? _recluster ? (*_recluster)(jet)
+                   : Recluster(cambridge_algorithm, JetDefinition::max_allowable_R)(jet)
+      : jet;
+  // issue a warning if the jet is not obtained through a C/A
+  // clustering
+  // if ((! j.has_associated_cluster_sequence()) ||
+  //     (j.validated_cs()->jet_def().jet_algorithm() != cambridge_algorithm))
+  //  _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.");
+  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 =
 …
     subjet = (*_subtractor)(subjet);
+  }
-  bool use_mu_cut = (_mu_cut != numeric_limits<double>::infinity());
   // variables for tracking what will happen
 …
   std::vector<double> dropped_symmetry;
   std::vector<double> dropped_mu;
+  double sym, mu2;
   // now recurse into the jet's structure
+  while (subjet.has_parents(piece1, piece2)) {
+    // 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) return PseudoJet();
+    if (_subtractor) {
+      piece1 = (*_subtractor)(piece1);
+      piece2 = (*_subtractor)(piece2);
+    }
+    // determine the symmetry parameter
+    double sym;
+    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");
+        return _result_no_substructure(subjet); //TBC: do we return the hardest parent? A NULL PseudoJet?
+      }
+      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) return PseudoJet();
+      sym = min(pt1, pt2) / 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);
+    // 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.
+    double mu2 = max(piece1.m2(), piece2.m2())/subjet.m2();
+    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");
+        return PseudoJet();
+      }
+      if (mu2 > 1) _mu2_gt1_warning.warn("RecursiveSymmetryCutBase encountered mu^2 value > 1");
+      if (mu2 > pow(_mu_cut,2)) tagged = false;
+    }
+    // if we've tagged the splitting, return the jet with its substructure
+    if (tagged) {
+      // record relevant information
+      StructureType * structure = new StructureType(subjet);
+      structure->_symmetry = sym;
+      structure->_mu       = (mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2);
+      structure->_delta_R  = piece1.delta_R(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<PseudoJetStructureBase>(structure));
+      return subjet;
+    }
+  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<PseudoJetStructureBase>(structure));
+      }
+      return result;
+    }
+    assert(status == recursion_dropped);
     // if desired, store information about dropped branches before recursing
 …
       dropped_mu.push_back((mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2));
+    }
+    // otherwise continue unclustering, allowing for the different
+    // ways of choosing which parent to look into
+    int choice;
+    if        (_recursion_choice == larger_mt) {
+      choice = piece1.mt2() > piece2.mt2() ? 1 : 2;
+    } else if (_recursion_choice == larger_pt) {
+      choice = piece1.pt2() > piece2.pt2() ? 1 : 2;
+    } else if (_recursion_choice == larger_m)  {
+      choice = piece1.m2()  > piece2.m2()  ? 1 : 2;
+    } else {
+      throw Error ("Unrecognized value for recursion_choice");
+    }
+    if (_verbose) cout << "choice is " << choice << endl;;
+    subjet = (choice == 1) ? piece1 : piece2;
+  } // (subjet.has_parents(...))
+  if (_verbose) cout << "reached end; returning null jet " << endl;
+  // decide on tagging versus grooming mode here
+  PseudoJet result = _result_no_substructure(subjet);
+  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<PseudoJetStructureBase>(structure));
+  }
+  return result;
+}
+    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<PseudoJetStructureBase>(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<double>::infinity());
+  mu2 = max(piece1.m2(), piece2.m2())/subjet.m2();
+  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 {
 …
   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);
 …
   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 << " and subtractor: " << _subtractor->description();
+    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){
 …
+//========================================================================
+// 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<const RecursiveSymmetryCutBase::StructureType*> to_parse;
+  to_parse.push_back(this);
+  unsigned int i_parse = 0;
+  while (i_parse<to_parse.size()){
+    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
+    count += current->_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<const CompositeJetStructure*>(current->_structure.get());
+    if (css == 0){ ++i_parse; continue; }
+    vector<PseudoJet> 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>()){
+        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<double> 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<double>();
+  // 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<double> all_dropped;
+  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
+  to_parse.push_back(this);
+  unsigned int i_parse = 0;
+  while (i_parse<to_parse.size()){
+    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
+    all_dropped.insert(all_dropped.end(), current->_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<const CompositeJetStructure*>(current->_structure.get());
+    if (css == 0){ ++i_parse; continue; }
+    vector<PseudoJet> 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>()){
+        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<double> 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<double>();
+  // 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<double> all_dropped;
+  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
+  to_parse.push_back(this);
+  unsigned int i_parse = 0;
+  while (i_parse<to_parse.size()){
+    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
+    all_dropped.insert(all_dropped.end(), current->_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<const CompositeJetStructure*>(current->_structure.get());
+    if (css == 0){ ++i_parse; continue; }
+    vector<PseudoJet> 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>()){
+        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<double> 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<double>();
+  // 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<double> all_dropped;
+  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
+  to_parse.push_back(this);
+  unsigned int i_parse = 0;
+  while (i_parse<to_parse.size()){
+    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
+    all_dropped.insert(all_dropped.end(), current->_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<const CompositeJetStructure*>(current->_structure.get());
+    if (css == 0){ ++i_parse; continue; }
+    vector<PseudoJet> 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>()){
+        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<const CompositeJetStructure*>(_structure.get());
+    if (css == 0) return local_max;
+    vector<PseudoJet> 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>()){
+        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<double, double> &p1, const pair<double, double> &p2) const{
+    return p1.second > p2.second;
+  }
+};
+//------------------------------------------------------------------------
+// the (zg,thetag) pairs of all the splitting that were found and passed the SD condition
+vector<pair<double,double> > 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<pair<double,double> >();
+  // otherwise fill a vector with all the prongs (no specific ordering)
+  vector<pair<double,double> > all;
+  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
+  to_parse.push_back(this);
+  unsigned int i_parse = 0;
+  while (i_parse<to_parse.size()){
+    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
+    all.push_back(pair<double,double>(current->_symmetry, current->_delta_R));
+    vector<PseudoJet> 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>()){
+        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

Note: See TracChangeset for help on using the changeset viewer.

Context Navigation

Changeset b7b836a in git for external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.cc

Legend:

external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.cc

Download in other formats: