Fork me on GitHub

source: git/external/fastjet/plugins/SISCone/SISConePlugin.cc@ 7b48085

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 7b48085 was 273e668, checked in by Pavel Demin <pavel.demin@…>, 10 years ago

upgrade FastJet to version 3.1.0

  • Property mode set to 100644
File size: 12.9 KB
Line 
1
2// fastjet stuff
3#include "fastjet/ClusterSequence.hh"
4#include "fastjet/SISConePlugin.hh"
5
6// sisocne stuff
7#include "momentum.h"
8#include "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;
17
18/// shortcut for converting siscone Cmomentum into PseudoJet
19template<> 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
26namespace 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);
63 j.set_structure_shared_ptr(SharedPtr<PseudoJetStructureBase>(
64 new SISConePlugin::UserScaleBaseStructureType<siscone::Cjet>(jet,_cs)));
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/////////////////////////////////////////////
80std::auto_ptr<SISConePlugin> SISConePlugin::stored_plugin;
81std::auto_ptr<std::vector<PseudoJet> > SISConePlugin::stored_particles;
82std::auto_ptr<Csiscone> SISConePlugin::stored_siscone;
83
84
85/////////////////////////////////////////////
86// now comes the implementation itself //
87/////////////////////////////////////////////
88string 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
139void 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
205 SharedPtr<siscone_plugin_internal::SISConeUserScale> internal_scale;
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
297void 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
351FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
Note: See TracBrowser for help on using the repository browser.