FastJet  3.1.0-beta.1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
12-boosted_higgs.cc
Go to the documentation of this file.
1 //----------------------------------------------------------------------
2 /// \file
3 /// \page Example12 12 - boosted Higgs tagging
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 /// - using a boosted higgs tagger
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 /// run it with : ./12-boosted_higgs < data/HZ-event-Hmass115.dat
20 ///
21 /// Source code: 12-boosted_higgs.cc
22 //----------------------------------------------------------------------
23 
24 
25 //STARTHEADER
26 // $Id: 12-boosted_higgs.cc 2684 2011-11-14 07:41:44Z soyez $
27 //
28 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
29 //
30 //----------------------------------------------------------------------
31 // This file is part of FastJet.
32 //
33 // FastJet is free software; you can redistribute it and/or modify
34 // it under the terms of the GNU General Public License as published by
35 // the Free Software Foundation; either version 2 of the License, or
36 // (at your option) any later version.
37 //
38 // The algorithms that underlie FastJet have required considerable
39 // development and are described in hep-ph/0512210. If you use
40 // FastJet as part of work towards a scientific publication, please
41 // include a citation to the FastJet paper.
42 //
43 // FastJet is distributed in the hope that it will be useful,
44 // but WITHOUT ANY WARRANTY; without even the implied warranty of
45 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
46 // GNU General Public License for more details.
47 //
48 // You should have received a copy of the GNU General Public License
49 // along with FastJet. If not, see <http://www.gnu.org/licenses/>.
50 //----------------------------------------------------------------------
51 //ENDHEADER
52 
53 #include <iostream> // needed for io
54 #include <sstream> // needed for internal io
55 #include <iomanip>
56 #include <cmath>
57 
58 #include <fastjet/ClusterSequence.hh>
59 #include <fastjet/tools/MassDropTagger.hh>
60 #include <fastjet/tools/Filter.hh>
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 &, const PseudoJet &);
108 
109 
110 //----------------------------------------------------------------------
111 // core of the program
112 //----------------------------------------------------------------------
113 int main(){
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 jet tagging using a mass drop tagger
157  //
158  // Note: if you prefer, you may as well use a CASubJetTagger
159  // CASubJetTagger ca_tagger;
160  // PseudoJet tagged = ca_tagger(jets[0]);
161  // This requires including fastjet/tools/CASubJetTagger.hh
162  // You also need to adapt the 2 lines below accessing
163  // the extra structural information provided by the tagger
164  //----------------------------------------------------------
165  MassDropTagger md_tagger(0.667, 0.09);
166  PseudoJet tagged = md_tagger(jets[0]);
167 
168  if (tagged == 0){
169  cout << "No substructure found" << endl;
170  return 0;
171  }
172 
173  PseudoJet parent1 = tagged.pieces()[0];
174  PseudoJet parent2 = tagged.pieces()[1];
175  cout << "Found suitable pair of subjets: " << endl;
176  cout << " " << parent1 << endl;
177  cout << " " << parent2 << endl;
178  cout << "Total = " << endl;
179  cout << " " << tagged << endl;
180  cout << "(mass drop = " << tagged.structure_of<MassDropTagger>().mu()
181  << ", y = " << tagged.structure_of<MassDropTagger>().y() << ")"
182  << endl << endl;
183 
184  // next we "filter" it, to remove UE & pileup contamination
185  //----------------------------------------------------------
186  //
187  // [there are two ways of doing this; here we directly use the
188  // existing cluster sequence and find the exclusive subjets of
189  // this_jet (i.e. work backwards within the cs starting from
190  // this_jet); alternatively one can recluster just the
191  // constituents of the jet]
192  //
193  // first get separation between the subjets (called Rbb -- assuming
194  // it's a Higgs!)
195  //
196  // See example 11-filter.cc for another way of implementing the dynamic
197  // Rfilt used below
198  double Rbb = parent1.delta_R(parent2);
199  double Rfilt = min(Rbb/2, 0.3); // somewhat arbitrary choice
200  unsigned nfilt = 3; // number of pieces we'll take
201  cout << "Subjet separation (Rbb) = " << Rbb << ", Rfilt = " << Rfilt << endl;
202 
203  Filter filter(JetDefinition(cambridge_algorithm, Rfilt, &flav_recombiner),
204  SelectorNHardest(nfilt));
205  PseudoJet filtered = filter(tagged);
206 
207  // now print out the filtered jets and reconstruct total
208  // at the same time
209  const vector<PseudoJet> & filtered_pieces = filtered.pieces();
210  cout << "Filtered pieces are " << endl;
211  for (unsigned i = 0; i < nfilt && i < filtered_pieces.size(); i++) {
212  cout << " " << filtered_pieces[i] << endl;
213  }
214  cout << "Filtered total is " << endl;
215  cout << " " << filtered << endl;
216 
217 }
218 
219 
220 //----------------------------------------------------------------------
221 // does the actual work for printing out a jet
222 //----------------------------------------------------------------------
223 ostream & operator<<(ostream & ostr, const PseudoJet & jet) {
224  ostr << "pt, y, phi ="
225  << " " << setw(10) << jet.perp()
226  << " " << setw(6) << jet.rap()
227  << " " << setw(6) << jet.phi()
228  << ", mass = " << setw(10) << jet.m()
229  << ", btag = " << jet.user_index();
230  return ostr;
231 }