00001 #ifndef __NNH_HH__
00002 #define __NNH_HH__
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
00033
00034 #include<fastjet/ClusterSequence.hh>
00035
00036
00037 FASTJET_BEGIN_NAMESPACE
00038
00040 class _NoInfo {};
00041
00043 template<class I> class NNHInfo {
00044 public:
00045 NNHInfo() : _info(NULL) {}
00046 NNHInfo(I * info) : _info(info) {}
00047 template<class NNBJ> void init_jet(NNBJ * briefjet, const fastjet::PseudoJet & jet, int index) { briefjet->init(jet, index, _info);}
00048 private:
00049 I * _info;
00050 };
00051
00053 template<> class NNHInfo<_NoInfo> {
00054 public:
00055 NNHInfo() {}
00056 NNHInfo(_NoInfo * info) {}
00057 template<class NNBJ> void init_jet(NNBJ * briefjet, const fastjet::PseudoJet & jet, int index) { briefjet->init(jet, index);}
00058 };
00059
00060
00061
00101
00102 template<class BJ, class I = _NoInfo> class NNH : public NNHInfo<I> {
00103 public:
00104
00107 NNH(const std::vector<PseudoJet> & jets) {start(jets);}
00108 NNH(const std::vector<PseudoJet> & jets, I * info) : NNHInfo<I>(info) {start(jets);}
00109
00110 void start(const std::vector<PseudoJet> & jets);
00111
00114 double dij_min(int & iA, int & iB);
00115
00117 void remove_jet(int iA);
00118
00121 void merge_jets(int iA, int iB, const PseudoJet & jet, int jet_index);
00122
00124 ~NNH() {
00125 delete[] briefjets;
00126 }
00127
00128 private:
00129 class NNBJ;
00130
00134 void set_NN_crosscheck(NNBJ * jet, NNBJ * begin, NNBJ * end);
00135
00138 void set_NN_nocross (NNBJ * jet, NNBJ * begin, NNBJ * end);
00139
00141 NNBJ * briefjets;
00142
00144 NNBJ * head, * tail;
00145
00147 int n;
00148
00150 std::vector<NNBJ *> where_is;
00151
00154 class NNBJ : public BJ {
00155 public:
00156 void init(const PseudoJet & jet, int index) {
00157 BJ::init(jet);
00158 other_init(index);
00159 }
00160 void init(const PseudoJet & jet, int index, I * info) {
00161 BJ::init(jet, info);
00162 other_init(index);
00163 }
00164 void other_init(int index) {
00165 _index = index;
00166 NN_dist = BJ::beam_distance();
00167 NN = NULL;
00168 }
00169 int index() const {return _index;}
00170
00171 double NN_dist;
00172 NNBJ * NN;
00173
00174 private:
00175 int _index;
00176 };
00177
00178 };
00179
00180
00181
00182
00183 template<class BJ, class I> void NNH<BJ,I>::start(const std::vector<PseudoJet> & jets) {
00184 n = jets.size();
00185 briefjets = new NNBJ[n];
00186 where_is.resize(2*n);
00187
00188 NNBJ * jetA = briefjets;
00189
00190
00191 for (int i = 0; i< n; i++) {
00192
00193 init_jet(jetA, jets[i], i);
00194 where_is[i] = jetA;
00195 jetA++;
00196 }
00197 tail = jetA;
00198 head = briefjets;
00199
00200
00201
00202 for (jetA = head + 1; jetA != tail; jetA++) {
00203
00204
00205
00206 set_NN_crosscheck(jetA, head, jetA);
00207 }
00208
00209 }
00210
00211
00212
00213 template<class BJ, class I> double NNH<BJ,I>::dij_min(int & iA, int & iB) {
00214
00215 double diJ_min = briefjets[0].NN_dist;
00216 int diJ_min_jet = 0;
00217 for (int i = 1; i < n; i++) {
00218 if (briefjets[i].NN_dist < diJ_min) {
00219 diJ_min_jet = i;
00220 diJ_min = briefjets[i].NN_dist;
00221 }
00222 }
00223
00224
00225 NNBJ * jetA = & briefjets[diJ_min_jet];
00226
00227 iA = jetA->index();
00228 iB = jetA->NN ? jetA->NN->index() : -1;
00229 return diJ_min;
00230 }
00231
00232
00233
00234
00235 template<class BJ, class I> void NNH<BJ,I>::remove_jet(int iA) {
00236 NNBJ * jetA = where_is[iA];
00237
00238
00239 tail--; n--;
00240
00241 *jetA = *tail;
00242
00243 where_is[jetA->index()] = jetA;
00244
00245 for (NNBJ * jetI = head; jetI != tail; jetI++) {
00246
00247 if (jetI->NN == jetA) set_NN_nocross(jetI, head, tail);
00248
00249
00250 if (jetI->NN == tail) {jetI->NN = jetA;}
00251 }
00252 }
00253
00254
00255
00256 template<class BJ, class I> void NNH<BJ,I>::merge_jets(int iA, int iB,
00257 const PseudoJet & jet, int index) {
00258
00259 NNBJ * jetA = where_is[iA];
00260 NNBJ * jetB = where_is[iB];
00261
00262
00263
00264
00265
00266 if (jetA < jetB) swap(jetA,jetB);
00267
00268
00269
00270 init_jet(jetB, jet, index);
00271
00272 if (index >= int(where_is.size())) where_is.resize(2*index);
00273 where_is[jetB->index()] = jetB;
00274
00275
00276
00277 tail--; n--;
00278
00279 *jetA = *tail;
00280
00281 where_is[jetA->index()] = jetA;
00282
00283 for (NNBJ * jetI = head; jetI != tail; jetI++) {
00284
00285 if (jetI->NN == jetA || jetI->NN == jetB) {
00286 set_NN_nocross(jetI, head, tail);
00287 }
00288
00289
00290
00291 double dist = jetI->distance(jetB);
00292 if (dist < jetI->NN_dist) {
00293 if (jetI != jetB) {
00294 jetI->NN_dist = dist;
00295 jetI->NN = jetB;
00296 }
00297 }
00298 if (dist < jetB->NN_dist) {
00299 if (jetI != jetB) {
00300 jetB->NN_dist = dist;
00301 jetB->NN = jetI;
00302 }
00303 }
00304
00305
00306 if (jetI->NN == tail) {jetI->NN = jetA;}
00307 }
00308 }
00309
00310
00311
00312
00313 template <class BJ, class I> void NNH<BJ,I>::set_NN_crosscheck(NNBJ * jet,
00314 NNBJ * begin, NNBJ * end) {
00315 double NN_dist = jet->beam_distance();
00316 NNBJ * NN = NULL;
00317 for (NNBJ * jetB = begin; jetB != end; jetB++) {
00318 double dist = jet->distance(jetB);
00319 if (dist < NN_dist) {
00320 NN_dist = dist;
00321 NN = jetB;
00322 }
00323 if (dist < jetB->NN_dist) {
00324 jetB->NN_dist = dist;
00325 jetB->NN = jet;
00326 }
00327 }
00328 jet->NN = NN;
00329 jet->NN_dist = NN_dist;
00330 }
00331
00332
00333
00334
00335
00336 template <class BJ, class I> void NNH<BJ,I>::set_NN_nocross(
00337 NNBJ * jet, NNBJ * begin, NNBJ * end) {
00338 double NN_dist = jet->beam_distance();
00339 NNBJ * NN = NULL;
00340 if (head < jet) {
00341 for (NNBJ * jetB = head; jetB != jet; jetB++) {
00342 double dist = jet->distance(jetB);
00343 if (dist < NN_dist) {
00344 NN_dist = dist;
00345 NN = jetB;
00346 }
00347 }
00348 }
00349 if (tail > jet) {
00350 for (NNBJ * jetB = jet+1; jetB != tail; jetB++) {
00351 double dist = jet->distance (jetB);
00352 if (dist < NN_dist) {
00353 NN_dist = dist;
00354 NN = jetB;
00355 }
00356 }
00357 }
00358 jet->NN = NN;
00359 jet->NN_dist = NN_dist;
00360 }
00361
00362
00363
00364
00365 FASTJET_END_NAMESPACE
00366
00367
00368 #endif // __NNH_HH__