FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
fastjet_subjets.cc
1 //STARTHEADER
2 // $Id: fastjet_subjets.cc 2577 2011-09-13 15:11:38Z salam $
3 //
4 // Copyright (c) 2005-2011, 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 and are described in hep-ph/0512210. If you use
16 // FastJet as part of work towards a scientific publication, please
17 // include a citation to the FastJet paper.
18 //
19 // FastJet is distributed in the hope that it will be useful,
20 // but WITHOUT ANY WARRANTY; without even the implied warranty of
21 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
22 // GNU General Public License for more details.
23 //
24 // You should have received a copy of the GNU General Public License
25 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
26 //----------------------------------------------------------------------
27 //ENDHEADER
28 
29 
30 //----------------------------------------------------------------------
31 // fastjet example program to show how to access subjets;
32 //
33 // See also fastjet_higgs_decomp.cc to see the use of subjets for
34 // identifying boosted higgs (and other objects)
35 //
36 // Compile it with: make fastjet_subjets
37 // run it with : ./fastjet_subjets < data/single-event.dat
38 //
39 //----------------------------------------------------------------------
40 #include "fastjet/PseudoJet.hh"
41 #include "fastjet/ClusterSequence.hh"
42 #include<iostream> // needed for io
43 #include<sstream> // needed for internal io
44 #include<vector>
45 #include <cstdio>
46 
47 using namespace std;
48 
49 // to avoid excessive typing, define an abbreviation for the
50 // fastjet namespace
51 namespace fj = fastjet;
52 
53 
54 // a declaration of a function that pretty prints a list of jets
55 void print_jets (const vector<fj::PseudoJet> &);
56 
57 // and this pretty prinst a single jet
58 void print_jet (const fj::PseudoJet & jet);
59 
60 // pretty print the jets and their subjets
61 void print_jets_and_sub (const vector<fj::PseudoJet> & jets,
62  double dcut);
63 
64 /// an example program showing how to use fastjet
65 int main (int argc, char ** argv) {
66 
67  vector<fj::PseudoJet> input_particles;
68 
69  // read in input particles
70  double px, py , pz, E;
71  while (cin >> px >> py >> pz >> E) {
72  // create a fj::PseudoJet with these components and put it onto
73  // back of the input_particles vector
74  input_particles.push_back(fj::PseudoJet(px,py,pz,E));
75  }
76 
77  // We'll do two subjet analyses, physically rather different
78  double R = 1.0;
79  //fj::JetDefinition kt_def(fj::kt_algorithm, R);
80  fj::JetDefinition cam_def(fj::cambridge_algorithm, R);
81 
82  // run the jet clustering with the above jet definition
83  //fj::ClusterSequence kt_seq(input_particles, kt_def);
84  fj::ClusterSequence cam_seq(input_particles, cam_def);
85 
86  // tell the user what was done
87  cout << "Ran " << cam_def.description() << endl;
88  cout << "Strategy adopted by FastJet was "<<
89  cam_seq.strategy_string()<<endl<<endl;
90 
91  // extract the inclusive jets with pt > 5 GeV
92  double ptmin = 5.0;
93  vector<fj::PseudoJet> inclusive_jets = cam_seq.inclusive_jets(ptmin);
94 
95  // for the subjets
96  double smallR = 0.4;
97  double dcut_cam = pow(smallR/R,2);
98 
99  // print them out
100  cout << "Printing inclusive jets (R = "<<R<<") with pt > "<< ptmin<<" GeV\n";
101  cout << "and their subjets with smallR = " << smallR << "\n";
102  cout << "---------------------------------------\n";
103  print_jets_and_sub(inclusive_jets, dcut_cam);
104  cout << endl;
105 
106 
107  // print them out
108  vector<fj::PseudoJet> exclusive_jets = cam_seq.exclusive_jets(dcut_cam);
109  cout << "Printing exclusive jets with dcut = "<< dcut_cam<<" \n";
110  cout << "--------------------------------------------\n";
111  print_jets(exclusive_jets);
112 
113 
114 }
115 
116 
117 //----------------------------------------------------------------------
118 /// a function that pretty prints a list of jets
119 void print_jets (const vector<fj::PseudoJet> & jets) {
120 
121  // sort jets into increasing pt
122  vector<fj::PseudoJet> sorted_jets = sorted_by_pt(jets);
123 
124  // label the columns
125  printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity",
126  "phi", "pt", "n constituents");
127 
128  // print out the details for each jet
129  for (unsigned int i = 0; i < sorted_jets.size(); i++) {
130  printf("%5u ",i);
131  print_jet(sorted_jets[i]);
132  }
133 
134 }
135 
136 
137 //----------------------------------------------------------------------
138 /// a function that pretty prints a list of jets
139 void print_jets_and_sub (const vector<fj::PseudoJet> & jets,
140  double dcut) {
141 
142  // sort jets into increasing pt
143  vector<fj::PseudoJet> sorted_jets = sorted_by_pt(jets);
144 
145  // label the columns
146  printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity",
147  "phi", "pt", "n constituents");
148 
149  // print out the details for each jet
150  for (unsigned int i = 0; i < sorted_jets.size(); i++) {
151  printf("%5u ",i);
152  print_jet(sorted_jets[i]);
153  vector<fj::PseudoJet> subjets = sorted_by_pt(sorted_jets[i].exclusive_subjets(dcut));
154  for (unsigned int j = 0; j < subjets.size(); j++) {
155  printf(" -sub-%02u ",j);
156  print_jet(subjets[j]);
157  }
158  }
159 
160 }
161 
162 
163 //----------------------------------------------------------------------
164 /// print a single jet
165 void print_jet (const fj::PseudoJet & jet) {
166  int n_constituents = jet.constituents().size();
167  printf("%15.8f %15.8f %15.8f %8u\n",
168  jet.rap(), jet.phi(), jet.perp(), n_constituents);
169 }
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:123
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2
Definition: PseudoJet.cc:799
deals with clustering
int main()
an example program showing how to use fastjet
Definition: 01-basic.cc:50
the FastJet namespace
void print_jets(const vector< fastjet::PseudoJet > &)
a function that pretty prints a list of jets
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:108
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:143
virtual std::vector< PseudoJet > constituents() const
retrieve the constituents.
Definition: PseudoJet.cc:580
Class to contain pseudojets, including minimal information of use to jet-clustering routines...
Definition: PseudoJet.hh:67
class that is intended to hold a full definition of the jet clusterer