FastJet  3.1.2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
ClusterSequence_Delaunay.cc
1 //FJSTARTHEADER
2 // $Id: ClusterSequence_Delaunay.cc 3475 2014-07-29 11:57:23Z salam $
3 //
4 // Copyright (c) 2005-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 
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 
49 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
50 
51 using 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.
60 void 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  auto_ptr<DynamicNearestNeighbours> DNN;
72 #ifndef DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
73  bool verbose = false;
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  // distance is immediately removed regardless of whether or not
129  // it is used.
130  // Some temporary testing code relating to problems with the gcc-3.2 compiler
131  //cout << "got here and size is "<< DijMap.size()<< " and it is "<<SmallestDij <<"\n";
132  //cout << jet_i << " "<< jet_j<<"\n";
133  DijMap.erase(DijMap.begin());
134  //cout << "got beyond here\n";
135 
136  // need to "prime" the validity of jet_j in such a way that
137  // if it corresponds to the beam then it is automatically valid.
138  recombine_with_beam = (jet_j == BeamJet);
139  if (!recombine_with_beam) {Valid2 = DNN->Valid(jet_j);}
140  else {Valid2 = true;}
141 
142  } while ( !DNN->Valid(jet_i) || !Valid2);
143 
144 
145  // The following part acts just on jet momenta and on the history.
146  // The action on the nearest-neighbour structures takes place
147  // later (only if at least 2 jets are around).
148  if (! recombine_with_beam) {
149  int nn; // will be index of new jet
150  _do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn);
151  //OBS // merge the two jets, add new jet, remove old ones
152  //OBS _jets.push_back(_jets[jet_i] + _jets[jet_j]);
153  //OBS
154  //OBS int nn = _jets.size()-1;
155  //OBS _jets[nn].set_cluster_hist_index(n+i);
156  //OBS
157  //OBS // get corresponding indices in history structure
158  //OBS int hist_i = _jets[jet_i].cluster_hist_index();
159  //OBS int hist_j = _jets[jet_j].cluster_hist_index();
160  //OBS
161  //OBS
162  //OBS _add_step_to_history(n+i,min(hist_i,hist_j), max(hist_i,hist_j),
163  //OBS _jets.size()-1, SmallestDij);
164 
165  // add new point to points vector
166  EtaPhi newpoint(_jets[nn].rap(), _jets[nn].phi_02pi());
167  newpoint.sanitize(); // make sure it is in correct range
168  points.push_back(newpoint);
169  } else {
170  // recombine the jet with the beam
171  _do_iB_recombination_step(jet_i, SmallestDij);
172  //OBS _add_step_to_history(n+i,_jets[jet_i].cluster_hist_index(),BeamJet,
173  //OBS Invalid, SmallestDij);
174  }
175 
176  // exit the loop because we do not want to look for nearest neighbours
177  // etc. of zero partons
178  if (i == n-1) {break;}
179 
180  vector<int> updated_neighbours;
181  if (! recombine_with_beam) {
182  // update DNN
183  int point3;
184  DNN->RemoveCombinedAddCombination(jet_i, jet_j,
185  points[points.size()-1], point3,
186  updated_neighbours);
187  // C++ beginners' comment: static_cast to unsigned int is necessary
188  // to do away with warnings about type mismatch between point3 (int)
189  // and points.size (unsigned int)
190  if (static_cast<unsigned int> (point3) != points.size()-1) {
191  throw Error("INTERNAL ERROR: point3 != points.size()-1");}
192  } else {
193  // update DNN
194  DNN->RemovePoint(jet_i, updated_neighbours);
195  }
196 
197  // update map
198  vector<int>::iterator it = updated_neighbours.begin();
199  for (; it != updated_neighbours.end(); ++it) {
200  int ii = *it;
201  _add_ktdistance_to_map(ii, DijMap, DNN.get());
202  }
203 
204  } // end clustering loop
205 
206 }
207 
208 
209 //----------------------------------------------------------------------
210 /// Add the current kt distance for particle ii to the map (DijMap)
211 /// using information from the DNN object. Work as follows:
212 ///
213 /// . if the kt is zero then it's nearest neighbour is taken to be the
214 /// the beam jet and the distance is zero.
215 ///
216 /// . if cylinder distance to nearest neighbour > _Rparam then it is
217 /// yiB that is smallest and this is added to map.
218 ///
219 /// . otherwise if the nearest neighbour jj has a larger kt then add
220 /// dij to the map.
221 ///
222 /// . otherwise do nothing
223 ///
224 void ClusterSequence::_add_ktdistance_to_map(
225  const int ii,
226  DistMap & DijMap,
227  const DynamicNearestNeighbours * DNN) {
228 
229  double yiB = jet_scale_for_algorithm(_jets[ii]);
230  if (yiB == 0.0) {
231  // in this case convention is that we do not worry about distances
232  // but directly state that nearest neighbour is beam
233  DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
234  } else {
235  double DeltaR2 = DNN->NearestNeighbourDistance(ii) * _invR2;
236  // Logic of following bit is: only add point to map if it has
237  // smaller kt2 than nearest neighbour j (if it has larger kt,
238  // then: either it is j's nearest neighbour and then we will
239  // include dij when we come to j; or it is not j's nearest
240  // neighbour and j will recombine with someone else).
241 
242  // If DeltaR2 > 1.0 then in any case it will recombine with beam rather
243  // than with any neighbours.
244  // (put general normalisation here at some point)
245  if (DeltaR2 > 1.0) {
246  DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
247  } else {
248  double kt2i = jet_scale_for_algorithm(_jets[ii]);
249  int jj = DNN->NearestNeighbourIndex(ii);
250  if (kt2i <= jet_scale_for_algorithm(_jets[jj])) {
251  double dij = DeltaR2 * kt2i;
252  DijMap.insert(DijEntry(dij, TwoVertices(ii,jj)));
253  }
254  }
255  }
256 }
257 
258 
259 FASTJET_END_NAMESPACE
260