fastjet 2.4.5
|
00001 //STARTHEADER 00002 // $Id: PxConePlugin.cc 1021 2008-01-15 20:32:25Z soyez $ 00003 // 00004 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam 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, write to the Free Software 00026 // Foundation, Inc.: 00027 // 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00028 //---------------------------------------------------------------------- 00029 //ENDHEADER 00030 00031 #include "fastjet/PxConePlugin.hh" 00032 00033 #include "fastjet/ClusterSequence.hh" 00034 #include <sstream> 00035 00036 // pxcone stuff 00037 #include "pxcone.h" 00038 00039 00040 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 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 // only have hh mode 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 // max number of allowed jets 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 // run pxcone 00081 pxcone( 00082 mode , // 1=>e+e-, 2=>hadron-hadron 00083 ntrak , // Number of particles 00084 itkdm , // First dimension of PTRAK array: 00085 ptrak , // Array of particle 4-momenta (Px,Py,Pz,E) 00086 cone_radius() , // Cone size (half angle) in radians 00087 min_jet_energy() , // Minimum Jet energy (GeV) 00088 overlap_threshold() , // Maximum fraction of overlap energy in a jet 00089 mxjet , // Maximum possible number of jets 00090 njet , // Number of jets found 00091 pjet , // 5-vectors of jets 00092 ipass, // Particle k belongs to jet number IPASS(k)-1 00093 // IPASS = -1 if not assosciated to a jet 00094 ijmul, // Jet i contains IJMUL[i] particles 00095 ierr // = 0 if all is OK ; = -1 otherwise 00096 ); 00097 00098 if (ierr != 0) throw Error("An error occurred while running PXCONE"); 00099 00100 // now transfer information back 00101 valarray<int> last_index_created(njet); 00102 00103 vector<vector<int> > jet_particle_content(njet); 00104 00105 // get a list of particles in each jet 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 // now transfer the jets back into our own structure -- we will 00112 // mimic the cone code with a sequential recombination sequence in 00113 // which the jets are built up by adding one particle at a time 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 // retrieve our misappropriated index for the jet 00121 int jet_j = jet_trak_list[ilist]; 00122 // do a fake recombination step with dij=0 00123 double dij = 0.0; 00124 //clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k); 00125 if (ilist != jet_trak_list.size()-1 || E_scheme_jets()) { 00126 // our E-scheme recombination in cases where it doesn't matter 00127 clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k); 00128 } else { 00129 // put in pxcone's momentum for the last recombination so that the 00130 // final inclusive jet corresponds exactly to PXCONE's 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 // NB: put a sensible looking d_iB just to be nice... 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 //cout << endl; 00146 //for (int ijet = 0; ijet < njet; ijet++) { 00147 // PseudoJet jet(pjet[ijet][0],pjet[ijet][1],pjet[ijet][2],pjet[ijet][3]); 00148 // cout << jet.perp() << " " << jet.rap() << endl; 00149 //} 00150 //cout << "-----------------------------------------------------\n"; 00151 //vector<PseudoJet> ourjets(clust_seq.inclusive_jets()); 00152 //for (vector<PseudoJet>::const_iterator ourjet = ourjets.begin(); 00153 // ourjet != ourjets.end(); ourjet++) { 00154 // cout << ourjet->perp() << " " << ourjet->rap() << endl; 00155 //} 00157 00158 delete[] ptrak; 00159 delete[] ipass; 00160 delete[] ijmul; 00161 delete[] pjet; 00162 } 00163 00164 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh