FastJet  3.2.1
EECambridgePlugin.cc
1 //FJSTARTHEADER
2 // $Id: EECambridgePlugin.cc 3433 2014-07-23 08:17:03Z salam $
3 //
4 // Copyright (c) 2007-2014, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
5 //
6 //----------------------------------------------------------------------
7 // This file is part of FastJet.
8 //
9 // FastJet is free software; you can redistribute it and/or modify
10 // it under the terms of the GNU General Public License as published by
11 // the Free Software Foundation; either version 2 of the License, or
12 // (at your option) any later version.
13 //
14 // The algorithms that underlie FastJet have required considerable
15 // development. They are described in the original FastJet paper,
16 // hep-ph/0512210 and in the manual, arXiv:1111.6097. If you use
17 // FastJet as part of work towards a scientific publication, please
18 // quote the version you use and include a citation to the manual and
19 // optionally also to hep-ph/0512210.
20 //
21 // FastJet is distributed in the hope that it will be useful,
22 // but WITHOUT ANY WARRANTY; without even the implied warranty of
23 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
24 // GNU General Public License for more details.
25 //
26 // You should have received a copy of the GNU General Public License
27 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
28 //----------------------------------------------------------------------
29 //FJENDHEADER
30 
31 // fastjet stuff
32 #include "fastjet/ClusterSequence.hh"
33 #include "fastjet/EECambridgePlugin.hh"
34 #include "fastjet/NNH.hh"
35 
36 // other stuff
37 #include <sstream>
38 #include <limits>
39 
40 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
41 
42 using namespace std;
43 
44 //----------------------------------------------------------------------
45 /// class to help run an e+e- Cambridge algorithm
46 class EECamBriefJet {
47 public:
48  void init(const PseudoJet & jet) {
49  double norm = 1.0/sqrt(jet.modp2());
50  nx = jet.px() * norm;
51  ny = jet.py() * norm;
52  nz = jet.pz() * norm;
53  }
54 
55  double distance(const EECamBriefJet * jet) const {
56  double dij = 1 - nx*jet->nx
57  - ny*jet->ny
58  - nz*jet->nz;
59  return dij;
60  }
61 
62  double beam_distance() const {
63  return numeric_limits<double>::max();
64  }
65 
66 private:
67  double nx, ny, nz;
68 };
69 
70 
71 string EECambridgePlugin::description () const {
72  ostringstream desc;
73  desc << "EECambridge plugin with ycut = " << ycut() ;
74  return desc.str();
75 }
76 
77 void EECambridgePlugin::run_clustering(ClusterSequence & cs) const {
78  int njets = cs.jets().size();
79  NNH<EECamBriefJet> nnh(cs.jets());
80 
81  double Q2 = cs.Q2();
82 
83  while (njets > 0) {
84  int i, j, k;
85  // here we get a minimum based on the purely angular variable from the NNH class
86  // (called dij there, but vij in the Cambridge article (which uses dij for
87  // a kt distance...)
88  double vij = nnh.dij_min(i, j); // i,j are return values...
89 
90  // next we work out the dij (ee kt distance), and based on its
91  // value decide whether we have soft-freezing (represented here by
92  // a "Beam" clustering) or not
93  double dij;
94  if (j >= 0) {
95  double scale = min(cs.jets()[i].E(), cs.jets()[j].E());
96  dij = 2 * vij * scale * scale;
97  if (dij > Q2 * ycut()) {
98  // we'll call the softer partner a "beam" jet
99  if (cs.jets()[i].E() > cs.jets()[j].E()) std::swap(i,j);
100  j = -1;
101  }
102  } else {
103  // for the last particle left, just use yij = 1
104  dij = Q2;
105  }
106 
107  if (j >= 0) {
108  cs.plugin_record_ij_recombination(i, j, dij, k);
109  nnh.merge_jets(i, j, cs.jets()[k], k);
110  } else {
112  nnh.remove_jet(i);
113  }
114  njets--;
115  }
116 
117 }
118 
119 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...
double norm(const VPoint p)
norm of a vector
Definition: Voronoi.hh:155
deals with clustering
double Q2() const
return Q()^2
Help solve closest pair problems with generic interparticle and beam distance (generic case) ...
Definition: NNH.hh:75
void plugin_record_ij_recombination(int jet_i, int jet_j, double dij, int &newjet_k)
record the fact that there has been a recombination between jets()[jet_i] and jets()[jet_k], with the specified dij, and return the index (newjet_k) allocated to the new jet, whose momentum is assumed to be the 4-vector sum of that of jet_i and jet_j
void plugin_record_iB_recombination(int jet_i, double diB)
record the fact that there has been a recombination between jets()[jet_i] and the beam...