FastJet 3.0.0
|
00001 //STARTHEADER 00002 // $Id: PxConePlugin.cc 2577 2011-09-13 15:11:38Z salam $ 00003 // 00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00005 // 00006 //---------------------------------------------------------------------- 00007 // This file is part of FastJet. 00008 // 00009 // FastJet is free software; you can redistribute it and/or modify 00010 // it under the terms of the GNU General Public License as published by 00011 // the Free Software Foundation; either version 2 of the License, or 00012 // (at your option) any later version. 00013 // 00014 // The algorithms that underlie FastJet have required considerable 00015 // development and are described in hep-ph/0512210. If you use 00016 // FastJet as part of work towards a scientific publication, please 00017 // include a citation to the FastJet paper. 00018 // 00019 // FastJet is distributed in the hope that it will be useful, 00020 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00021 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00022 // GNU General Public License for more details. 00023 // 00024 // You should have received a copy of the GNU General Public License 00025 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00026 //---------------------------------------------------------------------- 00027 //ENDHEADER 00028 00029 #include "fastjet/PxConePlugin.hh" 00030 00031 #include "fastjet/ClusterSequence.hh" 00032 #include <sstream> 00033 00034 // pxcone stuff 00035 #include "pxcone.h" 00036 00037 00038 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00039 00040 using namespace std; 00041 00042 string PxConePlugin::description () const { 00043 ostringstream desc; 00044 00045 desc << "PxCone jet algorithm with " 00046 << "cone_radius = " << cone_radius () << ", " 00047 << "min_jet_energy = " << min_jet_energy () << ", " 00048 << "overlap_threshold = " << overlap_threshold () << ", " 00049 << "E_scheme_jets = " << E_scheme_jets () 00050 << " (NB: non-standard version of PxCone, containing small bug fixes by Gavin Salam)"; 00051 00052 return desc.str(); 00053 } 00054 00055 00056 void PxConePlugin::run_clustering(ClusterSequence & clust_seq) const { 00057 00058 // only have hh mode 00059 int mode = 2; 00060 00061 int ntrak = clust_seq.jets().size(), itkdm = 4; 00062 double *ptrak = new double[ntrak*4+1]; 00063 for (int i = 0; i < ntrak; i++) { 00064 ptrak[4*i+0] = clust_seq.jets()[i].px(); 00065 ptrak[4*i+1] = clust_seq.jets()[i].py(); 00066 ptrak[4*i+2] = clust_seq.jets()[i].pz(); 00067 ptrak[4*i+3] = clust_seq.jets()[i].E(); 00068 } 00069 00070 // max number of allowed jets 00071 int mxjet = ntrak; 00072 int njet; 00073 double *pjet = new double[mxjet*5+1]; 00074 int *ipass = new int[ntrak+1]; 00075 int *ijmul = new int[mxjet+1]; 00076 int ierr; 00077 00078 // run pxcone 00079 pxcone( 00080 mode , // 1=>e+e-, 2=>hadron-hadron 00081 ntrak , // Number of particles 00082 itkdm , // First dimension of PTRAK array: 00083 ptrak , // Array of particle 4-momenta (Px,Py,Pz,E) 00084 cone_radius() , // Cone size (half angle) in radians 00085 min_jet_energy() , // Minimum Jet energy (GeV) 00086 overlap_threshold() , // Maximum fraction of overlap energy in a jet 00087 mxjet , // Maximum possible number of jets 00088 njet , // Number of jets found 00089 pjet , // 5-vectors of jets 00090 ipass, // Particle k belongs to jet number IPASS(k)-1 00091 // IPASS = -1 if not assosciated to a jet 00092 ijmul, // Jet i contains IJMUL[i] particles 00093 ierr // = 0 if all is OK ; = -1 otherwise 00094 ); 00095 00096 if (ierr != 0) throw Error("An error occurred while running PXCONE"); 00097 00098 // now transfer information back 00099 valarray<int> last_index_created(njet); 00100 00101 vector<vector<int> > jet_particle_content(njet); 00102 00103 // get a list of particles in each jet 00104 for (int itrak = 0; itrak < ntrak; itrak++) { 00105 int jet_i = ipass[itrak] - 1; 00106 if (jet_i >= 0) jet_particle_content[jet_i].push_back(itrak); 00107 } 00108 00109 // now transfer the jets back into our own structure -- we will 00110 // mimic the cone code with a sequential recombination sequence in 00111 // which the jets are built up by adding one particle at a time 00112 for(int ipxjet = njet-1; ipxjet >= 0; ipxjet--) { 00113 const vector<int> & jet_trak_list = jet_particle_content[ipxjet]; 00114 int jet_k = jet_trak_list[0]; 00115 00116 for (unsigned ilist = 1; ilist < jet_trak_list.size(); ilist++) { 00117 int jet_i = jet_k; 00118 // retrieve our misappropriated index for the jet 00119 int jet_j = jet_trak_list[ilist]; 00120 // do a fake recombination step with dij=0 00121 double dij = 0.0; 00122 //clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k); 00123 if (ilist != jet_trak_list.size()-1 || E_scheme_jets()) { 00124 // our E-scheme recombination in cases where it doesn't matter 00125 clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k); 00126 } else { 00127 // put in pxcone's momentum for the last recombination so that the 00128 // final inclusive jet corresponds exactly to PXCONE's 00129 clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, 00130 PseudoJet(pjet[5*ipxjet+0],pjet[5*ipxjet+1], 00131 pjet[5*ipxjet+2],pjet[5*ipxjet+3]), 00132 jet_k); 00133 } 00134 } 00135 00136 // NB: put a sensible looking d_iB just to be nice... 00137 double d_iB = clust_seq.jets()[jet_k].perp2(); 00138 clust_seq.plugin_record_iB_recombination(jet_k, d_iB); 00139 } 00140 00141 00142 //// following code is for testing only 00143 //cout << endl; 00144 //for (int ijet = 0; ijet < njet; ijet++) { 00145 // PseudoJet jet(pjet[ijet][0],pjet[ijet][1],pjet[ijet][2],pjet[ijet][3]); 00146 // cout << jet.perp() << " " << jet.rap() << endl; 00147 //} 00148 //cout << "-----------------------------------------------------\n"; 00149 //vector<PseudoJet> ourjets(clust_seq.inclusive_jets()); 00150 //for (vector<PseudoJet>::const_iterator ourjet = ourjets.begin(); 00151 // ourjet != ourjets.end(); ourjet++) { 00152 // cout << ourjet->perp() << " " << ourjet->rap() << endl; 00153 //} 00154 ////cout << endl; 00155 00156 delete[] ptrak; 00157 delete[] ipass; 00158 delete[] ijmul; 00159 delete[] pjet; 00160 } 00161 00162 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh