FastJet 3.0.2
CMSIterativeConePlugin.cc
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
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends