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 ActiveAreaSpec, Best, kt_algorithm, and print_jets().

                                  {
  
  // 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.

Definition at line 170 of file fastjet_subtraction.cc.

References fastjet::d0::inline_maths::phi(), and 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 -- you must decide the rapidity range
  // in which you measure the activity. (If doing a ghost-based active
  // area, make sure that your ghosts go up at least to \sim range+R).
  double range = 5.0;
  double median_pt_per_area = clust_seq.median_pt_per_unit_area(range);
  double median_pt_per_area4vector = clust_seq.median_pt_per_unit_area_4vector(range);

  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 sub_4vect = 
                       median_pt_per_area4vector*clust_seq.area_4vector(jets[j]);
    if (sub_4vect.perp2() >= jets[j].perp2() || 
        sub_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] - sub_4vect;
    }
    // NB We could also leave out the above "if": 
    // corrected_jets[j] = jets[j] - sub_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;
  cout << "median pt_over_area4vector = " << median_pt_per_area4vector << endl << endl;


}