#ifndef __FASTJET_FILTER_HH__
#define __FASTJET_FILTER_HH__

//STARTHEADER
// $Id: Filter.hh 1737 2010-07-08 16:19:56Z soyez $
//
// Copyright (c) 2009-2010, Matteo Cacciari, Gavin Salam and Gregory Soyez
//
//----------------------------------------------------------------------
// This file is part of FastJet.
//
//  FastJet 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.
//
//  The algorithms that underlie FastJet have required considerable
//  development and are described in hep-ph/0512210. If you use
//  FastJet as part of work towards a scientific publication, please
//  include a citation to the FastJet paper.
//
//  FastJet 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 FastJet; if not, write to the Free Software
//  Foundation, Inc.:
//      59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
//----------------------------------------------------------------------
//ENDHEADER

#include "fastjet/ClusterSequenceAreaBase.hh"
#include <cassert>
#include <cmath>
#include <iostream>
#include <limits>
#include <string>

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


//----------------------------------------------------------------------
/// class that helps perform filtering on Cambridge/Aachen jets, and
/// optionally subtraction (if rho > 0)
///
/// Use it to define a filter, and then use that filter with the
/// FilteredJet class. For example
///
///    // define a filter, on scale Rfilt, that keeps at least nfilt
///    // objects  and all those with pt > ptkeep
///    Filter filter(Rfilt, nfilt, ptkeep)
/// 
///    ....
///    ClusterSequence cs(...)
///    PseudoJet       jet(...)
///
///    // get filtered jet as follows (use mostly like normal jet)
///    FilteredJet fjet(cs, jet, filter);
///
/// Instead of specifying Rfilt (default to C/A algorithm), you can
/// also provide a JetDefinition of your choice.
///
/// Filtering with C/A (main jet alg) and C/A (subjet alg) can also
/// be done with subtraction by providing as extra parameter "rho"
/// to the constructor.
///
///
/// Credits
/// -------
///
/// Filtering with just "nfilt" was originally proposed for
/// boosted-object reconstruction in arXiv:0802.2470, and for normal
/// kinematic reconstruction (e.g. dijet mass peaks) in
/// arXiv:0810.1304.
///
/// Filtering with just "ptkeep" was proposed in arXiv:0912.1342 under
/// the name "trimming", though there "ptkeep" is expressed as a fraction
/// of a hard scale in the problem rather than in absolute terms.
///
///
class Filter {
public:
  /// trivial ctor
  /// Note: this is just for derived classes
  ///       a Filter initialised through this constructor will not work!
  Filter() :
    _forced_subdef(false){}

  /// default dtor
  virtual ~Filter(){}

  /// define a filter that decomposes the jet into subjets on angular
  /// scale Rfilt, extracts the nfilt hardest subjets, and any further
  /// ones that are above ptkeep. It uses the C/A algorithm for the
  /// subjets
  ///
  /// If rho is non-zero, and the original algorithm used 
  /// has areas with explicit ghosts, then by setting rho!=0, the 
  /// the filtering will be performed after noise subtraction on the subjets
  ///
  Filter(double Rfilt, int nfilt, 
	 double ptkeep = std::numeric_limits<double>::max(), double rho = 0.0) : 
    _Rfilt(Rfilt), _nfilt(nfilt), _ptkeep(ptkeep), _rho(rho), _forced_subdef(false) {}

  /// define a filter that decomposes the jet into subjets using the
  /// jet definition subjet_def, extracts the nfilt hardest subjets,
  /// and any further ones that are above ptkeep.
  ///
  /// Note: compared to the previous case, we have no support for
  ///       internal subtraction currently here
  Filter(const JetDefinition & subjet_def, int nfilt, 
	 double ptkeep = std::numeric_limits<double>::max(), double rho = 0.0) : 
    _subjet_def(subjet_def), _nfilt(nfilt), _ptkeep(ptkeep), _rho(rho), _forced_subdef(true) {}

  /// runs the filtering and sets kept and rejected to be the jets of interest
  /// (with non-zero rho, they will have been subtracted).
  /// We set it virtual so that derived classes can implement other 
  /// strategies to select the kept jets
  virtual void run_filter(const ClusterSequence & cs, const std::vector<PseudoJet> & jets, 
			  std::vector<PseudoJet> & kept, std::vector<PseudoJet> & rejected, 
			  ClusterSequence * &internal_cs) const;

  std::string description() const;

protected:
  /// return the dcut needed to see the requested _Rfilt
  double dcut (const ClusterSequence & cs) const {
    assert(cs.jet_def().jet_algorithm() == cambridge_algorithm 
	   || ((_Rfilt >= cs.jet_def().R()) && (_nfilt == 1)));
    double rtdcut = _Rfilt / cs.jet_def().R();
    return rtdcut*rtdcut;
  }

  /// sets filtered_elements to be all the subjets on which filtering will work
  /// [NB: this routine is work in progress as part of a transition to a Filter
  ///      that also works on jet collections]
  void set_filtered_elements(const ClusterSequence & cs, const std::vector<PseudoJet> & jets,
			     std::vector<PseudoJet> & filtered_elements,
                             ClusterSequence * &internal_cs) const;

  JetDefinition _subjet_def; ///< the jet def for subclustering (if not C/A)
  double _Rfilt;             ///< the C/A radius for subclustering
  int    _nfilt;             ///< the number of subjets to keeo
  double _ptkeep;            ///< also keep subjets above that pt
  double _rho;               ///< the background density (used for subtraction when possible)
  bool _forced_subdef;       ///< true if we have explicitly set _subjet_def
};



//----------------------------------------------------------------------
/// Class to contain a filtered jet.
///
/// It provides all the functionality of a PseudoJet (and is derived
/// from it), but the momentum is that of the jet after filtering
/// (including optional subtraction).
/// 
class FilteredJet : public PseudoJet {

public:
  /// constructor for a filtered jet
  /// if the filter requires subtraction, cs must be
  /// dynamic_castable to a ClusterSequenceAreaBase that uses
  /// explicit ghosts (and both original and filtering alg must be C/A)
  FilteredJet(const ClusterSequence & cs, const PseudoJet & jet, const Filter & filter) {
    _original_jets.push_back(jet);
    _internal_cs = NULL;
    run(cs, filter);
  }

  /// constructor for a filtered jet
  /// if the filter requires subtraction, cs must be
  /// dynamic_castable to a ClusterSequenceAreaBase that uses
  /// explicit ghosts.
  ///
  /// This provides filtering on the collection of input jets as a 
  /// __single entity__ (e.g. for boosted Higgs reconstruction,
  /// you can pass the two prongs of the Higgs and it will 
  /// produce a single jet that is the filtered Higgs)
  FilteredJet(const ClusterSequence & cs, const std::vector<PseudoJet> & jets, const Filter & filter) {
    _original_jets.clear();
    _internal_cs = NULL;
    copy(jets.begin(), jets.end(), back_inserter(_original_jets));
    run(cs, filter);
  }

  /// dtor
  ~FilteredJet(){
    if (_internal_cs != NULL) delete _internal_cs;
  }
  
  /// returns the original jet (the first of the original jets
  /// if you filtered a collection of jets)
  const PseudoJet & original() const {return _original_jets[0];}

  /// returns a vector of original jets
  const std::vector<PseudoJet> & original_jets() const {return _original_jets;}

  /// returns the subjets that were kept during the filtering procedure
  /// (subtracted if the filter requests it, and valid in the original cs)
  const std::vector<PseudoJet> & kept()     const {return _kept;}

  /// returns the subjets that were not kept during the filtering procedure
  /// (subtracted if the filter requests it, and valid in the original cs)
  const std::vector<PseudoJet> & rejected() const {return _rejected;}

  /// get the jet constituents
  std::vector<PseudoJet> constituents() const{
    std::vector<PseudoJet> constits;

    if (_internal_cs == NULL){
      // !!! watch out !!!
      /// this is going to crash if the mother cs goes out of scope
      for (std::vector<PseudoJet>::const_iterator sit = _kept.begin(); sit != _kept.end(); sit++){
	std::vector<PseudoJet> c = _original_cs->constituents(*sit);
	copy(c.begin(), c.end(), back_inserter(constits));
      }
    } else {
      for (std::vector<PseudoJet>::const_iterator sit = _kept.begin(); sit != _kept.end(); sit++){
	std::vector<PseudoJet> c = _internal_cs->constituents(*sit);
	copy(c.begin(), c.end(), back_inserter(constits));
      }
    }

    return constits;
  }
  
private:

  void run(const ClusterSequence & cs, const Filter & filter) {    
    _original_cs = &cs;

    // run the filter to get the subjets that are kept and/or rejected
    filter.run_filter(cs, _original_jets, _kept,_rejected, _internal_cs);

    // compute the final filtered jet
    reset(_kept[0]);
    for (unsigned i = 1; i < _kept.size(); i++) {
      cs.jet_def().recombiner()->plus_equal(*this, _kept[i]);
    }

    // make sure the subtracted jet has the same index (cluster and user)
    // (i.e. "looks like") the original jet
    if ((_internal_cs == NULL) && (_original_jets.size() == 1)) {
      set_cluster_hist_index(_original_jets[0].cluster_hist_index());
      set_user_index(_original_jets[0].user_index());
    } else {
      // give it an invalid cluster_hist_index and leave the user index to
      // whatever it is...
      set_cluster_hist_index(-1);
    }
  }

  std::vector<PseudoJet> _original_jets;
  std::vector<PseudoJet> _kept, _rejected;
  ClusterSequence *_internal_cs;
  const ClusterSequence *_original_cs; ///< yes, this is ugly!!
};


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

#endif