fastjet_subtraction.cc File Reference

#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequenceActiveArea.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::ClusterSequenceActiveArea &, const vector< fastjet::PseudoJet > &)
	a function that pretty prints a list of jets, and performs the subtraction in two different ways
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 59 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 finder and // the associated parameters double R = 0.7; fastjet::Strategy strategy = fastjet::Best; fastjet::JetDefinition jet_def(fastjet::kt_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); // run the jet clustering with the above jet definition. hard event first fastjet::ClusterSequenceActiveArea clust_seq(hard_event, jet_def, area_spec); // 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); // 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::ClusterSequenceActiveArea &  , const vector< fastjet::PseudoJet > &  )

a function that pretty prints a list of jets, and performs the subtraction in two different ways Definition at line 162 of file fastjet_subtraction.cc. References fastjet::ClusterSequenceActiveArea::area(), fastjet::ClusterSequenceActiveArea::area_4vector(), fastjet::PseudoJet::E(), fastjet::PseudoJet::perp2(), fastjet::ClusterSequenceActiveArea::pt_per_unit_area(), and fastjet::sorted_by_pt(). { // sort jets into increasing pt vector<fastjet::PseudoJet> jets = sorted_by_pt(unsorted_jets); // the corrected jets will go in here vector<fastjet::PseudoJet> corrected_jets(jets.size()); // get median pt per unit area // NB pt_per_unit_area exists only in ClusterSequenceActiveArea, and // not in the base class ClusterSequenceWithArea. // This might change in a future release double median_pt_per_area = clust_seq.pt_per_unit_area(); printf(" ijet rap phi Pt area Pt corr (rap corr phi corr Pt corr)ext\n"); for (size_t j = 0; j < jets.size(); j++) { // get area of each jet double area = clust_seq.area(jets[j]); // "standard" correction. Subtract only the Pt double pt_corr = jets[j].perp() - area*median_pt_per_area; // "extended" correction fastjet::PseudoJet area_4vect = median_pt_per_area*clust_seq.area_4vector(jets[j]); if (area_4vect.perp2() >= jets[j].perp2() || area_4vect.E() >= jets[j].E()) { // if the correction is too large, set the jet to zero corrected_jets[j] = 0.0 * jets[j]; } else { // otherwise do an E-scheme subtraction corrected_jets[j] = jets[j] - area_4vect; } // NB We could also leave out the above "if": // corrected_jets[j] = jets[j] - area_4vect; // but the result would be different, since we would not avoid // jets with negative Pt or energy printf("%5u %7.3f %7.3f %9.3f %7.3f %9.3f %7.3f %7.3f %9.3f\n", j,jets[j].rap(),jets[j].phi(),jets[j].perp(), area, pt_corr, corrected_jets[j].rap(),corrected_jets[j].phi(), corrected_jets[j].perp()); } cout << endl; cout << "median pt_over_area = " << median_pt_per_area << endl << endl; }

---