00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032 #include "fastjet/Error.hh"
00033 #include "fastjet/PseudoJet.hh"
00034 #include "fastjet/ClusterSequence.hh"
00035 #include<iostream>
00036 #include<sstream>
00037 #include<cmath>
00038 #include <cstdlib>
00039 #include<cassert>
00040
00041 #ifndef DROP_CGAL // in case we do not have the code for CGAL
00042 #include "fastjet/internal/Dnn4piCylinder.hh"
00043 #include "fastjet/internal/Dnn3piCylinder.hh"
00044 #include "fastjet/internal/Dnn2piCylinder.hh"
00045 #endif // DROP_CGAL
00046
00047 FASTJET_BEGIN_NAMESPACE
00048
00049 using namespace std;
00050
00051
00052
00058 void ClusterSequence::_delaunay_cluster () {
00059
00060 int n = _jets.size();
00061
00062 vector<EtaPhi> points(n);
00063 for (int i = 0; i < n; i++) {
00064 points[i] = EtaPhi(_jets[i].rap(),_jets[i].phi_02pi());
00065 points[i].sanitize();
00066 }
00067
00068
00069 DynamicNearestNeighbours * DNN;
00070 #ifndef DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
00071 bool verbose = false;
00072 bool ignore_nearest_is_mirror = (_Rparam < twopi);
00073 if (_strategy == NlnN4pi) {
00074 DNN = new Dnn4piCylinder(points,verbose);
00075 } else if (_strategy == NlnN3pi) {
00076 DNN = new Dnn3piCylinder(points,ignore_nearest_is_mirror,verbose);
00077 } else if (_strategy == NlnN) {
00078 DNN = new Dnn2piCylinder(points,ignore_nearest_is_mirror,verbose);
00079 } else
00080 #else
00081 if (_strategy == NlnN4pi || _strategy == NlnN3pi || _strategy == NlnN) {
00082 ostringstream err;
00083 err << "ERROR: Requested strategy "<<strategy_string()<<" but it is not"<<endl;
00084 err << " supported because FastJet was compiled without CGAL"<<endl;
00085 throw Error(err.str());
00086
00087 }
00088 #endif // DROP_CGAL
00089 {
00090 ostringstream err;
00091 err << "ERROR: Unrecognized value for strategy: "<<_strategy<<endl;
00092 assert(false);
00093 throw Error(err.str());
00094 }
00095
00096
00097
00098 DistMap DijMap;
00099
00100
00101
00102 for (int ii = 0; ii < n; ii++) {
00103 _add_ktdistance_to_map(ii, DijMap, DNN);
00104 }
00105
00106
00107
00108
00109 for (int i=0;i<n;i++) {
00110
00111
00112 TwoVertices SmallestDijPair;
00113 int jet_i, jet_j;
00114 double SmallestDij;
00115 bool Valid2;
00116 bool recombine_with_beam;
00117 do {
00118 SmallestDij = DijMap.begin()->first;
00119 SmallestDijPair = DijMap.begin()->second;
00120 jet_i = SmallestDijPair.first;
00121 jet_j = SmallestDijPair.second;
00122
00123
00124
00125
00126
00127 DijMap.erase(DijMap.begin());
00128
00129
00130
00131
00132 recombine_with_beam = (jet_j == BeamJet);
00133 if (!recombine_with_beam) {Valid2 = DNN->Valid(jet_j);}
00134 else {Valid2 = true;}
00135
00136 } while ( !DNN->Valid(jet_i) || !Valid2);
00137
00138
00139
00140
00141
00142 if (! recombine_with_beam) {
00143 int nn;
00144 _do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn);
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160 EtaPhi newpoint(_jets[nn].rap(), _jets[nn].phi_02pi());
00161 newpoint.sanitize();
00162 points.push_back(newpoint);
00163 } else {
00164
00165 _do_iB_recombination_step(jet_i, SmallestDij);
00166
00167
00168 }
00169
00170
00171
00172 if (i == n-1) {break;}
00173
00174 vector<int> updated_neighbours;
00175 if (! recombine_with_beam) {
00176
00177 int point3;
00178 DNN->RemoveCombinedAddCombination(jet_i, jet_j,
00179 points[points.size()-1], point3,
00180 updated_neighbours);
00181
00182
00183
00184 if (static_cast<unsigned int> (point3) != points.size()-1) {
00185 throw Error("INTERNAL ERROR: point3 != points.size()-1");}
00186 } else {
00187
00188 DNN->RemovePoint(jet_i, updated_neighbours);
00189 }
00190
00191
00192 vector<int>::iterator it = updated_neighbours.begin();
00193 for (; it != updated_neighbours.end(); ++it) {
00194 int ii = *it;
00195 _add_ktdistance_to_map(ii, DijMap, DNN);
00196 }
00197
00198 }
00199
00200
00201 delete DNN;
00202 }
00203
00204
00205
00220 void ClusterSequence::_add_ktdistance_to_map(
00221 const int & ii,
00222 DistMap & DijMap,
00223 const DynamicNearestNeighbours * DNN) {
00224
00225 double yiB = jet_scale_for_algorithm(_jets[ii]);
00226 if (yiB == 0.0) {
00227
00228
00229 DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
00230 } else {
00231 double DeltaR2 = DNN->NearestNeighbourDistance(ii) * _invR2;
00232
00233
00234
00235
00236
00237
00238
00239
00240
00241 if (DeltaR2 > 1.0) {
00242 DijMap.insert(DijEntry(yiB, TwoVertices(ii,-1)));
00243 } else {
00244 double kt2i = jet_scale_for_algorithm(_jets[ii]);
00245 int jj = DNN->NearestNeighbourIndex(ii);
00246 if (kt2i <= jet_scale_for_algorithm(_jets[jj])) {
00247 double dij = DeltaR2 * kt2i;
00248 DijMap.insert(DijEntry(dij, TwoVertices(ii,jj)));
00249 }
00250 }
00251 }
00252 }
00253
00254
00255 FASTJET_END_NAMESPACE
00256