fastjet 2.4.5
ClusterSequence_CP2DChan.cc
Go to the documentation of this file.
00001 //STARTHEADER
00002 // $Id: ClusterSequence_CP2DChan.cc 958 2007-11-28 17:43:48Z cacciari $
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 #include "fastjet/ClusterSequence.hh"
00032 #include "fastjet/internal/ClosestPair2D.hh"
00033 #include<limits>
00034 #include<vector>
00035 #include<cmath>
00036 #include<iostream>
00037 
00038 FASTJET_BEGIN_NAMESPACE      // defined in fastjet/internal/base.hh
00039 
00040 using namespace std;
00041 
00042 // place for things we don't want outside world to run into
00043 namespace Private {
00045   class MirrorInfo{
00046   public:
00047     int orig, mirror;
00048     MirrorInfo(int a, int b) : orig(a), mirror(b) {};
00049     MirrorInfo() {};
00050   };
00051 
00055   bool make_mirror(Coord2D & point, double Dlim) {
00056     if (point.y < Dlim)       {point.y += twopi; return true;}
00057     if (twopi-point.y < Dlim) {point.y -= twopi; return true;}
00058     return false;
00059   }
00060   
00061 }
00062 
00063 using namespace Private;
00064 
00065 
00066 //----------------------------------------------------------------------
00069 void ClusterSequence::_CP2DChan_limited_cluster (double Dlim) {
00070   
00071   unsigned int n = _initial_n;
00072 
00073   vector<MirrorInfo>   coordIDs(2*n); // coord IDs of a given jetID
00074   vector<int>          jetIDs(2*n);   // jet ID for a given coord ID
00075   vector<Coord2D>      coords(2*n);   // our coordinates (and copies)
00076 
00077   // start things off...
00078   double minrap = numeric_limits<double>::max();
00079   double maxrap = -minrap;
00080   int coord_index = -1;
00081   int n_active = 0;
00082   for (unsigned jet_i = 0; jet_i < _jets.size(); jet_i++) {
00083 
00084     // skip jets that already have children or that have infinite
00085     // rapidity
00086     if (_history[_jets[jet_i].cluster_hist_index()].child != Invalid ||
00087         (_jets[jet_i].E() == abs(_jets[jet_i].pz()) && 
00088          _jets[jet_i].perp2() == 0.0)
00089         ) {continue;}
00090 
00091     n_active++;
00092 
00093     coordIDs[jet_i].orig = ++coord_index;
00094     coords[coord_index]  = Coord2D(_jets[jet_i].rap(), _jets[jet_i].phi_02pi());
00095     jetIDs[coord_index]  = jet_i;
00096     minrap = min(coords[coord_index].x,minrap);
00097     maxrap = max(coords[coord_index].x,maxrap);
00098 
00099     Coord2D mirror_point(coords[coord_index]);
00100     if (make_mirror(mirror_point, Dlim)) {
00101       coordIDs[jet_i].mirror = ++coord_index;
00102       coords[coord_index] = mirror_point;
00103       jetIDs[coord_index] = jet_i;
00104     } else {
00105       coordIDs[jet_i].mirror = Invalid;
00106     }
00107   }
00108 
00109   // set them to their actual size...
00110   coords.resize(coord_index+1);
00111 
00112   // establish limits (with some leeway on rapidity)
00113   Coord2D left_edge(minrap-1.0, -pi);
00114   Coord2D right_edge(maxrap+1.0, 3*pi);
00115 
00116   //cerr << "minrap, maxrap = " << minrap << " " << maxrap << endl;
00117 
00118   // now create the closest pair search object
00119   ClosestPair2D cp(coords, left_edge, right_edge);
00120 
00121   // cross check to see what's going on...
00122   //cerr << "Tree depths before: ";
00123   //cp.print_tree_depths(cerr);
00124 
00125   // and start iterating...
00126   vector<Coord2D> new_points(2);
00127   vector<unsigned int> cIDs_to_remove(4);
00128   vector<unsigned int> new_cIDs(2);
00129 
00130   do {
00131     // find out which pair we recombine
00132     unsigned int cID1, cID2;
00133     double distance2;
00134     cp.closest_pair(cID1,cID2,distance2);
00135 
00136     // if the closest distance > Dlim, we exit and go to "inclusive" stage
00137     if (distance2 > Dlim*Dlim) {break;}
00138 
00139     // normalise distance as we like it
00140     distance2 *= _invR2;
00141 
00142     // do the recombination...
00143     int jet_i = jetIDs[cID1];
00144     int jet_j = jetIDs[cID2];
00145     int newjet_k;
00146     _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
00147 
00148     // don't bother with any further action if only one active particle
00149     // is left (also avoid closest-pair error [cannot remove last particle]).
00150     if (--n_active == 1) {break;}
00151 
00152     // now prepare operations on CP structure
00153     cIDs_to_remove.resize(0);
00154     cIDs_to_remove.push_back(coordIDs[jet_i].orig);
00155     cIDs_to_remove.push_back(coordIDs[jet_j].orig);
00156     if (coordIDs[jet_i].mirror != Invalid) 
00157       cIDs_to_remove.push_back(coordIDs[jet_i].mirror);
00158     if (coordIDs[jet_j].mirror != Invalid) 
00159       cIDs_to_remove.push_back(coordIDs[jet_j].mirror);
00160 
00161     Coord2D new_point(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
00162     new_points.resize(0);
00163     new_points.push_back(new_point);
00164     if (make_mirror(new_point, Dlim)) new_points.push_back(new_point);
00165     
00166     // carry out actions on search tree
00167     cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
00168 
00169     // now fill in info for new points...
00170     coordIDs[newjet_k].orig = new_cIDs[0];
00171     jetIDs[new_cIDs[0]]       = newjet_k;
00172     if (new_cIDs.size() == 2) {
00173       coordIDs[newjet_k].mirror = new_cIDs[1];
00174       jetIDs[new_cIDs[1]]         = newjet_k;
00175     } else {coordIDs[newjet_k].mirror = Invalid;}
00176     
00178     //n_active--;
00179     //if (n_active == 1) {break;}
00180 
00181   } while(true);
00182   
00183   // cross check to see what's going on...
00184   //cerr << "Max tree depths after: ";
00185   //cp.print_tree_depths(cerr);
00186 
00187 }
00188 
00189 
00190 //----------------------------------------------------------------------
00193 void ClusterSequence::_CP2DChan_cluster_2pi2R () {
00194 
00195   if (_jet_algorithm != cambridge_algorithm) throw Error("CP2DChan clustering method called for a jet-finder that is not the cambridge algorithm");
00196 
00197   // run the clustering with mirror copies kept such that only things
00198   // within _Rparam of a border are mirrored
00199   _CP2DChan_limited_cluster(_Rparam);
00200 
00201   //
00202   _do_Cambridge_inclusive_jets();
00203 }
00204 
00205 
00206 //----------------------------------------------------------------------
00209 void ClusterSequence::_CP2DChan_cluster_2piMultD () {
00210 
00211   // do a first run of clustering up to a small distance parameter,
00212   if (_Rparam >= 0.39) {
00213     _CP2DChan_limited_cluster(min(_Rparam/2,0.3));
00214   }
00215 
00216   // and then the final round of clustering
00217   _CP2DChan_cluster_2pi2R ();
00218 }
00219 
00220 
00221 //----------------------------------------------------------------------
00223 void ClusterSequence::_CP2DChan_cluster () {
00224 
00225   if (_jet_algorithm != cambridge_algorithm) throw Error("_CP2DChan_cluster called for a jet-finder that is not the cambridge algorithm");
00226 
00227   unsigned int n = _jets.size();
00228 
00229   vector<MirrorInfo>   coordIDs(2*n);  // link from original to mirror indices
00230   vector<int>          jetIDs(2*n);     // link from mirror to original indices
00231   vector<Coord2D>      coords(2*n);   // our coordinates (and copies)
00232 
00233   // start things off...
00234   double minrap = numeric_limits<double>::max();
00235   double maxrap = -minrap;
00236   int coord_index = 0;
00237   for (unsigned i = 0; i < n; i++) {
00238     // separate out points with infinite rapidity
00239     if (_jets[i].E() == abs(_jets[i].pz()) && _jets[i].perp2() == 0.0) {
00240       coordIDs[i] = MirrorInfo(BeamJet,BeamJet);
00241     } else {
00242       coordIDs[i].orig   = coord_index;
00243       coordIDs[i].mirror = coord_index+1;
00244       coords[coord_index]   = Coord2D(_jets[i].rap(), _jets[i].phi_02pi());
00245       coords[coord_index+1] = Coord2D(_jets[i].rap(), _jets[i].phi_02pi()+twopi);
00246       jetIDs[coord_index]   = i;
00247       jetIDs[coord_index+1] = i;
00248       minrap = min(coords[coord_index].x,minrap);
00249       maxrap = max(coords[coord_index].x,maxrap);
00250       coord_index += 2;
00251     }
00252   }
00253   // label remaining "mirror info" as saying that there are no jets
00254   for (unsigned i = n; i < 2*n; i++) {coordIDs[i].orig = Invalid;}
00255 
00256   // set them to their actual size...
00257   coords.resize(coord_index);
00258 
00259   // establish limits (with some leeway on rapidity)
00260   Coord2D left_edge(minrap-1.0, 0.0);
00261   Coord2D right_edge(maxrap+1.0, 2*twopi);
00262 
00263   // now create the closest pair search object
00264   ClosestPair2D cp(coords, left_edge, right_edge);
00265 
00266   // and start iterating...
00267   vector<Coord2D> new_points(2);
00268   vector<unsigned int> cIDs_to_remove(4);
00269   vector<unsigned int> new_cIDs(2);
00270   do {
00271     // find out which pair we recombine
00272     unsigned int cID1, cID2;
00273     double distance2;
00274     cp.closest_pair(cID1,cID2,distance2);
00275     distance2 *= _invR2;
00276 
00277     // if the closest distance > 1, we exit and go to "inclusive" stage
00278     if (distance2 > 1.0) {break;}
00279 
00280     // do the recombination...
00281     int jet_i = jetIDs[cID1];
00282     int jet_j = jetIDs[cID2];
00283     int newjet_k;
00284     _do_ij_recombination_step(jet_i, jet_j, distance2, newjet_k);
00285 
00286     // now prepare operations on CP structure
00287     cIDs_to_remove[0] = coordIDs[jet_i].orig;
00288     cIDs_to_remove[1] = coordIDs[jet_i].mirror;
00289     cIDs_to_remove[2] = coordIDs[jet_j].orig;
00290     cIDs_to_remove[3] = coordIDs[jet_j].mirror;
00291     new_points[0] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi());
00292     new_points[1] = Coord2D(_jets[newjet_k].rap(),_jets[newjet_k].phi_02pi()+twopi);
00293     // carry out the CP operation
00294     //cp.replace_many(cIDs_to_remove, new_points, new_cIDs);
00295     // remarkable the following is faster...
00296     new_cIDs[0] = cp.replace(cIDs_to_remove[0], cIDs_to_remove[2], new_points[0]);
00297     new_cIDs[1] = cp.replace(cIDs_to_remove[1], cIDs_to_remove[3], new_points[1]);
00298     // signal that the following jets no longer exist
00299     coordIDs[jet_i].orig = Invalid;
00300     coordIDs[jet_j].orig = Invalid;
00301     // and do bookkeeping for new jet
00302     coordIDs[newjet_k] = MirrorInfo(new_cIDs[0], new_cIDs[1]);
00303     jetIDs[new_cIDs[0]] = newjet_k;
00304     jetIDs[new_cIDs[1]] = newjet_k;
00305 
00306     // if we've reached one jet we should exit...
00307     n--;
00308     if (n == 1) {break;}
00309 
00310   } while(true);
00311   
00312 
00313   // now do "beam" recombinations 
00314   //for (unsigned int jet_i = 0; jet_i < coordIDs.size(); jet_i++) {
00315   //  // if we have a predeclared beam jet or a valid set of mirror
00316   //  // coordinates, recombine then recombine this jet with the beam
00317   //  if (coordIDs[jet_i].orig == BeamJet || coordIDs[jet_i].orig > 0) {
00318   //    // diB = 1 _always_ (for the cambridge alg.)
00319   //    _do_iB_recombination_step(jet_i, 1.0);
00320   //  }
00321   //}
00322 
00323   _do_Cambridge_inclusive_jets();
00324 
00325   //n = _history.size();
00326   //for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
00327   //  if (_history[hist_i].child == Invalid) {
00328   //    _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
00329   //  }
00330   //}
00331 
00332 }
00333 
00334 
00335 //----------------------------------------------------------------------
00336 void ClusterSequence::_do_Cambridge_inclusive_jets () {
00337   unsigned int n = _history.size();
00338   for (unsigned int hist_i = 0; hist_i < n; hist_i++) {
00339     if (_history[hist_i].child == Invalid) {
00340       _do_iB_recombination_step(_history[hist_i].jetp_index, 1.0);
00341     }
00342   }
00343 }
00344 
00345 FASTJET_END_NAMESPACE
00346 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines