52#include "fastjet/ClusterSequence.hh"
73typedef fj::JetDefinition::DefaultRecombiner
DefRecomb;
74class FlavourRecombiner :
public DefRecomb {
80 +
" (with user index addition)";}
98int main (
int argc,
char ** argv) {
100 vector<fj::PseudoJet> particles;
104 while (getline(cin,line)) {
105 if (line.substr(0,1) ==
"#") {
continue;}
106 istringstream linestream(line);
108 linestream >> px >> py >> pz >> E;
112 if (! (linestream >> btag)) btag = 0;
116 particle.set_user_index(btag);
117 particles.push_back(particle);
123 FlavourRecombiner flav_recombiner;
129 vector<fj::PseudoJet> jets =
sorted_by_pt(cs.inclusive_jets());
132 cout <<
"Ran: " << jet_def.description() << endl << endl;
133 cout <<
"Hardest jet: " << jets[0] << endl << endl;
142 double mass_drop_threshold = 0.667;
156 while ((had_parents = this_jet.
has_parents(parent1,parent2))) {
158 if (parent1.m() < parent2.m())
swap(parent1,parent2);
163 if (parent1.m() < mass_drop_threshold * this_jet.
m() &&
164 parent1.kt_distance(parent2) > pow(rtycut,2) * this_jet.
m2()) {
174 cout <<
"Found suitable pair of subjets: " << endl;
175 cout <<
" " << parent1 << endl;
176 cout <<
" " << parent2 << endl;
177 cout <<
"Total = " << endl;
178 cout <<
" " << this_jet << endl << endl;
189 double Rbb = sqrt(parent1.squared_distance(parent2));
190 double Rfilt = min(Rbb/2, 0.3);
192 cout <<
"Subjet separation (Rbb) = " << Rbb <<
", Rfilt = " << Rfilt << endl;
194 double dcut = pow(Rfilt/R,2);
200 cout <<
"Filtered pieces are " << endl;
201 cout <<
" " << filt_subjets[0] << endl;
203 for (
unsigned i = 1; i < nfilt && i < filt_subjets.size(); i++) {
204 cout <<
" " << filt_subjets[i] << endl;
205 flav_recombiner.plus_equal(filtered_total, filt_subjets[i]);
207 cout <<
"Filtered total is " << endl;
208 cout <<
" " << filtered_total << endl;
211 cout <<
"Did not find suitable hard substructure in this event." << endl;
218 ostr <<
"pt, y, phi ="
219 <<
" " << setw(10) << jet.
perp()
220 <<
" " << setw(6) << jet.
rap()
221 <<
" " << setw(6) << jet.
phi()
222 <<
", mass = " << setw(10) << jet.
m()
int main()
an example program showing how to use fastjet
A class that will provide the recombination scheme facilities and/or allow a user to extend these fac...
virtual void recombine(const PseudoJet &pa, const PseudoJet &pb, PseudoJet &pab) const override
recombine pa and pb and put result into pab
virtual std::string description() const override
return a textual description of the recombination scheme implemented here
class that is intended to hold a full definition of the jet clusterer
Class to contain pseudojets, including minimal information of use to jet-clustering routines.
double rap() const
returns the rapidity or some large value when the rapidity is infinite
void set_user_index(const int index)
set the user_index, intended to allow the user to add simple identifying information to a particle/je...
double phi() const
returns phi (in the range 0..2pi)
double perp() const
returns the scalar transverse momentum
double m2() const
returns the squared invariant mass // like CLHEP
virtual bool has_parents(PseudoJet &parent1, PseudoJet &parent2) const
check if it is the product of a recombination, in which case return the 2 parents through the 'parent...
std::vector< PseudoJet > exclusive_subjets(const double dcut) const
return a vector of all subjets of the current jet (in the sense of the exclusive algorithm) that woul...
double m() const
returns the invariant mass (If m2() is negative then -sqrt(-m2()) is returned, as in CLHEP)
int user_index() const
return the user_index,
RecombinationScheme
The various recombination schemes.
@ E_scheme
summing the 4-momenta
@ cambridge_algorithm
the longitudinally invariant variant of the cambridge algorithm (aka Aachen algoithm).
void swap(SharedPtr< T > &a, SharedPtr< T > &b)
swapping
vector< PseudoJet > sorted_by_pt(const vector< PseudoJet > &jets)
return a vector of jets sorted into decreasing kt2