FastJet  3.3.3
D0RunIIConePlugin.cc
1 //FJSTARTHEADER
2 // $Id: D0RunIIConePlugin.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/D0RunIIConePlugin.hh"
32 #include "fastjet/ClusterSequence.hh"
33 #include "fastjet/Error.hh"
34 #include <sstream>
35 
36 // D0 stuff
37 #include <list>
38 #include "ILConeAlgorithm.hpp"
39 #include "HepEntity.h"
40 
41 FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh
42 
43 using namespace std;
44 using namespace d0;
45 
46 const double D0RunIIConePlugin::_DEFAULT_split_ratio = 0.5 ; // overlap threshold
47 const double D0RunIIConePlugin::_DEFAULT_far_ratio = 0.5 ;
48 const double D0RunIIConePlugin::_DEFAULT_Et_min_ratio = 0.5 ;
49 const bool D0RunIIConePlugin::_DEFAULT_kill_duplicate = true ;
50 const double D0RunIIConePlugin::_DEFAULT_duplicate_dR = 0.005;
51 const double D0RunIIConePlugin::_DEFAULT_duplicate_dPT = 0.01 ;
52 const double D0RunIIConePlugin::_DEFAULT_search_factor = 1.0 ;
53 const double D0RunIIConePlugin::_DEFAULT_pT_min_leading_protojet = 0. ;
54 const double D0RunIIConePlugin::_DEFAULT_pT_min_second_protojet = 0. ;
55 const int D0RunIIConePlugin::_DEFAULT_merge_max = 10000;
56 const double D0RunIIConePlugin::_DEFAULT_pT_min_nomerge = 0. ;
57 
58 bool D0RunIIConePlugin::_first_time = true;
59 
60 string D0RunIIConePlugin::description () const {
61  ostringstream desc;
62 
63  desc << "D0 Run II Improved Legacy (midpoint) cone jet algorithm, with ";
64  desc << "cone_radius = " << cone_radius () << ", "
65  << "min_jet_Et = " << min_jet_Et () << ", "
66  << "split_ratio = " << split_ratio ();
67 
68  return desc.str();
69 }
70 
71 
72 void D0RunIIConePlugin::run_clustering(ClusterSequence & clust_seq) const {
73  // print a banner if we run this for the first time
74  _print_banner(clust_seq.fastjet_banner_stream());
75 
76  // create the entities needed by the D0 code
77  vector<HepEntity> entities(clust_seq.jets().size());
78  list<const HepEntity * > ensemble;
79  for (unsigned i = 0; i < clust_seq.jets().size(); i++) {
80  entities[i].Fill(clust_seq.jets()[i].E(),
81  clust_seq.jets()[i].px(),
82  clust_seq.jets()[i].py(),
83  clust_seq.jets()[i].pz(),
84  i);
85  // use only the particles that do not have infinite rapidity
86  if (abs(entities[i].pz) < entities[i].E) {
87  ensemble.push_back(& (entities[i]));
88  }
89  }
90 
91  // prepare the D0 algorithm
92  ILConeAlgorithm<HepEntity>
93  ilegac(cone_radius(),
94  min_jet_Et(),
95  split_ratio(),
96  far_ratio(),
97  Et_min_ratio(),
98  kill_duplicate(),
99  duplicate_dR(),
100  duplicate_dPT(),
101  search_factor(),
102  pT_min_leading_protojet(),
103  pT_min_second_protojet(),
104  merge_max(),
105  pT_min_nomerge());
106 
107  // run the algorithm
108  float Item_ET_Threshold = 0.;
109  list<HepEntity> jets;
110  ilegac.makeClusters(jets, ensemble, Item_ET_Threshold);
111 
112  // now transfer the information about the jets into the
113  // FastJet structure
114  for(int i = ilegac.ilcv.size()-1; i >= 0; i--) {
115 
116  std::list<const HepEntity*> tlist = ilegac.ilcv[i].LItems();
117  std::list<const HepEntity*>::iterator tk;
118 
119  // get first particle in list
120  tk = tlist.begin();
121 
122  // if there is no particle, just discard it
123  // Note: this unexpected behaviour has been observed when the
124  // min_jet_Et parameter was set to 0
125  if (tk==tlist.end())
126  continue;
127 
128  int jet_k = (*tk)->index;
129  // now merge with remaining particles in list
130  tk++;
131  for (; tk != tlist.end(); tk++) {
132  int jet_i = jet_k;
133  int jet_j = (*tk)->index;
134  // do a fake recombination step with dij=0
135  double dij = 0.0;
136  clust_seq.plugin_record_ij_recombination(jet_i, jet_j, dij, jet_k);
137  }
138 
139  // NB: put a sensible looking d_iB just to be nice...
140  double d_iB = clust_seq.jets()[jet_k].perp2();
141  clust_seq.plugin_record_iB_recombination(jet_k, d_iB);
142 
143  }
144 }
145 
146 // print a banner for reference to the 3rd-party code
147 void D0RunIIConePlugin::_print_banner(ostream *ostr) const{
148  if (! _first_time) return;
149  _first_time=false;
150 
151  // make sure the user has not set the banner stream to NULL
152  if (!ostr) return;
153 
154  (*ostr) << "#--------------------------------------------------------------------------" << endl;
155  (*ostr) << "# You are running the D0 Run II Cone plugin for FastJet " << endl;
156  (*ostr) << "# Original code by the D0 collaboration, provided by Lars Sonnenschein; " << endl;
157  (*ostr) << "# interface by FastJet authors " << endl;
158  (*ostr) << "# If you use this plugin, please cite " << endl;
159  (*ostr) << "# G. C. Blazey et al., hep-ex/0005012 " << endl;
160  (*ostr) << "# V. M. Abazov et al. [D0 Collaboration], arXiv:1110.3771 [hep-ex] " << endl;
161  (*ostr) << "# in addition to the usual FastJet reference. " << endl;
162  (*ostr) << "#--------------------------------------------------------------------------" << endl;
163 
164  // make sure we really have the output done.
165  ostr->flush();
166 }
167 
168 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...
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...