FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: CMSIterativeConePlugin.cc 1504 2009-04-10 13:39:48Z salam $ 00003 // 00004 // Copyright (c) 2007-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00005 // Copyright (c) ????-????, CMS [for the iterative-cone code itself] 00006 // 00007 //---------------------------------------------------------------------- 00008 // This file is part of FastJet. It contains code that has been 00009 // obtained from the CMS collaboration, revision 1.14 of the 00010 // CMSIterativeConeAlgorithm.cc file in CMSSW, see 00011 // http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/RecoJets/JetAlgorithms/src/CMSIterativeConeAlgorithm.cc?hideattic=0&revision=1.14&view=markup 00012 // 00013 // Permission has been granted by the CMS collaboration to release it 00014 // in FastJet under the terms of the GNU Public License(v2) (see the 00015 // COPYING file in the main FastJet directory for details). 00016 // Changes from the original file are listed below. 00017 // 00018 // FastJet is free software; you can redistribute it and/or modify 00019 // it under the terms of the GNU General Public License as published by 00020 // the Free Software Foundation; either version 2 of the License, or 00021 // (at your option) any later version. 00022 // 00023 // The algorithms that underlie FastJet have required considerable 00024 // development and are described in hep-ph/0512210. If you use 00025 // FastJet as part of work towards a scientific publication, please 00026 // include a citation to the FastJet paper. 00027 // 00028 // FastJet is distributed in the hope that it will be useful, 00029 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00030 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00031 // GNU General Public License for more details. 00032 // 00033 // You should have received a copy of the GNU General Public License 00034 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00035 //---------------------------------------------------------------------- 00036 //ENDHEADER 00037 00038 // List of changes compared to the original CMS code (revision 1.14 of 00039 // CMSIterativeConeAlgorithm.cc) 00040 // 00041 // 2009-05-10 Gavin Salam <salam@lpthe.jussieu.fr> 00042 // 00043 // * added radius and seed threshold information in the plugin 00044 // description 00045 // 00046 // 2009-01-06 Gregory Soyez <soyez@fastjet.fr> 00047 // 00048 // * Encapsulated the CMS code into a plugin for FastJet 00049 // * inserted the deltaPhi and deltaR2 codes from 00050 // DataFormats/Math/interface/deltaPhi.h (rev 1.1) 00051 // DataFormats/Math/interface/deltaR.h (rev 1.2) 00052 // * Adapted the code to use PseusoJet rather than 'InputItem' 00053 // and 'InputCollection' 00054 // * use the FastJet clustering history structures instead of 00055 // the ProtoJet one used by CMS. 00056 00057 00058 // fastjet stuff 00059 #include "fastjet/ClusterSequence.hh" 00060 #include "fastjet/CMSIterativeConePlugin.hh" 00061 00062 // other stuff 00063 #include <vector> 00064 #include <list> 00065 #include <sstream> 00066 #include "SortByEt.h" 00067 00068 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00069 00070 using namespace std; 00071 using namespace cms; 00072 00073 //------------------------------------------------------ 00074 // some tools 00075 //------------------------------------------------------ 00076 template <class T> 00077 T deltaPhi (T phi1, T phi2) { 00078 T result = phi1 - phi2; 00079 while (result > M_PI) result -= 2*M_PI; 00080 while (result <= -M_PI) result += 2*M_PI; 00081 return result; 00082 } 00083 00084 template <class T> 00085 T deltaR2 (T eta1, T phi1, T eta2, T phi2) { 00086 T deta = eta1 - eta2; 00087 T dphi = deltaPhi (phi1, phi2); 00088 return deta*deta + dphi*dphi; 00089 } 00090 00091 //------------------------------------------------------ 00092 bool CMSIterativeConePlugin::_first_time = true; 00093 00094 string CMSIterativeConePlugin::description () const { 00095 ostringstream desc; 00096 desc << "CMSIterativeCone plugin with R = " << theConeRadius << " and seed threshold = " << theSeedThreshold; 00097 return desc.str(); 00098 } 00099 00100 void CMSIterativeConePlugin::run_clustering(ClusterSequence & clust_seq) const { 00101 // print a banner if we run this for the first time 00102 _print_banner(clust_seq.fastjet_banner_stream()); 00103 00104 //make a list of input objects ordered by ET 00105 //cout << "copying the list of particles" << endl; 00106 list<PseudoJet> input; 00107 for (unsigned int i=0 ; i<clust_seq.jets().size() ; i++) { 00108 input.push_back(clust_seq.jets()[i]); 00109 } 00110 NumericSafeGreaterByEt<PseudoJet> compCandidate; 00111 //cout << "sorting" << endl; 00112 input.sort(compCandidate); 00113 00114 //find jets 00115 //cout << "launching the main loop" << endl; 00116 while( !input.empty() && (input.front().Et() > theSeedThreshold )) { 00117 //cone centre 00118 double eta0=input.front().eta(); 00119 double phi0=input.front().phi(); 00120 //protojet properties 00121 double eta=0; 00122 double phi=0; 00123 double et=0; 00124 //list of towers in cone 00125 list< list<PseudoJet>::iterator> cone; 00126 for(int iteration=0;iteration<100;iteration++){ 00127 //cout << "iterating" << endl; 00128 cone.clear(); 00129 eta=0; 00130 phi=0; 00131 et=0; 00132 for(list<PseudoJet>::iterator inp=input.begin(); 00133 inp!=input.end();inp++){ 00134 const PseudoJet tower = *inp; 00135 if( deltaR2(eta0,phi0,tower.eta(),tower.phi()) < 00136 theConeRadius*theConeRadius) { 00137 double tower_et = tower.Et(); 00138 cone.push_back(inp); 00139 eta+= tower_et*tower.eta(); 00140 double dphi=tower.phi()-phi0; 00141 if(dphi>M_PI) dphi-=2*M_PI; 00142 else if(dphi<=-M_PI) dphi+=2*M_PI; 00143 phi+=tower_et*dphi; 00144 et +=tower_et; 00145 } 00146 } 00147 eta=eta/et; 00148 phi=phi0+phi/et; 00149 if(phi>M_PI)phi-=2*M_PI; 00150 else if(phi<=-M_PI)phi+=2*M_PI; 00151 00152 if(fabs(eta-eta0)<.001 && fabs(phi-phi0)<.001) break;//stable cone found 00153 eta0=eta; 00154 phi0=phi; 00155 } 00156 00157 //cout << "make the jet final" << endl; 00158 00159 //make a final jet and remove the jet constituents from the input list 00160 // InputCollection jetConstituents; 00161 // list< list<InputItem>::iterator>::const_iterator inp; 00162 // for(inp=cone.begin();inp!=cone.end();inp++) { 00163 // jetConstituents.push_back(**inp); 00164 // input.erase(*inp); 00165 // } 00166 // fOutput->push_back (ProtoJet (jetConstituents)); 00167 // 00168 // IMPORTANT NOTE: 00169 // while the stability of the stable cone is tested using the Et 00170 // scheme recombination, the final jet uses E-scheme 00171 // recombination. 00172 // 00173 // The technique used here is the same as what we already used for 00174 // SISCone except that we act on the 'cone' list. 00175 // We successively merge the particles that make up the cone jet 00176 // until we have all particles in it. We start off with the zeroth 00177 // particle. 00178 list< list<PseudoJet>::iterator>::const_iterator inp; 00179 inp = cone.begin(); 00180 int jet_k = (*inp)->cluster_hist_index(); 00181 // gps tmp 00182 //float px=(*inp)->px(), py=(*inp)->py(), pz=(*inp)->pz(), E = (*inp)->E(); 00183 00184 // remove the particle from the list and jump to the next one 00185 input.erase(*inp); 00186 inp++; 00187 00188 // now loop over the remaining particles 00189 while (inp != cone.end()){ 00190 // take the last result of the merge 00191 int jet_i = jet_k; 00192 // and the next element of the jet 00193 int jet_j = (*inp)->cluster_hist_index(); 00194 // and merge them (with a fake dij) 00195 double dij = 0.0; 00196 00197 // create the new jet by hand so that we can adjust its user index 00198 // Note again the use of the E-scheme recombination here! 00199 PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j]; 00200 00201 // gps tmp to try out floating issues 00202 //px+=(*inp)->px(), py+=(*inp)->py(), pz+=(*inp)->pz(), E += (*inp)->E(); 00203 //PseudoJet newjet(px,py,pz,E); 00204 00205 clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, newjet, jet_k); 00206 00207 // remove the particle from the list and jump to the next one 00208 input.erase(*inp); 00209 inp++; 00210 } 00211 00212 // we have merged all the jet's particles into a single object, so now 00213 // "declare" it to be a beam (inclusive) jet. 00214 // [NB: put a sensible looking d_iB just to be nice...] 00215 double d_iB = clust_seq.jets()[jet_k].perp2(); 00216 clust_seq.plugin_record_iB_recombination(jet_k, d_iB); 00217 00218 00219 } //loop over seeds ended 00220 00221 } 00222 00223 // print a banner for reference to the 3rd-party code 00224 void CMSIterativeConePlugin::_print_banner(ostream *ostr) const{ 00225 if (! _first_time) return; 00226 _first_time=false; 00227 00228 // make sure the user has not set the banner stream to NULL 00229 if (!ostr) return; 00230 00231 (*ostr) << "#-------------------------------------------------------------------------" << endl; 00232 (*ostr) << "# You are running the CMS Iterative Cone plugin for FastJet " << endl; 00233 (*ostr) << "# Original code by the CMS collaboration adapted by the FastJet authors " << endl; 00234 (*ostr) << "# If you use this plugin, please cite " << endl; 00235 (*ostr) << "# G. L. Bayatian et al. [CMS Collaboration], " << endl; 00236 (*ostr) << "# CMS physics: Technical design report. " << endl; 00237 (*ostr) << "# in addition to the usual FastJet reference. " << endl; 00238 (*ostr) << "#-------------------------------------------------------------------------" << endl; 00239 00240 // make sure we really have the output done. 00241 ostr->flush(); 00242 } 00243 00244 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh