FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
14-groomers.cc
Go to the documentation of this file.
1 //----------------------------------------------------------------------
2 /// \file
3 /// \page Example14 14 - unified use of transformers
4 ///
5 /// fastjet example program to illustrate the use of the
6 /// fastjet::Filter and fastjet::Pruner classes in a unified way
7 /// through their derivation from fastjet::Transformer
8 ///
9 /// The two hardest jets in a boosted top event, clustered with an
10 /// (abnormally) large R, are then groomed using different tools. One
11 /// notes the reduction in the mass of the jets after grooming.
12 ///
13 /// run it with : ./14-groomers < data/boosted_top_event.dat
14 ///
15 /// Source code: 14-groomers.cc
16 //----------------------------------------------------------------------
17 
18 //STARTHEADER
19 // $Id: 14-groomers.cc 3302 2013-12-28 10:15:34Z salam $
20 //
21 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
22 //
23 //----------------------------------------------------------------------
24 // This file is part of FastJet.
25 //
26 // FastJet is free software; you can redistribute it and/or modify
27 // it under the terms of the GNU General Public License as published by
28 // the Free Software Foundation; either version 2 of the License, or
29 // (at your option) any later version.
30 //
31 // The algorithms that underlie FastJet have required considerable
32 // development and are described in hep-ph/0512210. If you use
33 // FastJet as part of work towards a scientific publication, please
34 // include a citation to the FastJet paper.
35 //
36 // FastJet is distributed in the hope that it will be useful,
37 // but WITHOUT ANY WARRANTY; without even the implied warranty of
38 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
39 // GNU General Public License for more details.
40 //
41 // You should have received a copy of the GNU General Public License
42 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
43 //----------------------------------------------------------------------
44 //ENDHEADER
45 
46 #include <fastjet/PseudoJet.hh>
47 #include <fastjet/ClusterSequence.hh>
48 #include <fastjet/Selector.hh>
49 #include <iostream>
50 #include "fastjet/tools/Filter.hh"
51 #include "fastjet/tools/Pruner.hh"
52 
53 #include <cstdio> // needed for io
54 
55 using namespace fastjet;
56 using namespace std;
57 
58 /// an example program showing how to use Filter and Pruner in FastJet
59 int main(){
60  // read in input particles
61  //----------------------------------------------------------
62  vector<PseudoJet> input_particles;
63 
64  double px, py , pz, E;
65  while (cin >> px >> py >> pz >> E) {
66  // create a PseudoJet with these components and put it onto
67  // back of the input_particles vector
68  input_particles.push_back(PseudoJet(px,py,pz,E));
69  }
70 
71  // get the resulting jets ordered in pt
72  //----------------------------------------------------------
73  JetDefinition jet_def(cambridge_algorithm, 1.5);
74  ClusterSequence clust_seq(input_particles, jet_def);
75  vector<PseudoJet> inclusive_jets =
76  sorted_by_pt(clust_seq.inclusive_jets(5.0));
77 
78  // label the columns
79  printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "mass");
80 
81  // print out the details for each jet
82  for (unsigned int i = 0; i < inclusive_jets.size(); i++) {
83  printf("%5u %15.8f %15.8f %15.8f %15.8f\n",
84  i, inclusive_jets[i].rap(), inclusive_jets[i].phi(),
85  inclusive_jets[i].perp(),inclusive_jets[i].m());
86  }
87 
88  // simple test to avoid that the example below crashes:
89  // make sure there is at least 3 jets above our 5 GeV
90  if (inclusive_jets.size()<2){
91  cout << "Please provide an event with at least 2 jets above 5 GeV" << endl;
92  return 1;
93  }
94 
95  // We will groom the two hardest jets of the event
96  //----------------------------------------------------------
97  vector<PseudoJet> candidates = SelectorNHardest(2)(inclusive_jets);
98 
99  // create 3 groomers
100  //----------------------------------------------------------
101  vector<Transformer *> groomers;
102 
103  // 1.
104  // the Cambridge/Aachen filter with Rfilt=0.3
105  // (simplified version of arXiv:0802.2470)
106  double Rfilt = 0.3;
107  unsigned int nfilt = 3;
108  groomers.push_back(new Filter(JetDefinition(cambridge_algorithm, Rfilt),
109  SelectorNHardest(nfilt) ) );
110 
111  // 2.
112  // Filtering with a pt cut as for trimming (arXiv:0912.1342)
113  double Rtrim = 0.2;
114  double ptfrac = 0.03;
115  groomers.push_back(new Filter(JetDefinition(kt_algorithm, Rtrim),
116  SelectorPtFractionMin(ptfrac) ) );
117 
118  // 3.
119  // Pruning (arXiv:0903.5081)
120  double zcut = 0.1;
121  double rcut_factor = 0.5;
122  groomers.push_back(new Pruner(cambridge_algorithm, zcut, rcut_factor));
123 
124  // apply the various groomers to the test PseudoJet's
125  // and show the result
126  //----------------------------------------------------------
127 
128  // print out original jet candidates
129  cout << "\nOriginal jets that will be grooomed: " << endl;
130  for (vector<PseudoJet>::iterator jit=candidates.begin(); jit!=candidates.end(); jit++){
131  const PseudoJet & c = *jit;
132  cout << " rap = " << c.rap() << ", phi = " << c.phi() << ", pt = " << c.perp()
133  << ", mass = " << c.m()
134  << " [" << c.description() << "]" << endl;
135  }
136 
137  // loop on groomers
138  for (unsigned int i=0; i < groomers.size(); i++){
139  const Transformer & f = *groomers[i];
140  cout << "\nUsing groomer: " << f.description() << endl;
141 
142  // loop on jet candidates
143  for (vector<PseudoJet>::iterator jit=candidates.begin(); jit!=candidates.end(); jit++){
144  const PseudoJet & c = *jit;
145 
146  // apply groomer f to jet c
147  PseudoJet j = f(c);
148 
149  // access properties specific to the given transformer
150  //
151  // We first make sure that the jet indeed has a structure
152  // compatible with the result of a Filter or Pruner (using
153  // has_structure_of()), and then retrieve the pieces rejected by the
154  // groomer (using structure_of())
155  int n_rejected;
156  if (j.has_structure_of<Filter>()) {
157  const Filter::StructureType & fj_struct = j.structure_of<Filter>();
158  n_rejected = fj_struct.rejected().size();
159  }else {
160  assert(j.has_structure_of<Pruner>()); // make sure
161  const Pruner::StructureType & fj_struct = j.structure_of<Pruner>();
162  n_rejected = fj_struct.rejected().size();
163  }
164 
165  // write out result
166  cout << " rap = " << j.rap() << ", phi = " << j.phi() << ", pt = " << j.perp()
167  << " mass = " << j.m()
168  << " [kept: " << j.pieces().size()
169  << ", rejected: " << n_rejected;
170  if (j.has_structure_of<Pruner>()) {
171  cout << ", Rcut: " << j.structure_of<Pruner>().Rcut();
172  }
173  cout << "]" << endl;
174  }
175  }
176 
177  // a bit of memory cleaning
178  for (unsigned int i=0; i < groomers.size(); i++) delete groomers[i];
179 
180  return 0;
181 }
double rap() const
returns the rapidity or some large value when the rapidity is infinite
Definition: PseudoJet.hh:123
std::string description() const
return a string describing what kind of PseudoJet we are dealing with
Definition: PseudoJet.cc:409
Selector SelectorPtFractionMin(double fraction)
select objects that carry at least a fraction "fraction" of the reference jet.
Definition: Selector.cc:1365
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
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:981
virtual std::string description() const =0
This should be overloaded to return a description of the Transformer.
Selector SelectorNHardest(unsigned int n)
select the n hardest objects
Definition: Selector.cc:1074
std::vector< PseudoJet > inclusive_jets(const double ptmin=0.0) const
return a vector of all jets (in the sense of the inclusive algorithm) with pt >= ptmin.
Class that helps perform filtering (Butterworth, Davison, Rubin and Salam, arXiv:0802.2470) and trimming (Krohn, Thaler and Wang, arXiv:0912.1342) on jets, optionally in conjunction with subtraction (Cacciari and Salam, arXiv:0707.1378).
Definition: Filter.hh:97
Base (abstract) class for a jet transformer.
Definition: Transformer.hh:71
the longitudinally invariant kt algorithm
std::vector< PseudoJet > rejected() const
return the constituents that have been rejected
Definition: Pruner.hh:200
virtual std::vector< PseudoJet > pieces() const
retrieve the pieces that make up the jet.
Definition: PseudoJet.cc:680
Transformer that prunes a jet.
Definition: Pruner.hh:107
int main()
an example program showing how to use Filter and Pruner in FastJet
Definition: 14-groomers.cc:59
The structure associated with a PseudoJet thas has gone through a Pruner transformer.
Definition: Pruner.hh:189
the FastJet namespace
Class to contain structure information for a filtered jet.
Definition: Filter.hh:199
const std::vector< PseudoJet > & rejected() const
returns the subjets that were not kept during the filtering procedure (subtracted if the filter reque...
Definition: Filter.hh:224
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
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
bool has_structure_of() const
check if the PseudoJet has the structure resulting from a Transformer (that is, its structure is comp...
Definition: PseudoJet.hh:971
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP) ...
Definition: PseudoJet.hh:920
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