fastjet_subtraction.cc

//STARTHEADER
// $Id: fastjet_subtraction.cc 2577 2011-09-13 15:11:38Z salam $
//
// Copyright (c) 2005-2011, Matteo Cacciari, Gavin P. Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
//  FastJet is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2 of the License, or
//  (at your option) any later version.
//
//  The algorithms that underlie FastJet have required considerable
//  development and are described in hep-ph/0512210. If you use
//  FastJet as part of work towards a scientific publication, please
//  include a citation to the FastJet paper.
//
//  FastJet is distributed in the hope that it will be useful,
//  but WITHOUT ANY WARRANTY; without even the implied warranty of
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//  GNU General Public License for more details.
//
//  You should have received a copy of the GNU General Public License
//  along with FastJet. If not, see <http://www.gnu.org/licenses/>.
//----------------------------------------------------------------------
//ENDHEADER







//----------------------------------------------------------------------
// fastjet subtraction example program. 
//
// Compile it with: make fastjet_subtraction
// run it with    : ./fastjet_subtraction < data/Pythia-dijet-ptmin100-lhc-pileup-1ev.dat
//
//----------------------------------------------------------------------
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include<iostream> // needed for io
#include<sstream>  // needed for internal io
#include<vector> 

using namespace std;

// A declaration of a function that pretty prints a list of jets
// The subtraction is also performed inside this function
void print_jets (const fastjet::ClusterSequenceAreaBase &, 
                 const vector<fastjet::PseudoJet> &);

/// an example program showing how to use fastjet
int main (int argc, char ** argv) {
  
  // 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::GhostedAreaSpec ghost_spec(ghost_etamax, active_area_repeats, 
                                    ghost_area);
  fastjet::AreaDefinition area_def(ghost_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::ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def);

  // 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;


}


//----------------------------------------------------------------------
/// a function that pretty prints a list of jets, and performs the subtraction
/// in two different ways, using a generic ClusterSequenceAreaBase
/// type object.
void print_jets (const fastjet::ClusterSequenceAreaBase & clust_seq, 
                 const vector<fastjet::PseudoJet> & unsorted_jets ) {

  // 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     = jets[j].area();

    // "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*jets[j].area_4vector();
    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;


}