Loading...
Searching...
No Matches
CMSIterativeConePlugin.cc
//FJSTARTHEADER
// $Id: CMSIterativeConePlugin.cc 1504 2009-04-10 13:39:48Z salam $
//
// Copyright (c) 2007-2021, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
// FastJet is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// The algorithms that underlie FastJet have required considerable
// development. They are described in the original FastJet paper,
// hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
// FastJet as part of work towards a scientific publication, please
// quote the version you use and include a citation to the manual and
// optionally also to hep-ph/0512210.
//
// FastJet is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//FJENDHEADER
// List of changes compared to the original CMS code (revision 1.14 of
// CMSIterativeConeAlgorithm.cc)
//
// 2009-05-10 Gavin Salam <salam@lpthe.jussieu.fr>
//
// * added radius and seed threshold information in the plugin
// description
//
// 2009-01-06 Gregory Soyez <soyez@fastjet.fr>
//
// * Encapsulated the CMS code into a plugin for FastJet
// * inserted the deltaPhi and deltaR2 codes from
// DataFormats/Math/interface/deltaPhi.h (rev 1.1)
// DataFormats/Math/interface/deltaR.h (rev 1.2)
// * Adapted the code to use PseusoJet rather than 'InputItem'
// and 'InputCollection'
// * use the FastJet clustering history structures instead of
// the ProtoJet one used by CMS.
// fastjet stuff
#include "fastjet/ClusterSequence.hh"
#include "fastjet/CMSIterativeConePlugin.hh"
// other stuff
#include <vector>
#include <list>
#include <sstream>
#include "SortByEt.h"
FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
using namespace std;
using namespace cms;
//------------------------------------------------------
// some tools
//------------------------------------------------------
template <class T>
T deltaPhi (T phi1, T phi2) {
T result = phi1 - phi2;
while (result > M_PI) result -= 2*M_PI;
while (result <= -M_PI) result += 2*M_PI;
return result;
}
template <class T>
T deltaR2 (T eta1, T phi1, T eta2, T phi2) {
T deta = eta1 - eta2;
T dphi = deltaPhi (phi1, phi2);
return deta*deta + dphi*dphi;
}
//------------------------------------------------------
thread_safety_helpers::FirstTimeTrue CMSIterativeConePlugin::_first_time;
string CMSIterativeConePlugin::description () const {
ostringstream desc;
desc << "CMSIterativeCone plugin with R = " << theConeRadius << " and seed threshold = " << theSeedThreshold;
return desc.str();
}
void CMSIterativeConePlugin::run_clustering(ClusterSequence & clust_seq) const {
// print a banner if we run this for the first time
_print_banner(clust_seq.fastjet_banner_stream());
//make a list of input objects ordered by ET
//cout << "copying the list of particles" << endl;
list<PseudoJet> input;
for (unsigned int i=0 ; i<clust_seq.jets().size() ; i++) {
input.push_back(clust_seq.jets()[i]);
}
NumericSafeGreaterByEt<PseudoJet> compCandidate;
//cout << "sorting" << endl;
input.sort(compCandidate);
//find jets
//cout << "launching the main loop" << endl;
while( !input.empty() && (input.front().Et() > theSeedThreshold )) {
//cone centre
double eta0=input.front().eta();
double phi0=input.front().phi();
//protojet properties
double eta=0;
double phi=0;
double et=0;
//list of towers in cone
list< list<PseudoJet>::iterator> cone;
for(int iteration=0;iteration<100;iteration++){
//cout << "iterating" << endl;
cone.clear();
eta=0;
phi=0;
et=0;
for(list<PseudoJet>::iterator inp=input.begin();
inp!=input.end();inp++){
const PseudoJet tower = *inp;
if( deltaR2(eta0,phi0,tower.eta(),tower.phi()) <
theConeRadius*theConeRadius) {
double tower_et = tower.Et();
cone.push_back(inp);
eta+= tower_et*tower.eta();
double dphi=tower.phi()-phi0;
if(dphi>M_PI) dphi-=2*M_PI;
else if(dphi<=-M_PI) dphi+=2*M_PI;
phi+=tower_et*dphi;
et +=tower_et;
}
}
eta=eta/et;
phi=phi0+phi/et;
if(phi>M_PI)phi-=2*M_PI;
else if(phi<=-M_PI)phi+=2*M_PI;
if(fabs(eta-eta0)<.001 && fabs(phi-phi0)<.001) break;//stable cone found
eta0=eta;
phi0=phi;
}
//cout << "make the jet final" << endl;
//make a final jet and remove the jet constituents from the input list
// InputCollection jetConstituents;
// list< list<InputItem>::iterator>::const_iterator inp;
// for(inp=cone.begin();inp!=cone.end();inp++) {
// jetConstituents.push_back(**inp);
// input.erase(*inp);
// }
// fOutput->push_back (ProtoJet (jetConstituents));
//
// IMPORTANT NOTE:
// while the stability of the stable cone is tested using the Et
// scheme recombination, the final jet uses E-scheme
// recombination.
//
// The technique used here is the same as what we already used for
// SISCone except that we act on the 'cone' list.
// We successively merge the particles that make up the cone jet
// until we have all particles in it. We start off with the zeroth
// particle.
list< list<PseudoJet>::iterator>::const_iterator inp;
inp = cone.begin();
int jet_k = (*inp)->cluster_hist_index();
// gps tmp
//float px=(*inp)->px(), py=(*inp)->py(), pz=(*inp)->pz(), E = (*inp)->E();
// remove the particle from the list and jump to the next one
input.erase(*inp);
inp++;
// now loop over the remaining particles
while (inp != cone.end()){
// take the last result of the merge
int jet_i = jet_k;
// and the next element of the jet
int jet_j = (*inp)->cluster_hist_index();
// and merge them (with a fake dij)
double dij = 0.0;
// create the new jet by hand so that we can adjust its user index
// Note again the use of the E-scheme recombination here!
PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j];
// gps tmp to try out floating issues
//px+=(*inp)->px(), py+=(*inp)->py(), pz+=(*inp)->pz(), E += (*inp)->E();
//PseudoJet newjet(px,py,pz,E);
clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, newjet, jet_k);
// remove the particle from the list and jump to the next one
input.erase(*inp);
inp++;
}
// we have merged all the jet's particles into a single object, so now
// "declare" it to be a beam (inclusive) jet.
// [NB: put a sensible looking d_iB just to be nice...]
double d_iB = clust_seq.jets()[jet_k].perp2();
clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
} //loop over seeds ended
}
// print a banner for reference to the 3rd-party code
void CMSIterativeConePlugin::_print_banner(ostream *ostr) const{
if (! _first_time()) return;
// make sure the user has not set the banner stream to NULL
if (!ostr) return;
(*ostr) << "#-------------------------------------------------------------------------" << endl;
(*ostr) << "# You are running the CMS Iterative Cone plugin for FastJet " << endl;
(*ostr) << "# Original code by the CMS collaboration adapted by the FastJet authors " << endl;
(*ostr) << "# If you use this plugin, please cite " << endl;
(*ostr) << "# G. L. Bayatian et al. [CMS Collaboration], " << endl;
(*ostr) << "# CMS physics: Technical design report. " << endl;
(*ostr) << "# in addition to the usual FastJet reference. " << endl;
(*ostr) << "#-------------------------------------------------------------------------" << endl;
// make sure we really have the output done.
ostr->flush();
}
FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
fastjet::CMSIterativeConePlugin::description
virtual std::string description() const
return a textual description of the jet-definition implemented in this plugin
Definition CMSIterativeConePlugin.cc:87
fastjet::CMSIterativeConePlugin::run_clustering
virtual void run_clustering(ClusterSequence &) const
given a ClusterSequence that has been filled up with initial particles, the following function should...
Definition CMSIterativeConePlugin.cc:93
fastjet::ClusterSequence
deals with clustering
Definition ClusterSequence.hh:63
fastjet::ClusterSequence::plugin_record_iB_recombination
void plugin_record_iB_recombination(int jet_i, double diB)
record the fact that there has been a recombination between jets()[jet_i] and the beam,...
Definition ClusterSequence.hh:343
fastjet::ClusterSequence::jets
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
Definition ClusterSequence.hh:1003
fastjet::ClusterSequence::fastjet_banner_stream
static std::ostream * fastjet_banner_stream()
returns a pointer to the stream to be used to print banners (cout by default).
Definition ClusterSequence.hh:579
fastjet::ClusterSequence::plugin_record_ij_recombination
void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, int &newjet_k)
record the fact that there has been a recombination between jets()[jet_i] and jets()[jet_k],...
Definition ClusterSequence.hh:326
fastjet::PseudoJet
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
Definition PseudoJet.hh:68
fastjet::PseudoJet::Et
double Et() const
return the transverse energy
Definition PseudoJet.hh:183
fastjet::PseudoJet::phi
double phi() const
returns phi (in the range 0..2pi)
Definition PseudoJet.hh:123