07-subtraction.cc

Go to the documentation of this file.

//----------------------------------------------------------------------
/// \file
/// \page Example07 07 - subtracting jet background contamination
///
/// fastjet subtraction example program. 
///
/// run it with    : ./07-subtraction < data/Pythia-Zp2jets-lhc-pileup-1ev.dat
///
/// Source code: 07-subtraction.cc
//----------------------------------------------------------------------

//STARTHEADER
// $Id: 07-subtraction.cc 2224 2011-06-03 14:22:31Z salam $
//
// Copyright (c) 2005-2011, Matteo Cacciari, Gavin 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, write to the Free Software
//  Foundation, Inc.:
//      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//----------------------------------------------------------------------
//ENDHEADER

#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include <iostream> // needed for io

using namespace std;

int main (int argc, char ** argv) {
  
  // read in input particles
  //
  // 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;
  
  // 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
  double particle_maxrap = 5.0;

  string line;
  int  nsub  = 0; // counter to keep track of which sub-event we're reading
  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 first one (the hard one) here
      if (nsub == 1) hard_event = full_event;
      nsub += 1;
    }
    if (line.substr(0,1) == "#") {continue;}
    double px,py,pz,E;
    linestream >> px >> py >> pz >> E;
    // you can construct 
    fastjet::PseudoJet particle(px,py,pz,E);

    // push event onto back of full_event vector
    if (abs(particle.rap()) <= particle_maxrap) full_event.push_back(particle);
  }

  // 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 a jet definition for the clustering
  // We use the anti-kt algorithm with a radius of 0.5
  //----------------------------------------------------------
  double R = 0.5;
  fastjet::JetDefinition jet_def(fastjet::antikt_algorithm, R);

  // create an area definition for the clustering
  //----------------------------------------------------------
  // ghosts should go up to the acceptance of the detector or
  // (with infinite acceptance) at least 2R beyond the region
  // where you plan to investigate jets.
  double ghost_maxrap = 6.0;
  fastjet::ActiveAreaSpec area_spec(ghost_maxrap);
  fastjet::AreaDefinition area_def(fastjet::active_area, area_spec);

  // run the jet clustering with the above jet and area definitions
  // for both the hard and full event
  //
  // We retrieve the jets above 7 GeV in both case (note that the
  // 7-GeV cut we be applied again later on after we subtract the jets
  // from the full event)
  // ----------------------------------------------------------
  fastjet::ClusterSequenceArea clust_seq_hard(hard_event, jet_def, area_def);
  fastjet::ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def);

  double ptmin = 7.0;
  vector<fastjet::PseudoJet> hard_jets = sorted_by_pt(clust_seq_hard.inclusive_jets(ptmin));
  vector<fastjet::PseudoJet> full_jets = sorted_by_pt(clust_seq_full.inclusive_jets(ptmin));

  // Now turn to the estimation of the background (for the full event)
  //
  // This also requires a ClusterSequenceArea.
  // In general, this ClusterSequenceArea does not need to be the same
  // as the one used (above) to cluster and extract the jets from the
  // event:
  //  - We strongly recommend using the kt or Cambridge/Aachen algorithm
  //    (a warning will be issued otherwise)
  //  - The choice of the radius is a bit more subtle. R=0.4 has been
  //    chosen to limit the impact of hard jets; in samples of
  //    dominantly sparse events it may cause the UE/pileup to be
  //    underestimated a little, a slightly larger value (0.5 or 0.6)
  //    may be better.
  //  - For the area definition, we recommend the use of explicit
  //    ghosts (i.e. active_area_explicit_ghosts)
  //    As mentionned in the area example (06-area.cc), ghosts should
  //    extend sufficiently far in rapidity to cover the jets used in
  //    the computation of the background (see also the comment below)
  //
  // ----------------------------------------------------------
  fastjet::JetDefinition jet_def_bkgd(fastjet::kt_algorithm, 0.4);
  fastjet::GhostedAreaSpec area_spec_bkgd(ghost_maxrap);
  fastjet::AreaDefinition area_def_bkgd(fastjet::active_area_explicit_ghosts, area_spec_bkgd);
  fastjet::ClusterSequenceArea clust_seq_bkgd(full_event, jet_def_bkgd, area_def_bkgd);

  // Once you have the ClusterSequenceArea, you can compute the
  // background. This is estimated over a given range
  // (RangeDefinition) i.e. only jets within that range will be used
  // to estimate the background. You shold thus make sure the ghosts
  // extend far enough in rapidity to cover the range, a warning will
  // be issued otherwise.
  //
  // The simplest way to define a RangeDefinition is through its
  // maximal |y| extent but other options are possible e.g. through a
  // minimal and maximal rapidity and minimal and maximal azimuthal
  // angle. If needed, you can even define your own ranges (a few are
  // provided with FastJet)
  //
  // Finally, the estimation of the background properties rho (the
  // average density per unit area) and sigma (the average
  // fluctuations per unit area) is done using
  // ClusterSequenceArea::get_median_rho_and_sigma(). This takes
  // 2 main parameters: the range discussed above and a boolean
  // controlling the use of 4-vector or scalar areas (we suggest using
  // 4-vector areas)
  //
  // ----------------------------------------------------------
  double range_maxrap = 4.5;  // we have a ghost_maxrap of 6.0, particles up to 5
  fastjet::RangeDefinition range(range_maxrap);

  bool use_4vector_area = true;
  
  double rho, sigma;
  clust_seq_bkgd.get_median_rho_and_sigma(range, use_4vector_area, rho, sigma);

  // show a summary of what was done so far
  //  - the description of the algorithms, areas and ranges used
  //  - the background properties
  //  - the jets in the hard event
  //----------------------------------------------------------
  cout << "Main clustering:" << endl;
  cout << "  Ran:   " << jet_def.description() << endl;
  cout << "  Area:  " << area_def.description() << endl;
  cout << "  Particles up to |y|=" << particle_maxrap << endl;
  cout << endl;

  cout << "Background estimation:" << endl;
  cout << "  Ran    " << jet_def_bkgd.description() << endl;
  cout << "  Area:  " << area_def_bkgd.description() << endl;
  cout << "  Range: " << range.description() << endl;
  cout << "  Giving, for the full event" << endl;
  cout << "    rho   = " << rho   << endl;
  cout << "    sigma = " << sigma << endl;
  cout << endl;

  cout << "Jets above " << ptmin << " GeV in the hard event (" << hard_event.size() << " particles)" << endl;
  cout << "---------------------------------------\n";
  printf("%5s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "area");
   for (unsigned int i = 0; i < hard_jets.size(); i++) {
    printf("%5u %15.8f %15.8f %15.8f %15.8f\n", i,
           hard_jets[i].rap(), hard_jets[i].phi(), hard_jets[i].perp(),
           hard_jets[i].area());
  }
  cout << endl;

  // Once the background properties have been computed, subtraction
  // can be applied on the jets
  //
  // This uses ClusterSequenceArea::subtracted_jet(jet, rho), with the
  // ClusterSequence used to cluster the jet and the background
  // density we have just computed
  // 
  // (Note that when using scalar areas, subtracted_pt should be used
  // instead of subtracted_jet)
  // 
  // We output the jets before and after subtraction
  // ----------------------------------------------------------
  cout << "Jets above " << ptmin << " GeV in the full event (" << full_event.size() << " particles)" << endl;
  cout << "---------------------------------------\n";
  printf("%5s %15s %15s %15s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt", "area", "rap_sub", "phi_sub", "pt_sub");
  unsigned int idx=0;
  for (unsigned int i=0; i<full_jets.size(); i++){
    // get the subtracted jet
    fastjet::PseudoJet subtracted_jet = clust_seq_full.subtracted_jet(full_jets[i], rho);

    // re-apply the pt cut
    if (subtracted_jet.perp2() >= ptmin*ptmin){
      printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n", idx,
             full_jets[i].rap(), full_jets[i].phi(), full_jets[i].perp(),
             full_jets[i].area(),
             subtracted_jet.rap(), subtracted_jet.phi(), 
             subtracted_jet.perp());
      idx++;
    }
  }

  return 0;
}