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