1 | //STARTHEADER |
---|
2 | // $Id: ClusterSequence_DumbN3.cc 859 2012-11-28 01:49:23Z pavel $ |
---|
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/PseudoJet.hh" |
---|
31 | #include "fastjet/ClusterSequence.hh" |
---|
32 | #include<iostream> |
---|
33 | #include<cmath> |
---|
34 | #include <cstdlib> |
---|
35 | #include<cassert> |
---|
36 | |
---|
37 | FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh |
---|
38 | |
---|
39 | |
---|
40 | using namespace std; |
---|
41 | |
---|
42 | |
---|
43 | //---------------------------------------------------------------------- |
---|
44 | /// Run the clustering in a very slow variant of the N^3 algorithm. |
---|
45 | /// |
---|
46 | /// The only thing this routine has going for it is that memory usage |
---|
47 | /// is O(N)! |
---|
48 | void ClusterSequence::_really_dumb_cluster () { |
---|
49 | |
---|
50 | // the array that will be overwritten here will be one |
---|
51 | // of pointers to jets. |
---|
52 | vector<PseudoJet *> jetsp(_jets.size()); |
---|
53 | vector<int> indices(_jets.size()); |
---|
54 | |
---|
55 | for (size_t i = 0; i<_jets.size(); i++) { |
---|
56 | jetsp[i] = & _jets[i]; |
---|
57 | indices[i] = i; |
---|
58 | } |
---|
59 | |
---|
60 | for (int n = jetsp.size(); n > 0; n--) { |
---|
61 | int ii, jj; |
---|
62 | // find smallest beam distance [remember jet_scale_for_algorithm |
---|
63 | // will return kt^2 for the kt-algorithm and 1 for the Cambridge/Aachen] |
---|
64 | double ymin = jet_scale_for_algorithm(*(jetsp[0])); |
---|
65 | ii = 0; jj = -2; |
---|
66 | for (int i = 0; i < n; i++) { |
---|
67 | double yiB = jet_scale_for_algorithm(*(jetsp[i])); |
---|
68 | if (yiB < ymin) { |
---|
69 | ymin = yiB; ii = i; jj = -2;} |
---|
70 | } |
---|
71 | |
---|
72 | // find smallest distance between pair of jetsp |
---|
73 | for (int i = 0; i < n-1; i++) { |
---|
74 | for (int j = i+1; j < n; j++) { |
---|
75 | //double y = jetsp[i]->kt_distance(*jetsp[j])*_invR2; |
---|
76 | double y = min(jet_scale_for_algorithm(*(jetsp[i])), |
---|
77 | jet_scale_for_algorithm(*(jetsp[j]))) |
---|
78 | * jetsp[i]->plain_distance(*jetsp[j])*_invR2; |
---|
79 | if (y < ymin) {ymin = y; ii = i; jj = j;} |
---|
80 | } |
---|
81 | } |
---|
82 | |
---|
83 | // output recombination sequence |
---|
84 | // old "ktclus" way of labelling |
---|
85 | //cout <<n<< " "<< ii+1 << " with " << jj+1 << "; y = "<< ymin<<endl; |
---|
86 | //OBS // new delaunay way of labelling |
---|
87 | //OBS int jjindex_or_beam, iiindex; |
---|
88 | //OBS if (jj < 0) {jjindex_or_beam = BeamJet; iiindex = indices[ii];} |
---|
89 | //OBS else { |
---|
90 | //OBS jjindex_or_beam = max(indices[ii],indices[jj]); |
---|
91 | //OBS iiindex = min(indices[ii],indices[jj]); |
---|
92 | //OBS } |
---|
93 | |
---|
94 | // now recombine |
---|
95 | int newn = 2*jetsp.size() - n; |
---|
96 | if (jj >= 0) { |
---|
97 | // combine pair |
---|
98 | int nn; // new jet index |
---|
99 | _do_ij_recombination_step(jetsp[ii]-&_jets[0], |
---|
100 | jetsp[jj]-&_jets[0], ymin, nn); |
---|
101 | |
---|
102 | // sort out internal bookkeeping |
---|
103 | jetsp[ii] = &_jets[nn]; |
---|
104 | // have jj point to jet that was pointed at by n-1 |
---|
105 | // (since original jj is no longer current, so put n-1 into jj) |
---|
106 | jetsp[jj] = jetsp[n-1]; |
---|
107 | indices[ii] = newn; |
---|
108 | indices[jj] = indices[n-1]; |
---|
109 | |
---|
110 | //OBS_jets.push_back(*jetsp[ii] + *jetsp[jj]); |
---|
111 | //OBSjetsp[ii] = &_jets[_jets.size()-1]; |
---|
112 | //OBS// have jj point to jet that was pointed at by n-1 |
---|
113 | //OBS// (since original jj is no longer current, so put n-1 into jj) |
---|
114 | //OBSjetsp[jj] = jetsp[n-1]; |
---|
115 | //OBS |
---|
116 | //OBSindices[ii] = newn; |
---|
117 | //OBSindices[jj] = indices[n-1]; |
---|
118 | //OBS_add_step_to_history(newn,iiindex, |
---|
119 | //OBS jjindex_or_beam,_jets.size()-1,ymin); |
---|
120 | } else { |
---|
121 | // combine ii with beam |
---|
122 | _do_iB_recombination_step(jetsp[ii]-&_jets[0], ymin); |
---|
123 | // put last jet (pointer) in place of ii (which has disappeared) |
---|
124 | jetsp[ii] = jetsp[n-1]; |
---|
125 | indices[ii] = indices[n-1]; |
---|
126 | //OBS_add_step_to_history(newn,iiindex,jjindex_or_beam,Invalid, ymin); |
---|
127 | } |
---|
128 | } |
---|
129 | |
---|
130 | } |
---|
131 | |
---|
132 | FASTJET_END_NAMESPACE |
---|
133 | |
---|