Fork me on GitHub

source: git/external/fastjet/plugins/Jade/JadePlugin.cc@ ca9e119

ImprovedOutputFile Timing dual_readout llp
Last change on this file since ca9e119 was 1d208a2, checked in by Pavel Demin <pavel.demin@…>, 8 years ago

update FastJet library to 3.2.1 and Nsubjettiness library to 2.2.4

  • Property mode set to 100644
File size: 6.5 KB
Line 
1//FJSTARTHEADER
2// $Id: JadePlugin.cc 4063 2016-03-04 10:31:40Z 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// fastjet stuff
32#include "fastjet/ClusterSequence.hh"
33#include "fastjet/JadePlugin.hh"
34#include <iostream>
35//#include "fastjet/internal/ClusterSequence_N2.icc"
36#include "fastjet/NNH.hh"
37#include "fastjet/NNFJN2Plain.hh"
38
39// other stuff
40#include <vector>
41#include <sstream>
42#include <limits>
43
44
45
46
47using namespace std;
48
49FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
50
51
52//----------------------------------------------------------------------
53/// class to help run a JADE algorithm
54///
55/// This class works both with NNH and NNFJN2Plain clustering
56/// helpers. They both use the same init(...) call, but for the
57/// clustering:
58///
59/// - NNH uses distance(...) and beam_distance()
60/// - NNFJPlainN2 uses geometrical_distance(...), momentum_factor()
61/// and geometrical_beam_distance()
62///
63/// For NNFJPlainN2 the 2 E_i E_j (1-cos theta_{ij}) factor
64/// gets broken up into
65///
66/// sqrt(2)*min(E_i,E_j) * [sqrt(2)*max(E_i,E_j) (1 - cos \theta_{ij})]
67///
68/// The second factor is what we call the "geometrical_distance" even
69/// though it isn't actually purely geometrical. But the fact that it
70/// gets multiplied by min(E_i,E_j) to get the full distance is
71/// sufficient for the validity of the FJ lemma, allowing for the use
72/// of NNFJN2Plain.
73class JadeBriefJet {
74public:
75 void init(const PseudoJet & jet) {
76 double norm = 1.0/sqrt(jet.modp2());
77 nx = jet.px() * norm;
78 ny = jet.py() * norm;
79 nz = jet.pz() * norm;
80 rt2E = sqrt(2.0)*jet.E();
81 }
82
83 double distance(const JadeBriefJet * jet) const {
84 double dij = 1 - nx*jet->nx
85 - ny*jet->ny
86 - nz*jet->nz;
87 dij *= rt2E*jet->rt2E;
88 return dij;
89 }
90
91 double geometrical_distance(const JadeBriefJet * jet) const {
92 double dij = 1 - nx*jet->nx
93 - ny*jet->ny
94 - nz*jet->nz;
95 dij *= max(rt2E,jet->rt2E);
96 return dij;
97 }
98
99 double momentum_factor() const {
100 return rt2E;
101 }
102
103 double beam_distance() const {
104 return numeric_limits<double>::max();
105 }
106
107 double geometrical_beam_distance() const {
108 // get a number that is almost the same as max(), just a little
109 // smaller so as to ensure that when we divide it by rt2E and then
110 // multiply it again, we won't get an overflow
111 const double almost_max = numeric_limits<double>::max() * (1 - 1e-13);
112 return almost_max / rt2E;
113 }
114
115private:
116 double rt2E, nx, ny, nz;
117};
118
119
120//----------------------------------------------------------------------
121string JadePlugin::description () const {
122 ostringstream desc;
123 desc << "e+e- JADE algorithm plugin";
124 switch(_strategy) {
125 case strategy_NNH:
126 desc << ", using NNH strategy"; break;
127 case strategy_NNFJN2Plain:
128 desc << ", using NNFJN2Plain strategy"; break;
129 default:
130 throw Error("Unrecognized strategy in JadePlugin");
131 }
132
133 return desc.str();
134}
135
136// //----------------------------------------------------------------------
137// void JadePlugin::run_clustering(ClusterSequence & cs) const {
138// int njets = cs.jets().size();
139//
140// //SharedPtr<NNBase<> > nn;
141// NNBase<> * nn;
142// switch(_strategy) {
143// case strategy_NNH:
144// //nn.reset(new NNH<JadeBriefJet>(cs.jets()));
145// nn = new NNH<JadeBriefJet>(cs.jets());
146// break;
147// case strategy_NNFJN2Plain:
148// //nn.reset(new NNFJN2Plain<JadeBriefJet>(cs.jets()));
149// nn = new NNFJN2Plain<JadeBriefJet>(cs.jets());
150// break;
151// default:
152// throw Error("Unrecognized strategy in JadePlugin");
153// }
154// //NNH<JadeBriefJet> nnh(cs.jets());
155// //NNFJN2Plain<JadeBriefJet> nnh(cs.jets());
156//
157// // if testing against Hoeth's implementation, need to rescale the
158// // dij by Q^2.
159// //double Q2 = cs.Q2();
160//
161// while (njets > 0) {
162// int i, j, k;
163// double dij = nn->dij_min(i, j);
164//
165// if (j >= 0) {
166// cs.plugin_record_ij_recombination(i, j, dij, k);
167// nn->merge_jets(i, j, cs.jets()[k], k);
168// } else {
169// double diB = cs.jets()[i].E()*cs.jets()[i].E(); // get new diB
170// cs.plugin_record_iB_recombination(i, diB);
171// nn->remove_jet(i);
172// }
173// njets--;
174// }
175// delete nn;
176// }
177
178
179template<class N> void JadePlugin::_actual_run_clustering(ClusterSequence & cs) const {
180
181 int njets = cs.jets().size();
182
183 N nn(cs.jets());
184
185 // if testing against Hoeth's implementation, need to rescale the
186 // dij by Q^2.
187 //double Q2 = cs.Q2();
188
189 while (njets > 0) {
190 int i, j, k;
191 double dij = nn.dij_min(i, j);
192
193 if (j >= 0) {
194 cs.plugin_record_ij_recombination(i, j, dij, k);
195 nn.merge_jets(i, j, cs.jets()[k], k);
196 } else {
197 double diB = cs.jets()[i].E()*cs.jets()[i].E(); // get new diB
198 cs.plugin_record_iB_recombination(i, diB);
199 nn.remove_jet(i);
200 }
201 njets--;
202 }
203
204}
205
206//----------------------------------------------------------------------
207void JadePlugin::run_clustering(ClusterSequence & cs) const {
208
209 switch(_strategy) {
210 case strategy_NNH:
211 _actual_run_clustering<NNH<JadeBriefJet> >(cs);
212 break;
213 case strategy_NNFJN2Plain:
214 _actual_run_clustering<NNFJN2Plain<JadeBriefJet> >(cs);
215 break;
216 default:
217 throw Error("Unrecognized strategy in JadePlugin");
218 }
219}
220
221
222FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
Note: See TracBrowser for help on using the repository browser.