Fork me on GitHub

Ignore:
File:
1 edited

Legend:

Unmodified
Added
Removed

  • external/fastjet/contribs/RecursiveTools/RecursiveSymmetryCutBase.cc

    r1f1f858 rb7b836a  
    1 // $Id: RecursiveSymmetryCutBase.cc 700 2014-07-07 12:50:05Z gsoyez $
     1// $Id: RecursiveSymmetryCutBase.cc 1080 2017-09-28 07:51:37Z gsoyez $
    22//
    33// Copyright (c) 2014-, Gavin P. Salam, Gregory Soyez, Jesse Thaler
    … …  
    2323#include "fastjet/JetDefinition.hh"
    2424#include "fastjet/ClusterSequenceAreaBase.hh"
     25#include <cassert>
    2526#include <algorithm>
    2627#include <cstdlib>
    … …  
    4243PseudoJet RecursiveSymmetryCutBase::result(const PseudoJet & jet) const {
    4344  // construct the input jet (by default, recluster with C/A)
    44 
    4545  if (! jet.has_constituents()){
    4646    throw Error("RecursiveSymmetryCutBase can only be applied to jets with constituents");
    4747  }
    4848
    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 
     49  PseudoJet j = _recluster_if_needed(jet);
     50
     51  // sanity check: the jet must have a valid CS
    6152  if (! j.has_valid_cluster_sequence()){
    6253    throw Error("RecursiveSymmetryCutBase can only be applied to jets associated to a (valid) cluster sequence");
    6354  }
    6455
     56  // check that area information is there in case we have a subtractor
     57  // GS: do we really need this since subtraction may not require areas?
    6558  if (_subtractor) {
    6659    const ClusterSequenceAreaBase * csab =
    … …  
    7568    subjet = (*_subtractor)(subjet);
    7669  }
    77 
    78   bool use_mu_cut = (_mu_cut != numeric_limits<double>::infinity());
    7970
    8071  // variables for tracking what will happen
    … …  
    8677  std::vector<double> dropped_symmetry;
    8778  std::vector<double> dropped_mu;
     79
     80  double sym, mu2;
    8881 
    8982  // 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     }
     83  RecursionStatus status;
     84  while ((status=recurse_one_step(subjet, piece1, piece2, sym, mu2)) != recursion_success) {
     85    // start with sanity checks:
     86    if ((status == recursion_issue) || (status == recursion_no_parents)) {
     87      // we should return piece1 by our convention for recurse_one_step
     88      PseudoJet result;
     89      if (status == recursion_issue){
     90        result = piece1;
     91        if (_verbose) cout << "reached end; returning null jet " << endl;
     92      } else {
     93        result = _result_no_substructure(piece1);
     94        if (_verbose) cout << "no parents found; returning last PJ or empty jet" << endl;
     95      }
     96     
     97      if (result != 0) {
     98        // if in grooming mode, add dummy structure information
     99        StructureType * structure = new StructureType(result);
     100        // structure->_symmetry = 0.0;
     101        // structure->_mu       = 0.0;
     102        // structure->_delta_R  = 0.0;
     103        if (_verbose_structure) { // still want to store verbose information about dropped branches
     104          structure->_has_verbose = true;
     105          structure->_dropped_symmetry = dropped_symmetry;
     106          structure->_dropped_mu = dropped_mu;
     107          structure->_dropped_delta_R = dropped_delta_R;
     108        }
     109        result.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
     110      }
     111 
     112      return result;
     113    }
     114
     115    assert(status == recursion_dropped);
    174116   
    175117    // if desired, store information about dropped branches before recursing
    … …  
    179121      dropped_mu.push_back((mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2));
    180122    }
    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 
     123
     124    subjet = piece1;
     125  }
     126 
     127
     128  // we've tagged the splitting, return the jet with its substructure
     129  StructureType * structure = new StructureType(subjet);
     130  structure->_symmetry = sym;
     131  structure->_mu       = (mu2 >= 0) ? sqrt(mu2) : -sqrt(-mu2);
     132  structure->_delta_R  = sqrt(squared_geometric_distance(piece1, piece2));
     133  if (_verbose_structure) {
     134    structure->_has_verbose = true;
     135    structure->_dropped_symmetry = dropped_symmetry;
     136    structure->_dropped_mu = dropped_mu;
     137    structure->_dropped_delta_R = dropped_delta_R;
     138  }
     139  subjet.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(structure));
     140  return subjet;
     141}
     142
     143
     144 
     145//----------------------------------------------------------------------
     146// the method below is the one actually performing one step of the
     147// recursion.
     148//
     149// It returns a status code (defined above)
     150//
     151// In case of success, all the information is filled
     152// In case of "no parents", piee1 is the same subjet
     153// In case of trouble, piece2 will be a 0 PJ and piece1 is the PJ we
     154//   should return (either 0 itself if the issue was critical, or
     155//   non-wero in case of a minor issue just causing the recursion to
     156//   stop)
     157RecursiveSymmetryCutBase::RecursionStatus
     158      RecursiveSymmetryCutBase::recurse_one_step(const PseudoJet & subjet,
     159                                                 PseudoJet &piece1, PseudoJet &piece2,
     160                                                 double &sym, double &mu2,
     161                                                 void *extra_parameters) const {
     162  if (!subjet.has_parents(piece1, piece2)){
     163    piece1 = subjet;   
     164    piece2 = PseudoJet();
     165    return recursion_no_parents;
     166  }
     167
     168  // first sanity check:
     169  // - zero or negative pts are not allowed for the input subjet
     170  // - zero or negative masses are not allowed for configurations
     171  //   in which the mass will effectively appear in a denominator
     172  //   (The masses will be checked later)
     173  if (subjet.pt2() <= 0){ // this is a critical problem, return an empty PJ
     174    piece1 = piece2 = PseudoJet();
     175    return recursion_issue;
     176  }
     177
     178  if (_subtractor) {
     179    piece1 = (*_subtractor)(piece1);
     180    piece2 = (*_subtractor)(piece2);
     181  }
     182   
     183  // determine the symmetry parameter
     184  if        (_symmetry_measure == y) {
     185    // the original d_{ij}/m^2 choice from MDT
     186    // first make sure the mass is sensible
     187    if (subjet.m2() <= 0) {
     188      _negative_mass_warning.warn("RecursiveSymmetryCutBase: cannot calculate y, because (sub)jet mass is negative; bailing out");
     189      // since rounding errors can give -ve masses, be a it more
     190      // tolerant and consider that no substructure has been found
     191      piece1 = _result_no_substructure(subjet);
     192      piece2 = PseudoJet();
     193      return recursion_issue;
     194    }
     195    sym = piece1.kt_distance(piece2) / subjet.m2();
     196   
     197  } else if (_symmetry_measure == vector_z) {
     198    // min(pt1, pt2)/(pt), where the denominator is a vector sum
     199    // of the two subjets
     200    sym = min(piece1.pt(), piece2.pt()) / subjet.pt();   
     201  } else if (_symmetry_measure == scalar_z) {
     202    // min(pt1, pt2)/(pt1+pt2), where the denominator is a scalar sum
     203    // of the two subjets
     204    double pt1 = piece1.pt();
     205    double pt2 = piece2.pt();
     206    // make sure denominator is non-zero
     207    sym = pt1 + pt2;
     208    if (sym == 0){ // this is a critical problem, return an empty PJ
     209      piece1 = piece2 = PseudoJet();
     210      return recursion_issue;
     211    }
     212    sym = min(pt1, pt2) / sym;
     213  } else if ((_symmetry_measure == theta_E) || (_symmetry_measure == cos_theta_E)){
     214    // min(E1, E2)/(E1+E2)
     215    double E1 = piece1.E();
     216    double E2 = piece2.E();
     217    // make sure denominator is non-zero
     218    sym = E1 + E2;
     219    if (sym == 0){ // this is a critical problem, return an empty PJ
     220      piece1 = piece2 = PseudoJet();
     221      return recursion_issue;
     222    }
     223    sym = min(E1, E2) / sym;
     224  } else {
     225    throw Error ("Unrecognized choice of symmetry_measure");
     226  }
     227 
     228  // determine the symmetry cut
     229  // (This function is specified in the derived classes)
     230  double this_symmetry_cut = symmetry_cut_fn(piece1, piece2, extra_parameters);
     231 
     232  // and make a first tagging decision based on symmetry cut
     233  bool tagged = (sym > this_symmetry_cut);
     234
     235  // if tagged based on symmetry cut, then check the mu cut (if relevant)
     236  // and update the tagging decision. Calculate mu^2 regardless, for cases
     237  // of users not cutting on mu2, but still interested in its value.
     238  bool use_mu_cut = (_mu_cut != numeric_limits<double>::infinity());
     239  mu2 = max(piece1.m2(), piece2.m2())/subjet.m2();
     240  if (tagged && use_mu_cut) {
     241    // first a sanity check -- mu2 won't be sensible if the subjet mass
     242    // is negative, so we can't then trust the mu cut - bail out
     243    if (subjet.m2() <= 0) {
     244      _negative_mass_warning.warn("RecursiveSymmetryCutBase: cannot trust mu, because (sub)jet mass is negative; bailing out");
     245      piece1 = piece2 = PseudoJet();
     246      return recursion_issue;
     247    }
     248    if (mu2 > 1) _mu2_gt1_warning.warn("RecursiveSymmetryCutBase encountered mu^2 value > 1");
     249    if (mu2 > pow(_mu_cut,2)) tagged = false;
     250  }
     251
     252  // we'll continue unclustering, allowing for the different
     253  // ways of choosing which parent to look into
     254  if        (_recursion_choice == larger_pt) {
     255    if (piece1.pt2() < piece2.pt2()) std::swap(piece1, piece2);   
     256  } else if (_recursion_choice == larger_mt) {
     257    if (piece1.mt2() < piece2.mt2()) std::swap(piece1, piece2);   
     258  } else if (_recursion_choice == larger_m)  {
     259    if (piece1.m2() < piece2.m2()) std::swap(piece1, piece2);   
     260  } else if (_recursion_choice == larger_E)  {
     261    if (piece1.E()  < piece2.E())  std::swap(piece1, piece2);   
     262  } else {
     263    throw Error ("Unrecognized value for recursion_choice");
     264  }   
     265
     266  return tagged ? recursion_success : recursion_dropped;
     267}
     268
     269 
    225270//----------------------------------------------------------------------
    226271string RecursiveSymmetryCutBase::description() const {
    … …  
    235280  case vector_z:
    236281    ostr << "vector_z"; break;
     282  case theta_E:
     283    ostr << "theta_E"; break;
     284  case cos_theta_E:
     285    ostr << "cos_theta_E"; break;
    237286  default:
    238287    cerr << "failed to interpret symmetry_measure" << endl; exit(-1);
    … …  
    254303  case larger_m:
    255304    ostr << "mass"; break;
     305  case larger_E:
     306    ostr << "energy"; break;
    256307  default:
    257308    cerr << "failed to interpret recursion_choice" << endl; exit(-1);
    … …  
    259310
    260311  if (_subtractor) {
    261     ostr << " and subtractor: " << _subtractor->description();
     312    ostr << ", subtractor: " << _subtractor->description();
    262313    if (_input_jet_is_subtracted) {ostr << " (input jet is assumed already subtracted)";}
    263314  }
     315
     316  if (_recluster) {
     317    ostr << " and reclustering using " << _recluster->description();
     318  }
     319 
    264320  return ostr.str();
    265321}
    266322
     323//----------------------------------------------------------------------
     324// helper for handling the reclustering
     325PseudoJet RecursiveSymmetryCutBase::_recluster_if_needed(const PseudoJet &jet) const{
     326  if (! _do_reclustering) return jet;
     327  if (_recluster) return (*_recluster)(jet);
     328  if (is_ee()){
     329#if FASTJET_VERSION_NUMBER >= 30100
     330    return Recluster(JetDefinition(ee_genkt_algorithm, JetDefinition::max_allowable_R, 0.0), true)(jet);
     331#else
     332    return Recluster(JetDefinition(ee_genkt_algorithm, JetDefinition::max_allowable_R, 0.0))(jet);
     333#endif
     334  }
     335
     336  return Recluster(cambridge_algorithm, JetDefinition::max_allowable_R)(jet);
     337} 
     338 
     339//----------------------------------------------------------------------
    267340// decide what to return when no substructure has been found
     341double RecursiveSymmetryCutBase::squared_geometric_distance(const PseudoJet &j1,
     342                                                            const PseudoJet &j2) const{
     343  if (_symmetry_measure == theta_E){
     344    double dot_3d = j1.px()*j2.px() + j1.py()*j2.py() + j1.pz()*j2.pz();
     345    double cos_theta = max(-1.0,min(1.0, dot_3d/sqrt(j1.modp2()*j2.modp2())));
     346    double theta = acos(cos_theta);
     347    return theta*theta;
     348  } else if (_symmetry_measure == cos_theta_E){
     349    double dot_3d = j1.px()*j2.px() + j1.py()*j2.py() + j1.pz()*j2.pz();
     350    return max(0.0, 2*(1-dot_3d/sqrt(j1.modp2()*j2.modp2())));
     351  }
     352
     353  return j1.squared_distance(j2);
     354}
     355
     356//----------------------------------------------------------------------
    268357PseudoJet RecursiveSymmetryCutBase::_result_no_substructure(const PseudoJet &last_parent) const{
    269358  if (_grooming_mode){
    … …  
    277366
    278367
     368//========================================================================
     369// implementation of the details of the structure
     370 
     371// the number of dropped subjets
     372int RecursiveSymmetryCutBase::StructureType::dropped_count(bool global) const {
     373  check_verbose("dropped_count()");
     374
     375  // if this jet has no substructure, just return an empty vector
     376  if (!has_substructure()) return _dropped_delta_R.size();
     377
     378  // deal with the non-global case
     379  if (!global) return _dropped_delta_R.size();
     380
     381  // for the global case, we've unfolded the recursion (likely more
     382  // efficient as it requires less copying)
     383  unsigned int count = 0;
     384  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
     385  to_parse.push_back(this);
     386
     387  unsigned int i_parse = 0;
     388  while (i_parse<to_parse.size()){
     389    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
     390    count += current->_dropped_delta_R.size();
     391
     392    // check if we need to recurse deeper in the substructure
     393    //
     394    // we can have 2 situations here for the underlying structure (the
     395    // one we've wrapped around):
     396    //  - it's of the clustering type
     397    //  - it's a composite jet
     398    // only in the 2nd case do we have to recurse deeper
     399    const CompositeJetStructure *css = dynamic_cast<const CompositeJetStructure*>(current->_structure.get());
     400    if (css == 0){ ++i_parse; continue; }
     401   
     402    vector<PseudoJet> prongs = css->pieces(PseudoJet()); // argument irrelevant
     403    assert(prongs.size() == 2);
     404    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
     405      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
     406        RecursiveSymmetryCutBase::StructureType* prong_structure
     407          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
     408        if (prong_structure->has_substructure())
     409          to_parse.push_back(prong_structure);
     410      }
     411    }
     412
     413    ++i_parse;
     414  }
     415  return count;
     416}
     417
     418// the delta_R of all the dropped subjets
     419vector<double> RecursiveSymmetryCutBase::StructureType::dropped_delta_R(bool global) const {
     420  check_verbose("dropped_delta_R()");
     421
     422  // if this jet has no substructure, just return an empty vector
     423  if (!has_substructure()) return vector<double>();
     424
     425  // deal with the non-global case
     426  if (!global) return _dropped_delta_R;
     427
     428  // for the global case, we've unfolded the recursion (likely more
     429  // efficient as it requires less copying)
     430  vector<double> all_dropped;
     431  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
     432  to_parse.push_back(this);
     433
     434  unsigned int i_parse = 0;
     435  while (i_parse<to_parse.size()){
     436    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
     437    all_dropped.insert(all_dropped.end(), current->_dropped_delta_R.begin(), current->_dropped_delta_R.end());
     438
     439    // check if we need to recurse deeper in the substructure
     440    //
     441    // we can have 2 situations here for the underlying structure (the
     442    // one we've wrapped around):
     443    //  - it's of the clustering type
     444    //  - it's a composite jet
     445    // only in the 2nd case do we have to recurse deeper
     446    const CompositeJetStructure *css = dynamic_cast<const CompositeJetStructure*>(current->_structure.get());
     447    if (css == 0){ ++i_parse; continue; }
     448   
     449    vector<PseudoJet> prongs = css->pieces(PseudoJet()); // argument irrelevant
     450    assert(prongs.size() == 2);
     451    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
     452      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
     453        RecursiveSymmetryCutBase::StructureType* prong_structure
     454          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
     455        if (prong_structure->has_substructure())
     456          to_parse.push_back(prong_structure);
     457      }
     458    }
     459
     460    ++i_parse;
     461  }
     462  return all_dropped;
     463}
     464
     465// the symmetry of all the dropped subjets
     466vector<double> RecursiveSymmetryCutBase::StructureType::dropped_symmetry(bool global) const {
     467  check_verbose("dropped_symmetry()");
     468
     469  // if this jet has no substructure, just return an empty vector
     470  if (!has_substructure()) return vector<double>();
     471
     472  // deal with the non-global case
     473  if (!global) return _dropped_symmetry;
     474
     475  // for the global case, we've unfolded the recursion (likely more
     476  // efficient as it requires less copying)
     477  vector<double> all_dropped;
     478  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
     479  to_parse.push_back(this);
     480
     481  unsigned int i_parse = 0;
     482  while (i_parse<to_parse.size()){
     483    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
     484    all_dropped.insert(all_dropped.end(), current->_dropped_symmetry.begin(), current->_dropped_symmetry.end());
     485
     486    // check if we need to recurse deeper in the substructure
     487    //
     488    // we can have 2 situations here for the underlying structure (the
     489    // one we've wrapped around):
     490    //  - it's of the clustering type
     491    //  - it's a composite jet
     492    // only in the 2nd case do we have to recurse deeper
     493    const CompositeJetStructure *css = dynamic_cast<const CompositeJetStructure*>(current->_structure.get());
     494    if (css == 0){ ++i_parse; continue; }
     495   
     496    vector<PseudoJet> prongs = css->pieces(PseudoJet()); // argument irrelevant
     497    assert(prongs.size() == 2);
     498    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
     499      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
     500        RecursiveSymmetryCutBase::StructureType* prong_structure
     501          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
     502        if (prong_structure->has_substructure())
     503          to_parse.push_back(prong_structure);
     504      }
     505    }
     506
     507    ++i_parse;
     508  }
     509  return all_dropped;
     510}
     511
     512// the mu of all the dropped subjets
     513vector<double> RecursiveSymmetryCutBase::StructureType::dropped_mu(bool global) const {
     514  check_verbose("dropped_mu()");
     515
     516  // if this jet has no substructure, just return an empty vector
     517  if (!has_substructure()) return vector<double>();
     518
     519  // deal with the non-global case
     520  if (!global) return _dropped_mu;
     521
     522  // for the global case, we've unfolded the recursion (likely more
     523  // efficient as it requires less copying)
     524  vector<double> all_dropped;
     525  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
     526  to_parse.push_back(this);
     527
     528  unsigned int i_parse = 0;
     529  while (i_parse<to_parse.size()){
     530    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
     531    all_dropped.insert(all_dropped.end(), current->_dropped_mu.begin(), current->_dropped_mu.end());
     532
     533    // check if we need to recurse deeper in the substructure
     534    //
     535    // we can have 2 situations here for the underlying structure (the
     536    // one we've wrapped around):
     537    //  - it's of the clustering type
     538    //  - it's a composite jet
     539    // only in the 2nd case do we have to recurse deeper
     540    const CompositeJetStructure *css = dynamic_cast<const CompositeJetStructure*>(current->_structure.get());
     541    if (css == 0){ ++i_parse; continue; }
     542   
     543    vector<PseudoJet> prongs = css->pieces(PseudoJet()); // argument irrelevant
     544    assert(prongs.size() == 2);
     545    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
     546      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
     547        RecursiveSymmetryCutBase::StructureType* prong_structure
     548          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
     549        if (prong_structure->has_substructure())
     550          to_parse.push_back(prong_structure);
     551      }
     552    }
     553
     554    ++i_parse;
     555  }
     556  return all_dropped;
     557}
     558
     559// the maximum of the symmetry over the dropped subjets
     560double RecursiveSymmetryCutBase::StructureType::max_dropped_symmetry(bool global) const {
     561  check_verbose("max_dropped_symmetry()");
     562
     563  // if there is no substructure, just exit
     564  if (!has_substructure()){ return 0.0; }
     565
     566  // local value of the max_dropped_symmetry
     567  double local_max = (_dropped_symmetry.size() == 0)
     568    ? 0.0 : *max_element(_dropped_symmetry.begin(),_dropped_symmetry.end());
     569
     570  // recurse down the structure if instructed to do so
     571  if (global){
     572    // we can have 2 situations here for the underlying structure (the
     573    // one we've wrapped around):
     574    //  - it's of the clustering type
     575    //  - it's a composite jet
     576    // only in the 2nd case do we have to recurse deeper
     577    const CompositeJetStructure *css = dynamic_cast<const CompositeJetStructure*>(_structure.get());
     578    if (css == 0) return local_max;
     579   
     580    vector<PseudoJet> prongs = css->pieces(PseudoJet()); // argument irrelevant
     581    assert(prongs.size() == 2);
     582    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
     583      // check if the prong has further substructure
     584      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
     585        RecursiveSymmetryCutBase::StructureType* prong_structure
     586          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
     587        local_max = max(local_max, prong_structure->max_dropped_symmetry(true));
     588      }
     589    }
     590  }
     591
     592  return local_max;
     593}
     594
     595//------------------------------------------------------------------------
     596// helper class to sort by decreasing thetag
     597class SortRecursiveSoftDropStructureZgThetagPair{
     598public:
     599  bool operator()(const pair<double, double> &p1, const pair<double, double> &p2) const{
     600    return p1.second > p2.second;
     601  }
     602};
     603//------------------------------------------------------------------------
     604
     605// the (zg,thetag) pairs of all the splitting that were found and passed the SD condition
     606vector<pair<double,double> > RecursiveSymmetryCutBase::StructureType::sorted_zg_and_thetag() const {
     607  //check_verbose("sorted_zg_and_thetag()");
     608
     609  // if this jet has no substructure, just return an empty vector
     610  if (!has_substructure()) return vector<pair<double,double> >();
     611 
     612  // otherwise fill a vector with all the prongs (no specific ordering)
     613  vector<pair<double,double> > all;
     614  vector<const RecursiveSymmetryCutBase::StructureType*> to_parse;
     615  to_parse.push_back(this);
     616
     617  unsigned int i_parse = 0;
     618  while (i_parse<to_parse.size()){
     619    const RecursiveSymmetryCutBase::StructureType *current = to_parse[i_parse];
     620    all.push_back(pair<double,double>(current->_symmetry, current->_delta_R));
     621
     622    vector<PseudoJet> prongs = current->pieces(PseudoJet());
     623    assert(prongs.size() == 2);
     624    for (unsigned int i_prong=0; i_prong<2; ++i_prong){
     625      if (prongs[i_prong].has_structure_of<RecursiveSymmetryCutBase>()){
     626        RecursiveSymmetryCutBase::StructureType* prong_structure
     627          = (RecursiveSymmetryCutBase::StructureType*) prongs[i_prong].structure_ptr();
     628        if (prong_structure->has_substructure())
     629          to_parse.push_back(prong_structure);
     630      }
     631    }
     632
     633    ++i_parse;
     634  }
     635
     636  sort(all.begin(), all.end(), SortRecursiveSoftDropStructureZgThetagPair());
     637  return all;
     638}
     639
    279640} // namespace contrib
    280641
Note: See TracChangeset for help on using the changeset viewer.