| [d7d2da3] | 1 |
|
|---|
| 2 | // fastjet stuff
|
|---|
| 3 | #include "fastjet/ClusterSequence.hh"
|
|---|
| 4 | #include "fastjet/SISConePlugin.hh"
|
|---|
| 5 |
|
|---|
| 6 | // sisocne stuff
|
|---|
| [fc39685] | 7 | #include "momentum.h"
|
|---|
| 8 | #include "siscone.h"
|
|---|
| [d7d2da3] | 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 |
|
|---|
| [273e668] | 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);
|
|---|
| 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 |
|
|---|
| [d7d2da3] | 76 |
|
|---|
| 77 | /////////////////////////////////////////////
|
|---|
| 78 | // static members declaration //
|
|---|
| 79 | /////////////////////////////////////////////
|
|---|
| [1d208a2] | 80 | SharedPtr<SISConePlugin> SISConePlugin::stored_plugin;
|
|---|
| 81 | SharedPtr<std::vector<PseudoJet> > SISConePlugin::stored_particles;
|
|---|
| 82 | SharedPtr<Csiscone> SISConePlugin::stored_siscone;
|
|---|
| [d7d2da3] | 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 () << ", ";
|
|---|
| [273e668] | 99 | if (_progressive_removal)
|
|---|
| 100 | desc << "progressive-removal mode, ";
|
|---|
| 101 | else
|
|---|
| 102 | desc << "overlap_threshold = " << overlap_threshold () << ", ";
|
|---|
| [d7d2da3] | 103 | desc << "n_pass_max = " << n_pass_max () << ", ";
|
|---|
| 104 | desc << "protojet_ptmin = " << protojet_ptmin() << ", ";
|
|---|
| [273e668] | 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 | }
|
|---|
| [d7d2da3] | 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 |
|
|---|
| [273e668] | 150 | if (caching() && !_progressive_removal) {
|
|---|
| [d7d2da3] | 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;
|
|---|
| [273e668] | 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 | }
|
|---|
| [d7d2da3] | 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 |
|
|---|
| [273e668] | 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.
|
|---|
| [d7d2da3] | 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.
|
|---|
| [35cdc46] | 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);
|
|---|
| [d7d2da3] | 294 | }
|
|---|
| 295 |
|
|---|
| 296 |
|
|---|
| 297 | void SISConePlugin::reset_stored_plugin() const{
|
|---|
| 298 | stored_plugin.reset( new SISConePlugin(*this));
|
|---|
| 299 | }
|
|---|
| 300 |
|
|---|
| [273e668] | 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 |
|
|---|
| [d7d2da3] | 351 | FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
|
|---|
