FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: PxConePlugin.cc 2777 2011-11-25 15:01:56Z soyez $ 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 bool PxConePlugin::_first_time = true; 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 // print a banner if we run this for the first time 00060 //_print_banner(clust_seq.fastjet_banner_stream()); 00061 00062 // only have hh mode 00063 int mode = 2; 00064 00065 int ntrak = clust_seq.jets().size(), itkdm = 4; 00066 double *ptrak = new double[ntrak*4+1]; 00067 for (int i = 0; i < ntrak; i++) { 00068 ptrak[4*i+0] = clust_seq.jets()[i].px(); 00069 ptrak[4*i+1] = clust_seq.jets()[i].py(); 00070 ptrak[4*i+2] = clust_seq.jets()[i].pz(); 00071 ptrak[4*i+3] = clust_seq.jets()[i].E(); 00072 } 00073 00074 // max number of allowed jets 00075 int mxjet = ntrak; 00076 int njet; 00077 double *pjet = new double[mxjet*5+1]; 00078 int *ipass = new int[ntrak+1]; 00079 int *ijmul = new int[mxjet+1]; 00080 int ierr; 00081 00082 // run pxcone 00083 pxcone( 00084 mode , // 1=>e+e-, 2=>hadron-hadron 00085 ntrak , // Number of particles 00086 itkdm , // First dimension of PTRAK array: 00087 ptrak , // Array of particle 4-momenta (Px,Py,Pz,E) 00088 cone_radius() , // Cone size (half angle) in radians 00089 min_jet_energy() , // Minimum Jet energy (GeV) 00090 overlap_threshold() , // Maximum fraction of overlap energy in a jet 00091 mxjet , // Maximum possible number of jets 00092 njet , // Number of jets found 00093 pjet , // 5-vectors of jets 00094 ipass, // Particle k belongs to jet number IPASS(k)-1 00095 // IPASS = -1 if not assosciated to a jet 00096 ijmul, // Jet i contains IJMUL[i] particles 00097 ierr // = 0 if all is OK ; = -1 otherwise 00098 ); 00099 00100 if (ierr != 0) throw Error("An error occurred while running PXCONE"); 00101 00102 // now transfer information back 00103 valarray<int> last_index_created(njet); 00104 00105 vector<vector<int> > jet_particle_content(njet); 00106 00107 // get a list of particles in each jet 00108 for (int itrak = 0; itrak < ntrak; itrak++) { 00109 int jet_i = ipass[itrak] - 1; 00110 if (jet_i >= 0) jet_particle_content[jet_i].push_back(itrak); 00111 } 00112 00113 // now transfer the jets back into our own structure -- we will 00114 // mimic the cone code with a sequential recombination sequence in 00115 // which the jets are built up by adding one particle at a time 00116 for(int ipxjet = njet-1; ipxjet >= 0; ipxjet--) { 00117 const vector<int> & jet_trak_list = jet_particle_content[ipxjet]; 00118 int jet_k = jet_trak_list[0]; 00119 00120 for (unsigned ilist = 1; ilist < jet_trak_list.size(); ilist++) { 00121 int jet_i = jet_k; 00122 // retrieve our misappropriated index for the jet 00123 int jet_j = jet_trak_list[ilist]; 00124 // do a fake recombination step with dij=0 00125 double dij = 0.0; 00126 //clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k); 00127 if (ilist != jet_trak_list.size()-1 || E_scheme_jets()) { 00128 // our E-scheme recombination in cases where it doesn't matter 00129 clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k); 00130 } else { 00131 // put in pxcone's momentum for the last recombination so that the 00132 // final inclusive jet corresponds exactly to PXCONE's 00133 clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, 00134 PseudoJet(pjet[5*ipxjet+0],pjet[5*ipxjet+1], 00135 pjet[5*ipxjet+2],pjet[5*ipxjet+3]), 00136 jet_k); 00137 } 00138 } 00139 00140 // NB: put a sensible looking d_iB just to be nice... 00141 double d_iB = clust_seq.jets()[jet_k].perp2(); 00142 clust_seq.plugin_record_iB_recombination(jet_k, d_iB); 00143 } 00144 00145 00146 //// following code is for testing only 00147 //cout << endl; 00148 //for (int ijet = 0; ijet < njet; ijet++) { 00149 // PseudoJet jet(pjet[ijet][0],pjet[ijet][1],pjet[ijet][2],pjet[ijet][3]); 00150 // cout << jet.perp() << " " << jet.rap() << endl; 00151 //} 00152 //cout << "-----------------------------------------------------\n"; 00153 //vector<PseudoJet> ourjets(clust_seq.inclusive_jets()); 00154 //for (vector<PseudoJet>::const_iterator ourjet = ourjets.begin(); 00155 // ourjet != ourjets.end(); ourjet++) { 00156 // cout << ourjet->perp() << " " << ourjet->rap() << endl; 00157 //} 00158 ////cout << endl; 00159 00160 delete[] ptrak; 00161 delete[] ipass; 00162 delete[] ijmul; 00163 delete[] pjet; 00164 } 00165 00166 // print a banner for reference to the 3rd-party code 00167 void PxConePlugin::_print_banner(ostream *ostr) const{ 00168 if (! _first_time) return; 00169 _first_time=false; 00170 00171 // make sure the user has not set the banner stream to NULL 00172 if (!ostr) return; 00173 00174 (*ostr) << "#-------------------------------------------------------------------------" << endl; 00175 (*ostr) << "# You are running the PxCone plugin for FastJet " << endl; 00176 (*ostr) << "# Original code by the Luis Del Pozo, David Ward and Michael H. Seymour " << endl; 00177 (*ostr) << "# If you use this plugin, please cite " << endl; 00178 (*ostr) << "# M. H. Seymour and C. Tevlin, JHEP 0611 (2006) 052 [hep-ph/0609100]. " << endl; 00179 (*ostr) << "# in addition to the usual FastJet reference. " << endl; 00180 (*ostr) << "#-------------------------------------------------------------------------" << endl; 00181 00182 // make sure we really have the output done. 00183 ostr->flush(); 00184 } 00185 00186 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh