FastJet 3.0.5
ClusterSequence_DumbN3.cc
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 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends