No OneTemporary
Actions

Size

74 KB

Subscribers

None

View Options

	Index: contrib/contribs/JetsWithoutJets/trunk/example.ref
	===================================================================
	--- contrib/contribs/JetsWithoutJets/trunk/example.ref (revision 347)
	+++ contrib/contribs/JetsWithoutJets/trunk/example.ref (revision 348)
	@@ -1,175 +1,182 @@
	#########
	## Read an event with 354 particles
	#########
	## Example of basic event shape analysis
	#########
	Jet parameters: R_jet=1, pTcut=200
	Trimming parameters: R_sub=0.2, fcut=0.05
	#####
	N_jet=1.99971
	N_jet_trim (using built-in)=1.97156
	N_jet_trim (using trimmer)=2
	#####
	H_T=1893.67
	H_T_trim (using built-in)=1867.45
	H_T_trim (using trimmer)=1867.45
	#####
	Missing H_T=76.1162
	Missing H_T_trim (using built-in)=91.8116
	Missing H_T_trim (using trimmer)=91.8116
	#####
	#--------------------------------------------------------------------------
	-# FastJet release 3.0.2
	+# FastJet release 3.0.3
	# M. Cacciari, G.P. Salam and G. Soyez
	# A software package for jet finding and analysis at colliders
	# http://fastjet.fr
	-#
	-# Please cite arXiv:1111.6097 if you use this package for scientific
	-# work and optionally also Phys. Lett. B641 (2006) [hep-ph/0512210].
	+#
	+# 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.
	#--------------------------------------------------------------------------
	#########
	## Example of different trimming methods
	#########
	Anti-kT jet parameters: R_jet=1, pTcut=200
	Trimming parameters: R_sub=0.2, fcut=0.05
	Untrimmed Mass = 39.9912
	Tree Trimmed Mass = 22.6675
	Jet Shape Trimmed Mass = 25.4948
	Event Shape Trimmed Mass = 25.4948
	#####
	#########
	## Example of variable pT analysis
	#########
	Jet parameters: R_jet=1, pTcut=200 (unless otherwise stated)
	Trimming parameters: R_sub=0.2, fcut=0.05
	#####
	N_jet from full pT_cut function=1.99971
	N_jet_trim from full pT_cut function=1.97156
	#####
	Variable pT_cut
	pT_cut=10, N_jet=4.57942
	pT_cut=20, N_jet=2.98989
	pT_cut=50, N_jet=2.96953
	pT_cut=100, N_jet=1.99971
	pT_cut=150, N_jet=1.99971
	#####
	pT of the hardest jet=983.637
	pT of the hardest trimmed jet=983.637
	#####
	#########
	## Example of variable R analysis
	#########
	Jet parameters: R_jet=1 (unless otherwise stated), pTcut=200
	Trimming parameters: R_sub=0.2, fcut=0.05
	#####
	N_jet from full R function=1.99971
	N_jet_trim from full R function=1.97156
	#####
	Variable Rjet
	R=0.2, N_jet=1.99997
	R=0.4, N_jet=2.00011
	R=0.6, N_jet=1.99988
	R=0.8, N_jet=2.00023
	R=1, N_jet=1.99971
	#########
	## Hardest jet axis and weights
	#########
	-Eta_hardest_jet=-0.868432
	+Without global consistency:
	+Rap_hardest_jet=-0.868431
	+Phi_hardest_jet=2.90578
	+pt_weight_hardest_jet=983.637
	+N_jet_weight_hardest_jet=0.999981
	+#
	+With global consistency:
	+Rap_hardest_jet=-0.868431
	Phi_hardest_jet=2.90578
	-pt_weight_hardest_jet=984.706
	+pt_weight_hardest_jet=983.637
	N_jet_weight_hardest_jet=0.999981
	#########
	## Examples showing more general shapes
	#########
	Jet parameters: R_jet=1, pTcut=200
	Trimming parameters: R_sub=0.2, fcut=0.05
	#####
	Summed Jet Count as event shape, R_jet=1, pT_cut=200
	Njet=1.99971
	(Storage array works? Yes)
	#
	Summed Jet Count as event shape, R_jet=1, pT_cut=200, trimming with R_sub=0.2, fcut=0.05
	Njet_trim=1.97156
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt as event shape, R_jet=1, pT_cut=200
	HT=1893.67
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt as event shape, R_jet=1, pT_cut=200, trimming with R_sub=0.2, fcut=0.05
	HT_trim=1867.45
	(Storage array works? Yes)
	#
	Missing Transverse Momentum as event shape, R_jet=1, pT_cut=200
	MissHT=76.1162
	(Storage array works? Yes)
	#
	Missing Transverse Momentum as event shape, R_jet=1, pT_cut=200, trimming with R_sub=0.2, fcut=0.05
	MissHT_trim=91.8116
	(Storage array works? Yes)
	#
	Summed Jet Mass as event shape, R_jet=1, pT_cut=200
	SigmaM=129.472
	(Storage array works? Yes)
	#
	Summed Jet Mass as event shape, R_jet=1, pT_cut=200, trimming with R_sub=0.2, fcut=0.05
	SigmaM_trim=125.566
	(Storage array works? Yes)
	#
	Summed Jet Mass^2 as event shape, R_jet=1, pT_cut=200
	SigmaM2=9852.27
	(Storage array works? Yes)
	#
	Summed Jet Mass^2 as event shape, R_jet=1, pT_cut=200, trimming with R_sub=0.2, fcut=0.05
	SigmaM2_trim=9114.19
	(Storage array works? Yes)
	#
	Trimmed Subjet multiplicity with R_jet=1, pT_cut=200, R_sub=0.2, and fcut=0.05
	SigmaNsub_trim=2.04888
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt^-2 as event shape, R_jet=1, pT_cut=200
	GenHT_-2=2.24003e-06
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt^-2 as event shape, R_jet=1, pT_cut=200, trimming with R_sub=0.2, fcut=0.05
	GenHT_-2_trim=2.2075e-06
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt^-1 as event shape, R_jet=1, pT_cut=200
	GenHT_-1=0.00211487
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt^-1 as event shape, R_jet=1, pT_cut=200, trimming with R_sub=0.2, fcut=0.05
	GenHT_-1_trim=0.00208462
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt^0 as event shape, R_jet=1, pT_cut=200
	GenHT_0=1.99971
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt^0 as event shape, R_jet=1, pT_cut=200, trimming with R_sub=0.2, fcut=0.05
	GenHT_0_trim=1.97156
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt^1 as event shape, R_jet=1, pT_cut=200
	GenHT_1=1893.67
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt^1 as event shape, R_jet=1, pT_cut=200, trimming with R_sub=0.2, fcut=0.05
	GenHT_1_trim=1867.45
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt^2 as event shape, R_jet=1, pT_cut=200
	GenHT_2=1.79595e+06
	(Storage array works? Yes)
	#
	Summed Jet Scalar Pt^2 as event shape, R_jet=1, pT_cut=200, trimming with R_sub=0.2, fcut=0.05
	GenHT_2_trim=1.77151e+06
	(Storage array works? Yes)
	#
	#####
	Index: contrib/contribs/JetsWithoutJets/trunk/JetsWithoutJets.hh
	===================================================================
	--- contrib/contribs/JetsWithoutJets/trunk/JetsWithoutJets.hh (revision 347)
	+++ contrib/contribs/JetsWithoutJets/trunk/JetsWithoutJets.hh (revision 348)
	@@ -1,672 +1,689 @@
	// JetsWithoutJets Package
	// Questions/Comments? danbert@mit.edu jthaler@mit.edu
	//
	// Copyright (c) 2013
	// Daniele Bertolini and Jesse Thaler
	//
	// $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_JETSWITHOUTJETS_HH__
	#define __FASTJET_CONTRIB_JETSWITHOUTJETS_HH__

	#include <fastjet/internal/base.hh>
	#include "fastjet/FunctionOfPseudoJet.hh"
	#include "fastjet/tools/Transformer.hh"
	#include "fastjet/JetDefinition.hh"
	#include "fastjet/PseudoJet.hh"
	#include "fastjet/ClusterSequence.hh"
	#include "fastjet/SharedPtr.hh"

	#include <string>
	#include <sstream>
	#include <algorithm>


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

	namespace contrib {


	//////
	//
	// Defining a function of a vector of PseudoJets
	//
	//////


	//----------------------------------------------------------------------
	//In current version of FastJet there is no standard interface for a function that takes a vector of PseudoJets as argument.
	//This class serves as a base class and it is the analog of FunctionOfPseudoJet, it only differs in taking a vector of PseudoJets as argument.
	//Should a standard interface become available, this class can be removed and the classes below
	//(JetMultiplicity,TransverseMomentum,MissingTransverseMomentum) should be derived from the standard base class.

	template<typename TOut>
	class MyFunctionOfVectorOfPseudoJets{
	public:

	MyFunctionOfVectorOfPseudoJets(){}

	MyFunctionOfVectorOfPseudoJets(const TOut &constant_value);

	virtual ~MyFunctionOfVectorOfPseudoJets(){}

	virtual std::string description() const{ return "";}

	virtual TOut result(const std::vector<PseudoJet> &pjs) const = 0;

	TOut operator()(const std::vector<PseudoJet> &pjs) const { return result(pjs);}
	};

	//----------------------------------------------------------------------
	// The following are example Functors to Use
	//----------------------------------------------------------------------


	//----------------------------------------------------------------------
	// Counts the number of jets (i.e. 1 per jet)
	class FunctorJetCount : public MyFunctionOfVectorOfPseudoJets<double> {
	double result(const std::vector<PseudoJet> & particles) const { return 1.0; }
	std::string description() const { return "Jet Count";}
	};


	//----------------------------------------------------------------------
	// Finds the summed pt of the jet constituents (NOT the pT of the jet)
	class FunctorJetScalarPt : public MyFunctionOfVectorOfPseudoJets<double> {
	double result(const std::vector<PseudoJet> & particles) const {
	double myPt = 0.0;
	for (unsigned int i = 0; i<particles.size(); i++) {
	myPt += particles[i].pt();
	}
	return myPt;
	}
	std::string description() const { return "Jet Scalar Pt";}
	};

	//----------------------------------------------------------------------
	// Raises FunctorJetScalarTransverseMomentum to an arbitrary Power
	class FunctorJetScalarPtToN : public MyFunctionOfVectorOfPseudoJets<double> {

	private:
	int _n;

	public:

	FunctorJetScalarPtToN(int n) : _n(n) {}
	double result(const std::vector<PseudoJet> & particles) const {
	FunctorJetScalarPt sumPt = FunctorJetScalarPt();
	return pow(sumPt(particles),_n);
	}
	std::string description() const {
	std::stringstream myStream;
	myStream << "Jet Scalar Pt^" << _n;
	return myStream.str();
	}
	};


	//----------------------------------------------------------------------
	// Finds the mass of the jet
	class FunctorJetMass : public MyFunctionOfVectorOfPseudoJets<double> {
	double result(const std::vector<PseudoJet> & particles) const {
	PseudoJet myJet(0,0,0,0);
	for (unsigned int i = 0; i<particles.size(); i++) {
	myJet += particles[i];
	}
	return myJet.m();
	}
	std::string description() const { return "Jet Mass";}
	};

	//----------------------------------------------------------------------
	// Finds the squared mass of the jet
	class FunctorJetMassSquared : public MyFunctionOfVectorOfPseudoJets<double> {
	double result(const std::vector<PseudoJet> & particles) const {
	PseudoJet myJet(0,0,0,0);
	for (unsigned int i = 0; i<particles.size(); i++) {
	myJet += particles[i];
	}
	return myJet.m2();
	}
	std::string description() const { return "Jet Mass^2";}
	};

	//////
	//
	// Storage Array (beta, to reduce computation time, only used in trimming at the moment)
	// Currently makes strips in rapidity, but could also do the same for azimuth if further speed up is needed.
	//
	//////

	class JWJStorageArray {

	private:

	double _Rjet;
	double _ptcut;

	std::vector<std::vector<std::vector<fastjet::PseudoJet> > > _regionStorage;
	std::vector<std::vector<bool> > _aboveCutBool;

	double _rapmax;

	int _maxRapIndex;
	double _rapSpread;
	int _maxPhiIndex;
	double _phiSpread;



	int getRapIndex(const fastjet::PseudoJet & currentJet) const;
	int getPhiIndex(const fastjet::PseudoJet & currentJet) const;

	public:


	JWJStorageArray() : _Rjet(0.0) {
	// at the moment, this value is hard coded
	_rapmax = 10.0;
	}

	void establishStorage(const std::vector<fastjet::PseudoJet> & my_jets, double Rjet, double ptcut);
	std::vector<fastjet::PseudoJet> & getStorageFor(const fastjet::PseudoJet & currentJet);
	bool aboveCutFor(const fastjet::PseudoJet & currentJet);

	};

	//////
	//
	// Extendable Jet-like Event Shapes (works on any function of PseudoJet that returns double)
	//
	//////


	//----------------------------------------------------------------------
	// Generic JetLikeEventShape which takes a Functor as an argument and then
	// sums over all "jets" in an event.
	// Optional trimming built in. Note that trimming first, then applying the event
	// shape gives a slightly different answer.
	// Only works with quantities that are doubles for each jet
	// (could templatize, but easier to deal with case by case)

	class JetLikeEventShape : public MyFunctionOfVectorOfPseudoJets<double> {

	protected:
	MyFunctionOfVectorOfPseudoJets<double> * _measurement;

	double _Rjet, _ptcut, _Rsub, _fcut;
	bool _trim;

	mutable bool _includeParticle;
	mutable double _weightParticle, _pt_in_Rjet, _pt_in_Rsub;
	mutable std::vector<PseudoJet> _nearbyParticles;

	void establishWeights(const PseudoJet myParticle, const std::vector<PseudoJet> & particles) const;

	bool _useStorageArray;
	mutable JWJStorageArray _myStorageArray;
	void establishStorage(const std::vector<PseudoJet> & particles) const;

	public:

	JetLikeEventShape(): _measurement(NULL), _useStorageArray(true), _myStorageArray() {}

	JetLikeEventShape(MyFunctionOfVectorOfPseudoJets<double> * measurement, double Rjet, double ptcut)
	: _measurement(measurement), _Rjet(Rjet), _ptcut(ptcut), _Rsub(Rjet), _fcut(1.0), _trim(false), _useStorageArray(true), _myStorageArray() {}
	JetLikeEventShape(MyFunctionOfVectorOfPseudoJets<double> * measurement, double Rjet, double ptcut, double Rsub, double fcut, bool useStorageArray = false)
	: _measurement(measurement), _Rjet(Rjet), _ptcut(ptcut), _Rsub(Rsub), _fcut(fcut), _trim(true), _useStorageArray(true), _myStorageArray() {}

	virtual ~JetLikeEventShape(){ if(_measurement) delete _measurement;}

	virtual double result(const std::vector<PseudoJet> & particles) const;

	void setUseStorageArray(bool useStorageArray) {_useStorageArray = useStorageArray;}

	std::string jetParameterString() const {
	std::stringstream stream;
	stream << "R_jet=" << _Rjet << ", pT_cut=" << _ptcut;
	if (_trim) stream << ", trimming with R_sub=" << _Rsub <<", fcut=" << _fcut;
	return stream.str();
	}

	virtual std::string description() const {
	return "Summed " + _measurement->description() + " as event shape, " + jetParameterString();
	}

	};

	//////
	//
	// Built-In Jet-like Event Shapes
	//
	//////

	//----------------------------------------------------------------------
	// JetMultiplicity defines event shape jet counting.
	class ShapeJetMultiplicity: public JetLikeEventShape {

	public:

	ShapeJetMultiplicity(double Rjet, double ptcut): JetLikeEventShape(new FunctorJetCount(),Rjet,ptcut) {}
	ShapeJetMultiplicity(double Rjet, double ptcut, double Rsub, double fcut) : JetLikeEventShape(new FunctorJetCount(),Rjet,ptcut,Rsub,fcut) {}
	virtual ~ShapeJetMultiplicity(){}
	};


	//----------------------------------------------------------------------
	// ScalarPT defines scalar pT sum.
	class ShapeScalarPt: public JetLikeEventShape {

	public:

	ShapeScalarPt(double Rjet, double ptcut): JetLikeEventShape(new FunctorJetScalarPt(),Rjet,ptcut) {}
	ShapeScalarPt(double Rjet, double ptcut, double Rsub, double fcut) : JetLikeEventShape(new FunctorJetScalarPt(),Rjet,ptcut,Rsub,fcut) {}
	virtual ~ShapeScalarPt(){}
	};

	//----------------------------------------------------------------------
	// TransverseMomentum defines scalar pT sum.
	class ShapeScalarPtToN: public JetLikeEventShape {

	public:

	ShapeScalarPtToN(double n, double Rjet, double ptcut): JetLikeEventShape(new FunctorJetScalarPtToN(n),Rjet,ptcut) {}
	ShapeScalarPtToN(double n, double Rjet, double ptcut, double Rsub, double fcut) :
	JetLikeEventShape(new FunctorJetScalarPtToN(n),Rjet,ptcut,Rsub,fcut) {}
	virtual ~ShapeScalarPtToN(){}
	};


	//----------------------------------------------------------------------
	// Sum over jet masses
	class ShapeSummedMass: public JetLikeEventShape {

	public:

	ShapeSummedMass(double Rjet, double ptcut): JetLikeEventShape(new FunctorJetMass(),Rjet,ptcut) {}
	ShapeSummedMass(double Rjet, double ptcut, double Rsub, double fcut) : JetLikeEventShape(new FunctorJetMass(),Rjet,ptcut,Rsub,fcut) {}
	virtual ~ShapeSummedMass(){}
	};

	//----------------------------------------------------------------------
	// Sum over jet mass^2
	class ShapeSummedMassSquared: public JetLikeEventShape {

	public:

	ShapeSummedMassSquared(double Rjet, double ptcut): JetLikeEventShape(new FunctorJetMassSquared(),Rjet,ptcut) {}
	ShapeSummedMassSquared(double Rjet, double ptcut, double Rsub, double fcut) : JetLikeEventShape(new FunctorJetMassSquared(),Rjet,ptcut,Rsub,fcut) {}
	virtual ~ShapeSummedMassSquared(){}
	};


	//////
	//
	// More complicated Jet-like Event Shapes that need special coding
	//
	//////

	//----------------------------------------------------------------------
	// ShapeMissingPt is pT of the vector sum of jets momenta.
	// Needs special coding

	class ShapeMissingPt : public JetLikeEventShape {

	public:

	ShapeMissingPt(double Rjet, double ptcut): JetLikeEventShape(NULL,Rjet,ptcut) {}
	ShapeMissingPt(double Rjet, double ptcut, double Rsub, double fcut) : JetLikeEventShape(NULL,Rjet,ptcut,Rsub,fcut) {}
	virtual ~ShapeMissingPt(){}

	double result(const std::vector<PseudoJet> & particles) const;

	std::string description() const;

	};

	//----------------------------------------------------------------------
	// ShapeTrimmedSubjetMultiplicity is a counter for subjets within trimmed fat jets.

	class ShapeTrimmedSubjetMultiplicity: public JetLikeEventShape {

	private:

	double _ptsubcut; // additional subjet pT cut, if desired

	public:

	ShapeTrimmedSubjetMultiplicity(double Rjet, double ptcut, double Rsub, double fcut, double ptsubcut = 0.0) : JetLikeEventShape(NULL,Rjet,ptcut,Rsub,fcut), _ptsubcut(ptsubcut) {}
	virtual ~ShapeTrimmedSubjetMultiplicity(){}

	double result(const std::vector<PseudoJet> & particles) const;

	std::string description() const;

	};


	//////
	//
	// Shape Trimming
	//
	//////

	//----------------------------------------------------------------------
	//Shape Trimming. It is implemented as a selector.
	Selector SelectorShapeTrimming(double Rjet, double ptcut, double Rsub, double fcut);


	//----------------------------------------------------------------------
	//Jet Shape Trimming. It is implemented as a selector, and takes the jet pT as the sum of the given four-vectors
	Selector SelectorJetShapeTrimming(double Rsub, double fcut);

	//----------------------------------------------------------------------
	//Also implement a transformer version to act on individual jets
	class JetShapeTrimmer : public Transformer {

	private:

	double _Rsub, _fcut;
	Selector _selector;

	public:

	JetShapeTrimmer(double Rsub, double fcut) : _Rsub(Rsub), _fcut(fcut){
	_selector = SelectorJetShapeTrimming(_Rsub,_fcut);
	}

	virtual PseudoJet result(const PseudoJet & original) const {
	return join(_selector(original.constituents()));
	}

	std::string jetParameterString() const {
	std::stringstream stream;
	stream << "R_sub=" << _Rsub <<", fcut=" << _fcut;
	return stream.str();
	}

	virtual std::string description() const {
	return "Jet shape trimmer, " + jetParameterString();
	}
	};


	//////
	//
	// Variable pT_cut
	//
	//////

	class JetLikeEventShape_VariablePtCut {

	protected:

	MyFunctionOfVectorOfPseudoJets<double> * _measurement;
	double _Rjet, _Rsub, _fcut, _offset;
	bool _trim;

	void _buildStepFunction(const std::vector<PseudoJet> particles) const;

	mutable std::vector<double> _ptR_values, _eventShape_values;
	mutable std::vector<PseudoJet> _nearbyParticles;
	mutable double _pt_in_Rjet, _pt_in_Rsub;

	public:

	JetLikeEventShape_VariablePtCut(): _measurement(NULL) {}

	JetLikeEventShape_VariablePtCut(MyFunctionOfVectorOfPseudoJets<double> * measurement, double Rjet, double offset = 0.0):
	_measurement(measurement), _Rjet(Rjet), _Rsub(Rjet), _fcut(1.0), _offset(offset), _trim(false) {}

	JetLikeEventShape_VariablePtCut(MyFunctionOfVectorOfPseudoJets<double> * measurement, double Rjet, double Rsub, double fcut, double offset = 0.0):
	_measurement(measurement), _Rjet(Rjet), _Rsub(Rsub), _fcut(fcut), _offset(offset), _trim(true) {}

	virtual ~JetLikeEventShape_VariablePtCut() {if(_measurement) delete _measurement;}

	// Initialization: input particles and build step function.
	void set_input(const std::vector<PseudoJet> & particles) const {_buildStepFunction(particles);}

	// Get the event shape for any value of ptcut.
	virtual double eventShapeFor(const double ptcut_0) const;

	// Pseudo-inverse: get ptcut for a given event shape value. If initialized it uses an offset, default offset is zero.
	virtual double ptCutFor(const double eventShape_0) const;

	// Get the full arrays.
	virtual std::vector<double> ptCutArray() const {return _ptR_values;}
	virtual std::vector<double> eventShapeArray() const {return _eventShape_values;}


	std::string ParameterString() const {
	std::stringstream stream;
	stream << "R_jet=" << _Rjet;
	if (_trim) stream << ", trimming with R_sub=" << _Rsub <<", fcut=" << _fcut;
	stream << ", offset for inverse function=" << _offset;
	return stream.str();
	}

	virtual std::string description() const {
	return _measurement->description() + "as function of pT_cut, " + ParameterString();
	}

	};

	//----------------------------------------------------------------------
	// Njet with variable pT_cut

	class ShapeJetMultiplicity_VariablePtCut: public JetLikeEventShape_VariablePtCut {

	public:

	// Default offset is 0.5
	ShapeJetMultiplicity_VariablePtCut(double Rjet, double offset = 0.5):
	JetLikeEventShape_VariablePtCut(new FunctorJetCount(), Rjet, offset) {};

	ShapeJetMultiplicity_VariablePtCut(double Rjet, double Rsub, double fcut, double offset = 0.5):
	JetLikeEventShape_VariablePtCut(new FunctorJetCount(), Rjet, Rsub, fcut, offset) {};

	virtual ~ShapeJetMultiplicity_VariablePtCut(){}

	};



	//////
	//
	// Njet with variable R_jet
	//
	//////

	class ShapeJetMultiplicity_VariableR {

	protected:

	double _ptcut, _Rsub, _fcut;
	bool _trim;

	void _buildStepFunction(const std::vector<PseudoJet> particles) const;

	mutable std::vector<double> _dR_values, _eventShape_values;
	mutable std::vector<PseudoJet> _nearbyParticles;
	mutable double _pt_in_Rjet;

	public:

	ShapeJetMultiplicity_VariableR() {}

	ShapeJetMultiplicity_VariableR(double ptcut): _ptcut(ptcut), _Rsub(0.0), _fcut(1.0), _trim(false) {}

	ShapeJetMultiplicity_VariableR(double ptcut, double Rsub, double fcut): _ptcut(ptcut), _Rsub(Rsub), _fcut(fcut), _trim(true) {}

	virtual ~ShapeJetMultiplicity_VariableR() {}

	// Initialization: input particles and build step function.
	void set_input(const std::vector<PseudoJet> & particles) const {_buildStepFunction(particles);}

	// Get the event shape for any value of R.
	virtual double eventShapeFor(const double Rjet_0) const;

	// Get the full arrays.
	virtual std::vector<double> RArray() const {return _dR_values;}
	virtual std::vector<double> eventShapeArray() const {return _eventShape_values;}

	std::string ParameterString() const {
	std::stringstream stream;
	stream << "pT_cut=" << _ptcut;
	if (_trim) stream << ", trimming with R_sub=" << _Rsub <<", fcut=" << _fcut;
	return stream.str();
	}

	virtual std::string description() const {
	return "Jet multiplicity as function of Rjet, " + ParameterString();
	}

	};

	//////
	//
	-// Jet Axes
	+// Finding Jet Axes with the Event Shape Density
	//
	//////


	//--------------------------------------------------------
	//Functor to find jet axis according to a given jet definition

	-class FunctorJetAxis : public MyFunctionOfVectorOfPseudoJets< std::vector<PseudoJet> > {
	+class FunctorJetAxis : public MyFunctionOfVectorOfPseudoJets< PseudoJet > {

	public:

	FunctorJetAxis(): _jetDef(NULL) {}
	- FunctorJetAxis(fastjet::JetDefinition *jetDef) : _jetDef(jetDef) {}
	+ FunctorJetAxis(fastjet::JetDefinition *jetDef) : _jetDef(jetDef) {
	+ _jetDef->delete_recombiner_when_unused();
	+ }
	~FunctorJetAxis(){ if(_jetDef) delete _jetDef;}

	- std::vector<PseudoJet> result(const std::vector<PseudoJet> & particles) const {
	-
	- fastjet::ClusterSequence clustSeq(particles, *_jetDef);
	- std::vector <fastjet::PseudoJet> myAxes = clustSeq.inclusive_jets(0.0);
	-
	- return(myAxes);
	+ PseudoJet result(const std::vector<PseudoJet> & particles) const {
	+ fastjet::ClusterSequence clustSeq(particles, *_jetDef);
	+ fastjet::PseudoJet myAxis = clustSeq.exclusive_jets(1)[0]; // should cluster everything
	+ return(myAxis);
	}
	- std::string description() const { return "Jet axes with " + _jetDef->description();}
	-
	+
	+ std::string description() const { return "Jet axis with " + _jetDef->description();}

	private:

	fastjet::JetDefinition *_jetDef;

	};

	//--------------------------------------------------------
	+// FunctorJetAxis
	+
	+class FunctorLightLikeVersion : public FunctionOfPseudoJet< PseudoJet > {
	+
	+public:
	+
	+ FunctorLightLikeVersion() {}
	+
	+ // Convert to light-like object
	+ PseudoJet result(const PseudoJet & jet) const {
	+ PseudoJet p;
	+ p.reset_PtYPhiM(1.0, jet.rap(), jet.phi(),0.0);
	+ p.set_user_index(jet.user_index()); // carry user index infomation
	+ return p;
	+ }
	+
	+ std::string description() const { return "Light-like version";}
	+
	+};
	+
	+//--------------------------------------------------------
	//Winner-take-all recombiner

	-class WinnerTakeAll : public fastjet::JetDefinition::Recombiner {
	+class WinnerTakeAllRecombiner : public fastjet::JetDefinition::Recombiner {

	public:

	- WinnerTakeAll() {}
	+ WinnerTakeAllRecombiner() {}

	/// return a textual description of the recombination scheme implemented here
	virtual std::string description() const {
	return "WinnerTakeAll Recombination Scheme";
	}

	/// recombine pa and pb and put result into pab
	virtual void recombine(const PseudoJet & pa, const PseudoJet & pb, PseudoJet & pab) const {
	PseudoJet harder = pa;
	PseudoJet softer = pb;
	if (harder.pt() < softer.pt()) std::swap(harder,softer);

	PseudoJet direction = lightLikeVersion(harder);
	double newpT = harder.pt() + softer.pt();
	pab = newpT * direction;
	}

	protected:

	- // Convert to light-like object
	- PseudoJet lightLikeVersion(const PseudoJet& jet) const {
	- double eta = jet.eta();
	- double phi = jet.phi();
	- PseudoJet p(cos(phi),sin(phi),sinh(eta),cosh(eta));
	- p.set_user_index(jet.user_index()); // carry user index infomation
	- return p;
	- }
	+ FunctorLightLikeVersion lightLikeVersion;

	};


	//--------------------------------------------------------
	//Hybrid probability density shape for jet axes position

	class EventShapeDensity_JetAxes {

	public:

	EventShapeDensity_JetAxes(){}

	- EventShapeDensity_JetAxes(double Rjet, double ptcut) : _Rjet(Rjet), _ptcut(ptcut), _applyGlobalConsistency(false) {}
	+ EventShapeDensity_JetAxes(double Rjet, double ptcut, bool applyGlobalConsistency = false)
	+ : _Rjet(Rjet), _ptcut(ptcut), _applyGlobalConsistency(applyGlobalConsistency), _useStorageArray(true) {}

	~EventShapeDensity_JetAxes(){}
	-

	//Set input and find local axes
	- void set_input(const std::vector<PseudoJet> & particles) const {_find_local_axes(particles);}
	+ void set_input(const std::vector<PseudoJet> & particles) const {
	+ _find_local_axes(particles);
	+ find_axes_and_weights();
	+ }

	//Turn on/off global consistency requirement
	void setGlobalConsistencyCheck(bool applyGlobalConsistency) {_applyGlobalConsistency = applyGlobalConsistency;}

	- //Find axes and weights
	+ //(re)find axes and weights
	void find_axes_and_weights() const;

	//Return a vector of PseudoJets sorted by pT. pT of each axis is its corresponding pT weight.
	std::vector<PseudoJet> axes() const { return _distinctAxes; }
	+
	//Return the corresponding vector of Njet weights
	std::vector<double> Njet_weights() const { return _tot_Njet_weights; }

	+ void setUseStorageArray(bool useStorageArray) {_useStorageArray = useStorageArray;}

	private:

	double _Rjet, _ptcut;
	bool _applyGlobalConsistency;

	- void _find_local_axes(const std::vector<PseudoJet> particles) const;
	+ void _find_local_axes(const std::vector<PseudoJet> & particles) const;
	bool _isStable(const int thisAxis) const;

	mutable unsigned int _N;
	mutable std::vector<double> _Njet_weights, _pt_weights, _tot_Njet_weights;
	mutable std::vector<int> _axes;
	+
	mutable std::vector<PseudoJet> _myParticles, _distinctAxes;

	- // Convert to light-like object with pt=1
	- PseudoJet _lightLikeVersion(const PseudoJet& jet) const {
	- double eta = jet.eta();
	- double phi = jet.phi();
	- PseudoJet p(cos(phi),sin(phi),sinh(eta),cosh(eta));
	- p.set_user_index(jet.user_index()); // carry user index infomation
	- return p;
	- }
	+ FunctorLightLikeVersion _lightLikeVersion;
	+
	+ bool _useStorageArray;
	+ mutable JWJStorageArray _myStorageArray;
	+
	};


	} // namespace contrib

	FASTJET_END_NAMESPACE

	#endif // __FASTJET_CONTRIB_JETSWITHOUTJETS_HH__
	Index: contrib/contribs/JetsWithoutJets/trunk/example.cc
	===================================================================
	--- contrib/contribs/JetsWithoutJets/trunk/example.cc (revision 347)
	+++ contrib/contribs/JetsWithoutJets/trunk/example.cc (revision 348)
	@@ -1,387 +1,390 @@
	// Example showing usage of JetsWithoutJets classes.
	//
	// Compile it with "make example" and run it with
	//
	// ./example < ../data/single-event.dat
	//
	// Copyright (c) 2013
	// Daniele Bertolini and Jesse Thaler
	//
	// $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 <iomanip>
	#include <iostream>
	#include <sstream>
	#include <ctime>


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

	#include "JetsWithoutJets.hh" // In external code, this should be fastjet/contrib/JetsWithoutJets.hh

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

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

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

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

	//----------------------------------------------------------
	// illustrate how this JetsWithoutJets contrib works

	analyze(event);
	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);
	}
	}

	// Main code
	void analyze(const vector<PseudoJet> & input_particles) {

	// Jet parameters.
	double Rjet = 1.0;
	double pTcut = 200.0;

	// Subjet parameters. Need for trimming.
	double Rsub = 0.2;
	double fcut = 0.05;
	double ptsubcut = 50.0; // additional ptsubcut for subjet multiplicity

	//////////
	// Basic event shape analysis.
	//////////

	// Selectors of the type SelectorEventTrimmer are trimmers of the event.
	Selector trimmer=SelectorShapeTrimming(Rjet,pTcut,Rsub,fcut);

	// Classes derived from JetLikeEventShape take as input a vector<PseudoJet>
	// representing the event, and return a number.
	// They have two constructors: you can specify Rjet and pTcut solely or
	// in addition to those also Rsub and fcut (if you want to implement
	// built-in trimming at the same time).

	// N.B.: trimming an event and then calculating the event shape is NOT
	// the same as calculating the trimmed event shape. See the README for
	// more information.

	// ShapeJetMultiplicity = jet counting
	ShapeJetMultiplicity Nj(Rjet,pTcut);
	ShapeJetMultiplicity Nj_trim(Rjet,pTcut,Rsub,fcut); // with trimming

	// ShapeScalarPt = HT (summed jet pT)
	ShapeScalarPt HT(Rjet,pTcut);
	ShapeScalarPt HT_trim(Rjet,pTcut,Rsub,fcut); // with trimming

	// ShapeMissingPt = missing pT
	ShapeMissingPt HTmiss(Rjet,pTcut);
	ShapeMissingPt HTmiss_trim(Rjet,pTcut,Rsub,fcut); // with trimming

	cout << setprecision(6);
	cout << "#########" << endl;
	cout << "## Example of basic event shape analysis" << endl;
	cout << "#########" << endl;
	cout << "Jet parameters: R_jet=" << Rjet << ", pTcut=" << pTcut <<endl;
	cout << "Trimming parameters: R_sub=" << Rsub << ", fcut=" << fcut <<endl;
	cout << "#####" << endl;
	cout << "N_jet=" << Nj(input_particles) << endl;
	cout << "N_jet_trim (using built-in)=" << Nj_trim(input_particles) << endl;
	cout << "N_jet_trim (using trimmer)=" << Nj(trimmer(input_particles)) << endl; // N.B.: Not exactly equivalent
	cout << "#####" << endl;
	cout << "H_T=" << HT(input_particles) << endl;
	cout << "H_T_trim (using built-in)=" << HT_trim(input_particles) << endl;
	cout << "H_T_trim (using trimmer)=" << HT(trimmer(input_particles)) << endl; // N.B.: Not exactly equivalent
	cout << "#####" << endl;
	cout << "Missing H_T=" << HTmiss(input_particles) << endl;
	cout << "Missing H_T_trim (using built-in)=" << HTmiss_trim(input_particles) << endl;
	cout << "Missing H_T_trim (using trimmer)=" << HTmiss(trimmer(input_particles)) << endl; // N.B.: Not exactly equivalent
	cout << "#####" << endl;

	//////////
	// Testing different trimming methods
	//////////

	// Different trimming mechanisms
	Selector eventShapeTrimmer=SelectorShapeTrimming(Rjet,pTcut,Rsub,fcut);
	JetShapeTrimmer jetShapeTrimmer(Rsub,fcut);
	Filter treeTrimmer(Rsub, SelectorPtFractionMin(fcut));

	// Clustering with anti-kt algorithm
	JetAlgorithm algorithm = antikt_algorithm;
	JetDefinition jetDef = JetDefinition(algorithm,Rjet);
	ClusterSequence clust_seq(input_particles,jetDef);
	PseudoJet hardestJet = sorted_by_pt(clust_seq.inclusive_jets(pTcut))[0];

	// Clustering after applying event shape trimming
	ClusterSequence clust_seq_estrim(eventShapeTrimmer(input_particles),jetDef);
	PseudoJet hardestJet_estrim = sorted_by_pt(clust_seq_estrim.inclusive_jets(pTcut))[0];

	double untrimmedMass = hardestJet.m();
	double treeTrimmedMass = treeTrimmer(hardestJet).m();
	double jetShapeTrimmedMass = jetShapeTrimmer(hardestJet).m();
	double eventShapeTrimmedMass = hardestJet_estrim.m();

	cout << setprecision(6);
	cout << "#########" << endl;
	cout << "## Example of different trimming methods" << endl;
	cout << "#########" << endl;
	cout << "Anti-kT jet parameters: R_jet=" << Rjet << ", pTcut=" << pTcut <<endl;
	cout << "Trimming parameters: R_sub=" << Rsub << ", fcut=" << fcut <<endl;
	cout << "Untrimmed Mass = " << untrimmedMass << endl;
	cout << "Tree Trimmed Mass = " << treeTrimmedMass << endl;
	cout << "Jet Shape Trimmed Mass = " << jetShapeTrimmedMass << endl;
	cout << "Event Shape Trimmed Mass = " << eventShapeTrimmedMass << endl;
	cout << "#####" << endl;

	//////////
	// Variable pT analysis.
	//////////

	// These use classes derived from JetLikeEventShape_VariablePtCut, where
	// the full pT information is kept in calculating the event shape.
	// One can calculate the value of the event shape at a given pT cut, or
	// one can calculate the inverse function to find the value of the pT cut
	// needed to return a given value of the event shape.

	ShapeJetMultiplicity_VariablePtCut Nj_allPt(Rjet);
	ShapeJetMultiplicity_VariablePtCut Nj_allPt_trim(Rjet,Rsub,fcut); // with trimming
	Nj_allPt.set_input(input_particles);
	Nj_allPt_trim.set_input(input_particles);

	cout << "#########" << endl;
	cout << "## Example of variable pT analysis" << endl;
	cout << "#########" << endl;
	cout << "Jet parameters: R_jet=" << Rjet << ", pTcut=" << pTcut << " (unless otherwise stated)" <<endl;
	cout << "Trimming parameters: R_sub=" << Rsub << ", fcut=" << fcut <<endl;
	cout << "#####" << endl;
	cout << "N_jet from full pT_cut function=" << Nj_allPt.eventShapeFor(pTcut) << endl;
	cout << "N_jet_trim from full pT_cut function=" << Nj_allPt_trim.eventShapeFor(pTcut) << endl;
	cout << "#####" << endl;
	cout << "Variable pT_cut" << endl;
	double pTcut_values[] = {10,20,50,100,150};
	for (unsigned int i=0; i<5; i++){
	cout << "pT_cut=" << pTcut_values[i] << ", N_jet=" << Nj_allPt.eventShapeFor(pTcut_values[i]) << endl;
	}
	cout << "#####" << endl;
	cout << "pT of the hardest jet=" << Nj_allPt.ptCutFor(1) << endl;
	cout << "pT of the hardest trimmed jet=" << Nj_allPt_trim.ptCutFor(1) << endl;
	cout << "#####" << endl;

	// Can get the full arrays too, though these are not output on the example.ref file
	std::vector<double> ptval = Nj_allPt.ptCutArray();
	std::vector<double> NjPtval = Nj_allPt.eventShapeArray();

	//////////
	// Variable R analysis.
	//////////

	// These use the class ShapeJetMultiplicity_VariableR (at present there is no
	// general class for variable R observables because of the high computational time).
	// Here, the full R information is kept.

	ShapeJetMultiplicity_VariableR Nj_allR(pTcut);
	ShapeJetMultiplicity_VariableR Nj_allR_trim(pTcut,Rsub,fcut); // with trimming
	Nj_allR.set_input(input_particles);
	Nj_allR_trim.set_input(input_particles);

	cout << "#########" << endl;
	cout << "## Example of variable R analysis" << endl;
	cout << "#########" << endl;
	cout << "Jet parameters: R_jet=" << Rjet << " (unless otherwise stated), pTcut=" << pTcut <<endl;
	cout << "Trimming parameters: R_sub=" << Rsub << ", fcut=" << fcut <<endl;
	cout << "#####" << endl;
	cout << "N_jet from full R function=" << Nj_allR.eventShapeFor(Rjet) << endl;
	cout << "N_jet_trim from full R function=" << Nj_allR_trim.eventShapeFor(Rjet) << endl;
	cout << "#####" << endl;
	cout << "Variable Rjet" << endl;
	double Rjet_values[] = {0.2,0.4,0.6,0.8,1};
	for(unsigned int i=0; i<5; i++){
	cout << "R=" << Rjet_values[i] << ", N_jet=" << Nj_allR.eventShapeFor(Rjet_values[i]) << endl;
	}

	// Can get the full arrays too, though these are not output on the example.ref file
	std::vector<double> Rval = Nj_allR.RArray();
	std::vector<double> NjRval = Nj_allR.eventShapeArray();


	//////////
	- // Jet Axes
	+ // Finding Jet Axes
	//////////

	- EventShapeDensity_JetAxes axes_eventShape(Rjet,0);
	- axes_eventShape.set_input(input_particles);
	-
	- //Find axes and weights
	- axes_eventShape.find_axes_and_weights();
	+ EventShapeDensity_JetAxes axes_eventShape(Rjet,pTcut);
	+ axes_eventShape.set_input(input_particles); //Loads particles and finds axes and weights
	+
	//Get jet axes (by default sorted by pt). Pt of each axis corresponds to its pt weight.
	vector<PseudoJet> axes = axes_eventShape.axes();
	//Get N_jet weights
	vector<double> Njw = axes_eventShape.Njet_weights();

	- //If you want turn on global consistency and recalculate axes and weights. Default is no.
	+ //If you want turn on global consistency and recalculate axes and weights
	axes_eventShape.setGlobalConsistencyCheck(true);
	- axes_eventShape.find_axes_and_weights();
	-
	- axes = axes_eventShape.axes();
	- Njw = axes_eventShape.Njet_weights();
	+ axes_eventShape.find_axes_and_weights(); // uses stored information to save time; don't need to set_input again
	+ vector<PseudoJet> axes_gc = axes_eventShape.axes();
	+ vector<double> Njw_gc = axes_eventShape.Njet_weights();

	//Output hardest jet eta/phi/pT/N_jet weight
	cout << "#########" << endl;
	cout << "## Hardest jet axis and weights" << endl;
	cout << "#########" << endl;
	- cout << "Eta_hardest_jet=" << axes[0].eta() <<endl;
	+ cout << "Without global consistency:" << endl;
	+ cout << "Rap_hardest_jet=" << axes[0].rap() <<endl;
	cout << "Phi_hardest_jet=" << axes[0].phi() <<endl;
	cout << "pt_weight_hardest_jet=" << axes[0].pt() <<endl;
	cout << "N_jet_weight_hardest_jet=" << Njw[0] <<endl;
	-
	+ cout << "#" << endl;
	+ cout << "With global consistency:" << endl;
	+ cout << "Rap_hardest_jet=" << axes_gc[0].rap() <<endl;
	+ cout << "Phi_hardest_jet=" << axes_gc[0].phi() <<endl;
	+ cout << "pt_weight_hardest_jet=" << axes_gc[0].pt() <<endl;
	+ cout << "N_jet_weight_hardest_jet=" << Njw_gc[0] <<endl;

	//////////
	// Showing a wide range of event shapes
	//////////

	// This emphasizes that the jet-like event shapes are all derived from
	// MyFunctionOfVectorOfPseudoJets
	vector<JetLikeEventShape*> myShapes;
	vector<string > myNames;

	myShapes.push_back( new ShapeJetMultiplicity(Rjet,pTcut) );
	myNames.push_back( "Njet");
	myShapes.push_back( new ShapeJetMultiplicity(Rjet,pTcut,Rsub,fcut) ); //with trimming
	myNames.push_back( "Njet_trim");

	myShapes.push_back( new ShapeScalarPt(Rjet,pTcut) );
	myNames.push_back( "HT");
	myShapes.push_back( new ShapeScalarPt(Rjet,pTcut,Rsub,fcut) ); //with trimming
	myNames.push_back( "HT_trim");

	myShapes.push_back( new ShapeMissingPt(Rjet,pTcut) );
	myNames.push_back( "MissHT");
	myShapes.push_back( new ShapeMissingPt(Rjet,pTcut,Rsub,fcut) ); //with trimming
	myNames.push_back( "MissHT_trim");

	myShapes.push_back( new ShapeSummedMass(Rjet,pTcut) );
	myNames.push_back( "SigmaM");
	myShapes.push_back( new ShapeSummedMass(Rjet,pTcut,Rsub,fcut) ); //with trimming
	myNames.push_back( "SigmaM_trim");

	myShapes.push_back( new ShapeSummedMassSquared(Rjet,pTcut) );
	myNames.push_back( "SigmaM2");
	myShapes.push_back( new ShapeSummedMassSquared(Rjet,pTcut,Rsub,fcut) ); //with trimming
	myNames.push_back( "SigmaM2_trim");

	myShapes.push_back( new ShapeTrimmedSubjetMultiplicity(Rjet,pTcut,Rsub,fcut,ptsubcut) ); //with trimming
	myNames.push_back( "SigmaNsub_trim");

	for (int n = -2; n <= 2; n++) {
	myShapes.push_back( new ShapeScalarPtToN(n,Rjet,pTcut) );
	myShapes.push_back( new ShapeScalarPtToN(n,Rjet,pTcut,Rsub,fcut) ); //with trimming

	stringstream myStream;
	myStream << "GenHT_" << n;
	myNames.push_back( myStream.str());
	myNames.push_back( myStream.str() + "_trim");
	}

	// outputing the variables
	cout << "#########" << endl;
	cout << "## Examples showing more general shapes" << endl;
	cout << "#########" << endl;
	cout << "Jet parameters: R_jet=" << Rjet << ", pTcut=" << pTcut << endl;
	cout << "Trimming parameters: R_sub=" << Rsub << ", fcut=" << fcut << endl;
	cout << "#####" << endl;
	for (unsigned int i = 0; i < myShapes.size(); i++) {
	cout << myShapes.at(i)->description() << endl; // ideally each even shape should have a description

	myShapes.at(i)->setUseStorageArray(false);
	double valueNoStorage = myShapes.at(i)->result(input_particles);
	myShapes.at(i)->setUseStorageArray(true);
	double valueYesStorage = myShapes.at(i)->result(input_particles);

	cout << myNames.at(i) << "="
	<< valueNoStorage << endl;
	cout << "(Storage array works? " << (valueNoStorage == valueYesStorage ? "Yes" : "No!!") << ")" << endl;
	cout << "#" << endl;

	delete myShapes.at(i);
	}
	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;

	cout << setprecision(6);

	num_iter = 1000;
	clock_begin = clock();
	ShapeScalarPt NjetA(Rjet,pTcut);
	cout << "timing " << NjetA.description() << endl;
	for (int t = 0; t < num_iter; t++) {
	NjetA(input_particles);
	}
	clock_end = clock();
	cout << "Method A: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per Njet"<< endl;

	num_iter = 1000;
	clock_begin = clock();
	JetLikeEventShape NjetB(new FunctorJetScalarPt(),Rjet,pTcut);
	NjetB.setUseStorageArray(false);
	cout << "timing " << NjetB.description() << endl;
	for (int t = 0; t < num_iter; t++) {
	NjetB(input_particles);
	}
	clock_end = clock();
	cout << "Method B: " << (clock_end-clock_begin)/double(CLOCKS_PER_SECnum_iter)1000 << " ms per Njet"<< endl;


	}

	}
	Index: contrib/contribs/JetsWithoutJets/trunk/JetsWithoutJets.cc
	===================================================================
	--- contrib/contribs/JetsWithoutJets/trunk/JetsWithoutJets.cc (revision 347)
	+++ contrib/contribs/JetsWithoutJets/trunk/JetsWithoutJets.cc (revision 348)
	@@ -1,788 +1,803 @@
	// JetsWithoutJets Package
	// Questions/Comments? danbert@mit.edu jthaler@mit.edu
	//
	// Copyright (c) 2013
	// Daniele Bertolini and Jesse Thaler
	//
	// $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 "JetsWithoutJets.hh"
	#include <sstream>

	using namespace std;


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

	namespace contrib {

	//////
	//
	// Storage Array (beta, to reduce computation time, only used in trimming at the moment)
	// Currently makes strips in rapidity, but could also do the same for azimuth if further speed up is needed.
	//
	//////

	void JWJStorageArray::establishStorage(const vector<fastjet::PseudoJet> & my_jets, double Rjet, double ptcut) {

	// set radius and ptcut
	_Rjet = Rjet;
	_ptcut = ptcut;
	FunctorJetScalarPt sumPt = FunctorJetScalarPt();

	// Dividing up Phase Space (need to be synced with getRapIndex/getPhiIndex)
	// Make sure that there are equal spacing of regions
	_maxRapIndex = (int) floor((_rapmax)/_Rjet);// 2 rapmax / 2 Rjet
	_rapSpread = 2.0*_rapmax/(floor((_rapmax)/_Rjet));
	_maxPhiIndex = (int) floor((M_PI)/_Rjet); // 2 pi / 2 Rjet
	_phiSpread = 2.0*M_PI/(floor((M_PI)/_Rjet));

	// Initializing storage
	_regionStorage.resize(_maxRapIndex);
	_aboveCutBool.resize(_maxRapIndex);
	for (int rapIndex = 0 ; rapIndex < _maxRapIndex ; rapIndex++) {
	_regionStorage[rapIndex].resize(_maxPhiIndex);
	_aboveCutBool[rapIndex].resize(_maxPhiIndex);
	for (int phiIndex = 0; phiIndex < _maxPhiIndex; phiIndex++) {
	_regionStorage[rapIndex][phiIndex].clear();
	}
	}

	// looping through particles and assigning to bins of size 2R by 2R (approximately)
	for (unsigned int i = 0; i < my_jets.size(); i++) {
	PseudoJet currentJet = my_jets[i];
	double rap = currentJet.rap();
	double phi = currentJet.phi();

	int lowRap = (int) floor((rap + _rapmax)/_rapSpread);
	int highRap = (int) ceil((rap + _rapmax)/_rapSpread);
	int lowPhi = (int) floor(phi/_phiSpread);
	int highPhi = (int) ceil(phi/_phiSpread);
	if (highPhi >= _maxPhiIndex) highPhi = highPhi - _maxPhiIndex; // loop around

	if (lowRap < 0) lowRap = 0;
	if (lowRap >= _maxRapIndex) lowRap = _maxRapIndex-1;
	if (highRap < 0) highRap = 0;
	if (highRap >= _maxRapIndex) highRap = _maxRapIndex-1;

	// Fill in storage. For particles at periphery in rapidity, only fill two bins instead of four. If phi overlaps, do the same thing.
	_regionStorage[lowRap][lowPhi].push_back(currentJet);
	if (lowPhi != highPhi) _regionStorage[lowRap][highPhi].push_back(currentJet);
	if (lowRap != highRap) _regionStorage[highRap][lowPhi].push_back(currentJet);
	if (lowRap != highRap && lowPhi != highPhi) _regionStorage[highRap][highPhi].push_back(currentJet);

	}

	// store whether above ptCut values
	for (int rapIndex = 0; rapIndex < _maxRapIndex; rapIndex++) {
	for (int phiIndex = 0; phiIndex < _maxPhiIndex; phiIndex++) {
	_aboveCutBool[rapIndex][phiIndex] = (sumPt(_regionStorage[rapIndex][phiIndex]) >= ptcut);
	}
	}

	}

	int JWJStorageArray::getRapIndex(const fastjet::PseudoJet & currentJet) const {
	double rap = currentJet.rap();
	int rapIndex = round((rap + _rapmax)/_rapSpread);

	if (rapIndex < 0) rapIndex = 0;
	if (rapIndex >= _maxRapIndex) rapIndex = _maxRapIndex - 1;

	return rapIndex;
	}

	int JWJStorageArray::getPhiIndex(const fastjet::PseudoJet & currentJet) const {
	double phi = currentJet.phi();
	int phiIndex = round(phi/_phiSpread);

	if (phiIndex >= _maxPhiIndex) phiIndex = phiIndex - _maxPhiIndex;

	return phiIndex;
	}

	vector<fastjet::PseudoJet> & JWJStorageArray::getStorageFor(const fastjet::PseudoJet & currentJet) {
	return _regionStorage[getRapIndex(currentJet)][getPhiIndex(currentJet)];
	}

	bool JWJStorageArray::aboveCutFor(const fastjet::PseudoJet & currentJet) {
	return _aboveCutBool[getRapIndex(currentJet)][getPhiIndex(currentJet)];
	}


	//////
	//
	// Extendable Jet-like Event Shape
	//
	//////

	// If called, this helper function establishes the storage array for speed up
	void JetLikeEventShape::establishStorage(const std::vector<PseudoJet> & particles) const {
	if (_useStorageArray) {
	_myStorageArray.establishStorage(particles,_Rjet, _ptcut);
	} else {
	// do nothing
	}
	}



	// This helper function calculates the weights assigned to each particle.
	// In most cases, the user does not need to access the information set by this function
	void JetLikeEventShape::establishWeights(const PseudoJet myParticle, const std::vector<PseudoJet> & particles) const {

	FunctorJetScalarPt sumPt = FunctorJetScalarPt();
	_includeParticle = false;
	_weightParticle = 0.0;
	_pt_in_Rjet = 0.0;
	_pt_in_Rsub = 0.0;

	// Find particles in my _Rjet neighborhood
	fastjet::Selector sel = fastjet::SelectorCircle(_Rjet);
	sel.set_reference(myParticle);

	- if (_useStorageArray) { // start from rapidity strips
	+ if (_useStorageArray) { // start from rapidity/phi blocks
	if (!_myStorageArray.aboveCutFor(myParticle)) {
	return; //don't do any further analysis if not above jet pTcut
	} else {
	_nearbyParticles = sel(_myStorageArray.getStorageFor(myParticle));
	}
	} else { // use all particles
	_nearbyParticles = sel(particles);
	}
	// See if there is enough pt in the neighborhood to pass _ptcut
	_pt_in_Rjet = sumPt(_nearbyParticles);
	if (_pt_in_Rjet < _ptcut) return;

	// If trimming is on, do check of _fcut as well
	if (_trim) {
	assert(_Rsub <= _Rjet);
	fastjet::Selector sel_sub = fastjet::SelectorCircle(_Rsub);
	sel_sub.set_reference(myParticle);
	std::vector<fastjet::PseudoJet> near_particles_sub = sel_sub(_nearbyParticles); // Save time and only look at near particles
	_pt_in_Rsub = sumPt(near_particles_sub);
	if (_pt_in_Rsub/_pt_in_Rjet < _fcut) return;
	}

	// If enough pt, find weight and apply measurement
	_includeParticle = true;
	_weightParticle = myParticle.pt()/_pt_in_Rjet;

	}

	// For a generic event shape, the result is obtained by simply multiplying the weight
	// of a given particle by the _measurement acting on the _nearbyParticles
	double JetLikeEventShape::result(const std::vector<PseudoJet> & particles) const {

	double answer = 0.0;
	establishStorage(particles);

	for (unsigned int i = 0; i < particles.size(); i++) {
	establishWeights(particles[i],particles);
	if (_includeParticle) {
	double measurement = _measurement->result(_nearbyParticles);
	answer += measurement*_weightParticle;
	}
	}

	return(answer);
	}



	//////
	//
	// MissingTransverseMomentum
	//
	//////

	// To calculate missing pT, we cannot use a standard _measurement, so we code this by hand.
	double ShapeMissingPt::result(const std::vector<PseudoJet> & particles) const {

	PseudoJet p_tot(0,0,0,0);
	establishStorage(particles);


	for (unsigned int i = 0; i < particles.size(); i++) {
	establishWeights(particles[i],particles);
	if (_includeParticle) {
	p_tot += particles[i];
	}
	}

	return(p_tot.pt());
	}


	string ShapeMissingPt::description() const {
	return "Missing Transverse Momentum as event shape, " + jetParameterString();
	}


	//////
	//
	// TrimmedSubjetMultiplicity
	//
	//////

	// While it would be possible to calculate trimmed subjet multiplicity using a
	// _measurement functor, it is computationally more efficient to use the information
	// already captured by establishWeights().
	double ShapeTrimmedSubjetMultiplicity::result(const vector<PseudoJet> & particles) const {

	double answer = 0.0;
	establishStorage(particles);

	for (unsigned int i = 0; i < particles.size(); i++) {
	establishWeights(particles[i],particles);
	if (_includeParticle && _pt_in_Rsub > _ptsubcut) {
	answer += particles[i].pt()/_pt_in_Rsub;
	}
	}
	return(answer);

	}

	string ShapeTrimmedSubjetMultiplicity::description() const {
	ostringstream oss;
	oss << "Trimmed Subjet multiplicity with R_jet=" << _Rjet << ", pT_cut=" << _ptcut << ", R_sub=" << _Rsub << ", and fcut=" << _fcut;
	return oss.str();
	}

	//////
	//
	// Shape Trimming
	//
	//////


	//----------------------------------------------------------------------
	// Helper for SelectorShapeTrimming.
	// Class derived from SelectorWorker: pass, terminator, applies_jet_by_jet, and description are overloaded.
	// It can't be applied to an individual jet since it requires information about the neighbourhood of the jet.
	class SW_ShapeTrimming: public SelectorWorker {


	private:

	double _Rjet, _ptcut, _Rsub, _fcut;

	bool _useTempStorage;

	public:

	SW_ShapeTrimming(double Rjet, double ptcut, double Rsub, double fcut, bool useTempStorage = true) : _Rjet(Rjet), _ptcut(ptcut), _Rsub(Rsub), _fcut(fcut), _useTempStorage(useTempStorage) {};

	virtual bool pass(const PseudoJet &) const {
	if (!applies_jet_by_jet())
	throw Error("Cannot apply this selector worker to an individual jet");
	return false;
	}

	virtual void terminator(vector<const PseudoJet *> & jets) const {

	// Copy pointers content to a vector of PseudoJets. Check if the pointer has been already nullified.
	vector<PseudoJet> my_jets;
	vector<unsigned int> indices;
	for (unsigned int i=0; i < jets.size(); i++){
	if (jets[i]){
	indices.push_back(i);
	my_jets.push_back(*jets[i]);
	}
	}

	FunctorJetScalarPt sumPt = FunctorJetScalarPt();
	JWJStorageArray myStorageArray = JWJStorageArray();

	-
	- // to increase speed in very high multiplicity events, make a grid of rapidity strips
	+ // to increase speed in very high multiplicity events, make a grid of rapidity/phi blocks
	if (_useTempStorage) {
	myStorageArray.establishStorage(my_jets,_Rjet,_ptcut);
	}

	for (unsigned int i=0; i < my_jets.size(); i++){

	// Select particles in radius _Rjet
	fastjet::Selector sel = fastjet::SelectorCircle(_Rjet);
	sel.set_reference(my_jets[i]);
	vector<PseudoJet> near_particles;

	- if (_useTempStorage) { // start from rapidity strips
	+ if (_useTempStorage) { // start from rapidity/phi blocks
	if (!myStorageArray.aboveCutFor(my_jets[i])) {
	jets[indices[i]] = NULL;
	continue; // no need to consider that particle further
	} else {
	near_particles = sel(myStorageArray.getStorageFor(my_jets[i]));
	}
	} else { // use all particles
	near_particles = sel(my_jets);
	}

	// Calculate enclosed pT
	double pt_Rjet = sumPt(near_particles);

	// Check if pT is enough
	// Need to have < instead of <= in order to have consistent description with shape variables
	if (pt_Rjet < _ptcut) {
	jets[indices[i]] = NULL;
	continue;
	}

	// If so, select particles in radius
	fastjet::Selector sel_sub = fastjet::SelectorCircle(_Rsub);
	sel_sub.set_reference(my_jets[i]);
	vector<PseudoJet> near_particles_sub = sel_sub(near_particles);


	// Calculate enclosed pT
	double pt_Rsub = sumPt(near_particles_sub);

	// Need to have < instead of <= in order to have consistent description with shape variables
	if (pt_Rsub/pt_Rjet < _fcut) {
	jets[indices[i]] = NULL;
	continue;
	}
	}
	}

	virtual bool applies_jet_by_jet() const {return false;}

	std::string jetParameterString() const {
	std::stringstream stream;
	stream << "R_jet=" << _Rjet << ", pT_cut=" << _ptcut << ", R_sub=" << _Rsub <<", fcut=" << _fcut;
	return stream.str();
	}

	virtual string description() const {
	return "Shape trimmer, " + jetParameterString();
	}
	};

	//----------------------------------------------------------------------
	// Helper for SelectorJetShapeTrimming.
	// Class derived from SelectorWorker: pass, terminator, applies_jet_by_jet, and description are overloaded.
	// This cannot be applied to an individual jet.
	class SW_JetShapeTrimming: public SelectorWorker {

	private:

	double _Rsub, _fcut;

	public:

	SW_JetShapeTrimming(double Rsub, double fcut) : _Rsub(Rsub), _fcut(fcut) {};

	virtual bool pass(const PseudoJet &) const {
	if (!applies_jet_by_jet())
	throw Error("Cannot apply this selector worker to an individual jet");
	return false;
	}

	virtual void terminator(vector<const PseudoJet *> & jets) const {

	// Copy pointers content to a vector of PseudoJets. Check if the pointer has been already nullified.
	vector<PseudoJet> my_jets;
	vector<unsigned int> indices;
	for (unsigned int i=0; i < jets.size(); i++){
	if (jets[i]){
	indices.push_back(i);
	my_jets.push_back(*jets[i]);
	}
	}

	FunctorJetScalarPt sumPt = FunctorJetScalarPt();

	// take sum pt of given jet
	double pt_Rjet = sumPt(my_jets);

	for (unsigned int i=0; i < my_jets.size(); i++){

	// Select particles in sub radius
	fastjet::Selector sel_sub = fastjet::SelectorCircle(_Rsub);
	sel_sub.set_reference(my_jets[i]);
	vector<PseudoJet> near_particles_sub = sel_sub(my_jets);

	// Calculate enclosed pT
	double pt_Rsub = sumPt(near_particles_sub);

	// Need to have < instead of <= in order to have consistent description with shape variables
	if (pt_Rsub/pt_Rjet < _fcut) {
	jets[indices[i]] = NULL;
	continue;
	}
	}
	}

	virtual bool applies_jet_by_jet() const {return false;}

	std::string jetParameterString() const {
	std::stringstream stream;
	stream << "R_sub=" << _Rsub <<", fcut=" << _fcut;
	return stream.str();
	}

	virtual string description() const {
	return "Jet shape trimmer, " + jetParameterString();
	}
	};


	Selector SelectorShapeTrimming(double Rjet, double ptcut, double Rsub, double fcut){
	return Selector(new SW_ShapeTrimming(Rjet,ptcut,Rsub,fcut));
	}

	Selector SelectorJetShapeTrimming(double Rsub, double fcut){
	return Selector(new SW_JetShapeTrimming(Rsub,fcut));
	}

	//////
	//
	// Variable pT_cut
	//
	//////

	-
	// helper to sort a vector of vector<double> by comparing just the first entry.
	bool _mySortFunction (std::vector<double> v_0, std::vector<double> v_1) { return (v_0[0] > v_1[0]); }

	// As described in the appendix of the physics paper, one can construct the inverse
	// of an event shape with respect to pT by calculating all of the possible pTi,R values and sorting them.
	// This helper function accomplishes that task.
	void JetLikeEventShape_VariablePtCut::_buildStepFunction(const std::vector<PseudoJet> particles) const {

	FunctorJetScalarPt sumPt = FunctorJetScalarPt();
	fastjet::Selector sel = fastjet::SelectorCircle(_Rjet);
	_pt_in_Rjet = 0.0;
	_pt_in_Rsub = 0.0;

	std::vector< std::vector<double> > myValues;
	myValues.resize(0);

	// this acts much like establishWeights did for JetLikeEventShape,
	// except each value of _pt_in_Rjet is associated with the corresponding
	// weighted measurement.
	for (unsigned int i = 0; i < particles.size(); i++){
	sel.set_reference(particles[i]);
	_nearbyParticles = sel(particles);
	_pt_in_Rjet = sumPt(_nearbyParticles);

	if(_trim) {
	assert(_Rsub <= _Rjet);
	fastjet::Selector sel_sub = fastjet::SelectorCircle(_Rsub);
	sel_sub.set_reference(particles[i]);
	std::vector<PseudoJet> near_particles_sub = sel_sub(_nearbyParticles);
	_pt_in_Rsub = sumPt(near_particles_sub);
	if (_pt_in_Rsub / _pt_in_Rjet >= _fcut) {
	double measurement = _measurement->result(_nearbyParticles);
	std::vector<double> point(2);
	point[0] = _pt_in_Rjet;
	point[1] = measurement*particles[i].pt()/_pt_in_Rjet;
	myValues.push_back(point);
	}
	}

	else {

	double measurement = _measurement->result(_nearbyParticles);
	std::vector<double> point(2);
	point[0] = _pt_in_Rjet;
	point[1] = measurement*particles[i].pt()/_pt_in_Rjet;
	myValues.push_back(point);
	}
	}



	std::sort (myValues.begin(), myValues.end(), _mySortFunction);

	// Calculating the cumulative sum of the measurements up to a given pt value
	_ptR_values.resize(0);
	_eventShape_values.resize(0);
	if (!myValues.empty()){

	_ptR_values.push_back(myValues[0][0]);
	_eventShape_values.push_back(myValues[0][1]);

	for (unsigned int i = 1; i < myValues.size(); i++){
	_ptR_values.push_back(myValues[i][0]);
	_eventShape_values.push_back(_eventShape_values[i-1] + myValues[i][1]);
	}
	}

	myValues.clear();
	}

	double JetLikeEventShape_VariablePtCut::eventShapeFor(const double ptcut_0) const {

	// TODO: find a more efficient data structure to do this searching

	double eventShape = 0.0;

	if ( ptcut_0 <= _ptR_values.back() ) eventShape = _eventShape_values.back();

	else if (ptcut_0 > _ptR_values.back() && ptcut_0 <= _ptR_values.front()) {
	for (unsigned int i = 0; i < _ptR_values.size()-1; i++){
	if (ptcut_0 <= _ptR_values[i] && ptcut_0 > _ptR_values[i+1]){
	eventShape = _eventShape_values[i];
	break;
	}
	}
	}
	return (eventShape);
	}

	double JetLikeEventShape_VariablePtCut::ptCutFor(const double eventShape_0) const {

	// TODO: find a more efficient data structure to do this searching

	double ptcut = 0.0;

	double new_eventShape_0 = eventShape_0 - _offset;

	if ( new_eventShape_0 < 0 \|\| new_eventShape_0 > _eventShape_values.back() ) {
	cout << ">>> ERROR: "<< _measurement->description() <<" min allowed value is 0, max allowed value is " << _eventShape_values.back() << endl;

	} else if ( new_eventShape_0 <= _eventShape_values.front() ) ptcut = _ptR_values.front();

	else {
	for (unsigned int i = 0; i < _eventShape_values.size()-1; i++){
	if (new_eventShape_0 > _eventShape_values[i] && new_eventShape_0 <= _eventShape_values[i+1]){
	ptcut = _ptR_values[i+1];
	break;
	}
	}
	}

	return(ptcut);
	}


	//////
	//
	// Njet with variable R_jet
	//
	//////

	// Similar to the same named function in JetLikeEventShape_VariablePtCut, except now we have
	// to store a lot more information.
	void ShapeJetMultiplicity_VariableR::_buildStepFunction(const std::vector<PseudoJet> particles) const {

	unsigned int N = particles.size();

	FunctorJetScalarPt sumPt = FunctorJetScalarPt();

	std::vector< std::vector<double> > myValues;
	myValues.resize(0);

	for(unsigned int i = 0; i < N; i++) {
	unsigned int j = i+1;
	while(j < N) {
	vector<double> myPair(3);
	myPair[0] = particles[i].delta_R(particles[j]);
	myPair[1] = i;
	myPair[2] = j;
	myValues.push_back(myPair);
	j++;
	}
	}

	std::sort (myValues.begin(), myValues.end(), _mySortFunction);

	_eventShape_values.resize(0);
	_dR_values.resize(0);

	// Initialize
	double _tot_pt = sumPt(particles);
	if(_tot_pt >= _ptcut) _eventShape_values.push_back(1.0);
	else _eventShape_values.push_back(0.0);

	std::vector<double> pTR(N);
	std::fill(pTR.begin(), pTR.end(), _tot_pt);

	if (_trim) {
	std::vector<double> pTRsub;
	fastjet::Selector sel_sub = fastjet::SelectorCircle(_Rsub);
	for (unsigned int i = 0; i < N; i++) {
	sel_sub.set_reference(particles[i]);
	std::vector<PseudoJet> near_particles_sub = sel_sub(particles);
	pTRsub.push_back(sumPt(near_particles_sub));
	}


	for (unsigned int i = 0; i < myValues.size(); i++) {

	_dR_values.push_back(myValues[i][0]);

	int id1 = (int)myValues[i][1];
	int id2 = (int)myValues[i][2];

	pTR[id1] -= particles[id2].pt();
	pTR[id2] -= particles[id1].pt();

	double myEventShape = 0;

	for(unsigned int j = 0; j < N; j++) {
	if(pTR[j] >= _ptcut && pTRsub[j] / pTR[j] >= _fcut) myEventShape += particles[j].pt() / pTR[j];
	}
	_eventShape_values.push_back(myEventShape);
	}
	}

	else {

	for (unsigned int i = 0; i < myValues.size(); i++) {

	_dR_values.push_back(myValues[i][0]);

	int id1 = (int)myValues[i][1];
	int id2 = (int)myValues[i][2];

	pTR[id1] -= particles[id2].pt();
	pTR[id2] -= particles[id1].pt();

	double myEventShape = 0;

	for(unsigned int j = 0; j < N; j++) {
	if(pTR[j] >= _ptcut) myEventShape += particles[j].pt() / pTR[j];
	}
	_eventShape_values.push_back(myEventShape);
	}
	}

	_dR_values.push_back(0.0);
	myValues.clear();

	}


	double ShapeJetMultiplicity_VariableR::eventShapeFor(const double Rjet_0) const {

	double eventShape = 0.0;

	if( Rjet_0 < _Rsub ) cout <<">>> WARNING: R_jet < R_sub= " << _Rsub << endl;

	if( Rjet_0 >= _dR_values.front() ) eventShape = _eventShape_values.front();
	else if ( Rjet_0 < _dR_values.front() && Rjet_0 >= _dR_values.back() ){

	for(unsigned int i = 0; i < _dR_values.size()-1; i++) {
	if( Rjet_0 < _dR_values[i] && Rjet_0 >= _dR_values[i+1]) {
	eventShape = _eventShape_values[i+1];
	break;
	}
	}
	}
	else cout <<">>> ERROR: R_jet should be >= 0" << endl;
	return(eventShape);
	}


	//////
	//
	-// Jet Axes
	+// Finding Jet Axes with the Event Shape Density
	//
	//////

	-void EventShapeDensity_JetAxes::_find_local_axes(const std::vector<PseudoJet> particles) const {
	-
	+void EventShapeDensity_JetAxes::_find_local_axes(const std::vector<PseudoJet> & particles) const {
	+
	_myParticles = particles;
	_N = particles.size();

	_axes.resize(0);
	_Njet_weights.resize(0);
	_pt_weights.resize(0);
	+
	+ // Indexing particles (note that this is happening internally, so shouldn't modify user input)
	+ for (unsigned int i=0; i<_N; i++) _myParticles[i].set_user_index(i);

	- //Indexing particles.
	- for(unsigned int i=0; i<_N; i++) _myParticles[i].set_user_index(i);
	-
	- //Find axis associated with each particle.
	- for(unsigned int i=0; i<_N; i++) {
	+ // For speed up/
	+ if (_useStorageArray) {
	+ _myStorageArray.establishStorage(_myParticles,_Rjet, _ptcut);
	+ }
	+
	+ // Find axis associated with each particle.
	+ for (unsigned int i=0; i<_N; i++) {
	+
	+ PseudoJet myParticle = _myParticles[i];
	+ FunctorJetScalarPt sumPt = FunctorJetScalarPt();
	+
	+ //Recombiner has to carry user_index information.
	+ FunctorJetAxis axis(new fastjet::JetDefinition(fastjet::cambridge_algorithm, 2.0*_Rjet, new WinnerTakeAllRecombiner(), fastjet::Best));

	- FunctorJetScalarPt sumPt = FunctorJetScalarPt();
	- //Recombiner has to carry user_index information.
	- FunctorJetAxis axis(new fastjet::JetDefinition(fastjet::cambridge_algorithm, 2*_Rjet, new WinnerTakeAll(), fastjet::Best));
	-
	- fastjet::Selector sel = SelectorCircle(_Rjet);
	- sel.set_reference(_myParticles[i]);
	- vector<PseudoJet> nearbyParticles = sel(_myParticles);
	-
	- int myAxis = -1;
	- double pt_w = 0;
	- double Njet_w = 0;
	-
	- double pt_in_Rjet = sumPt(nearbyParticles);
	- if(pt_in_Rjet >= _ptcut) {
	- myAxis = axis(nearbyParticles)[0].user_index();
	- pt_w = _myParticles[i].pt();
	- Njet_w = pt_w / pt_in_Rjet;
	- }
	-
	- _axes.push_back(myAxis);
	- _pt_weights.push_back(pt_w);
	- _Njet_weights.push_back(Njet_w);
	+ fastjet::Selector sel = SelectorCircle(_Rjet);
	+ sel.set_reference(myParticle);
	+ vector<PseudoJet> nearbyParticles;

	+ if (_useStorageArray) { // start from phi/rapidity strips
	+ if (!_myStorageArray.aboveCutFor(myParticle)) {
	+ nearbyParticles.resize(0); //don't do any further analysis if not above jet pTcut
	+ } else {
	+ nearbyParticles = sel(_myStorageArray.getStorageFor(myParticle));
	+ }
	+ } else { // use all particles
	+ nearbyParticles = sel(_myParticles);
	+ }
	+
	+ int myAxis = -1;
	+ double pt_w = 0;
	+ double Njet_w = 0;
	+
	+ double pt_in_Rjet = sumPt(nearbyParticles);
	+ if(pt_in_Rjet >= _ptcut) {
	+ myAxis = axis(nearbyParticles).user_index();
	+ pt_w = _myParticles[i].pt();
	+ Njet_w = pt_w / pt_in_Rjet;
	+ }
	+
	+ _axes.push_back(myAxis);
	+ _pt_weights.push_back(pt_w);
	+ _Njet_weights.push_back(Njet_w);
	+
	}
	}
	-
	+
	void EventShapeDensity_JetAxes::find_axes_and_weights() const {


	//If requested establish global consistency.
	- if(_applyGlobalConsistency) {
	+ if (_applyGlobalConsistency) {
	for(unsigned int i=0; i < _N; i++) {
	- while(_axes[i]!= -1 && !_isStable(_axes[i])) _axes[i] = _axes[_axes[i]];
	+ while(_axes[i]!= -1 && !_isStable(_axes[i])) _axes[i] = _axes[_axes[i]];
	}
	}


	//Find distinct axes.
	//Output a vector of Pseudojets (sorted by pT), with pT of each axis corresponding to the pT weight.
	vector<double> tot_Njet_weights(_N,0), tot_pt_weights(_N,0);

	- for(unsigned int i=0; i<_N; i++) {
	- if(_axes[i] != -1) {
	- tot_pt_weights[_axes[i]] += _pt_weights[i];
	- tot_Njet_weights[_axes[i]] += _Njet_weights[i];
	- }
	+ for (unsigned int i=0; i<_N; i++) {
	+ if (_axes[i] != -1) {
	+ tot_pt_weights[_axes[i]] += _pt_weights[i];
	+ tot_Njet_weights[_axes[i]] += _Njet_weights[i];
	+ }
	}

	- for(unsigned int i=0; i<_N; i++){
	- if(tot_pt_weights[i] > 0) {
	- PseudoJet myAxis = _lightLikeVersion(_myParticles[i]) * tot_pt_weights[i];
	- _distinctAxes.push_back(myAxis);
	- }
	+ for (unsigned int i=0; i<_N; i++){
	+ if(tot_pt_weights[i] > 0) {
	+ PseudoJet myAxis = _lightLikeVersion(_myParticles[i]) * tot_pt_weights[i];
	+ _distinctAxes.push_back(myAxis);
	+ }
	}

	_distinctAxes = sorted_by_pt(_distinctAxes);

	//Find the corresponding N_jet weights.
	for(unsigned int i=0; i<_distinctAxes.size(); i++) _tot_Njet_weights.push_back(tot_Njet_weights[_distinctAxes[i].user_index()]);

	}

	bool EventShapeDensity_JetAxes::_isStable(const int thisAxis) const {

	bool answer = false;
	if(_axes[thisAxis] == thisAxis) answer = true;
	return(answer);
	}


	} // namespace contrib

	FASTJET_END_NAMESPACE

File Metadata

Mime Type: text/x-diff
Expires: Wed, May 14, 11:38 AM (10 h, 3 m)
Storage Engine: blob
Storage Format: Raw Data
Storage Handle: 5111473
Default Alt Text: (74 KB)

No OneTemporaryActions

View Options

File Metadata

Event Timeline

No OneTemporary
Actions