FastJet 3.0.2
12-boosted_higgs.cc
Go to the documentation of this file.
00001 //----------------------------------------------------------------------
00002 /// \file
00003 /// \page Example12 12 - boosted Higgs tagging
00004 ///
00005 /// fastjet example program, illustration of carrying out boosted
00006 /// Higgs subjet ID analysis
00007 ///
00008 /// It illustrates two kinds of functionality: 
00009 ///
00010 ///  - using a boosted higgs tagger
00011 ///  - following information on a b-tag through the jet
00012 ///
00013 /// This kind of functionality was used in arXiv:0802.2470
00014 /// (Butterworth, Davison, Rubin & Salam) for boosted Higgs searches,
00015 /// and related functionality was used in arXiv:0806.0848 (Kaplan,
00016 /// Rehermann, Schwartz & Tweedie) in searching for boosted tops
00017 /// (without b-tag assumptions).
00018 ///
00019 /// run it with    : ./12-boosted_higgs < data/HZ-event-Hmass115.dat
00020 ///
00021 /// Source code: 12-boosted_higgs.cc
00022 //----------------------------------------------------------------------
00023 
00024 
00025 //STARTHEADER
00026 // $Id: 12-boosted_higgs.cc 2684 2011-11-14 07:41:44Z soyez $
00027 //
00028 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00029 //
00030 //----------------------------------------------------------------------
00031 // This file is part of FastJet.
00032 //
00033 //  FastJet is free software; you can redistribute it and/or modify
00034 //  it under the terms of the GNU General Public License as published by
00035 //  the Free Software Foundation; either version 2 of the License, or
00036 //  (at your option) any later version.
00037 //
00038 //  The algorithms that underlie FastJet have required considerable
00039 //  development and are described in hep-ph/0512210. If you use
00040 //  FastJet as part of work towards a scientific publication, please
00041 //  include a citation to the FastJet paper.
00042 //
00043 //  FastJet is distributed in the hope that it will be useful,
00044 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00045 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00046 //  GNU General Public License for more details.
00047 //
00048 //  You should have received a copy of the GNU General Public License
00049 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00050 //----------------------------------------------------------------------
00051 //ENDHEADER
00052 
00053 #include <iostream> // needed for io
00054 #include <sstream>  // needed for internal io
00055 #include <iomanip>  
00056 #include <cmath>
00057 
00058 #include <fastjet/ClusterSequence.hh>
00059 #include <fastjet/tools/MassDropTagger.hh>
00060 #include <fastjet/tools/Filter.hh>
00061 
00062 using namespace std;
00063 using namespace fastjet;
00064 
00065 
00066 //----------------------------------------------------------------------
00067 // set up a class to give standard (by default E-scheme)
00068 // recombination, with additional tracking of flavour information in
00069 // the user_index. 
00070 //
00071 // b-tagged particles are assumed to have their user_index set to 1,
00072 // and other particles should have user_index to 0.
00073 //
00074 // Watch out however that, by default, the user_index of a particle is
00075 // set to -1 and you may not have control over that (e.g. if you
00076 // compute the jet area using explicit ghosts, the ghosts will have a
00077 // default user_index of -1). For that reason, if one of the particle
00078 // being combined has a user index of -1, we assume it is not b-tagged
00079 // (i.e. we count it as 0 in the recombination)
00080 //
00081 // This will work for native algorithms, but not for all plugins
00082 //----------------------------------------------------------------------
00083 typedef JetDefinition::DefaultRecombiner DefRecomb;
00084 
00085 class FlavourRecombiner : public  DefRecomb {
00086 public:
00087   FlavourRecombiner(RecombinationScheme recomb_scheme = E_scheme) : 
00088     DefRecomb(recomb_scheme) {};
00089 
00090   virtual std::string description() const {
00091     return DefRecomb::description()+" (with user index addition)";}
00092 
00093   /// recombine pa and pb and put result into pab
00094   virtual void recombine(const PseudoJet & pa, const PseudoJet & pb, 
00095                          PseudoJet & pab) const {
00096     DefRecomb::recombine(pa,pb,pab);
00097     // Note: see the above discussion for the fact that we consider
00098     // negative user indices as "0"
00099     pab.set_user_index(max(pa.user_index(),0) + max(pb.user_index(),0));
00100   }
00101 };
00102 
00103 
00104 //----------------------------------------------------------------------
00105 // forward declaration for printing out info about a jet
00106 //----------------------------------------------------------------------
00107 ostream & operator<<(ostream &, const PseudoJet &);
00108 
00109 
00110 //----------------------------------------------------------------------
00111 // core of the program
00112 //----------------------------------------------------------------------
00113 int main(){
00114 
00115   vector<PseudoJet> particles;
00116 
00117   // read in data in format px py pz E b-tag [last of these is optional]
00118   // lines starting with "#" are considered as comments and discarded
00119   //----------------------------------------------------------
00120 
00121   string line;
00122   while (getline(cin,line)) {
00123     if (line.substr(0,1) == "#") {continue;}
00124     istringstream linestream(line);
00125     double px,py,pz,E;
00126     linestream >> px >> py >> pz >> E;
00127 
00128     // optionally read in btag information
00129     int    btag;
00130     if (! (linestream >> btag)) btag = 0;
00131 
00132     // construct the particle
00133     PseudoJet particle(px,py,pz,E);
00134     particle.set_user_index(btag); // btag info goes in user index, for flavour tracking
00135     particles.push_back(particle);
00136   }
00137 
00138 
00139   // set up the jet finding
00140   //
00141   // This also shows how to use the "FlavourRecombiner" user-defined
00142   // recombiner 
00143   // ----------------------------------------------------------
00144   double R = 1.2;
00145   FlavourRecombiner flav_recombiner; // for tracking flavour
00146   JetDefinition jet_def(cambridge_algorithm, R, &flav_recombiner);
00147 
00148   
00149   // run the jet finding; find the hardest jet
00150   ClusterSequence cs(particles, jet_def);
00151   vector<PseudoJet> jets = sorted_by_pt(cs.inclusive_jets());
00152 
00153   cout << "Ran: " << jet_def.description() << endl << endl;
00154   cout << "Hardest jet: " << jets[0] << endl << endl;
00155 
00156   // now do jet tagging using a mass drop tagger
00157   //
00158   // Note: if you prefer, you may as well use a CASubJetTagger
00159   //    CASubJetTagger ca_tagger;
00160   //    PseudoJet tagged = ca_tagger(jets[0]);
00161   // This requires including fastjet/tools/CASubJetTagger.hh 
00162   // You also need to adapt the 2 lines below accessing
00163   // the extra structural information provided by the tagger
00164   //----------------------------------------------------------
00165   MassDropTagger md_tagger(0.667, 0.09);
00166   PseudoJet tagged = md_tagger(jets[0]);
00167 
00168   if (tagged == 0){
00169     cout << "No substructure found" << endl;
00170     return 0;
00171   }
00172 
00173   PseudoJet parent1 = tagged.pieces()[0];
00174   PseudoJet parent2 = tagged.pieces()[1];
00175   cout << "Found suitable pair of subjets: " << endl;
00176   cout << " " << parent1 << endl;
00177   cout << " " << parent2 << endl;
00178   cout << "Total = " << endl;
00179   cout << " " << tagged << endl;
00180   cout << "(mass drop = " << tagged.structure_of<MassDropTagger>().mu()
00181        << ", y = " << tagged.structure_of<MassDropTagger>().y() << ")" 
00182        << endl << endl;
00183 
00184   // next we "filter" it, to remove UE & pileup contamination
00185   //----------------------------------------------------------
00186   //
00187   // [there are two ways of doing this; here we directly use the
00188   // existing cluster sequence and find the exclusive subjets of
00189   // this_jet (i.e. work backwards within the cs starting from
00190   // this_jet); alternatively one can recluster just the
00191   // constituents of the jet]
00192   //
00193   // first get separation between the subjets (called Rbb -- assuming
00194   // it's a Higgs!)
00195   //
00196   // See example 11-filter.cc for another way of implementing the dynamic
00197   // Rfilt used below
00198   double   Rbb = parent1.delta_R(parent2);
00199   double   Rfilt = min(Rbb/2, 0.3); // somewhat arbitrary choice
00200   unsigned nfilt = 3;               // number of pieces we'll take
00201   cout << "Subjet separation (Rbb) = " << Rbb << ", Rfilt = " << Rfilt << endl;
00202 
00203   Filter filter(JetDefinition(cambridge_algorithm, Rfilt, &flav_recombiner),
00204                 SelectorNHardest(nfilt));
00205   PseudoJet filtered = filter(tagged);
00206 
00207   // now print out the filtered jets and reconstruct total 
00208   // at the same time
00209   const vector<PseudoJet> & filtered_pieces = filtered.pieces();
00210   cout << "Filtered pieces are " << endl;
00211   for (unsigned i = 0; i < nfilt && i < filtered_pieces.size(); i++) {
00212     cout << " " << filtered_pieces[i] << endl;
00213   }
00214   cout << "Filtered total is " << endl;
00215   cout << " " << filtered << endl;
00216 
00217 }
00218 
00219 
00220 //----------------------------------------------------------------------
00221 // does the actual work for printing out a jet
00222 //----------------------------------------------------------------------
00223 ostream & operator<<(ostream & ostr, const PseudoJet & jet) {
00224   ostr << "pt, y, phi =" 
00225        << " " << setw(10) << jet.perp() 
00226        << " " << setw(6) <<  jet.rap()  
00227        << " " << setw(6) <<  jet.phi()  
00228        << ", mass = " << setw(10) << jet.m()
00229        << ", btag = " << jet.user_index();
00230   return ostr;
00231 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends