example.cc
No OneTemporary
Actions

Size

12 KB

Subscribers

None

example.cc
View Options

	// Nsubjettiness Package
	// Questions/Comments? jthaler@jthaler.net
	//
	// Copyright (c) 2011-13
	// Jesse Thaler, Ken Van Tilburg, and Christopher K. Vermilion
	//
	//----------------------------------------------------------------------
	// This file is part of FastJet contrib.
	//
	// It 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.
	//
	// It 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 this code. If not, see <http://www.gnu.org/licenses/>.
	//----------------------------------------------------------------------


	#include <iomanip>
	#include <stdlib.h>
	#include <stdio.h>
	#include <time.h>
	#include <iostream>
	#include <istream>
	#include <fstream>
	#include <sstream>
	#include <string>

	#include "fastjet/PseudoJet.hh"
	#include <sstream>
	#include "Nsubjettiness.hh" // In external code, this should be fastjet/contrib/Nsubjettiness.hh
	#include "Njettiness.hh"
	#include "NjettinessPlugin.hh"

	using namespace std;
	using namespace fastjet;
	using namespace fastjet::contrib;

	// forward declaration to make things clearer
	void read_event(vector<PseudoJet> &event);

	////////
	//
	// Helper function for output
	//
	///////

	void PrintJets(const vector <PseudoJet>& jets) {

	if (jets.size() == 0) return;
	const NjettinessExtras * extras = njettiness_extras(jets[0]);

	if (jets[0].has_area()) {
	if (extras == NULL) {
	printf("%5s %15s %15s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt","m","e","area"); // label the columns
	for (unsigned int i = 0; i < jets.size(); i++) {
	printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n",i, jets[i].rap(),jets[i].phi(),jets[i].perp(),jets[i].m(),jets[i].e(),(jets[i].has_area() ? jets[i].area() : 0.0 ));
	}
	}
	else {
	fastjet::PseudoJet total(0,0,0,0);
	printf("%5s %15s %15s %15s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt","m","e","subTau","area"); // label the columns
	for (unsigned int i = 0; i < jets.size(); i++) {
	printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n",i, jets[i].rap(),jets[i].phi(),jets[i].perp(),jets[i].m(),jets[i].e(),extras->subTau(jets[i]),(jets[i].has_area() ? jets[i].area() : 0.0 ));
	total += jets[i];
	}
	printf("%5s %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n","total", total.rap(), total.phi(), total.perp(),total.m(),total.e(),extras->totalTau(),total.area());
	}
	} else {
	if (extras == NULL) {
	printf("%5s %15s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt","m","e"); // label the columns
	for (unsigned int i = 0; i < jets.size(); i++) {
	printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f\n",i, jets[i].rap(),jets[i].phi(),jets[i].perp(),jets[i].m(),jets[i].e());
	}
	}
	else {
	fastjet::PseudoJet total(0,0,0,0);
	printf("%5s %15s %15s %15s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt","m","e","subTau"); // label the columns
	for (unsigned int i = 0; i < jets.size(); i++) {
	printf("%5u %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n",i, jets[i].rap(),jets[i].phi(),jets[i].perp(),jets[i].m(),jets[i].e(),extras->subTau(jets[i]));
	total += jets[i];
	}
	printf("%5s %15.8f %15.8f %15.8f %15.8f %15.8f %15.8f\n","total", total.rap(), total.phi(), total.perp(),total.m(),total.e(),extras->totalTau());
	}
	}

	}


	void analyze(const vector<PseudoJet> & input_particles, int i_event, int n_event) {

	/////// N-subjettiness /////////////////////////////

	// Initial clustering with anti-kt algorithm
	JetAlgorithm algorithm = antikt_algorithm;
	double jet_rad = 1.00; // jet radius for anti-kt algorithm
	JetDefinition jetDef = JetDefinition(algorithm,jet_rad,E_scheme,Best);
	ClusterSequence clust_seq(input_particles,jetDef);
	vector<PseudoJet> antikt_jets = sorted_by_pt(clust_seq.inclusive_jets());

	// Defining Nsubjettiness parameters
	double beta = 1.0; // power for angular dependence, e.g. beta = 1 --> linear k-means, beta = 2 --> quadratic/classic k-means
	double R0 = 1.0; // Characteristic jet radius for normalization
	double Rcut = 1.0; // maximum R particles can be from axis to be included in jet

	for (int j = 0; j < 2; j++) { // Two hardest jets per event
	if (antikt_jets[j].perp() > 200) {
	vector<PseudoJet> jet_constituents = clust_seq.constituents(antikt_jets[j]);

	//
	// If you don't want subjets, you can use the simple functor Nsubjettiness:
	//

	// 1-subjettiness
	Nsubjettiness nSub1KT(1, Njettiness::kt_axes, beta, R0, Rcut);
	double tau1 = nSub1KT(antikt_jets[j]);
	Nsubjettiness nSub1Min(1, Njettiness::min_axes, beta, R0, Rcut);
	double tau1min = nSub1Min(antikt_jets[j]);
	Nsubjettiness nSub1OnePass(1, Njettiness::onepass_kt_axes, beta, R0, Rcut);
	double tau1onepass = nSub1OnePass(antikt_jets[j]);

	// 2-subjettiness
	Nsubjettiness nSub2KT(2, Njettiness::kt_axes, beta, R0, Rcut);
	double tau2 = nSub2KT(antikt_jets[j]);
	Nsubjettiness nSub2Min(2, Njettiness::min_axes, beta, R0, Rcut);
	double tau2min = nSub2Min(antikt_jets[j]);
	Nsubjettiness nSub2OnePass(2, Njettiness::onepass_kt_axes, beta, R0, Rcut);
	double tau2onepass = nSub2OnePass(antikt_jets[j]);

	// 3-subjettiness
	Nsubjettiness nSub3KT(3, Njettiness::kt_axes, beta, R0, Rcut);
	double tau3 = nSub3KT(antikt_jets[j]);
	Nsubjettiness nSub3Min(3, Njettiness::min_axes, beta, R0, Rcut);
	double tau3min = nSub3Min(antikt_jets[j]);
	Nsubjettiness nSub3OnePass(3, Njettiness::onepass_kt_axes, beta, R0, Rcut);
	double tau3onepass = nSub3OnePass(antikt_jets[j]);

	//
	// Or, if you want subjets, use the FastJet plugin on a jet's constituents
	//


	JetDefinition nsub_jetDef1(new NjettinessPlugin(1, Njettiness::kt_axes, 1.0, 1.0, 1.0));
	ClusterSequence nsub_seq1(antikt_jets[j].constituents(), nsub_jetDef1);
	vector<PseudoJet> kt1jets = nsub_seq1.inclusive_jets();

	JetDefinition nsub_jetDef2(new NjettinessPlugin(2, Njettiness::kt_axes, 1.0, 1.0, 1.0));
	ClusterSequence nsub_seq2(antikt_jets[j].constituents(), nsub_jetDef2);
	vector<PseudoJet> kt2jets = nsub_seq2.inclusive_jets();

	JetDefinition nsub_jetDef3(new NjettinessPlugin(3, Njettiness::kt_axes, 1.0, 1.0, 1.0));
	ClusterSequence nsub_seq3(antikt_jets[j].constituents(), nsub_jetDef3);
	vector<PseudoJet> kt3jets = nsub_seq3.inclusive_jets();

	// JetDefinition nsubMin_jetDef(new NjettinessPlugin(3, Njettiness::min_axes, 1.0, 1.0, 1.0));
	// ClusterSequence nsubMin_seq(antikt_jets[j].constituents(), nsubMin_jetDef);
	// vector<PseudoJet> min3jets = nsubMin_seq.inclusive_jets();

	JetDefinition nsubOnePass_jetDef1(new NjettinessPlugin(1, Njettiness::onepass_kt_axes, 1.0, 1.0, 1.0));
	ClusterSequence nsubOnePass_seq1(antikt_jets[j].constituents(), nsubOnePass_jetDef1);
	vector<PseudoJet> onepass1jets = nsubOnePass_seq1.inclusive_jets();

	JetDefinition nsubOnePass_jetDef2(new NjettinessPlugin(2, Njettiness::onepass_kt_axes, 1.0, 1.0, 1.0));
	ClusterSequence nsubOnePass_seq2(antikt_jets[j].constituents(), nsubOnePass_jetDef2);
	vector<PseudoJet> onepass2jets = nsubOnePass_seq2.inclusive_jets();

	JetDefinition nsubOnePass_jetDef3(new NjettinessPlugin(3, Njettiness::onepass_kt_axes, 1.0, 1.0, 1.0));
	ClusterSequence nsubOnePass_seq3(antikt_jets[j].constituents(), nsubOnePass_jetDef3);
	vector<PseudoJet> onepass3jets = nsubOnePass_seq3.inclusive_jets();

	// etc...

	// Print data to screen
	if (true) {
	printf("-------------------------------------------------------------------------------------"); printf("\n");
	printf("-------------------------------------------------------------------------------------"); printf("\n");
	cout << "Beta = " << beta << endl;
	cout << "kT Axes:" << endl;
	PrintJets(kt1jets);
	PrintJets(kt2jets);
	PrintJets(kt3jets);
	// cout << "Minimization Axes:" << endl;
	// PrintJets(min1jets);
	// PrintJets(min2jets);
	// PrintJets(min3jets);
	cout << "One Pass Minimization Axes from kT" << endl;
	PrintJets(onepass1jets);
	PrintJets(onepass2jets);
	PrintJets(onepass3jets);
	printf("-------------------------------------------------------------------------------------"); printf("\n");
	cout << "Beta = " << beta << endl;
	cout << " kT: " << "tau1: " << tau1 << " tau2: " << tau2 << " tau3: " << tau3 << " tau2/tau1: " << tau2/tau1 << " tau3/tau2: " << tau3/tau2 << endl;
	cout << " Min: " << "tau1: " << tau1min << " tau2: " << tau2min << " tau3: " << tau3min << " tau2/tau1: " << tau2min/tau1min << " tau3/tau2: " << tau3min/tau2min << endl;
	cout << "OnePass: " << "tau1: " << tau1onepass << " tau2: " << tau2onepass << " tau3: " << tau3onepass << " tau2/tau1: " << tau2onepass/tau1onepass << " tau3/tau2: " << tau3onepass/tau2onepass << endl;
	cout << endl;
	printf("-------------------------------------------------------------------------------------"); printf("\n");
	printf("-------------------------------------------------------------------------------------"); printf("\n");

	}
	}
	}


	////////// N-jettiness as a jet algorithm ///////////////////////////

	// WARNING: This is extremely preliminary
	// You can also find jets with Njettiness:

	// double ghost_maxrap = 5.0; // e.g. if particles go up to y=5
	// AreaDefinition area_def(active_area_explicit_ghosts, GhostedAreaSpec(ghost_maxrap));

	NjettinessPlugin njet_plugin(3, Njettiness::onepass_kt_axes, 1.0, 1.0, 1.0);
	JetDefinition njet_jetDef(&njet_plugin);
	// ClusterSequenceArea njet_seq(input_particles, njet_jetDef,area_def);
	ClusterSequence njet_seq(input_particles, njet_jetDef);
	vector<PseudoJet> njet_jets = njet_seq.inclusive_jets();

	NjettinessPlugin geo_plugin(3, NsubGeometricParameters(1.0));
	JetDefinition geo_jetDef(&geo_plugin);
	// ClusterSequenceArea geo_seq(input_particles, geo_jetDef,area_def);
	ClusterSequence geo_seq(input_particles, geo_jetDef);
	vector<PseudoJet> geo_jets = geo_seq.inclusive_jets();

	// The axes might point in a different direction than the jets
	// Using the NjettinessExtras pointer (ClusterSequence::Extras) to access that information
	vector<PseudoJet> njet_axes;
	// const NjettinessExtras * extras = njettiness_extras(njet_jets[0]);
	const NjettinessExtras * extras = njettiness_extras(njet_seq);
	if (extras != NULL) {
	njet_axes = extras->axes();
	}

	if (false) {
	printf("-------------------------------------------------------------------------------------"); printf("\n");
	cout << "Event-wide Jets from One-Pass Minimization (beta = 1.0)" << endl;
	PrintJets(njet_jets);
	cout << "Event-wide Axis Location for Above Jets" << endl;
	PrintJets(njet_axes);
	cout << "Event-wide Jets from Geometric Measure" << endl;
	PrintJets(geo_jets);
	printf("-------------------------------------------------------------------------------------"); printf("\n");
	}
	}




	//----------------------------------------------------------------------
	int main(){

	//----------------------------------------------------------
	// read in input particles
	vector<PseudoJet> event;
	read_event(event);
	cout << "# read an event with " << event.size() << " particles" << endl;

	//----------------------------------------------------------
	// illustrate how this Nsubjettiness contrib works

	analyze(event,0,0);

	return 0;
	}

	// read in input particles
	void read_event(vector<PseudoJet> &event){
	string line;
	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") {return;}
	if (line.substr(0,1) == "#") {continue;}
	double px,py,pz,E;
	linestream >> px >> py >> pz >> E;
	PseudoJet particle(px,py,pz,E);

	// push event onto back of full_event vector
	event.push_back(particle);
	}
	}

File Metadata

Mime Type: text/x-c
Expires: Sat, Dec 21, 5:12 PM (14 h, 49 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4023592
Default Alt Text: example.cc (12 KB)

example.ccNo OneTemporaryActions

example.ccView Options

File Metadata

Event Timeline

example.cc
No OneTemporary
Actions

example.cc
View Options