AxesFinder.hh
No OneTemporary
Actions

Size

15 KB

Subscribers

None

AxesFinder.hh
View Options

	// Nsubjettiness Package
	// Questions/Comments? jthaler@jthaler.net
	//
	// Copyright (c) 2011-14
	// Jesse Thaler, Ken Van Tilburg, Christopher K. Vermilion, and TJ Wilkason
	//
	//----------------------------------------------------------------------
	// 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_AXESFINDER_HH__
	#define __FASTJET_CONTRIB_AXESFINDER_HH__

	#include "WinnerTakeAllRecombiner.hh"
	#include "MeasureFunction.hh"

	#include "fastjet/PseudoJet.hh"
	#include "fastjet/ClusterSequence.hh"
	#include "fastjet/JetDefinition.hh"

	#include <cmath>
	#include <vector>
	#include <list>

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

	namespace contrib{

	///////
	//
	// Axes Finder Options
	//
	///////

	//------------------------------------------------------------------------
	/// \class AxesFinder
	// This is the base class for all axes finders. These axes are used along with the MeasureFunctions to calculate
	// tau_N. There are different implementations of axes finding that are defined in derived classes below.
	class AxesFinder {

	protected:
	AxesFinder() {}

	public:
	virtual ~AxesFinder(){}

	virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector<fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) = 0;

	};

	//------------------------------------------------------------------------
	/// \class AxesFinderFromExclusiveJetDefinition
	// This class finds axes by clustering the particles and then finding the exclusive jets. This can be implemented
	// with different jet algorithms.
	class AxesFinderFromExclusiveJetDefinition : public AxesFinder {

	private:
	fastjet::JetDefinition _def;

	public:
	AxesFinderFromExclusiveJetDefinition(fastjet::JetDefinition def) : _def(def) {}

	virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
	fastjet::ClusterSequence jet_clust_seq(inputs, _def);
	return jet_clust_seq.exclusive_jets(n_jets);
	}
	};

	//------------------------------------------------------------------------
	/// \class AxesFinderFromWTA_KT
	// This class finds axes by finding the exlusive jets after clustering according to a kT algorithm and a
	// winner take all recombination scheme.
	class AxesFinderFromWTA_KT : public AxesFinderFromExclusiveJetDefinition {
	private:
	const WinnerTakeAllRecombiner *recomb;
	public:
	AxesFinderFromWTA_KT() : AxesFinderFromExclusiveJetDefinition(
	fastjet::JetDefinition(fastjet::kt_algorithm,
	fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
	recomb = new WinnerTakeAllRecombiner(),
	fastjet::Best)) {}
	~AxesFinderFromWTA_KT() {delete recomb;}
	};

	//------------------------------------------------------------------------
	/// \class AxesFinderFromWTA_CA
	// This class finds axes by finding the exlusive jets after clustering according to a CA algorithm and a
	// winner take all recombination scheme.
	class AxesFinderFromWTA_CA : public AxesFinderFromExclusiveJetDefinition {
	private:
	const WinnerTakeAllRecombiner *recomb;
	public:
	AxesFinderFromWTA_CA() : AxesFinderFromExclusiveJetDefinition(
	fastjet::JetDefinition(fastjet::cambridge_algorithm,
	fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
	recomb = new WinnerTakeAllRecombiner(),
	fastjet::Best)) {}
	~AxesFinderFromWTA_CA() {delete recomb;}
	};

	// class added by TJW -- 1/27
	//------------------------------------------------------------------------
	/// \class AxesFinderFromWTA2_KT
	// This class finds axes by finding the exlusive jets after clustering according to a kT algorithm and a
	// winner take all recombination scheme with alpha = 2.
	class AxesFinderFromWTA2_KT : public AxesFinderFromExclusiveJetDefinition {
	private:
	const WinnerTakeAllRecombiner *recomb;
	public:
	AxesFinderFromWTA2_KT() : AxesFinderFromExclusiveJetDefinition(
	fastjet::JetDefinition(fastjet::kt_algorithm,
	fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
	recomb = new WinnerTakeAllRecombiner(2), // uses alpha = 2 here -- TJW 1/27
	fastjet::Best)) {}
	~AxesFinderFromWTA2_KT() {delete recomb;}
	};

	// class added by TJW -- 1/27
	//------------------------------------------------------------------------
	/// \class AxesFinderFromWTA2_CA
	// This class finds axes by finding the exlusive jets after clustering according to a CA algorithm and a
	// winner take all recombination scheme with alpha = 2.
	class AxesFinderFromWTA2_CA : public AxesFinderFromExclusiveJetDefinition {
	private:
	const WinnerTakeAllRecombiner *recomb;
	public:
	AxesFinderFromWTA2_CA() : AxesFinderFromExclusiveJetDefinition(
	fastjet::JetDefinition(fastjet::cambridge_algorithm,
	fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
	recomb = new WinnerTakeAllRecombiner(2), //uses alpha = 2 here -- TJW 1/27
	fastjet::Best)) {}
	~AxesFinderFromWTA2_CA() {delete recomb;}
	};

	//------------------------------------------------------------------------
	/// \class AxesFinderFromKT
	// This class finds axes by finding the exlusive jets after clustering according to a kT algorithm and a
	// E_scheme recombination.
	class AxesFinderFromKT : public AxesFinderFromExclusiveJetDefinition {
	public:
	AxesFinderFromKT() : AxesFinderFromExclusiveJetDefinition(
	fastjet::JetDefinition(fastjet::kt_algorithm,
	fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
	fastjet::E_scheme,
	fastjet::Best)) {}
	};

	//------------------------------------------------------------------------
	/// \class AxesFinderFromCA
	// This class finds axes by finding the exlusive jets after clustering according to a CA algorithm and a
	// E_scheme recombination.
	class AxesFinderFromCA : public AxesFinderFromExclusiveJetDefinition {
	public:
	AxesFinderFromCA() : AxesFinderFromExclusiveJetDefinition(
	fastjet::JetDefinition(fastjet::cambridge_algorithm,
	fastjet::JetDefinition::max_allowable_R, //maximum jet radius constant
	fastjet::E_scheme,
	fastjet::Best)) {}
	};


	//------------------------------------------------------------------------
	/// \class AxesFinderFromHardestJetDefinition
	// This class finds axes by clustering the particles and then finding the n hardest inclusive jets.
	// This can be implemented with different jet algorithms.
	class AxesFinderFromHardestJetDefinition : public AxesFinder {

	private:
	fastjet::JetDefinition _def;

	public:
	AxesFinderFromHardestJetDefinition(fastjet::JetDefinition def) : _def(def) {}

	virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
	fastjet::ClusterSequence jet_clust_seq(inputs, _def);
	std::vector<fastjet::PseudoJet> myJets = sorted_by_pt(jet_clust_seq.inclusive_jets());
	myJets.resize(n_jets); // only keep n hardest
	return myJets;
	}
	};

	//------------------------------------------------------------------------
	/// \class AxesFinderFromAntiKT
	// This class finds axes by finding the n hardest jets after clustering the particles according
	// to an anti kT algorithm and E_scheme.
	class AxesFinderFromAntiKT : public AxesFinderFromHardestJetDefinition {
	public:
	AxesFinderFromAntiKT(double R0) : AxesFinderFromHardestJetDefinition(fastjet::JetDefinition(fastjet::antikt_algorithm,R0,fastjet::E_scheme,fastjet::Best)) {}
	};


	//------------------------------------------------------------------------
	/// \class AxesFinderFromUserInput
	// This class allows the user to manually define the axes.
	class AxesFinderFromUserInput : public AxesFinder {

	public:
	AxesFinderFromUserInput() {}
	~AxesFinderFromUserInput() {}

	virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputs, const std::vector<fastjet::PseudoJet>& currentAxes) {
	assert(currentAxes.size() == (unsigned int) n_jets);
	return currentAxes;
	}
	};

	//This is a helper class for the Minimum Axes Finders. It is defined later.
	class LightLikeAxis;


	//------------------------------------------------------------------------
	/// \class AxesFinderFromOnePassMinimization
	// This class defines an AxesFinder that uses Kmeans minimization, but only on a single pass.
	class AxesFinderFromOnePassMinimization : public AxesFinder {

	private:
	AxesFinder* _startingFinder;
	double _precision; // Desired precision in axes alignment
	int _halt; // maximum number of steps per iteration

	double _beta;
	double _Rcutoff;

	DefaultUnnormalizedMeasure _measureFunction;

	public:

	// From a startingFinder, try to minimize the unnormalized_measure
	AxesFinderFromOnePassMinimization(AxesFinder* startingFinder, double beta, double Rcutoff)
	: _startingFinder(startingFinder),
	_precision(0.0001), //hard coded for now
	_halt(1000), //hard coded for now
	_beta(beta),
	_Rcutoff(Rcutoff),
	_measureFunction(beta, Rcutoff)
	{}

	~AxesFinderFromOnePassMinimization() {
	delete _startingFinder; //TODO: Convert to smart pointers to avoid this.
	}

	virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& currentAxes);

	template <int N> std::vector<LightLikeAxis> UpdateAxesFast(const std::vector <LightLikeAxis> & old_axes,
	const std::vector <fastjet::PseudoJet> & inputJets);

	std::vector<LightLikeAxis> UpdateAxes(const std::vector <LightLikeAxis> & old_axes,
	const std::vector <fastjet::PseudoJet> & inputJets);

	};


	//------------------------------------------------------------------------
	/// \class AxesFinderFromKmeansMinimization
	// This class finds finds axes by using Kmeans clustering to minimizaiton N-jettiness. Given a first set of
	// starting axes, it updates n times to get as close to the global minimum as possible. This class calls OnePass many times,
	// added noise to the axes.
	class AxesFinderFromKmeansMinimization : public AxesFinder{

	private:
	int _n_iterations; // Number of iterations to run (0 for no minimization, 1 for one-pass, >>1 for global minimum)
	double _noise_range; // noise range for random initialization

	DefaultUnnormalizedMeasure _measureFunction; //function to test whether minimum is reached
	AxesFinderFromOnePassMinimization _onePassFinder; //one pass finder for minimization

	PseudoJet jiggle(const PseudoJet& axis);

	public:
	AxesFinderFromKmeansMinimization(AxesFinder *startingFinder, double beta, double Rcutoff, int n_iterations) :
	_n_iterations(n_iterations),
	_noise_range(1.0), // hard coded for the time being
	_measureFunction(beta, Rcutoff),
	_onePassFinder(startingFinder, beta, Rcutoff)
	{}

	virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & inputJets, const std::vector<fastjet::PseudoJet>& currentAxes);

	};

	//------------------------------------------------------------------------
	/// \class AxesFinderFromGeometricMinimization
	// This class finds axes by minimizing the Lorentz dot product distance between axes and particles. Given a first set of starting axes,
	// it essentially does stable cone finxing.
	class AxesFinderFromGeometricMinimization : public AxesFinder {

	private:
	AxesFinder* _startingFinder;
	MeasureFunction* _function;
	double _Rcutoff;
	double _nAttempts;
	double _accuracy;


	public:
	AxesFinderFromGeometricMinimization(AxesFinder* startingFinder, double Rcutoff) : _startingFinder(startingFinder), _Rcutoff(Rcutoff) {
	_nAttempts = 100;
	_accuracy = 0.000000001;
	_function = new GeometricMeasure(_Rcutoff);
	}

	~AxesFinderFromGeometricMinimization() {
	delete _startingFinder; //TODO: Convert to smart pointers to avoid this.
	delete _function;
	}
	virtual std::vector<fastjet::PseudoJet> getAxes(int n_jets, const std::vector <fastjet::PseudoJet> & particles, const std::vector<fastjet::PseudoJet>& currentAxes);
	};

	//------------------------------------------------------------------------
	/// \class LightLikeAxis
	// This is a helper class for the minimum Axes Finders classes above. It creates a convenient way of defining axes
	// in order to better facilitate calculations.
	class LightLikeAxis {
	private:
	double _rap, _phi, _weight, _mom;

	double DistanceSq(double rap2, double phi2) const {
	double rap1 = _rap;
	double phi1 = _phi;

	double distRap = rap1-rap2;
	double distPhi = std::fabs(phi1-phi2);
	if (distPhi > M_PI) {distPhi = 2.0*M_PI - distPhi;}
	return sq(distRap) + sq(distPhi);
	}

	double Distance(double rap2, double phi2) const {
	return std::sqrt(DistanceSq(rap2,phi2));
	}


	public:
	LightLikeAxis() : _rap(0.0), _phi(0.0), _weight(0.0), _mom(0.0) {}
	LightLikeAxis(double my_rap, double my_phi, double my_weight, double my_mom) :
	_rap(my_rap), _phi(my_phi), _weight(my_weight), _mom(my_mom) {}
	double rap() const {return _rap;}
	double phi() const {return _phi;}
	double weight() const {return _weight;}
	double mom() const {return _mom;}
	void set_rap(double my_set_rap) {_rap = my_set_rap;}
	void set_phi(double my_set_phi) {_phi = my_set_phi;}
	void set_weight(double my_set_weight) {_weight = my_set_weight;}
	void set_mom(double my_set_mom) {_mom = my_set_mom;}
	void reset(double my_rap, double my_phi, double my_weight, double my_mom) {_rap=my_rap; _phi=my_phi; _weight=my_weight; _mom=my_mom;}

	// return PseudoJet with information
	fastjet::PseudoJet ConvertToPseudoJet();

	double DistanceSq(const fastjet::PseudoJet& input) const {
	return DistanceSq(input.rap(),input.phi());
	}

	double Distance(const fastjet::PseudoJet& input) const {
	return std::sqrt(DistanceSq(input));
	}

	double DistanceSq(const LightLikeAxis& input) const {
	return DistanceSq(input.rap(),input.phi());
	}

	double Distance(const LightLikeAxis& input) const {
	return std::sqrt(DistanceSq(input));
	}

	};

	} //namespace contrib

	FASTJET_END_NAMESPACE

	#endif // __FASTJET_CONTRIB_AXESFINDER_HH__

File Metadata

Mime Type: text/x-c++
Expires: Fri, Apr 4, 9:37 PM (1 d, 4 h)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 4737421
Default Alt Text: AxesFinder.hh (15 KB)

AxesFinder.hhNo OneTemporaryActions

AxesFinder.hhView Options

File Metadata

Event Timeline

AxesFinder.hh
No OneTemporary
Actions

AxesFinder.hh
View Options