FastJet  3.3.0
SISConeSphericalPlugin.cc
1 
2 // fastjet stuff
3 #include "fastjet/ClusterSequence.hh"
4 #include "fastjet/SISConeSphericalPlugin.hh"
5 
6 //// siscone stuff
7 #include "siscone/spherical/momentum.h"
8 #include "siscone/spherical/siscone.h"
9 
10 // other stuff
11 #include<sstream>
12 
13 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
14 
15 using namespace std;
16 using namespace siscone_spherical;
17 
18 /// shortcut for converting siscone CSphmomentum into PseudoJet
19 template<> PseudoJet::PseudoJet(const siscone_spherical::CSphmomentum & four_vector) {
20  (*this) = PseudoJet(four_vector.px,four_vector.py,four_vector.pz,
21  four_vector.E);
22 }
23 
24 
25 //======================================================================
26 // wrap-up around siscone's user-defined scales
27 namespace siscone_plugin_internal{
28  /// @ingroup internal
29  /// \class SISConeUserScale
30  /// class that makes the transition between the internal SISCone
31  /// user-defined scale choice (using SISCone's Cjet) and
32  /// user-defined scale choices in the plugn above (using FastJet's
33  /// PseudoJets)
34  class SISConeSphericalUserScale : public siscone_spherical::CSphsplit_merge::Cuser_scale_base{
35  public:
36  /// ctor takes the "fastjet-style" user-defined scale as well as a
37  /// reference to the current cluster sequence (to access the
38  /// particles if needed)
39  SISConeSphericalUserScale(const SISConeSphericalPlugin::UserScaleBase *user_scale,
40  const ClusterSequence &cs)
41  : _user_scale(user_scale), _cs(cs){}
42 
43  /// returns the scale associated to a given jet
44  virtual double operator()(const siscone_spherical::CSphjet &jet) const{
45  return _user_scale->result(_build_PJ_from_Cjet(jet));
46  }
47 
48  /// returns true id the scasle associated to jet a is larger than
49  /// the scale associated to jet b
50  virtual bool is_larger(const siscone_spherical::CSphjet &a, const siscone_spherical::CSphjet &b) const{
51  return _user_scale->is_larger(_build_PJ_from_Cjet(a), _build_PJ_from_Cjet(b));
52  }
53 
54  private:
55  /// constructs a PseudoJet from a siscone::Cjet
56  ///
57  /// Note that it is tempting to overload the PseudoJet ctor. This
58  /// would not work because down the line we need to access the
59  /// original PseudoJet through the ClusterSequence and therefore
60  /// the PseudoJet structure needs to be aware of the
61  /// ClusterSequence.
62  PseudoJet _build_PJ_from_Cjet(const siscone_spherical::CSphjet &jet) const{
63  PseudoJet j(jet.v.px, jet.v.py, jet.v.pz, jet.v.E);
66  return j;
67  }
68 
69  const SISConeSphericalPlugin::UserScaleBase *_user_scale;
70  const ClusterSequence & _cs;
71  };
72 }
73 
74 // end of the internal material
75 //======================================================================
76 
77 
78 /////////////////////////////////////////////
79 // static members declaration //
80 /////////////////////////////////////////////
81 SharedPtr<SISConeSphericalPlugin> SISConeSphericalPlugin::stored_plugin;
82 SharedPtr<std::vector<PseudoJet> > SISConeSphericalPlugin::stored_particles;
83 SharedPtr<CSphsiscone> SISConeSphericalPlugin::stored_siscone;
84 
85 
86 /////////////////////////////////////////////
87 // now comes the implementation itself //
88 /////////////////////////////////////////////
89 string SISConeSphericalPlugin::description () const {
90  ostringstream desc;
91 
92  const string on = "on";
93  const string off = "off";
94 
95  string sm_scale_string = "split-merge uses " +
96  split_merge_scale_name(Esplit_merge_scale(split_merge_scale()));
97 
98  desc << "Spherical SISCone jet algorithm with " ;
99  desc << "cone_radius = " << cone_radius () << ", ";
100  if (_progressive_removal)
101  desc << "progressive-removal mode, ";
102  else
103  desc << "overlap_threshold = " << overlap_threshold () << ", ";
104  desc << "n_pass_max = " << n_pass_max () << ", ";
105  desc << "protojet_Emin = " << protojet_Emin() << ", ";
106  if (_progressive_removal && _user_scale) {
107  desc << "using a user-defined scale for ordering of stable cones";
108  string user_scale_desc = _user_scale->description();
109  if (user_scale_desc != "") desc << " (" << user_scale_desc << ")";
110  } else {
111  desc << sm_scale_string;
112  }
113  if (!_progressive_removal){
114  desc << "caching turned " << (caching() ? on : off);
115  desc << ", SM stop scale = " << _split_merge_stopping_scale;
116  }
117 
118  // add a note to the description if we use the pt-weighted splitting
119  if (_use_E_weighted_splitting){
120  desc << ", using E-weighted splitting";
121  }
122 
123  if (_use_jet_def_recombiner){
124  desc << ", using jet-definition's own recombiner";
125  }
126 
127  // create a fake siscone object so that we can find out more about it
128  CSphsiscone siscone;
129  if (siscone.merge_identical_protocones) {
130  desc << ", and (IR unsafe) merge_indentical_protocones=true" ;
131  }
132 
133  desc << ", SISCone code v" << siscone_version();
134 
135  return desc.str();
136 }
137 
138 
139 // overloading the base class implementation
140 void SISConeSphericalPlugin::run_clustering(ClusterSequence & clust_seq) const {
141 
142  CSphsiscone::set_banner_stream(clust_seq.fastjet_banner_stream());
143 
144  CSphsiscone local_siscone;
145  CSphsiscone * siscone;
146 
147  unsigned n = clust_seq.jets().size();
148 
149  bool new_siscone = true; // by default we'll be running it
150 
151  if (caching() && !_progressive_removal) {
152 
153  // Establish if we have a cached run with the same R, npass and
154  // particles. If not then do any tidying up / reallocation that's
155  // necessary for the next round of caching, otherwise just set
156  // relevant pointers so that we can reuse and old run.
157  if (stored_siscone.get() != 0) {
158  new_siscone = !(stored_plugin->cone_radius() == cone_radius()
159  && stored_plugin->n_pass_max() == n_pass_max()
160  && stored_particles->size() == n);
161  if (!new_siscone) {
162  for(unsigned i = 0; i < n; i++) {
163  // only check momentum because indices will be correctly dealt
164  // with anyway when extracting the clustering order.
165  new_siscone |= !have_same_momentum(clust_seq.jets()[i],
166  (*stored_particles)[i]);
167  }
168  }
169  }
170 
171  // allocate the new siscone, etc., if need be
172  if (new_siscone) {
173  stored_siscone .reset( new CSphsiscone );
174  stored_particles.reset( new std::vector<PseudoJet>(clust_seq.jets()));
175  reset_stored_plugin();
176  }
177 
178  siscone = stored_siscone.get();
179  } else {
180  siscone = &local_siscone;
181  }
182 
183  // make sure stopping scale is set in siscone
184  siscone->SM_var2_hardest_cut_off = _split_merge_stopping_scale*_split_merge_stopping_scale;
185 
186  // set the specific parameters
187  // when running with ghosts for passive areas, do not put the
188  // ghosts into the stable-cone search (not relevant)
189  siscone->stable_cone_soft_E2_cutoff = ghost_separation_scale()
190  * ghost_separation_scale();
191  // set the type of splitting we want (default=std one, true->pt-weighted split)
192  siscone->set_E_weighted_splitting(_use_E_weighted_splitting);
193 
194 
195  if (new_siscone) {
196  // transfer fastjet initial particles into the siscone type
197  std::vector<CSphmomentum> siscone_momenta(n);
198  for(unsigned i = 0; i < n; i++) {
199  const PseudoJet & p = clust_seq.jets()[i]; // shorthand
200  siscone_momenta[i] = CSphmomentum(p.px(), p.py(), p.pz(), p.E());
201  }
202 
203  // run the jet finding
204  //cout << "plg sms: " << split_merge_scale() << endl;
205  if (_progressive_removal){
206  // handle the optional user-defined scale choice
208  if (_user_scale){
209  internal_scale.reset(new siscone_plugin_internal::SISConeSphericalUserScale(_user_scale, clust_seq));
210  siscone->set_user_scale(internal_scale.get());
211  }
212  siscone->compute_jets_progressive_removal(siscone_momenta, cone_radius(),
213  n_pass_max(), protojet_or_ghost_Emin(),
214  Esplit_merge_scale(split_merge_scale()));
215  } else {
216  siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(),
217  n_pass_max(), protojet_or_ghost_Emin(),
218  Esplit_merge_scale(split_merge_scale()));
219  }
220  } else {
221  // rerun the jet finding
222  // just run the overlap part of the jets.
223  //cout << "plg rcmp sms: " << split_merge_scale() << endl;
224  siscone->recompute_jets(overlap_threshold(), protojet_or_ghost_Emin(),
225  Esplit_merge_scale(split_merge_scale()));
226  }
227 
228  // extract the jets [in reverse order -- to get nice ordering in pt at end]
229  int njet = siscone->jets.size();
230 
231  // allocate space for the extras object
233 
234  // the ordering in which the inclusive jets are transfered here is
235  // deliberate and ensures that when a user asks for
236  // inclusive_jets(), they are provided in the order in which SISCone
237  // created them.
238  for (int ijet = njet-1; ijet >= 0; ijet--) {
239  const CSphjet & jet = siscone->jets[ijet]; // shorthand
240 
241  // Successively merge the particles that make up the cone jet
242  // until we have all particles in it. Start off with the zeroth
243  // particle.
244  int jet_k = jet.contents[0];
245  for (unsigned ipart = 1; ipart < jet.contents.size(); ipart++) {
246  // take the last result of the merge
247  int jet_i = jet_k;
248  // and the next element of the jet
249  int jet_j = jet.contents[ipart];
250  // and merge them (with a fake dij)
251  double dij = 0.0;
252 
253  if (_use_jet_def_recombiner) {
254  clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
255  } else {
256  // create the new jet by hand so that we can adjust its user index
257  PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j];
258  // set the user index to be the pass in which the jet was discovered
259  // *** The following line was commented for 3.0.1 ***
260  //newjet.set_user_index(jet.pass);
261  clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, newjet, jet_k);
262  }
263 
264  }
265 
266  // we have merged all the jet's particles into a single object, so now
267  // "declare" it to be a beam (inclusive) jet.
268  // [NB: put a sensible looking d_iB just to be nice...]
269  double d_iB = clust_seq.jets()[jet_k].perp2();
270  clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
271 
272  // now record the pass of the jet in the extras object
273  extras->_pass[clust_seq.jets()[jet_k].cluster_hist_index()] = jet.pass;
274  }
275 
276  // now copy the list of protocones into an "extras" objects
277  for (unsigned ipass = 0; ipass < siscone->protocones_list.size(); ipass++) {
278  for (unsigned ipc = 0; ipc < siscone->protocones_list[ipass].size(); ipc++) {
279  PseudoJet protocone(siscone->protocones_list[ipass][ipc]);
280  protocone.set_user_index(ipass);
281  extras->_protocones.push_back(protocone);
282  }
283  }
284  extras->_most_ambiguous_split = siscone->most_ambiguous_split;
285 
286  // tell it what the jet definition was
287  extras->_jet_def_plugin = this;
288 
289  // give the extras object to the cluster sequence.
290  //
291  // As of v3.1 of FastJet, extras are automatically owned (as
292  // SharedPtr) by the ClusterSequence and auto_ptr is deprecated. So
293  // we can use a simple pointer here
294  clust_seq.plugin_associate_extras(extras);
295 }
296 
297 void SISConeSphericalPlugin::reset_stored_plugin() const{
298  stored_plugin.reset( new SISConeSphericalPlugin(*this));
299 }
300 
301 
302 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
void plugin_associate_extras(Extras *extras_in)
the plugin can associate some extra information with the ClusterSequence object by calling this funct...
void set_user_index(const int index)
set the user_index, intended to allow the user to add simple identifying information to a particle/je...
Definition: PseudoJet.hh:338
deals with clustering
T * get() const
get the stored pointer
Definition: SharedPtr.hh:250
Implementation of the spherical version of the SISCone algorithm (plugin for fastjet v2...
void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, int &newjet_k)
record the fact that there has been a recombination between jets()[jet_i] and jets()[jet_k], with the specified dij, and return the index (newjet_k) allocated to the new jet, whose momentum is assumed to be the 4-vector sum of that of jet_i and jet_j
template class derived from UserScaleBase::StryctureType that works for both SISCone jet classes impl...
void plugin_record_iB_recombination(int jet_i, double diB)
record the fact that there has been a recombination between jets()[jet_i] and the beam...
an implementation of C++0x shared pointers (or boost&#39;s)
Definition: SharedPtr.hh:121
bool have_same_momentum(const PseudoJet &jeta, const PseudoJet &jetb)
returns true if the momenta of the two input jets are identical
Definition: PseudoJet.cc:336
base class for user-defined ordering of stable cones (used for prorgessive removal) ...
void set_structure_shared_ptr(const SharedPtr< PseudoJetStructureBase > &structure)
set the associated structure
Definition: PseudoJet.cc:468
static std::ostream * fastjet_banner_stream()
returns a pointer to the stream to be used to print banners (cout by default).
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
Class that provides extra information about a SISCone clustering.
void reset()
reset the pointer to default value (NULL)
Definition: SharedPtr.hh:161