30 #include "fastjet/Error.hh"
31 #include "fastjet/PseudoJet.hh"
32 #include "fastjet/ClusterSequence.hh"
33 #include "fastjet/internal/DynamicNearestNeighbours.hh"
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"
47 FASTJET_BEGIN_NAMESPACE
58 void ClusterSequence::_delaunay_cluster () {
62 vector<EtaPhi> points(n);
63 for (
int i = 0; i < n; i++) {
64 points[i] = EtaPhi(_jets[i].rap(),_jets[i].phi_02pi());
69 auto_ptr<DynamicNearestNeighbours> DNN;
70 #ifndef DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
72 bool ignore_nearest_is_mirror = (_Rparam < twopi);
73 if (_strategy == NlnN4pi) {
74 DNN.reset(
new Dnn4piCylinder(points,verbose));
75 }
else if (_strategy == NlnN3pi) {
76 DNN.reset(
new Dnn3piCylinder(points,ignore_nearest_is_mirror,verbose));
77 }
else if (_strategy == NlnN) {
78 DNN.reset(
new Dnn2piCylinder(points,ignore_nearest_is_mirror,verbose));
81 if (_strategy == NlnN4pi || _strategy == NlnN3pi || _strategy == NlnN) {
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());
91 err <<
"ERROR: Unrecognized value for strategy: "<<_strategy<<endl;
93 throw Error(err.str());
102 for (
int ii = 0; ii < n; ii++) {
103 _add_ktdistance_to_map(ii, DijMap, DNN.get());
109 for (
int i=0;i<n;i++) {
112 TwoVertices SmallestDijPair;
116 bool recombine_with_beam;
118 SmallestDij = DijMap.begin()->first;
119 SmallestDijPair = DijMap.begin()->second;
120 jet_i = SmallestDijPair.first;
121 jet_j = SmallestDijPair.second;
127 DijMap.erase(DijMap.begin());
132 recombine_with_beam = (jet_j == BeamJet);
133 if (!recombine_with_beam) {Valid2 = DNN->Valid(jet_j);}
134 else {Valid2 =
true;}
136 }
while ( !DNN->Valid(jet_i) || !Valid2);
142 if (! recombine_with_beam) {
144 _do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn);
160 EtaPhi newpoint(_jets[nn].rap(), _jets[nn].phi_02pi());
162 points.push_back(newpoint);
165 _do_iB_recombination_step(jet_i, SmallestDij);
172 if (i == n-1) {
break;}
174 vector<int> updated_neighbours;
175 if (! recombine_with_beam) {
178 DNN->RemoveCombinedAddCombination(jet_i, jet_j,
179 points[points.size()-1], point3,
184 if (static_cast<unsigned int> (point3) != points.size()-1) {
185 throw Error(
"INTERNAL ERROR: point3 != points.size()-1");}
188 DNN->RemovePoint(jet_i, updated_neighbours);
192 vector<int>::iterator it = updated_neighbours.begin();
193 for (; it != updated_neighbours.end(); ++it) {
195 _add_ktdistance_to_map(ii, DijMap, DNN.get());
218 void ClusterSequence::_add_ktdistance_to_map(
221 const DynamicNearestNeighbours * DNN) {
223 double yiB = jet_scale_for_algorithm(_jets[ii]);
227 DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
229 double DeltaR2 = DNN->NearestNeighbourDistance(ii) * _invR2;
240 DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
242 double kt2i = jet_scale_for_algorithm(_jets[ii]);
243 int jj = DNN->NearestNeighbourIndex(ii);
244 if (kt2i <= jet_scale_for_algorithm(_jets[jj])) {
245 double dij = DeltaR2 * kt2i;
246 DijMap.insert(DijEntry(dij, TwoVertices(ii,jj)));
253 FASTJET_END_NAMESPACE