FastJet 3.0.2
|
00001 //STARTHEADER 00002 // $Id: ClusterSequence_DumbN3.cc 2690 2011-11-14 14:57:54Z soyez $ 00003 // 00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez 00005 // 00006 //---------------------------------------------------------------------- 00007 // This file is part of FastJet. 00008 // 00009 // FastJet is free software; you can redistribute it and/or modify 00010 // it under the terms of the GNU General Public License as published by 00011 // the Free Software Foundation; either version 2 of the License, or 00012 // (at your option) any later version. 00013 // 00014 // The algorithms that underlie FastJet have required considerable 00015 // development and are described in hep-ph/0512210. If you use 00016 // FastJet as part of work towards a scientific publication, please 00017 // include a citation to the FastJet paper. 00018 // 00019 // FastJet is distributed in the hope that it will be useful, 00020 // but WITHOUT ANY WARRANTY; without even the implied warranty of 00021 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00022 // GNU General Public License for more details. 00023 // 00024 // You should have received a copy of the GNU General Public License 00025 // along with FastJet. If not, see <http://www.gnu.org/licenses/>. 00026 //---------------------------------------------------------------------- 00027 //ENDHEADER 00028 00029 00030 #include "fastjet/PseudoJet.hh" 00031 #include "fastjet/ClusterSequence.hh" 00032 #include<iostream> 00033 #include<cmath> 00034 #include <cstdlib> 00035 #include<cassert> 00036 00037 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00038 00039 00040 using namespace std; 00041 00042 00043 //---------------------------------------------------------------------- 00044 /// Run the clustering in a very slow variant of the N^3 algorithm. 00045 /// 00046 /// The only thing this routine has going for it is that memory usage 00047 /// is O(N)! 00048 void ClusterSequence::_really_dumb_cluster () { 00049 00050 // the array that will be overwritten here will be one 00051 // of pointers to jets. 00052 vector<PseudoJet *> jetsp(_jets.size()); 00053 vector<int> indices(_jets.size()); 00054 00055 for (size_t i = 0; i<_jets.size(); i++) { 00056 jetsp[i] = & _jets[i]; 00057 indices[i] = i; 00058 } 00059 00060 for (int n = jetsp.size(); n > 0; n--) { 00061 int ii, jj; 00062 // find smallest beam distance [remember jet_scale_for_algorithm 00063 // will return kt^2 for the kt-algorithm and 1 for the Cambridge/Aachen] 00064 double ymin = jet_scale_for_algorithm(*(jetsp[0])); 00065 ii = 0; jj = -2; 00066 for (int i = 0; i < n; i++) { 00067 double yiB = jet_scale_for_algorithm(*(jetsp[i])); 00068 if (yiB < ymin) { 00069 ymin = yiB; ii = i; jj = -2;} 00070 } 00071 00072 // find smallest distance between pair of jetsp 00073 for (int i = 0; i < n-1; i++) { 00074 for (int j = i+1; j < n; j++) { 00075 //double y = jetsp[i]->kt_distance(*jetsp[j])*_invR2; 00076 double y = min(jet_scale_for_algorithm(*(jetsp[i])), 00077 jet_scale_for_algorithm(*(jetsp[j]))) 00078 * jetsp[i]->plain_distance(*jetsp[j])*_invR2; 00079 if (y < ymin) {ymin = y; ii = i; jj = j;} 00080 } 00081 } 00082 00083 // output recombination sequence 00084 // old "ktclus" way of labelling 00085 //cout <<n<< " "<< ii+1 << " with " << jj+1 << "; y = "<< ymin<<endl; 00086 //OBS // new delaunay way of labelling 00087 //OBS int jjindex_or_beam, iiindex; 00088 //OBS if (jj < 0) {jjindex_or_beam = BeamJet; iiindex = indices[ii];} 00089 //OBS else { 00090 //OBS jjindex_or_beam = max(indices[ii],indices[jj]); 00091 //OBS iiindex = min(indices[ii],indices[jj]); 00092 //OBS } 00093 00094 // now recombine 00095 int newn = 2*jetsp.size() - n; 00096 if (jj >= 0) { 00097 // combine pair 00098 int nn; // new jet index 00099 _do_ij_recombination_step(jetsp[ii]-&_jets[0], 00100 jetsp[jj]-&_jets[0], ymin, nn); 00101 00102 // sort out internal bookkeeping 00103 jetsp[ii] = &_jets[nn]; 00104 // have jj point to jet that was pointed at by n-1 00105 // (since original jj is no longer current, so put n-1 into jj) 00106 jetsp[jj] = jetsp[n-1]; 00107 indices[ii] = newn; 00108 indices[jj] = indices[n-1]; 00109 00110 //OBS_jets.push_back(*jetsp[ii] + *jetsp[jj]); 00111 //OBSjetsp[ii] = &_jets[_jets.size()-1]; 00112 //OBS// have jj point to jet that was pointed at by n-1 00113 //OBS// (since original jj is no longer current, so put n-1 into jj) 00114 //OBSjetsp[jj] = jetsp[n-1]; 00115 //OBS 00116 //OBSindices[ii] = newn; 00117 //OBSindices[jj] = indices[n-1]; 00118 //OBS_add_step_to_history(newn,iiindex, 00119 //OBS jjindex_or_beam,_jets.size()-1,ymin); 00120 } else { 00121 // combine ii with beam 00122 _do_iB_recombination_step(jetsp[ii]-&_jets[0], ymin); 00123 // put last jet (pointer) in place of ii (which has disappeared) 00124 jetsp[ii] = jetsp[n-1]; 00125 indices[ii] = indices[n-1]; 00126 //OBS_add_step_to_history(newn,iiindex,jjindex_or_beam,Invalid, ymin); 00127 } 00128 } 00129 00130 } 00131 00132 FASTJET_END_NAMESPACE 00133