FastJet 3.0.4
fastjet_boosted_higgs.cc
00001 //STARTHEADER
00002 // $Id: fastjet_boosted_higgs.cc 2577 2011-09-13 15:11:38Z salam $
00003 //
00004 // Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
00005 //
00006 //----------------------------------------------------------------------
00007 // This file is part of FastJet.
00008 //
00009 //  FastJet is free software; you can redistribute it and/or modify
00010 //  it under the terms of the GNU General Public License as published by
00011 //  the Free Software Foundation; either version 2 of the License, or
00012 //  (at your option) any later version.
00013 //
00014 //  The algorithms that underlie FastJet have required considerable
00015 //  development and are described in hep-ph/0512210. If you use
00016 //  FastJet as part of work towards a scientific publication, please
00017 //  include a citation to the FastJet paper.
00018 //
00019 //  FastJet is distributed in the hope that it will be useful,
00020 //  but WITHOUT ANY WARRANTY; without even the implied warranty of
00021 //  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00022 //  GNU General Public License for more details.
00023 //
00024 //  You should have received a copy of the GNU General Public License
00025 //  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
00026 //----------------------------------------------------------------------
00027 //ENDHEADER
00028 
00029 
00030 //----------------------------------------------------------------------
00031 // fastjet example program, illustration of carrying out boosted
00032 // Higgs subjet ID analysis
00033 //
00034 // It illustrates two kinds of functionality: 
00035 //
00036 //  - following the decomposition of a jet into pieces
00037 //  - following information on a b-tag through the jet
00038 //
00039 // This kind of functionality was used in arXiv:0802.2470
00040 // (Butterworth, Davison, Rubin & Salam) for boosted Higgs searches,
00041 // and related functionality was used in arXiv:0806.0848 (Kaplan,
00042 // Rehermann, Schwartz & Tweedie) in searching for boosted tops
00043 // (without b-tag assumptions).
00044 
00045 // Compile it with: make fastjet_boosted_higgs
00046 // run it with    : ./fastjet_boosted_higgs < data/HZ-event-Hmass115.dat
00047 //
00048 //----------------------------------------------------------------------
00049 
00050 #include "fastjet/ClusterSequence.hh"
00051 #include<iostream> // needed for io
00052 #include<sstream>  // needed for internal io
00053 #include<iomanip>  
00054 #include<cmath>
00055 
00056 using namespace std;
00057 namespace fj = fastjet;
00058 
00059 
00060 
00061 
00062 //----------------------------------------------------------------------
00063 /// set up a class to give standard (by default E-scheme)
00064 /// recombination, with additional tracking of flavour information in
00065 /// the user_index. 
00066 ///
00067 /// If you use this, you must explicitly set the user index to 0 for
00068 /// non-flavoured particles (the default value is -1);
00069 ///
00070 /// This will work for native algorithms, but not for all plugins
00071 typedef fj::JetDefinition::DefaultRecombiner DefRecomb;
00072 class FlavourRecombiner : public  DefRecomb {
00073 public:
00074   FlavourRecombiner(fj::RecombinationScheme recomb_scheme = fj::E_scheme) : 
00075     DefRecomb(recomb_scheme) {};
00076 
00077   virtual std::string description() const {return DefRecomb::description()
00078       +" (with user index addition)";}
00079 
00080   /// recombine pa and pb and put result into pab
00081   virtual void recombine(const fj::PseudoJet & pa, 
00082                          const fj::PseudoJet & pb, 
00083                          fj::PseudoJet & pab) const {
00084     DefRecomb::recombine(pa,pb,pab);
00085     pab.set_user_index(pa.user_index() + pb.user_index());
00086   }
00087 
00088 };
00089 
00090 
00091 /// forward declaration for printing out info about a jet
00092 ostream & operator<<(ostream &, fj::PseudoJet &);
00093 
00094 
00095 //----------------------------------------------------------------------
00096 int main (int argc, char ** argv) {
00097   
00098   vector<fj::PseudoJet> particles;
00099 
00100   // read in data in format px py pz E b-tag [last of these is optional]
00101   string line;
00102   while (getline(cin,line)) {
00103     if (line.substr(0,1) == "#") {continue;}
00104     istringstream linestream(line);
00105     double px,py,pz,E;
00106     linestream >> px >> py >> pz >> E;
00107 
00108     // optionally read in btag information
00109     int    btag;
00110     if (! (linestream >> btag)) btag = 0;
00111 
00112     // construct the particle
00113     fj::PseudoJet particle(px,py,pz,E);
00114     particle.set_user_index(btag); // btag info goes in user index, for flavour tracking
00115     particles.push_back(particle);
00116   }
00117 
00118 
00119   // set up the jet finding
00120   double R = 1.2;
00121   FlavourRecombiner flav_recombiner; // for tracking flavour
00122   fj::JetDefinition jet_def(fj::cambridge_algorithm, R, &flav_recombiner);
00123 
00124   
00125   // run the jet finding; find the hardest jet
00126   fj::ClusterSequence cs(particles, jet_def);
00127   vector<fj::PseudoJet> jets = sorted_by_pt(cs.inclusive_jets());
00128 
00129 
00130   cout << "Ran: " << jet_def.description() << endl << endl;
00131   cout << "Hardest jet: " << jets[0] << endl << endl;
00132 
00133 
00134   /// now do the subjet decomposition;
00135   //
00136   /// when unpeeling a C/A jet, often only a very soft piece may break off;
00137   /// the mass_drop_threshold indicates how much "lighter" the heavier of the two
00138   /// resulting pieces must be in order for us to consider that we've really
00139   /// seen some form of substructure
00140   double mass_drop_threshold = 0.667; 
00141   /// QCD backgrounds that give larger jet masses have a component
00142   /// where a quite soft gluon is emitted; to eliminate part of this
00143   /// one can place a cut on the asymmetry of the branching; 
00144   ///
00145   /// Here the cut is expressed in terms of y, the kt-distance scaled
00146   /// to the squared jet mass; an easier way to see it is in terms of
00147   /// a requirement on the momentum fraction in the splitting: z/(1-z)
00148   /// and (1-z)/z > rtycut^2 [the correspondence holds only at LO]
00149   double rtycut              = 0.3;
00150 
00151   fj::PseudoJet this_jet = jets[0], parent1, parent2;
00152   bool had_parents;
00153 
00154   while ((had_parents = this_jet.has_parents(parent1,parent2))) {
00155     // make parent1 the more massive jet
00156     if (parent1.m() < parent2.m()) swap(parent1,parent2);
00157     //
00158     // if we pass the conditions on the mass drop and its degree of
00159     // asymmetry (z/(1-z) \sim kt_dist/m^2 > rtycut), then we've found
00160     // something interesting, so exit the loop
00161     if (parent1.m() < mass_drop_threshold * this_jet.m() &&
00162         parent1.kt_distance(parent2) > pow(rtycut,2) * this_jet.m2()) {
00163       break;
00164     } else {
00165       // otherwise try a futher decomposition on the more massive jet
00166       this_jet = parent1;
00167     }
00168   }
00169 
00170   // look to see what we found
00171   if (had_parents) {
00172     cout << "Found suitable pair of subjets: " << endl;
00173     cout << " " << parent1 << endl;
00174     cout << " " << parent2 << endl;
00175     cout << "Total = " << endl;
00176     cout << " " << this_jet << endl << endl;
00177 
00178     // next we "filter" it, to remove UE & pileup contamination
00179     //
00180     // [there are two ways of doing this; here we directly use the
00181     // exsiting cluster sequence and find the exclusive subjets of
00182     // this_jet (i.e. work backwards within the cs starting from
00183     // this_jet); alternatively one can recluster just the
00184     // constituents of the jet]
00185     //
00186     // first get separation between the subjets (called Rbb -- assuming it's a Higgs!)
00187     double   Rbb = sqrt(parent1.squared_distance(parent2));
00188     double   Rfilt = min(Rbb/2, 0.3); // somewhat arbitrary choice
00189     unsigned nfilt = 3;               // number of pieces we'll take
00190     cout << "Subjet separation (Rbb) = " << Rbb << ", Rfilt = " << Rfilt << endl;
00191 
00192     double   dcut  = pow(Rfilt/R,2);  // for C/A get a view at Rfilt by
00193                                     // using a dcut=(Rfilt/R)^2
00194     vector<fj::PseudoJet> filt_subjets = sorted_by_pt(this_jet.exclusive_subjets(dcut));
00195 
00196     // now print out the filtered jets and reconstruct total 
00197     // at the same time
00198     cout << "Filtered pieces are " << endl;
00199     cout << " " << filt_subjets[0] << endl;
00200     fj::PseudoJet filtered_total = filt_subjets[0];
00201     for (unsigned i = 1; i < nfilt && i < filt_subjets.size(); i++) {
00202       cout << " " << filt_subjets[i] << endl;
00203       flav_recombiner.plus_equal(filtered_total, filt_subjets[i]);
00204     }
00205     cout << "Filtered total is " << endl;
00206     cout << " " << filtered_total << endl;
00207 
00208   } else {
00209     cout << "Did not find suitable hard substructure in this event." << endl;
00210   }
00211 }
00212 
00213 
00214 /// does the actual work for printing out a jet
00215 ostream & operator<<(ostream & ostr, fj::PseudoJet & jet) {
00216   ostr << "pt, y, phi =" 
00217        << " " << setw(10) << jet.perp() 
00218        << " " << setw(6) <<  jet.rap()  
00219        << " " << setw(6) <<  jet.phi()  
00220        << ", mass = " << setw(10) << jet.m()
00221        << ", btag = " << jet.user_index();
00222   return ostr;
00223 }
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends