fastjet_boosted_higgs.cc
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052 #include "fastjet/ClusterSequence.hh"
00053 #include<iostream>
00054 #include<sstream>
00055 #include<iomanip>
00056 #include<cmath>
00057
00058 using namespace std;
00059 namespace fj = fastjet;
00060
00061
00062
00063
00064
00073 typedef fj::JetDefinition::DefaultRecombiner DefRecomb;
00074 class FlavourRecombiner : public DefRecomb {
00075 public:
00076 FlavourRecombiner(fj::RecombinationScheme recomb_scheme = fj::E_scheme) :
00077 DefRecomb(recomb_scheme) {};
00078
00079 virtual std::string description() const {return DefRecomb::description()
00080 +" (with user index addition)";}
00081
00083 virtual void recombine(const fj::PseudoJet & pa,
00084 const fj::PseudoJet & pb,
00085 fj::PseudoJet & pab) const {
00086 DefRecomb::recombine(pa,pb,pab);
00087 pab.set_user_index(pa.user_index() + pb.user_index());
00088 }
00089
00090 };
00091
00092
00094 ostream & operator<<(ostream &, fj::PseudoJet &);
00095
00096
00097
00098 int main (int argc, char ** argv) {
00099
00100 vector<fj::PseudoJet> particles;
00101
00102
00103 string line;
00104 while (getline(cin,line)) {
00105 if (line.substr(0,1) == "#") {continue;}
00106 istringstream linestream(line);
00107 double px,py,pz,E;
00108 linestream >> px >> py >> pz >> E;
00109
00110
00111 int btag;
00112 if (! (linestream >> btag)) btag = 0;
00113
00114
00115 fj::PseudoJet particle(px,py,pz,E);
00116 particle.set_user_index(btag);
00117 particles.push_back(particle);
00118 }
00119
00120
00121
00122 double R = 1.2;
00123 FlavourRecombiner flav_recombiner;
00124 fj::JetDefinition jet_def(fj::cambridge_algorithm, R, &flav_recombiner);
00125
00126
00127
00128 fj::ClusterSequence cs(particles, jet_def);
00129 vector<fj::PseudoJet> jets = sorted_by_pt(cs.inclusive_jets());
00130
00131
00132 cout << "Ran: " << jet_def.description() << endl << endl;
00133 cout << "Hardest jet: " << jets[0] << endl << endl;
00134
00135
00137
00142 double mass_drop_threshold = 0.667;
00151 double rtycut = 0.3;
00152
00153 fj::PseudoJet this_jet = jets[0], parent1, parent2;
00154 bool had_parents;
00155
00156 while ((had_parents = cs.has_parents(this_jet,parent1,parent2))) {
00157
00158 if (parent1.m() < parent2.m()) swap(parent1,parent2);
00159
00160
00161
00162
00163 if (parent1.m() < mass_drop_threshold * this_jet.m() &&
00164 parent1.kt_distance(parent2) > pow(rtycut,2) * this_jet.m2()) {
00165 break;
00166 } else {
00167
00168 this_jet = parent1;
00169 }
00170 }
00171
00172
00173 if (had_parents) {
00174 cout << "Found suitable pair of subjets: " << endl;
00175 cout << " " << parent1 << endl;
00176 cout << " " << parent2 << endl;
00177 cout << "Total = " << endl;
00178 cout << " " << this_jet << endl << endl;
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189 double Rbb = sqrt(parent1.squared_distance(parent2));
00190 double Rfilt = min(Rbb/2, 0.3);
00191 unsigned nfilt = 3;
00192 cout << "Subjet separation (Rbb) = " << Rbb << ", Rfilt = " << Rfilt << endl;
00193
00194 double dcut = pow(Rfilt/R,2);
00195
00196 vector<fj::PseudoJet> filt_subjets = sorted_by_pt(cs.exclusive_subjets(this_jet, dcut));
00197
00198
00199
00200 cout << "Filtered pieces are " << endl;
00201 cout << " " << filt_subjets[0] << endl;
00202 fj::PseudoJet filtered_total = filt_subjets[0];
00203 for (unsigned i = 1; i < nfilt && i < filt_subjets.size(); i++) {
00204 cout << " " << filt_subjets[i] << endl;
00205 flav_recombiner.plus_equal(filtered_total, filt_subjets[i]);
00206 }
00207 cout << "Filtered total is " << endl;
00208 cout << " " << filtered_total << endl;
00209
00210 } else {
00211 cout << "Did not find suitable hard substructure in this event." << endl;
00212 }
00213 }
00214
00215
00217 ostream & operator<<(ostream & ostr, fj::PseudoJet & jet) {
00218 ostr << "pt, y, phi ="
00219 << " " << setw(10) << jet.perp()
00220 << " " << setw(6) << jet.rap()
00221 << " " << setw(6) << jet.phi()
00222 << ", mass = " << setw(10) << jet.m()
00223 << ", btag = " << jet.user_index();
00224 return ostr;
00225 }