FastJet 3.0.0
|
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. 00009 // 00010 // FastJet is free software; you can redistribute it and/or modify 00011 // it under the terms of the GNU General Public License as published by 00012 // the Free Software Foundation; either version 2 of the License, or 00013 // (at your option) any later version. 00014 // 00015 // The algorithms that underlie FastJet have required considerable 00016 // development and are described in hep-ph/0512210. If you use 00017 // FastJet as part of work towards a scientific publication, please 00018 // include a citation to the FastJet paper. 00019 // 00020 // FastJet is distributed in the hope that it will be useful, 00021 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00022 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00023 // GNU General Public License for more details. 00024 // 00025 // You should have received a copy of the GNU General Public License 00026 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00027 //---------------------------------------------------------------------- 00028 //ENDHEADER 00029 00030 // fastjet stuff 00031 #include "fastjet/ClusterSequence.hh" 00032 #include "fastjet/CMSIterativeConePlugin.hh" 00033 00034 // other stuff 00035 #include <vector> 00036 #include <list> 00037 #include <sstream> 00038 #include "SortByEt.h" 00039 00040 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00041 00042 using namespace std; 00043 using namespace cms; 00044 00045 //------------------------------------------------------ 00046 // some tools 00047 //------------------------------------------------------ 00048 template <class T> 00049 T deltaPhi (T phi1, T phi2) { 00050 T result = phi1 - phi2; 00051 while (result > M_PI) result -= 2*M_PI; 00052 while (result <= -M_PI) result += 2*M_PI; 00053 return result; 00054 } 00055 00056 template <class T> 00057 T deltaR2 (T eta1, T phi1, T eta2, T phi2) { 00058 T deta = eta1 - eta2; 00059 T dphi = deltaPhi (phi1, phi2); 00060 return deta*deta + dphi*dphi; 00061 } 00062 00063 //------------------------------------------------------ 00064 00065 00066 string CMSIterativeConePlugin::description () const { 00067 ostringstream desc; 00068 desc << "CMSIterativeCone plugin with R = " << theConeRadius << " and seed threshold = " << theSeedThreshold; 00069 return desc.str(); 00070 } 00071 00072 void CMSIterativeConePlugin::run_clustering(ClusterSequence & clust_seq) const { 00073 00074 // This code is adapted from CMSIterativeConeAlgorithms.cc from the 00075 // CMSSW software. 00076 // The adaptation is just meant to use 00077 // - the FastJet 4-vectors instead of the CMS ones 00078 // - the FastJet clustering history structures instead of the 00079 // ProtoJet one used by CMS. 00080 00081 //make a list of input objects ordered by ET 00082 //cout << "copying the list of particles" << endl; 00083 list<PseudoJet> input; 00084 for (unsigned int i=0 ; i<clust_seq.jets().size() ; i++) { 00085 input.push_back(clust_seq.jets()[i]); 00086 } 00087 NumericSafeGreaterByEt<PseudoJet> compCandidate; 00088 //cout << "sorting" << endl; 00089 input.sort(compCandidate); 00090 00091 //find jets 00092 //cout << "launching the main loop" << endl; 00093 while( !input.empty() && (input.front().Et() > theSeedThreshold )) { 00094 //cone centre 00095 double eta0=input.front().eta(); 00096 double phi0=input.front().phi(); 00097 //protojet properties 00098 double eta=0; 00099 double phi=0; 00100 double et=0; 00101 //list of towers in cone 00102 list< list<PseudoJet>::iterator> cone; 00103 for(int iteration=0;iteration<100;iteration++){ 00104 //cout << "iterating" << endl; 00105 cone.clear(); 00106 eta=0; 00107 phi=0; 00108 et=0; 00109 for(list<PseudoJet>::iterator inp=input.begin(); 00110 inp!=input.end();inp++){ 00111 const PseudoJet tower = *inp; 00112 if( deltaR2(eta0,phi0,tower.eta(),tower.phi()) < 00113 theConeRadius*theConeRadius) { 00114 double tower_et = tower.Et(); 00115 cone.push_back(inp); 00116 eta+= tower_et*tower.eta(); 00117 double dphi=tower.phi()-phi0; 00118 if(dphi>M_PI) dphi-=2*M_PI; 00119 else if(dphi<=-M_PI) dphi+=2*M_PI; 00120 phi+=tower_et*dphi; 00121 et +=tower_et; 00122 } 00123 } 00124 eta=eta/et; 00125 phi=phi0+phi/et; 00126 if(phi>M_PI)phi-=2*M_PI; 00127 else if(phi<=-M_PI)phi+=2*M_PI; 00128 00129 if(fabs(eta-eta0)<.001 && fabs(phi-phi0)<.001) break;//stable cone found 00130 eta0=eta; 00131 phi0=phi; 00132 } 00133 00134 //cout << "make the jet final" << endl; 00135 00136 //make a final jet and remove the jet constituents from the input list 00137 // InputCollection jetConstituents; 00138 // list< list<InputItem>::iterator>::const_iterator inp; 00139 // for(inp=cone.begin();inp!=cone.end();inp++) { 00140 // jetConstituents.push_back(**inp); 00141 // input.erase(*inp); 00142 // } 00143 // fOutput->push_back (ProtoJet (jetConstituents)); 00144 // 00145 // IMPORTANT NOTE: 00146 // while the stability of the stable cone is tested using the Et 00147 // scheme recombination, the final jet uses E-scheme 00148 // recombination. 00149 // 00150 // The technique used here is the same as what we already used for 00151 // SISCone except that we act on the 'cone' list. 00152 // We successively merge the particles that make up the cone jet 00153 // until we have all particles in it. We start off with the zeroth 00154 // particle. 00155 list< list<PseudoJet>::iterator>::const_iterator inp; 00156 inp = cone.begin(); 00157 int jet_k = (*inp)->cluster_hist_index(); 00158 // gps tmp 00159 //float px=(*inp)->px(), py=(*inp)->py(), pz=(*inp)->pz(), E = (*inp)->E(); 00160 00161 // remove the particle from the list and jump to the next one 00162 input.erase(*inp); 00163 inp++; 00164 00165 // now loop over the remaining particles 00166 while (inp != cone.end()){ 00167 // take the last result of the merge 00168 int jet_i = jet_k; 00169 // and the next element of the jet 00170 int jet_j = (*inp)->cluster_hist_index(); 00171 // and merge them (with a fake dij) 00172 double dij = 0.0; 00173 00174 // create the new jet by hand so that we can adjust its user index 00175 // Note again the use of the E-scheme recombination here! 00176 PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j]; 00177 00178 // gps tmp to try out floating issues 00179 //px+=(*inp)->px(), py+=(*inp)->py(), pz+=(*inp)->pz(), E += (*inp)->E(); 00180 //PseudoJet newjet(px,py,pz,E); 00181 00182 clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, newjet, jet_k); 00183 00184 // remove the particle from the list and jump to the next one 00185 input.erase(*inp); 00186 inp++; 00187 } 00188 00189 // we have merged all the jet's particles into a single object, so now 00190 // "declare" it to be a beam (inclusive) jet. 00191 // [NB: put a sensible looking d_iB just to be nice...] 00192 double d_iB = clust_seq.jets()[jet_k].perp2(); 00193 clust_seq.plugin_record_iB_recombination(jet_k, d_iB); 00194 00195 00196 } //loop over seeds ended 00197 00198 00199 } 00200 00201 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh