00001
00002
00003 #include "fastjet/ClusterSequence.hh"
00004 #include "SISConePlugin.hh"
00005
00006
00007 #include "momentum.h"
00008 #include "siscone.h"
00009
00010
00011 #include<sstream>
00012
00013 FASTJET_BEGIN_NAMESPACE
00014
00015 using namespace std;
00016 using namespace siscone;
00017
00018
00019 auto_ptr<SISConePlugin > SISConePlugin::stored_plugin ;
00020 auto_ptr<vector<PseudoJet> > SISConePlugin::stored_particles ;
00021 auto_ptr<Csiscone > SISConePlugin::stored_siscone ;
00022
00023 string SISConePlugin::description () const {
00024 ostringstream desc;
00025
00026 const string on = "on";
00027 const string off = "off";
00028
00029 string sm_scale_string = "split-merge uses " +
00030 split_merge_scale_name(Esplit_merge_scale(split_merge_scale()));
00031
00032 desc << "SISCone jet algorithm with " ;
00033 desc << "cone_radius = " << cone_radius () << ", ";
00034 desc << "overlap_threshold = " << overlap_threshold () << ", ";
00035 desc << "n_pass_max = " << n_pass_max () << ", ";
00036 desc << "protojet_ptmin = " << protojet_ptmin() << ", ";
00037 desc << sm_scale_string << ", ";
00038 desc << "caching turned " << (caching() ? on : off);
00039
00040
00041 Csiscone siscone;
00042 if (siscone.merge_identical_protocones) {
00043 desc << ", and (IR unsafe) merge_indentical_protocones=true" ;
00044 }
00045
00046 return desc.str();
00047 }
00048
00049
00051 template<> PseudoJet::PseudoJet(const Cmomentum & four_vector) {
00052 (*this) = PseudoJet(four_vector.px,four_vector.py,four_vector.pz,
00053 four_vector.E);
00054 }
00055
00056
00057 void SISConePlugin::run_clustering(ClusterSequence & clust_seq) const {
00058
00059
00060 Csiscone * siscone;
00061 Csiscone local_siscone;
00062
00063 unsigned n = clust_seq.jets().size();
00064
00065 bool new_siscone = true;
00066
00067 if (caching()) {
00068
00069
00070
00071
00072
00073 if (stored_siscone.get() != 0) {
00074 new_siscone = !(stored_plugin->cone_radius() == cone_radius()
00075 && stored_plugin->n_pass_max() == n_pass_max()
00076 && stored_particles->size() == n);
00077 if (!new_siscone) {
00078 for(unsigned i = 0; i < n; i++) {
00079
00080
00081 new_siscone |= !have_same_momentum(clust_seq.jets()[i],
00082 (*stored_particles)[i]);
00083 }
00084 }
00085 }
00086
00087
00088 if (new_siscone) {
00089 stored_siscone .reset( new Csiscone );
00090 stored_particles.reset( new vector<PseudoJet>(clust_seq.jets()));
00091 stored_plugin .reset( new SISConePlugin(*this) );
00092 }
00093
00094 siscone = stored_siscone.get();
00095 } else {
00096 siscone = &local_siscone;
00097 }
00098
00099 if (new_siscone) {
00100
00101 vector<Cmomentum> siscone_momenta(n);
00102 for(unsigned i = 0; i < n; i++) {
00103 const PseudoJet & p = clust_seq.jets()[i];
00104 siscone_momenta[i] = Cmomentum(p.px(), p.py(), p.pz(), p.E());
00105 }
00106
00107
00108
00109 siscone->compute_jets(siscone_momenta, cone_radius(), overlap_threshold(),
00110 n_pass_max(), protojet_ptmin(),
00111 Esplit_merge_scale(split_merge_scale()));
00112 } else {
00113
00114
00115 siscone->recompute_jets(overlap_threshold(), protojet_ptmin(),
00116 Esplit_merge_scale(split_merge_scale()));
00117 }
00118
00119
00120 int njet = siscone->jets.size();
00121
00122 for (int ijet = njet-1; ijet >= 0; ijet--) {
00123 const Cjet & jet = siscone->jets[ijet];
00124
00125
00126
00127
00128 int jet_k = jet.contents[0];
00129 for (unsigned ipart = 1; ipart < jet.contents.size(); ipart++) {
00130
00131 int jet_i = jet_k;
00132
00133 int jet_j = jet.contents[ipart];
00134
00135 double dij = 0.0;
00136
00137
00138 PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j];
00139
00140
00141 newjet.set_user_index(jet.pass);
00142
00143 clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, newjet, jet_k);
00144 }
00145
00146
00147
00148 double d_iB = clust_seq.jets()[jet_k].perp2();
00149 clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
00150 }
00151
00152
00153 SISConeExtras * extras = new SISConeExtras;
00154 for (unsigned ipass = 0; ipass < siscone->protocones_list.size(); ipass++) {
00155 for (unsigned ipc = 0; ipc < siscone->protocones_list[ipass].size(); ipc++) {
00156 double rap = siscone->protocones_list[ipass][ipc].eta;
00157 double phi = siscone->protocones_list[ipass][ipc].phi;
00158 PseudoJet protocone(cos(phi),sin(phi),sinh(rap),cosh(rap));
00159
00160 protocone.set_user_index(ipass);
00161 extras->_protocones.push_back(protocone);
00162 }
00163 }
00164 extras->_most_ambiguous_split = siscone->most_ambiguous_split;
00165
00166
00167 extras->_jet_def_plugin = this;
00168
00169
00170 clust_seq.plugin_associate_extras(auto_ptr<ClusterSequence::Extras>(extras));
00171 }
00172
00173
00175 string SISConeExtras::description() const {
00176 ostringstream ostr;
00177 ostr << "This SISCone clustering found " << protocones().size()
00178 << " stable protocones";
00179 return ostr.str();
00180 }
00181
00182
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 FASTJET_END_NAMESPACE