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