FastJet  3.0.6
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ClusterSequence_DumbN3.cc
1 //STARTHEADER
2 // $Id: ClusterSequence_DumbN3.cc 2690 2011-11-14 14:57:54Z soyez $
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