FastJet  3.3.0
SISConePlugin.cc
1 
2 // fastjet stuff
3 #include "fastjet/ClusterSequence.hh"
4 #include "fastjet/SISConePlugin.hh"
5 
6 // sisocne stuff
7 #include "siscone/momentum.h"
8 #include "siscone/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;
17 
18 /// shortcut for converting siscone Cmomentum into PseudoJet
19 template<> PseudoJet::PseudoJet(const siscone::Cmomentum & four_vector) {
20  (*this) = PseudoJet(four_vector.px,four_vector.py,four_vector.pz,
21  four_vector.E);
22 }
23 
24 //======================================================================
25 // wrap-up around siscone's user-defined scales
26 namespace siscone_plugin_internal{
27  /// @ingroup internal
28  /// \class SISConeUserScale
29  /// class that makes the transition between the internal SISCone
30  /// user-defined scale choice (using SISCone's Cjet) and
31  /// user-defined scale choices in the plugn above (using FastJet's
32  /// PseudoJets)
33  class SISConeUserScale : public siscone::Csplit_merge::Cuser_scale_base{
34  public:
35  /// ctor takes the "fastjet-style" user-defined scale as well as a
36  /// reference to the current cluster sequence (to access the
37  /// particles if needed)
38  SISConeUserScale(const SISConePlugin::UserScaleBase *user_scale,
39  const ClusterSequence &cs)
40  : _user_scale(user_scale), _cs(cs){}
41 
42  /// returns the scale associated to a given jet
43  virtual double operator()(const siscone::Cjet &jet) const{
44  return _user_scale->result(_build_PJ_from_Cjet(jet));
45  }
46 
47  /// returns true id the scasle associated to jet a is larger than
48  /// the scale associated to jet b
49  virtual bool is_larger(const siscone::Cjet &a, const siscone::Cjet &b) const{
50  return _user_scale->is_larger(_build_PJ_from_Cjet(a), _build_PJ_from_Cjet(b));
51  }
52 
53  private:
54  /// constructs a PseudoJet from a siscone::Cjet
55  ///
56  /// Note that it is tempting to overload the PseudoJet ctor. This
57  /// would not work because down the line we need to access the
58  /// original PseudoJet through the ClusterSequence and therefore
59  /// the PseudoJet structure needs to be aware of the
60  /// ClusterSequence.
61  PseudoJet _build_PJ_from_Cjet(const siscone::Cjet &jet) const{
62  PseudoJet j(jet.v.px, jet.v.py, jet.v.pz, jet.v.E);
65  return j;
66  }
67 
68  const SISConePlugin::UserScaleBase *_user_scale;
69  const ClusterSequence & _cs;
70  };
71 }
72 
73 // end of the internal material
74 //======================================================================
75 
76 
77 /////////////////////////////////////////////
78 // static members declaration //
79 /////////////////////////////////////////////
80 SharedPtr<SISConePlugin> SISConePlugin::stored_plugin;
81 SharedPtr<std::vector<PseudoJet> > SISConePlugin::stored_particles;
82 SharedPtr<Csiscone> SISConePlugin::stored_siscone;
83 
84 
85 /////////////////////////////////////////////
86 // now comes the implementation itself //
87 /////////////////////////////////////////////
88 string SISConePlugin::description () const {
89  ostringstream desc;
90 
91  const string on = "on";
92  const string off = "off";
93 
94  string sm_scale_string = "split-merge uses " +
95  split_merge_scale_name(Esplit_merge_scale(split_merge_scale()));
96 
97  desc << "SISCone jet algorithm with " ;
98  desc << "cone_radius = " << cone_radius () << ", ";
99  if (_progressive_removal)
100  desc << "progressive-removal mode, ";
101  else
102  desc << "overlap_threshold = " << overlap_threshold () << ", ";
103  desc << "n_pass_max = " << n_pass_max () << ", ";
104  desc << "protojet_ptmin = " << protojet_ptmin() << ", ";
105  if (_progressive_removal && _user_scale) {
106  desc << "using a user-defined scale for ordering of stable cones";
107  string user_scale_desc = _user_scale->description();
108  if (user_scale_desc != "") desc << " (" << user_scale_desc << ")";
109  } else {
110  desc << sm_scale_string;
111  }
112  if (!_progressive_removal){
113  desc << ", caching turned " << (caching() ? on : off);
114  desc << ", SM stop scale = " << _split_merge_stopping_scale;
115  }
116 
117  // add a note to the description if we use the pt-weighted splitting
118  if (_use_pt_weighted_splitting){
119  desc << ", using pt-weighted splitting";
120  }
121 
122  if (_use_jet_def_recombiner){
123  desc << ", using jet-definition's own recombiner";
124  }
125 
126  // create a fake siscone object so that we can find out more about it
127  Csiscone siscone;
128  if (siscone.merge_identical_protocones) {
129  desc << ", and (IR unsafe) merge_indentical_protocones=true" ;
130  }
131 
132  desc << ", SISCone code v" << siscone_version();
133 
134  return desc.str();
135 }
136 
137 
138 // overloading the base class implementation
139 void SISConePlugin::run_clustering(ClusterSequence & clust_seq) const {
140 
141  Csiscone::set_banner_stream(clust_seq.fastjet_banner_stream());
142 
143  Csiscone local_siscone;
144  Csiscone * siscone;
145 
146  unsigned n = clust_seq.jets().size();
147 
148  bool new_siscone = true; // by default we'll be running it
149 
150  if (caching() && !_progressive_removal) {
151 
152  // Establish if we have a cached run with the same R, npass and
153  // particles. If not then do any tidying up / reallocation that's
154  // necessary for the next round of caching, otherwise just set
155  // relevant pointers so that we can reuse and old run.
156  if (stored_siscone.get() != 0) {
157  new_siscone = !(stored_plugin->cone_radius() == cone_radius()
158  && stored_plugin->n_pass_max() == n_pass_max()
159  && stored_particles->size() == n);
160  if (!new_siscone) {
161  for(unsigned i = 0; i < n; i++) {
162  // only check momentum because indices will be correctly dealt
163  // with anyway when extracting the clustering order.
164  new_siscone |= !have_same_momentum(clust_seq.jets()[i],
165  (*stored_particles)[i]);
166  }
167  }
168  }
169 
170  // allocate the new siscone, etc., if need be
171  if (new_siscone) {
172  stored_siscone .reset( new Csiscone );
173  stored_particles.reset( new std::vector<PseudoJet>(clust_seq.jets()));
174  reset_stored_plugin();
175  }
176 
177  siscone = stored_siscone.get();
178  } else {
179  siscone = &local_siscone;
180  }
181 
182  // make sure stopping scale is set in siscone
183  siscone->SM_var2_hardest_cut_off = _split_merge_stopping_scale*_split_merge_stopping_scale;
184 
185  // set the specific parameters
186  // when running with ghosts for passive areas, do not put the
187  // ghosts into the stable-cone search (not relevant)
188  siscone->stable_cone_soft_pt2_cutoff = ghost_separation_scale()
189  * ghost_separation_scale();
190  // set the type of splitting we want (default=std one, true->pt-weighted split)
191  siscone->set_pt_weighted_splitting(_use_pt_weighted_splitting);
192 
193  if (new_siscone) {
194  // transfer fastjet initial particles into the siscone type
195  std::vector<Cmomentum> siscone_momenta(n);
196  for(unsigned i = 0; i < n; i++) {
197  const PseudoJet & p = clust_seq.jets()[i]; // shorthand
198  siscone_momenta[i] = Cmomentum(p.px(), p.py(), p.pz(), p.E());
199  }
200 
201  // run the jet finding
202  //cout << "plg sms: " << split_merge_scale() << endl;
203  if (_progressive_removal){
204  // handle the optional user-defined scale choice
206  if (_user_scale){
207  internal_scale.reset(new siscone_plugin_internal::SISConeUserScale(_user_scale, clust_seq));
208  siscone->set_user_scale(internal_scale.get());
209  }
210  siscone->compute_jets_progressive_removal(siscone_momenta, cone_radius(),
211  n_pass_max(), protojet_or_ghost_ptmin(),
212  Esplit_merge_scale(split_merge_scale()));
213  } else {
214  siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(),
215  n_pass_max(), protojet_or_ghost_ptmin(),
216  Esplit_merge_scale(split_merge_scale()));
217  }
218  } else {
219  // rerun the jet finding
220  // just run the overlap part of the jets.
221  //cout << "plg rcmp sms: " << split_merge_scale() << endl;
222  siscone->recompute_jets(overlap_threshold(), protojet_or_ghost_ptmin(),
223  Esplit_merge_scale(split_merge_scale()));
224  }
225 
226  // extract the jets [in reverse order -- to get nice ordering in pt at end]
227  int njet = siscone->jets.size();
228 
229  // allocate space for the extras object
230  SISConeExtras * extras = new SISConeExtras(n);
231 
232  // the ordering in which the inclusive jets are transfered here is
233  // deliberate and ensures that when a user asks for
234  // inclusive_jets(), they are provided in the order in which SISCone
235  // created them.
236  for (int ijet = njet-1; ijet >= 0; ijet--) {
237  const Cjet & jet = siscone->jets[ijet]; // shorthand
238 
239  // Successively merge the particles that make up the cone jet
240  // until we have all particles in it. Start off with the zeroth
241  // particle.
242  int jet_k = jet.contents[0];
243  for (unsigned ipart = 1; ipart < jet.contents.size(); ipart++) {
244  // take the last result of the merge
245  int jet_i = jet_k;
246  // and the next element of the jet
247  int jet_j = jet.contents[ipart];
248  // and merge them (with a fake dij)
249  double dij = 0.0;
250 
251  if (_use_jet_def_recombiner) {
252  clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
253  } else {
254  // create the new jet by hand so that we can adjust its user index
255  PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j];
256  // set the user index to be the pass in which the jet was discovered
257  // *** The following line was commented for 3.0.1 ***
258  //newjet.set_user_index(jet.pass);
259  clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, newjet, jet_k);
260  }
261 
262  }
263 
264  // we have merged all the jet's particles into a single object, so now
265  // "declare" it to be a beam (inclusive) jet.
266  // [NB: put a sensible looking d_iB just to be nice...]
267  double d_iB = clust_seq.jets()[jet_k].perp2();
268  clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
269 
270  // now record the pass of the jet in the extras object
271  extras->_pass[clust_seq.jets()[jet_k].cluster_hist_index()] = jet.pass;
272  }
273 
274  // now copy the list of protocones into an "extras" objects
275  for (unsigned ipass = 0; ipass < siscone->protocones_list.size(); ipass++) {
276  for (unsigned ipc = 0; ipc < siscone->protocones_list[ipass].size(); ipc++) {
277  PseudoJet protocone(siscone->protocones_list[ipass][ipc]);
278  protocone.set_user_index(ipass);
279  extras->_protocones.push_back(protocone);
280  }
281  }
282  extras->_most_ambiguous_split = siscone->most_ambiguous_split;
283 
284  // tell it what the jet definition was
285  extras->_jet_def_plugin = this;
286 
287  // give the extras object to the cluster sequence.
288  //
289  // As of v3.1 of FastJet, extras are automatically owned (as
290  // SharedPtr) by the ClusterSequence and auto_ptr is deprecated. So
291  // we can use a simple pointer here
292  //clust_seq.plugin_associate_extras(std::auto_ptr<ClusterSequence::Extras>(extras));
293  clust_seq.plugin_associate_extras(extras);
294 }
295 
296 
297 void SISConePlugin::reset_stored_plugin() const{
298  stored_plugin.reset( new SISConePlugin(*this));
299 }
300 
301 // //======================================================================
302 // // material to handle user-defined scales
303 //
304 // //--------------------------------------------------
305 // // SISCone structure type
306 //
307 // // the textual descripotion
308 // std::string SISConePlugin::UserScaleBase::StructureType::description() const{
309 // return "PseudoJet wrapping a siscone::Cjet stable cone";
310 // }
311 //
312 // // retrieve the constituents
313 // //
314 // // if you simply need to iterate over the constituents, it will be
315 // // faster to access them via constituent(i)
316 // vector<PseudoJet> SISConePlugin::UserScaleBase::StructureType::constituents(const PseudoJet &) const{
317 // vector<PseudoJet> constits;
318 // constits.reserve(_jet.n);
319 // for (unsigned int i=0; i<(unsigned int)_jet.n;i++)
320 // constits.push_back(constituent(i));
321 // return constits;
322 // }
323 //
324 // // returns the number of constituents
325 // unsigned int SISConePlugin::UserScaleBase::StructureType::size() const{
326 // return _jet.n;
327 // }
328 //
329 // // returns the index (in the original particle list) of the ith
330 // // constituent
331 // int SISConePlugin::UserScaleBase::StructureType::constituent_index(unsigned int i) const{
332 // return _jet.contents[i];
333 // }
334 //
335 // // returns the ith constituent (as a PseusoJet)
336 // const PseudoJet & SISConePlugin::UserScaleBase::StructureType::constituent(unsigned int i) const{
337 // return _cs.jets()[_jet.contents[i]];
338 // }
339 //
340 // // returns the scalar pt of this stable cone
341 // double SISConePlugin::UserScaleBase::StructureType::pt_tilde() const{
342 // return _jet.pt_tilde;
343 // }
344 //
345 // // returns the sm_var2 (signed ordering variable squared) for this stable cone
346 // double SISConePlugin::UserScaleBase::StructureType::ordering_var2() const{
347 // return _jet.sm_var2;
348 // }
349 
350 
351 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
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 that provides extra information about a SISCone clustering.
Implementation of the SISCone algorithm (plugin for fastjet v2.1 upwards)
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...
void reset()
reset the pointer to default value (NULL)
Definition: SharedPtr.hh:161