FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
12-boosted_higgs-old.cc
Go to the documentation of this file.
1 //----------------------------------------------------------------------
2 /// \file
3 /// \page Example12old 12 - boosted Higgs tagging (old version)
4 ///
5 /// fastjet example program, illustration of carrying out boosted
6 /// Higgs subjet ID analysis
7 ///
8 /// It illustrates two kinds of functionality:
9 ///
10 /// - following the decomposition of a jet into pieces
11 /// - following information on a b-tag through the jet
12 ///
13 /// This kind of functionality was used in arXiv:0802.2470
14 /// (Butterworth, Davison, Rubin & Salam) for boosted Higgs searches,
15 /// and related functionality was used in arXiv:0806.0848 (Kaplan,
16 /// Rehermann, Schwartz & Tweedie) in searching for boosted tops
17 /// (without b-tag assumptions).
18 ///
19 /// Note that this example is deprecated --- see 12-boosted_higgs.cc
20 /// for the newest version --- so it is not built by default
21 ///
22 /// run it with : ./12-boosted_higgs-old < data/HZ-event-Hmass115.dat
23 ///
24 /// Source code: 12-boosted_higgs-old.cc
25 //----------------------------------------------------------------------
26 
27 
28 //STARTHEADER
29 // $Id: 12-boosted_higgs-old.cc 2577 2011-09-13 15:11:38Z salam $
30 //
31 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
32 //
33 //----------------------------------------------------------------------
34 // This file is part of FastJet.
35 //
36 // FastJet is free software; you can redistribute it and/or modify
37 // it under the terms of the GNU General Public License as published by
38 // the Free Software Foundation; either version 2 of the License, or
39 // (at your option) any later version.
40 //
41 // The algorithms that underlie FastJet have required considerable
42 // development and are described in hep-ph/0512210. If you use
43 // FastJet as part of work towards a scientific publication, please
44 // include a citation to the FastJet paper.
45 //
46 // FastJet is distributed in the hope that it will be useful,
47 // but WITHOUT ANY WARRANTY; without even the implied warranty of
48 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
49 // GNU General Public License for more details.
50 //
51 // You should have received a copy of the GNU General Public License
52 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
53 //----------------------------------------------------------------------
54 //ENDHEADER
55 
56 #include "fastjet/ClusterSequence.hh"
57 #include <iostream> // needed for io
58 #include <sstream> // needed for internal io
59 #include <iomanip>
60 #include <cmath>
61 
62 using namespace std;
63 using namespace fastjet;
64 
65 
66 //----------------------------------------------------------------------
67 // set up a class to give standard (by default E-scheme)
68 // recombination, with additional tracking of flavour information in
69 // the user_index.
70 //
71 // b-tagged particles are assumed to have their user_index set to 1,
72 // and other particles should have user_index to 0.
73 //
74 // Watch out however that, by default, the user_index of a particle is
75 // set to -1 and you may not have control over that (e.g. if you
76 // compute the jet area using explicit ghosts, the ghosts will have a
77 // default user_index of -1). For that reason, if one of the particle
78 // being combined has a user index of -1, we assume it is not b-tagged
79 // (i.e. we count it as 0 in the recombination)
80 //
81 // This will work for native algorithms, but not for all plugins
82 //----------------------------------------------------------------------
84 
85 class FlavourRecombiner : public DefRecomb {
86 public:
87  FlavourRecombiner(RecombinationScheme recomb_scheme = E_scheme) :
88  DefRecomb(recomb_scheme) {};
89 
90  virtual std::string description() const {
91  return DefRecomb::description()+" (with user index addition)";}
92 
93  /// recombine pa and pb and put result into pab
94  virtual void recombine(const PseudoJet & pa, const PseudoJet & pb,
95  PseudoJet & pab) const {
96  DefRecomb::recombine(pa,pb,pab);
97  // Note: see the above discussion for the fact that we consider
98  // negative user indices as "0"
99  pab.set_user_index(max(pa.user_index(),0) + max(pb.user_index(),0));
100  }
101 };
102 
103 
104 //----------------------------------------------------------------------
105 // forward declaration for printing out info about a jet
106 //----------------------------------------------------------------------
107 ostream & operator<<(ostream &, PseudoJet &);
108 
109 
110 //----------------------------------------------------------------------
111 // core of the program
112 //----------------------------------------------------------------------
113 int main (int argc, char ** argv) {
114 
115  vector<PseudoJet> particles;
116 
117  // read in data in format px py pz E b-tag [last of these is optional]
118  // lines starting with "#" are considered as comments and discarded
119  //----------------------------------------------------------
120 
121  string line;
122  while (getline(cin,line)) {
123  if (line.substr(0,1) == "#") {continue;}
124  istringstream linestream(line);
125  double px,py,pz,E;
126  linestream >> px >> py >> pz >> E;
127 
128  // optionally read in btag information
129  int btag;
130  if (! (linestream >> btag)) btag = 0;
131 
132  // construct the particle
133  PseudoJet particle(px,py,pz,E);
134  particle.set_user_index(btag); // btag info goes in user index, for flavour tracking
135  particles.push_back(particle);
136  }
137 
138 
139  // set up the jet finding
140  //
141  // This also shows how to use the "FlavourRecombiner" user-defined
142  // recombiner
143  // ----------------------------------------------------------
144  double R = 1.2;
145  FlavourRecombiner flav_recombiner; // for tracking flavour
146  JetDefinition jet_def(cambridge_algorithm, R, &flav_recombiner);
147 
148 
149  // run the jet finding; find the hardest jet
150  ClusterSequence cs(particles, jet_def);
151  vector<PseudoJet> jets = sorted_by_pt(cs.inclusive_jets());
152 
153  cout << "Ran: " << jet_def.description() << endl << endl;
154  cout << "Hardest jet: " << jets[0] << endl << endl;
155 
156  // now do the subjet decomposition
157  //----------------------------------------------------------
158  //
159  // when unpeeling a C/A jet, often only a very soft piece may break off;
160  // the mass_drop_threshold indicates how much "lighter" the heavier of the two
161  // resulting pieces must be in order for us to consider that we've really
162  // seen some form of substructure
163  double mass_drop_threshold = 0.667;
164  // QCD backgrounds that give larger jet masses have a component
165  // where a quite soft gluon is emitted; to eliminate part of this
166  // one can place a cut on the asymmetry of the branching;
167  //
168  // Here the cut is expressed in terms of y, the kt-distance scaled
169  // to the squared jet mass; an easier way to see it is in terms of
170  // a requirement on the momentum fraction in the splitting: z/(1-z)
171  // and (1-z)/z > rtycut^2 [the correspondence holds only at LO]
172  double rtycut = 0.3;
173 
174  PseudoJet this_jet = jets[0];
175  PseudoJet parent1, parent2;
176  bool had_parents;
177 
178  while ((had_parents = this_jet.has_parents(parent1,parent2))) {
179  // make parent1 the more massive jet
180  if (parent1.m() < parent2.m()) swap(parent1,parent2);
181 
182  // if we pass the conditions on the mass drop and its degree of
183  // asymmetry (z/(1-z) \sim kt_dist/m^2 > rtycut), then we've found
184  // something interesting, so exit the loop
185  if (parent1.m() < mass_drop_threshold * this_jet.m() &&
186  parent1.kt_distance(parent2) > pow(rtycut,2) * this_jet.m2()) {
187  break;
188  } else {
189  // otherwise try a futher decomposition on the more massive jet
190  this_jet = parent1;
191  }
192  }
193 
194  // look to see what we found
195  if (!had_parents) {
196  cout << "Did not find suitable hard substructure in this event." << endl;
197  return 0;
198  }
199 
200  cout << "Found suitable pair of subjets: " << endl;
201  cout << " " << parent1 << endl;
202  cout << " " << parent2 << endl;
203  cout << "Total = " << endl;
204  cout << " " << this_jet << endl << endl;
205 
206  // next we "filter" it, to remove UE & pileup contamination
207  //----------------------------------------------------------
208  //
209  // [there are two ways of doing this; here we directly use the
210  // exsiting cluster sequence and find the exclusive subjets of
211  // this_jet (i.e. work backwards within the cs starting from
212  // this_jet); alternatively one can recluster just the
213  // constituents of the jet]
214  //
215  // first get separation between the subjets (called Rbb -- assuming it's a Higgs!)
216  double Rbb = sqrt(parent1.squared_distance(parent2));
217  double Rfilt = min(Rbb/2, 0.3); // somewhat arbitrary choice
218  unsigned nfilt = 3; // number of pieces we'll take
219  cout << "Subjet separation (Rbb) = " << Rbb << ", Rfilt = " << Rfilt << endl;
220 
221  double dcut = pow(Rfilt/R,2); // for C/A get a view at Rfilt by
222  // using a dcut=(Rfilt/R)^2
223  vector<PseudoJet> filt_subjets = sorted_by_pt(this_jet.exclusive_subjets(dcut));
224 
225  // now print out the filtered jets and reconstruct total
226  // at the same time
227  cout << "Filtered pieces are " << endl;
228  cout << " " << filt_subjets[0] << endl;
229  PseudoJet filtered_total = filt_subjets[0];
230  for (unsigned i = 1; i < nfilt && i < filt_subjets.size(); i++) {
231  cout << " " << filt_subjets[i] << endl;
232  flav_recombiner.plus_equal(filtered_total, filt_subjets[i]);
233  }
234  cout << "Filtered total is " << endl;
235  cout << " " << filtered_total << endl;
236 
237 }
238 
239 
240 //----------------------------------------------------------------------
241 // does the actual work for printing out a jet
242 //----------------------------------------------------------------------
243 ostream & operator<<(ostream & ostr, PseudoJet & jet) {
244  ostr << "pt, y, phi ="
245  << " " << setw(10) << jet.perp()
246  << " " << setw(6) << jet.rap()
247  << " " << setw(6) << jet.phi()
248  << ", mass = " << setw(10) << jet.m()
249  << ", btag = " << jet.user_index();
250  return ostr;
251 }
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
void set_user_index(const int index)
set the user_index, intended to allow the user to add simple identifying information to a particle/je...
Definition: PseudoJet.hh:334
deals with clustering
std::vector< PseudoJet > exclusive_subjets(const double dcut) const
return a vector of all subjets of the current jet (in the sense of the exclusive algorithm) that woul...
Definition: PseudoJet.cc:603
virtual bool has_parents(PseudoJet &parent1, PseudoJet &parent2) const
check if it is the product of a recombination, in which case return the 2 parents through the 'parent...
Definition: PseudoJet.cc:547
double m2() const
returns the squared invariant mass // like CLHEP
Definition: PseudoJet.hh:148
ostream & operator<<(ostream &, PseudoJet &)
does the actual work for printing out a jet
A class that will provide the recombination scheme facilities and/or allow a user to extend these fac...
virtual std::string description() const
return a textual description of the recombination scheme implemented here
int main()
an example program showing how to use fastjet
Definition: 01-basic.cc:50
the FastJet namespace
double squared_distance(const PseudoJet &other) const
returns squared cylinder (rap-phi) distance between this jet and another
Definition: PseudoJet.hh:186
RecombinationScheme
The various recombination schemes.
double kt_distance(const PseudoJet &other) const
returns kt distance (R=1) between this jet and another
Definition: PseudoJet.cc:378
void swap(SharedPtr< T > &a, SharedPtr< T > &b)
swapping
Definition: SharedPtr.hh:384
double phi() const
returns phi (in the range 0..2pi)
Definition: PseudoJet.hh:108
virtual void recombine(const PseudoJet &pa, const PseudoJet &pb, PseudoJet &pab) const
recombine pa and pb and put result into pab
int user_index() const
return the user_index,
Definition: PseudoJet.hh:331
double perp() const
returns the scalar transverse momentum
Definition: PseudoJet.hh:143
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