FastJet 3.0.2
|
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