SISConePlugin.cc
Go to the documentation of this file.00001
00002
00003 #include "fastjet/ClusterSequence.hh"
00004 #include "fastjet/SISConePlugin.hh"
00005
00006
00007 #include "siscone/momentum.h"
00008 #include "siscone/siscone.h"
00009
00010
00011 #include<sstream>
00012
00013 FASTJET_BEGIN_NAMESPACE
00014
00015 using namespace std;
00016 using namespace siscone;
00017
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
00026
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
00034
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
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
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
00076 void SISConePlugin::run_clustering(ClusterSequence & clust_seq) const {
00077
00078 Csiscone local_siscone;
00079 Csiscone * siscone;
00080
00081 unsigned n = clust_seq.jets().size();
00082
00083 bool new_siscone = true;
00084
00085 if (caching()) {
00086
00087
00088
00089
00090
00091 if (stored_siscone.get() != 0) {
00092 new_siscone = !(stored_plugin->cone_radius() == cone_radius()
00093 && stored_plugin->n_pass_max() == n_pass_max()
00094 && stored_particles->size() == n);
00095 if (!new_siscone) {
00096 for(unsigned i = 0; i < n; i++) {
00097
00098
00099 new_siscone |= !have_same_momentum(clust_seq.jets()[i],
00100 (*stored_particles)[i]);
00101 }
00102 }
00103 }
00104
00105
00106 if (new_siscone) {
00107 stored_siscone .reset( new Csiscone );
00108 stored_particles.reset( new std::vector<PseudoJet>(clust_seq.jets()));
00109 reset_stored_plugin();
00110 }
00111
00112 siscone = stored_siscone.get();
00113 } else {
00114 siscone = &local_siscone;
00115 }
00116
00117
00118 siscone->SM_var2_hardest_cut_off = _split_merge_stopping_scale*_split_merge_stopping_scale;
00119
00120
00121
00122
00123 siscone->stable_cone_soft_pt2_cutoff = ghost_separation_scale()
00124 * ghost_separation_scale();
00125
00126 siscone->set_pt_weighted_splitting(_use_pt_weighted_splitting);
00127
00128 if (new_siscone) {
00129
00130 std::vector<Cmomentum> siscone_momenta(n);
00131 for(unsigned i = 0; i < n; i++) {
00132 const PseudoJet & p = clust_seq.jets()[i];
00133 siscone_momenta[i] = Cmomentum(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_ptmin(),
00140 Esplit_merge_scale(split_merge_scale()));
00141 } else {
00142
00143
00144
00145 siscone->recompute_jets(overlap_threshold(), protojet_or_ghost_ptmin(),
00146 Esplit_merge_scale(split_merge_scale()));
00147 }
00148
00149
00150 int njet = siscone->jets.size();
00151
00152
00153 SISConeExtras * extras = new SISConeExtras(n);
00154
00155 for (int ijet = njet-1; ijet >= 0; ijet--) {
00156 const Cjet & 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
00210 void SISConePlugin::reset_stored_plugin() const{
00211 stored_plugin.reset( new SISConePlugin(*this));
00212 }
00213
00214 FASTJET_END_NAMESPACE