test.cc

#include <iostream>
#include <fastjet/ClusterSequence.hh>
#include <fastjet/tools/Recluster.hh>

using namespace std;
using namespace fastjet;

//----------------------------------------------------------------------
// set up a class to give standard (by default E-scheme)
// recombination, with additional tracking of flavour information in
// the user_index. 
//
// b-tagged particles are assumed to have their user_index set to 1,
// and other particles should have user_index to 0.
//
// Watch out however that, by default, the user_index of a particle is
// set to -1 and you may not have control over that (e.g. if you
// compute the jet area using explicit ghosts, the ghosts will have a
// default user_index of -1). For that reason, if one of the particle
// being combined has a user index of -1, we assume it is not b-tagged
// (i.e. we count it as 0 in the recombination)
//
// This will work for native algorithms, but not for all plugins
//----------------------------------------------------------------------
typedef JetDefinition::DefaultRecombiner DefRecomb;

class FlavourRecombiner : public  DefRecomb {
public:
  FlavourRecombiner(RecombinationScheme recomb_scheme = E_scheme) : 
    DefRecomb(recomb_scheme) {};

  virtual std::string description() const {
    return DefRecomb::description()+" (with user index addition)";}

  /// recombine pa and pb and put result into pab
  virtual void recombine(const PseudoJet & pa, const PseudoJet & pb, 
                         PseudoJet & pab) const {
    DefRecomb::recombine(pa,pb,pab);
    // Note: see the above discussion for the fact that we consider
    // negative user indices as "0"
    pab.set_user_index(max(pa.user_index(),0) + max(pb.user_index(),0));
  }
};

//----------------------------------------------------------------------
// forward declaration for printing out info about a jet
//----------------------------------------------------------------------
ostream & operator<<(ostream &, const PseudoJet &);

//----------------------------------------------------------------------
// get the hardest jet
//----------------------------------------------------------------------
PseudoJet hardest_jet(vector<PseudoJet> &particles){
  double R = 1.2;
  FlavourRecombiner *flav_recombiner=new FlavourRecombiner; // for tracking flavour
  JetDefinition jet_def(cambridge_algorithm, R, flav_recombiner);
  jet_def.delete_recombiner_when_unused();

  ClusterSequence *cs = new ClusterSequence(particles, jet_def);
  vector<PseudoJet> jets = sorted_by_pt(cs->inclusive_jets());
  cs->delete_self_when_unused();

  cout << "Ran: " << jet_def.description() << endl << endl;
  cout << "Hardest jet: " << jets[0] << endl << endl;
  
  return jets[0];
}

//----------------------------------------------------------------------
// core of the program
//----------------------------------------------------------------------
int main(){
  vector<PseudoJet> particles;

  // read in data in format px py pz E b-tag [last of these is optional]
  // lines starting with "#" are considered as comments and discarded
  //----------------------------------------------------------

  string line;
  while (getline(cin,line)) {
    if (line.substr(0,1) == "#") {continue;}
    istringstream linestream(line);
    double px,py,pz,E;
    linestream >> px >> py >> pz >> E;

    // optionally read in btag information
    int    btag;
    if (! (linestream >> btag)) btag = 0;

    // construct the particle
    PseudoJet particle(px,py,pz,E);
    particle.set_user_index(btag); // btag info goes in user index, for flavour tracking
    particles.push_back(particle);
  }


  // set up the jet finding
  //
  // This also shows how to use the "FlavourRecombiner" user-defined
  // recombiner 
  // ----------------------------------------------------------
  PseudoJet j = hardest_jet(particles);

  Recluster f1(kt_algorithm, 0.2);
  PseudoJet fj1 = f1(j);
  cout << "Recluster1 jet: " << fj1 << endl << endl;

  //delete the original CS
  j=PseudoJet();
  
  Recluster f2(cambridge_algorithm, 1.0);
  PseudoJet fj2 = f2(fj1);
  cout << "Recluster1 jet: " << fj2 << endl << endl;

  return 0;
}

//----------------------------------------------------------------------
// does the actual work for printing out a jet
//----------------------------------------------------------------------
ostream & operator<<(ostream & ostr, const PseudoJet & jet) {
  ostr << "pt, y, phi =" 
       << " " << jet.perp() 
       << " "  <<  jet.rap()  
       << " "  <<  jet.phi()  
       << ", mass = " << jet.m()
       << ", btag = " << jet.user_index();
  return ostr;
}