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