56 #include "fastjet/ClusterSequence.hh"
63 using namespace fastjet;
85 class FlavourRecombiner :
public DefRecomb {
90 virtual std::string description()
const {
113 int main (
int argc,
char ** argv) {
115 vector<PseudoJet> particles;
122 while (getline(cin,line)) {
123 if (line.substr(0,1) ==
"#") {
continue;}
124 istringstream linestream(line);
126 linestream >> px >> py >> pz >> E;
130 if (! (linestream >> btag)) btag = 0;
134 particle.set_user_index(btag);
135 particles.push_back(particle);
145 FlavourRecombiner flav_recombiner;
146 JetDefinition jet_def(cambridge_algorithm, R, &flav_recombiner);
151 vector<PseudoJet> jets =
sorted_by_pt(cs.inclusive_jets());
153 cout <<
"Ran: " << jet_def.
description() << endl << endl;
154 cout <<
"Hardest jet: " << jets[0] << endl << endl;
163 double mass_drop_threshold = 0.667;
178 while ((had_parents = this_jet.
has_parents(parent1,parent2))) {
180 if (parent1.
m() < parent2.
m())
swap(parent1,parent2);
185 if (parent1.
m() < mass_drop_threshold * this_jet.
m() &&
186 parent1.
kt_distance(parent2) > pow(rtycut,2) * this_jet.
m2()) {
196 cout <<
"Did not find suitable hard substructure in this event." << endl;
200 cout <<
"Found suitable pair of subjets: " << endl;
201 cout <<
" " << parent1 << endl;
202 cout <<
" " << parent2 << endl;
203 cout <<
"Total = " << endl;
204 cout <<
" " << this_jet << endl << endl;
217 double Rfilt = min(Rbb/2, 0.3);
219 cout <<
"Subjet separation (Rbb) = " << Rbb <<
", Rfilt = " << Rfilt << endl;
221 double dcut = pow(Rfilt/R,2);
227 cout <<
"Filtered pieces are " << endl;
228 cout <<
" " << filt_subjets[0] << endl;
229 PseudoJet filtered_total = filt_subjets[0];
230 for (
unsigned i = 1; i < nfilt && i < filt_subjets.size(); i++) {
231 cout <<
" " << filt_subjets[i] << endl;
232 flav_recombiner.plus_equal(filtered_total, filt_subjets[i]);
234 cout <<
"Filtered total is " << endl;
235 cout <<
" " << filtered_total << endl;
244 ostr <<
"pt, y, phi ="
245 <<
" " << setw(10) << jet.
perp()
246 <<
" " << setw(6) << jet.
rap()
247 <<
" " << setw(6) << jet.
phi()
248 <<
", mass = " << setw(10) << jet.
m()