gSystem->Load("libfastjet");
Here we convert an ASCII file in ROOT format to then fill the particles necessary to Fastjet. Notice the inclusion of the header file. Given that we instructed ROOT about its location and the Fastjet library has been loaded, the function can be smoothly JIT-ted.
%%cpp -d
#include "fastjet/ClusterSequence.hh"
void fillInputParticles(vector<fastjet::PseudoJet>& input_particles, const char* inputFileName)
{
auto ntupleFormat = "px:py:pz:E";
TNtuple input_particles_ntuple("InputParticles","Input Particles",ntupleFormat);
auto n_particles = input_particles_ntuple.ReadFile(inputFileName,ntupleFormat);
input_particles.reserve(n_particles);
for (auto i : ROOT::TSeqI(n_particles)) {
input_particles_ntuple.GetEntry(i);
auto v = input_particles_ntuple.GetArgs();
fastjet::PseudoJet particle(v[0],v[1],v[2],v[3]);
particle.set_user_index(i);
input_particles.emplace_back(particle);
}
}
We now need to:
This function is flexible: the jet algorithm, its radius and the jets' minimum transverse momentum can be set by the user. This is useful to visualise the behaviour of the different algorithms.
%%cpp -d
THStack* getJetsComponents(const vector<fastjet::PseudoJet>& input_particles,
fastjet::JetAlgorithm jetAlgorithm = fastjet::antikt_algorithm,
double R = .6,
double ptMin = 14.)
{
// create a jet definition:
// a jet algorithm with a given radius parameter
fastjet::JetDefinition jet_def(jetAlgorithm, R);
fastjet::ClusterSequence clust_seq(input_particles, jet_def);
// get the resulting jets ordered in pt
auto inclusive_jets = sorted_by_pt(clust_seq.inclusive_jets(ptMin));
auto hs = new THStack("hs","Jets;rapidity;#phi");
int i=0;
for (auto&& jet : inclusive_jets) {
// get the constituents of the jet
auto constituents = jet.constituents();
auto h2Name = TString::Format("JetHist_%d",i);
if (auto oldH2 = (TH2F*) gDirectory->GetObjectChecked(h2Name,"TH2F")) {
delete oldH2;
}
auto h2 = new TH2F(h2Name,"JetHist",48, -5, 5, 48, 0, TMath::TwoPi());
h2->SetFillColor(i+1);
h2->SetLineWidth(1);
h2->SetLineColor(kBlack);
for (auto&& constituent : constituents){
h2->Fill(constituent.rap(), constituent.phi_02pi(), constituent.perp());
}
hs->Add(h2);
i++;
}
hs->Draw();// This creates the axes
hs->GetXaxis()->SetTitleOffset(1.7);
hs->GetYaxis()->SetTitleOffset(1.7);
return hs;
}
We now prepare all the elements needed to build our plot: a canvas to hold the graphics primitives, the containers and the stack of histograms containing the jets' components. First of all we uncompress the input file if needed.
%%python
inputFileName = 'Pythia-Z2jets-lhc-pileup-1ev.dat'
import os
if not os.path.exists(inputFileName):
import urllib2
response = urllib2.urlopen('https://raw.githubusercontent.com/dpiparo/swanExamples/master/notebooks/Pythia-Z2jets-lhc-pileup-1ev.dat')
filecontent = response.read()
with open(inputFileName,"w") as f_out:
f_out.write(filecontent)
TCanvas c;
c.SetTheta(90);
c.SetPhi(37);
unique_ptr<THStack> jetComponents(nullptr);
vector<fastjet::PseudoJet> input_particles;
fillInputParticles(input_particles, "Pythia-Z2jets-lhc-pileup-1ev.dat");
We are ready to go now: we can play with the parameters of the clustering and visualise the result. Each jet is visualised with a different colour in the Rapidity-Phi plane.
jetComponents.reset(getJetsComponents(input_particles, fastjet::antikt_algorithm, .8, 14.));
jetComponents->Draw("LEGO1 0");
c.Draw();
#-------------------------------------------------------------------------- # FastJet release 3.1.1 # M. Cacciari, G.P. Salam and G. Soyez # A software package for jet finding and analysis at colliders # http://fastjet.fr # # Please cite EPJC72(2012)1896 [arXiv:1111.6097] if you use this package # for scientific work and optionally PLB641(2006)57 [hep-ph/0512210]. # # FastJet is provided without warranty under the terms of the GNU GPLv2. # It uses T. Chan's closest pair algorithm, S. Fortune's Voronoi code # and 3rd party plugin jet algorithms. See COPYING file for details. #--------------------------------------------------------------------------