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