00001 //STARTHEADER 00002 // $Id: ClusterSequence_DumbN3.cc 293 2006-08-17 19:38:38Z salam $ 00003 // 00004 // Copyright (c) 2005-2006, Matteo Cacciari and Gavin Salam 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, write to the Free Software 00026 // Foundation, Inc.: 00027 // 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA 00028 //---------------------------------------------------------------------- 00029 //ENDHEADER 00030 00031 00032 #include "fastjet/PseudoJet.hh" 00033 #include "fastjet/ClusterSequence.hh" 00034 #include<iostream> 00035 #include<cmath> 00036 #include <cstdlib> 00037 #include<cassert> 00038 00039 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh 00040 00041 00042 using namespace std; 00043 00044 00045 //---------------------------------------------------------------------- 00050 void ClusterSequence::_really_dumb_cluster () { 00051 00052 // the array that will be overwritten here will be one 00053 // of pointers to jets. 00054 vector<PseudoJet *> jetsp(_jets.size()); 00055 vector<int> indices(_jets.size()); 00056 00057 for (size_t i = 0; i<_jets.size(); i++) { 00058 jetsp[i] = & _jets[i]; 00059 indices[i] = i; 00060 } 00061 00062 for (int n = jetsp.size(); n > 0; n--) { 00063 int ii, jj; 00064 // find smallest beam distance [remember jet_scale_for_algorithm 00065 // will return kt^2 for the kt-algorithm and 1 for the Cambridge/Aachen] 00066 double ymin = jet_scale_for_algorithm(*(jetsp[0])); 00067 ii = 0; jj = -2; 00068 for (int i = 0; i < n; i++) { 00069 double yiB = jet_scale_for_algorithm(*(jetsp[i])); 00070 if (yiB < ymin) { 00071 ymin = yiB; ii = i; jj = -2;} 00072 } 00073 00074 // find smallest distance between pair of jetsp 00075 for (int i = 0; i < n-1; i++) { 00076 for (int j = i+1; j < n; j++) { 00077 //double y = jetsp[i]->kt_distance(*jetsp[j])*_invR2; 00078 double y = min(jet_scale_for_algorithm(*(jetsp[i])), 00079 jet_scale_for_algorithm(*(jetsp[j]))) 00080 * jetsp[i]->plain_distance(*jetsp[j])*_invR2; 00081 if (y < ymin) {ymin = y; ii = i; jj = j;} 00082 } 00083 } 00084 00085 // output recombination sequence 00086 // old "ktclus" way of labelling 00087 //cout <<n<< " "<< ii+1 << " with " << jj+1 << "; y = "<< ymin<<endl; 00088 // new delaunay way of labelling 00089 int jjindex_or_beam, iiindex; 00090 if (jj < 0) {jjindex_or_beam = BeamJet; iiindex = indices[ii];} 00091 else { 00092 jjindex_or_beam = max(indices[ii],indices[jj]); 00093 iiindex = min(indices[ii],indices[jj]); 00094 } 00095 00096 // now recombine 00097 int newn = 2*jetsp.size() - n; 00098 if (jj >= 0) { 00099 // combine pair 00100 int nn; // new jet index 00101 _do_ij_recombination_step(jetsp[ii]-&_jets[0], 00102 jetsp[jj]-&_jets[0], ymin, nn); 00103 00104 // sort out internal bookkeeping 00105 jetsp[ii] = &_jets[nn]; 00106 // have jj point to jet that was pointed at by n-1 00107 // (since original jj is no longer current, so put n-1 into jj) 00108 jetsp[jj] = jetsp[n-1]; 00109 indices[ii] = newn; 00110 indices[jj] = indices[n-1]; 00111 00112 //OBS_jets.push_back(*jetsp[ii] + *jetsp[jj]); 00113 //OBSjetsp[ii] = &_jets[_jets.size()-1]; 00114 //OBS// have jj point to jet that was pointed at by n-1 00115 //OBS// (since original jj is no longer current, so put n-1 into jj) 00116 //OBSjetsp[jj] = jetsp[n-1]; 00117 //OBS 00118 //OBSindices[ii] = newn; 00119 //OBSindices[jj] = indices[n-1]; 00120 //OBS_add_step_to_history(newn,iiindex, 00121 //OBS jjindex_or_beam,_jets.size()-1,ymin); 00122 } else { 00123 // combine ii with beam 00124 _do_iB_recombination_step(jetsp[ii]-&_jets[0], ymin); 00125 // put last jet (pointer) in place of ii (which has disappeared) 00126 jetsp[ii] = jetsp[n-1]; 00127 indices[ii] = indices[n-1]; 00128 //OBS_add_step_to_history(newn,iiindex,jjindex_or_beam,Invalid, ymin); 00129 } 00130 } 00131 00132 } 00133 00134 FASTJET_END_NAMESPACE 00135