Fork me on GitHub

source: svn/trunk/external/fastjet/ClusterSequence_DumbN3.cc@ 1378

Last change on this file since 1378 was 859, checked in by Pavel Demin, 12 years ago

update fastjet to version 3.0.3

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id Revision Date
File size: 4.5 KB
Line 
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
37FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
38
39
40using 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)!
48void 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
132FASTJET_END_NAMESPACE
133
Note: See TracBrowser for help on using the repository browser.