FastJet 3.4.1
ClusterSequence_Delaunay.cc
1//FJSTARTHEADER
2// $Id$
3//
4// Copyright (c) 2005-2023, 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
32#include "fastjet/Error.hh"
33#include "fastjet/PseudoJet.hh"
34#include "fastjet/ClusterSequence.hh"
35#include "fastjet/internal/DynamicNearestNeighbours.hh"
36#include<iostream>
37#include<sstream>
38#include<cmath>
39#include <cstdlib>
40#include<cassert>
41#include<memory>
42//
43#ifndef DROP_CGAL // in case we do not have the code for CGAL
44#include "fastjet/internal/Dnn4piCylinder.hh"
45#include "fastjet/internal/Dnn3piCylinder.hh"
46#include "fastjet/internal/Dnn2piCylinder.hh"
47#endif // DROP_CGAL
48
49FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
50
51using namespace std;
52
53
54//----------------------------------------------------------------------
55/// Run the clustering using a Hierarchical Delaunay triangulation and
56/// STL maps to achieve O(N*ln N) behaviour.
57///
58/// There may be internally asserted assumptions about absence of
59/// points with coincident eta-phi coordinates.
60void ClusterSequence::_delaunay_cluster () {
61
62 int n = _jets.size();
63
64 vector<EtaPhi> points(n); // recall EtaPhi is just a typedef'd pair<double>
65 for (int i = 0; i < n; i++) {
66 points[i] = EtaPhi(_jets[i].rap(),_jets[i].phi_02pi());
67 points[i].sanitize(); // make sure things are in the right range
68 }
69
70 // initialise our DNN structure with the set of points
71 const bool verbose = false;
72 SharedPtr<DynamicNearestNeighbours> DNN;
73#ifndef DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
74 bool ignore_nearest_is_mirror = (_Rparam < twopi);
75 if (_strategy == NlnN4pi) {
76 DNN.reset(new Dnn4piCylinder(points,verbose));
77 } else if (_strategy == NlnN3pi) {
78 DNN.reset(new Dnn3piCylinder(points,ignore_nearest_is_mirror,verbose));
79 } else if (_strategy == NlnN) {
80 DNN.reset(new Dnn2piCylinder(points,ignore_nearest_is_mirror,verbose));
81 } else
82#else
83 if (_strategy == NlnN4pi || _strategy == NlnN3pi || _strategy == NlnN) {
84 ostringstream err;
85 err << "ERROR: Requested strategy "<<strategy_string()<<" but it is not"<<endl;
86 err << " supported because FastJet was compiled without CGAL"<<endl;
87 throw Error(err.str());
88 //assert(false);
89 } else
90#endif // DROP_CGAL
91 {
92 //ostringstream err;
93 //err << "ERROR: Unrecognized value for strategy: "<<_strategy<<endl;
94 //throw Error(err.str());
95 //-----------------------------------------------------------------
96 // The code should never reach this point, because the checks above
97 // should always handle all _strategy values for which
98 // _delaunay_cluster() is called
99 assert(false);
100 }
101
102 // We will find nearest neighbour for each vertex, and include
103 // distance in map (NB DistMap is a typedef given in the .h file)
104 DistMap DijMap;
105
106 // fill the map with the minimal (as far as we know) subset of Dij
107 // distances (i.e. nearest neighbour ones).
108 for (int ii = 0; ii < n; ii++) {
109 _add_ktdistance_to_map(ii, DijMap, DNN.get());
110 }
111
112 // run the clustering (go up to i=n-1, but then will stop half-way down,
113 // when we reach that point -- it will be the final beam jet and there
114 // will be no nearest neighbours to find).
115 for (int i=0;i<n;i++) {
116 // find nearest vertices
117 // NB: skip cases where the point is not there anymore!
118 TwoVertices SmallestDijPair;
119 int jet_i, jet_j;
120 double SmallestDij;
121 bool Valid2;
122 bool recombine_with_beam;
123 do {
124 SmallestDij = DijMap.begin()->first;
125 SmallestDijPair = DijMap.begin()->second;
126 jet_i = SmallestDijPair.first;
127 jet_j = SmallestDijPair.second;
128 if (verbose) cout << "CS_Delaunay found recombination candidate: " << jet_i << " " << jet_j << " " << SmallestDij << endl; // GPS debugging
129 // distance is immediately removed regardless of whether or not
130 // it is used.
131 // Some temporary testing code relating to problems with the gcc-3.2 compiler
132 //cout << "got here and size is "<< DijMap.size()<< " and it is "<<SmallestDij <<"\n";
133 //cout << jet_i << " "<< jet_j<<"\n";
134 DijMap.erase(DijMap.begin());
135 //cout << "got beyond here\n";
136
137 // need to "prime" the validity of jet_j in such a way that
138 // if it corresponds to the beam then it is automatically valid.
139 recombine_with_beam = (jet_j == BeamJet);
140 if (!recombine_with_beam) {Valid2 = DNN->Valid(jet_j);}
141 else {Valid2 = true;}
142 if (verbose) cout << "CS_Delaunay validities i & j: " << DNN->Valid(jet_i) << " " << Valid2 << endl;
143 } while ( !DNN->Valid(jet_i) || !Valid2);
144
145
146 // The following part acts just on jet momenta and on the history.
147 // The action on the nearest-neighbour structures takes place
148 // later (only if at least 2 jets are around).
149 if (! recombine_with_beam) {
150 int nn; // will be index of new jet
151 if (verbose) cout << "CS_Delaunay call _do_ij_recomb: " << jet_i << " " << jet_j << " " << SmallestDij << endl; // GPS debug
152 _do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn);
153 //OBS // merge the two jets, add new jet, remove old ones
154 //OBS _jets.push_back(_jets[jet_i] + _jets[jet_j]);
155 //OBS
156 //OBS int nn = _jets.size()-1;
157 //OBS _jets[nn].set_cluster_hist_index(n+i);
158 //OBS
159 //OBS // get corresponding indices in history structure
160 //OBS int hist_i = _jets[jet_i].cluster_hist_index();
161 //OBS int hist_j = _jets[jet_j].cluster_hist_index();
162 //OBS
163 //OBS
164 //OBS _add_step_to_history(n+i,min(hist_i,hist_j), max(hist_i,hist_j),
165 //OBS _jets.size()-1, SmallestDij);
166
167 // add new point to points vector
168 EtaPhi newpoint(_jets[nn].rap(), _jets[nn].phi_02pi());
169 newpoint.sanitize(); // make sure it is in correct range
170 points.push_back(newpoint);
171 } else {
172 // recombine the jet with the beam
173 if (verbose) cout << "CS_Delaunay call _do_iB_recomb: " << jet_i << " " << SmallestDij << endl; // GPS debug
174 _do_iB_recombination_step(jet_i, SmallestDij);
175 //OBS _add_step_to_history(n+i,_jets[jet_i].cluster_hist_index(),BeamJet,
176 //OBS Invalid, SmallestDij);
177 }
178
179 // exit the loop because we do not want to look for nearest neighbours
180 // etc. of zero partons
181 if (i == n-1) {break;}
182
183 vector<int> updated_neighbours;
184 if (! recombine_with_beam) {
185 // update DNN
186 int point3;
187 DNN->RemoveCombinedAddCombination(jet_i, jet_j,
188 points[points.size()-1], point3,
189 updated_neighbours);
190 // C++ beginners' comment: static_cast to unsigned int is necessary
191 // to do away with warnings about type mismatch between point3 (int)
192 // and points.size (unsigned int)
193 if (static_cast<unsigned int> (point3) != points.size()-1) {
194 throw Error("INTERNAL ERROR: point3 != points.size()-1");}
195 } else {
196 // update DNN
197 DNN->RemovePoint(jet_i, updated_neighbours);
198 }
199
200 // update map
201 vector<int>::iterator it = updated_neighbours.begin();
202 for (; it != updated_neighbours.end(); ++it) {
203 int ii = *it;
204 _add_ktdistance_to_map(ii, DijMap, DNN.get());
205 }
206
207 } // end clustering loop
208
209}
210
211
212//----------------------------------------------------------------------
213/// Add the current kt distance for particle ii to the map (DijMap)
214/// using information from the DNN object. Work as follows:
215///
216/// . if the kt is zero then it's nearest neighbour is taken to be the
217/// the beam jet and the distance is zero.
218///
219/// . if cylinder distance to nearest neighbour > _Rparam then it is
220/// yiB that is smallest and this is added to map.
221///
222/// . otherwise if the nearest neighbour jj has a larger kt then add
223/// dij to the map.
224///
225/// . otherwise do nothing
226///
227void ClusterSequence::_add_ktdistance_to_map(
228 const int ii,
229 DistMap & DijMap,
230 const DynamicNearestNeighbours * DNN) {
231
232 double yiB = jet_scale_for_algorithm(_jets[ii]);
233 if (yiB == 0.0) {
234 // in this case convention is that we do not worry about distances
235 // but directly state that nearest neighbour is beam
236 DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
237 } else {
238 double DeltaR2 = DNN->NearestNeighbourDistance(ii) * _invR2;
239 // Logic of following bit is: only add point to map if it has
240 // smaller kt2 than nearest neighbour j (if it has larger kt,
241 // then: either it is j's nearest neighbour and then we will
242 // include dij when we come to j; or it is not j's nearest
243 // neighbour and j will recombine with someone else).
244
245 // If DeltaR2 > 1.0 then in any case it will recombine with beam rather
246 // than with any neighbours.
247 // (put general normalisation here at some point)
248 if (DeltaR2 > 1.0) {
249 DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
250 } else {
251 double kt2i = jet_scale_for_algorithm(_jets[ii]);
252 int jj = DNN->NearestNeighbourIndex(ii);
253 if (kt2i <= jet_scale_for_algorithm(_jets[jj])) {
254 double dij = DeltaR2 * kt2i;
255 DijMap.insert(DijEntry(dij, TwoVertices(ii,jj)));
256 }
257 }
258 }
259}
260
261
262FASTJET_END_NAMESPACE
263