SISConeSphericalPlugin.cc
Go to the documentation of this file.00001
00002
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
00011 #include<sstream>
00012
00013 FASTJET_BEGIN_NAMESPACE
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
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
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
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
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
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;
00083
00084 if (caching()) {
00085
00086
00087
00088
00089
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
00097
00098 new_siscone |= !have_same_momentum(clust_seq.jets()[i],
00099 (*stored_particles)[i]);
00100 }
00101 }
00102 }
00103
00104
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
00117 siscone->SM_var2_hardest_cut_off = _split_merge_stopping_scale*_split_merge_stopping_scale;
00118
00119
00120
00121
00122 siscone->stable_cone_soft_E2_cutoff = ghost_separation_scale()
00123 * ghost_separation_scale();
00124
00125 siscone->set_E_weighted_splitting(_use_E_weighted_splitting);
00126
00127
00128 if (new_siscone) {
00129
00130 std::vector<CSphmomentum> siscone_momenta(n);
00131 for(unsigned i = 0; i < n; i++) {
00132 const PseudoJet & p = clust_seq.jets()[i];
00133 siscone_momenta[i] = CSphmomentum(p.px(), p.py(), p.pz(), p.E());
00134 }
00135
00136
00137
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
00143
00144
00145 siscone->recompute_jets(overlap_threshold(), protojet_or_ghost_Emin(),
00146 Esplit_merge_scale(split_merge_scale()));
00147 }
00148
00149
00150 int njet = siscone->jets.size();
00151
00152
00153 SISConeSphericalExtras * extras = new SISConeSphericalExtras(n);
00154
00155 for (int ijet = njet-1; ijet >= 0; ijet--) {
00156 const CSphjet & jet = siscone->jets[ijet];
00157
00158
00159
00160
00161 int jet_k = jet.contents[0];
00162 for (unsigned ipart = 1; ipart < jet.contents.size(); ipart++) {
00163
00164 int jet_i = jet_k;
00165
00166 int jet_j = jet.contents[ipart];
00167
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
00174 PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j];
00175
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
00183
00184
00185 double d_iB = clust_seq.jets()[jet_k].perp2();
00186 clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
00187
00188
00189 extras->_pass[clust_seq.jets()[jet_k].cluster_hist_index()] = jet.pass;
00190 }
00191
00192
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
00203 extras->_jet_def_plugin = this;
00204
00205
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