FastJet 3.0.2
12-boosted_higgs-old.cc
Go to the documentation of this file.
00001 //----------------------------------------------------------------------
00002 /// \file
00003 /// \page Example12old 12 - boosted Higgs tagging (old version)
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 ///  - following the decomposition of a jet into pieces
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 /// Note that this example is deprecated --- see 12-boosted_higgs.cc
00020 /// for the newest version --- so it is not built by default
00021 ///
00022 /// run it with    : ./12-boosted_higgs-old < data/HZ-event-Hmass115.dat
00023 ///
00024 /// Source code: 12-boosted_higgs-old.cc
00025 //----------------------------------------------------------------------
00026 
00027 
00028 //STARTHEADER
00029 // $Id: 12-boosted_higgs-old.cc 2577 2011-09-13 15:11:38Z salam $
00030 //
00031 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00032 //
00033 //----------------------------------------------------------------------
00034 // This file is part of FastJet.
00035 //
00036 //  FastJet is free software; you can redistribute it and/or modify
00037 //  it under the terms of the GNU General Public License as published by
00038 //  the Free Software Foundation; either version 2 of the License, or
00039 //  (at your option) any later version.
00040 //
00041 //  The algorithms that underlie FastJet have required considerable
00042 //  development and are described in hep-ph/0512210. If you use
00043 //  FastJet as part of work towards a scientific publication, please
00044 //  include a citation to the FastJet paper.
00045 //
00046 //  FastJet is distributed in the hope that it will be useful,
00047 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00048 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00049 //  GNU General Public License for more details.
00050 //
00051 //  You should have received a copy of the GNU General Public License
00052 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00053 //----------------------------------------------------------------------
00054 //ENDHEADER
00055 
00056 #include "fastjet/ClusterSequence.hh"
00057 #include <iostream> // needed for io
00058 #include <sstream>  // needed for internal io
00059 #include <iomanip>  
00060 #include <cmath>
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 &, PseudoJet &);
00108 
00109 
00110 //----------------------------------------------------------------------
00111 // core of the program
00112 //----------------------------------------------------------------------
00113 int main (int argc, char ** argv) {
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 the subjet decomposition
00157   //----------------------------------------------------------
00158   //
00159   // when unpeeling a C/A jet, often only a very soft piece may break off;
00160   // the mass_drop_threshold indicates how much "lighter" the heavier of the two
00161   // resulting pieces must be in order for us to consider that we've really
00162   // seen some form of substructure
00163   double mass_drop_threshold = 0.667; 
00164   // QCD backgrounds that give larger jet masses have a component
00165   // where a quite soft gluon is emitted; to eliminate part of this
00166   // one can place a cut on the asymmetry of the branching; 
00167   //
00168   // Here the cut is expressed in terms of y, the kt-distance scaled
00169   // to the squared jet mass; an easier way to see it is in terms of
00170   // a requirement on the momentum fraction in the splitting: z/(1-z)
00171   // and (1-z)/z > rtycut^2 [the correspondence holds only at LO]
00172   double rtycut              = 0.3;
00173 
00174   PseudoJet this_jet = jets[0];
00175   PseudoJet parent1, parent2;
00176   bool had_parents;
00177 
00178   while ((had_parents = this_jet.has_parents(parent1,parent2))) {
00179     // make parent1 the more massive jet
00180     if (parent1.m() < parent2.m()) swap(parent1,parent2);
00181 
00182     // if we pass the conditions on the mass drop and its degree of
00183     // asymmetry (z/(1-z) \sim kt_dist/m^2 > rtycut), then we've found
00184     // something interesting, so exit the loop
00185     if (parent1.m() < mass_drop_threshold * this_jet.m() &&
00186         parent1.kt_distance(parent2) > pow(rtycut,2) * this_jet.m2()) {
00187       break;
00188     } else {
00189       // otherwise try a futher decomposition on the more massive jet
00190       this_jet = parent1;
00191     }
00192   }
00193 
00194   // look to see what we found
00195   if (!had_parents) {
00196     cout << "Did not find suitable hard substructure in this event." << endl;
00197     return 0;
00198   }
00199 
00200   cout << "Found suitable pair of subjets: " << endl;
00201   cout << " " << parent1 << endl;
00202   cout << " " << parent2 << endl;
00203   cout << "Total = " << endl;
00204   cout << " " << this_jet << endl << endl;
00205 
00206   // next we "filter" it, to remove UE & pileup contamination
00207   //----------------------------------------------------------
00208   //
00209   // [there are two ways of doing this; here we directly use the
00210   // exsiting cluster sequence and find the exclusive subjets of
00211   // this_jet (i.e. work backwards within the cs starting from
00212   // this_jet); alternatively one can recluster just the
00213   // constituents of the jet]
00214   //
00215   // first get separation between the subjets (called Rbb -- assuming it's a Higgs!)
00216   double   Rbb = sqrt(parent1.squared_distance(parent2));
00217   double   Rfilt = min(Rbb/2, 0.3); // somewhat arbitrary choice
00218   unsigned nfilt = 3;               // number of pieces we'll take
00219   cout << "Subjet separation (Rbb) = " << Rbb << ", Rfilt = " << Rfilt << endl;
00220 
00221   double   dcut  = pow(Rfilt/R,2);  // for C/A get a view at Rfilt by
00222   // using a dcut=(Rfilt/R)^2
00223   vector<PseudoJet> filt_subjets = sorted_by_pt(this_jet.exclusive_subjets(dcut));
00224 
00225   // now print out the filtered jets and reconstruct total 
00226   // at the same time
00227   cout << "Filtered pieces are " << endl;
00228   cout << " " << filt_subjets[0] << endl;
00229   PseudoJet filtered_total = filt_subjets[0];
00230   for (unsigned i = 1; i < nfilt && i < filt_subjets.size(); i++) {
00231     cout << " " << filt_subjets[i] << endl;
00232     flav_recombiner.plus_equal(filtered_total, filt_subjets[i]);
00233   }
00234   cout << "Filtered total is " << endl;
00235   cout << " " << filtered_total << endl;
00236 
00237 }
00238 
00239 
00240 //----------------------------------------------------------------------
00241 // does the actual work for printing out a jet
00242 //----------------------------------------------------------------------
00243 ostream & operator<<(ostream & ostr, PseudoJet & jet) {
00244   ostr << "pt, y, phi =" 
00245        << " " << setw(10) << jet.perp() 
00246        << " " << setw(6) <<  jet.rap()  
00247        << " " << setw(6) <<  jet.phi()  
00248        << ", mass = " << setw(10) << jet.m()
00249        << ", btag = " << jet.user_index();
00250   return ostr;
00251 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends