FastJet  3.3.0
10-subjets.cc
Go to the documentation of this file.
1 //----------------------------------------------------------------------
2 /// \file
3 /// \page Example10 10 - extracting subjets
4 ///
5 /// fastjet example program to show how to access subjets;
6 ///
7 /// See also 12-boosted_higgs.cc to see the use of subjets for
8 /// identifying boosted higgs (and other objects)
9 ///
10 /// run it with : ./10-subjets < data/single-event.dat
11 ///
12 /// Source code: 10-subjets.cc
13 //----------------------------------------------------------------------
14 
15 //STARTHEADER
16 // $Id: 10-subjets.cc 2742 2011-11-22 22:27:08Z salam $
17 //
18 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
19 //
20 //----------------------------------------------------------------------
21 // This file is part of FastJet.
22 //
23 // FastJet is free software; you can redistribute it and/or modify
24 // it under the terms of the GNU General Public License as published by
25 // the Free Software Foundation; either version 2 of the License, or
26 // (at your option) any later version.
27 //
28 // The algorithms that underlie FastJet have required considerable
29 // development and are described in hep-ph/0512210. If you use
30 // FastJet as part of work towards a scientific publication, please
31 // include a citation to the FastJet paper.
32 //
33 // FastJet is distributed in the hope that it will be useful,
34 // but WITHOUT ANY WARRANTY; without even the implied warranty of
35 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
36 // GNU General Public License for more details.
37 //
38 // You should have received a copy of the GNU General Public License
39 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
40 //----------------------------------------------------------------------
41 //ENDHEADER
42 
43 #include "fastjet/ClusterSequence.hh"
44 #include <iostream> // needed for io
45 #include <cstdio> // needed for io
46 
47 using namespace std;
48 using namespace fastjet;
49 
50 int main(){
51 
52  // read in input particles
53  //----------------------------------------------------------
54  vector<PseudoJet> input_particles;
55 
56  double px, py , pz, E;
57  while (cin >> px >> py >> pz >> E) {
58  // create a PseudoJet with these components and put it onto
59  // back of the input_particles vector
60  input_particles.push_back(PseudoJet(px,py,pz,E));
61  }
62 
63 
64  // create a jet definition:
65  // for subjet studies, Cambridge/Aachen is the natural algorithm
66  //----------------------------------------------------------
67  double R = 1.0;
68  JetDefinition jet_def(cambridge_algorithm, R);
69 
70 
71  // run the jet clustering with the above jet definition
72  // and get the jets above 5 GeV
73  //----------------------------------------------------------
74  ClusterSequence clust_seq(input_particles, jet_def);
75  double ptmin = 6.0;
76  vector<PseudoJet> inclusive_jets = sorted_by_pt(clust_seq.inclusive_jets(ptmin));
77 
78  // extract the subjets at a smaller angular scale (Rsub=0.5)
79  //
80  // This is done by using ClusterSequence::exclusive_subjets(dcut):
81  // for the Cambridge/Aachen algorithm, running with R and then
82  // asking for exclusive subjets with dcut should give the same
83  // subjets as rerunning the algorithm with R'=R*sqrt(dcut) on the
84  // jet's constituents.
85  //
86  // At the same time we output a summary of what has been done and the
87  // resulting subjets
88  //----------------------------------------------------------
89  double Rsub = 0.5;
90  double dcut = pow(Rsub/R,2);
91 
92  // a "header" for the output
93  cout << "Ran " << jet_def.description() << endl;
94  cout << "Showing the jets above " << ptmin << " GeV" << endl;
95  cout << "And their subjets for Rsub = " << Rsub << endl;
96  printf("%10s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "n constituents");
97 
98  // show the jets and their subjets
99  for (unsigned int i = 0; i < inclusive_jets.size(); i++) {
100  // get the subjets
101  vector<PseudoJet> subjets = sorted_by_pt(inclusive_jets[i].exclusive_subjets(dcut));
102 
103  cout << endl;
104  // print the jet and its subjets
105  printf("%5u %15.8f %15.8f %15.8f %8d\n", i,
106  inclusive_jets[i].rap(), inclusive_jets[i].phi(),
107  inclusive_jets[i].perp(), int(inclusive_jets[i].constituents().size()));
108 
109  for (unsigned int j=0; j<subjets.size(); j++)
110  printf(" sub%4u %15.8f %15.8f %15.8f %8u\n", j,
111  subjets[j].rap(), subjets[j].phi(),
112  subjets[j].perp(),
113  (unsigned int) subjets[j].constituents().size());
114  }
115 
116  return 0;
117 }
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
int main()
an example program showing how to use fastjet
Definition: 01-basic.cc:50
the FastJet namespace
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