Fork me on GitHub

source: svn/trunk/Utilities/Fastjet/src/ClusterSequence_Delaunay.cc@ 38

Last change on this file since 38 was 11, checked in by severine ovyn, 16 years ago

Fastjet added; CDFCones directory has been changed

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