50 #include "fastjet/ClusterSequence.hh"
57 namespace fj = fastjet;
71 typedef fj::JetDefinition::DefaultRecombiner
DefRecomb;
72 class FlavourRecombiner :
public DefRecomb {
78 +
" (with user index addition)";}
96 int main (
int argc,
char ** argv) {
98 vector<fj::PseudoJet> particles;
102 while (getline(cin,line)) {
103 if (line.substr(0,1) ==
"#") {
continue;}
104 istringstream linestream(line);
106 linestream >> px >> py >> pz >> E;
110 if (! (linestream >> btag)) btag = 0;
114 particle.set_user_index(btag);
115 particles.push_back(particle);
121 FlavourRecombiner flav_recombiner;
127 vector<fj::PseudoJet> jets =
sorted_by_pt(cs.inclusive_jets());
130 cout <<
"Ran: " << jet_def.
description() << endl << endl;
131 cout <<
"Hardest jet: " << jets[0] << endl << endl;
140 double mass_drop_threshold = 0.667;
154 while ((had_parents = this_jet.
has_parents(parent1,parent2))) {
156 if (parent1.m() < parent2.m())
swap(parent1,parent2);
161 if (parent1.m() < mass_drop_threshold * this_jet.
m() &&
162 parent1.kt_distance(parent2) > pow(rtycut,2) * this_jet.
m2()) {
172 cout <<
"Found suitable pair of subjets: " << endl;
173 cout <<
" " << parent1 << endl;
174 cout <<
" " << parent2 << endl;
175 cout <<
"Total = " << endl;
176 cout <<
" " << this_jet << endl << endl;
187 double Rbb = sqrt(parent1.squared_distance(parent2));
188 double Rfilt = min(Rbb/2, 0.3);
190 cout <<
"Subjet separation (Rbb) = " << Rbb <<
", Rfilt = " << Rfilt << endl;
192 double dcut = pow(Rfilt/R,2);
198 cout <<
"Filtered pieces are " << endl;
199 cout <<
" " << filt_subjets[0] << endl;
201 for (
unsigned i = 1; i < nfilt && i < filt_subjets.size(); i++) {
202 cout <<
" " << filt_subjets[i] << endl;
203 flav_recombiner.plus_equal(filtered_total, filt_subjets[i]);
205 cout <<
"Filtered total is " << endl;
206 cout <<
" " << filtered_total << endl;
209 cout <<
"Did not find suitable hard substructure in this event." << endl;
216 ostr <<
"pt, y, phi ="
217 <<
" " << setw(10) << jet.
perp()
218 <<
" " << setw(6) << jet.
rap()
219 <<
" " << setw(6) << jet.
phi()
220 <<
", mass = " << setw(10) << jet.
m()