No OneTemporary
Actions

Size

63 KB

Subscribers

None

View Options

	Index: contrib/contribs/VariableR/trunk/example.cc
	===================================================================
	--- contrib/contribs/VariableR/trunk/example.cc (revision 907)
	+++ contrib/contribs/VariableR/trunk/example.cc (revision 908)
	@@ -1,283 +1,290 @@
	// VariableR Package
	// Questions/Comments? jthaler@jthaler.net
	//
	// Copyright (c) 2009-2016
	// David Krohn, Gregory Soyez, Jesse Thaler, and Lian-Tao Wang
	//
	// $Id$
	//----------------------------------------------------------------------
	// 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 <iostream>
	#include <sstream>
	#include <stdio.h>
	#include <ctime>


	#include "fastjet/PseudoJet.hh"
	#include <sstream>
	#include "VariableRPlugin.hh" // In external code, this should be fastjet/contrib/VariableRPlugin.hh
	//#include "VariableR.hh" // This header is still available for backwards compatibility.

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

	void print_jets (const fastjet::ClusterSequence &,
	const vector<fastjet::PseudoJet> &);

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

	//----------------------------------------------------------------------
	int main(){
	-
	- //----------------------------------------------------------
	- // read in input particles
	- vector<PseudoJet> event;
	- read_event(event);
	- cout << "# read an event with " << event.size() << " particles" << endl;
	-
	- //----------------------------------------------------------
	- // illustrate how VariableR contrib works
	- // anti-kT variable R
	- //----------------------------------------------------------
	-
	- // defining parameters
	- double rho(2000.), min_r(0.0), max_r(2.0),ptmin(5.0);
	-
	- VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE);
	- fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	- fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	-
	- // tell the user what was done
	- cout << "# Ran " << jet_defAKT.description() << endl;
	-
	- // extract the inclusive jets with pt > 5 GeV
	- vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	-
	- // print them out
	- cout << "Printing inclusive jets with pt > "<< ptmin <<" GeV\n";
	- cout << "---------------------------------------\n";
	- print_jets(clust_seqAKT, inclusive_jetsAKT);
	- cout << endl;
	-
	- //----------------------------------------------------------
	- // Same for Cambridge-Aachen variable R
	- //----------------------------------------------------------
	-
	- VariableRPlugin lvjet_pluginCA(rho, min_r, max_r, VariableRPlugin::CALIKE);
	- fastjet::JetDefinition jet_defCA(&lvjet_pluginCA);
	- fastjet::ClusterSequence clust_seqCA(event, jet_defCA);
	-
	- // tell the user what was done
	- cout << "# Ran " << jet_defCA.description() << endl;
	-
	- // extract the inclusive jets with pt > 5 GeV
	- vector<fastjet::PseudoJet> inclusive_jetsCA = clust_seqCA.inclusive_jets(ptmin);
	-
	- // print them out
	- cout << "Printing inclusive jets with pt > "<< ptmin <<" GeV\n";
	- cout << "---------------------------------------\n";
	- print_jets(clust_seqCA, inclusive_jetsCA);
	- cout << endl;
	-
	- //----------------------------------------------------------
	- // Illustrating preclustering feature new in v1.1
	- // Need to have minimum jet radius for preclustering to make sense
	- //----------------------------------------------------------
	-
	- min_r = 0.4; // change small radius to allow for preclustering
	- VariableRPlugin lvjet_pluginAKT_precluster(rho, min_r, max_r, VariableRPlugin::AKTLIKE,true);
	- fastjet::JetDefinition jet_defAKT_precluster(&lvjet_pluginAKT_precluster);
	- fastjet::ClusterSequence clust_seqAKT_precluster(event, jet_defAKT_precluster);
	-
	- // tell the user what was done
	- cout << "# Ran " << jet_defAKT_precluster.description() << endl;
	-
	- // extract the inclusive jets with pt > 5 GeV
	- vector<fastjet::PseudoJet> inclusive_jetsAKT_precluster = clust_seqAKT_precluster.inclusive_jets(ptmin);
	-
	- // print them out
	- cout << "Printing inclusive jets with pt > "<< ptmin <<" GeV\n";
	- cout << "---------------------------------------\n";
	- print_jets(clust_seqAKT_precluster, inclusive_jetsAKT_precluster);
	- cout << endl;
	-
	- //----------------------------------------------------------
	- // As a cross check, same as above, but with no preclustering
	- // (Results should be nearly identical)
	- //----------------------------------------------------------
	-
	- VariableRPlugin lvjet_pluginAKT_noprecluster(rho, min_r, max_r, VariableRPlugin::AKTLIKE,false);
	- fastjet::JetDefinition jet_defAKT_noprecluster(&lvjet_pluginAKT_noprecluster);
	- fastjet::ClusterSequence clust_seqAKT_noprecluster(event, jet_defAKT_noprecluster);
	-
	- // tell the user what was done
	- cout << "# Ran " << jet_defAKT_noprecluster.description() << endl;
	-
	- // extract the inclusive jets with pt > 5 GeV
	- vector<fastjet::PseudoJet> inclusive_jetsAKT_noprecluster = clust_seqAKT_noprecluster.inclusive_jets(ptmin);
	-
	- // print them out
	- cout << "Printing inclusive jets with pt > "<< ptmin <<" GeV\n";
	- cout << "---------------------------------------\n";
	- print_jets(clust_seqAKT_noprecluster, inclusive_jetsAKT_noprecluster);
	- cout << endl;
	-
	- //----------------------------------------------------------
	- // An example with generalised kt
	- //----------------------------------------------------------
	- VariableRPlugin lvjet_pluginGenKT(rho, min_r, max_r, 0.5);
	- fastjet::JetDefinition jet_defGenKT(&lvjet_pluginGenKT);
	- fastjet::ClusterSequence clust_seqGenKT(event, jet_defGenKT);
	-
	- // tell the user what was done
	- cout << "# Ran " << jet_defGenKT.description() << endl;
	-
	- // extract the inclusive jets with pt > 5 GeV
	- vector<fastjet::PseudoJet> inclusive_jetsGenKT = clust_seqGenKT.inclusive_jets(ptmin);
	-
	- // print them out
	- cout << "Printing inclusive jets with pt > "<< ptmin <<" GeV\n";
	- cout << "---------------------------------------\n";
	- print_jets(clust_seqGenKT, inclusive_jetsGenKT);
	- cout << endl;
	-
	-
	-
	- // timing tests for the developers
	- double do_timing_test = true;
	- if (do_timing_test) {
	-
	- clock_t clock_begin, clock_end;
	- double num_iter;
	-
	- num_iter = 10;
	- min_r = 0.4;
	-
	- // Testing Native
	-
	- clock_begin = clock();
	- for (int t = 0; t < num_iter; t++) {
	- VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE,false,VariableRPlugin::Native);
	- fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	- fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	- vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	- }
	- clock_end = clock();
	- cout << "# Native: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per event"<< endl;
	-
	-
	- // Testing NNH
	-
	- clock_begin = clock();
	- for (int t = 0; t < num_iter; t++) {
	- VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE,false,VariableRPlugin::NNH);
	- fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	- fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	- vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	- }
	- clock_end = clock();
	- cout << "# NNH: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per event"<< endl;
	-
	+
	+ //----------------------------------------------------------
	+ // read in input particles
	+ vector<PseudoJet> event;
	+ read_event(event);
	+ cout << "# read an event with " << event.size() << " particles" << endl;
	+
	+ //----------------------------------------------------------
	+ // illustrate how VariableR contrib works
	+ // anti-kT variable R
	+ //----------------------------------------------------------
	+
	+ // defining parameters
	+ double rho = 2000.0;
	+ double min_r = 0.0;
	+ double max_r = 2.0;
	+ double ptmin = 5.0;
	+
	+ VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE);
	+ fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	+ fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	+
	+ // tell the user what was done
	+ cout << "# Ran " << jet_defAKT.description() << endl;
	+
	+ // extract the inclusive jets with pt > 5 GeV
	+ vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	+
	+ // print them out
	+ cout << "Printing inclusive jets with pt > "<< ptmin <<" GeV\n";
	+ cout << "---------------------------------------\n";
	+ print_jets(clust_seqAKT, inclusive_jetsAKT);
	+ cout << endl;
	+
	+ //----------------------------------------------------------
	+ // Same for Cambridge-Aachen variable R
	+ //----------------------------------------------------------
	+
	+ VariableRPlugin lvjet_pluginCA(rho, min_r, max_r, VariableRPlugin::CALIKE);
	+ fastjet::JetDefinition jet_defCA(&lvjet_pluginCA);
	+ fastjet::ClusterSequence clust_seqCA(event, jet_defCA);
	+
	+ // tell the user what was done
	+ cout << "# Ran " << jet_defCA.description() << endl;
	+
	+ // extract the inclusive jets with pt > 5 GeV
	+ vector<fastjet::PseudoJet> inclusive_jetsCA = clust_seqCA.inclusive_jets(ptmin);
	+
	+ // print them out
	+ cout << "Printing inclusive jets with pt > "<< ptmin <<" GeV\n";
	+ cout << "---------------------------------------\n";
	+ print_jets(clust_seqCA, inclusive_jetsCA);
	+ cout << endl;
	+
	+ //----------------------------------------------------------
	+ // Illustrating preclustering feature new in v1.1
	+ // Need to have minimum jet radius for preclustering to make sense
	+ //----------------------------------------------------------
	+
	+ min_r = 0.4; // change small radius to allow for preclustering
	+ bool use_preclustering = true;
	+ VariableRPlugin lvjet_pluginAKT_precluster(rho, min_r, max_r, VariableRPlugin::AKTLIKE,use_preclustering);
	+ fastjet::JetDefinition jet_defAKT_precluster(&lvjet_pluginAKT_precluster);
	+ fastjet::ClusterSequence clust_seqAKT_precluster(event, jet_defAKT_precluster);
	+
	+ // tell the user what was done
	+ cout << "# Ran " << jet_defAKT_precluster.description() << endl;
	+
	+ // extract the inclusive jets with pt > 5 GeV
	+ vector<fastjet::PseudoJet> inclusive_jetsAKT_precluster = clust_seqAKT_precluster.inclusive_jets(ptmin);
	+
	+ // print them out
	+ cout << "Printing inclusive jets with pt > "<< ptmin <<" GeV\n";
	+ cout << "---------------------------------------\n";
	+ print_jets(clust_seqAKT_precluster, inclusive_jetsAKT_precluster);
	+ cout << endl;
	+
	+ //----------------------------------------------------------
	+ // As a cross check, same as above, but with no preclustering
	+ // (Results should be nearly identical)
	+ //----------------------------------------------------------
	+
	+ use_preclustering = false;
	+ VariableRPlugin lvjet_pluginAKT_noprecluster(rho, min_r, max_r, VariableRPlugin::AKTLIKE,use_preclustering);
	+ fastjet::JetDefinition jet_defAKT_noprecluster(&lvjet_pluginAKT_noprecluster);
	+ fastjet::ClusterSequence clust_seqAKT_noprecluster(event, jet_defAKT_noprecluster);
	+
	+ // tell the user what was done
	+ cout << "# Ran " << jet_defAKT_noprecluster.description() << endl;
	+
	+ // extract the inclusive jets with pt > 5 GeV
	+ vector<fastjet::PseudoJet> inclusive_jetsAKT_noprecluster = clust_seqAKT_noprecluster.inclusive_jets(ptmin);
	+
	+ // print them out
	+ cout << "Printing inclusive jets with pt > "<< ptmin <<" GeV\n";
	+ cout << "---------------------------------------\n";
	+ print_jets(clust_seqAKT_noprecluster, inclusive_jetsAKT_noprecluster);
	+ cout << endl;
	+
	+ //----------------------------------------------------------
	+ // An example with generalised kt
	+ //----------------------------------------------------------
	+
	+ VariableRPlugin lvjet_pluginGenKT(rho, min_r, max_r, -0.5); // p = -0.5 (halfway between CA and anti-kt)
	+ fastjet::JetDefinition jet_defGenKT(&lvjet_pluginGenKT);
	+ fastjet::ClusterSequence clust_seqGenKT(event, jet_defGenKT);
	+
	+ // tell the user what was done
	+ cout << "# Ran " << jet_defGenKT.description() << endl;
	+
	+ // extract the inclusive jets with pt > 5 GeV
	+ vector<fastjet::PseudoJet> inclusive_jetsGenKT = clust_seqGenKT.inclusive_jets(ptmin);
	+
	+ // print them out
	+ cout << "Printing inclusive jets with pt > "<< ptmin <<" GeV\n";
	+ cout << "---------------------------------------\n";
	+ print_jets(clust_seqGenKT, inclusive_jetsGenKT);
	+ cout << endl;
	+
	+
	+
	+ // timing tests for the developers
	+ double do_timing_test = false;
	+ if (do_timing_test) {
	+
	+ clock_t clock_begin, clock_end;
	+ double num_iter;
	+
	+ num_iter = 10;
	+ min_r = 0.4;
	+ use_preclustering = false;
	+
	+ // Testing Native
	+
	+ clock_begin = clock();
	+ for (int t = 0; t < num_iter; t++) {
	+ VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE,use_preclustering,VariableRPlugin::Native);
	+ fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	+ fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	+ vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	+ }
	+ clock_end = clock();
	+ cout << "# Native: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per event"<< endl;
	+
	+
	+ // Testing NNH
	+
	+ clock_begin = clock();
	+ for (int t = 0; t < num_iter; t++) {
	+ VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE,use_preclustering,VariableRPlugin::NNH);
	+ fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	+ fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	+ vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	+ }
	+ clock_end = clock();
	+ cout << "# NNH: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per event"<< endl;
	+
	#if FASTJET_VERSION_NUMBER >= 30200
	-
	- // Testing N2Tiled
	-
	- clock_begin = clock();
	- for (int t = 0; t < num_iter; t++) {
	- VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE,false,VariableRPlugin::N2Tiled);
	- fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	- fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	- vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	- }
	- clock_end = clock();
	- cout << "# N2Tiled: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per event"<< endl;
	-
	-
	- // Testing N2Plain
	-
	- clock_begin = clock();
	- for (int t = 0; t < num_iter; t++) {
	- VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE,false,VariableRPlugin::N2Plain);
	- fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	- fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	- vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	- }
	- clock_end = clock();
	- cout << "# N2Plain: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per event"<< endl;
	-
	-
	+
	+ // Testing N2Tiled
	+
	+ clock_begin = clock();
	+ for (int t = 0; t < num_iter; t++) {
	+ VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE,use_preclustering,VariableRPlugin::N2Tiled);
	+ fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	+ fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	+ vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	+ }
	+ clock_end = clock();
	+ cout << "# N2Tiled: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per event"<< endl;
	+
	+
	+ // Testing N2Plain
	+
	+ clock_begin = clock();
	+ for (int t = 0; t < num_iter; t++) {
	+ VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE,use_preclustering,VariableRPlugin::N2Plain);
	+ fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	+ fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	+ vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	+ }
	+ clock_end = clock();
	+ cout << "# N2Plain: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per event"<< endl;
	+
	+
	#endif
	-
	- // Testing Best
	-
	- clock_begin = clock();
	- for (int t = 0; t < num_iter; t++) {
	- VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE,false,VariableRPlugin::Best);
	- fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	- fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	- vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	- }
	- clock_end = clock();
	- cout << "# Best: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per event"<< endl;
	-
	- }
	-
	-
	- return 0;
	+
	+ // Testing Best
	+
	+ clock_begin = clock();
	+ for (int t = 0; t < num_iter; t++) {
	+ VariableRPlugin lvjet_pluginAKT(rho, min_r, max_r, VariableRPlugin::AKTLIKE,use_preclustering,VariableRPlugin::Best);
	+ fastjet::JetDefinition jet_defAKT(&lvjet_pluginAKT);
	+ fastjet::ClusterSequence clust_seqAKT(event, jet_defAKT);
	+ vector<fastjet::PseudoJet> inclusive_jetsAKT = clust_seqAKT.inclusive_jets(ptmin);
	+ }
	+ clock_end = clock();
	+ cout << "# Best: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per event"<< endl;
	+
	+ }
	+
	+
	+ 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);
	- }
	+ 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);
	+ }
	}

	//----------------------------------------------------------------------
	/// a function that pretty prints a list of jets
	void print_jets (const fastjet::ClusterSequence & clust_seq,
	const vector<fastjet::PseudoJet> & jets) {
	-
	- // sort jets into increasing pt
	- vector<fastjet::PseudoJet> sorted_jets = sorted_by_pt(jets);
	-
	- // label the columns
	- printf("%5s %10s %10s %10s %10s %10s %10s\n","jet #", "rapidity",
	- "phi", "pt","m","e", "n constituents");
	-
	- // print out the details for each jet
	- for (unsigned int i = 0; i < sorted_jets.size(); i++) {
	- int n_constituents = clust_seq.constituents(sorted_jets[i]).size();
	- printf("%5u %10.3f %10.3f %10.3f %10.3f %10.3f %8u\n",
	- i, sorted_jets[i].rap(), sorted_jets[i].phi(),
	- sorted_jets[i].perp(),sorted_jets[i].m(),sorted_jets[i].e(), n_constituents);
	- }
	+
	+ // sort jets into increasing pt
	+ vector<fastjet::PseudoJet> sorted_jets = sorted_by_pt(jets);
	+
	+ // label the columns
	+ printf("%5s %10s %10s %10s %10s %10s %10s\n","jet #", "rapidity",
	+ "phi", "pt","m","e", "n constituents");
	+
	+ // print out the details for each jet
	+ for (unsigned int i = 0; i < sorted_jets.size(); i++) {
	+ int n_constituents = clust_seq.constituents(sorted_jets[i]).size();
	+ printf("%5u %10.3f %10.3f %10.3f %10.3f %10.3f %8u\n",
	+ i, sorted_jets[i].rap(), sorted_jets[i].phi(),
	+ sorted_jets[i].perp(),sorted_jets[i].m(),sorted_jets[i].e(), n_constituents);
	+ }
	}
	Index: contrib/contribs/VariableR/trunk/VariableRPlugin.hh
	===================================================================
	--- contrib/contribs/VariableR/trunk/VariableRPlugin.hh (revision 907)
	+++ contrib/contribs/VariableR/trunk/VariableRPlugin.hh (revision 908)
	@@ -1,193 +1,193 @@
	// VariableR Package
	// Questions/Comments? jthaler@jthaler.net
	//
	// Copyright (c) 2009-2016
	// David Krohn, Gregory Soyez, Jesse Thaler, and Lian-Tao Wang
	//
	// $Id$
	//----------------------------------------------------------------------
	// 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/>.
	//----------------------------------------------------------------------

	#ifndef __FASTJET_CONTRIB_VARIABLERPLUGIN_HH__
	#define __FASTJET_CONTRIB_VARIABLERPLUGIN_HH__

	#include <fastjet/internal/base.hh>
	#include <fastjet/config.h>

	#include "fastjet/JetDefinition.hh"
	#include "fastjet/PseudoJet.hh"
	#include "fastjet/ClusterSequence.hh"
	#include <fastjet/LimitedWarning.hh>

	#include <map>
	#include <queue>

	using namespace std;

	FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh

	namespace contrib {

	////////
	//
	// Core VR Code
	//
	////////

	// In version 1.1, this replaces CoreJetAlgorithm.
	// This acts like any fastjet plugin since it implements run_clustering
	class VariableRPlugin : public JetDefinition::Plugin {

	public:
	/// Type of clustering
	///
	/// Since version 1.2.0 of VariableR, the clustering is treated as
	/// a generalised-kt algorithm and the previous "ClusterType"
	/// argument is replaced by the "p" parameter of the
	/// generalised-kt algorithm. The definitions below are shorthand
	/// for the antikt, C/A and kt algorithm which also allow for
	/// backwards compatibility.
	/// (These are initialized in VariableRPlugin.cc.)
	static const double CALIKE; // = 0.0;
	static const double KTLIKE; // = 1.0;
	static const double AKTLIKE; // = -1.0;

	- /// for backwards compatibility reason we also define ClusterType
	+ /// for backwards compatibility reasons, we also define ClusterType
	/// as "double"
	typedef double ClusterType;

	/// The strategy to be used with the clustering
	enum Strategy{
	Best, ///< currently N2Tiled or N2Plain for FJ>3.2.0, Native for FastJet<3.2.0
	N2Tiled, ///< the default (faster in most cases) [requires FastJet>=3.2.0]
	N2Plain, ///< [requires FastJet>=3.2.0]
	NNH, ///< slower but already available for FastJet<3.2.0
	Native ///< original local implemtation of the clustering [the default for FastJet<3.2.0]
	};


	- /// constructor that sets VR algorithm parameters
	+ /// Constructor that sets VR algorithm parameters
	///
	/// \param rho mass scale for effective radius (i.e. R ~ rho/pT)
	/// \param min_r minimum jet radius
	/// \param max_r maximum jet radius
	/// \param clust_type whether to use CA-like, kT-like, or anti-kT-like distance measure
	/// (this value is the same as the p exponent in generalized-kt, with
	/// anti-kt = -1.0, CA = 0.0, and kT = 1.0)
	/// \param precluster whether to use optional kT subjets (of size min_r) for preclustering
	/// (true is much faster, default=false). At the moment, the only option
	/// for preclustering is kT (use fastjet::NestedDefsPjugin otherwise)
	/// \param strategy decodes which algorithm to apply for the clustering
	///
	/// Note that pre-clustering is deprecated and will likely be
	/// removed in a future releasse of this contrib. You can get the
	/// same behaviour using the NestedDefs fastjet plugin:
	///
	/// \code
	/// std::list<JetDefinition> jet_defs;
	/// jet_defs.push_back(JetDefinition(kt_algorithm, min_r));
	/// contrib::VariableRPlugin vr_plugin(rho,min_r,max_r,clust_type);
	/// jet_defs.push_back(JetDefinition(&vr_plugin));
	///
	/// NestedDefsPlugin nd_plugin(jet_defs);
	/// JetDefinition jet_def(&nd_plugin);
	/// \endcode
	///
	/// For FastJet>=3.2.0, we will use by default the N2Tiled
	/// strategy unless pre-clustering is requested. Otherwise, we use
	/// the "old" native strategy.
	VariableRPlugin(double rho, double min_r, double max_r, ClusterType clust_type, bool precluster = false,
	Strategy requested_strategy = Best);

	// virtual function from JetDefinition::Plugin that implements the actual VR algorithm
	void run_clustering(ClusterSequence & cs) const;

	// information string
	virtual string description() const;

	// TODO: have this return a non-trivial answer.
	virtual double R() const{ return _max_r; }

	private:

	// parameters to define VR algorithm
	double _rho2, _min_r2, _max_r, _max_r2;
	ClusterType _clust_type;
	Strategy _requested_strategy;

	// For preclustering, can use kT algorithm to make subclusters of size min_r
	bool _precluster;
	JetDefinition _pre_jet_def;

	// warn about pre-clustering being deprecated
	static LimitedWarning _preclustering_deprecated_warning;

	// helper function to decide what strategy is best
	//
	// The optimal strategy will depend on the multiplicity and _max_r
	Strategy _best_strategy(unsigned int N) const;

	// helper function to apply kT preclustering if desired
	void _preclustering(ClusterSequence & cs, set<int>& unmerged_jets) const;

	// native implementation of the clustering
	void _native_clustering(ClusterSequence &cs) const;

	// implementation of the clustering using FastJet NN*** classes
	template<typename NN>
	void _NN_clustering(ClusterSequence &cs, NN &nn) const;

	// Helper struct to store two jets and a distance measure
	struct JetDistancePair{
	int j1,j2;
	double distance;
	};

	// Helper comparitor class for comparing JetDistancePairs
	class CompareJetDistancePair {
	public:
	CompareJetDistancePair(){};
	bool operator() (const JetDistancePair & lhs, const JetDistancePair &rhs) const {
	return (lhs.distance > rhs.distance);
	}
	};

	// helper function to merge jet with beam
	void _merge_jet_with_beam(ClusterSequence & clust_seq, JetDistancePair & jdp, set<int>& unmerged_jets) const;

	// helper function to merge two jets into a new pseudojet
	void _merge_jets(ClusterSequence & clust_seq,
	JetDistancePair & jdp,
	priority_queue< JetDistancePair, vector<JetDistancePair>, CompareJetDistancePair > &jet_queue,
	set<int>& unmerged_jets) const;

	// use ClusterMode to determine jet-jet and jet-beam distance
	inline double _get_JJ_distance_measure(const PseudoJet& j1, const PseudoJet& j2) const;
	inline double _get_JB_distance_measure(const PseudoJet& jet) const;

	// helper function to establish measures
	void _setup_distance_measures(ClusterSequence & clust_seq,
	vector<JetDistancePair> &jet_vec,
	set<int>& unmerged_jets) const;

	};

	} // namespace contrib

	FASTJET_END_NAMESPACE

	#endif // __FASTJET_CONTRIB_VARIABLERPLUGIN_HH__
	Index: contrib/contribs/VariableR/trunk/example.ref
	===================================================================
	--- contrib/contribs/VariableR/trunk/example.ref (revision 907)
	+++ contrib/contribs/VariableR/trunk/example.ref (revision 908)
	@@ -1,70 +1,71 @@
	# read an event with 354 particles
	#--------------------------------------------------------------------------
	-# FastJet release 3.0.6
	+# FastJet release 3.2.0-devel
	# 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.
	#--------------------------------------------------------------------------
	-# Ran Variable R (0903.0392), AKT, rho=2000.0, min_r=0.0, max_r=2.0
	+# Ran Variable R (0903.0392), AKT, rho=2000.0, min_r=0.0, max_r=2.0, strategy=Best
	Printing inclusive jets with pt > 5 GeV
	---------------------------------------
	jet # rapidity phi pt m e n constituents
	0 -0.881 2.908 991.422 222.818 1436.367 65
	1 0.095 6.046 986.931 461.292 1094.380 121
	2 -3.376 0.863 18.189 32.359 543.629 45
	3 3.487 5.797 8.210 17.192 311.682 34
	4 -4.586 3.971 7.090 15.267 825.609 28
	5 5.824 3.590 7.032 12.730 2459.371 22

	-# Ran Variable R (0903.0392), CA, rho=2000.0, min_r=0.0, max_r=2.0
	+# Ran Variable R (0903.0392), CA, rho=2000.0, min_r=0.0, max_r=2.0, strategy=Best
	Printing inclusive jets with pt > 5 GeV
	---------------------------------------
	jet # rapidity phi pt m e n constituents
	0 -0.876 2.915 990.767 187.190 1420.188 58
	1 0.086 6.046 984.618 434.657 1080.312 103
	2 -2.301 0.993 17.150 19.523 131.047 39
	3 2.421 5.294 9.652 20.340 127.684 35
	4 -4.548 1.813 7.785 9.217 570.017 22
	5 5.726 3.628 6.924 12.983 2256.767 21
	6 -4.505 4.802 6.596 16.142 788.616 28

	-# Ran Variable R (0903.0392), AKT, rho=2000.0, min_r=0.4, max_r=2.0, with precluster
	+# Ran Variable R (0903.0392), AKT, rho=2000.0, min_r=0.4, max_r=2.0, with precluster, strategy=Best
	Printing inclusive jets with pt > 5 GeV
	---------------------------------------
	jet # rapidity phi pt m e n constituents
	0 -0.881 2.907 991.435 225.793 1437.746 68
	1 0.096 6.045 986.968 460.724 1094.229 121
	2 -3.391 0.854 17.828 33.757 567.489 46
	3 3.536 5.855 8.089 15.905 306.543 31
	4 -4.651 3.961 6.808 13.699 800.780 26
	5 5.904 3.612 6.735 11.543 2449.325 21

	-# Ran Variable R (0903.0392), AKT, rho=2000.0, min_r=0.4, max_r=2.0
	+# Ran Variable R (0903.0392), AKT, rho=2000.0, min_r=0.4, max_r=2.0, strategy=Best
	Printing inclusive jets with pt > 5 GeV
	---------------------------------------
	jet # rapidity phi pt m e n constituents
	0 -0.881 2.908 991.422 222.818 1436.367 65
	1 0.095 6.046 986.931 461.292 1094.380 121
	2 -3.376 0.863 18.189 32.359 543.629 45
	3 3.487 5.797 8.210 17.192 311.682 34
	4 -4.586 3.971 7.090 15.267 825.609 28
	5 5.824 3.590 7.032 12.730 2459.371 22

	-# Ran Variable R (0903.0392), GenKT(p=0.5), rho=2000.0, min_r=0.4, max_r=2.0
	+# Ran Variable R (0903.0392), GenKT(p=-0.5), rho=2000.0, min_r=0.4, max_r=2.0, strategy=Best
	Printing inclusive jets with pt > 5 GeV
	---------------------------------------
	jet # rapidity phi pt m e n constituents
	- 0 0.005 6.059 991.515 728.701 1230.506 172
	- 1 -0.876 2.915 990.891 181.336 1419.597 55
	- 2 -4.402 1.319 9.257 14.931 716.854 32
	- 3 5.584 3.574 7.527 18.137 2611.993 30
	- 4 3.457 5.876 7.302 13.688 246.297 29
	- 5 -4.723 4.292 6.509 9.647 654.821 21
	+ 0 -0.881 2.908 991.422 222.818 1436.367 65
	+ 1 0.095 6.046 986.931 461.292 1094.380 121
	+ 2 -3.208 0.702 16.658 25.928 381.653 40
	+ 3 3.407 5.842 7.882 14.301 246.597 30
	+ 4 5.627 3.643 7.650 16.535 2531.455 26
	+ 5 -4.681 4.291 6.710 10.160 656.977 22
	+ 6 -4.708 2.153 5.233 5.719 429.710 13

	Index: contrib/contribs/VariableR/trunk/ChangeLog
	===================================================================
	--- contrib/contribs/VariableR/trunk/ChangeLog (revision 907)
	+++ contrib/contribs/VariableR/trunk/ChangeLog (revision 908)
	@@ -1,80 +1,85 @@
	2016-03-09 <jthaler>
	VariableRPlugin.{hh,cc}:
	- Fixed -Wgnu-static-float-init warning. Small typos in comments fixed.
	+ Fixed -Wgnu-static-float-init warning. Small typos in comments fixed.
	+ Some commented out code removed in preparation for release.
	+ example.{cc,ref}:
	+ Change gen-kt example to p = -0.5 (since that is more physical)
	+ README:
	+ A bit more documentation added.
	2016-03-09 Gregory Soyez <soyez@fastjet.fr>
	VariableRPlugin.{hh,cc}:
	implemented as a gen-kt type of clustering. ClusterType is now
	defined as a double so one can pass any value with AKTLIKE,
	CALIKE and KTLIKE defined as -1, 0 and 1 respectively
	example.{cc,ref}: added an example using gen-kt(p=1/2)
	2016-03-09 <jthaler>
	AUTHORS,example.cc,VariableR.hh,VariableRPlugin.{cc,hh}:
	Updated authorship/copyright information
	2016-03-09 Gregory Soyez <soyez@fastjet.fr>
	VariableRPlugin.{hh,cc}:
	issued a warning for deprecated pre-clustering
	best strategy now has a transition between N2Plain and N2Tiled
	2016-03-08 <jthaler>
	example.cc,VariableR.hh,VariableRPlugin.{hh,cc}:
	Added fastjet authors to copyright statement
	AUTHORS: Added fastjet authors to the list
	NEWS, README: Fixed some typos
	VariableRPlugin.cc:
	Added "#" to timing test to avoid "make check" failing.
	Fixed error message to say that NNH is indeed supported in FJ<3.2.0
	example.cc:
	Adding in separate timing tests for the various strategies
	2016-03-08 Gregory Soyez <soyez@fastjet.fr>
	VariableRPlugin.{hh,cc}: adjusted a few details for compilation w FJ>3.2.0
	VariableRPlugin.{hh,cc}: added a "Best" strategy which selects a decent default
	2016-03-07 Gregory Soyez <soyez@fastjet.fr>
	example.{cc,ref}: prefixed the description line by "#" to avoid its
	inclusion in make check
	VariableRPlugin.{hh,cc}: For FastJet>=3.2.0, used the NNFJN2Tiled,
	NNFJN2Plain (and NNH) classes. The NNFJN2Tiled interface provides
	a much faster clustering than the native implementation.
	In practice, we introduced a "Strategy" class in VariableR that
	can take the following values: N2Tiled (the default if available),
	N2Plain, NNH or Native (the default for FastJet<3.2.0.
	VERSION: switched to 1.2.0-devel
	2014-06-02 <jthaler>
	Added dates to NEWS, wrapped to 80 characters.
	2014-05-30 <jthaler>
	Addressed suggestions of FJ authors (thanks!)
	Reorganized VariableRPlugin.hh for readability
	Wrapped README file to 80 characters
	Added Error("...")s for invalid parameter choices
	2014-05-18 <jthaler>
	Added VariableRPlugin.hh/.cc to stick to naming conventions.
	VariableR.hh maintained for backwards compatability
	Updated README to use new naming convention
	2014-05-17 <jthaler>
	Clarified NEWS about backwards compatibility.
	2014-05-16 <jthaler>
	Added additional comments for readability.
	First v1.1.0 release candidate.
	Added extra test in example.cc
	Changed the name VR -> VariableRPlugin
	Put classes JetDistancePair and CompareJetDistancePair into VariableRPlugin class.
	Renamed VRClustType -> ClusterType and put in VariableRPlugin class
	2014-04-14 <jthaler>
	Fix compile issue with example (needed stdio.h)
	Removed old CoreJetAlgorithm
	Changed CAVR and AKTVR to point to the new code.
	Added KTVR class.
	Removed example_new since this debugging is no longer needed.
	Added test with preclustering to example
	2014-03-31 <jthaler>
	Included David Krohn's new version of the code, which should have same behavior but faster.
	Old version of the code temporarily kept for debugging purposes.
	Moved print_jets out of VariableR.hh/.cc
	Added example_new file for debugging purposes
	Reindented the files
	2013-02-23 <jthaler>
	Fixed dependency issue (now using make depend)
	2013-02-22 <jthaler>
	Fixed memory management and failed make check issues.
	2013-02-21 <jthaler>
	First version submitted to fjcontrib
	2013-02-21 <jthaler>
	Initial creation based on previous plugin hosted at http://www.jthaler.net/jets/
	Index: contrib/contribs/VariableR/trunk/VariableRPlugin.cc
	===================================================================
	--- contrib/contribs/VariableR/trunk/VariableRPlugin.cc (revision 907)
	+++ contrib/contribs/VariableR/trunk/VariableRPlugin.cc (revision 908)
	@@ -1,478 +1,467 @@
	// VariableR Package
	// Questions/Comments? jthaler@jthaler.net
	//
	// Copyright (c) 2009-2016
	// David Krohn, Gregory Soyez, Jesse Thaler, and Lian-Tao Wang
	//
	// $Id: VariableR.cc 596 2014-04-16 22:15:15Z jthaler $
	//----------------------------------------------------------------------
	// 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 "VariableRPlugin.hh"

	#include <cstdio>
	#include "math.h"
	#include <iomanip>
	#include <cmath>
	#include <map>
	#include <sstream>
	#include <queue>

	#include <fastjet/NNH.hh>
	#if FASTJET_VERSION_NUMBER >= 30200
	#include <fastjet/NNFJN2Plain.hh>
	#include <fastjet/NNFJN2Tiled.hh>
	#endif

	FASTJET_BEGIN_NAMESPACE // defined in fastjet/internal/base.hh

	namespace contrib {

	//----------------------------------------------------------------------
	- // Defining static const.
	+ // Defining static constants
	const double VariableRPlugin::CALIKE = 0.0;
	const double VariableRPlugin::KTLIKE = 1.0;
	const double VariableRPlugin::AKTLIKE = -1.0;

	//----------------------------------------------------------------------
	// classes to help run a Variable R algorithm using NN-type classes

	// class carrying particle-independent information
	class VariableRNNInfo {
	public:
	VariableRNNInfo(double rho2_in, double min_r2_in, double max_r2_in,
	VariableRPlugin::ClusterType clust_type_in)
	: _rho2(rho2_in), _min_r2(min_r2_in), _max_r2(max_r2_in),
	_clust_type(clust_type_in) {}

	double rho2() const {return _rho2; }
	double min_r2() const {return _min_r2;}
	double max_r2() const {return _max_r2;}
	double momentum_scale_of_pt2(double pt2) const {
	+
	+ // The commented out code below is only marginally faster and
	+ // doesn't always play nice with numerical precision
	// if (_clust_type == VariableRPlugin::AKTLIKE){ return 1.0/pt2;}
	// else if (_clust_type == VariableRPlugin::KTLIKE){ return pt2;}
	// else { return 1.0;} // VariableRPlugin::CALIKE
	+
	return pow(pt2,_clust_type);
	}

	private:
	double _rho2; ///< constant that controls the overall R magnitude
	double _min_r2; ///< minimal allowed radius
	double _max_r2; ///< maximal allowed radius
	VariableRPlugin::ClusterType _clust_type;
	};

	// class carrying the minimal info required for the clustering
	class VariableRBriefJet {
	public:
	void init(const PseudoJet & jet, VariableRNNInfo *info) {
	_rap = jet.rap();
	_phi = jet.phi();

	double pt2 = jet.pt2();

	// get the appropriate "radius"
	_beam_R2 = info->rho2()/pt2;
	if (_beam_R2 > info->max_r2()){ _beam_R2 = info->max_r2();}
	else if (_beam_R2 < info->min_r2()){ _beam_R2 = info->min_r2();}

	// get the appropriate momentum scale
	_mom_factor2 = info->momentum_scale_of_pt2(pt2);
	}

	double geometrical_distance(const VariableRBriefJet * jet) const {
	double dphi = std::abs(_phi - jet->_phi);
	double deta = (_rap - jet->_rap);
	if (dphi > pi) {dphi = twopi - dphi;}
	return dphidphi + detadeta;
	}

	double geometrical_beam_distance() const { return _beam_R2; }

	double momentum_factor() const{ return _mom_factor2; }

	/// make this BJ class compatible with the use of NNH
	double distance(const VariableRBriefJet * other_bj_jet){
	double mom1 = momentum_factor();
	double mom2 = other_bj_jet->momentum_factor();
	return (mom1<mom2 ? mom1 : mom2) * geometrical_distance(other_bj_jet);
	}
	double beam_distance(){
	return momentum_factor() * geometrical_beam_distance();
	}

	// the following are required by N2Tiled
	inline double rap() const{ return _rap;}
	inline double phi() const{ return _phi;}

	private:
	double _rap, _phi, _mom_factor2, _beam_R2;
	};

	//----------------------------------------------------------------------
	// now the implementation of VariableR itself
	//----------------------------------------------------------------------
	LimitedWarning VariableRPlugin::_preclustering_deprecated_warning;

	// Constructor that sets VR algorithm parameters
	// - rho mass scale for effective radius (i.e. R ~ rho/pT)
	// - min_r minimum jet radius
	// - max_r maximum jet radius
	// - clust_type whether to use CA-like, kT-like, or anti-kT-like distance measure
	// - strategy one of Best (the default), N2Tiled , N2Plain or NNH
	VariableRPlugin::VariableRPlugin(double rho, double min_r, double max_r,
	ClusterType clust_type, bool precluster,
	Strategy requested_strategy)
	:_rho2(rhorho), _min_r2(min_rmin_r), _max_r(max_r), _max_r2(max_r*max_r),
	_clust_type(clust_type), _requested_strategy(requested_strategy),
	_precluster(precluster), _pre_jet_def(kt_algorithm, min_r){
	// Added errors for user input.
	if (min_r < 0.0) throw Error("VariableRPlugin: Minimum radius must be positive.");
	if (precluster && min_r==0.0) throw Error("VariableRPlugin: To apply preclustering, minimum radius must be non-zero.");
	if (max_r < 0.0) throw Error("VariableRPlugin: Maximum radius must be positive.");
	if (min_r > max_r) throw Error("VariableRPlugin: Minimum radius must be bigger than or equal to maximum radius.");

	// decide the strategy

	#if FASTJET_VERSION_NUMBER < 30200
	// this is only supported for the Best and Native strategies
	if ((requested_strategy!=Best) && (requested_strategy!=Native) && (requested_strategy!=NNH))
	throw Error("VariableRPlugin: with FastJet<3.2.0, Native, Best, and NNH are the only supported strategies");
	#endif

	// with pre-clustering we're forced into the native clustering
	// Q: do we issue a warning? throw?
	if (precluster){
	// this is only supported for the Best and Native strategies
	if ((requested_strategy!=Best) && (requested_strategy!=Native))
	throw Error("VariableRPlugin: pre-clustering is only supported for the Native and Best strategies");

	- _preclustering_deprecated_warning.warn("VariableRPlugin: internal pre-clustering is deprecated, use the NestedDefs FastJet plugin instead.");
	+ _preclustering_deprecated_warning.warn("VariableRPlugin: internal pre-clustering is deprecated; use the NestedDefs FastJet plugin instead.");
	}
	}

	//----------------------------------------------------------------------
	// Description of algorithm, including parameters
	string VariableRPlugin::description() const{
	stringstream myStream("");

	myStream << "Variable R (0903.0392), ";

	if (_clust_type == AKTLIKE){
	myStream << "AKT";
	} else if (_clust_type == CALIKE){
	myStream << "CA";
	} else if (_clust_type == KTLIKE){
	myStream << "KT";
	} else {
	myStream << "GenKT(p=" << _clust_type << ")";
	}

	myStream << fixed << setprecision(1) << ", rho=" << sqrt(_rho2);
	myStream << ", min_r=" << sqrt(_min_r2);
	myStream << ", max_r=" << sqrt(_max_r2);
	myStream << (_precluster ? ", with precluster" : "");
	switch (_requested_strategy){
	case Best: myStream << ", strategy=Best"; break;
	case N2Tiled: myStream << ", strategy=N2Tiled"; break;
	case N2Plain: myStream << ", strategy=N2Plain"; break;
	case NNH: myStream << ", strategy=NNH"; break;
	case Native: myStream << ", strategy=Native"; break;
	};

	return myStream.str();
	}

	//----------------------------------------------------------------------
	// do the clustering
	void VariableRPlugin::run_clustering(ClusterSequence & cs) const {
	Strategy strategy = _requested_strategy;

	// decide the best option upon request
	if (_requested_strategy==Best){
	strategy = _best_strategy(cs.jets().size());
	}

	// handle the case of native clustering
	if (strategy==Native){
	_native_clustering(cs);
	return;
	}

	// handle the NN-type clustering
	VariableRNNInfo info(_rho2, _min_r2, _max_r2, _clust_type);

	#if FASTJET_VERSION_NUMBER >= 30200
	if (strategy==N2Tiled){
	NNFJN2Tiled<VariableRBriefJet,VariableRNNInfo> nnt(cs.jets(), _max_r, &info);
	_NN_clustering(cs, nnt);
	} else if (strategy==N2Plain){
	NNFJN2Plain<VariableRBriefJet,VariableRNNInfo> nnp(cs.jets(), &info);
	_NN_clustering(cs, nnp);
	} else { // NNH is the only option left
	#endif
	fastjet::NNH<VariableRBriefJet,VariableRNNInfo> nnh(cs.jets(), &info);
	_NN_clustering(cs, nnh);
	#if FASTJET_VERSION_NUMBER >= 30200
	}
	#endif
	}

	//---------------------------------------------------------------------
	// decide the optimal strategy
	VariableRPlugin::Strategy VariableRPlugin::_best_strategy(unsigned int N) const{
	#if FASTJET_VERSION_NUMBER >= 30200
	// pre-clustering requires the native implementation
	if (_precluster) return Native;
	// use the FastJet (v>/3.1) transition between N2Plain and N2Tiled
	if (N <= 30 \|\| N <= 39.0/(max(_max_r, 0.1) + 0.6)) return N2Plain;
	return N2Tiled;
	#else
	return Native;
	#endif
	}

	//---------------------------------------------------------------------
	// the "new" NN-style clustering
	template<typename NN>
	void VariableRPlugin::_NN_clustering(ClusterSequence &cs, NN &nn) const{
	int njets = cs.jets().size();
	while (njets > 0) {
	int i, j, k;
	double dij = nn.dij_min(i, j);

	if (j >= 0) {
	cs.plugin_record_ij_recombination(i, j, dij, k);
	nn.merge_jets(i, j, cs.jets()[k], k);
	} else {
	cs.plugin_record_iB_recombination(i, dij);
	nn.remove_jet(i);
	}
	njets--;
	}
	}



	//----------------------------------------------------------------------
	- // the code below is for the "native" clustering
	+ // the code below is for the "native" clustering only

	// precluster into kT subjets if desired.
	void VariableRPlugin::_preclustering(ClusterSequence & cs, set<int>& unmerged_jets) const {
	int cntr(0);
	for(vector<PseudoJet>::const_iterator it = cs.jets().begin(); it != cs.jets().end(); it++)
	unmerged_jets.insert(unmerged_jets.end(), cntr++);

	// Make preclusters
	ClusterSequence pre_cs(cs.jets(), _pre_jet_def);
	vector<PseudoJet> preclustered_jets = pre_cs.inclusive_jets();
	vector<int> particle_jet_indices = pre_cs.particle_jet_indices(preclustered_jets);

	// This code take preclustered objects and puts them into the ClusterSequence tree of VR
	// TODO: Figure out if there is a better way to do this step.
	for(int i = 0 ; i < (int)preclustered_jets.size(); i++){
	queue<int> constit_indices;
	for(int j = 0 ; j < (int)particle_jet_indices.size(); j++)
	if(particle_jet_indices[j] == i)
	constit_indices.push(j);

	int final_jet;
	while(constit_indices.size() > 1){
	int indx1 = constit_indices.front();
	unmerged_jets.erase(indx1);
	constit_indices.pop();
	int indx2 = constit_indices.front();
	unmerged_jets.erase(indx2);
	constit_indices.pop();
	cs.plugin_record_ij_recombination(indx1, indx2, 0., final_jet);
	constit_indices.push(final_jet);
	unmerged_jets.insert(unmerged_jets.end(), final_jet);
	}
	}
	}

	// Implements VR alorithm (virtual function from JetDefinition::Plugin)
	void VariableRPlugin::_native_clustering(ClusterSequence & cs) const {
	set<int> unmerged_jets;

	if(_precluster){ // do kT preclustering
	assert(_min_r2 > 0.); // since this is (min_r)^2, this should never happen
	_preclustering(cs, unmerged_jets);
	} else // make a list of the unmerged jets
	for(int i = 0 ; i < (int)cs.jets().size(); i++)
	unmerged_jets.insert(unmerged_jets.end(), i);

	// priority_queue is part of the C++ standard library.
	// The objects being sorted are JetDistancePairs
	// The comparison function is just asking who has the smallest distance
	// The distances are set initially by _setup_distance_measures and then updated by the while loop below.
	vector<JetDistancePair> jet_vec;
	_setup_distance_measures(cs, jet_vec, unmerged_jets);
	priority_queue< JetDistancePair, vector<JetDistancePair>, CompareJetDistancePair > jet_queue(jet_vec.begin(),jet_vec.end());

	// go through the jet_queue until empty
	while(!jet_queue.empty()){

	// find the closest pair
	JetDistancePair jdpair = jet_queue.top();
	jet_queue.pop();

	// Rebuild the jet_queue
	// DK - the 1.5 below was just found empirically
	// JDT - this code is safe, but I'm not 100% sure why it is needed.
	// It rebuilds the jet_queue instead of letting the below functions do the hard work.
	if(jet_queue.size() > 50 && jet_queue.size() > 1.5unmerged_jets.size()unmerged_jets.size()){
	jet_vec.clear();
	_setup_distance_measures(cs, jet_vec, unmerged_jets);
	priority_queue< JetDistancePair, vector<JetDistancePair>, CompareJetDistancePair > tmp_jet_queue(jet_vec.begin(),jet_vec.end());
	swap(jet_queue,tmp_jet_queue);
	}

	// make sure not merged
	if((unmerged_jets.find(jdpair.j1) == unmerged_jets.end()) \|\| (jdpair.j2 != -1 && unmerged_jets.find(jdpair.j2) == unmerged_jets.end()))
	continue;


	if(jdpair.j2 == -1) // If closest distance is to beam, then merge with beam
	_merge_jet_with_beam(cs, jdpair, unmerged_jets);
	else // Otherwise, merge jets back together
	_merge_jets(cs, jdpair, jet_queue, unmerged_jets);
	}
	}


	// Add final jet to clust_seq. No need to update jet_queue, since this jet has already been deleted.
	void VariableRPlugin::_merge_jet_with_beam(ClusterSequence & clust_seq, JetDistancePair & jdp, set<int>& unmerged_jets) const{
	clust_seq.plugin_record_iB_recombination(jdp.j1, jdp.distance);
	unmerged_jets.erase(jdp.j1);
	}


	// Add jet merging to clust_seq. Here, we need to update the priority_queue since some of the measures have changes.
	void VariableRPlugin::_merge_jets(ClusterSequence & clust_seq,
	JetDistancePair & jdp,
	priority_queue< JetDistancePair, vector<JetDistancePair>, CompareJetDistancePair > &jet_queue,
	set<int>& unmerged_jets) const{


	int new_jet_num;
	clust_seq.plugin_record_ij_recombination(jdp.j1,jdp.j2,jdp.distance,new_jet_num);

	unmerged_jets.erase(jdp.j1);
	unmerged_jets.erase(jdp.j2);

	// take the resulting jet and recompute all distances with it
	for(set<int>::iterator it = unmerged_jets.begin(); it != unmerged_jets.end(); it++){
	JetDistancePair jpair;
	jpair.j1 = new_jet_num;
	jpair.j2 = (*it);
	jpair.distance = _get_JJ_distance_measure(clust_seq.jets()[*it],clust_seq.jets()[new_jet_num]);
	jet_queue.push(jpair);
	}
	unmerged_jets.insert(unmerged_jets.end(), new_jet_num);

	// also add the new distance to beam
	JetDistancePair jpair;
	jpair.j1 = new_jet_num;
	jpair.j2 = -1; // -1 is for the beam
	jpair.distance = _get_JB_distance_measure(clust_seq.jets()[new_jet_num]);
	jet_queue.push(jpair);
	}

	// Initial distance setup
	void VariableRPlugin::_setup_distance_measures(ClusterSequence & clust_seq,
	vector<JetDistancePair> &jet_vec,
	set<int> & unmerged_jets) const{
	JetDistancePair jpair;

	// Add the jet-jet distances
	for(set<int>::iterator it1 = unmerged_jets.begin(); it1 != unmerged_jets.end(); it1++){
	for(set<int>::iterator it2 = it1; it2 != unmerged_jets.end(); it2++){
	if((it1) != (it2)){
	jpair.j1 = (*it1);
	jpair.j2 = (*it2);
	jpair.distance = _get_JJ_distance_measure(clust_seq.jets()[it1],clust_seq.jets()[it2]);
	jet_vec.push_back(jpair);
	}
	}
	// Add the jet-beam distances, and set initial merge info
	jpair.j1 = (*it1);
	jpair.j2 = -1; // -1 is for the beam
	jpair.distance = _get_JB_distance_measure(clust_seq.jets()[*it1]);
	jet_vec.push_back(jpair);
	}
	}

	// get the dij between two jets
	- // Different measures for AKTLIKE, CALIKE, and KTLIKE
	- // TODO: add arbitrary p weighting
	+ // Different measures for AKTLIKE, CALIKE, and KTLIKE, as well as for GenKT
	double VariableRPlugin::_get_JJ_distance_measure(const PseudoJet& j1, const PseudoJet& j2) const {
	- // Old code with ClusterType as an enum
	- // double ret;
	- // switch(_clust_type){
	- // case AKTLIKE:
	- // ret = min(1./j1.perp2(), 1./j2.perp2());
	- // break;
	- // case CALIKE :
	- // ret = 1.;
	- // break;
	- // case KTLIKE:
	- // ret = min(j1.perp2(), j2.perp2());
	- // break;
	- // default:
	- // assert(false);
	- // }

	- // new code with generic p (note that we avoid having to take the
	- // "pow" twice)
	+ // new code since version 1.2 using generic p for GenKT
	+ // Note that this is written to avoid having to take the "pow" twice.
	// We also keep the +-1 cases separate (for the same reason)
	double ret;
	if (_clust_type == AKTLIKE){
	ret = min(1.0/j1.perp2(), 1.0/j2.perp2());
	} else if (_clust_type == CALIKE){
	ret = 1.0;
	} else if (_clust_type == KTLIKE){
	ret = min(j1.perp2(), j2.perp2());
	} else if (_clust_type>=0){ // (include = just to be sure)
	ret = pow(min(j1.perp2(), j2.perp2()), _clust_type);
	} else {
	ret = pow(min(1.0/j1.perp2(), 1.0/j2.perp2()), -_clust_type);
	}

	ret *= j1.squared_distance(j2);
	return ret;
	}

	// jet diB between jet and beam
	- // Different measures for AKTLIKE, CALIKE, and KTLIKE
	- // TODO: add arbitrary p weighting
	+ // Different measures for AKTLIKE, CALIKE, and KTLIKE, as well as for GenKT
	double VariableRPlugin::_get_JB_distance_measure(const PseudoJet& jet) const{
	double pre_factor = pow(jet.perp2(), _clust_type);
	double geom_factor = _rho2 / jet.perp2();

	// Q: do we want to implement the simplification of the pt2
	// factors explicitly for KTLIKE?
	+ // A: Since the "Native" code is being replaced eventually
	+ // (and since KTVR is not recommended), I don't think this is necessary.
	if(geom_factor < _min_r2) return _min_r2 * pre_factor;
	if(geom_factor > _max_r2) return _max_r2 * pre_factor;
	return geom_factor * pre_factor;
	}



	} // namespace contrib

	FASTJET_END_NAMESPACE
	Index: contrib/contribs/VariableR/trunk/NEWS
	===================================================================
	--- contrib/contribs/VariableR/trunk/NEWS (revision 907)
	+++ contrib/contribs/VariableR/trunk/NEWS (revision 908)
	@@ -1,30 +1,30 @@
	1.2.0 (2016-03-XX): Added an implementation based on the NNFJN2Tiled
	(and similar) classes in FastJet, bringing
	sizable speed improvements. It is introduced
	through a "Strategy" argument and requires
	FastJet>=3.2.0 ["Native" clustering used for older
	- versions or to recover the previous behaviour].
	+ versions or to recover the previous behavior].
	Please use fastjet::NestedDefsPlugin whenever
	pre-clustering is needed.
	1.1.1 (2014-06-02): Added VariableRPlugin.cc/.hh to conform to recommended
	FastJet naming structure. Old VariableR.hh header preserved
	for backwards compatibility.
	1.1.0 (2014-04-18): New version incorporates optimized code from David Krohn,
	which should be faster than the old version. The AKTVR and
	CAVR interfaces are the same as v1.0 (with identical
	results). An additional VariableRPlugin interface is
	available, which is now recommended. This new interface
	offers the possibility of setting a min_r > 0 as well as
	using kT preclustering for speed gains.

	1.0.1 (2013-02-23): Fixed memory management and potential dependency issue
	(thanks FJ authors)
	1.0.0 (2013-02-22): New release using FastJet contrib framework. All of the
	code has been merged into the VariableR.hh and VariableR.cc
	files (whereas previously there were separate files for AKT,
	CA, and the CoreJetAlgorithm classes). The example program
	runs both AKTVR and CAVR.

	Previous version available from:
	http://jthaler.net/jets/VRPlugins.tar.gz

	\ No newline at end of file
	Index: contrib/contribs/VariableR/trunk/README
	===================================================================
	--- contrib/contribs/VariableR/trunk/README (revision 907)
	+++ contrib/contribs/VariableR/trunk/README (revision 908)
	@@ -1,68 +1,71 @@
	The VariableR package is based on the physics described in:

	Jets with Variable R.
	David Krohn, Jesse Thaler, and Lian-Tao Wang.
	JHEP 0906:059 (2009), arXiv:0903.0392.

	---

	Since version 1.2, one uses the following constructor for creating a
	VariableRPlugin (see the VariableRplugin.hh header file for more
	details)

	VariableRPlugin(double rho, double min_r, double max_r,
	ClustType clust_type,
	bool precluster=false,
	Strategy strategy=VARIABLE_R_DEFAULT_STRATEGY);

	with

	rho : mass parameter of VR algorithm, where the effective radius is R~rho/pT.
	min_r: minimum allowed jet radius
	max_r: maximum allowed jet radius

	- clust_type: AKTLIKE, CALIKE, KTLIKE, for determining the cluster ordering
	- Recommended usage is AKTLIKE.
	+ clust_type: AKTLIKE, CALIKE, KTLIKE, for determining the cluster ordering.
	+ For generalized-kt-like behavior, you can enter an arbitrary double
	+ which is matched to the p value of generalized kt.
	+ Recommended usage is AKTLIKE = -1.0
	precluster: true to build small kT jets of size min_r before applying VR
	strategy : specify the internal clustering strategy. NNFJN2Tiled, the default
	for FastJet>=3.2.0 is expected to be faster. Other options are
	NNFJN2Plain, NNH, and Native. For FastJet<3.2.0, or when
	pre-clustering is requested, the Native strategy is used by default.

	+
	The VariableRPlugin can be used via the standard FastJet plugin mechanism.

	Note that pre-clustering is deprecated: one should rather use the
	FastJet NestedDefs plugin along the following lines:

	#include <fastjet/NestedDefs.hh>
	#include <fastjet/contrib/VariableRPlugin.hh>

	std::list<JetDefinition> jet_defs;
	jet_defs.push_back(JetDefinition(kt_algorithm, min_r));
	contrib::VariableRPlugin vr_plugin(rho,min_r,max_r,clust_type);
	jet_defs.push_back(JetDefinition(&vr_plugin));

	NestedDefsPlugin nd_plugin(jet_defs);
	JetDefinition jet_def(&nd_plugin);

	For backwards compatibility with Version 1.0, we define the AKTVR and CAVR
	classes with fewer options (plus a KTVR for consistency). These are available
	in the VariableR.hh header.

	AKTVR(double rho, double max_r): variable-R anti-kT algorithm
	CAVR(double rho, double max_r): variable-R Cambridge-Aachen algorithm
	KTVR(double rho, double max_r): variable-R kT algorithm

	---

	Note:
	- Version 1.2 is based on the NNFJN2Tiled and NNFJN2Plain pairwise
	clustering classes available in FastJet>=3.2.0. For earlier versions
	of FastJet (or when pre-clustering is requested), one reverts back to
	David Krohn's implementation ("Native" strategy)
	- Version 1.1 is based on a new clustering code written by David
	Krohn, which should be faster than the previous version.
	- Version 1.0 had a CoreJetAlgorithm, which allowed the definition of
	generic pairwise recombination scheme. This functionality is now
	available in the FastJet NNH class, so CoreJetAlgorithm is no longer
	needed.
	\ No newline at end of file

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 10:17 AM (1 d, 17 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5111084
Default Alt Text: (63 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions