FastJet  3.3.0
ClusterSequence_DumbN3.cc
1 //FJSTARTHEADER
2 // $Id: ClusterSequence_DumbN3.cc 3433 2014-07-23 08:17:03Z salam $
3 //
4 // Copyright (c) 2005-2014, 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. They are described in the original FastJet paper,
16 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17 // FastJet as part of work towards a scientific publication, please
18 // quote the version you use and include a citation to the manual and
19 // optionally also to hep-ph/0512210.
20 //
21 // FastJet is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28 //----------------------------------------------------------------------
29 //FJENDHEADER
30 
31 
32 #include "fastjet/PseudoJet.hh"
33 #include "fastjet/ClusterSequence.hh"
34 #include<iostream>
35 #include<cmath>
36 #include <cstdlib>
37 #include<cassert>
38 
39 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
40 
41 
42 using namespace std;
43 
44 
45 //----------------------------------------------------------------------
46 /// Run the clustering in a very slow variant of the N^3 algorithm.
47 ///
48 /// The only thing this routine has going for it is that memory usage
49 /// is O(N)!
50 void ClusterSequence::_really_dumb_cluster () {
51 
52  // the array that will be overwritten here will be one
53  // of pointers to jets.
54  vector<PseudoJet *> jetsp(_jets.size());
55  vector<int> indices(_jets.size());
56 
57  for (size_t i = 0; i<_jets.size(); i++) {
58  jetsp[i] = & _jets[i];
59  indices[i] = i;
60  }
61 
62  for (int n = jetsp.size(); n > 0; n--) {
63  int ii, jj;
64  // find smallest beam distance [remember jet_scale_for_algorithm
65  // will return kt^2 for the kt-algorithm and 1 for the Cambridge/Aachen]
66  double ymin = jet_scale_for_algorithm(*(jetsp[0]));
67  ii = 0; jj = -2;
68  for (int i = 0; i < n; i++) {
69  double yiB = jet_scale_for_algorithm(*(jetsp[i]));
70  if (yiB < ymin) {
71  ymin = yiB; ii = i; jj = -2;}
72  }
73 
74  // find smallest distance between pair of jetsp
75  for (int i = 0; i < n-1; i++) {
76  for (int j = i+1; j < n; j++) {
77  //double y = jetsp[i]->kt_distance(*jetsp[j])*_invR2;
78  double y = min(jet_scale_for_algorithm(*(jetsp[i])),
79  jet_scale_for_algorithm(*(jetsp[j])))
80  * jetsp[i]->plain_distance(*jetsp[j])*_invR2;
81  if (y < ymin) {ymin = y; ii = i; jj = j;}
82  }
83  }
84 
85  // output recombination sequence
86  // old "ktclus" way of labelling
87  //cout <<n<< " "<< ii+1 << " with " << jj+1 << "; y = "<< ymin<<endl;
88  //OBS // new delaunay way of labelling
89  //OBS int jjindex_or_beam, iiindex;
90  //OBS if (jj < 0) {jjindex_or_beam = BeamJet; iiindex = indices[ii];}
91  //OBS else {
92  //OBS jjindex_or_beam = max(indices[ii],indices[jj]);
93  //OBS iiindex = min(indices[ii],indices[jj]);
94  //OBS }
95 
96  // now recombine
97  int newn = 2*jetsp.size() - n;
98  if (jj >= 0) {
99  // combine pair
100  int nn; // new jet index
101  _do_ij_recombination_step(jetsp[ii]-&_jets[0],
102  jetsp[jj]-&_jets[0], ymin, nn);
103 
104  // sort out internal bookkeeping
105  jetsp[ii] = &_jets[nn];
106  // have jj point to jet that was pointed at by n-1
107  // (since original jj is no longer current, so put n-1 into jj)
108  jetsp[jj] = jetsp[n-1];
109  indices[ii] = newn;
110  indices[jj] = indices[n-1];
111 
112  //OBS_jets.push_back(*jetsp[ii] + *jetsp[jj]);
113  //OBSjetsp[ii] = &_jets[_jets.size()-1];
114  //OBS// have jj point to jet that was pointed at by n-1
115  //OBS// (since original jj is no longer current, so put n-1 into jj)
116  //OBSjetsp[jj] = jetsp[n-1];
117  //OBS
118  //OBSindices[ii] = newn;
119  //OBSindices[jj] = indices[n-1];
120  //OBS_add_step_to_history(newn,iiindex,
121  //OBS jjindex_or_beam,_jets.size()-1,ymin);
122  } else {
123  // combine ii with beam
124  _do_iB_recombination_step(jetsp[ii]-&_jets[0], ymin);
125  // put last jet (pointer) in place of ii (which has disappeared)
126  jetsp[ii] = jetsp[n-1];
127  indices[ii] = indices[n-1];
128  //OBS_add_step_to_history(newn,iiindex,jjindex_or_beam,Invalid, ymin);
129  }
130  }
131 
132 }
133 
134 FASTJET_END_NAMESPACE
135