fastjet_subtraction.cc File Reference

#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include <iostream>
#include <sstream>
#include <vector>

Include dependency graph for fastjet_subtraction.cc:

Go to the source code of this file.

Functions
void	print_jets (const fastjet::ClusterSequenceAreaBase &clust_seq, const vector< fastjet::PseudoJet > &unsorted_jets)
	a function that pretty prints a list of jets, and performs the subtraction in two different ways, using a generic ClusterSequenceAreaBase type object.
int	main (int argc, char **argv)
	an example program showing how to use fastjet

---

Function Documentation

int main	(	int	argc,
		char **	argv
	)

an example program showing how to use fastjet

Definition at line 61 of file fastjet_subtraction.cc.

References fastjet::Best, fastjet::ClusterSequence::inclusive_jets(), fastjet::kt_algorithm, print_jets(), fastjet::PseudoJet::rap(), fastjet::PseudoJet::set_user_index(), and fastjet::ClusterSequence::strategy_string().

                                  {
  
  // since we use here simulated data we can split the hard event
  // from the full (i.e. with pileup added) one
  vector<fastjet::PseudoJet> hard_event, full_event;
  
  // maximum rapidity that we accept for the input particles
  double etamax = 6.0;
  
  // read in input particles. Keep the hard event generated by PYTHIA
  // separated from the full event, so as to be able to gauge the
  // "goodness" of the subtraction from the full event which also
  // includes pileup
  string line;
  int  nsub  = 0;
  while (getline(cin, line)) {
    istringstream linestream(line);
    // take substrings to avoid problems when there are extra "pollution"
    // characters (e.g. line-feed).
    if (line.substr(0,4) == "#END") {break;}
    if (line.substr(0,9) == "#SUBSTART") {
      // if more sub events follow, make copy of hard one here
      if (nsub == 1) hard_event = full_event;
      nsub += 1;
    }
    if (line.substr(0,1) == "#") {continue;}
    valarray<double> fourvec(4);
    linestream >> fourvec[0] >> fourvec[1] >> fourvec[2] >> fourvec[3];
    fastjet::PseudoJet psjet(fourvec);
    psjet.set_user_index(0);
    // push event onto back of full_event vector
    if (abs(psjet.rap()) < etamax) {full_event.push_back(psjet);}
  }

  // if we have read in only one event, copy it across here...
  if (nsub == 1) hard_event = full_event;

  // if there was nothing in the event 
  if (nsub == 0) {
    cerr << "Error: read empty event\n";
    exit(-1);
  }
  
  
  // create an object that represents your choice of jet algorithm and 
  // the associated parameters
  double R = 0.7;
  fastjet::Strategy strategy = fastjet::Best;
  fastjet::JetDefinition jet_def(fastjet::kt_algorithm, R, strategy);
  //fastjet::JetDefinition jet_def(fastjet::cambridge_algorithm, R, strategy);
  //fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, R, strategy);

  // create an object that specifies how we to define the (active) area
  double ghost_etamax = 6.0;
  int    active_area_repeats = 1;
  double ghost_area    = 0.01;
  fastjet::ActiveAreaSpec area_spec(ghost_etamax, active_area_repeats, 
                                    ghost_area);
  fastjet::AreaDefinition area_def(area_spec);
  //fastjet::AreaDefinition area_def(fastjet::VoronoiAreaSpec(1.0));

  // run the jet clustering with the above jet definition. hard event first
  fastjet::ClusterSequenceArea clust_seq(hard_event, jet_def, area_def);
  //fastjet::ClusterSequencePassiveArea clust_seq(hard_event, jet_def);


  // extract the inclusive jets with pt > 5 GeV, sorted by pt
  double ptmin = 5.0;
  vector<fastjet::PseudoJet> inclusive_jets = clust_seq.inclusive_jets(ptmin);

  // print them out
  cout << "" <<endl;
  cout << "Hard event only"<<endl;
  cout << "Number of input particles: "<<hard_event.size()<<endl;
  cout << "Strategy used: "<<clust_seq.strategy_string()<<endl;
  cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV\n";
  cout << "---------------------------------------\n";
  print_jets(clust_seq, inclusive_jets);
  cout << endl;


  // repeat everything on the full event

  // run the jet clustering with the above jet definition
  //fastjet::ClusterSequenceActiveArea clust_seq_full(full_event, jet_def, area_spec);
  fastjet::ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def);
  //fastjet::ClusterSequencePassiveArea clust_seq_full(full_event, jet_def,0.9);

  // extract the inclusive jets with pt > 20 GeV, sorted by pt
  ptmin = 20.0;
  inclusive_jets = clust_seq_full.inclusive_jets(ptmin);

  // print them out
  cout << "Full event, with pileup, and its subtraction"<<endl;
  cout << "Number of input particles: "<<full_event.size()<<endl;
  cout << "Strategy used: "<<clust_seq_full.strategy_string()<<endl;
  cout << "Printing inclusive jets with pt > "<< ptmin<<" GeV (before subtraction)\n";
  cout << "---------------------------------------\n";
  print_jets(clust_seq_full, inclusive_jets);
  cout << endl;


}

void print_jets	(	const fastjet::ClusterSequenceAreaBase &	clust_seq,
		const vector< fastjet::PseudoJet > &	unsorted_jets
	)

a function that pretty prints a list of jets, and performs the subtraction in two different ways, using a generic ClusterSequenceAreaBase type object.

---