FastJet 3.0.2
|
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 }