Fork me on GitHub

source: svn/trunk/external/fastjet/plugins/CMSIterativeCone/CMSIterativeConePlugin.cc@ 1373

Last change on this file since 1373 was 859, checked in by Pavel Demin, 12 years ago

update fastjet to version 3.0.3

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 8.9 KB
Line 
1//STARTHEADER
2// $Id: CMSIterativeConePlugin.cc 859 2012-11-28 01:49:23Z pavel $
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
68FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
69
70using namespace std;
71using namespace cms;
72
73//------------------------------------------------------
74// some tools
75//------------------------------------------------------
76template <class T>
77T 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
84template <class T>
85T 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//------------------------------------------------------
92bool CMSIterativeConePlugin::_first_time = true;
93
94string CMSIterativeConePlugin::description () const {
95 ostringstream desc;
96 desc << "CMSIterativeCone plugin with R = " << theConeRadius << " and seed threshold = " << theSeedThreshold;
97 return desc.str();
98}
99
100void 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
224void 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
244FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
Note: See TracBrowser for help on using the repository browser.