FastJet  3.3.3
CDFMidPointPlugin.cc
1 //FJSTARTHEADER
2 // $Id: CDFMidPointPlugin.cc 4420 2019-11-29 09:28:20Z soyez $
3 //
4 // Copyright (c) 2005-2019, 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 #include "fastjet/CDFMidPointPlugin.hh"
32 #include "fastjet/ClusterSequence.hh"
33 #include "fastjet/Error.hh"
34 #include <sstream>
35 
36 // CDF stuff
37 #include "MidPointAlgorithm.hh"
38 #include "PhysicsTower.hh"
39 #include "Cluster.hh"
40 
41 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
42 
43 using namespace std;
44 using namespace cdf;
45 
46 bool CDFMidPointPlugin::_first_time = true;
47 
48 string CDFMidPointPlugin::description () const {
49  ostringstream desc;
50 
51  string sm_scale_string = "split-merge uses ";
52  switch(_sm_scale) {
53  case SM_pt:
54  sm_scale_string += "pt";
55  break;
56  case SM_Et:
57  sm_scale_string += "Et";
58  break;
59  case SM_mt:
60  sm_scale_string += "mt";
61  break;
62  case SM_pttilde:
63  sm_scale_string += "pttilde (scalar sum of pts)";
64  break;
65  default:
66  ostringstream err;
67  err << "Unrecognized split-merge scale choice = " << _sm_scale;
68  throw Error(err.str());
69  }
70 
71 
72  if (cone_area_fraction() == 1) {
73  desc << "CDF MidPoint jet algorithm, with " ;
74  } else {
75  desc << "CDF MidPoint+Searchcone jet algorithm, with ";
76  }
77  desc << "seed_threshold = " << seed_threshold () << ", "
78  << "cone_radius = " << cone_radius () << ", "
79  << "cone_area_fraction = " << cone_area_fraction () << ", "
80  << "max_pair_size = " << max_pair_size () << ", "
81  << "max_iterations = " << max_iterations () << ", "
82  << "overlap_threshold = " << overlap_threshold () << ", "
83  << sm_scale_string ;
84 
85  return desc.str();
86 }
87 
88 
89 void CDFMidPointPlugin::run_clustering(ClusterSequence & clust_seq) const {
90  // print a banner if we run this for the first time
91  _print_banner(clust_seq.fastjet_banner_stream());
92 
93  // create the physics towers needed by the CDF code
94  vector<PhysicsTower> towers;
95  towers.reserve(clust_seq.jets().size());
96  for (unsigned i = 0; i < clust_seq.jets().size(); i++) {
97  LorentzVector fourvect(clust_seq.jets()[i].px(),
98  clust_seq.jets()[i].py(),
99  clust_seq.jets()[i].pz(),
100  clust_seq.jets()[i].E());
101  PhysicsTower tower(fourvect);
102  // misuse one of the indices for tracking, since the MidPoint
103  // implementation doesn't seem to make use of these indices
104  tower.calTower.iEta = i;
105  towers.push_back(tower);
106  }
107 
108  // prepare the CDF algorithm
109  MidPointAlgorithm m(_seed_threshold,_cone_radius,_cone_area_fraction,
110  _max_pair_size,_max_iterations,_overlap_threshold,
111  MidPointAlgorithm::SplitMergeScale(_sm_scale));
112 
113  // run the CDF algorithm
114  std::vector<Cluster> jets;
115  m.run(towers,jets);
116 
117 
118  // now transfer the jets back into our own structure -- we will
119  // mimic the cone code with a sequential recombination sequence in
120  // which the jets are built up by adding one particle at a time
121  for(vector<Cluster>::const_iterator jetIter = jets.begin();
122  jetIter != jets.end(); jetIter++) {
123  const vector<PhysicsTower> & tower_list = jetIter->towerList;
124  int jet_k = tower_list[0].calTower.iEta;
125 
126  int ntow = int(jetIter->towerList.size());
127  for (int itow = 1; itow < ntow; itow++) {
128  int jet_i = jet_k;
129  // retrieve our misappropriated index for the jet
130  int jet_j = tower_list[itow].calTower.iEta;
131  // do a fake recombination step with dij=0
132  double dij = 0.0;
133  clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
134  }
135 
136  // NB: put a sensible looking d_iB just to be nice...
137  double d_iB = clust_seq.jets()[jet_k].perp2();
138  clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
139  }
140 
141 
142  // following code is for testing only
143  //cout << endl;
144  //for(vector<Cluster>::const_iterator jetIter = jets.begin();
145  // jetIter != jets.end(); jetIter++) {
146  // cout << jetIter->fourVector.pt() << " " << jetIter->fourVector.y() << endl;
147  //}
148  //cout << "-----------------------------------------------------\n";
149  //vector<PseudoJet> ourjets(clust_seq.inclusive_jets());
150  //for (vector<PseudoJet>::const_reverse_iterator ourjet = ourjets.rbegin();
151  // ourjet != ourjets.rend(); ourjet++) {
152  // cout << ourjet->perp() << " " << ourjet->rap() << endl;
153  //}
154  //cout << endl;
155 }
156 
157 // print a banner for reference to the 3rd-party code
158 void CDFMidPointPlugin::_print_banner(ostream *ostr) const{
159  if (! _first_time) return;
160  _first_time=false;
161 
162  // make sure the user has not set the banner stream to NULL
163  if (!ostr) return;
164 
165  (*ostr) << "#-------------------------------------------------------------------------" << endl;
166  (*ostr) << "# You are running the CDF MidPoint plugin for FastJet " << endl;
167  (*ostr) << "# This is based on an implementation provided by Joey Huston. " << endl;
168  (*ostr) << "# If you use this plugin, please cite " << endl;
169  (*ostr) << "# G. C. Blazey et al., hep-ex/0005012. " << endl;
170  (*ostr) << "# in addition to the usual FastJet reference. " << endl;
171  (*ostr) << "#-------------------------------------------------------------------------" << endl;
172 
173  // make sure we really have the output done.
174  ostr->flush();
175 }
176 
177 FASTJET_END_NAMESPACE // defined in fastjet/internal/base.hh
deals with clustering
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...
base class corresponding to errors that can be thrown by FastJet
Definition: Error.hh:47
static std::ostream * fastjet_banner_stream()
returns a pointer to the stream to be used to print banners (cout by default).
const std::vector< PseudoJet > & jets() const
allow the user to access the internally stored _jets() array, which contains both the initial particl...