FastJet 3.0.2
SISConeSphericalPlugin.cc
00001 
00002 // fastjet stuff
00003 #include "fastjet/ClusterSequence.hh"
00004 #include "fastjet/SISConeSphericalPlugin.hh"
00005 
00006 //// siscone stuff
00007 #include "siscone/spherical/momentum.h"
00008 #include "siscone/spherical/siscone.h"
00009 
00010 // other stuff
00011 #include<sstream>
00012 
00013 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00014 
00015 using namespace std;
00016 using namespace siscone_spherical;
00017 
00018 /// shortcut for converting siscone CSphmomentum into PseudoJet
00019 template<> PseudoJet::PseudoJet(const siscone_spherical::CSphmomentum & four_vector) {
00020   (*this) = PseudoJet(four_vector.px,four_vector.py,four_vector.pz,
00021                       four_vector.E);
00022 }
00023 
00024 /////////////////////////////////////////////
00025 // static members declaration              //
00026 /////////////////////////////////////////////
00027 std::auto_ptr<SISConeSphericalPlugin>  SISConeSphericalPlugin::stored_plugin;
00028 std::auto_ptr<std::vector<PseudoJet> > SISConeSphericalPlugin::stored_particles;
00029 std::auto_ptr<CSphsiscone>             SISConeSphericalPlugin::stored_siscone;
00030 
00031 
00032 /////////////////////////////////////////////
00033 // now comes the implementation itself     //
00034 /////////////////////////////////////////////
00035 string SISConeSphericalPlugin::description () const {
00036   ostringstream desc;
00037   
00038   const string on = "on";
00039   const string off = "off";
00040 
00041   string sm_scale_string = "split-merge uses " + 
00042     split_merge_scale_name(Esplit_merge_scale(split_merge_scale()));
00043 
00044   desc << "Spherical SISCone jet algorithm with " ;
00045   desc << "cone_radius = "       << cone_radius        () << ", ";
00046   desc << "overlap_threshold = " << overlap_threshold  () << ", ";
00047   desc << "n_pass_max = "        << n_pass_max         () << ", ";
00048   desc << "protojet_Emin = "     << protojet_Emin()      << ", ";
00049   desc <<  sm_scale_string                                << ", ";
00050   desc << "caching turned "      << (caching() ? on : off);
00051   desc << ", SM stop scale = "     << _split_merge_stopping_scale;
00052 
00053   // add a note to the description if we use the pt-weighted splitting
00054   if (_use_E_weighted_splitting){
00055     desc << ", using E-weighted splitting";
00056   }
00057 
00058   if (_use_jet_def_recombiner){
00059     desc << ", using jet-definition's own recombiner";
00060   }
00061 
00062   // create a fake siscone object so that we can find out more about it
00063   CSphsiscone siscone;
00064   if (siscone.merge_identical_protocones) {
00065     desc << ", and (IR unsafe) merge_indentical_protocones=true" ;
00066   }
00067 
00068   desc << ", SISCone code v" << siscone_version();
00069 
00070   return desc.str();
00071 }
00072 
00073 
00074 // overloading the base class implementation
00075 void SISConeSphericalPlugin::run_clustering(ClusterSequence & clust_seq) const {
00076 
00077   CSphsiscone::set_banner_stream(clust_seq.fastjet_banner_stream());
00078 
00079   CSphsiscone   local_siscone;
00080   CSphsiscone * siscone;
00081 
00082   unsigned n = clust_seq.jets().size();
00083 
00084   bool new_siscone = true; // by default we'll be running it
00085 
00086   if (caching()) {
00087 
00088     // Establish if we have a cached run with the same R, npass and
00089     // particles. If not then do any tidying up / reallocation that's
00090     // necessary for the next round of caching, otherwise just set
00091     // relevant pointers so that we can reuse and old run.
00092     if (stored_siscone.get() != 0) {
00093       new_siscone = !(stored_plugin->cone_radius()   == cone_radius()
00094                       && stored_plugin->n_pass_max() == n_pass_max()
00095                       && stored_particles->size()    == n);
00096       if (!new_siscone) {
00097         for(unsigned i = 0; i < n; i++) {
00098           // only check momentum because indices will be correctly dealt
00099           // with anyway when extracting the clustering order.
00100           new_siscone |= !have_same_momentum(clust_seq.jets()[i],
00101                                              (*stored_particles)[i]);
00102         }
00103       }
00104     }
00105 
00106     // allocate the new siscone, etc., if need be
00107     if (new_siscone) {
00108       stored_siscone  .reset( new CSphsiscone );
00109       stored_particles.reset( new std::vector<PseudoJet>(clust_seq.jets()));
00110       reset_stored_plugin();
00111     }
00112 
00113     siscone = stored_siscone.get();
00114   } else {
00115     siscone = &local_siscone;
00116   }
00117 
00118   // make sure stopping scale is set in siscone
00119   siscone->SM_var2_hardest_cut_off = _split_merge_stopping_scale*_split_merge_stopping_scale;
00120 
00121   // set the specific parameters
00122   // when running with ghosts for passive areas, do not put the
00123   // ghosts into the stable-cone search (not relevant)
00124   siscone->stable_cone_soft_E2_cutoff = ghost_separation_scale()
00125                                       * ghost_separation_scale();
00126   // set the type of splitting we want (default=std one, true->pt-weighted split)
00127   siscone->set_E_weighted_splitting(_use_E_weighted_splitting);
00128 
00129 
00130   if (new_siscone) {
00131     // transfer fastjet initial particles into the siscone type
00132     std::vector<CSphmomentum> siscone_momenta(n);
00133     for(unsigned i = 0; i < n; i++) {
00134       const PseudoJet & p = clust_seq.jets()[i]; // shorthand
00135       siscone_momenta[i] = CSphmomentum(p.px(), p.py(), p.pz(), p.E());
00136     }
00137 
00138     // run the jet finding
00139     //cout << "plg sms: " << split_merge_scale() << endl;
00140     siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(),
00141                           n_pass_max(), protojet_or_ghost_Emin(), 
00142                           Esplit_merge_scale(split_merge_scale()));
00143   } else {
00144     // rerun the jet finding
00145     // just run the overlap part of the jets.
00146     //cout << "plg rcmp sms: " << split_merge_scale() << endl;
00147     siscone->recompute_jets(overlap_threshold(), protojet_or_ghost_Emin(), 
00148                             Esplit_merge_scale(split_merge_scale()));
00149   }
00150 
00151   // extract the jets [in reverse order -- to get nice ordering in pt at end]
00152   int njet = siscone->jets.size();
00153 
00154   // allocate space for the extras object
00155   SISConeSphericalExtras * extras = new SISConeSphericalExtras(n);
00156 
00157   for (int ijet = njet-1; ijet >= 0; ijet--) {
00158     const CSphjet & jet = siscone->jets[ijet]; // shorthand
00159 
00160     // Successively merge the particles that make up the cone jet
00161     // until we have all particles in it.  Start off with the zeroth
00162     // particle.
00163     int jet_k = jet.contents[0];
00164     for (unsigned ipart = 1; ipart < jet.contents.size(); ipart++) {
00165       // take the last result of the merge
00166       int jet_i = jet_k;
00167       // and the next element of the jet
00168       int jet_j = jet.contents[ipart];
00169       // and merge them (with a fake dij)
00170       double dij = 0.0;
00171 
00172       if (_use_jet_def_recombiner) {
00173        clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
00174       } else {
00175        // create the new jet by hand so that we can adjust its user index
00176        PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j];
00177        // set the user index to be the pass in which the jet was discovered
00178        // *** The following line was commented for 3.0.1 ***
00179        //newjet.set_user_index(jet.pass);
00180        clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, newjet, jet_k);
00181       }
00182 
00183     }
00184 
00185     // we have merged all the jet's particles into a single object, so now
00186     // "declare" it to be a beam (inclusive) jet.
00187     // [NB: put a sensible looking d_iB just to be nice...]
00188     double d_iB = clust_seq.jets()[jet_k].perp2();
00189     clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
00190 
00191     // now record the pass of the jet in the extras object
00192     extras->_pass[clust_seq.jets()[jet_k].cluster_hist_index()] = jet.pass;
00193   }
00194 
00195   // now copy the list of protocones into an "extras" objects
00196   for (unsigned ipass = 0; ipass < siscone->protocones_list.size(); ipass++) {
00197     for (unsigned ipc = 0; ipc < siscone->protocones_list[ipass].size(); ipc++) {
00198       PseudoJet protocone(siscone->protocones_list[ipass][ipc]);
00199       protocone.set_user_index(ipass);
00200       extras->_protocones.push_back(protocone);
00201     }
00202   }
00203   extras->_most_ambiguous_split = siscone->most_ambiguous_split;
00204 
00205   // tell it what the jet definition was
00206   extras->_jet_def_plugin = this;
00207 
00208   // give the extras object to the cluster sequence.
00209   clust_seq.plugin_associate_extras(std::auto_ptr<ClusterSequence::Extras>(extras));
00210 }
00211 
00212 void SISConeSphericalPlugin::reset_stored_plugin() const{
00213   stored_plugin.reset( new SISConeSphericalPlugin(*this));
00214 }
00215 
00216 
00217 FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends