// 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
//
//////


//----------------------------------------------------------------------
// MyFunctionOfVectorOfPseudoJets
// 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 could be removed and the classes below (like JetLikeEventShape,
// FunctorJetCount, etc.) could derive 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
//
//-GPS: is there any reason to call these Functors? A Functor has a
//      quite specific mathematical meaning, and not all aspects are
//      relevant here, e.g. "associates to each morphism
//      f:X\rightarrow Y \in C a morphism F(f):F(X) \rightarrow F(Y)
//      \in D "
//----------------------------------------------------------------------

//////
//
// Example Functors to Use
//
//////

//----------------------------------------------------------------------
// FunctorJetCount
// Counts the number of jets (i.e. 1 per jet)
//-GPS: how does this count particles? It always returns 1. This seems a bit too trivial...
class FunctorJetCount : public MyFunctionOfVectorOfPseudoJets<double> {

   double result(const std::vector<PseudoJet> & particles) const { return 1.0; }
   
   std::string description() const { return "Jet Count";}
};


//----------------------------------------------------------------------
// FunctorJetScalarPt
// Finds the summed pt of the jet constituents (NOT the pT of the jet)
//-GPS: should this be called ScalarPtSum? It doesn't necessarily involve 
//      any jets at all, but it does do a sum.
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";}
};

//----------------------------------------------------------------------
// FunctorJetScalarPtToN
// Raises FunctorJetScalarTransverseMomentum to an arbitrary Power
//-GPS: similar comment as for FunctorJetScalarPt
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();
      }
};


//----------------------------------------------------------------------
// FunctorJetMass
// Finds the mass of the jet
//-GPS: similar comment as for FunctorJetScalarPt
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";}
};

//----------------------------------------------------------------------
// FunctorJetMassSquared
// Finds the squared mass of the jet
//-GPS: similar comment as for FunctorJetScalarPt
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 (to reduce computation time)
//
//////


//----------------------------------------------------------------------
// JWJStorageArray
// This class divides an event into (overlapping) regions of size 2R x 2R
// and caches the partition for later use.  Partitions that are below
// the ptcut can be ignored entirely.  This saves some computational time by
// avoiding having to calculate irrelevant distances.
class JWJStorageArray {

private:

   double _Rjet;
   double _ptcut;

   // storage for partitions of particles.
   std::vector<std::vector<std::vector<fastjet::PseudoJet> > > _regionStorage; 
   
   // bool for whether region needs to be considered.
   std::vector<std::vector<bool> > _aboveCutBool;

   double _rapmax;

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

   // helper functions to determine which region a pseudo jet is in   
   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;
   } 

   // make the storage array
   void establishStorage(const std::vector<fastjet::PseudoJet> & my_jets, double Rjet, double ptcut);
   
   // give the 2R x 2R region relevant for a given particle
   std::vector<fastjet::PseudoJet> & getStorageFor(const fastjet::PseudoJet & currentJet);
   
   // give whether the 2R x 2R region has enough pT to merit further examination 
   bool aboveCutFor(const fastjet::PseudoJet & currentJet);

};

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


//----------------------------------------------------------------------
// JetLikeEventShape
// Generic JetLikeEventShape which takes a Functor (derived from 
// MyFunctionOfVectorOfPseudoJets) as an argument and then sums over all
// "jets" (i.e. regions of radius R around each particle) 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)
//-GPS: why is this "JetLikeEventShape"? The result isn't a jet. You
//      mean an object that acts like an event shape as applied to jets?
class JetLikeEventShape : public MyFunctionOfVectorOfPseudoJets<double> {

   protected:
      MyFunctionOfVectorOfPseudoJets<double> * _measurement; //jet measurement to use

      double _Rjet, _ptcut, _Rsub, _fcut;
      bool _trim;
      
      //-GPS if at all possible, it's good to avoid mutables for temporary storage like
      //-GPS this, because it makes it impossible to run the code in parallel on
      //-GPS two separate threads. You may not think you need to worry about that now,
      //-GPS but it can create debugging nightmares later.

      // used for temporary storage of information when establishing weights
      mutable bool _includeParticle;
      mutable double _weightParticle, _pt_in_Rjet, _pt_in_Rsub;
      mutable std::vector<PseudoJet> _nearbyParticles;
      
      // establishes the weights for each jet region to avoid double-counting.
      void establishWeights(const PseudoJet myParticle, const std::vector<PseudoJet> & particles) const;
      
      // optional storage array to cache the 2R x 2R partitions.
      bool _useStorageArray;
      mutable JWJStorageArray _myStorageArray;
      void establishStorage(const std::vector<PseudoJet> & particles) const;

   public:
   
      JetLikeEventShape(): _measurement(NULL), _useStorageArray(true), _myStorageArray() {}
      
      // constructor without trimming
  //-GPS: comments here are a bit thin! In particular, it's useful to mention that measurement
  //      will be automatically deleted when the class is destroyed.
      JetLikeEventShape(MyFunctionOfVectorOfPseudoJets<double> * measurement, double Rjet, double ptcut)
         : _measurement(measurement), _Rjet(Rjet), _ptcut(ptcut), _Rsub(Rjet), _fcut(1.0), _trim(false), _useStorageArray(true), _myStorageArray() {}
      
      // constructor with trimming
      JetLikeEventShape(MyFunctionOfVectorOfPseudoJets<double> * measurement, double Rjet, double ptcut, double Rsub, double fcut)
         : _measurement(measurement), _Rjet(Rjet), _ptcut(ptcut), _Rsub(Rsub), _fcut(fcut), _trim(true), _useStorageArray(true), _myStorageArray() {}

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

      // the event shape value
      virtual double result(const std::vector<PseudoJet> & particles) const;
      
      // Default is to use storage array.  Can use this to turn off.
      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
//
//////

//----------------------------------------------------------------------
// ShapeJetMultiplicity
// 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(){}
};


//----------------------------------------------------------------------
// ShapeScalarPt
// Defines event shape 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(){}
};

//----------------------------------------------------------------------
// ShapeScalarPtToN
// Sum of pT^n of jets 
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(){}
};


//----------------------------------------------------------------------
// ShapeSummedMass
// 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(){}
};

//----------------------------------------------------------------------
// ShapeSummedMassSquared
// 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
// ShapeMissingPt is pT of the vector sum of jets momenta.
// Needs special coding because it requires measuring the magnitude of the
// transverse momentum vector
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
// A counter for subjets within trimmed fat jets.  Requires special coding
// to avoid too many repeated calculations
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
//
//////

//----------------------------------------------------------------------
// SelectorShapeTrimming
// Event-wide trimming is implemented as a selector.
Selector SelectorShapeTrimming(double Rjet, double ptcut, double Rsub, double fcut);


//----------------------------------------------------------------------
// SelectorJetShapeTrimming
// Jet-based trimming, implimented as a selector.
// It takes the jet pT as the sum of the given four-vectors
Selector SelectorJetShapeTrimming(double Rsub, double fcut);

//----------------------------------------------------------------------
// JetShapeTrimmer
// For convenience, a Transformer version of jet-shape trimming which 
// acts on individual jets.  Simply a wrapper for SelectorJetShapeTrimming
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
//
//////

//----------------------------------------------------------------------
// JetLikeEventShape_VariablePtCut
// Generic class to get the whole range of event shape values as ptcut 
// is changed.  Currently, this does not use the JWJStorageArray for
// speed up, but that can be added upon request.
class JetLikeEventShape_VariablePtCut {

protected:

   // The jet measurement of interest
   MyFunctionOfVectorOfPseudoJets<double> * _measurement;

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

   // Build up the event shape as a function of ptcut.
   void _buildStepFunction(const std::vector<PseudoJet> particles) const; 
  
   // temporary storage of information
   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) {}

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

   // constructor with trimming
   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();
      }
      
};

//----------------------------------------------------------------------
// ShapeJetMultiplicity_VariablePtCut
// Example of an event shape with a variable ptcut.  Here, we are just doing
// jet multiplicity, but the same could be done for any of the other event shapes
class ShapeJetMultiplicity_VariablePtCut: public JetLikeEventShape_VariablePtCut {

public:

   // As per the physics paper, the 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
//
//////

//----------------------------------------------------------------------
// ShapeJetMultiplicity_VariableR
// This is the only event shape we have coded up with the ability to do
// a variable R. Currently, this does not use the JWJStorageArray for
// speed up, but that can be added upon request.
//
//-GPS it's not clear from this what the class really does - is it
//-Variable R in the sense of the Variable R paper from a couple of
//-years ago, or is it variable R in the sense of being able to handle
//-multiple R values, any R value?
//
//-You should also signal that this uses N^2 memory (16 N^2 bytes as
//-far as I can tell), because that can become problematic - e.g. if
//-someone in HI tries this, with 30k particles, they'll get 14GB usage. 
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() {}

   // constructor without trimming
   ShapeJetMultiplicity_VariableR(double ptcut): _ptcut(ptcut), _Rsub(0.0), _fcut(1.0), _trim(false) {}
   
   // constructor with trimming
   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.
   //
   //-GPS: these are both huge arrays (N^2); it's much better to return 
   //-a constant reference to the internal array, so that time isn't
   //-wasted copying them. Also, do they really need to be virtual?
   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();
      }
      
};

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


//--------------------------------------------------------
// FunctorJetAxis
// Functor to find jet axis according to a given jet definition.  Needed
// for the event shape densities
class FunctorJetAxis : public MyFunctionOfVectorOfPseudoJets< PseudoJet > {

public:

   FunctorJetAxis(){}
   FunctorJetAxis(fastjet::JetDefinition &jetDef) : _jetDef(jetDef) {} 
   ~FunctorJetAxis(){}

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

private:

   fastjet::JetDefinition _jetDef;

};

//--------------------------------------------------------
// FunctorLightLikeVersion
// Takes a pseudojet and returns a light-like version with no mass
// energy set to 1, but maintaining the input rapidity and azimuth
//-GPS: above name doesn't correspond to the actual class name!
//-GPS: it's not entirely clear what this class does from the name.
//-GPS: It's actually doing two things: making it light-like and giving it 
//      a unit pt. Should it really be LightLikeUnitPtVersion?
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);
      //-GPS: why take only user index? It might be better to 
      //      say p = jet and then use reset_momentum_PtYPhiM(...)
      p.set_user_index(jet.user_index()); // carry user index infomation 
      return p;
   }
  
  //-GPS: this might also be more explicit about the unit pt
   std::string description() const { return "Light-like version";}
   
};

//--------------------------------------------------------
// WinnerTakeAllRecombiner
// Winner-take-all recombiner.  Returns a recombined pseudojet whose pt
// is the sum of the input pts, but direction is the harder jet direction.
class WinnerTakeAllRecombiner : public  fastjet::JetDefinition::Recombiner {

public:

   WinnerTakeAllRecombiner() {}

   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;
     //-GPS: pt costs a sqrt in FJ. Better to cache the pt's, since they'll be reused below
      if (harder.pt() < softer.pt()) std::swap(harder,softer);

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

protected:

   FunctorLightLikeVersion lightLikeVersion;
  
};


//--------------------------------------------------------
// EventShapeDensity_JetAxes
// This is the hybrid probability density shape for jet axes position
// It calculates axis directions using the winner-take-all recombination,
// and also returns weights.
class EventShapeDensity_JetAxes {

public:

   EventShapeDensity_JetAxes(){}
  
   // using default cluster measure (anti-kT with winner-take-all recombination)
   EventShapeDensity_JetAxes(double Rjet, double ptcut, bool applyGlobalConsistency = false): 
      _Rjet(Rjet), _ptcut(ptcut),
      _jetDef(fastjet::JetDefinition(fastjet::antikt_algorithm, 2.0*Rjet, new WinnerTakeAllRecombiner(), fastjet::Best)),
      _applyGlobalConsistency(applyGlobalConsistency), _useStorageArray(true) {
         _jetDef.delete_recombiner_when_unused();
      }

   // Using arbitrary jet clustering procedure
   EventShapeDensity_JetAxes(double Rjet, double ptcut, const fastjet::JetAlgorithm &jetAlgo, bool applyGlobalConsistency = false):
	_Rjet(Rjet), _ptcut(ptcut), _jetDef(fastjet::JetDefinition(jetAlgo, 2.0*Rjet, new WinnerTakeAllRecombiner(), fastjet::Best)),
	_applyGlobalConsistency(applyGlobalConsistency), _useStorageArray(true) { _jetDef.delete_recombiner_when_unused();}
  
   ~EventShapeDensity_JetAxes(){}
   
   // Set input and find local axes	
   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;}

   //(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;}


   //Description
   std::string ParameterString() const {
         std::stringstream stream;
         stream << "R_jet=" << _Rjet << ", pT_cut=" << _ptcut;
	 stream << ". Local clustering with " << _jetDef.description()<<".";
         if (_applyGlobalConsistency) stream << " Global consistency condition on.";
         return stream.str();
   }


   std::string description() const { return "Hybrid event shape density for finding jet axes." + ParameterString();}


private:
 
   double _Rjet, _ptcut;
   mutable fastjet::JetDefinition _jetDef;
   bool _applyGlobalConsistency;
   
   void _find_local_axes(const std::vector<PseudoJet> & particles) const;
   bool _isStable(const int thisAxis) const;

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

   FunctorLightLikeVersion _lightLikeVersion;

   bool _useStorageArray;
   mutable JWJStorageArray _myStorageArray;

};


} // namespace contrib

FASTJET_END_NAMESPACE

#endif  // __FASTJET_CONTRIB_JETSWITHOUTJETS_HH__