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 | |
---|
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 | ///////////////////////////////////////////// |
---|
26 | // static members declaration // |
---|
27 | ///////////////////////////////////////////// |
---|
28 | std::auto_ptr<SISConePlugin> SISConePlugin::stored_plugin; |
---|
29 | std::auto_ptr<std::vector<PseudoJet> > SISConePlugin::stored_particles; |
---|
30 | std::auto_ptr<Csiscone> SISConePlugin::stored_siscone; |
---|
31 | |
---|
32 | |
---|
33 | ///////////////////////////////////////////// |
---|
34 | // now comes the implementation itself // |
---|
35 | ///////////////////////////////////////////// |
---|
36 | string SISConePlugin::description () const { |
---|
37 | ostringstream desc; |
---|
38 | |
---|
39 | const string on = "on"; |
---|
40 | const string off = "off"; |
---|
41 | |
---|
42 | string sm_scale_string = "split-merge uses " + |
---|
43 | split_merge_scale_name(Esplit_merge_scale(split_merge_scale())); |
---|
44 | |
---|
45 | desc << "SISCone jet algorithm with " ; |
---|
46 | desc << "cone_radius = " << cone_radius () << ", "; |
---|
47 | desc << "overlap_threshold = " << overlap_threshold () << ", "; |
---|
48 | desc << "n_pass_max = " << n_pass_max () << ", "; |
---|
49 | desc << "protojet_ptmin = " << protojet_ptmin() << ", "; |
---|
50 | desc << sm_scale_string << ", "; |
---|
51 | desc << "caching turned " << (caching() ? on : off); |
---|
52 | desc << ", SM stop scale = " << _split_merge_stopping_scale; |
---|
53 | |
---|
54 | // add a note to the description if we use the pt-weighted splitting |
---|
55 | if (_use_pt_weighted_splitting){ |
---|
56 | desc << ", using pt-weighted splitting"; |
---|
57 | } |
---|
58 | |
---|
59 | if (_use_jet_def_recombiner){ |
---|
60 | desc << ", using jet-definition's own recombiner"; |
---|
61 | } |
---|
62 | |
---|
63 | // create a fake siscone object so that we can find out more about it |
---|
64 | Csiscone siscone; |
---|
65 | if (siscone.merge_identical_protocones) { |
---|
66 | desc << ", and (IR unsafe) merge_indentical_protocones=true" ; |
---|
67 | } |
---|
68 | |
---|
69 | desc << ", SISCone code v" << siscone_version(); |
---|
70 | |
---|
71 | return desc.str(); |
---|
72 | } |
---|
73 | |
---|
74 | |
---|
75 | // overloading the base class implementation |
---|
76 | void SISConePlugin::run_clustering(ClusterSequence & clust_seq) const { |
---|
77 | |
---|
78 | Csiscone::set_banner_stream(clust_seq.fastjet_banner_stream()); |
---|
79 | |
---|
80 | Csiscone local_siscone; |
---|
81 | Csiscone * siscone; |
---|
82 | |
---|
83 | unsigned n = clust_seq.jets().size(); |
---|
84 | |
---|
85 | bool new_siscone = true; // by default we'll be running it |
---|
86 | |
---|
87 | if (caching()) { |
---|
88 | |
---|
89 | // Establish if we have a cached run with the same R, npass and |
---|
90 | // particles. If not then do any tidying up / reallocation that's |
---|
91 | // necessary for the next round of caching, otherwise just set |
---|
92 | // relevant pointers so that we can reuse and old run. |
---|
93 | if (stored_siscone.get() != 0) { |
---|
94 | new_siscone = !(stored_plugin->cone_radius() == cone_radius() |
---|
95 | && stored_plugin->n_pass_max() == n_pass_max() |
---|
96 | && stored_particles->size() == n); |
---|
97 | if (!new_siscone) { |
---|
98 | for(unsigned i = 0; i < n; i++) { |
---|
99 | // only check momentum because indices will be correctly dealt |
---|
100 | // with anyway when extracting the clustering order. |
---|
101 | new_siscone |= !have_same_momentum(clust_seq.jets()[i], |
---|
102 | (*stored_particles)[i]); |
---|
103 | } |
---|
104 | } |
---|
105 | } |
---|
106 | |
---|
107 | // allocate the new siscone, etc., if need be |
---|
108 | if (new_siscone) { |
---|
109 | stored_siscone .reset( new Csiscone ); |
---|
110 | stored_particles.reset( new std::vector<PseudoJet>(clust_seq.jets())); |
---|
111 | reset_stored_plugin(); |
---|
112 | } |
---|
113 | |
---|
114 | siscone = stored_siscone.get(); |
---|
115 | } else { |
---|
116 | siscone = &local_siscone; |
---|
117 | } |
---|
118 | |
---|
119 | // make sure stopping scale is set in siscone |
---|
120 | siscone->SM_var2_hardest_cut_off = _split_merge_stopping_scale*_split_merge_stopping_scale; |
---|
121 | |
---|
122 | // set the specific parameters |
---|
123 | // when running with ghosts for passive areas, do not put the |
---|
124 | // ghosts into the stable-cone search (not relevant) |
---|
125 | siscone->stable_cone_soft_pt2_cutoff = ghost_separation_scale() |
---|
126 | * ghost_separation_scale(); |
---|
127 | // set the type of splitting we want (default=std one, true->pt-weighted split) |
---|
128 | siscone->set_pt_weighted_splitting(_use_pt_weighted_splitting); |
---|
129 | |
---|
130 | if (new_siscone) { |
---|
131 | // transfer fastjet initial particles into the siscone type |
---|
132 | std::vector<Cmomentum> siscone_momenta(n); |
---|
133 | for(unsigned i = 0; i < n; i++) { |
---|
134 | const PseudoJet & p = clust_seq.jets()[i]; // shorthand |
---|
135 | siscone_momenta[i] = Cmomentum(p.px(), p.py(), p.pz(), p.E()); |
---|
136 | } |
---|
137 | |
---|
138 | // run the jet finding |
---|
139 | //cout << "plg sms: " << split_merge_scale() << endl; |
---|
140 | siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(), |
---|
141 | n_pass_max(), protojet_or_ghost_ptmin(), |
---|
142 | Esplit_merge_scale(split_merge_scale())); |
---|
143 | } else { |
---|
144 | // rerun the jet finding |
---|
145 | // just run the overlap part of the jets. |
---|
146 | //cout << "plg rcmp sms: " << split_merge_scale() << endl; |
---|
147 | siscone->recompute_jets(overlap_threshold(), protojet_or_ghost_ptmin(), |
---|
148 | Esplit_merge_scale(split_merge_scale())); |
---|
149 | } |
---|
150 | |
---|
151 | // extract the jets [in reverse order -- to get nice ordering in pt at end] |
---|
152 | int njet = siscone->jets.size(); |
---|
153 | |
---|
154 | // allocate space for the extras object |
---|
155 | SISConeExtras * extras = new SISConeExtras(n); |
---|
156 | |
---|
157 | for (int ijet = njet-1; ijet >= 0; ijet--) { |
---|
158 | const Cjet & jet = siscone->jets[ijet]; // shorthand |
---|
159 | |
---|
160 | // Successively merge the particles that make up the cone jet |
---|
161 | // until we have all particles in it. Start off with the zeroth |
---|
162 | // particle. |
---|
163 | int jet_k = jet.contents[0]; |
---|
164 | for (unsigned ipart = 1; ipart < jet.contents.size(); ipart++) { |
---|
165 | // take the last result of the merge |
---|
166 | int jet_i = jet_k; |
---|
167 | // and the next element of the jet |
---|
168 | int jet_j = jet.contents[ipart]; |
---|
169 | // and merge them (with a fake dij) |
---|
170 | double dij = 0.0; |
---|
171 | |
---|
172 | if (_use_jet_def_recombiner) { |
---|
173 | clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k); |
---|
174 | } else { |
---|
175 | // create the new jet by hand so that we can adjust its user index |
---|
176 | PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j]; |
---|
177 | // set the user index to be the pass in which the jet was discovered |
---|
178 | // *** The following line was commented for 3.0.1 *** |
---|
179 | //newjet.set_user_index(jet.pass); |
---|
180 | clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, newjet, jet_k); |
---|
181 | } |
---|
182 | |
---|
183 | } |
---|
184 | |
---|
185 | // we have merged all the jet's particles into a single object, so now |
---|
186 | // "declare" it to be a beam (inclusive) jet. |
---|
187 | // [NB: put a sensible looking d_iB just to be nice...] |
---|
188 | double d_iB = clust_seq.jets()[jet_k].perp2(); |
---|
189 | clust_seq.plugin_record_iB_recombination(jet_k, d_iB); |
---|
190 | |
---|
191 | // now record the pass of the jet in the extras object |
---|
192 | extras->_pass[clust_seq.jets()[jet_k].cluster_hist_index()] = jet.pass; |
---|
193 | } |
---|
194 | |
---|
195 | // now copy the list of protocones into an "extras" objects |
---|
196 | for (unsigned ipass = 0; ipass < siscone->protocones_list.size(); ipass++) { |
---|
197 | for (unsigned ipc = 0; ipc < siscone->protocones_list[ipass].size(); ipc++) { |
---|
198 | PseudoJet protocone(siscone->protocones_list[ipass][ipc]); |
---|
199 | protocone.set_user_index(ipass); |
---|
200 | extras->_protocones.push_back(protocone); |
---|
201 | } |
---|
202 | } |
---|
203 | extras->_most_ambiguous_split = siscone->most_ambiguous_split; |
---|
204 | |
---|
205 | // tell it what the jet definition was |
---|
206 | extras->_jet_def_plugin = this; |
---|
207 | |
---|
208 | // give the extras object to the cluster sequence. |
---|
209 | clust_seq.plugin_associate_extras(std::auto_ptr<ClusterSequence::Extras>(extras)); |
---|
210 | } |
---|
211 | |
---|
212 | |
---|
213 | void SISConePlugin::reset_stored_plugin() const{ |
---|
214 | stored_plugin.reset( new SISConePlugin(*this)); |
---|
215 | } |
---|
216 | |
---|
217 | FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh |
---|