FastJet  3.0.6
CMSIterativeConePlugin.cc
1 //STARTHEADER
2 // $Id: CMSIterativeConePlugin.cc 1504 2009-04-10 13:39:48Z salam $
3 //
4 // Copyright (c) 2007-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 // Copyright (c) ????-????, CMS [for the iterative-cone code itself]
6 //
7 //----------------------------------------------------------------------
8 // This file is part of FastJet. It contains code that has been
9 // obtained from the CMS collaboration, revision 1.14 of the
10 // CMSIterativeConeAlgorithm.cc file in CMSSW, see
11 // http://cmssw.cvs.cern.ch/cgi-bin/cmssw.cgi/CMSSW/RecoJets/JetAlgorithms/src/CMSIterativeConeAlgorithm.cc?hideattic=0&revision=1.14&view=markup
12 //
13 // Permission has been granted by the CMS collaboration to release it
14 // in FastJet under the terms of the GNU Public License(v2) (see the
15 // COPYING file in the main FastJet directory for details).
16 // Changes from the original file are listed below.
17 //
18 // FastJet is free software; you can redistribute it and/or modify
19 // it under the terms of the GNU General Public License as published by
20 // the Free Software Foundation; either version 2 of the License, or
21 // (at your option) any later version.
22 //
23 // The algorithms that underlie FastJet have required considerable
24 // development and are described in hep-ph/0512210. If you use
25 // FastJet as part of work towards a scientific publication, please
26 // include a citation to the FastJet paper.
27 //
28 // FastJet is distributed in the hope that it will be useful,
29 // but WITHOUT ANY WARRANTY; without even the implied warranty of
30 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
31 // GNU General Public License for more details.
32 //
33 // You should have received a copy of the GNU General Public License
34 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
35 //----------------------------------------------------------------------
36 //ENDHEADER
37 
38 // List of changes compared to the original CMS code (revision 1.14 of
39 // CMSIterativeConeAlgorithm.cc)
40 //
41 // 2009-05-10 Gavin Salam <salam@lpthe.jussieu.fr>
42 //
43 // * added radius and seed threshold information in the plugin
44 // description
45 //
46 // 2009-01-06 Gregory Soyez <soyez@fastjet.fr>
47 //
48 // * Encapsulated the CMS code into a plugin for FastJet
49 // * inserted the deltaPhi and deltaR2 codes from
50 // DataFormats/Math/interface/deltaPhi.h (rev 1.1)
51 // DataFormats/Math/interface/deltaR.h (rev 1.2)
52 // * Adapted the code to use PseusoJet rather than 'InputItem'
53 // and 'InputCollection'
54 // * use the FastJet clustering history structures instead of
55 // the ProtoJet one used by CMS.
56 
57 
58 // fastjet stuff
59 #include "fastjet/ClusterSequence.hh"
60 #include "fastjet/CMSIterativeConePlugin.hh"
61 
62 // other stuff
63 #include <vector>
64 #include <list>
65 #include <sstream>
66 #include "SortByEt.h"
67 
68 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
69 
70 using namespace std;
71 using namespace cms;
72 
73 //------------------------------------------------------
74 // some tools
75 //------------------------------------------------------
76 template <class T>
77 T deltaPhi (T phi1, T phi2) {
78  T result = phi1 - phi2;
79  while (result > M_PI) result -= 2*M_PI;
80  while (result <= -M_PI) result += 2*M_PI;
81  return result;
82 }
83 
84 template <class T>
85 T deltaR2 (T eta1, T phi1, T eta2, T phi2) {
86  T deta = eta1 - eta2;
87  T dphi = deltaPhi (phi1, phi2);
88  return deta*deta + dphi*dphi;
89 }
90 
91 //------------------------------------------------------
92 bool CMSIterativeConePlugin::_first_time = true;
93 
94 string CMSIterativeConePlugin::description () const {
95  ostringstream desc;
96  desc << "CMSIterativeCone plugin with R = " << theConeRadius << " and seed threshold = " << theSeedThreshold;
97  return desc.str();
98 }
99 
100 void CMSIterativeConePlugin::run_clustering(ClusterSequence & clust_seq) const {
101  // print a banner if we run this for the first time
102  _print_banner(clust_seq.fastjet_banner_stream());
103 
104  //make a list of input objects ordered by ET
105  //cout << "copying the list of particles" << endl;
106  list<PseudoJet> input;
107  for (unsigned int i=0 ; i<clust_seq.jets().size() ; i++) {
108  input.push_back(clust_seq.jets()[i]);
109  }
110  NumericSafeGreaterByEt<PseudoJet> compCandidate;
111  //cout << "sorting" << endl;
112  input.sort(compCandidate);
113 
114  //find jets
115  //cout << "launching the main loop" << endl;
116  while( !input.empty() && (input.front().Et() > theSeedThreshold )) {
117  //cone centre
118  double eta0=input.front().eta();
119  double phi0=input.front().phi();
120  //protojet properties
121  double eta=0;
122  double phi=0;
123  double et=0;
124  //list of towers in cone
125  list< list<PseudoJet>::iterator> cone;
126  for(int iteration=0;iteration<100;iteration++){
127  //cout << "iterating" << endl;
128  cone.clear();
129  eta=0;
130  phi=0;
131  et=0;
132  for(list<PseudoJet>::iterator inp=input.begin();
133  inp!=input.end();inp++){
134  const PseudoJet tower = *inp;
135  if( deltaR2(eta0,phi0,tower.eta(),tower.phi()) <
136  theConeRadius*theConeRadius) {
137  double tower_et = tower.Et();
138  cone.push_back(inp);
139  eta+= tower_et*tower.eta();
140  double dphi=tower.phi()-phi0;
141  if(dphi>M_PI) dphi-=2*M_PI;
142  else if(dphi<=-M_PI) dphi+=2*M_PI;
143  phi+=tower_et*dphi;
144  et +=tower_et;
145  }
146  }
147  eta=eta/et;
148  phi=phi0+phi/et;
149  if(phi>M_PI)phi-=2*M_PI;
150  else if(phi<=-M_PI)phi+=2*M_PI;
151 
152  if(fabs(eta-eta0)<.001 && fabs(phi-phi0)<.001) break;//stable cone found
153  eta0=eta;
154  phi0=phi;
155  }
156 
157  //cout << "make the jet final" << endl;
158 
159  //make a final jet and remove the jet constituents from the input list
160  // InputCollection jetConstituents;
161  // list< list<InputItem>::iterator>::const_iterator inp;
162  // for(inp=cone.begin();inp!=cone.end();inp++) {
163  // jetConstituents.push_back(**inp);
164  // input.erase(*inp);
165  // }
166  // fOutput->push_back (ProtoJet (jetConstituents));
167  //
168  // IMPORTANT NOTE:
169  // while the stability of the stable cone is tested using the Et
170  // scheme recombination, the final jet uses E-scheme
171  // recombination.
172  //
173  // The technique used here is the same as what we already used for
174  // SISCone except that we act on the 'cone' list.
175  // We successively merge the particles that make up the cone jet
176  // until we have all particles in it. We start off with the zeroth
177  // particle.
178  list< list<PseudoJet>::iterator>::const_iterator inp;
179  inp = cone.begin();
180  int jet_k = (*inp)->cluster_hist_index();
181  // gps tmp
182  //float px=(*inp)->px(), py=(*inp)->py(), pz=(*inp)->pz(), E = (*inp)->E();
183 
184  // remove the particle from the list and jump to the next one
185  input.erase(*inp);
186  inp++;
187 
188  // now loop over the remaining particles
189  while (inp != cone.end()){
190  // take the last result of the merge
191  int jet_i = jet_k;
192  // and the next element of the jet
193  int jet_j = (*inp)->cluster_hist_index();
194  // and merge them (with a fake dij)
195  double dij = 0.0;
196 
197  // create the new jet by hand so that we can adjust its user index
198  // Note again the use of the E-scheme recombination here!
199  PseudoJet newjet = clust_seq.jets()[jet_i] + clust_seq.jets()[jet_j];
200 
201  // gps tmp to try out floating issues
202  //px+=(*inp)->px(), py+=(*inp)->py(), pz+=(*inp)->pz(), E += (*inp)->E();
203  //PseudoJet newjet(px,py,pz,E);
204 
205  clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, newjet, jet_k);
206 
207  // remove the particle from the list and jump to the next one
208  input.erase(*inp);
209  inp++;
210  }
211 
212  // we have merged all the jet's particles into a single object, so now
213  // "declare" it to be a beam (inclusive) jet.
214  // [NB: put a sensible looking d_iB just to be nice...]
215  double d_iB = clust_seq.jets()[jet_k].perp2();
216  clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
217 
218 
219  } //loop over seeds ended
220 
221 }
222 
223 // print a banner for reference to the 3rd-party code
224 void CMSIterativeConePlugin::_print_banner(ostream *ostr) const{
225  if (! _first_time) return;
226  _first_time=false;
227 
228  // make sure the user has not set the banner stream to NULL
229  if (!ostr) return;
230 
231  (*ostr) << "#-------------------------------------------------------------------------" << endl;
232  (*ostr) << "# You are running the CMS Iterative Cone plugin for FastJet " << endl;
233  (*ostr) << "# Original code by the CMS collaboration adapted by the FastJet authors " << endl;
234  (*ostr) << "# If you use this plugin, please cite " << endl;
235  (*ostr) << "# G. L. Bayatian et al. [CMS Collaboration], " << endl;
236  (*ostr) << "# CMS physics: Technical design report. " << endl;
237  (*ostr) << "# in addition to the usual FastJet reference. " << endl;
238  (*ostr) << "#-------------------------------------------------------------------------" << endl;
239 
240  // make sure we really have the output done.
241  ostr->flush();
242 }
243 
244 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
fastjet::PseudoJet::Et
double Et() const
return the transverse energy
Definition: PseudoJet.hh:166
fastjet::ClusterSequence
Definition: ClusterSequence.hh:59
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:346
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:959
fastjet::PseudoJet
Definition: PseudoJet.hh:65
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:567
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:329
fastjet::PseudoJet::phi
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:106