//----------------------------------------------------------------------
/// \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 2684 2011-11-14 07:41:44Z soyez $
//
// 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 .
//----------------------------------------------------------------------
//ENDHEADER
#include "fastjet/PseudoJet.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/Selector.hh"
#include "fastjet/tools/JetMedianBackgroundEstimator.hh"
#include "fastjet/tools/Subtractor.hh"
#include // needed for io
using namespace std;
using namespace fastjet;
int main(){
// 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 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
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;
JetDefinition jet_def(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;
GhostedAreaSpec area_spec(ghost_maxrap);
AreaDefinition area_def(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)
// ----------------------------------------------------------
ClusterSequenceArea clust_seq_hard(hard_event, jet_def, area_def);
ClusterSequenceArea clust_seq_full(full_event, jet_def, area_def);
double ptmin = 7.0;
vector hard_jets = sorted_by_pt(clust_seq_hard.inclusive_jets(ptmin));
vector full_jets = sorted_by_pt(clust_seq_full.inclusive_jets(ptmin));
// Now turn to the estimation of the background (for the full event)
//
// There are different ways to do that. In general, this also
// requires clustering the particles that will be handled internally
// in FastJet.
//
// The suggested way to proceed is to use a BackgroundEstimator
// constructed from the following 3 arguments:
// - a jet definition used to cluster the particles.
// . 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.
// - An area definition for which 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)
// - A Selector specifying the range over which we will keep the
// jets entering the estimation of the background (you should
// thus make sure the ghosts extend far enough in rapidity to
// cover the range, a warning will be issued otherwise).
// In this particular example, the two hardest jets in the event
// are removed from the background estimation
// ----------------------------------------------------------
JetDefinition jet_def_bkgd(kt_algorithm, 0.4);
AreaDefinition area_def_bkgd(active_area_explicit_ghosts,
GhostedAreaSpec(ghost_maxrap));
Selector selector = SelectorAbsRapMax(4.5) * (!SelectorNHardest(2));
JetMedianBackgroundEstimator bkgd_estimator(selector, jet_def_bkgd, area_def_bkgd);
// To help manipulate the background estimator, we also provide a
// transformer that allows to apply directly the background
// subtraction on the jets. This will use the background estimator
// to compute rho for the jets to be subtracted.
// ----------------------------------------------------------
Subtractor subtractor(&bkgd_estimator);
// Finally, once we have an event, we can just tell the background
// estimator to use that list of particles
// This could be done directly when declaring the background
// estimator but the usage below can more easily be accomodated to a
// loop over a set of events.
// ----------------------------------------------------------
bkgd_estimator.set_particles(full_event);
// 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 << " " << bkgd_estimator.description() << endl << endl;;
cout << " Giving, for the full event" << endl;
cout << " rho = " << bkgd_estimator.rho() << endl;
cout << " sigma = " << bkgd_estimator.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. Subtraction is performed on the
// full 4-vector
//
// 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;
// get the subtracted jets
vector subtracted_jets = subtractor(full_jets);
for (unsigned int i=0; i= 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_jets[i].rap(), subtracted_jets[i].phi(),
subtracted_jets[i].perp());
idx++;
}
}
return 0;
}