FastJet  3.3.0
13-boosted_top.cc
Go to the documentation of this file.
1 //----------------------------------------------------------------------
2 /// \file
3 /// \page Example13 13 - boosted top tagging
4 ///
5 /// fastjet example program, illustration of carrying out boosted
6 /// top subjet ID analysis using the Johns Hopkins top tagger as
7 /// introduced in arXiv:0806.0848 (Kaplan, Rehermann, Schwartz
8 /// and Tweedie)
9 ///
10 /// run it with : ./13-boosted_top < data/boosted_top_event.dat
11 ///
12 /// Source code: 13-boosted_top.cc
13 //----------------------------------------------------------------------
14 
15 
16 //STARTHEADER
17 // $Id: 13-boosted_top.cc 2684 2011-11-14 07:41:44Z soyez $
18 //
19 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
20 //
21 //----------------------------------------------------------------------
22 // This file is part of FastJet.
23 //
24 // FastJet is free software; you can redistribute it and/or modify
25 // it under the terms of the GNU General Public License as published by
26 // the Free Software Foundation; either version 2 of the License, or
27 // (at your option) any later version.
28 //
29 // The algorithms that underlie FastJet have required considerable
30 // development and are described in hep-ph/0512210. If you use
31 // FastJet as part of work towards a scientific publication, please
32 // include a citation to the FastJet paper.
33 //
34 // FastJet is distributed in the hope that it will be useful,
35 // but WITHOUT ANY WARRANTY; without even the implied warranty of
36 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
37 // GNU General Public License for more details.
38 //
39 // You should have received a copy of the GNU General Public License
40 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
41 //----------------------------------------------------------------------
42 //ENDHEADER
43 
44 #include <iostream> // needed for io
45 #include <sstream> // needed for internal io
46 #include <iomanip>
47 #include <cmath>
48 
49 #include <fastjet/ClusterSequence.hh>
50 #include <fastjet/Selector.hh>
51 #include <fastjet/tools/JHTopTagger.hh>
52 
53 using namespace std;
54 using namespace fastjet;
55 
56 
57 //----------------------------------------------------------------------
58 // forward declaration for printing out info about a jet
59 //----------------------------------------------------------------------
60 ostream & operator<<(ostream &, const PseudoJet &);
61 
62 //----------------------------------------------------------------------
63 // core of the program
64 //----------------------------------------------------------------------
65 int main(){
66 
67  vector<PseudoJet> particles;
68 
69  // read in data in format px py pz E b-tag [last of these is optional]
70  // lines starting with "#" are considered as comments and discarded
71  //----------------------------------------------------------
72  string line;
73  while (getline(cin,line)) {
74  if (line.substr(0,1) == "#") {continue;}
75  istringstream linestream(line);
76  double px,py,pz,E;
77  linestream >> px >> py >> pz >> E;
78 
79  // construct the particle
80  particles.push_back(PseudoJet(px,py,pz,E));
81  }
82 
83  // compute the parameters to be used through the analysis
84  // ----------------------------------------------------------
85  double Et=0;
86  for (unsigned int i=0; i<particles.size(); i++)
87  Et += particles[i].perp();
88 
89  double R, delta_p, delta_r;
90  if (Et>2600){ R=0.4; delta_p=0.05; delta_r=0.19;}
91  else if (Et>1600){ R=0.6; delta_p=0.05; delta_r=0.19;}
92  else if (Et>1000){ R=0.8; delta_p=0.10; delta_r=0.19;}
93  else{ cerr << "Et has to be at least 1 TeV"<< endl; return 1;}
94 
95  double ptmin = min(500.0, 0.7*Et/2);
96 
97  // find the jets
98  // ----------------------------------------------------------
99  JetDefinition jet_def(cambridge_algorithm, R);
100  ClusterSequence cs(particles, jet_def);
101  vector<PseudoJet> jets = sorted_by_pt(cs.inclusive_jets());
102 
103  cout << "Ran: " << jet_def.description() << endl << endl;
104  cout << "2 Hardest jets: " << jets[0] << endl
105  << " " << jets[1] << endl << endl;
106 
107  if (jets[0].perp()<ptmin){
108  cout << "No jet above the ptmin threshold" << endl;
109  return 2;
110  }
111 
112  // now do jet tagging using the Johns Hopkins top tagger
113  // For simplicity, we just apply it to the hardest jet.
114  //
115  // In addition to delta_p and delta_r, note that there are two
116  // further parameters to the JH top tagger that here are implicitly
117  // set to their defaults:
118  //
119  // - cos_theta_W_max (defaults to 0.7)
120  // - mW (defaults to 80.4).
121  //
122  // The value for mW implicitly assumes that momenta are passed in
123  // GeV.
124  // ----------------------------------------------------------
125  JHTopTagger top_tagger(delta_p, delta_r);
126  top_tagger.set_top_selector(SelectorMassRange(150,200));
127  top_tagger.set_W_selector (SelectorMassRange( 65, 95));
128 
129  PseudoJet tagged = top_tagger(jets[0]);
130 
131  cout << "Ran the following top tagger: " << top_tagger.description() << endl;
132 
133  if (tagged == 0){
134  cout << "No top substructure found" << endl;
135  return 0;
136  }
137 
138  cout << "Found top substructure from the hardest jet:" << endl;
139  cout << " top candidate: " << tagged << endl;
140  cout << " |_ W candidate: " << tagged.structure_of<JHTopTagger>().W() << endl;
141  cout << " | |_ W subjet 1: " << tagged.structure_of<JHTopTagger>().W1() << endl;
142  cout << " | |_ W subjet 2: " << tagged.structure_of<JHTopTagger>().W2() << endl;
143  cout << " | cos(theta_W) = " << tagged.structure_of<JHTopTagger>().cos_theta_W() << endl;
144  cout << " |_ non-W subjet: " << tagged.structure_of<JHTopTagger>().non_W() << endl;
145 }
146 
147 
148 //----------------------------------------------------------------------
149 // does the actual work for printing out a jet
150 //----------------------------------------------------------------------
151 ostream & operator<<(ostream & ostr, const PseudoJet & jet) {
152  ostr << "pt, y, phi =" << setprecision(6)
153  << " " << setw(9) << jet.perp()
154  << " " << setw(9) << jet.rap()
155  << " " << setw(9) << jet.phi()
156  << ", mass = " << setw(9) << jet.m();
157  return ostr;
158 }
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:774
deals with clustering
Class that helps perform boosted top tagging using the "Johns Hopkins" method from arXiv:0806...
Definition: JHTopTagger.hh:120
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:145
int main()
an example program showing how to use fastjet
Definition: 01-basic.cc:50
the FastJet namespace
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) ...
Definition: PseudoJet.hh:951
const TransformerType::StructureType & structure_of() const
this is a helper to access any structure created by a Transformer (that is, of type Transformer::Stru...
Definition: PseudoJet.hh:1012
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:110
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:125
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
Selector SelectorMassRange(double mmin, double mmax)
select objects with Mmin <= Mass <= Mmax
Definition: Selector.cc:766
class that is intended to hold a full definition of the jet clusterer