FastJet 3.4.1
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
13FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
14
15using namespace std;
16using namespace siscone_spherical;
17
18/// shortcut for converting siscone CSphmomentum into PseudoJet
19template<> 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
27namespace 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);
64 j.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(
65 new SISConeSphericalPlugin::UserScaleBaseStructureType<siscone_spherical::CSphjet>(jet,_cs)));
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/////////////////////////////////////////////
81SharedPtr<SISConeSphericalPlugin> SISConeSphericalPlugin::stored_plugin;
82SharedPtr<std::vector<PseudoJet> > SISConeSphericalPlugin::stored_particles;
83SharedPtr<CSphsiscone> SISConeSphericalPlugin::stored_siscone;
84
85
86/////////////////////////////////////////////
87// now comes the implementation itself //
88/////////////////////////////////////////////
89string 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
140void 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
297void SISConeSphericalPlugin::reset_stored_plugin() const{
298 stored_plugin.reset( new SISConeSphericalPlugin(*this));
299}
300
301
302FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
deals with clustering
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,...
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
void plugin_associate_extras(Extras *extras_in)
the plugin can associate some extra information with the ClusterSequence object by calling this funct...
static std::ostream * fastjet_banner_stream()
returns a pointer to the stream to be used to print banners (cout by default).
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],...
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition: PseudoJet.hh:68
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:392
base class for user-defined ordering of stable cones (used for prorgessive removal)
Class that provides extra information about a SISCone clustering.
Implementation of the spherical version of the SISCone algorithm (plugin for fastjet v2....
An implementation of shared pointers that is broadly similar to C++11 shared_ptr (https://en....
Definition: SharedPtr.hh:341
T * get() const
get the stored pointer
Definition: SharedPtr.hh:473
void reset()
reset the pointer to default value (NULL)
Definition: SharedPtr.hh:381
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:426