ClusterSequence_N2.cc

Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: ClusterSequence_N2.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 // The plain N^2 part of the ClusterSequence class -- separated out
00033 // from the rest of the class implementation so as to speed up
00034 // compilation of this particular part while it is under test.
00035 
00036 #include "fastjet/PseudoJet.hh"
00037 #include "fastjet/ClusterSequence.hh"
00038 #include<iostream>
00039 
00040 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00041 
00042 
00043 using namespace std;
00044 
00045 //----------------------------------------------------------------------
00048 void ClusterSequence::_simple_N2_cluster() {
00049   int n = _jets.size();
00050   BriefJet * briefjets = new BriefJet[n];
00051   BriefJet * jetA = briefjets, * jetB;
00052   
00053   // initialise the basic jet info 
00054   for (int i = 0; i< n; i++) {
00055     _bj_set_jetinfo(jetA, i);
00056     jetA++; // move on to next entry of briefjets
00057   }
00058   BriefJet * tail = jetA; // a semaphore for the end of briefjets
00059   BriefJet * head = briefjets; // a nicer way of naming start
00060 
00061   // now initialise the NN distances: jetA will run from 1..n-1; and
00062   // jetB from 0..jetA-1
00063   for (jetA = head + 1; jetA != tail; jetA++) {
00064     // set NN info for jetA based on jets running from head..jetA-1,
00065     // checking in the process whether jetA itself is an undiscovered
00066     // NN of one of those jets.
00067     _bj_set_NN_crosscheck(jetA, head, jetA);
00068   }
00069 
00070 
00071   // now create the diJ (where J is i's NN) table -- remember that 
00072   // we differ from standard normalisation here by a factor of R2
00073   double * diJ = new double[n];
00074   jetA = head;
00075   for (int i = 0; i < n; i++) {
00076     diJ[i] = _bj_diJ(jetA);
00077     jetA++; // have jetA follow i
00078   }
00079 
00080   // now run the recombination loop
00081   int history_location = n-1;
00082   while (tail != head) {
00083 
00084     // find the minimum of the diJ on this round
00085     double diJ_min = diJ[0];
00086     int diJ_min_jet = 0;
00087     for (int i = 1; i < n; i++) {
00088       if (diJ[i] < diJ_min) {diJ_min_jet = i; diJ_min  = diJ[i];}
00089     }
00090     
00091     // do the recombination between A and B
00092     history_location++;
00093     jetA = & briefjets[diJ_min_jet];
00094     jetB = jetA->NN;
00095     // put the normalisation back in
00096     diJ_min *= _invR2; 
00097     if (jetB != NULL) {
00098       // jet-jet recombination
00099       // If necessary relabel A & B to ensure jetB < jetA, that way if
00100       // the larger of them == newtail then that ends up being jetA and 
00101       // the new jet that is added as jetB is inserted in a position that
00102       // has a future!
00103       if (jetA < jetB) {swap(jetA,jetB);}
00104 
00105       int nn; // new jet index
00106       _do_ij_recombination_step(jetA->_jets_index, jetB->_jets_index, diJ_min, nn);
00107 
00108       //OBS // get the two history indices
00109       //OBS int hist_a = _jets[jetA->_jets_index].cluster_hist_index();
00110       //OBS int hist_b = _jets[jetB->_jets_index].cluster_hist_index();
00111       //OBS // create the recombined jet
00112       //OBS _jets.push_back(_jets[jetA->_jets_index] + _jets[jetB->_jets_index]);
00113       //OBS int nn = _jets.size() - 1;
00114       //OBS _jets[nn].set_cluster_hist_index(history_location);
00115       //OBS // update history
00116       //OBS _add_step_to_history(history_location, 
00117       //OBS                min(hist_a,hist_b),max(hist_a,hist_b),
00118       //OBS                        nn, diJ_min);
00119       // what was jetB will now become the new jet
00120       _bj_set_jetinfo(jetB, nn);
00121     } else {
00122       // jet-beam recombination
00123       _do_iB_recombination_step(jetA->_jets_index, diJ_min);
00124       //OBS // get the hist_index
00125       //OBS int hist_a = _jets[jetA->_jets_index].cluster_hist_index();
00126       //OBS _add_step_to_history(history_location,hist_a,BeamJet,Invalid,diJ_min); 
00127     }
00128 
00129     // now update our nearest neighbour info and diJ table
00130     // first reduce size of table
00131     tail--; n--;
00132     // Copy last jet contents and diJ info into position of jetA
00133     *jetA = *tail;
00134     diJ[jetA - head] = diJ[tail-head];
00135 
00136     // Initialise jetB's NN distance as well as updating it for 
00137     // other particles.
00138     // NB: by having different loops for jetB == or != NULL we could
00139     //     perhaps save a few percent (usually avoid one if inside loop),
00140     //     but will not do it for now because on laptop fluctuations are
00141     //     too large to reliably measure a few percent difference...
00142     for (BriefJet * jetI = head; jetI != tail; jetI++) {
00143       // see if jetI had jetA or jetB as a NN -- if so recalculate the NN
00144       if (jetI->NN == jetA || jetI->NN == jetB) {
00145         _bj_set_NN_nocross(jetI, head, tail);
00146         diJ[jetI-head] = _bj_diJ(jetI); // update diJ 
00147       } 
00148       // check whether new jetB is closer than jetI's current NN and
00149       // if need be update things
00150       if (jetB != NULL) {
00151         double dist = _bj_dist(jetI,jetB);
00152         if (dist < jetI->NN_dist) {
00153           if (jetI != jetB) {
00154             jetI->NN_dist = dist;
00155             jetI->NN = jetB;
00156             diJ[jetI-head] = _bj_diJ(jetI); // update diJ...
00157           }
00158         }
00159         if (dist < jetB->NN_dist) {
00160           if (jetI != jetB) {
00161             jetB->NN_dist = dist;
00162             jetB->NN      = jetI;}
00163         }
00164       }
00165       // if jetI's NN is the new tail then relabel it so that it becomes jetA
00166       if (jetI->NN == tail) {jetI->NN = jetA;}
00167     }
00168 
00169     
00170     if (jetB != NULL) {diJ[jetB-head] = _bj_diJ(jetB);}
00171 
00172 
00173   }
00174   
00175 
00176   // final cleaning up;
00177   delete[] diJ;
00178   delete[] briefjets;
00179 }
00180 
00181 FASTJET_END_NAMESPACE
00182 

Generated on Tue Dec 18 17:05:03 2007 for fastjet by  doxygen 1.5.2