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