FastJet  3.2.2
JadePlugin.cc
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 
47 using namespace std;
48 
49 FASTJET_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.
73 class JadeBriefJet {
74 public:
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 
115 private:
116  double rt2E, nx, ny, nz;
117 };
118 
119 
120 //----------------------------------------------------------------------
121 string 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 
179 template<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
199  nn.remove_jet(i);
200  }
201  njets--;
202  }
203 
204 }
205 
206 //----------------------------------------------------------------------
207 void 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 
222 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
double norm(const VPoint p)
norm of a vector
Definition: Voronoi.hh:155
deals with clustering
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...
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...