//STARTHEADER
// $Id: Filter.cc 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 "Filter.hh"
#include <cassert>
#include <algorithm>
#include <fastjet/ClusterSequenceActiveAreaExplicitGhosts.hh>
#include <sstream>

using namespace std;


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

//----------------------------------------------------------------------
string Filter::description() const {
  ostringstream ostr;
  ostr << "Filter with subjet_def = ";
  if (_forced_subdef) {ostr << _subjet_def.description() << ", ";}
  else { ostr << "C/A algorithm with R = " << _Rfilt << ", " ;}
  ostr << "nfilt = " << _nfilt << ", ptkeep = " << _ptkeep << ", rho = " << _rho;
  return ostr.str();
}
 
//----------------------------------------------------------------------
void Filter::set_filtered_elements(const ClusterSequence & cs, 
                                   const vector<PseudoJet> & jets,
				   std::vector<PseudoJet> & filtered_elements,
                                   ClusterSequence* & internal_cs) const {

  assert(jets.size() > 0);

  // there are 4 (2x2) possibilities here according to whether C/A is
  // used for the native CS and the CS used to compute the subjets.
  // Each of this 4 cases is treated below.
  //
  // Only when both algorithms are C/A can we easily get the subjets
  // using exclusive jets. In all other cases, an explicit
  // reclustering of the constituents is necessary.
  if ( (cs.jet_def().jet_algorithm() != cambridge_algorithm) || (_forced_subdef) ){
    // create the jet definition for the subclustering
    // If we're to used C/A, use the same recombiner as parent alg
    JetDefinition subjet_def;
    if (_forced_subdef) {
      subjet_def = _subjet_def;
    } else {
      subjet_def = JetDefinition(cambridge_algorithm, _Rfilt);
      subjet_def.set_recombiner(cs.jet_def().recombiner());
    }

    // for a generic case, the fact that the native CS has explicit
    // ghosts should be sufficient?  one should be able to know if a
    // given particle is a ghost, then create a csaeg for the
    // subclustering and hence know the area of the subjets
    if (_rho != 0.0) {
      const ClusterSequenceAreaBase * csab = dynamic_cast<const ClusterSequenceAreaBase *>(&cs);

      // make sure we have all relevant information
      if (!csab) throw Error("Attempt to filter and subtract (non-zero rho) without area info in cs");
      if (!csab->has_explicit_ghosts())
	throw Error("Attempt to filter and subtract (non-zero rho) without explicit ghosts");

      // get the constituents of all the (sub)jets
      vector<PseudoJet> constituents;
      vector<PseudoJet> ghosts;
      vector<PseudoJet> local_constituents;
      for (unsigned i = 0; i < jets.size(); i++){
	local_constituents = cs.constituents(jets[i]);
	for (vector<PseudoJet>::iterator cit = local_constituents.begin(); cit != local_constituents.end(); cit++){
	  if (csab->is_pure_ghost(*cit))
	    ghosts.push_back(*cit);
	  else 
	    constituents.push_back(*cit);
	}
      }

      double ghost_area=(ghosts.size()==0) ? 0.0 : csab->area(ghosts[0]);

      // recluster them
      ClusterSequenceActiveAreaExplicitGhosts * subcs = new ClusterSequenceActiveAreaExplicitGhosts(constituents, subjet_def, ghosts, ghost_area);
      filtered_elements = subcs->subtracted_jets(_rho);
      internal_cs = subcs;
    } else {
      // get the constituents of all the (sub)jets
      vector<PseudoJet> constituents;
      for (unsigned i = 0; i < jets.size(); i++) cs.add_constituents(jets[i], constituents);

      // recluster them
      ClusterSequence *subcs = new ClusterSequence(constituents, subjet_def);
      filtered_elements = subcs->inclusive_jets();
      internal_cs = subcs;      
    }
    //used_original_cs = false;

  } else { // default case of C/A jets reclustered with C/A

    // get the subjets from the first jet
    filtered_elements = cs.exclusive_subjets(jets[0], dcut(cs));
    // and then copy across the subjets from the others
    for (unsigned i = 1; i < jets.size(); i++) {
      vector<PseudoJet> excl = cs.exclusive_subjets(jets[i], dcut(cs));
      copy(excl.begin(), excl.end(), back_inserter(filtered_elements));
    }
    internal_cs = NULL;
    //used_original_cs = true;

    // perform subtraction if requested and possible
    if (_rho != 0.0) {
      const ClusterSequenceAreaBase * csab = dynamic_cast<const ClusterSequenceAreaBase *>(&cs);
      
      // make sure we have all relevant information
      if (!csab) throw Error("Attempt to filter and subtract (non-zero rho) without area info in cs");
      if (!csab->has_explicit_ghosts())
	throw Error("Attempt to filter and subtract (non-zero rho) without explicit ghosts");
      
      // now subtract the jets
      for (unsigned i = 0; i < filtered_elements.size(); i++) {
	//std::cout << filtered_elements[i].perp() << " " << csab->area(filtered_elements[i]) << " ";
	filtered_elements[i] = csab->subtracted_jet(filtered_elements[i], _rho);
	//std::cout << filtered_elements[i].perp() << std::endl;
      }
    }

  }

  // order the filtered elements in pt
  filtered_elements = sorted_by_pt(filtered_elements);
}

//----------------------------------------------------------------------
/// return a vector of subjets, which are the ones that would be kept by the filtering
void Filter::run_filter(const ClusterSequence & cs, const vector<PseudoJet> & jets,
			std::vector<PseudoJet> & kept, std::vector<PseudoJet> & rejected, 
                        ClusterSequence * &internal_cs) const {

  set_filtered_elements(cs, jets, kept, internal_cs);
  
  rejected.clear(); // make sure it's clean

  unsigned nkept = _nfilt;
  for (unsigned i = _nfilt; i < kept.size(); i++) {
    if (kept[i].perp() > _ptkeep) {nkept++;}
    else {rejected.push_back(kept[i]);}
  }
  if (kept.size() > nkept) kept.resize(nkept);

}


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