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 | |
---|
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 |
---|