FastJet 3.0.2
|
00001 00002 // fastjet stuff 00003 #include "fastjet/ClusterSequence.hh" 00004 #include "fastjet/SISConePlugin.hh" 00005 00006 // sisocne stuff 00007 #include "siscone/momentum.h" 00008 #include "siscone/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; 00017 00018 /// shortcut for converting siscone Cmomentum into PseudoJet 00019 template<> PseudoJet::PseudoJet(const siscone::Cmomentum & four_vector) { 00020 (*this) = PseudoJet(four_vector.px,four_vector.py,four_vector.pz, 00021 four_vector.E); 00022 } 00023 00024 00025 ///////////////////////////////////////////// 00026 // static members declaration // 00027 ///////////////////////////////////////////// 00028 std::auto_ptr<SISConePlugin> SISConePlugin::stored_plugin; 00029 std::auto_ptr<std::vector<PseudoJet> > SISConePlugin::stored_particles; 00030 std::auto_ptr<Csiscone> SISConePlugin::stored_siscone; 00031 00032 00033 ///////////////////////////////////////////// 00034 // now comes the implementation itself // 00035 ///////////////////////////////////////////// 00036 string SISConePlugin::description () const { 00037 ostringstream desc; 00038 00039 const string on = "on"; 00040 const string off = "off"; 00041 00042 string sm_scale_string = "split-merge uses " + 00043 split_merge_scale_name(Esplit_merge_scale(split_merge_scale())); 00044 00045 desc << "SISCone jet algorithm with " ; 00046 desc << "cone_radius = " << cone_radius () << ", "; 00047 desc << "overlap_threshold = " << overlap_threshold () << ", "; 00048 desc << "n_pass_max = " << n_pass_max () << ", "; 00049 desc << "protojet_ptmin = " << protojet_ptmin() << ", "; 00050 desc << sm_scale_string << ", "; 00051 desc << "caching turned " << (caching() ? on : off); 00052 desc << ", SM stop scale = " << _split_merge_stopping_scale; 00053 00054 // add a note to the description if we use the pt-weighted splitting 00055 if (_use_pt_weighted_splitting){ 00056 desc << ", using pt-weighted splitting"; 00057 } 00058 00059 if (_use_jet_def_recombiner){ 00060 desc << ", using jet-definition's own recombiner"; 00061 } 00062 00063 // create a fake siscone object so that we can find out more about it 00064 Csiscone siscone; 00065 if (siscone.merge_identical_protocones) { 00066 desc << ", and (IR unsafe) merge_indentical_protocones=true" ; 00067 } 00068 00069 desc << ", SISCone code v" << siscone_version(); 00070 00071 return desc.str(); 00072 } 00073 00074 00075 // overloading the base class implementation 00076 void SISConePlugin::run_clustering(ClusterSequence & clust_seq) const { 00077 00078 Csiscone::set_banner_stream(clust_seq.fastjet_banner_stream()); 00079 00080 Csiscone local_siscone; 00081 Csiscone * siscone; 00082 00083 unsigned n = clust_seq.jets().size(); 00084 00085 bool new_siscone = true; // by default we'll be running it 00086 00087 if (caching()) { 00088 00089 // Establish if we have a cached run with the same R, npass and 00090 // particles. If not then do any tidying up / reallocation that's 00091 // necessary for the next round of caching, otherwise just set 00092 // relevant pointers so that we can reuse and old run. 00093 if (stored_siscone.get() != 0) { 00094 new_siscone = !(stored_plugin->cone_radius() == cone_radius() 00095 && stored_plugin->n_pass_max() == n_pass_max() 00096 && stored_particles->size() == n); 00097 if (!new_siscone) { 00098 for(unsigned i = 0; i < n; i++) { 00099 // only check momentum because indices will be correctly dealt 00100 // with anyway when extracting the clustering order. 00101 new_siscone |= !have_same_momentum(clust_seq.jets()[i], 00102 (*stored_particles)[i]); 00103 } 00104 } 00105 } 00106 00107 // allocate the new siscone, etc., if need be 00108 if (new_siscone) { 00109 stored_siscone .reset( new Csiscone ); 00110 stored_particles.reset( new std::vector<PseudoJet>(clust_seq.jets())); 00111 reset_stored_plugin(); 00112 } 00113 00114 siscone = stored_siscone.get(); 00115 } else { 00116 siscone = &local_siscone; 00117 } 00118 00119 // make sure stopping scale is set in siscone 00120 siscone->SM_var2_hardest_cut_off = _split_merge_stopping_scale*_split_merge_stopping_scale; 00121 00122 // set the specific parameters 00123 // when running with ghosts for passive areas, do not put the 00124 // ghosts into the stable-cone search (not relevant) 00125 siscone->stable_cone_soft_pt2_cutoff = ghost_separation_scale() 00126 * ghost_separation_scale(); 00127 // set the type of splitting we want (default=std one, true->pt-weighted split) 00128 siscone->set_pt_weighted_splitting(_use_pt_weighted_splitting); 00129 00130 if (new_siscone) { 00131 // transfer fastjet initial particles into the siscone type 00132 std::vector<Cmomentum> siscone_momenta(n); 00133 for(unsigned i = 0; i < n; i++) { 00134 const PseudoJet & p = clust_seq.jets()[i]; // shorthand 00135 siscone_momenta[i] = Cmomentum(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_ptmin(), 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_ptmin(), 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 SISConeExtras * extras = new SISConeExtras(n); 00156 00157 for (int ijet = njet-1; ijet >= 0; ijet--) { 00158 const Cjet & 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 00213 void SISConePlugin::reset_stored_plugin() const{ 00214 stored_plugin.reset( new SISConePlugin(*this)); 00215 } 00216 00217 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh