PxConePlugin.cc
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "fastjet/PxConePlugin.hh"
00032
00033 #include "fastjet/ClusterSequence.hh"
00034 #include <sstream>
00035
00036
00037 #include "pxcone.h"
00038
00039
00040 FASTJET_BEGIN_NAMESPACE
00041
00042 using namespace std;
00043
00044 string PxConePlugin::description () const {
00045 ostringstream desc;
00046
00047 desc << "PxCone jet algorithm with "
00048 << "cone_radius = " << cone_radius () << ", "
00049 << "min_jet_energy = " << min_jet_energy () << ", "
00050 << "overlap_threshold = " << overlap_threshold () << ", "
00051 << "E_scheme_jets = " << E_scheme_jets ()
00052 << " (NB: non-standard version of PxCone, containing small bug fixes by Gavin Salam)";
00053
00054 return desc.str();
00055 }
00056
00057
00058 void PxConePlugin::run_clustering(ClusterSequence & clust_seq) const {
00059
00060
00061 int mode = 2;
00062
00063 int ntrak = clust_seq.jets().size(), itkdm = 4;
00064 double *ptrak = new double[ntrak*4+1];
00065 for (int i = 0; i < ntrak; i++) {
00066 ptrak[4*i+0] = clust_seq.jets()[i].px();
00067 ptrak[4*i+1] = clust_seq.jets()[i].py();
00068 ptrak[4*i+2] = clust_seq.jets()[i].pz();
00069 ptrak[4*i+3] = clust_seq.jets()[i].E();
00070 }
00071
00072
00073 int mxjet = ntrak;
00074 int njet;
00075 double *pjet = new double[mxjet*5+1];
00076 int *ipass = new int[ntrak+1];
00077 int *ijmul = new int[mxjet+1];
00078 int ierr;
00079
00080
00081 pxcone(
00082 mode ,
00083 ntrak ,
00084 itkdm ,
00085 ptrak ,
00086 cone_radius() ,
00087 min_jet_energy() ,
00088 overlap_threshold() ,
00089 mxjet ,
00090 njet ,
00091 pjet ,
00092 ipass,
00093
00094 ijmul,
00095 ierr
00096 );
00097
00098 if (ierr != 0) throw Error("An error occurred while running PXCONE");
00099
00100
00101 valarray<int> last_index_created(njet);
00102
00103 vector<vector<int> > jet_particle_content(njet);
00104
00105
00106 for (int itrak = 0; itrak < ntrak; itrak++) {
00107 int jet_i = ipass[itrak] - 1;
00108 if (jet_i >= 0) jet_particle_content[jet_i].push_back(itrak);
00109 }
00110
00111
00112
00113
00114 for(int ipxjet = njet-1; ipxjet >= 0; ipxjet--) {
00115 const vector<int> & jet_trak_list = jet_particle_content[ipxjet];
00116 int jet_k = jet_trak_list[0];
00117
00118 for (unsigned ilist = 1; ilist < jet_trak_list.size(); ilist++) {
00119 int jet_i = jet_k;
00120
00121 int jet_j = jet_trak_list[ilist];
00122
00123 double dij = 0.0;
00124
00125 if (ilist != jet_trak_list.size()-1 || E_scheme_jets()) {
00126
00127 clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
00128 } else {
00129
00130
00131 clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij,
00132 PseudoJet(pjet[5*ipxjet+0],pjet[5*ipxjet+1],
00133 pjet[5*ipxjet+2],pjet[5*ipxjet+3]),
00134 jet_k);
00135 }
00136 }
00137
00138
00139 double d_iB = clust_seq.jets()[jet_k].perp2();
00140 clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
00141 }
00142
00143
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00157
00158 delete[] ptrak;
00159 delete[] ipass;
00160 delete[] ijmul;
00161 delete[] pjet;
00162 }
00163
00164 FASTJET_END_NAMESPACE