fastjet 2.4.5
SISConeSphericalPlugin.cc
Go to the documentation of this file.
00001 
00002 // fastjet stuff
00003 #include "fastjet/ClusterSequence.hh"
00004 #include "fastjet/SISConeSphericalPlugin.hh"
00005 
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 
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 
00025 // static members declaration              //
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 
00033 // now comes the implementation itself     //
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   local_siscone;
00078   CSphsiscone * siscone;
00079 
00080   unsigned n = clust_seq.jets().size();
00081 
00082   bool new_siscone = true; // by default we'll be running it
00083 
00084   if (caching()) {
00085 
00086     // Establish if we have a cached run with the same R, npass and
00087     // particles. If not then do any tidying up / reallocation that's
00088     // necessary for the next round of caching, otherwise just set
00089     // relevant pointers so that we can reuse and old run.
00090     if (stored_siscone.get() != 0) {
00091       new_siscone = !(stored_plugin->cone_radius()   == cone_radius()
00092                       && stored_plugin->n_pass_max() == n_pass_max()
00093                       && stored_particles->size()    == n);
00094       if (!new_siscone) {
00095         for(unsigned i = 0; i < n; i++) {
00096           // only check momentum because indices will be correctly dealt
00097           // with anyway when extracting the clustering order.
00098           new_siscone |= !have_same_momentum(clust_seq.jets()[i],
00099                                              (*stored_particles)[i]);
00100         }
00101       }
00102     }
00103 
00104     // allocate the new siscone, etc., if need be
00105     if (new_siscone) {
00106       stored_siscone  .reset( new CSphsiscone );
00107       stored_particles.reset( new std::vector<PseudoJet>(clust_seq.jets()));
00108       reset_stored_plugin();
00109     }
00110 
00111     siscone = stored_siscone.get();
00112   } else {
00113     siscone = &local_siscone;
00114   }
00115 
00116   // make sure stopping scale is set in siscone
00117   siscone->SM_var2_hardest_cut_off = _split_merge_stopping_scale*_split_merge_stopping_scale;
00118 
00119   // set the specific parameters
00120   // when running with ghosts for passive areas, do not put the
00121   // ghosts into the stable-cone search (not relevant)
00122   siscone->stable_cone_soft_E2_cutoff = ghost_separation_scale()
00123                                       * ghost_separation_scale();
00124   // set the type of splitting we want (default=std one, true->pt-weighted split)
00125   siscone->set_E_weighted_splitting(_use_E_weighted_splitting);
00126 
00127 
00128   if (new_siscone) {
00129     // transfer fastjet initial particles into the siscone type
00130     std::vector<CSphmomentum> siscone_momenta(n);
00131     for(unsigned i = 0; i < n; i++) {
00132       const PseudoJet & p = clust_seq.jets()[i]; // shorthand
00133       siscone_momenta[i] = CSphmomentum(p.px(), p.py(), p.pz(), p.E());
00134     }
00135 
00136     // run the jet finding
00137     //cout << "plg sms: " << split_merge_scale() << endl;
00138     siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(),
00139                           n_pass_max(), protojet_or_ghost_Emin(), 
00140                           Esplit_merge_scale(split_merge_scale()));
00141   } else {
00142     // rerun the jet finding
00143     // just run the overlap part of the jets.
00144     //cout << "plg rcmp sms: " << split_merge_scale() << endl;
00145     siscone->recompute_jets(overlap_threshold(), protojet_or_ghost_Emin(), 
00146                             Esplit_merge_scale(split_merge_scale()));
00147   }
00148 
00149   // extract the jets [in reverse order -- to get nice ordering in pt at end]
00150   int njet = siscone->jets.size();
00151 
00152   // allocate space for the extras object
00153   SISConeSphericalExtras * extras = new SISConeSphericalExtras(n);
00154 
00155   for (int ijet = njet-1; ijet >= 0; ijet--) {
00156     const CSphjet & jet = siscone->jets[ijet]; // shorthand
00157 
00158     // Successively merge the particles that make up the cone jet
00159     // until we have all particles in it.  Start off with the zeroth
00160     // particle.
00161     int jet_k = jet.contents[0];
00162     for (unsigned ipart = 1; ipart < jet.contents.size(); ipart++) {
00163       // take the last result of the merge
00164       int jet_i = jet_k;
00165       // and the next element of the jet
00166       int jet_j = jet.contents[ipart];
00167       // and merge them (with a fake dij)
00168       double dij = 0.0;
00169 
00170       if (_use_jet_def_recombiner) {
00171        clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
00172       } else {
00173        // create the new jet by hand so that we can adjust its user index
00174        PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j];
00175        // set the user index to be the pass in which the jet was discovered
00176        newjet.set_user_index(jet.pass);
00177        clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, newjet, jet_k);
00178       }
00179 
00180     }
00181 
00182     // we have merged all the jet's particles into a single object, so now
00183     // "declare" it to be a beam (inclusive) jet.
00184     // [NB: put a sensible looking d_iB just to be nice...]
00185     double d_iB = clust_seq.jets()[jet_k].perp2();
00186     clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
00187 
00188     // now record the pass of the jet in the extras object
00189     extras->_pass[clust_seq.jets()[jet_k].cluster_hist_index()] = jet.pass;
00190   }
00191 
00192   // now copy the list of protocones into an "extras" objects
00193   for (unsigned ipass = 0; ipass < siscone->protocones_list.size(); ipass++) {
00194     for (unsigned ipc = 0; ipc < siscone->protocones_list[ipass].size(); ipc++) {
00195       PseudoJet protocone(siscone->protocones_list[ipass][ipc]);
00196       protocone.set_user_index(ipass);
00197       extras->_protocones.push_back(protocone);
00198     }
00199   }
00200   extras->_most_ambiguous_split = siscone->most_ambiguous_split;
00201 
00202   // tell it what the jet definition was
00203   extras->_jet_def_plugin = this;
00204 
00205   // give the extras object to the cluster sequence.
00206   clust_seq.plugin_associate_extras(std::auto_ptr<ClusterSequence::Extras>(extras));
00207 }
00208 
00209 void SISConeSphericalPlugin::reset_stored_plugin() const{
00210   stored_plugin.reset( new SISConeSphericalPlugin(*this));
00211 }
00212 
00213 
00214 FASTJET_END_NAMESPACE      // defined in fastjet/internal/base.hh
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines