Fork me on GitHub

source: git/external/fastjet/plugins/CMSIterativeCone/CMSIterativeConePlugin.cc@ 4b9a2dc

ImprovedOutputFile Timing dual_readout llp
Last change on this file since 4b9a2dc was 35cdc46, checked in by Pavel Demin <demin@…>, 10 years ago

upgrade FastJet to version 3.1.0-beta.1, upgrade Nsubjettiness to version 2.1.0, add SoftKiller version 1.0.0

  • Property mode set to 100644
File size: 8.5 KB
Line 
1//FJSTARTHEADER
2// $Id: CMSIterativeConePlugin.cc 1504 2009-04-10 13:39:48Z salam $
3//
4// Copyright (c) 2007-2014, 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
61FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
62
63using namespace std;
64using namespace cms;
65
66//------------------------------------------------------
67// some tools
68//------------------------------------------------------
69template <class T>
70T 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
77template <class T>
78T 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//------------------------------------------------------
85bool CMSIterativeConePlugin::_first_time = true;
86
87string CMSIterativeConePlugin::description () const {
88 ostringstream desc;
89 desc << "CMSIterativeCone plugin with R = " << theConeRadius << " and seed threshold = " << theSeedThreshold;
90 return desc.str();
91}
92
93void 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
217void 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
237FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
Note: See TracBrowser for help on using the repository browser.