FastJet 3.0.2
ClusterSequence_Delaunay.cc
00001 //STARTHEADER
00002 // $Id: ClusterSequence_Delaunay.cc 2577 2011-09-13 15:11:38Z salam $
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/Error.hh"
00031 #include "fastjet/PseudoJet.hh"
00032 #include "fastjet/ClusterSequence.hh"
00033 #include<iostream>
00034 #include<sstream>
00035 #include<cmath>
00036 #include <cstdlib>
00037 #include<cassert>
00038 #include<memory>
00039 //
00040 #ifndef DROP_CGAL // in case we do not have the code for CGAL
00041 #include "fastjet/internal/Dnn4piCylinder.hh"
00042 #include "fastjet/internal/Dnn3piCylinder.hh"
00043 #include "fastjet/internal/Dnn2piCylinder.hh"
00044 #endif //  DROP_CGAL 
00045 
00046 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00047 
00048 using namespace std;
00049 
00050 
00051 //----------------------------------------------------------------------
00052 /// Run the clustering using a Hierarchical Delaunay triangulation and
00053 /// STL maps to achieve O(N*ln N) behaviour.
00054 ///
00055 /// There may be internally asserted assumptions about absence of
00056 /// points with coincident eta-phi coordinates.
00057 void ClusterSequence::_delaunay_cluster () {
00058 
00059   int n = _jets.size();
00060 
00061   vector<EtaPhi> points(n); // recall EtaPhi is just a typedef'd pair<double>
00062   for (int i = 0; i < n; i++) {
00063     points[i] = EtaPhi(_jets[i].rap(),_jets[i].phi_02pi());
00064     points[i].sanitize(); // make sure things are in the right range
00065   }
00066 
00067   // initialise our DNN structure with the set of points
00068   auto_ptr<DynamicNearestNeighbours> DNN;
00069 #ifndef DROP_CGAL // strategy = NlnN* are not supported if we drop CGAL...
00070   bool verbose = false;
00071   bool ignore_nearest_is_mirror = (_Rparam < twopi);
00072   if (_strategy == NlnN4pi) {
00073     DNN.reset(new Dnn4piCylinder(points,verbose));
00074   } else if (_strategy == NlnN3pi) {
00075     DNN.reset(new Dnn3piCylinder(points,ignore_nearest_is_mirror,verbose));
00076   } else if (_strategy == NlnN) {
00077     DNN.reset(new Dnn2piCylinder(points,ignore_nearest_is_mirror,verbose));
00078   } else 
00079 #else
00080   if (_strategy == NlnN4pi || _strategy == NlnN3pi || _strategy == NlnN) {
00081     ostringstream err;
00082     err << "ERROR: Requested strategy "<<strategy_string()<<" but it is not"<<endl;
00083     err << "       supported because FastJet was compiled without CGAL"<<endl;
00084     throw Error(err.str());
00085     //assert(false);
00086   }
00087 #endif // DROP_CGAL
00088   {
00089     ostringstream err;
00090     err << "ERROR: Unrecognized value for strategy: "<<_strategy<<endl;
00091     assert(false);
00092     throw Error(err.str());
00093   }
00094 
00095   // We will find nearest neighbour for each vertex, and include
00096   // distance in map (NB DistMap is a typedef given in the .h file)
00097   DistMap DijMap;
00098 
00099   // fill the map with the minimal (as far as we know) subset of Dij
00100   // distances (i.e. nearest neighbour ones).
00101   for (int ii = 0; ii < n; ii++) {
00102     _add_ktdistance_to_map(ii, DijMap, DNN.get());
00103   }
00104 
00105   // run the clustering (go up to i=n-1, but then will stop half-way down,
00106   // when we reach that point -- it will be the final beam jet and there
00107   // will be no nearest neighbours to find).
00108   for (int i=0;i<n;i++) {
00109     // find nearest vertices
00110     // NB: skip cases where the point is not there anymore!
00111     TwoVertices SmallestDijPair;
00112     int jet_i, jet_j;
00113     double SmallestDij;
00114     bool Valid2;
00115     bool recombine_with_beam;
00116     do { 
00117       SmallestDij = DijMap.begin()->first;
00118       SmallestDijPair = DijMap.begin()->second;
00119       jet_i = SmallestDijPair.first;
00120       jet_j = SmallestDijPair.second;
00121       // distance is immediately removed regardless of whether or not
00122       // it is used.
00123       // Some temporary testing code relating to problems with the gcc-3.2 compiler
00124       //cout << "got here and size is "<< DijMap.size()<< " and it is "<<SmallestDij <<"\n";
00125       //cout <<  jet_i << " "<< jet_j<<"\n";
00126       DijMap.erase(DijMap.begin());
00127       //cout << "got beyond here\n";
00128 
00129       // need to "prime" the validity of jet_j in such a way that 
00130       // if it corresponds to the beam then it is automatically valid.
00131       recombine_with_beam = (jet_j == BeamJet);
00132       if (!recombine_with_beam) {Valid2 = DNN->Valid(jet_j);} 
00133       else {Valid2 = true;}
00134 
00135     } while ( !DNN->Valid(jet_i) || !Valid2);
00136 
00137 
00138     // The following part acts just on jet momenta and on the history.
00139     // The action on the nearest-neighbour structures takes place
00140     // later (only if at least 2 jets are around).
00141     if (! recombine_with_beam) {
00142       int nn; // will be index of new jet
00143       _do_ij_recombination_step(jet_i, jet_j, SmallestDij, nn);
00144       //OBS // merge the two jets, add new jet, remove old ones
00145       //OBS _jets.push_back(_jets[jet_i] + _jets[jet_j]);
00146       //OBS 
00147       //OBS int nn = _jets.size()-1;
00148       //OBS _jets[nn].set_cluster_hist_index(n+i);
00149       //OBS 
00150       //OBS // get corresponding indices in history structure
00151       //OBS int hist_i = _jets[jet_i].cluster_hist_index();
00152       //OBS int hist_j = _jets[jet_j].cluster_hist_index();
00153       //OBS 
00154       //OBS 
00155       //OBS _add_step_to_history(n+i,min(hist_i,hist_j), max(hist_i,hist_j),
00156       //OBS                   _jets.size()-1, SmallestDij);
00157 
00158       // add new point to points vector
00159       EtaPhi newpoint(_jets[nn].rap(), _jets[nn].phi_02pi());
00160       newpoint.sanitize(); // make sure it is in correct range
00161       points.push_back(newpoint);
00162     } else {
00163       // recombine the jet with the beam
00164       _do_iB_recombination_step(jet_i, SmallestDij);
00165       //OBS _add_step_to_history(n+i,_jets[jet_i].cluster_hist_index(),BeamJet,
00166       //OBS                        Invalid, SmallestDij);
00167     }
00168 
00169     // exit the loop because we do not want to look for nearest neighbours
00170     // etc. of zero partons
00171     if (i == n-1) {break;}
00172 
00173     vector<int> updated_neighbours;
00174     if (! recombine_with_beam) {
00175       // update DNN
00176       int point3;
00177       DNN->RemoveCombinedAddCombination(jet_i, jet_j, 
00178                                        points[points.size()-1], point3,
00179                                        updated_neighbours);
00180       // C++ beginners' comment: static_cast to unsigned int is necessary
00181       // to do away with warnings about type mismatch between point3 (int) 
00182       // and points.size (unsigned int)
00183       if (static_cast<unsigned int> (point3) != points.size()-1) {
00184         throw Error("INTERNAL ERROR: point3 != points.size()-1");}
00185     } else {
00186       // update DNN
00187       DNN->RemovePoint(jet_i, updated_neighbours);
00188     }
00189 
00190     // update map
00191     vector<int>::iterator it = updated_neighbours.begin();
00192     for (; it != updated_neighbours.end(); ++it) {
00193       int ii = *it;
00194       _add_ktdistance_to_map(ii, DijMap, DNN.get());
00195     }
00196       
00197   } // end clustering loop 
00198   
00199 }
00200 
00201 
00202 //----------------------------------------------------------------------
00203 /// Add the current kt distance for particle ii to the map (DijMap)
00204 /// using information from the DNN object. Work as follows:
00205 /// 
00206 /// . if the kt is zero then it's nearest neighbour is taken to be the
00207 ///   the beam jet and the distance is zero.
00208 ///
00209 /// . if cylinder distance to nearest neighbour > _Rparam then it is
00210 ///   yiB that is smallest and this is added to map.
00211 ///
00212 /// . otherwise if the nearest neighbour jj has a larger kt then add
00213 ///   dij to the map.
00214 ///
00215 /// . otherwise do nothing
00216 ///
00217 void ClusterSequence::_add_ktdistance_to_map(
00218                           const int & ii, 
00219                           DistMap & DijMap,
00220                           const DynamicNearestNeighbours * DNN) {
00221   
00222   double yiB = jet_scale_for_algorithm(_jets[ii]);
00223   if (yiB == 0.0) {
00224     // in this case convention is that we do not worry about distances
00225     // but directly state that nearest neighbour is beam
00226     DijMap.insert(DijEntry(yiB,  TwoVertices(ii,-1)));
00227   } else {
00228     double DeltaR2 = DNN->NearestNeighbourDistance(ii) * _invR2;
00229     // Logic of following bit is: only add point to map if it has
00230     // smaller kt2 than nearest neighbour j (if it has larger kt,
00231     // then: either it is j's nearest neighbour and then we will
00232     // include dij when we come to j; or it is not j's nearest
00233     // neighbour and j will recombine with someone else).
00234     
00235     // If DeltaR2 > 1.0 then in any case it will recombine with beam rather
00236     // than with any neighbours.
00237     // (put general normalisation here at some point)
00238     if (DeltaR2 > 1.0) {
00239       DijMap.insert(DijEntry(yiB,  TwoVertices(ii,-1)));
00240     } else {
00241       double kt2i = jet_scale_for_algorithm(_jets[ii]);
00242       int jj = DNN->NearestNeighbourIndex(ii);
00243       if (kt2i <= jet_scale_for_algorithm(_jets[jj])) {
00244         double dij = DeltaR2 * kt2i;
00245         DijMap.insert(DijEntry(dij, TwoVertices(ii,jj)));
00246       }
00247     }
00248   }
00249 }
00250 
00251 
00252 FASTJET_END_NAMESPACE
00253 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends