FastJet  3.1.3
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
fastjet_boosted_higgs.cc
1 //STARTHEADER
2 // $Id: fastjet_boosted_higgs.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, illustration of carrying out boosted
32 // Higgs subjet ID analysis
33 //
34 // It illustrates two kinds of functionality:
35 //
36 // - following the decomposition of a jet into pieces
37 // - following information on a b-tag through the jet
38 //
39 // This kind of functionality was used in arXiv:0802.2470
40 // (Butterworth, Davison, Rubin & Salam) for boosted Higgs searches,
41 // and related functionality was used in arXiv:0806.0848 (Kaplan,
42 // Rehermann, Schwartz & Tweedie) in searching for boosted tops
43 // (without b-tag assumptions).
44 
45 // Compile it with: make fastjet_boosted_higgs
46 // run it with : ./fastjet_boosted_higgs < data/HZ-event-Hmass115.dat
47 //
48 //----------------------------------------------------------------------
49 
50 #include "fastjet/ClusterSequence.hh"
51 #include<iostream> // needed for io
52 #include<sstream> // needed for internal io
53 #include<iomanip>
54 #include<cmath>
55 
56 using namespace std;
57 namespace fj = fastjet;
58 
59 
60 
61 
62 //----------------------------------------------------------------------
63 /// set up a class to give standard (by default E-scheme)
64 /// recombination, with additional tracking of flavour information in
65 /// the user_index.
66 ///
67 /// If you use this, you must explicitly set the user index to 0 for
68 /// non-flavoured particles (the default value is -1);
69 ///
70 /// This will work for native algorithms, but not for all plugins
71 typedef fj::JetDefinition::DefaultRecombiner DefRecomb;
72 class FlavourRecombiner : public DefRecomb {
73 public:
74  FlavourRecombiner(fj::RecombinationScheme recomb_scheme = fj::E_scheme) :
75  DefRecomb(recomb_scheme) {};
76 
77  virtual std::string description() const {return DefRecomb::description()
78  +" (with user index addition)";}
79 
80  /// recombine pa and pb and put result into pab
81  virtual void recombine(const fj::PseudoJet & pa,
82  const fj::PseudoJet & pb,
83  fj::PseudoJet & pab) const {
84  DefRecomb::recombine(pa,pb,pab);
85  pab.set_user_index(pa.user_index() + pb.user_index());
86  }
87 
88 };
89 
90 
91 /// forward declaration for printing out info about a jet
92 ostream & operator<<(ostream &, fj::PseudoJet &);
93 
94 
95 //----------------------------------------------------------------------
96 int main (int argc, char ** argv) {
97 
98  vector<fj::PseudoJet> particles;
99 
100  // read in data in format px py pz E b-tag [last of these is optional]
101  string line;
102  while (getline(cin,line)) {
103  if (line.substr(0,1) == "#") {continue;}
104  istringstream linestream(line);
105  double px,py,pz,E;
106  linestream >> px >> py >> pz >> E;
107 
108  // optionally read in btag information
109  int btag;
110  if (! (linestream >> btag)) btag = 0;
111 
112  // construct the particle
113  fj::PseudoJet particle(px,py,pz,E);
114  particle.set_user_index(btag); // btag info goes in user index, for flavour tracking
115  particles.push_back(particle);
116  }
117 
118 
119  // set up the jet finding
120  double R = 1.2;
121  FlavourRecombiner flav_recombiner; // for tracking flavour
122  fj::JetDefinition jet_def(fj::cambridge_algorithm, R, &flav_recombiner);
123 
124 
125  // run the jet finding; find the hardest jet
126  fj::ClusterSequence cs(particles, jet_def);
127  vector<fj::PseudoJet> jets = sorted_by_pt(cs.inclusive_jets());
128 
129 
130  cout << "Ran: " << jet_def.description() << endl << endl;
131  cout << "Hardest jet: " << jets[0] << endl << endl;
132 
133 
134  /// now do the subjet decomposition;
135  //
136  /// when unpeeling a C/A jet, often only a very soft piece may break off;
137  /// the mass_drop_threshold indicates how much "lighter" the heavier of the two
138  /// resulting pieces must be in order for us to consider that we've really
139  /// seen some form of substructure
140  double mass_drop_threshold = 0.667;
141  /// QCD backgrounds that give larger jet masses have a component
142  /// where a quite soft gluon is emitted; to eliminate part of this
143  /// one can place a cut on the asymmetry of the branching;
144  ///
145  /// Here the cut is expressed in terms of y, the kt-distance scaled
146  /// to the squared jet mass; an easier way to see it is in terms of
147  /// a requirement on the momentum fraction in the splitting: z/(1-z)
148  /// and (1-z)/z > rtycut^2 [the correspondence holds only at LO]
149  double rtycut = 0.3;
150 
151  fj::PseudoJet this_jet = jets[0], parent1, parent2;
152  bool had_parents;
153 
154  while ((had_parents = this_jet.has_parents(parent1,parent2))) {
155  // make parent1 the more massive jet
156  if (parent1.m() < parent2.m()) swap(parent1,parent2);
157  //
158  // if we pass the conditions on the mass drop and its degree of
159  // asymmetry (z/(1-z) \sim kt_dist/m^2 > rtycut), then we've found
160  // something interesting, so exit the loop
161  if (parent1.m() < mass_drop_threshold * this_jet.m() &&
162  parent1.kt_distance(parent2) > pow(rtycut,2) * this_jet.m2()) {
163  break;
164  } else {
165  // otherwise try a futher decomposition on the more massive jet
166  this_jet = parent1;
167  }
168  }
169 
170  // look to see what we found
171  if (had_parents) {
172  cout << "Found suitable pair of subjets: " << endl;
173  cout << " " << parent1 << endl;
174  cout << " " << parent2 << endl;
175  cout << "Total = " << endl;
176  cout << " " << this_jet << endl << endl;
177 
178  // next we "filter" it, to remove UE & pileup contamination
179  //
180  // [there are two ways of doing this; here we directly use the
181  // exsiting cluster sequence and find the exclusive subjets of
182  // this_jet (i.e. work backwards within the cs starting from
183  // this_jet); alternatively one can recluster just the
184  // constituents of the jet]
185  //
186  // first get separation between the subjets (called Rbb -- assuming it's a Higgs!)
187  double Rbb = sqrt(parent1.squared_distance(parent2));
188  double Rfilt = min(Rbb/2, 0.3); // somewhat arbitrary choice
189  unsigned nfilt = 3; // number of pieces we'll take
190  cout << "Subjet separation (Rbb) = " << Rbb << ", Rfilt = " << Rfilt << endl;
191 
192  double dcut = pow(Rfilt/R,2); // for C/A get a view at Rfilt by
193  // using a dcut=(Rfilt/R)^2
194  vector<fj::PseudoJet> filt_subjets = sorted_by_pt(this_jet.exclusive_subjets(dcut));
195 
196  // now print out the filtered jets and reconstruct total
197  // at the same time
198  cout << "Filtered pieces are " << endl;
199  cout << " " << filt_subjets[0] << endl;
200  fj::PseudoJet filtered_total = filt_subjets[0];
201  for (unsigned i = 1; i < nfilt && i < filt_subjets.size(); i++) {
202  cout << " " << filt_subjets[i] << endl;
203  flav_recombiner.plus_equal(filtered_total, filt_subjets[i]);
204  }
205  cout << "Filtered total is " << endl;
206  cout << " " << filtered_total << endl;
207 
208  } else {
209  cout << "Did not find suitable hard substructure in this event." << endl;
210  }
211 }
212 
213 
214 /// does the actual work for printing out a jet
215 ostream & operator<<(ostream & ostr, fj::PseudoJet & jet) {
216  ostr << "pt, y, phi ="
217  << " " << setw(10) << jet.perp()
218  << " " << setw(6) << jet.rap()
219  << " " << setw(6) << jet.phi()
220  << ", mass = " << setw(10) << jet.m()
221  << ", btag = " << jet.user_index();
222  return ostr;
223 }
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
RecombinationScheme
The various recombination schemes.
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