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